Variational inequality approach to enforcing the

arXiv:1611.08758v2 [cs.CE] 4 Dec 2016

non-negative constraint for advection-diffusion
equations
An e-print of the paper is available on arXiv: 1611.08758.
Authored by

J. Chang
Graduate Student, University of Houston.

K. B. Nakshatrala
Department of Civil & Environmental Engineering
University of Houston, Houston, Texas 77204–4003.
phone: +1-713-743-4418, e-mail: knakshatrala@uh.edu
website: http://www.cive.uh.edu/faculty/nakshatrala

These figures depict the concentration profiles of the advection-diffusion equation under the Discontinuous Galerkin formulation (left) and the Discontinuous
Galerkin formulation with bounded constraints enforced through a variational inequality (right). Only the regions where concentrations meet the threshold of 0.5 and
above are shown.
2016
Computational & Applied Mechanics Laboratory

Variational inequality approach to enforcing the non-negative
constraint for advection-diffusion equations
J. Chang, and K. B. Nakshatrala
Department of Civil and Environmental Engineering, University of Houston.
Correspondence to: knakshatrala@uh.edu
Abstract. Predictive simulations are crucial for the success of many subsurface applications, and
it is highly desirable to obtain accurate non-negative solutions for transport equations in these
numerical simulations. To this end, optimization-based methodologies based on quadratic programming (QP) have been shown to be a viable approach to ensuring discrete maximum principles
and the non-negative constraint for anisotropic diffusion equations. In this paper, we propose a
computational framework based on the variational inequality (VI) which can also be used to enforce
important mathematical properties (e.g., maximum principles) and physical constraints (e.g., the
non-negative constraint). We demonstrate that this framework is not only applicable to diffusion
equations but also to non-symmetric advection-diffusion equations. An attractive feature of the
proposed framework is that it works with with any weak formulation for the advection-diffusion
equations, including single-field formulations, which are computationally attractive. A particular
emphasis is placed on the parallel and algorithmic performance of the VI approach across large-scale
and heterogeneous problems. It is also shown that QP and VI are equivalent under certain conditions. State-of-the-art QP and VI solvers available from the PETSc library are used on a variety
of steady-state 2D and 3D benchmarks, and a comparative study on the scalability between the
QP and VI solvers is presented. We then extend the proposed framework to transient problems by
simulating the miscible displacement of fluids in a heterogeneous porous medium and illustrate the
importance of enforcing maximum principles for these types of coupled problems. Our numerical
experiments indicate that VIs are indeed a viable approach for enforcing the maximum principles
and the non-negative constraint in a large-scale computing environment. Also provided are Firedrake project files as well as a discussion on the computer implementation to help facilitate readers
in understanding the proposed framework.

A list of abbreviations
ABC

CG
DG
GAL
GMRES
KSP

Arnold-Beltrami-Childress
Conjugate Gradient method
Discontinuous Galerkin
(Continuous) Galerkin
Generalized Minimal Residual method
Krylov subspace iterative solver

Key words and phrases. anisotropy; variational inequalities; quadratic programming; non-negative solutions;
maximum principles; parallel computing; advection-diffusion equations; miscible displacement.
1

MCP
MP
MPI
PETSc
QP
QP - TRON
SP
SUPG
TAO
VI
VI - SS
VI - RS
WF

Mixed complementarity problem
Minimization problem
Message Passing Interface
Portable Extensible Toolkit for Scientific Computation [Balay et al., 2016]
Quadratic programming
Trust region Newton method
Strong problem
Streamlined Upwind Petrov-Galerkin
Toolkit for Advanced Optimization [Munson et al., 2014]
Variational Inequality
Semi-smooth method
Reduced-space active-set method
Weak formulation
1. INTRODUCTION AND MOTIVATION

This paper presents a numerical methodology based on variational inequalities for anisotropic
diffusion and advection-diffusion equations that satisfies discrete maximum principles, meets the
non-negative constraint, and is well-suited for solving large-scale problems using parallel computing.
We now provide a motivation behind the current work, a discussion on prior works, and an outline
of our approach highlighting the significance of the work presented herein.
The diffusion and advection-diffusion equations are important partial differential equations
which are commonly used to model flow and transport of chemical species in porous media. Some
of the applications include subsurface remediation [EPA, 2004; Harp and Vesselinov, 2013; Heikoop
et al., 2014] and transport of radionuclides [Hammond and Lichtner, 2010; Genty and Potier, 2011].
Since these important problems are not analytically tractable, one needs to rely on predictive
numerical simulations. An important aspect in a predictive simulation of these equations is to
satisfy the non-negative constraint of concentration of chemical species.
Research efforts over the years have successfully created numerical models and discretization for
these equations, but they are not without their setbacks. For example, non-monotone discretizations
like the finite element method may result in spurious oscillations with high Péclet numbers. Other
common issues that may occur within highly heterogeneous and anisotropic diffusion-type equations
are violations of the maximum principle and the non-negative constraint [Ciarlet and Raviart, 1973;
Lipnikov et al., 2007; Liska and Shashkov, 2008; Genty and Potier, 2011]. Such numerical setbacks
can result in algorithmic failures or sharp fronts that may result in erroneous approximations of
reactive transport. Moreover, several important applications which require accurate predictive
capabilities of transport solvers are often large-scale and cannot be solved on a single computer. It
is important for numerical algorithms to not only ensure maximum principle but scale well with
respect to both problem size and computing concurrency. Obtaining numerical solutions within
a reasonable amount of time is the ultimate goal when selecting or designing algorithms that
are robust and can ensure non-negative concentrations for a wide range of subsurface transport
applications.
1.1. Prior works on non-negative formulations. The prior non-negative formulations can
be broadly classified into the following five categories:
2

(a) Reporting the violations: In [Payette et al., 2012], several cases of violations of the maximum
principle and the non-negative constraint have been showcased for different anisotropic diffusivity tensors. This paper also demonstrates that h- and p-refinements do not eliminate these
violations. The adverse effects due to violations of the non-negative constraint for non-linear
ecological models and chemically reacting flows have been illustrated in [Mudunuru et al.,
2015]. Neither of these papers have provided any fix to overcome these violations.
(b) Mesh restrictions: The first work on maximum principles under the finite element method can
be traced back to the seminal paper by [Ciarlet and Raviart, 1973]. This paper considered
isotropic diffusion, and has shown that an acute-angled triangular mesh (which is a restriction
on the mesh) will satisfy the maximum principle under the finite element method. Anisotropic
diffusion equations have been addressed in [Huang and Wang, 2015], wherein they developed an
algorithm to generate metric-based meshes to satisfy the maximum principle for such equations.
[Mudunuru and Nakshatrala, 2016b] have addressed various versions of maximum principles
for diffusion and advection-diffusion equations and studied the performance of metric-based
meshes for these equations. This paper highlighted the main deficiency of metric-based meshes,
which is the need to alter the mesh for different diffusivity tensors. A comprehensive list and
discussion of other prior works related to enforcing mesh restrictions to meet the maximum
principle and the non-negative constraint can also be found in [Mudunuru and Nakshatrala,
2016b].
(c) Developing or altering formulations in the continuum setting: Two works that fall under this
category are [Harari, 2004; Pal et al., 2016], both of which addressed transient transport
problems. [Harari, 2004] utilized a stabilized method that is available for Helmholtz-type
equations to construct a stabilized formulation for transient isotropic diffusion equations to
meet the maximum principle. This approach, as presented in [Harari, 2004], is applicable to
one-dimensional setting. [Pal et al., 2016] meets maximum principles for transient transport
equations by employing two techniques. They rewrote transient transport equations, which
are parabolic in nature, into modified Maxwell-Cattaneo equations, which are hyperbolic in
nature, and employed the space-time Discontinuous Galerkin approach.
(d) Non-finite element approaches: A finite volume-based approach, to enforce the non-negative
constraint, as proposed in [Potier, 2005], involves a non-linear iterative procedure to select
appropriate collocation points for cell concentrations. This technique has been refined by
several others including [Lipnikov et al., 2007] and [Sheng and Yuan, 2016]. Other similar
approaches include the mimetic finite difference method [da Veiga et al., 2014], which ensures
monotonicity and positivity. Since neither the finite difference nor finite volume methods are
based on weak formulations, these mentioned works cannot be trivially extended to the finite
element method.
(e) Optimization-based techniques at the discrete level : Several studies over the years [Nakshatrala and Valocchi, 2009; Nagarajan and Nakshatrala, 2011; Nakshatrala et al., 2013, 2016]
have focused on the development of optimization-based methodologies that enforce the maximum principle and the non-negative constraint for diffusion problems. An optimization-based
methodology based on the work of the aforementioned studies has been applied to enforce
maximum principles advection-diffusion equations [Mudunuru and Nakshatrala, 2016a]. By
reformulating the advection-diffusion problem as a mixed finite element formulation under the

3

least-squares formalism, one introduces flux variables into the problem. The discrete formulation is also symmetric and positive-definite, so one can easily apply both bounded constraints
and equality constraints to ensure non-negative solutions and local mass conservation respectively. It should be noted that one may also employ normal equations or the least-squares
approach to ensure that the minimization problem for non-symmetric problems is convex
[Demmel, 1997; Burdakov et al., 2012; Pal et al., 2016; Chang et al., 2016b]. All of these
studies have employed quadratic programming (QP) techniques to enforce the maximum principle on 2D academic problems, but the problems studied are small-scale and do not require
state-of-the-art Krylov subspace (KSP) iterative solvers and preconditioners. Moreover, it is extremely difficult to find solvers for least-squares or penalty-type problems. To this end, we are
interested in numerical formulations and solvers that are suitable for large-scale applications.
It has been shown recently that parallel optimization-based solvers can handle large-scale heterogeneous and anisotropic diffusion in [Chang et al., 2016a], so we want to extend the work
presented in that study to advection-diffusion equations.
1.2. Our approach and its salient features. The main contribution of this paper is to
present a finite element computational framework applicable to both diffusion and advectiondiffusion equations that meets the maximum principle and satisfies the non-negative constraint.
The framework is built by rewriting the weak formulation (WF) as a variational inequality (VI) [Ulbrich, 2011].
The field of VIs grew from a problem posed by Antonio Signorini [Signorini, 1933, 1959]. This
problem was later coined as “Signorini problem” by Gaetano Fichera, who was a student of Signorini. Fichera posed the problem more precisely and obtained a variational inequality corresponding to the problem using which he established existence and uniqueness of solutions [Fichera, 1964,
1965]. VIs have been employed to study contact problems [Kikuchi and Oden, 1988; Hlavacek
et al., 2012], obstacle problems [Rodrigues, 1987], elastoplastic problems [Hlavacek et al., 2012;
Han and Reddy, 2012] and other problems arising in mechanics and mathematics [Kinderlehrer
and Stampacchia, 2000]. If the bilinear form under the WF is symmetric, one can rewrite the WF
as a QP, which is a special case of VIs [Chipot, 2009]. Most of the existing single-field WFs for
advection-diffusion equations do not have symmetric bilinear forms, and hence one cannot construct
equivalent problems under QP. To the best of the authors’ knowledge, VIs have not be employed
to develop numerical formulations to satisfy maximum principles and the non-negative constraint
for anisotropic advection-diffusion equations.
The framework is particularly suited for large-scale problems, which is the case with many
practical subsurface applications. The proposed framework enjoys the following salient features:
(I) One can enforce bounded constraints for any transport problems that may be non-symmetric
or nonlinear.
(II) One can employ any numerical formulation, even a single-field formulation, for solving
advection-diffusion equations.
(III) One can leverage on existing high performance computing libraries and toolkits (e.g.,
solvers and preconditioners).
(IV) The framework is amenable for parallel computing, which will be illustrated using both
strong and weak scaling studies.
The rest of the paper is organized as follows. In Section 2, we present the boundary value
problem for steady-state diffusion and advection-diffusion equations. In Section 3, we present the
4

variational inequality (VI) and the various single-field weak formulations (WF) in the continuous
setting. In Section 4, we propose the computational framework in a discrete setting and discuss
in detail the specific solvers and implementation procedure. In Section 5, numerical results for
the steady-state governing equations under the proposed framework are shown, and we conduct
a thorough strong and weak-scaling study to demonstrate the parallel performance. In Section
6, we provide an extension of the proposed framework to transient problems and illustrate the
performance of miscible displacement in porous media, which is a coupled non-linear phenomenon.
Conclusions are drawn in Section 7. To facilitate the reader to be able to reproduce the results
given in this paper, sample Firedrake project files along with the discussion on the solution strategy
for large-scale Darcy equations are provided in Appendices A and B.
On the notational front, we denote all the continuum vectors by lowercase boldface unitalicized
letters, and the vectors in the discrete setting are denoted by lowercase boldface italic letters. We
denote all the continuum second-order tensors by boldface uppercase unitalicized letters, and all the
finite element matrices are denoted by uppercase boldface italicized letters. In this paper, repeated
indices do not imply summation. Other notational conventions are introduced as needed.
2. GOVERNING EQUATIONS IN THE CONTINUOUS SETTING
Let Ω ⊂ Rd be a bounded open domain, where ‘d’ denotes the number of spatial dimensions.
The boundary is denoted by ∂Ω = Ω − Ω, where a superposed bar denotes the set closure. We
denote the set of all k-times continuously differentiable functions on Ω by C k (Ω). We denote the set
of all functions in C 0 (Ω) that are continuous to the boundary by C 0 (Ω). A spatial point is denoted
by x ∈ Ω. The gradient and divergence operators with respect to x are, respectively, denoted by
b (x).
grad[·] and div[·]. The unit outward normal to boundary is denoted by n
Let c(x) denote the concentration field. The boundary is divided into two parts: ΓD and ΓN ,
such that ΓD ∪ ΓN = ∂Ω and ΓD ∩ ΓN = ∅. ΓD is that part of the boundary on which Dirichlet
boundary conditions are enforced (i.e., concentration is prescribed). ΓN is the part of the boundary
on which Neumann boundary conditions are enforced (i.e., flux is prescribed). When advection
is present, the Neumann boundary is further divided into inflow and outflow regions, which are
defined as follows:
n
o
N
b
:=
x
∈
Γ
v(x)
·
n
(x)
<
0
(2.1a)
ΓN
inflow
n
o
N
b (x) ≥ 0
ΓN
v(x) · n
(2.1b)
outflow := x ∈ Γ
For uniqueness of the solution under a steady-state response,
 we assume that concentration is
D
prescribed on a non-zero part of the boundary (i.e., meas Γ > 0).
2.1. Strong problems (SP). The strong problem (SP) for steady-state diffusion reads: Find
c(x) ∈ C 2 (Ω) ∩ C 0 (Ω) such that we have
−div[D(x)grad[c(x)]] = f (x)
c(x) = cp (x)
−b
n(x) · D(x)grad[c(x)] = q p (x)

in Ω

(2.2a)

on ΓD

(2.2b)

on ΓN

(2.2c)

and the SP for steady-state advection-diffusion reads: Find c(x) ∈ C 2 (Ω) ∩ C 0 (Ω) such that we
have
v(x) · grad[c(x)] − div[D(x)grad[c(x)]] = f (x)
5

in Ω

(2.3a)

c(x) = cp (x)

on ΓD

(2.3b)

b (x) · (v(x)c(x) − D(x)grad[c(x)]) = q p (x)
n

on ΓN
inflow

(2.3c)

on ΓN
outflow

(2.3d)

−b
n(x) · D(x)grad[c(x)] = q p (x)

where v(x) is the advective velocity, f (x) is the prescribed volumetric source/sink, cp (x) is the
prescribed concentration on the boundary, q p (x) is the prescribed flux on the boundary, and D(x)
is the second-order diffusivity tensor. The diffusivity tensor is assumed to be bounded and uniformly
elliptic. That is, there exist two constants 0 < ξ1 ≤ ξ2 < +∞ such that
ξ1 y · y ≤ y · D(x)y ≤ ξ2 y · y

∀y ∈ Rd

(2.4)

Moreover, the diffusivity tensor is assumed to be symmetric. That is,
D(x) = DT (x) ∀x ∈ Ω

(2.5)

A solution to SP is commonly referred to as a classical solution.
2.2. Maximum principle and the non-negative constraint. From the theory of partial
differential equations, it is well-known that a classical solution for the above mentioned SPs satisfies
maximum principles. For completeness, we provide below the statement of the classical maximum
principle of second-order elliptic partial differential equations with Dirichlet boundary conditions
on the entire boundary.
Theorem 1. (Classical maximum principle) If ΓD = ∂Ω, c(x) ∈ C 2 (Ω) ∩ C 0 (Ω) and f (x) ≤ 0,
then
max c(x) ≤ max cp (x)
x∈Ω

x∈∂Ω

Proof. A proof can be found in [Gilbarg and Trudinger, 2001].

(2.6)


A generalization of the classical maximum principle that is relevant to this paper is provided
in [Mudunuru and Nakshatrala, 2016a]. Specifically, they have extended the classical maximum
principle on four fronts: the regularity of the solution is relaxed to C 1 (Ω) ∩ C 0 (Ω), the regularity of
the volumetric source f (x) is relaxed to the space of square integrable functions, the boundary can
have both Dirichlet and Neumann boundary conditions, and the Neumann boundary conditions
are further divided into inflow and outflow (i.e., similar to equations (2.1a)–(2.1b)). For the sake of
brevity, we defer all interested readers to the suggested reference. Another property that is relevant
to this paper is non-negative solutions, which can be shown to be a special case of maximum
principles. In particular, the above maximum principle implies the following result:
Corollary 1. (Non-negative solutions) If ΓD = ∂Ω, c(x) ∈ C 2 (Ω) ∩ C 0 (Ω), cp (x) ≥ 0, and
f (x) ≥ 0, then
0 ≤ c(x)

∀x ∈ Ω

(2.7)

The central aim of this paper is to obtain numerical solutions to the above governing equations
(i.e., equations (2.2) and (2.3)) that respect maximum principles and the non-negative constraint.
The main task will then be to find an appropriate setting for numerical solutions. The finite
difference method directly discretizes the SP. However, under the finite element method, the SP is
rewritten as a WF, which is equivalent to the SP under some regularity assumptions. A solution to
a WF is referred to as a weak solution. As mentioned in Section 1 and will be shown using several
6

Regularity
of solution

Weak form
Solution space
is a closed
convex set

Strong form

Symmetry of
bilinear form

Minimization
problem

of
ry rm
t
e fo
m
m ear
y
S lin
bi

Variational
inequality

Figure 1. Relationships between the strong problem (SP), a weak formulation
(WF), a variational inequality (VI), and the minimization problem (MP).
examples later in this paper, a WF does not guarantee non-negative solutions in the discrete
setting. To overcome this deficiency, some non-negative formulations, especially for diffusion-type
equations, have rewritten the WF as an equivalent minimization problem (MP) and augmented with
bound constraints. However, it needs to be emphasized that such conversion is not always possible,
which is the case with the typical WF for advection-diffusion equations, as these formulations
have non-symmetric bilinear forms. In order to handle non-self-adjoint differential operators (e.g.,
advection-diffusion equation) and WFs with non-symmetric bilinear forms, we rewrite a given WF
as a VI. In order to satisfy maximum principles and the non-negative constraint, we restrict the
feasible solution space of the VI formulation using bound constraints. It needs to be mentioned
that one can pose the VI as an equivalent MP only if the bilinear form is symmetric. Figure 1
illustrates the various ways of rewriting the governing equations, and the conditions under which one
form is equivalent to the other. We now present various WFs for diffusion and advection-diffusion
equations, which will form the basis for our proposed VI-based formulations.
3. VARIATIONAL INEQUALITIES AND WEAK FORMULATIONS
The non-negative constraint and maximum principles restrict the feasible solution space to a
closed convex set. A variational inequality (VI) is basically a variational problem on a convex set,
which need not be a vector space. To this end, let C denote the solution space for the concentration
field, and K be a closed convex subset of C. The subset K is defined by the underlying maximum
principles and the non-negative constraint. The formulation based on VIs corresponding to the
mention SPs can be compactly written as: Find c(x) ∈ K such that we have
B(w − c; c) ≥ L(w − c)

∀w(x) ∈ K

(3.1)

where B(·; ·) is a bilinear form and L(·) is a linear functional, whose specific choices are provided
by the associated weak formulation. A WF can be abstractly written as: Find c(x) ∈ C such that
we have
B(w; c) = L(w)

∀w(x) ∈ W

(3.2)

where C and W are appropriate function spaces for a given WF. Our intention is to illustrate
the applicability of the proposed VI framework to a variety of WFs. To this end, we employ the
continuous Galerkin (GAL), Streamlined Upwind Petrov-Galerkin (SUPG), and Discontinuous Galerkin
7

(DG) formulations, which are documented below. For convenience, the standard L2 inner-product
over K is denoted as follows:
Z
(a; b)K =
a(x) · b(x) dK
(3.3)
K

3.1. Continuous Galerkin. The relevant function spaces are:
n
o
C := c(x) ∈ H 1 (Ω) c(x) = cp (x) on ΓD
n
o
W := w(x) ∈ H 1 (Ω) w(x) = 0 on ΓD

(3.4)
(3.5)

where H 1 (Ω) is a Sobolev space [Brezis, 2010]. We assume that f (x) ∈ H −1 (Ω), which is a dual
space corresponding to H 1 (Ω). We employ the GAL formulation for the diffusion problem, for
which the bilinear form and linear functional are:


(3.6)
BGAL (w; c) := grad[w(x)]; D(x)grad[c(x)]



Ω
LGAL (w) := w(x); f (x) − w(x); q p (x) N
(3.7)
Ω

Γ

For the advection-diffusion problem, spurious oscillations may arise under the GAL formulation for
high Péclet numbers. Herein, we employ the SUPG formulation [Brooks and Hughes, 1982], and
the corresponding bilinear form and linear functional are:




BSUPG (w; c) := BRES (w; c) + w(x); v(x) · grad[c(x)] + grad[w(x)]; D(x)grad[c(x)]
(3.8)
Ω
Ω




LSUPG (w) := LRES (w) + w(x); f (x) − w(x); q p (x) N
(3.9)
Ω

Γ

where the residual terms


h
BRES (w; c) :=
v(x) · grad[w(x)]; v(x) · grad[c(x)] − div [D(x)grad[c(x)]]
2kv(x)k
Ω


h
LRES (w) :=
v(x) · grad[w(x)]; f (x)
2kv(x)k
Ω

(3.10)
(3.11)

and h denotes the element-wise diameter.
3.2. Discontinuous Galerkin. For several transport applications, it is highly desirable to
possess element-wise mass balance property, as it is an important fundamental physical law [Turner
et al., 2011]. This is particularly true when the transport is coupled with chemical reactions and
biofilm growth [Werth and Valocchi, 2005; von der Schulenburg et al., 2009]. The GAL and SUPG
formulations do not possess this property without any further modification or enrichment to their
formulations. One way to ensure this property under a single-field finite element framework is
through the use of the DG formulations (see [Arnold et al., 2002; Riviére and Wheeler, 2002;
Cockburn, 2003; Li and Riviére, 2015b,a, 2016; Pal et al., 2016] and the references within for
further details). To present the DG formulation employed in the paper, we now introduce relevant
notation.
The domain Ω is divided into S subdomains:
Ω=

S
[
i=1
8

ωi

(3.12)

The boundary of the subdomain ωi is denoted by ∂ωi . The interior face between ωi and ωj is
denoted by Γij . That is,
Γij = ∂ωi ∩ ∂ωj

(3.13)

The set of all points on the interior faces is denoted by Γint . Mathematically,
S
[

Γint =

Γij

(3.14)

i=1,i<j

For an interior face, we denote the subdomains shared by this face by ω + and ω − . The outward
b + and n
b − . Employing
normals on this face for these subdomains are, respectively, denoted by n
Brezzi’s notation [Arnold et al., 2002], the average and jump operators on an interior face are
defined as follows:

c+ + c−
c :=
2

 
b + + c− n
b−
c := c+ n

and

(3.15)

where
c+ = c|∂ω+

and c− = c|∂ω−

(3.16)

One of the most popular DG formulations is the Interior Penalty method, which for equation
(2.2) is written as:




 
BDG (w; c) := grad[w(x)]; D(x)grad[c(x)] − w(x) ; D(x)grad[c(x)]
Ω
Γint



 

γ 
+  D(x)grad[w(x)] ; c(x)
+
(3.17)
w(x) ; c(x)
h
Γint
Γint




LDG (w) := w(x); f (x) − w(x); q p (x) N
(3.18)
Ω

Γ

where the penalty term γ = 2 (d+1)
[Shahbazi, 2005] for first-order elements and  ∈ [−1, 0, 1]
d
denotes the Symmetric, Incomplete, and Non-symmetric Interior Penalty methods respectively.
For equation (2.3), the DG formulation can be written as:




 
BDG (w; c) := grad[w(x)]; D(x)grad[c(x)] − w(x) ; D(x)grad[c(x)]
Ω
Γint



 

γ 
w(x) ; c(x)
+  D(x)grad[w(x)] ; c(x)
+
h
Γint
Γint




 up
− w(x); v(x) · grad[c(x)] − w(x) ; c (x)v(x)
(3.19)
Ω
Γint




LDG (w) := w(x); f (x) − w(x); q p (x) N
(3.20)
Ω

Γ

where the upwinding term cup (x) is defined as:
 +
b + (x) > 0
c (x) if v(x) · n
up
c (x) =
−
c (x) otherwise

(3.21)

For the remainder of this paper, we shall consider only the Symmetric Interior Penalty method
where  = −1.
9

3.3. A theoretical discussion. The bilinear form is assumed to be continuous (i.e., bounded
above). That is, there exists a constant κ1 > 0 such that
B(w; c) ≤ κ1 kckkwk

∀c(x), w(x) ∈ C

(3.22)

In addition, the bilinear form is assumed to be coercive. That is, there exists a constant κ2 > 0
such that
κ2 kck2 ≤ B(c; c)

∀c(x) ∈ C

(3.23)

Recall that L(·) is assumed to be a linear continuous functional on C. Then, from the Lax-Milgram
theorem [Brenner and Scott, 1994], it is known that a unique solution exists under the WF. Under
the same conditions on the bilinear form and the linear functional, a unique solution exists for the
associated VI if K ⊂ C is a closed convex subset [Lions and Stampacchia, 1967]. A solution of the
VI is a solution of the WF if C = K. Moreover, if the bilinear form is symmetric, that is,
B(w; c) = B(c; w)

(3.24)

then the WF and the VI are equivalent to the following MP: Find c(x) ∈ C such that
minimize
c(x)∈C

1
B(c; c) − L(c)
2

(3.25)

These relations are pictorially described in Figure 1. From a theoretical point of view, it is important
to note that the VIs that we will be dealing with for steady-state problems will be elliptic of first
kind. For further details on infinite-dimensional VIs, see [Glowinski, 1984; Duvaut and Lions,
1976].
4. PROPOSED COMPUTATIONAL FRAMEWORK IN A DISCRETE SETTING
We denote the total number of degrees-of-freedom by “ndof s”. We also denote the vector of
ones by 1, whose size will be apparent from the context in which it is used. The component-wise
inequalities are denoted by  and . That is,
ab

implies that

ai ≤ bi ∀i

(4.1a)

ab

implies that

ai ≥ bi ∀i

(4.1b)

The vector of unknown nodal concentrations is denoted by c, and the corresponding nodal source
vector is denoted by f . The coefficient matrix after a finite element discretization is denoted by K.
Note that the vectors c and f are of size ndof s × 1, and the matrix K is of size ndof s × ndof s.
We denote the standard inner-product in Euclidean spaces by h·; ·i. That is,
ha; bi =

ndof
Xs

ai bi

∀a, b ∈ Rndof s

(4.2)

i

The formulation in the discrete setting will be posed as a mixed complementarity problem (MCP)
[Kinderlehrer and Stampacchia, 2000]. For convenience, we define h ∈ Rndof s as:
h := Kc − f

(4.3)

The corresponding MCP reads: Find cmin 1  c  cmax 1 such that for each i ∈ {1, ..., ndof s}
hi (c) ≥ 0

if cmin = ci

(4.4a)

hi (c) = 0

if cmin ≤ ci ≤ cmax

(4.4b)

10

c − c∗

c

∗

−h(c∗)

normal cone

feasible region c
closed convex
set K
h(c∗)

Figure 2. This figure illustrates the condition under which a solution exists for
a variational inequality of the form hh(c), c̃ − ci ≥ 0 ∀c̃ ∈ K. Here, c∗ de∗
∗
notes
a solution of the VI. The normal
n
o cone of K at c is defined as N (c ) :=
w ∈ Rndof s hw; c − c∗ i ≤ 0 ∀c ∈ K .
hi (c) ≤ 0

if ci = cmax

(4.4c)

where cmin and cmax , respectively, denote the minimum and maximum concentrations, which are
provided by the maximum principle or the non-negative constraint. Simple complementarity conditions arise from the first-order optimality conditions in optimization. For bound-constrained
optimization, h corresponds to the gradient of the objective functional. If one has only the nonnegative constraints (i.e., cmin = 0 and cmax = +∞), then the problem reduces to a non-linear
complementarity problem, which is a special case of MCP. For details on non-linear complementarity problems, see [Facchinei and Pang, 2003]. Note that the feasible region, which is restricted
by the bound constraints, form a parallelepiped, which is a convex set [Boyd and Vandenberghe,
2004].
Let the feasible region K be a convex subset of Rndof s . In our case, the feasible region is
restricted by constraints which are in the form of finite number of linear equalities and inequalities.
This makes the feasible region to be a polyhedron, which is a convex set [Boyd and Vandenberghe,
2004]. It should be noted that bound constraints are a special case of linear inequalities. With this
machinery at our disposal, one can pose the second formulation based on variational inequalities,
which reads: Find c ∈ K such that we have
hKc; v − ci ≥ hf ; v − ci

∀v ∈ K

(4.5)

Note that MCP is a special case of VIs in which the feasible region is a parallelepiped (i.e., one has
only bound constraints). The conditions under which a solution exists for the finite-dimensional
VI given in equation (4.5) is pictorially described in Figure 2.
If the coefficient matrix K is symmetric, one can alternatively enforce maximum principles and
the non-negative constraint using QP, which has been illustrated in [Nagarajan and Nakshatrala,
2011; Nakshatrala et al., 2016] for small-scale problems, and in [Chang et al., 2016a] for largescale problems in parallel environments. Therefore, this approach is only applicable for formally
self-adjoint differential operators. The formulation can be posed as follows:
1
minimize
hc; Kci − hc; f i
(4.6a)
2
c∈Rndof s
11

c2
K

R
c1

0
BR (0)

KR = K ∩ BR (0)

Figure 3. A pictorial description of BR (0) and KR . These sets are used in Theorem 4.
subject to cmin 1  c  cmax 1

(4.6b)

In addition, if K is positive-definite the objective function becomes convex. The resulting optimization problem then belongs to the special case of convex quadratic programming for which
sophisticated solvers exist.
Remark 1. It should be mentioned that a quick fix to eliminate negative violations is through
the so-called clipping procedure. However, this procedure is rather ad hoc and, more importantly, it
is not variationally consistent. On the other hand, the proposed VI-based computational framework
not only ensures non-negative solutions but also has a firm variational basis. We will also illustrate
that the solutions under the proposed framework need not necessarily match the solution under the
clipping procedure.
4.1. Theoretical results in the discrete setting. In this paper, we are interested in problems with two different cases of bound constraints. In the first case, we have both lower and upper
bounds. In the second case, we have only the lower bound. The lower bound typically comes the
non-negative constraint, and the upper bound comes from maximum principles. We now discuss
existence results for finite-dimensional VIs under the mentioned two cases of bound constraints.
We begin by noting that the feasible set K will be convex and closed for both the sets of bound
constraints. In the first case, the feasible set will also be bounded, which makes the feasible set to
be compact (which, in the context of Euclidean spaces, is equivalent to closed and bounded). We
therefore deal with both the cases separately.
Theorem 2. (Existence based on compactness of K) If K is compact and convex, then a solution
exists to the finite-dimensional VI (4.5).
Proof. A proof can be constructed using the Brouwer’s fixed point theorem and can be found
in [Facchinei and Pang, 2003].

Theorem 3. (Existence based on positive-definiteness of sym[K]) If the symmetric part of the
coefficient matrix K (i.e., sym[K]) is positive-definite, a solution to the finite-dimensional VI (4.5)
exists. Note that the feasible set K need not be compact.
12

Proof. Let
g(c) := Kc − f

(4.7)

The VI then becomes
hg(c); e
c − ci

∀e
c∈K

(4.8)

Clearly, the function g(c) is continuous. Moreover, the function g(c) satisfies the following coercive
condition:
hg(c) − g(e
c); c − e
ci
→∞
(4.9)
e
kc − ck
as kck → ∞. To wit, since sym[K] is positive-definite and symmetric, the minimum eigenvalue
λmin is real and positive. One can then write:
λmin kc − e
ck2 ≤ (c − e
c) · sym[K] (c − e
c) = (c − e
c) · K (c − e
c) = hg(c) − g(e
c); c − e
ci

(4.10)

That is,
λmin kc − e
ck ≤

hg(c) − g(e
c); c − e
ci
kc − e
ck

(4.11)

We thus have shown that the function g(c) is continuous and coercive. Under such conditions, a
solution exists to the VI (4.8) (e.g., see [Nagurney, 2002; Nagurney and Zhang, 2012]).

We next present another existence theorem which is particularly useful when the feasible set is
unbounded (for example, when we have only one of the bounds – either lower or upper bounds).
Let BR (0) is a hypersphere of radius R centered at 0, and let KR = K ∩ BR (0) (see Figure 3).
Clearly, KR is bounded.
Theorem 4. (Existence based on KR ) A solution exists to the VI (4.5) on K (which need not
be bounded) if and only if there exists R > 0 and a solution c∗ ∈ KR that satisfies the following VI:
hKc∗ ; e
c − c∗ i

∀e
c ∈ KR

(4.12)

which is defined on a bounded set.
Proof. See [Nagurney, 2002].



Theorem 5. (Uniqueness) If the symmetric part of the coefficient matrix K is positive-definite,
then the finite-dimensional VI (4.5) has a unique solution if it exists.
Proof. On the contrary, assume that c1 and c2 are two different solutions of the VI (4.5).
This implies that
hKc1 ; v − c1 i ≥ hf ; v − c1 i

∀v ∈ K

(4.13)

hKc2 ; v − c2 i ≥ hf ; v − c2 i

∀v ∈ K

(4.14)

Since c1 , c2 ∈ K, choose v = c2 in (4.13) and v = c1 in (4.14). This results in
hKc1 ; c2 − c1 i ≥ hf ; c2 − c1 i

(4.15)

hKc2 ; c1 − c2 i ≥ hf ; c1 − c2 i

(4.16)

Summing the above two inequalities and invoking the linearity in the second slot, we obtain
hK(c1 − c2 ); c1 − c2 i ≤ 0
13

(4.17)

which further implies that
hsym[K](c1 − c2 ); c1 − c2 i ≤ 0

(4.18)

On the other hand, the positive-definiteness of sym[K] and our assumption c1 − c2 6= 0 imply that
hsym[K](c1 − c2 ); c1 − c2 i > 0
which contradicts the inequality given by equation (4.18). Hence, c1 = c2 .

(4.19)


Using the aforementioned general existence and uniqueness results for VIs, we now establish
the existence and uniqueness of solutions under the proposed framework in the discrete setting.
Theorem 6. (Well-posedness of the proposed framework) Unique solutions exist for the VIs
from the GAL, SUPG and DG WFs under lower bounds (which arise from the non-negative constraint) and under both lower and upper bounds (which arise from maximum principles).
Proof. First, it should be noted that the symmetric part of the coefficient matrices under the
GAL and DG formulations are positive-definite. On the other hand, the stabilization term under
the SUPG formulation does not guarantee that the symmetric part of the coefficient matrix to be
positive-definite. It should also be noted that the stabilization term in equation (3.10) is O(h),
where h is the characteristic mesh size. This means that there exist a critical mesh size, hcrit , such
that if h < hcrit then the contribution from the residual terms to the coefficient matrix will be
small, and the resulting symmetric part of the coefficient matrix will be positive-definite.
(Existence.) If both the lower and upper bounds are present, the feasible region will be compact.
For this case of bound constraints, Theorem 2 establishes the existence of solutions for the VIs
arising from all the three WFs (i.e., GAL, SUPG and DG). If only the lower bounds are present,
Theorem 3 will ensure the existence of solutions for VIs arising from the GAL and DG formulations,
and Theorem 4 will ensure the existence of solutions for the VIs arising from the SUPG formulation
on a general mesh. Of course, if the mesh is adequately refined (i.e., h < hcrit ) then Theorem 3 can
also ensure the existence of a solution for the VIs arising under the SUPG formulation.
(Uniqueness.) Theorem 5 provides the uniqueness of solution for the VIs arising from the GAL
and SUPG formulations. As discussed above, upon an adequate mesh refinement, sym[K] will
be positive-definite under the SUPG formulation. On those meshes, Theorem 5 will provide the
uniqueness of solutions for the VIs arising from the SUPG formulation.

4.2. Computer implementation details. In this paper, the proposed QP and VI-based
formulations for the GAL, SUPG, and DG formulations are implemented through the Firedrake
project (see Appendix A for further details), but one can employ any other finite element library.
The primary advantage of the Firedrake project is that it provides easy access to parallel solvers,
specifically the PETSc and TAO libraries [Balay et al., 2016, 1997; Munson et al., 2014] which are
built on top of Message Passing Interface (MPI). Appropriate iterative solvers and preconditioners
are needed for large-scale problems, and the PETSc library provides the necessary data structures.
The Conjugate Gradient (CG) method is used for symmetric problems like the diffusion equation
whereas the Generalized Minimal Residual (GMRES) method is used for the non-symmetric advectiondiffusion equation. HYPRE ’s algebraic multi-grid package [Law] is used as the preconditioner.
14

4.2.1. Solvers. The main ingredient of the proposed computational framework is to solve finitedimensional VIs. There are several solvers available for solving these type of inequalities in a
large-scale parallel environment. However, the performance of these solvers is problem-specific. It
is, therefore, necessary to identify the best performing VI solver for our case, which is primarily
to enforcing maximum principles and the non-negative constraint. To this end, we consider the
following VI and QP solvers available through the PETSc and TAO libraries:
(i) Semi-smooth (VI - SS): TAO’s implementation of the semi-smooth algorithm [Luca et al.,
1996; Munson et al., 2001] reformulates the MCP as a non-smooth system of equations using
the Fischer-Burmeister function [Fischer, 1992]. This function, φ : R2 → R, is defined as
p
φ(a, b) := a2 + b2 − a − b
(4.20a)
φ(a, b) = 0

⇔

a ≥ 0, b ≥ 0, ab = 0

(4.20b)

The reformulation of the MCP is handled component-wise, and the system of equations
Φ(c) = 0 where Φ : Rndof s → Rndof s is expressed as:


φ (ci − cmin , hi (c))
if − ∞ < cmin < cmax = ∞




φ
(c
−
c
,
−h
(c))
if
− ∞ = cmin < cmax < ∞

max
i
i
Φi (c) :=
(4.21)
φ (ci − cmin , φ (cmax − ci , −hi (c))) if − ∞ < cmin < cmax < ∞



−h
(c)
if
−
∞
=
c
<
c
=
∞
i
min
max


 c
if − ∞ < cmin = cmax < ∞
min − ci
It should be noted that Φ(c) is not differentiable everywhere but it still satisfies a semismoothness property [Mifflin, 1977; Qi, 1993; Qi and Sun, 1993]. The above system of equations is used to compute a descent direction, and the solver finishes as soon as the natural
merit function Ψ(c) := 21 kΦ(c)k22 meets some level of tolerance. We also employ TAO’s feasible line-search algorithm which ensures that the solution is within the bounds by using a
projected Armijo line search [Armijo, 1966].
(ii) Reduced-space active-set (VI - RS): The reduced-space active-set method selects an activeset and solves a reduced linear system of equations to calculate a direction of the gradient
descent. The active and inactive sets are, respectively, defined as:
A(c) := {i ∈ {1, ..., ndof s}

ci = 0 and hi (c) > 0}

(4.22a)

I(c) := {i ∈ {1, ..., ndof s}

ci > 0 or hi (c) ≤ 0}

(4.22b)

The active set A(c) represents regions where the lower bound is active thus the function value
can be ignored, and the inactive set I(c) contains everything else. The descent direction of
the active set is set to zero whereas the descent direction of the inactive set is approximated,
and the solution is updated using a projected line search. As far as we know, there is little
documentation on the theoretical and mathematical convergence properties for this particular
algorithm, but the computational results from [Benson and Munson, 2006] demonstrate that
this solver is robust and can handle a wide range of applications. For further implementation
details of these two VI solvers, we defer all interested readers to [Benson and Munson, 2006;
Munson et al., 2014] and the references within.
(iii) Trust region Newton (QP - TRON): Unlike the previous two solvers, the trust region Newton
method [Lin and Moré, 1999] is an active-set solver designed for large-scale minimization
problems. It uses the gradient projection to generate a Cauchy step and the preconditioned
15

CG with an incomplete Cholesky factorization to generate a descent direction. Each iteration
of the TRON algorithm solves a reduced linear system containing variables that lie between
the lower and upper bounds. The algorithm then applies a trust region to the conjugate
gradients to ensure convergence. The algorithmic scalability and hardware performance of
this solver has been thoroughly documented in [Chang et al., 2016a], so the computational
results arising from QP - TRON serves primarily as a benchmark for comparison with the
VI solvers.
We acknowledge that there may be several other QP and VI solvers which are not covered in this
paper. Nonetheless, the computational framework that we propose is algorithm-independent and
platform-agnostic, so one is free to either employ different solvers or modify the above implementations of the QP and VI algorithms to cater to specific needs and applications.
4.2.2. An outline of the algorithm. The performance of non-linear and optimization-based
solvers depends on accurate initial guesses. To this end, we propose the following steps for the
overall implementation of the proposed computational framework:
Step 1: Assemble K and f
Step 2: Solve for c0 .
Step 3: Clip c0 and obtain cCLIP . Formally:
1
cCLIP = arg min kc − c0 k2
2

subject to: cmin 1  c  cmax 1

(4.23)

Step 4: Solve the bounded constraint problem under the QP or VI framework with cCLIP as the
initial guess.
It should be emphasized that one need not solve equation (4.23) to implement clipping procedure.
Instead, one trims nodal values to meet the desired bounds. Since the governing equations are linear,
K and f only need to be assembled once and are reused for the various QP and VI evaluation
routines for Step 4. Python implementations of the VI - SS, VI - RS, and QP - TRON solvers
leveraging petsc4py [Dalcin et al., 2011] capabilities are shown in Listings 5, 6, and 7 respectively of
Appendix A. For the steady-state 3D benchmarks in the next section, the KSP relative tolerance
is set to 10−7 for the solver in Step 2, whereas the KSP relative tolerance for approximating the
gradient descents in Step 4 is set to 10−3 . It was shown in [Chang et al., 2016a] that relaxing
the tolerance requires more non-linear iterations but lessens the overall solve time. Relaxing the
tolerance also lessens the arithmetic intensity where the performance is governed by the memory
bandwidth thus making it less likely to achieve good speedup on a shared compute node. In other
words, the parallel efficiency of the QP and VI solvers are likely to be worse than solving the
WF with standard KSP convergence tolerances. The absolute convergence tolerances for the QP
and VI solvers are set to 10−8 although it should be mentioned that the optimal values depends
on the application at hand. All large-scale computations are conducted on Intel Xeon E5-2680v2
processors where each MPI process is restricted to a single core.
5. NUMERICAL RESULTS FOR STEADY-STATE RESPONSE
5.1. 2D benchmarks. We now examine 2D problems in order to demonstrate the effectiveness
of the proposed computational algorithms for ensuring discrete maximum principles and the nonnegative constraint. Only the GAL and SUPG formulations are employed in this numerical study.
First, let us consider the pure diffusion equation on a bi-unit square: Ω := (0, 1) × (0, 1) as shown
16

cP(x; y = 1) = 0

cP(x; y = 1) = 0

cP (x = 0; y) = 0

D(x)
f (x; y)

cPhole (x; y) = 1

cP (x = 0; y) = 0

cP (x = 0; y) = 0

cP (x = 0; y) = 0

D(x)
vx; vy
f (x; y) = 0

cP(x; y = 0) = 0

cP(x; y = 0) = 0

(a) Diffusion problem

(b) Advection-diffusion problem

Figure 4. 2D benchmarks: Pictorial description of the boundary value problems
for the diffusion and advection-diffusion examples. The finite element mesh for the
diffusion problem contains 40,000 quadrilateral elements and 40,401 nodes. The
mesh for the advection-diffusion problem is unstructured containing 96,430 triangular elements and 48,663 nodes.

in Figure 4a. The following heterogeneous and anisotropic diffusivity tensor similar to the one
considered in [Potier, 2005] is used:

D(x) =

y 2 + x2
−(1 − )xy
−(1 − )xy
x2 + y 2


(5.1)


 

where  = 10−4 . The forcing function is defined as f (x, y) = 1 if (x, y) ∈ 83 , 58 × 35 , 58 and
zero elsewhere. Homogeneous boundary conditions are applied on all four sides of the domain.
Numerical solutions under the GAL, VI - SS, VI - RS, and QP - TRON methods with uniform
quadrilateral elements of h-size = 1/200 are shown in Figure 5. All three non-negative solvers
successfully eliminate negative concentrations, and the absolute difference plots in Figure 6 show
that their results are quite different than from the one arising from the standard clipping procedure.
Moreover, the absolute differences between the various QP and VI solvers, as seen from Figure 7,
are extremely small and suggest that the QP and VI solvers have similar numerical accuracy.
Next we consider the advection-diffusion problem under the SUPG formulation where only VI SS and VI - RS methods are applicable for enforcing the maximum principle and the non-negative


constraint.
Consider a bi-unit square: Ω := (0, 1) × (0, 1) with a square hole of dimension 49 , 59 ×
4 5
9 , 9 as shown in Figure 4b. The mesh is discretized into 96,430 unstructured triangular elements
and 48,663 vertices. Homogeneous boundary conditions are applied on the outside boundary, and
a Dirichlet boundary value cp (x) = 1 is applied on the interior boundary Γhole . The velocity vector
17

(a) cGAL

(b) cTRON

(c) cSS

(d) cRS

Figure 5. 2D diffusion: concentrations under the Galerkin (cGAL ), TRON
(cTRON ), semi-smooth (cSS ), and reduced-space active-set (cRS ) methods where the
white regions represent negative concentrations.
field v(x) is characterized by the following:
vx = cos(2πy 2 )

(5.2a)
2

vy = sin(2πx) + cos(2πx )
18

(5.2b)

(a) kcCLIP − cTRON k

(b) kcCLIP − cSS k

(c) kcCLIP − cRS k

Figure 6. 2D diffusion: absolute difference in concentrations between the clipped
solution and the non-negative solution.
and the diffusivity tensor D(x) for this problem is the dispersion tensor:
D(x) = (αT kvk + DM ) I + (αL − αT )

v⊗v
kvk

(5.3)

where αL = 10−1 , αT = 10−5 , and DM = 10−9 denote the longitudinal dispersivity, transverse
dispersivity and molecular diffusivity, respectively. Figure 8 depicts the numerical solutions under
19

(a) kcTRON − cSS k

(b) kcTRON − cRS k

(c) kcSS − cRS k

Figure 7. 2D diffusion: absolute difference in concentrations between the various
non-negative methodologies

the SUPG, VI - SS, and VI - RS formulations. We see that the SUPG formulation results in negative
concentrations as well as concentrations greater than the maximum prescribed boundary condition
whereas the two VI solvers successfully correct these concentrations. The absolute difference plots,
as seen in Figure 9, indicate that the VI are also similar to one another and differ from the clipping
procedure. These 2D benchmarks suggest that the QP and VI solvers are accurate alternatives to the
20

(a) cSUPG

(b) cSS

(c) cRS

Figure 8. 2D advection-diffusion: concentrations under the SUPG (cSUPG ), semismooth (cSS ), and reduced-space active-set (cRS ) methods where the white regions
represent negative concentrations.
clipping procedure for satisfying the discrete maximum principle and the non-negative constraint.
Listings 1 and 2 contain the Firedrake project files for solving the GAL and SUPG formulations,
respectively.
21

(a) kcCLIP − cSS k

(b) kcCLIP − cRS k

(c) kcSS − cRS k

Figure 9. 2D advection-diffusion: absolute difference in concentrations between
the clipped and non-negative solutions.

5.2. 3D benchmark. We now consider a 3D problem designed to capture two particular
aspects that may arise in large-scale applications: 1) chaotic advection, which is pervasive in many
porous media applications [Lester et al., 2013], and 2) random point sources, which in subsurface
remediation problems are the sites where potential contaminant leaks occur. Predictive modeling
involving such important aspects require numerical methodologies that are not accurate but fast
22

cP(x; y; z = 1) = 0
cP(x = 0; y; z) = 0

cP(x; y = 0; z) = 0

D(x)
vx ; v y ; v z
f (x; y; z)

cP(x; y = 1; z) = 0

cP(x = 1; y; z) = 0
cP(x; y; z = 0) = 0

Figure 10. 3D benchmarks: Left figure contains a pictorial description of the
boundary value problem. Right figure contains the corresponding velocity contour
and vector field for the ABC flow (right) for the steady-state diffusion and advectiondiffusion problems.
and scalable in a parallel environment. Herein, our goal is to study the computational performance
of the various QP and VI solvers under the GAL, SUPG, and DG formulations.
Consider a unit cube domain as shown in Figure 10 with chaotic advection flow characterized
by the Arnold-Beltrami-Childress (ABC) flow [Zhao et al., 1993; Dombre et al.]:
vx = 0.3 sin(2πz) + cos(3πy)

(5.4a)

vy = 0.65 sin(2πx) + 0.3 cos(5πz)

(5.4b)

vz = sin(4πy) + 0.65 cos(6πy)

(5.4c)

For this problem, we shall also let D(x) denote the dispersion tensor as shown in equation
(5.3) where αL = 10−1 , αT = 10−5 , and DM = 10−9 . All six faces of the cube have homogeneous
boundary conditions, and the following forcing function consisting of 8 randomly located point
sources is used throughout the domain:

1 if (x, y, z) ∈ [0.4, 0.2, 0.1] × [0.5, 0.3, 0.2]




1 if (x, y, z) ∈ [0.8, 0.4, 0.2] × [0.9, 0.5, 0.3]





1 if (x, y, z) ∈ [0.5, 0.7, 0.3] × [0.6, 0.8, 0.4]




 1 if (x, y, z) ∈ [0.3, 0.5, 0.2] × [0.4, 0.6, 0.3]
f (x, y, z) =
(5.5)
1 if (x, y, z) ∈ [0.5, 0.2, 0.6] × [0.6, 0.3, 0.7]


 1 if (x, y, z) ∈ [0.6, 0.5, 0.7] × [0.7, 0.6, 0.8]




1 if (x, y, z) ∈ [0.4, 0.7, 0.8] × [0.5, 0.8, 0.9]




1 if (x, y, z) ∈ [0.1, 0.4, 0.7] × [0.2, 0.5, 0.8]



0 otherwise
23

(a) GAL

(b) GAL with VI - SS

(c) DG

(d) DG with VI - SS

Figure 11. 3D diffusion: 3D contours of the concentrations for the GAL and DG
formulations with and without and VI - SS for h-size = 1/80, where the purple
contours represent regions with negative concentrations (see online version for color
figures.)

To understand the parallel and algorithmic scalability of the QP and VI solvers, various levels
of mesh refinement are considered, ranging from 1,331 to 1,030,301 degrees-of-freedom for the
24

Table 1. 3D diffusion: minimum and maximum concentrations for various level of
mesh refinement under the GAL formulation.
h-size
1/10
1/20
1/30
1/40
1/50
1/60
1/70
1/80

Min. concentration
-0.0224497
-0.0071611
-0.0083804
-0.0062918
-0.0067679
-0.0072030
-0.0066007
-0.0059264

Max. concentration
0.368322
0.339679
0.481598
0.378390
0.477119
0.518469
0.498127
0.484484

% degrees-of-freedom violated
280/1,331 → 21.0%
2,462/9,261 → 26.6%
8,449/29,791 → 28.4%
20,195/68,921 → 29.3%
39,500/132,651 → 29.8%
68,161/226,981 → 30.0%
109,554/357,911 → 30.6%
160,925/531,441 → 30.3%

Table 2. 3D diffusion: minimum and maximum concentrations for various level of
mesh refinement under the DG formulation.
h-size
1/10
1/20
1/30
1/40
1/50
1/60
1/70
1/80

Min. concentration
-0.0226040
-0.0071913
-0.0082811
-0.0062341
-0.0067168
-0.0071682
-0.0065727
-0.0058998

Max. concentration
0.372831
0.341955
0.483264
0.379389
0.478146
0.519338
0.498775
0.485012

% degrees-of-freedom violated
3,704/8,000 → 46.3%
27,496/64,000 → 43.0%
91,176/216,000 → 42.2%
213,000/512,000 → 41.6%
410,976/1,000,000 → 41.1%
702,504/1,728,000 → 40.7%
1,114,856/2,744,000 → 40.6%
1,624,496/4,096,000 → 39.7%

GAL/SUPG formulations and ranging from 8,000 to 4,096,000 degrees-of-freedom for the DG formulations. Up to 16 MPI processes are used to study the weak-scaling and strong-scaling potential
of these solvers.
Figure 11 depicts the GAL and DG solutions for the diffusion equation with and without VI
- SS. It can be seen from the figures that negative concentrations are present regardless which
finite element formulation is used. Tables 1 and 2 indicate that negative concentrations arise for
the GAL and DG formulations, respectively, even as h-size is refined. It is interesting to note
that the DG formulation not only has more degrees-of-freedom but has more regions with negative
concentrations than its GAL counterpart. Using the initial guess solver from Step 2 of the proposed
framework in 4.2.2 as a baseline for comparison, Tables 3 and 4 demonstrate how the wall-clock time
and number of KSP/VI/QP solver iterations vary with h-refinement under a single MPI process. It
should be noted that the timings for the QP and VI solvers consider both the assembly of the data
structures as well as the actual solver. The heterogeneous nature of the problem causes the number
of solvers iterations to increase with problem size, but the iteration counts begin to stabilize as the
problem gets bigger. The VI - RS method outperforms VI - SS in both wall-clock time and VI
iterations but has similar performance to QP - TRON.
Next we perform weak-scaling studies to investigate how increasing both problem size and
number of MPI processes affects the performance of the VI and QP solvers. Each MPI process
will handle approximately 100k degrees-of-freedom so the h-sizes for the GAL case are 1/46, 1/58,
25

Table 3. 3D diffusion: wall-clock time and number of solver iterations (KSP, VI,
or QP) for various levels of mesh refinement under the GAL formulation.
h-size
1/10
1/20
1/30
1/40
1/50
1/60
1/70
1/80

GAL
time (s)
0.003
0.036
0.165
0.525
1.293
2.556
4.747
7.962

iters
9
15
20
24
28
31
35
39

VI - SS
time (s)
iters
0.027
5
0.477
12
2.624
18
7.576
20
21.49
27
43.72
30
76.21
31
140.7
37

VI - RS
time (s)
iters
0.008
2
0.147
5
0.765
7
2.246
8
5.381
9
12.01
11
18.27
10
36.40
13

QP - TRON
time (s)
iters
0.007
2
0.135
5
0.650
6
1.758
6
5.330
9
12.20
11
17.81
9
38.06
13

Table 4. 3D diffusion: wall-clock time and number of solver iterations (KSP, VI,
or QP) for various levels of mesh refinement under the DG formulation.
h-size
1/10
1/20
1/30
1/40
1/50
1/60
1/70
1/80

DG
time (s)
iters
0.030
10
0.446
14
2.148
18
6.278
21
14.29
24
28.25
27
51.95
31
85.53
34

VI - SS
time (s)
iters
0.748
12
1.410
20
73.52
27
251.5
34
623.6
39
1290
45
2560
51
5049
54

VI - RS
time (s)
iters
0.221
5
4.528
8
23.85
10
69.33
11
171.7
13
360.1
15
620.2
16
1107
19

QP - TRON
time (s)
iters
0.186
5
3.715
7
22.64
10
70.77
11
170.8
12
388.6
15
639.0
15
1291
17

1/73, 1/92, and 1/116 for 1, 2, 4, 8, and 16 processes respectively whereas the h-sizes for the DG
case are 1/23, 1/29, 1/37, 1/46, and 1/58 for 1, 2, 4, 8, and 16 processes respectively. Figure 12
contains the scaling plots as well as the parallel efficiencies in the weak sense under the GAL and DG
formulations. Generally speaking, the non-negative methodologies do not scale as well in the weak
sense for two reasons: 1) the wall-clock time and solver iteration for the QP and VI methodologies
are not linearly proportional to problem size (increasing solver iterations with h-size as seen from
Tables 3 and 4), and 2) the lower KSP relative tolerance for the gradient descent computations
make the overall solver more sensitive to the memory-bandwidth, meaning that speedup is reduced
as the compute nodes become populated with more MPI processes (see Sections 4 and 5 of [Chang
et al., 2016a] for a more thorough discussion).
However, the weak-scaling plots alone make it difficult to distinguish whether parallel performance deteriorates due to communication overhead or suboptimal algorithmic convergence. To
better understand why parallel performance degrades as the number of MPI processes increases,
we conduct strong-scaling studies by setting the h-size to 1/80 and 1/40 for the GAL and DG
formulations respectively (roughly 500k degrees-of-freedom) and study how increasing the number
of MPI processes (hence communication overhead) affects the parallel performance. Figure 13 contains the strong-scaling plots, and we see that the QP and VI solvers still do not scale as well.
26

100
GAL
VI - SS
VI - RS
QP - TRON

90
80

Efficiency (%)

Wall-clock time (s)

102

101

70
60
50
40
GAL
VI - SS
VI - RS
QP - TRON

30
20
100
1

2

4

8

10

16

1

2

# of cores

4

8

16

# of cores

(a) GAL - solve time

(b) GAL - parallel efficiency
100
90
80

Efficiency (%)

Wall-clock time (s)

10

DG
VI - SS
VI - RS
QP - TRON

2

101

70
60
50
40
DG
VI - SS
VI - RS
QP - TRON

30
100

20
1

2

4

8

10

16

# of cores

1

2

4

8

16

# of cores

(c) DG - solve time

(d) DG - parallel efficiency

Figure 12. 3D diffusion: weak-scaling plots with approximately 100k degrees-offreedom per core and the corresponding parallel efficiencies.

Regardless of the finite element formulation used, the QP and VI - RS methods have roughly the
same strong-scaling performance whereas the VI - SS method has slightly better strong-scaling.
For the advection-diffusion equation, the same problem is considering but advection due to the
ABC flow is now taken into account. A Firedrake project implementation of the DG formulation
can be found in Listing 3. Like the diffusion equation, the advection-diffusion equation also exhibits
negative concentrations as seen from Figure 14. Table 5 depicts violations under the SUPG formulation to be no greater than 30% whereas the DG formulation exhibits huge violations as seen from
Table 6. Moreover, the single MPI process metrics shown in Tables 7 and 8 clearly indicate that the
advection-diffusion equations are generally more expensive to solve than its diffusion counterpart.
These metrics tell us that the VI iterations also begin to stabilize as the problem size increases and
that the VI - RS method is faster than the VI - SS method for the advection-diffusion equation.
27

100
90
80

Efficiency (%)

Wall-clock time (s)

10

2

101

GAL
VI - SS
VI - RS
QP - TRON

100
1

2

70
60
50

GAL
VI - SS
VI - RS
QP - TRON

40

4

8

30

16

1

2

# of cores

4

8

16

# of cores

(a) GAL - solver time

(b) GAL - parallel efficiency
100
90

Efficiency (%)

Wall-clock time (s)

102

101

DG
VI - SS
VI - RS
QP - TRON

100

1

2

80
70
60
DG
VI - SS
VI - RS
QP - TRON

50

4

8

40

16

# of cores

1

2

4

8

16

# of cores

(c) DG - solver time

(d) DG - parallel efficiency

Figure 13. 3D diffusion: strong-scaling plots for approximately 500k degrees of
freedom (h-size = 1/80 and 1/40 for GAL and DG respectively) and the corresponding parallel efficiencies.
The weak-scaling plots, as seen from Figure 15, indicate that the VI solvers are identical in
performance, but the strong-scaling plots from Figure 16 suggest that VI - SS has slightly better
scaling than VI - RS. These steady-state numerical experiments suggest that the VI - RS is the
preferred methodology for solving large-scale advection-diffusion equations. However, if advection
becomes negligible thus reducing the system to a symmetric diffusion problem, one could use either
QP - TRON or VI - RS as these solvers are equal in both parallel and algorithmic performance.
Remark 2. These parallel performance studies do not indicate how truly efficient the PETSc
and TAO implementations of the QP and VI solvers are in the context of high performance computing. It should be noted that a serially efficient algorithm will likely have poor parallel efficiency due
to dominating effects from communication overhead and memory latencies. An approximation of
28

(a) SUPG

(b) SUPG with VI - SS

(c) DG

(d) DG with VI - SS

Figure 14. 3D advection-diffusion: 3D contours of the concentrations for the
SUPG and DG formulations with and without and VI - SS for h-size = 1/80, where
the purple contours represent regions with negative concentrations (see online version for color figures.)

the hardware performance for these solvers can be made through the perfect cache model described
in Section 4 of [Chang et al., 2016a]. However, a more thorough and detailed performance model
29

Table 5. 3D advection-diffusion: minimum and maximum concentrations for various level of mesh refinement under the SUPG formulation.
h-size
1/10
1/20
1/30
1/40
1/50
1/60
1/70
1/80

Min. concentration
-0.0135676
-0.0068733
-0.0091657
-0.0055686
-0.0064795
-0.0063168
-0.0053682
-0.0045065

Max. concentration
0.187489
0.180922
0.210942
0.171690
0.185440
0.189047
0.179675
0.172049

% degrees-of-freedom violated
212/1,331 → 15.9%
2,323/9,261 → 25.1%
7,964/29,791 → 26.7%
18,235/68,921 → 26.5%
35,221/132,651 → 26.6%
61,171/226,981 → 26.9%
99,668/357,911 → 27.8%
147,462/531,441 → 27.7%

Table 6. 3D advection-diffusion: minimum and maximum concentrations for various level of mesh refinement under the DG formulation.
h-size
1/10
1/20
1/30
1/40
1/50
1/60
1/70
1/80

Min. concentration
-0.0151514
-0.0162537
-0.0137824
-0.0067079
-0.0057574
-0.0065627
-0.0069389
-0.0066445

Max. concentration
0.259127
0.211295
0.237722
0.186956
0.203852
0.203093
0.193418
0.199912

% degrees-of-freedom violated
4,976/8,000 → 62.2%
37,464/64,000 → 58.5%
120,296/216,000 → 55.7%
276,832/512,000 → 54.1%
526,080/1,000,000 → 52.6%
891,768/1,728,000 → 51.6%
1,410,208/2,744,000 → 51.4%
2,069,752/4,096,000 → 50.5%

Table 7. 3D advection-diffusion: wall-clock time and number of linear (KSP) or
nonlinear (VI) solve iterations for various levels of mesh refinement under the SUPG
formulation.
h-size
1/10
1/20
1/30
1/40
1/50
1/60
1/70
1/80

SUPG
time (s)
0.003
0.045
0.221
0.768
2.019
4.530
8.178
14.70

iters
9
16
22
29
35
43
47
55

VI - SS
time (s)
iters
0.028
5
0.504
11
2.525
14
10.21
20
27.27
23
65.61
30
117.3
28
229.4
34

VI - RS
time (s)
iters
0.009
2
0.146
4
0.710
5
2.592
6
8.842
10
19.20
10
34.54
11
66.07
12

is needed to fully understand the performance of important memory-bandwidth limited algorithms
such as the ones studied in this paper and is part of our future work.
30

Table 8. 3D advection-diffusion: wall-clock time and number of linear (KSP) or
nonlinear (VI) solve iterations for various levels of mesh refinement under the DG
formulation.
h-size
1/10
1/20
1/30
1/40
1/50
1/60
1/70
1/80

DG
time (s)
0.031
0.459
2.352
6.819
15.49
30.65
57.56
97.15

iters
10
14
19
22
25
28
32
36

VI - SS
time (s)
iters
0.807
13
14.46
21
73.84
27
251.3
35
629.1
40
1387
47
2267
51
7105
60

VI - RS
time (s)
iters
0.202
5
3.890
8
19.56
10
56.94
11
137.1
13
292.3
15
570.0
18
917.2
18

6. EXTENSION TO TRANSIENT ANALYSIS AND COUPLED PROBLEMS
We now illustrate that the proposed framework, which is based on variational inequalities,
can be extended to perform a transient analysis. The resulting governing equations will then be
parabolic variational inequalities. This extension will be illustrated by considering the displacement
of miscible fluids in porous media wherein a fluid displaces a fluid with higher viscosity [Stalkup,
1983]. Some of the applications of miscible displacement include oil recovery and carbon-dioxide
sequestration [Chen and Meiburg, 1998a,b]. The phenomenon is commonly modeled using coupled
flow and transport equations, which will be presented below. In this section, we will also show
how negative concentrations can have serious ramifications when simulating non-linear transport
phenomenon like the displacement of miscible fluids.
6.1. Governing equations and temporal discretization. We denote the time by t ∈
[0, T ], where T denotes the length of the time interval of interest. For the Darcy equation, the
boundary is divided into two parts: Γp and Γv , such that Γp ∪ Γv = ∂Ω and Γp ∩ Γv = ∅. Γp
and Γv denote the parts of the boundary on which pressure and velocity boundary conditions are
enforced respectively. We shall denote time-dependent pressure by p(x, t), time-dependent velocity
by v(x, t), concentration-dependent viscosity by µ(c(x, t)), permeability by k(x), density by ρ, timedependent specific body force by b(x), time-dependent concentration by c(x, t), prescribed initial
concentration by c0 (x), time dependent volumetric source by f (x, t), and time-dependent diffusivity
tensor by D(x, t). For the boundary conditions, the prescribed time-dependent concentration is
denoted by cp (x, t), prescribed time-dependent pressure by pp (x, t), prescribed time-dependent
normal component of the velocity by vn (x, t), and prescribed time-dependent flux by q p (x, t). The
initial boundary value problem for the coupled flow and advective-diffusive equations can be written
as follows:
µ(c(x, t))
v(x, t) + grad[p(x, t)] = ρb(x, t)
in Ω × (0, T )
(6.1a)
k(x)
div[v(x, t)] = 0
in Ω × (0, T )
(6.1b)
p(x, t) = pp (x, t)

on Γp × (0, T )
v

b = vn (x, t)
v(x, t) · n

on Γ × (0, T )
31

(6.1c)
(6.1d)

100
SUPG
VI - SS
VI - RS

90
80

Efficiency (%)

Wall-clock time (s)

102

101

70
60
50
40
30

SUPG
VI - SS
VI - RS

20
100

1

2

4

8

10

16

1

2

# of cores

4

8

16

# of cores

(a) SUPG - solve time

(b) SUPG - parallel efficiency
100
90
80

Efficiency (%)

Wall-clock time (s)

10

DG
VI - SS
VI - RS

2

101

70
60
50
40
30

10

DG
VI - SS
VI - RS

20

0

1

2

4

8

10

16

1

# of cores

2

4

8

16

# of cores

(c) DG - solve time

(d) DG - parallel efficiency

Figure 15. 3D advection-diffusion: weak-scaling plots with approximately 100k
degrees-of-freedom per core and the corresponding parallel efficiencies.
∂c(x, t)
+ v(x, t) · grad[c(x, t)] − div[D(x, t)grad[c(x, t)]] = f (x, t) in Ω × (0, T )
∂t
c(x, t) = cp (x, t)
on ΓD × (0, T )
p

b (x) · (v(x, t)c(x, t) − D(x, t)grad[c(x, t)]) = q (x, t)
n
−b
n(x) · D(x, t)grad[c(x, t)] = q p (x, t)

N

on Γinflow × (0, T )

(6.1e)
(6.1f)
(6.1g)

on ΓN
outflow × (0, T ) (6.1h)

c(x, 0) = c0 (x)

in Ω

(6.1i)

where equations (6.1a) through (6.1d) represent the Darcy equation, and equations (6.1e) through
(6.1i) represent the transient advection-diffusion equation. To complete the coupled problem, the
viscosity is assumed to depend exponentially on concentration:


µ(c(x, t)) = µ0 exp Rc c(x, t)
(6.2a)
32

100
90

Efficiency (%)

Wall-clock time (s)

102

101

SUPG
VI - SS
VI - RS

100
1

2

80
70
60
SUPG
VI - SS
VI - RS

50

4

8

40

16

1

2

# of cores

4

8

16

# of cores

(a) SUPG - solve time

(b) SUPG - parallel efficiency
100
90

Efficiency (%)

Wall-clock time (s)

102

101

DG
VI - SS
VI - RS

100

1

2

80
70
60
DG
VI - SS
VI - RS

50

4

8

40

16

1

2

# of cores

4

8

16

# of cores

(c) DG - solve time

(d) DG - parallel efficiency

Figure 16. 3D advection-diffusion: strong-scaling plots for approximately 500k
degrees of freedom (h-size = 1/80 and 1/40 for GAL and DG respectively) and the
corresponding parallel efficiencies.



µ(c(x, t)) = µ0 exp Rc (1 − c(x, t))

(6.2b)

where µ0 is the base viscosity of the less viscous fluid and Rc is the log-mobility ratio in an isothermal
miscible displacement.
To solve the transient advection-diffusion equation, we employ the method of horizontal lines
[Grossmann et al., 2007], which first discretizes the time derivatives, thereby giving rise to timeindependent equations. The time interval of interest is divided into N sub-intervals. That is,
[0, T ] :=

N
[

[tn , tn+1 ]

n=0
33

(6.3)

pp (x, t) = 202650 Pa
cp (x, t) = 0.0

c0 (x) = 1.0

pp (x, t) = 101325 Pa
q p (x, t) = 0.0

25 m

q p (x, t) = vn (x, t) = 0.0

q p (x, t) = vn (x, t) = 0.0
50 m

(b) Log scale permeability field (m2 )

(a) Problem description

Figure 17. 2D miscible displacement: Pictorial description of the boundary value
problems for the coupled Darcy and advection-diffusion equations and the corresponding random permeability.
where tn denotes the n-th time-level. We assume that the time-step is uniform, which can be
written as:
∆t = tn+1 − tn

(6.4)

One can then employ the finite-dimensional VI solvers for these resulting equations, which were
described earlier in this paper. This implies that we will still be solving elliptic VIs of first kind
but at each time level. This procedure will be illustrated below using the backward Euler method.
However, a detailed discussion on the effect of time-stepping schemes in meeting maximum principles can be found in [Nakshatrala et al., 2016]. For a transient analysis, the proposed framework
outlined in Section 4.2.2 is modified as follows:
Step 1: Set t = 0.0, n = 0, and c(n) = c0 .
Step 2: Solve Darcy equation:
(a) Compute µ(c(n) ).
Table 9. 2D miscible displacement: problem parameters
Parameter
b(x)
µ(c(x, t))
µ0
Rc
k(x)
f (x, t)
ρ
αL
αT
αD
Number of elements
Darcy degrees-of-freedom
Advection-diffusion degrees-of-freedom
34

Value
{0, 0} m/s2
(6.2a)
3.95 · 10−5 Pa s
3
varies
0
479 kg/m3
10−1 m
10−5 m
10−9 m2 /s
31,250
94,125
125,000

(a) cDG at t = 0.5 years

(b) cDG at t = 1.0 years

(c) cRS at t = 0.5 years

(d) cRS at t = 1.0 years

Figure 18. 2D miscible displacement: Concentration fields under the DG (cDG )
formulation and VI - RS (cRS ) method at various time levels. White regions denote
violations of the maximum principle and the non-negative constraint.
(b) Assemble K vv , K vp , K pv , K pp , f v and f p .
(c) Solve for v (n) .
Step 3: Solve advection-diffusion equation:
(a) Compute D(n) .
(b) Assemble K c and f c using c(n) .
(n+1)
(c) Solve for cDG .
(n+1)
(n+1)
(d) Clip cDG and obtain cCLIP .
(n+1)
(n+1)
(e) Solve the bounded constraint problem for cRS with cCLIP as the initial guess.
(n+1)
Step 4: Set c(n+1) ←− cRS , t ←− t + ∆t, and n ←− n + 1.
Step 5: If n < N go to Step2.
where K vv , K vp , K pv , K pp , f v , and f p are the assembled matrices and vectors for the Darcy
equation, and K c and f c are for the transient advection-diffusion equation. The finite element
discretization and solution strategy for the steady-state Darcy equations can be found in Appendix
B.
6.2. Numerical results. Consider a 50m×25m rectangular domain with heterogeneous permeability, as shown in Figure 17. The flow will be modeled using Darcy equations, in which the
viscosity depends on the concentration of the attendant chemical species, and the transport of the
chemical species will be modeled using advection-diffusion equations. For the flow subproblem, we
prescribe the pressure boundary conditions on the left and right sides of the domain and no flow
35

(a) cCLIP - cRS at t = 0.5 years

(b) cCLIP - cRS at t = 1.0 years

(c) cDG - cRS at t = 0.5 years

(d) cDG - cRS at t = 1.0 years

Figure 19. 2D miscible displacement: Differences between concentration fields
under the DG formulation (cDG ), VI - RS method (cRS ), and clipping procedure
(cCLIP ).
boundary conditions on the top and bottom. For the transport subproblem, an initial concentration
of unity is prescribed everywhere in the domain, and a Dirichlet boundary condition of zero along
the left side and zero flux boundary conditions on the remaining sides are prescribed. A time-step
∆t = 1 day is used to simulate the miscible displacement over a time interval T = 1 year. All
other problem parameters and material properties can be found in Table 9. Figure 18 depicts the
concentration profiles under the DG and VI - RS methods at t = 0.5 and t = 1.0 years. It can be
seen that violations in the maximum principles occur even under the coupled flow and transport
computational framework. Furthermore, the violations do not go away as the simulation progresses
in time. As it may be difficult to distinguish between the VI - RS and DG or clipping procedures
by directly plotting the solutions, we show the differences in the solutions in Figure 19. It can be
seen that there are significant discrepancies in the development of the plumes.
To understand the performance of our VI-based solvers for large-scale versions of this problem,
we now consider a 50m×25m×5m box domain with heterogeneous permeability as shown in Figure
20. Table 10 lists all the necessary problem parameters, and the same time-step and time interval
from the previous problem is used. This problem is now solved in parallel across 40 MPI processes,
and Figure 21 depicts the numerical results under the DG formulation without VI - RS and DG
formulation with VI - RS. The exact regions where violations in the maximum principle and the
non-negative constraint occur are shown in Figures 21c and 21d. First, we note that this proposed
framework can successfully eliminate the violations that occur in a large-scale miscible displacement
simulation. We also note that the development and displacement of the plumes is significant affected
36

vn (x, t) = 0.0
q P (x, t) = 0.0

vn (x, t) = 0.0
q P (x, t) = 0.0

pP (x, t) = 202,650 Pa
cP (x, t) = 1.0

c0 (x) = 0.0

pP (x, t) = 101,325 Pa
q P (x, t) = 0.0
vn (x, t) = 0.0
q P (x, t) = 0.0
vn (x, t) = 0.0
q P (x, t) = 0.0

(b) Log scale permeability field (m2 )

(a) Problem description

Figure 20. 3D miscible displacement: Pictorial description of the boundary value
problems (50m×25m×5m domain with 781,250 structured hexahedrons) for the coupled Darcy and advection-diffusion equations and the corresponding random permeability.
by whether VI - RS is applied or not; the differences between Figures 21b and 21f are quite evident
unlike the 2D example. Next, we note from Figure 22 that enforcing the bounded constraints under
the VI - RS method will not drastically increase the total number of KSP iterations needed for
either the Darcy or advection-diffusion equations. However, the wall-clock time shown in Figure 23
indicates that the VI - RS method is very expensive. The oscillatory behavior of both the solver
iterations and wall-clock time at each time level is largely attributed to the heterogeneous nature of
the problem as well as the number of maximum principle violating degrees-of-freedom that naturally
Table 10. 3D miscible displacement: problem parameters
Parameter
b(x)
µ(c(x, t))
µ0
Rc
k(x)
f (x, t)
ρ
αL
αT
αD
Number of elements
Darcy degrees-of-freedom
Advection-diffusion degrees-of-freedom
37

Value
{0, 0, −9.81} m/s2
(6.2b)
3.95 · 10−5 Pa s
3
varies
0
479 kg/m3
10−1 m
10−5 m
10−9 m2 /s
781,250
3,165,625
6,250,000

(a) cDG at t = 0.4 years

(b) cDG at t = 1.0 years

(c) cDG at t = 0.4 years

(d) cDG at t = 1.0 years

(e) cRS at t = 0.4 years

(f) cRS at t = 1.0 years

Figure 21. 3D miscible displacement: Top (DG) and bottom (VI) show regions
with concentrations above 0.5. Middle figures show regions with concentrations
above 1.0 (green) and below 0.0 (purple) (see online version for color figures.)
38

35
VI - RS
KSP Darcy (no VI)
KSP Darcy (with VI)
KSP advection-diffusion (no VI)
KSP advection-diffusion (with VI)

30

Solver iterations

25

20

15

10

5

0

0

0.1

0.2

0.3

0.4
0.5
0.6
Time level (years)

0.7

0.8

0.9

1

Figure 22. 3D miscible displacement: Total number of KSP iterations at each
time level for the Darcy and advection-diffusion equation with and without VI. Also
shown is the total number of VI iterations at each time level.

400
No VI
With VI

350

Wall-clock time (s)

300
250

Total wall-clock time = 10.8 hours

200
150
100
50
Total wall-clock time = 33.3 minutes
0

0

0.1

0.2

0.3

0.4
0.5
0.6
Time level (years)

0.7

0.8

0.9

1

Figure 23. 3D miscible displacement: The wall-clock time at each time level with
and without VI across 40 cores. Also shown is the sum of the wall-clock time across
all time levels.
arise out of the DG formulation. Although VI-based solvers like VI - RS can enforce maximum
principles and the non-negative constraint, we have observed that applying such methodologies can
39

make the overall advection-diffusion finite element simulation up to 20 times as expensive even in
a parallel environment.
Before we draw any further conclusions of this paper, we acknowledge that we did not perform a numerical h-convergence study. This is due to, despite our best efforts, failure to find
an advection-diffusion boundary value problem that considers anisotropy, has an analytical solution, and violates discrete maximum principles. We, therefore, illustrated the performance of the
proposed computational framework through other means, as presented in the previous sections.

7. CONCLUDING REMARKS
We presented a robust computational framework based on VIs for diffusion and advectiondiffusion equations that satisfies the discrete maximum principles and the non-negative constraint.
The framework is applicable to large-scale and transient problems, and can be solved in a parallel
setting. The main contributions of this paper and the salient features of the proposed formulation
can be summarized as follows.
(C1) Realizing and posing the advection-diffusion problem as a variational inequality (VI) to
meet the discrete maximum principles and the non-negative constraint.
(C2) For large-scale problems, we have demonstrated that QP solvers, which is a special case
of VIs, are just as good as VI solvers for symmetric and positive-definite problems like the
diffusion equation. On the other hand, the proposed VI-based framework can also handle
non-self-adjoint operators.
(C3) Unlike the non-negative framework proposed in [Mudunuru and Nakshatrala, 2016a], which
is based on a mixed least-squares WF, the proposed framework can utilize any finite
element formulation including single-field formulations, and these formulations need not
result in symmetric and positive definite coefficient matrices.
(C4) The proposed framework allows one to leverage on existing state-of-the-art computational
frameworks for solving VIs. In particular, the Firedrake project, which provides access to
parallel solvers in PETSc and TAO libraries, can serve as a suitable platform for implementing the proposed framework, as illustrated in this paper.
(C5) This framework is suitable for many important applications like miscible displacement,
subsurface remediation, and transport of radionuclides. In these applications, one encounters not only highly anisotropic medium properties but also highly non-linear phenomena
due to aqueous complexation and kinetic reactions.
A logical next-step of our work is to extend this computational framework to advective-diffusivereactive equations.

ACKNOWLEDGMENTS
The authors acknowledge the financial support from the Department of Energy Office of Science
Graduate Student Research (SCGSR) award. The authors also acknowledge the use of the Opuntia
Cluster from the Center of Advanced Computing and Data Systems (CACDS) at the University of
Houston to carry out the research presented here. The opinions expressed in this paper are those
of the authors and do not necessarily reflect that of the sponsors.
40

Appendix A. Firedrake Project
The Firedrake project [Rathgeber et al., 2015; Luporini et al., 2016, 2015] is a python-based
library that provides an automated system for the solution of partial differential equations using
the finite element method. Like the FEniCS Project [Logg et al., 2012; Alnæs et al., 2015], it is
also built upon several scientific packages and can employ parallel computing tools across either
CPUs or GPUs to obtain the solution. Two of its main leveraged components are the Unified Form
Language (UFL) [Alnæs et al., 2014], used to declare finite element discretizations of variational
forms, and the PyOP2 system [Rathgeber et al., 2012; Markall et al., 2013], used for the parallel
assembly of the finite element discrete formulations. The main difference between the FEniCS and
Firedrake project is that all data structures, linear solvers, non-linear solvers, and optimization
solvers for the latter are provided entirely by the PETSc and TAO libraries. The mesh can either
be generated internally or imported from third party mesh generators like GMSH [Geuzaine and
Remacle, 2009], and the parallel partitioning of the mesh is achieved through packages like Chaco
[Hendrickson and Leland, 1995]. Another important feature utilized in this paper is extruded
meshes. The internal mesh algorithm generates and partitions a 2D quadrilateral base mesh and is
extruded into a hexahedron mesh using the algorithms listed in [Homolya and Ham, 2016; McRae
et al., 2014]. To facilitate the readers to be able to reproduce the results presented in this paper,
we provided some useful Firedrake-related files below.

Listing 1. 2D GAL diffusion example
1
2

# Load firedrake environment
from firedrake import *

3
4
5

# Create mesh
mesh = UnitSquareMesh(200,200,quadrilateral=True)

6
7
8
9
10
11

#
P
Q
u
v

Function spaces
= FunctionSpace(mesh, 'Lagrange', 1)
= TensorFunctionSpace(mesh, 'Lagrange', 1)
= TrialFunction(P)
= TestFunction(P)

12
13
14
15
16
17
18
19

# Bounds
cmin = 0.0
cmax = PETSc.INFINITY
lb = Function(P)
lb.assign(cmin)
ub = Function(P)
ub.assign(cmax)

20
21
22
23
24

# Diffusion tensor
eps = 1e−4
D = interpolate(Expression(('eps*x[0]*x[0]+x[1]*x[1]','−(1−eps)*x[0]*x[1]',
'−(1−eps)*x[0]*x[1]','eps*x[1]*x[1]+x[0]*x[0]'),eps=eps),Q)

25
26
27

# Forcing function
lb = 0.375
41

28
29
30

ub = 0.625
f = interpolate(Expression(('x[0] ≥ lb && x[0]≤ ub && x[1]
x[1] ≤ ub ? 1.0 : 0.0'),lb=lb,ub=ub),P)

≥

lb &&

31
32
33
34

# GAL formulation
a = dot(D * grad(u), grad(v)) * dx
L = v * f * dx

35
36
37

# Homogeneous boundary conditions
bcs = DirichletBC(P, Constant(0.0), (1,2,3,4))

38
39
40

# Assemble coefficient matrix
A = assemble(a, bcs=bcs)

41
42
43
44
45
46
47
48
49

# Assemble forcing vector
tmp = Function(P)
bcs.apply(tmp)
b = Function(P)
bfree = assemble(L)
rhs bcs = assemble(action(a,tmp))
b.assign(bfree − rhs bcs)
bcs.apply(b)

50
51
52
53
54
55
56

# Create PETSc solver
initial solver = PETSc.KSP().create(PETSc.COMM WORLD)
initial solver.setOptionsPrefix("initial ")
initial solver.setOperators(A.M.handle)
initial solver.setFromOptions()
initial solver.setUp()

57
58
59
60
61

# Solve problem
solution = Function(P)
with b.dat.vec ro as b vec, solution.dat.vec as sol vec:
initial solver.solve(b vec,sol vec)

Listing 2. 2D SUPG advection-diffusion example
1
2

# Load firedrake environment
from firedrake import *

3
4
5

# Load GMSH file
mesh = Mesh('square hole.msh')

6
7
8
9
10
11

#
P
V
u
v

Function spaces
= FunctionSpace(mesh, 'Lagrange', 1)
= VectorFunctionSpace(mesh, 'Lagrange', 1)
= TrialFunction(P)
= TestFunction(P)

12
13
14

# Bounds
cmin = 0.0
42

15
16
17
18
19

cmax = 1.0
lb = Function(P)
lb.assign(cmin)
ub = Function(P)
ub.assign(cmax)

20
21
22
23

# Velocity field
velocity = interpolate(Expression(('cos(2*pi*x[1]*x[1])',
'sin(2*pi*x[0])+cos(2*pi*x[0]*x[0])')),V)

24
25
26
27
28
29
30
31
32

# Diffusion tensor
alphaT = Constant(1e−5)
alphaL = Constant(1e−1)
alphaD = Constant(1e−9)
normv = sqrt(dot(velocity,velocity))
Id = Identity(mesh.geometric dimension())
D = (alphaD + alphaT*normv)*Id +
(alphaL − alphaT)*outer(velocity,velocity)/normv

33
34
35

# Forcing function
f = Constant(0.0)

36
37
38
39
40
41
42
43

# SUPG weak form
h = CellSize(mesh)
Pe = h/(2*normv)*dot(velocity,grad(v))
ar = Pe*(dot(velocity,grad(u)) − div(D*grad(u)))*dx
a = ar + v*dot(velocity,grad(u))*dx + dot(grad(v),D*grad(u))*dx
a = dot(D*grad(u), grad(v))*dx
L = (v + Pe)*f*dx

44
45
46
47
48

# Boundary conditions
bc1 = DirichletBC(P, Constant(0.0), (12,13,14,15)) # Outer square
bc2 = DirichletBC(P, Constant(1.0), (16,17,18,19)) # Inner square
bcs = [bc1,bcs2]

49
50
51

# Homogeneous boundary conditions
bcs = DirichletBC(P, Constant(0.0), (1,2,3,4))

52
53
54

# Assemble coefficient matrix
A = assemble(a, bcs=bcs)

55
56
57
58
59
60
61
62
63
64
65

# Assemble forcing vector
tmp = Function(P)
for bc in bcs:
bc.apply(tmp)
b = Function(P)
bfree = assemble(L)
rhs bcs = assemble(action(a,tmp))
b.assign(bfree − rhs bcs)
for bc in bcs:
bc.apply(b)

66

43

67
68
69
70
71
72

# Create PETSc solver
initial solver = PETSc.KSP().create(PETSc.COMM WORLD)
initial solver.setOptionsPrefix("initial ")
initial solver.setOperators(A.M.handle)
initial solver.setFromOptions()
initial solver.setUp()

73
74
75
76
77

# Solve problem
solution = Function(P)
with b.dat.vec ro as b vec, solution.dat.vec as sol vec:
initial solver.solve(b vec,sol vec)

Listing 3. 3D DG advection-diffusion example
1
2

# Load firedrake environment
from firedrake import *

3
4
5

# Number of elements in each spatial dimension
seed = 40

6
7
8
9
10

# 2D base mesh
mesh = UnitSquareMesh(seed,seed,quadrilateral=True)
# Extruded mesh
mesh = ExtrudedMesh(meshbase,seed)

11
12
13
14
15
16

#
P
V
u
v

Function spaces
= FunctionSpace(mesh, 'DG', 1)
= VectorFunctionSpace(mesh, 'DG', 1)
= TrialFunction(P)
= TestFunction(P)

17
18
19
20
21
22
23
24

# Bounds
cmin = 0.0
cmax = PETSc.INFINITY
lb = Function(P)
lb.assign(cmin)
ub = Function(P)
ub.assign(cmax)

25
26
27
28
29

# Velocity field
velocity = interpolate(Expression(('0.3*sin(2*pi*x[2])+cos(3*pi*x[1])',
'0.65*sin(2*pi*x[0])+0.3*cos(5*pi*x[2])',
'sin(4*pi*x[1])+0.65*cos(6*pi*x[0])')),V)

30
31
32
33
34
35
36
37

# Diffusion tensor
alphaT = Constant(1e−5)
alphaL = Constant(1e−1)
alphaD = Constant(1e−9)
normv = sqrt(dot(velocity,velocity))
Id = Identity(mesh.geometric dimension())
D = (alphaD + alphaT*normv)*Id +
44

38

(alphaL − alphaT)*outer(velocity,velocity)/normv

39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

# Forcing function
f1 = interpolate(Expression(('x[0] ≥
x[1] ≤ 0.3 && x[2] ≥ 0.1 && x[2]
f2 = interpolate(Expression(('x[0] ≥
x[1] ≤ 0.5 && x[2] ≥ 0.2 && x[2]
f3 = interpolate(Expression(('x[0] ≥
x[1] ≤ 0.8 && x[2] ≥ 0.3 && x[2]
f4 = interpolate(Expression(('x[0] ≥
x[1] ≤ 0.6 && x[2] ≥ 0.2 && x[2]
f5 = interpolate(Expression(('x[0] ≥
x[1] ≤ 0.3 && x[2] ≥ 0.6 && x[2]
f6 = interpolate(Expression(('x[0] ≥
x[1] ≤ 0.6 && x[2] ≥ 0.7 && x[2]
f7 = interpolate(Expression(('x[0] ≥
x[1] ≤ 0.8 && x[2] ≥ 0.8 && x[2]
f8 = interpolate(Expression(('x[0] ≥
x[1] ≤ 0.5 && x[2] ≥ 0.7 && x[2]
f = f1 + f2 + f3 + f4 + f5 + f6 + f7

0.4 && x[0]
0.2 ? 1.0
0.8 && x[0]
≤ 0.3 ? 1.0
0.5 && x[0]
≤ 0.4 ? 1.0
0.3 && x[0]
≤ 0.3 ? 1.0
0.5 && x[0]
≤ 0.7 ? 1.0
0.6 && x[0]
≤ 0.8 ? 1.0
0.4 && x[0]
≤ 0.9 ? 1.0
0.1 && x[0]
≤ 0.8 ? 1.0
+ f8
≤

0.5 && x[1]
: 0.0')),P)
≤ 0.9 && x[1]
: 0.0')),P)
≤ 0.6 && x[1]
: 0.0')),P)
≤ 0.4 && x[1]
: 0.0')),P)
≤ 0.6 && x[1]
: 0.0')),P)
≤ 0.7 && x[1]
: 0.0')),P)
≤ 0.5 && x[1]
: 0.0')),P)
≤ 0.2 && x[1]
: 0.0')),P)
≤

≥

0.2 &&

≥

0.4 &&

≥

0.7 &&

≥

0.5 &&

≥

0.2 &&

≥

0.5 &&

≥

0.7 &&

≥

0.4 &&

58
59
60
61
62
63

# Parameters
h = Constant(1/float(seed)) # h−size
gamma = Constant(8/3) # Penalty term
n = FacetNormal(mesh) # Unit outward normal
vn = 0.5*(dot(velocity,n) + abs(dot(velocity,n))) # Upwinding term

64
65
66
67
68
69
70
71
72
73

# DG weak formulation
a = inner(D * grad(u), grad(v)) * dx(degree=(3,3)) \
− dot(jump(v,n),avg(D*grad(u)))*(dS h + dS v) \
− dot(avg(D*grad(v)),jump(u,n))*(dS h + dS v) \
+ gamma/h*dot(jump(v,n),jump(u,n))*(dS h + dS v) \
− dot(grad(v),velocity*u)*dx(degree=(3,3)) \
+ dot(jump(v),vn('+')*u('+')−vn('−')*u('−'))*(dS h+dS v) \
+ dot(v, vn*u)*(ds v+ds t+ds b)
L = v * f * dx(degree=(3,3))

74
75
76
77
78
79

# Boundary conditions
bc1 = DirichletBC(V, Constant(0.0), (1,2,3,4), method="geometric")
bc2 = DirichletBC(V, Constant(0.0), "bottom", method="geometric")
bc3 = DirichletBC(V, Constant(0.0), "top", method="geometric")
bcs = [bc1, bc2, bc3]

80
81
82

# Homogeneous boundary conditions
bcs = DirichletBC(P, Constant(0.0), (1,2,3,4))

83
84
85

# Assemble coefficient matrix
A = assemble(a, bcs=bcs)

86
87
88
89

# Assemble forcing vector
tmp = Function(P)
for bc in bcs:
45

90
91
92
93
94
95
96

bc.apply(tmp)
b = Function(P)
bfree = assemble(L)
rhs bcs = assemble(action(a,tmp))
b.assign(bfree − rhs bcs)
for bc in bcs:
bc.apply(b)

97
98
99
100
101
102
103

# Create PETSc solver
initial solver = PETSc.KSP().create(PETSc.COMM WORLD)
initial solver.setOptionsPrefix("initial ")
initial solver.setOperators(A.M.handle)
initial solver.setFromOptions()
initial solver.setUp()

104
105
106
107
108

# Solve problem
solution = Function(P)
with b.dat.vec ro as b vec, solution.dat.vec as sol vec:
initial solver.solve(b vec,sol vec)

Listing 4. GMSH geometry file for Listing 2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

Point(1) = {0, 0, 0, 1.0};
Point(2) = {1, 0, 0, 1.0};
Point(3) = {1, 1, 0, 1.0};
Point(4) = {0, 1, 0, 1.0};
Point(5) = {4/9, 4/9, 0, 1.0};
Point(6) = {5/9, 4/9, 0, 1.0};
Point(7) = {5/9, 5/9, 0, 1.0};
Point(8) = {4/9, 5/9, 0, 1.0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line Loop(9) = {4, 1, 2, 3};
Line Loop(10) = {8, 5, 6, 7};
Plane Surface(11) = {9, 10};
Physical Line(12) = {4};
Physical Line(13) = {1};
Physical Line(14) = {2};
Physical Line(15) = {3};
Physical Line(16) = {7};
Physical Line(17) = {6};
Physical Line(18) = {5};
Physical Line(19) = {8};
Physical Surface(20) = {11};
46

Listing 5. Semi-smooth (VI - SS) method
1
2
3

# Create TAO object
ss solver = PETSc.TAO().create(PETSc.COMM WORLD)
ss solver.setOptionsPrefix("ss ")

4
5
6
7
8
9
10
11
12
13

# Semi−smooth call−backs
def ss formJac(tao, petsc x, petsc J, petsc JP, A=None, a=None, bcs=None):
A = assemble(a, bcs=bcs, tensor=A)
A.M. force evaluation()
def ss formFunc(tao, petsc x, petsc g, A=None, a=None, b=None, bcs=None):
A = assemble(a, bcs=bcs, tensor=A)
with b.dat.vec as b vec:
A.M.handle.mult(petsc x, petsc g)
petsc g.axpy(−1.0,b vec)

14
15
16
17
18
19
20
21
22

# Setup
ss con = Function(solution.function space())
with ss con.dat.vec as con vec, lb.dat.vec as lb vec, ub.dat.vec as ub vec:
ss solver.setConstraints(ss formFunc, con vec, kargs={'A':A,'a':a,'b':b,'bcs':bcs})
ss solver.setJacobian(ss formJac,A.M.handle, kargs={'A':A,'a':a,'bcs':bcs})
ss solver.setType(PETSc.TAO.Type.SSFLS) # can also be ASFLS/SSILS/ASILS
ss solver.setVariableBounds(lb vec,ub vec)
ss solver.setFromOptions()

23
24
25
26
27

# Solve problem
def viss():
with solution.dat.vec as sol vec:
ss solver.solve(sol vec)

Listing 6. Reduced-space active set (VI - RS) method
1
2
3

# Create SNES object
rs solver = PETSc.SNES().create(PETSc.COMM WORLD)
rs solver.setOptionsPrefix("rs ")

4
5
6
7
8
9
10
11
12

# Reduced−space active−set call−backs
def rs formJac(snes, petsc x, petsc J, petsc JP):
pass
def rs formFunc(snes, petsc x, petsc g, A=None, b=None):
with b.dat.vec as b vec:
A.M.handle.mult(petsc x, petsc g)
petsc g.axpy(−1.0,b vec)
rs con = Function(solution.function space())

13
14
15
16
17
18
19
20

# Solve problem
def virs():
with solution.dat.vec as sol vec, lb.dat.vec as lb vec, ub.dat.vec as ub vec:
with rs con.dat.vec as con vec:
rs solver.setFunction(rs formFunc, con vec, kargs={'A':A,'b':b})
rs solver.setJacobian(rs formJac,A.M.handle)
rs solver.setType(PETSc.SNES.Type.VINEWTONRSLS)
47

21
22
23
24

rs
rs
rs
rs

solver.setVariableBounds(lb vec,ub vec)
solver.setFromOptions()
solver.solve(None,sol vec)
solver.reset()

Listing 7. Trust region Newton (QP - TRON) method
1
2
3

# Create TAO object
tron solver = PETSc.TAO().create(PETSc.COMM WORLD)
tron solver.setOptionsPrefix("tron ")

4
5
6
7
8
9
10
11
12
13
14

# TRON call−backs
def tron formHess(tao, petsc x, petsc H, petsc HP):
pass
def tron formObjGrad(tao, petsc x, petsc g, A=None, b=None):
with b.dat.vec ro as b vec:
A.M.handle.mult(petsc x, petsc g)
xtHx = petsc x.dot(petsc g)
xtf = petsc x.dot(b vec)
petsc g.axpy(−1.0,b vec)
return 0.5*xtHx − xtf

15
16
17
18
19
20
21
22

# Setup
with lb.dat.vec as lb vec, ub.dat.vec as ub vec:
tron solver.setVariableBounds(lb vec,ub vec)
tron solver.setObjectiveGradient(tron formObjGrad, kargs={'A':A,'b':b})
tron solver.setHessian(tron formHess,A.M.handle)
tron solver.setType(PETSc.TAO.Type.TRON)
tron solver.setFromOptions()

23
24
25
26
27

# Solve problem
def qptron():
with solution.dat.vec as sol vec:
tron solver.solve(sol vec)

Appendix B. Solution strategy for the Darcy equation
B.1. Weak formulation. Let w(x) and q(x) represent the weighting functions for velocity
and pressure respectively. The relevant function spaces read as follows:
n
o
b = vn on Γv
V := v ∈ (L2 (Ω))d div[v] ∈ L2 (Ω), v · n
(B.1a)
n
o
b = 0 on Γv
W := w ∈ (L2 (Ω))d div[w] ∈ L2 (Ω), w · n
(B.1b)
P := L2 (Ω)

(B.1c)

where L2 (Ω) is the space of square integrable functions. The WF under the classical mixed formulation for the Darcy equations (6.1a) through (6.1d) reads: Find v(x) ∈ V and p(x) ∈ P such that
we have:
B(w, q; v, p) = L(w, q)

∀w(x) ∈ W, q(x) ∈ P
48

(B.2)

where the bilinear form and linear functional are:


µ(c(x))
B(w, q; v, p) := w(x);
v(x)
− (div[w(x)]; p(x))Ω − (q(x); div[w(x)])Ω
k(x)
Ω

b (x); pP Γp
L(w, q) := (w(x); ρb(x))Ω − w(x) · n

(B.3)
(B.4)

The lowest order Raviart-Thomas (RT0) space [Raviart and Thomas, 1977] is employed because it
ensures element-wise mass conservation. To map the RT0 element onto quadrilateral and extruded
hexahedrons, contravariant Piola mapping is used (see [Rognes et al., 2009; Bercea et al., 2016] for
further details). The discrete formulations may be assembled into the following block format:
   

v
fv
K vv K vp
=
(B.5)
p
fp
K pv K pp
where the terms in equation (B.3) respectively correspond to K vv , K vp , and K pv , and equation
(B.4) corresponds to f v . It should be noted that K pp and f p are a zero matrix and zero vector,
respectively.
B.2. Preconditioning methodology. Equation (B.5) is a saddle-point system which is tricky
to precondition effectively for large-scale problems. Several classes of iterative solvers and preconditioning strategies exist for these types of problems [Benzi et al., 2005; Elman et al., 2006; Murphy
et al., 2000]. One could alternatively employ hybridization techniques [Cockburn et al., 2009] which
introduces Lagrange multipliers to reduce the difficulty of solving such problems. In this paper, we
employed a Schur complement approach to precondition the saddle-point system. Conceptually,
the problem at hand is a 2×2 block matrix:


K vv K vp
K=
(B.6)
K pv
0
which admits a full factorization of




I
0
K vv 0
I K −1
vv K vp
K=
K pv K −1
I
0
S
0
I
vv

(B.7)

where I is the identity matrix and
S = −K pv K −1
vv K vp
is the Schur complement. The inverse can therefore be written as:

  −1


I −K −1
K vv
0
I
0
−1
vv K vp
K =
0
I
0
S −1
−K pv K −1
I
vv

(B.8)

(B.9)

−1
The task at hand is to approximate K −1
vv and S . Since K vv is a mass matrix for the Darcy
equation, we can invert it using the ILU(0) (incomplete lower upper) solver. We note that the
Schur complement is spectrally a Laplacian, so we can employ a diagonal mass-lumping of K vv to
give a good approximation to K −1
vv . That is, we can use

S p = −K pv diag (K vv )−1 K vp

(B.10)

to precondition the inner solver inverting S. For this block we employ the multi-grid V-cycle on
S p using the Trilinos ML package ([Sala et al., 2004]). These blocks are symmetric and positivedefinite so one could employ the CG solvers to obtain the inverses. When the inverses are obtained,
only a single sweep of flexible GMRES is needed to obtain the full solution. However, instead of
individually solving for K −1
vv and S p , we could alternatively apply a single sweep of ILU(0) and
49

multi-grid, respectively, and rely on GMRES to solve the entire block system. By providing less
accurate approximations of the inner individual blocks, the number of GMRES iterations for the
overall system increases but the numerical accuracy remains the same. We have found that this
methodology is computationally less expensive and more practical for large-scale computations.
One could alternatively employ one of the factorizations (either lower or upper) to decrease the
computational cost associated with setting up the preconditioner. Below are some necessary PETSc
command-line options for the described Schur complement approach.

Listing 8. PETSc solver options for the Schur complement approach
1
2
3

parameters = {
# Outer solver
"ksp type": "gmres",

4

# Schur complement with full factorization
"pc type": "fieldsplit",
"pc fieldsplit type": "schur",
"pc fieldsplit schur fact type": "full",

5
6
7
8
9

# Diagonal mass lumping
"pc fieldsplit schur precondition": "selfp",

10
11
12

# Single sweep of ILU(0) for the mass matrix
"fieldsplit 0 ksp type": "preonly",
"fieldsplit 0 pc type": "ilu",

13
14
15
16

# Single sweep of multi−grid for the Schur complement
"fieldsplit 1 ksp type": "preonly",
"fieldsplit 1 pc type": "ml"

17
18
19
20

}

References
M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E.
Rognes, and G. N. Wells. The fenics project version 1.5. Archive of Numerical Software, 3(100):
9–23, 2015.
M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells. Unified form language: A
domain-specific language for weak formulations of partial differential equations. ACM Transactions on Mathematical Software, 40(2):9, 2014.
L. Armijo. Minimization of functions having Lipschitz-continuous first partial derivatives. Pacific
Journal of Mathematics, 16:1–3, 1966.
D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified analysis of discontinuous Galerkin
methods for elliptic problems. SIAM Journal on Numerical Analysis, 39:1749–1779, 2002.
S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith. Efficient management of parallelism in
object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen,
editors, Modern Software Tools in Scientific Computing, pages 163–202. Birkhäuser Press, 1997.
S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout,
W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini,
50

H. Zhang, and H. Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 3.7,
Argonne National Laboratory, 2016.
S. Benson and T. S. Munson. Flexible complementarity solvers for large-scale applications. Pacific
Journal of Mathematics, 21:155–168, 2006.
M. Benzi, G. H. Golub, and J. Liesen. Numerical solution of saddle point problems. Acta Numerica,
14:1–137, 2005.
G. T. Bercea, A. T. T. McRae, D. A. Ham, L. Mitchell, F. Rathgeber, L. Nardi, F. Luporini, and
P. H. J. Kelly. A numbering algorithm for finite element on extruded meshes which avoids the
unstructured mesh penalty. Submitted to Geoscientific Model Development (GMD), 2016. URL
http://arxiv.org/abs/1604.05937.
S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, Cambridge,
UK, 2004.
S. C. Brenner and L. R. Scott. The Mathematical Theory of Finite Element Methods. SpringerVerlag, New York, USA, 1994.
H. Brezis. Functional Analysis, Sobolev Spaces and Partial Differential Equations. Springer Science
& Business Media, 2010.
A. N. Brooks and T. J. R. Hughes. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations.
Computer Methods in Applied Mechanics and Engineering, 32:199–259, 1982.
O. Burdakov, I. Kapyrin, and Y. Vassilevski. Monotonicity recovering and accuracy preserving
optimization methods for postprocessing finite element solutions. Journal of Computational
Physics, 231:3126–3142, 2012.
J. Chang, S. Karra, and K. B. Nakshatrala. Large-scale optimization-based non-negative computational framework for diffusion equations: Parallel implementation and performance studies.
Journal of Scientific Computing, DOI: 10.1007/s10915-016-0250-5, 2016a.
J. Chang, K. B. Nakshatrala, and J. N. Reddy. Modification to Darcy-Forchheimer model due to
pressure-dependent viscosity: consequences and numerical solutions. To appear in Journal of
Porous Media, 2016b.
C. Y. Chen and E. Meiburg. Miscible porous media displacements in the quarter five-spot configuration. part 1. the homogeneous case. Journal of Fluid Mechanics, 371:233–268, 1998a.
C. Y. Chen and E. Meiburg. Miscible porous media displacements in the quarter five-spot configuration. part 2. effect of heterogeneities. Journal of Fluid Mechanics, 371:269–299, 1998b.
M. Chipot. Elliptic Equations: An Introductory Course. Birkhäuser, Basel, Switzerland, 2009.
P. G. Ciarlet and P. A. Raviart. Maximum principle and uniform convergence for the finite element
method. Computer Methods in Applied Methods and Engineering, 2:17–31, 1973.
B. Cockburn. Discontinuous Galerkin methods. Zeitschrift für Angewandte Mathematik und
Mechanik, 83(11):731–754, 2003. doi: 10.1002/zamm.200310088.
B. Cockburn, J. Gopalakrishnan, and R. Lazarov. Unified hybridization of discontinuous Galerkin,
mixed, and continuous Galerkin methods for second order elliptic problems. SIAM Journal on
Numerical Analysis, 47(2):1319–1365, 2009.
L. B. da Veiga, K. Lipnikov, and G. Manzini. The mimetic finite difference method for elliptic
problems, volume 11. Springer, 2014.
L. D. Dalcin, R. R. Paz, P. A. Kler, and A. Cosimo. Parallel distributed computing using Python.
Advances in Water Resources, 34(9):1124–1139, 2011.
51

J. W. Demmel. Applied Numerical Linear Algebra. Society for Industrial and Applied Mathematics,
Philadelphia, PA, USA, 1997.
T. Dombre, U. Frisch, J. M. Greene, M. Hénon, A. Mehr, and A. M. Soward. Chaotic streamlines
in the ABC flows. Journal of Fluid Mechanics, 167:353–391.
G. Duvaut and J. L. Lions. Inequalities in Mechanics and Physics. Springer-Verlag, Berlin, Germany, 1976.
H. C. Elman, D. J. Silvester, and A. J. Wathen. Finite elements and fast iterative solvers. Journal
of Fluid Mechanics, 557(1):474–475, 2006.
US EPA. Cleaning up the nation’s waste sites: Markets and technology trends. Technical Report
EPA 542-R-04-015, 2004.
F. Facchinei and J.-S. Pang. Finite-Dimensional Variational Inequalities and Complementarity
Problems. Volume I. Springer-Verlag, New York, USA, 2003.
G. Fichera. Problemi elastostatici con vincoli unilaterali: il problema di Signorini con ambigue
condizioni al contorno. Accademia nazionale dei Lincei, 1964.
G. Fichera. Linear Elliptic Differential Systems and Eigenvalue Problems. Number 8. Springer,
1965.
A. Fischer. A special Newton-type optimization method. Optimization, 24:269–284, 1992.
A. Genty and C. Le Potier. Maximum and minimum principles for radionuclide transport calculations in geological radioactive waste repository: comparison between a mixed hybrid finite
element method and finite volume element discretizations. Transport in Porous Media, 88:65–86,
2011.
C. Geuzaine and J. F. Remacle. Gmsh: a three-dimensional finite element mesh generator with
built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, 79:1309–1331, 2009.
D. Gilbarg and N. S. Trudinger. Elliptic Partial Differential Equations of Second Order. Springer,
New York, USA, 2001.
R. Glowinski. Numerical Methods for Nonlinear Variational Problems. Springer-Verlag, Berlin,
Germany, 1984.
C. Grossmann, H. G. Roos, and M. Stynes. Numerical Treatment of Partial Differential Equations.
Springer, Berlin, Germany, 2007.
G. E. Hammond and P. C. Lichtner. Field-scale model for the natural attenuation of uranium at
the hanford 300 area using high-performance computing. Water Resources Research, 46:W09602,
2010.
W. Han and B. D. Reddy. Plasticity: Mathematical Theory and Numerical Analysis, volume 9.
Springer Science & Business Media, 2012.
I. Harari. Stability of semidiscrete formulations for parabolic problems at small time steps. Computer Methods in Applied Mechanics and Engineering, 193(15):1491–1516, 2004.
D. R. Harp and V. V. Vesselinov. Contaminant remediation decision analysis using information
gap theory. Stochastic Environmental Research and Risk Assessment, 27:159–168, 2013.
J. M. Heikoop, T. M. Johnson, K. H. Birdsell, P. Longmire, D. D. Hickmott, E. P. Jacobs, D. E.
Broxton, D. Katzman, V. V. Vesselinov, M. Ding, D. T. Vanimana, S. L. Reneaua, T. J. Goering,
J. Glessnerb, and A. Basu. Isotopic evidence for reduction of anthropogenic hexavalent chromium
in Los Alamos National Laboratory groundwater. Chemical Geology, 373:1–9, 2014.

52

B. Hendrickson and R. Leland. A multilevel algorithm for partitioning graphs. In Supercomputing
’95: Proceedings of the 1995 ACM/IEEE Conference on Supercomputing (CDROM), page 28,
New York, 1995. ACM Press. ISBN 0-89791-816-9. doi: http://doi.acm.org/10.1145/224170.
224228.
I. Hlavacek, J. Haslinger, J. Necas, and J. Lovisek. Solution of Variational Inequalities in Mechanics,
volume 66. Springer Science & Business Media, 2012.
M. Homolya and D. A. Ham. A parallel edge orientation algorithm for quadrilateral meshes. To
appear in SIAM Journal on Scientific Computing, 2016. URL http://arxiv.org/abs/1505.
03357.
W. Huang and Y. Wang. Discrete maximum principle for the weak Galerkin method for anisotropic
diffusion problems. Communications in Computational Physics, 18(1):65–90, 2015.
N. Kikuchi and J. T. Oden. Contact Problems in Elasticity: A Study of Variational Inequalities
and Finite Element Methods, volume 8. SIAM, 1988.
D. Kinderlehrer and G. Stampacchia. An Introduction to Variational Inequalities and Their Applications. SIAM Classics in Applied Mathematics, New York, USA, 2000.
HYPRE: High Performance Preconditioners. Lawrence Livermore National Laboratory. http:
//www.llnl.gov/CASC/hypre/.
D. R. Lester, G. Metcalfe, and M. G. Trefry. Is chaotic advection inherent to porous media flow?
Physical Review Letters, 111(17):174101, 2013.
J. Li and B. Riviére. High order discontinuous Galerkin method for simulating miscible flooding in
porous media. Computational Geoscience, 19:1251–1268, 2015a.
J. Li and B. Riviére. Numerical solutions of the incompressible miscible displacement equations in
heterogeneous media. Computer Methods in Applied Mechanics and Engineering, 292:107–121,
2015b.
J. Li and B. Riviére. Numerical modeling of miscible viscous fingering instabilities by high-order
methods. Transport in Porous Media, 113:607628, 2016.
C. J. Lin and J. Moré. Newton’s method for large bound-constrained optimization problems. SIAM
Journal on Optimization, 9:1100–1127, 1999.
J. L. Lions and G. Stampacchia. Variational inequalities. Communications on Pure and Applied
Mathematics, 20(3):493–519, 1967.
K. Lipnikov, M. Shashkov, D. Svyatskiy, and Y. Yassilevski. Monotone finite volume schemes for
diffusion equations on unstructured triangular and shape-regular polygonal meshes. Journal of
Computational Physics, 227:492–512, 2007.
R. Liska and M. Shashkov. Enforcing the discrete maximum principle for linear finite element
solutions for elliptic problems. Communications in Computational Physics, 3:852–877, 2008.
A. Logg, K. A. Mardal, and G. N. Wells. Automated solution of differential equations by the finite
element method: The FEniCS book, volume 84. Springer Science & Business Media, 2012.
T. De Luca, F. Facchinei, and C. Kanzow. A semismooth equation approach to the solution of
nonlinear complementarity problems. Mathematical Programming, 75:407–439, 1996.
F. Luporini, A. L. Varbanescu, F. Rathgeber, G. T. Bercea, J. Ramanujam, D. A. Ham, and P. H. J.
Kelly. Cross-loop optimization of arithmetic intensity for finite element local assembly. ACM
Transactions on Architecture and Code Optimization, 11(4):57:1–57:25, 2015.
F. Luporini, D. A. Ham, and P. H. J. Kelly. An algorithm for the optimization of finite element
integration loops. Submitted to ACM TOMS, 2016. URL http://arxiv.org/abs/1604.05872.
53

G. R. Markall, F. Rathgeber, L. Mitchell, N. Loriant, C. Bertolli, D. A. Ham, and P. H. J.
Kelly. Performance-portable finite element assembly using pyop2 and fenics. In Julian Martin
Kunkel, Thomas Ludwig, and Hans Werner Meuer, editors, 28th International Supercomputing
Conference, ISC, Proceedings, volume 7905 of Lecture Notes in Computer Science, pages 279–
289. Springer, 2013. doi: 10.1007/978-3-642-38750-0 21. URL http://dx.doi.org/10.1007/
978-3-642-38750-0_21.
A. T. T. McRae, G. T. Bercea, L. Mitchell, D. A. Ham, and C. J. Cotter. Automated generation
and symbolic manipulation of tensor product finite elements. Submitted to SIAM Journal on
Scientific Computing, 2014. URL http://arxiv.org/abs/1411.2940.
R. Mifflin. Semismooth and semiconvex functions in constrained optimization. SIAM Journal on
Control and Optimization, 15:957–972, 1977.
M. K. Mudunuru and K. B. Nakshatrala. On enforcing maximum principles and achieving elementwise species balance for advection-diffusion-reaction equations under the finite element method.
Journal of Computational Physics, 305:448–493, 2016a.
M. K. Mudunuru and K. B. Nakshatrala. On mesh restrictions to satisfy comparison principles,
maximum principles, and the non-negative constraint: Recent developments and new results.
Mechanics of Advanced Materials and Structures, DOI: 10.1080/15502287.2016.1166160, 2016b.
M. K. Mudunuru, M. Shabouei, and K. B. Nakshatrala. On local and global species conservation
errors for nonlinear ecological models and chemical reacting flows. In ASME 2015 International
Mechanical Engineering Congress and Exposition, pages V009T12A018–V009T12A018. American
Society of Mechanical Engineers, 2015.
T. Munson, J. Sarich, S. Wild, S. Benson, and L. C. McInnes. Toolkit for advanced optimization
(tao) users manual. Technical Report ANL/MCS-TM-322 - Revision 3.5, Argonne National
Laboratory, 2014. URL http://www.mcs.anl.gov/tao.
T. S. Munson, F. Facchinei, M. C. Ferris, A. Fischer, and C. Kanzow. The semismooth algorithm
for large scale complementarity problems. INFORMS Journal on Computing, 13:294–311, 2001.
M. F. Murphy, G. H. Golub, and A. J. Wathen. A note on preconditioning for indefinite linear
systems. SIAM Journal on Scientific Computing, 21(6):1969–1972, 2000.
H. Nagarajan and K. B. Nakshatrala. Enforcing the non-negativity constraint and maximum
principles for diffusion with decay on general computational grids. International Journal for
Numerical Methods in Fluids, 67:820–847, 2011.
A. Nagurney. Variational Inequalities. Technical report, Isenberg School of Management, University of Massachusetts Amherst, 2002. URL https://supernet.isenberg.umass.edu/Austria_
Lectures/fvisli.pdf.
A. Nagurney and D. Zhang. Projected Dynamical Systems and Variational Inequalities with Applications, volume 2. Springer Science & Business Media, 2012.
K. B. Nakshatrala and A. J. Valocchi. Non-negative mixed finite element formulations for a tensorial
diffusion equation. Journal of Computational Physics, 228:6726–6752, 2009.
K. B. Nakshatrala, M. K. Mudunuru, and A. J. Valocchi. A numerical framework for diffusioncontrolled bimolecular-reactive systems to enforce maximum principles and the non-negative
constraint. Journal of Computational Physics, 253:278–307, 2013.
K. B. Nakshatrala, H. Nagarajan, and M. Shabouei. A numerical methodology for enforcing maximum principles and the non-negative constraint for transient diffusion equations. Communications in Computational Physics, 19:53–93, 2016.
54

R. K. Pal, R. Abedi, A. Madhukar, and R. B. Haber. Adaptive spacetime discontinuous Galerkin
method for hyperbolic advection-diffusion with a non-negativity constraint. International Journal
for Numerical Methods in Engineering, 105:963–989, 2016.
G. S. Payette, K. B. Nakshatrala, and J. N. Reddy. On the performance of high-order finite elements
with respect to maximum principles and the non-negative constraint for diffusion-type equations.
International Journal for Numerical Methods in Engineering, 91:742–771, 2012.
C. Le Potier. Finite volume monotone scheme for highly anisotropic diffusion operators on unstructured triangular meshes. Comptes Rendus Mathematique, 341:787–792, 2005.
L. Qi. Convergence analysis of some algorithms for solving nonsmooth equations. Mathematics of
Operations Research, 18:227–244, 1993.
L. Qi and J. Sun. A nonsmooth version of Newton’s method. Mathematical Programming, 58:
353–368, 1993.
F. Rathgeber, G. R. Markall, L. Mitchell, N. Loriant, D. A. Ham, C. Bertolli, and P. H. J. Kelly.
Pyop2: A high-level framework for performance-portable simulations on unstructured meshes. In
High Performance Computing, Networking Storage and Analysis, SC Companion:, pages 1116–
1123, Los Alamitos, CA, USA, 2012. IEEE Computer Society. ISBN 978-1-4673-3049-7. doi: 10.
1109/SC.Companion.2012.134. URL http://dx.doi.org/10.1109/SC.Companion.2012.134.
F. Rathgeber, D. A. Ham, L. Mitchell, M. Lange, F. Luporini, A. T. T. McRae, G. T. Bercea, G. R.
Markall, and P. H. J. Kelly. Firedrake: automating the finite element method by composing
abstractions. Submitted to ACM TOMS, 2015. URL http://arxiv.org/abs/1501.01809.
P. A. Raviart and J. M. Thomas. A mixed finite element method for 2-nd order elliptic problems.
In Mathematical aspects of finite element methods, pages 292–315. Springer, 1977.
B. Riviére and M. Wheeler. Discontinuous Galerkin methods for and transport in porous media.
Communications in Numerical Methods in Engineering, 18:63–68, 2002.
J.F. Rodrigues. Obstacle Problems in Mathematical Physics, volume 134. Elsevier, 1987.
M. E. Rognes, R. C. Kirby, and A. Logg. Efficient assembly of h(div) and h(curl) conforming finite
elements. SIAM Journal on Scientific Computing, 31(6):4130–4151, 2009.
M. Sala, J. J. Hu, and R. S. Tuminaro. ML3.1 Smoothed Aggregation User’s Guide. Technical
Report SAND2004-4821, Sandia National Laboratories, 2004.
K. Shahbazi. Short note: An explicit expression for the penalty parameter of the interior penalty
method. Journal of Computational Physics, 205(2):401–407, May 2005. ISSN 0021-9991. doi:
10.1016/j.jcp.2004.11.017. URL http://dx.doi.org/10.1016/j.jcp.2004.11.017.
Z. Sheng and G. Yuan. A new nonlinear finite volume scheme preserving positivity for diffusion
equations. Journal of Computational Physics, 315:182–193, 2016.
A. Signorini. Sopra alcune questioni di statica dei sistemi continui. Annali della Scuola Normale
Superiore di Pisa-Classe di Scienze, 2(2):231–251, 1933.
A. Signorini. Questioni di elasticità non linearizzata e semilinearizzata. Rendiconti Di Matematica
e Delle Sue Applicazioni, V. Serie, 18:95–139, 1959. ISSN 0034-4427.
F. I. Stalkup. Miscible displacement. Society of Petroleum Engineers, 1983.
D. Z. Turner, K. B. Nakshatrala, M. J. Martinez, and P. K. Notz. Modeling subsurface water
resource systems involving heterogeneous porous media using the variational multiscale formulation. Journal of Hydrology, 428-429:1–14, 2011.
M. Ulbrich. Semismooth Newton Methods for Variational Inequalities and Constrained Optimization
Problems in Function Spaces. Society for Industrial and Applied Mathematics, Philadelphia, PA,
55

USA, 2011.
D. A. von der Schulenburg, T. R. R. Pintelon, C. Picioreanu, M. C.M. Van Loosdrecht, and M. L.
Johns. Three-dimensional simulations of biofilm growth in porous media. AIChE Journal, 55(2):
494–504, 2009.
C. E. Knutson C. J. Werth and A. J. Valocchi. Pore-scale simulation of biomass growth along the
transverse mixing zone of a model two-dimensional porous medium. Water Resources Research,
41(7), 2005.
X. H. Zhao, K. H. Kwek, J. B. Li, and K. L. Huang. Chaotic and resonant streamlines in the ABC
flow. SIAM Journal on Applied Mathematics, 53:71–77, 1993.

56

