# An automatic-adaptivity stabilized finite element method via residual
minimization for heterogeneous, anisotropic advection-diffusion-reaction
problems

Roberto J. Cier<EMAIL_ADDRESS>Sergio Rojas  Victor M. Calo  School of
Civil and Mechanical Engineering, Curtin University, Kent Street, Bentley,
Perth, WA 6102, Australia  School of Earth and Planetary Sciences, Curtin
University, Kent Street, Bentley, Perth, WA 6102, Australia  Mineral
Resources, Commonwealth Scientific and Industrial Research Organisation
(CSIRO), Kensington, Perth, WA 6152, Australia

###### Abstract

In this paper, we describe a stable finite element formulation for advection-
diffusion-reaction problems that allows for robust automatic adaptive
strategies to be easily implemented. We consider locally vanishing,
heterogeneous, and anisotropic diffusivities, as well as advection-dominated
diffusion problems. The general stabilized finite element framework was
described and analyzed in [1] for linear problems in general, and tested for
pure advection problems. The method seeks for the discrete solution through a
residual minimization process on a proper stable discontinuous Galerkin (dG)
dual norm. This technique leads to a saddle-point problem that delivers a
stable discrete solution and a robust error estimate that can drive mesh
adaptivity. In this work, we demonstrate the efficiency of the method in
extreme scenarios, delivering stable solutions. The quality and performance of
the solutions are comparable to classical discontinuous Galerkin formulations
in the respective discrete space norm on each mesh. Meanwhile, this technique
allows us to solve on coarse meshes and adapt the solution to achieve a user-
specified solution quality.

###### keywords:

stabilized finite elements , residual minimization , inf-sup stability ,
advection-diffusion-reaction , adaptive mesh refinement

††journal: Finite Elements in Analysis and Design

##  1 Introduction

Advection-diffusion-reaction problems arise in a wide range of phenomena
relevant to many areas of applied physics and engineering, such as flow in
porous media (e.g., reservoir engineering [2, 3] and groundwater flow [3, 4])
and drug delivery [5, 6, 7]. These processes generally involve heterogeneous
and highly anisotropic diffusion tensors, representing varying material
properties (e.g., permeability, porosity) in the domain [8, 9, 10]. Thus, the
accuracy and stability of the numerical approximation have been the focus of
intense research for several decades. Moreover, in advection-dominated
regimes, this problem behaves as a hyperbolic partial differential equation
(PDE), implying that an inaccurate numerical approximation could produce non-
physical oscillatory discrete solutions on coarse meshes.

Traditionally, stabilization terms modify standard finite elements to improve
the properties of the discrete solution. A plethora of conforming numerical
schemes exist which add extra scaled residual terms [11, 12, 13], and the
residual-free bubble approximation (RFB) [14]. Despite these efforts, the
stabilization process still represents a significant challenge for the
scientific community due to its relevance in the outcome of numerical
simulations. For instance, the discrete solution quality and solver
performance strongly depend on the appropriate selection of the penalization
(scaling) parameters. Alternatively, non-conforming schemes build stability
differently; for example, the discontinuous Galerkin (dG) class of methods,
achieve stability by enforcing an element-by-element discretization and
introducing a suitable choice of numerical traces [15, 16, 17, 18, 19, 20].
Other conforming stabilized formulations are the minimal residual methods,
such as the Least Squares Finite Element Method (LS-FEM) [21] or the
Discontinuous Petrov-Galerkin (DPG) method [22, 23, 24, 25, 26, 27, 28, 29],
which minimize the discrete residual with respect to an artificial energy norm
and, thus, achieve the sought stability.

Recently, [1] introduced, analyzed, and numerically tested a stabilized finite
element formulation in abstract form for any linear partial differential
equation system. The procedure minimizes the discrete residual of a conforming
trial space in the dual norm of a suitable dG test space in which the trial
space is a proper subspace. Hence, this formulation endows the discrete,
continuous solution with the stability properties of the dG formulation that
defines the dual norm. The method expresses the constrained residual
minimization as a saddle-point problem that delivers a conforming
approximation and an on-the-fly error estimate that drives automatic mesh
refinement. Besides, since the formulation builds on the non-conformity of the
underlying dG method, this attribute, in turn, allows us to consider strong
norms for the test space when the trial space has high regularity. Measuring
the error in stronger norms is an outstanding feature of this technology when
compared against alternative ones such as the LS-FEM and DPG methods. This
non-conforming residual minimization successfully tackled the analysis of
compaction banding in geomaterials [30], the constraint enforcement in
advection-dominated problems [31], and incompressible flows [32, 33], as well
as goal-oriented adaptive mesh refinements [34].

In this paper, we apply this stabilized finite element method to advection-
diffusion-reaction problems and show the impact of the norm selection on the
quality of the resulting solutions. We consider different possible scenarios
for this class of elliptic problems, such as advection-dominated diffusion,
and heterogeneous, with locally vanishing or highly-anisotropic diffusivities.
Section 2 describes the advection-diffusion-reaction problem and a
discontinuous Galerkin formulation that allows us to solve the problem via
residual minimization. Section 3 presents the method and the key ingredients
for its formulation, and applies it to the model problem. Finally, Section 4
discusses several numerical examples that show the efficiency and robustness
of this formulation.

##  2 The advection-diffusion-reaction problem

Let ${\Omega\subset\mathbb{R}^{d}}$, with dimension ${d=2,3}$, be an open and
bounded Lipschitz domain with boundary ${\Gamma:=\partial\Omega}$, and outward
unit normal vector $\mathbf{n}$. Using the standard notation of Hilbert and
Banach spaces, let ${K\in\left[L^{\infty}(\Omega)\right]^{d,d}}$ be a
diffusion tensor, to be symmetric and positive definite in $\Omega$. Let
${\mathbf{b}\in\left[L^{\infty}(\Omega)\right]^{d}}$ denote a divergence-free
(almost everywhere) advection coefficient, and ${\sigma\in
L^{\infty}(\Omega)}$ be a reactive coefficient. We write the advection-
diffusion-reaction problem as follows:

$\left\\{\begin{array}[]{l}\text{Find }u\text{ such that:}\vskip 3.0pt plus
1.0pt minus 1.0pt\\\ \begin{aligned} -\mathrm{div}\left(K\nabla
u\right)-\mathbf{b}\cdot\nabla u+\sigma\,u&=f,&&\text{ in }\Omega,\vskip 3.0pt
plus 1.0pt minus 1.0pt\\\ u&=0,&&\text{ on
}\Gamma,\end{aligned}\end{array}\right.$ (1)

where ${f\in L^{2}(\Omega)}$ denotes an spatial source.

In the present work, we focus on two different types of advection-diffusion-
reaction problems:

  * •

Advection-dominated problems, that is, problems where
${{0<\|K\|_{\infty},{\|\sigma\|_{\infty}}<<\|\mathbf{b}\|_{\infty}}}$. This
kind of problems lead unstable solutions when using the standard FEM on coarse
meshes.

  * •

Highly heterogeneous and anisotropic diffusion problems, that is, problems
where the diffusion takes locally small values, leading to advection-dominated
regimes.

These two scenarios lead to sharp inner and boundary layers, which are
difficult to capture with standard FEM formulation as they induce spurious
oscillations (see [35, 36]).

###  2.1 Continuous weak variational formulation

The weak formulation of (1) reads:

$\left\\{\begin{array}[]{l}\text{Find }u\in H^{1}_{0}(\Omega),\text{ such
that:}\vskip 3.0pt plus 1.0pt minus 1.0pt\\\ b(u,v)=\ell(v),\quad\forall v\in
H^{1}_{0}(\Omega),\end{array}\right.$ (2)

with bilinear form $b(u,v)=(K\nabla u,\nabla
v)_{0,\Omega}+(\mathbf{b}\cdot\nabla u,v)_{0,\Omega}+(\sigma u,v)_{0,\Omega}$,
and linear form $\ell(v)=(f,v)_{0,\Omega}$, where ${(\cdot,\cdot)_{0,\Omega}}$
denotes the $L^{2}$-scalar product in $\Omega$. In what follows, we assume
that there is a real number ${\sigma_{0}>0}$ such that

$\displaystyle\sigma-\frac{1}{2}\nabla\cdot\mathbf{b}\geq\sigma_{0}\quad\textrm{a.e.
in }\Omega.$ (3)

Furthermore, we assume that the smallest eigenvalue of $K$ is bounded from
below by a positive constant $K_{0}$. Then, owing to the Lax-Milgram Lemma,
problem (2) is well-posed (see e.g., [4]).

###  2.2 Discrete setting

Let $\\{\mathscr{P}_{h}\\}$ be a family of simplicial meshes of $\Omega$ and,
for simplicity, we assume that any mesh exactly represents $\Omega$ in
$\mathscr{P}_{h}$, that is, $\Omega$ is a polygon or a polyhedron. Let $T$ be
a generic element in $\mathscr{P}_{h}$, and denote by $\partial T$ its
boundary, by $h_{T}$ its diameter, and by $\mathbf{n}_{T}$ its outward unit
normal. We set ${h=\max_{T\in\mathscr{P}_{h}}h_{T}}$. We define the classical
dG approximation space as

$\displaystyle V_{h}:=\\{v_{h}\in L^{2}(\Omega)\ |\ \forall
T\in\mathscr{P}_{h},v_{h}|_{T}\in\mathds{P}^{p}(T)\\},$ (4)

where $\mathds{P}^{p}(T)$ denotes the space of polynomial functions with
degree smaller or equal than $p$. We set the extended space
${V_{h,\\#}=H^{2}(\mathscr{P}_{h})+V_{h}}$ for convenience. We say that $F$ is
an “interior face” if there exist two elements
$\\{T^{-}(F),\,{T^{+}(F)\\}\in\mathscr{P}_{h}}$, such that ${T^{-}(F)\cap
T^{+}(F)=F}$ and $F$ has nonzero measure. We collect all faces $F$ of
$\mathscr{P}_{h}$ into the set
${\mathscr{S}_{h}=\bigcup_{T\in\mathscr{P}_{h}}F}$. We denote by
$\mathscr{S}^{\partial}_{h}$ the boundary skeleton
${\mathscr{S}^{\partial}_{h}=\mathscr{S}_{h}\cap\Gamma}$, and by
$\mathscr{S}^{0}_{h}$ the internal skeleton
${\mathscr{S}^{0}_{h}=\mathscr{S}_{h}\backslash\Gamma}$. Over each
$F\in\mathscr{S}_{h}$, we set $\mathbf{n}_{F}$ as a predefined unit normal,
oriented from $T^{-}(F)$ to $T^{+}(F)$, being coincident with $\mathbf{n}$
when ${F\in\mathscr{S}_{h}^{\partial}}$, and $h_{F}$ as the diameter of the
face $F$. On interior faces, any function ${v_{h}\in V_{h}}$ is two-valued,
with values $v_{h}^{+}$ and $v_{h}^{-}$, defined with respect to the
predefined normal $\mathbf{n}_{F}$. Thus, the jump $[\\![v_{h}]\\!]_{F}$ and
the weighted average $\\{\\{v_{h}\\}\\}_{\omega}$ functions are defined as:

$\displaystyle[\\![v_{h}]\\!]_{F}:=v_{h}^{+}-v_{h}^{-},$
$\displaystyle\\{\\{v_{h}\\}\\}_{\omega}:=\omega^{-}v_{h}^{-}+\omega^{+}v_{h}^{+},$

where the weights satisfy ${\omega^{-}+\omega^{+}=1}$, with
${\omega^{-},\omega^{+}\geq 0}$. In particular, when considering heterogeneous
tensorial diffusivities, we choose the weights accounting the diffusivity
structure:

$\displaystyle\omega^{-}=\frac{\delta_{Kn}^{+}}{\delta_{Kn}^{+}+\delta_{Kn}^{-}},\quad\quad\omega^{+}=\frac{\delta_{Kn}^{-}}{\delta_{Kn}^{+}+\delta_{Kn}^{-}},$
(5)

with $\delta^{\mp}_{Kn}=\mathbf{n}_{F}\cdot K^{\mp}\mathbf{n}_{F}$ if
$F\in\mathscr{S}_{h}^{0}$, and $\delta_{Kn}=\mathbf{n}_{F}\cdot
K\mathbf{n}_{F}$ if $F\in\mathscr{S}_{h}^{\partial}$. When $K$ is a continuous
tensor (homogeneous diffusion), the weights reduce to
${\omega^{-}=\omega^{+}=1/2}$. Finally, on a boundary face
${F\in\mathscr{S}^{\partial}_{h}}$, we set
${[\\![v_{h}]\\!]_{F}=\\{\\{v_{h}\\}\\}_{F}=v_{h}|_{F}}$. We omit the
subscript $F$ in the jump and weighted average functions when there is no
ambiguity.

###  2.3 Discontinuous Galerkin variational formulation

In this section, we briefly discuss a stable discontinuous Galerkin (dG)
formulation for problem (1). It combines the Symmetry Weighted Interior
Penalty (SWIP) scheme [37, 4] that handles general diffusivities, combined
with the Upwinding (UPW) method [19, 38] that handles the advection-reaction
contribution.

Considering the discrete setting described in Section 2.2, the dG formulation
of problem (2) reads:

$\left\\{\begin{array}[]{l}\text{Find }u_{h}\in V_{h},\text{ such that:}\vskip
3.0pt plus 1.0pt minus 1.0pt\\\ \begin{aligned}
b_{h}(u^{\text{dG}}_{h},v_{h})&:=b_{h}^{\textrm{diff}}(u^{\text{dG}}_{h},v_{h})+b_{h}^{\textrm{adv}}(u^{\text{dG}}_{h},v_{h})=\ell_{h}(v_{h}),&&\forall
v_{h}\in V_{h},\end{aligned}\end{array}\right.$ (6)

where the bilinear and linear forms are

$\displaystyle b_{h}^{\textrm{diff}}(u_{h},v_{h})$
$\displaystyle:=\displaystyle\sum_{T\in\mathscr{P}_{h}}(K\nabla
u_{h}\,,\,\nabla v_{h})_{0,T}$
$\displaystyle+\displaystyle\sum_{F\in\mathscr{S}_{h}}\left[\displaystyle\,-\left([\\![u_{h}]\\!]\,,\,\mathbf{n}_{F}\cdot\\{\\{K\nabla
v_{h}\\}\\}_{\omega}\right)_{0,F}-\displaystyle\left(\mathbf{n}_{F}\cdot\\{\\{K\nabla
u_{h}\\}\\}_{\omega}\,,\,[\\![v_{h}]\\!]\right)_{0,F}+\gamma_{F}\left([\\![u_{h}]\\!]\,,\,[\\![v_{h}]\\!]\right)_{0,F}\right],\vskip
3.0pt plus 1.0pt minus 1.0pt$ $\displaystyle b_{h}^{\textrm{adv}}(w,v)$
$\displaystyle:=\displaystyle\sum_{T\in\mathscr{P}_{h}}(\,\mathbf{b}\cdot\nabla
u_{h}+\sigma\,u_{h}\,,\,v_{h})_{0,T}+\displaystyle\sum_{F\in\mathscr{S}^{\partial}_{h}}\left((\,\mathbf{b}\cdot\mathbf{n}_{F})^{\ominus}u_{h},v_{h}\right)_{0,F}$
$\displaystyle+\sum_{F\in\mathscr{S}^{0}_{h}}\left[\dfrac{1}{2}\left(\,\left|\,\mathbf{b}\cdot\mathbf{n}_{F}\right|\,[\\![u_{h}]\\!]\,,\,[\\![v_{h}]\\!]\right)_{0,F}-\left(\,\mathbf{b}\cdot\mathbf{n}_{F}\,[\\![u_{h}]\\!]\,,\,\\{\\{v_{h}\\}\\}\right)_{0,F}\right],\vskip
3.0pt plus 1.0pt minus 1.0pt$

and

$\displaystyle\ell_{h}(v_{h})$
$\displaystyle:=\sum_{T\in\mathscr{P}_{h}}(f,v_{h})_{0,T}$

In the above, $(\cdot)^{\ominus}$ denotes the negative part of $x$ (i.e.,
${x^{\ominus}:=\frac{1}{2}\left(|x|-x\right)}$ for any real number $x$). For
problems with diffusion, there exist many suitable choices for the penalty
parameter $\gamma_{F}$ (e.g., [39, 40] analyzed the penalty parameters and
their dependence on the polynomial order of approximation; [4] defined and
analyzed the impact of $\gamma_{K}$ as the harmonic average of the “normal”
permeabilities; [41], introduced a mesh-dependent penalty parameter). In this
work, following [4, 39, 42], we set the penalty parameter $\gamma_{F}$ as
${\gamma_{F}=\displaystyle\eta\,{\gamma_{K}}}$. Here, $\eta>0$ represents an
element-wise parameter defined as:

$\displaystyle\forall F\in\mathscr{S}^{0}_{h},$ $\displaystyle\eta$
$\displaystyle=\frac{1}{2}\frac{(p+1)(p+d)}{d}\left(\frac{\mathscr{A}(\partial
T^{+})}{\mathscr{V}(T^{+})}+\frac{\mathscr{A}(\partial
T^{-})}{\mathscr{V}(T^{-})}\right),$ (7) $\displaystyle\forall
F\in\mathscr{S}^{\partial}_{h},$ $\displaystyle\eta$
$\displaystyle=\frac{(p+1)(p+d)}{d}\frac{\mathscr{A}(\partial
T)}{\mathscr{V}(T)},$ (8)

where $p$ is the polynomial order of the test space $V_{h}$, $d$ is the
dimension, and $\mathscr{A}$ and $\mathscr{V}$ denote area and volume,
respectively, for $d=3$, and length and area, respectively, for $d=2$. We
define ${\gamma_{K}}$ as follows:

$\displaystyle\forall F\in\mathscr{S}^{0}_{h},$ $\displaystyle\gamma_{K}$
$\displaystyle=(\omega^{-})^{2}\delta^{-}_{Kn}+(\omega^{+})^{2}\delta^{+}_{Kn},$
(9) $\displaystyle\forall F\in\mathscr{S}^{\partial}_{h},$
$\displaystyle\gamma_{K}$ $\displaystyle=\delta_{Kn},$ (10)

thus, recalling weights definition (5), we derive that

$\displaystyle\forall F\in\mathscr{S}_{h}^{0},$ $\displaystyle\gamma_{K}$
$\displaystyle=\frac{\delta_{Kn}^{+}\delta_{Kn}^{-}}{\delta_{Kn}^{+}+\delta_{Kn}^{-}}.$
(11)

When we consider scalar diffusivities, that is, $K=\varkappa I$ for some
scalar function ${\varkappa:\Omega\rightarrow\mathbb{R}}$, we recover the
symmetric interior penalty (IP or SIP) method [43, 44, 45] given that the
penalty parameters reduces to $\gamma_{F}=\eta\varkappa$.

######  Remark 1.

The dG formulation allows us to weakly impose non-homogeneous Dirichlet
boundary conditions through modifying the right-hand side of (6) as follows:
if we look for a solution of problem (1) satisfying $u=g_{D}$ on $\Gamma$,
being ${g_{D}\in H^{1/2}(\Gamma)}$ a boundary source, then we rewrite the
linear form $\ell_{h}(v_{h})$ as:

$\displaystyle\displaystyle\ell_{h}(v_{h})$
$\displaystyle:=\displaystyle\sum_{T\in\mathscr{P}_{h}}(f,v_{h})_{0,T}+\sum_{F\in\mathscr{S}_{h}^{\partial}}\left[\,-\left(g_{D},\mathbf{n}_{F}\cdot
K\nabla
v_{h}\right)_{0,F}+\gamma_{F}\left(g_{D},v_{h}\right)_{0,F}+\left((\,\mathbf{b}\cdot\mathbf{n}_{F})^{\ominus}g_{D},v_{h}\right)_{0,F}\right].\vskip
3.0pt plus 1.0pt minus 1.0pt$

###  2.4 Discrete norms, well-posedness and a priori error estimate

We provide the discrete space $V_{h}$ with the following norm:

$\displaystyle\|w\|^{2}_{V_{h}}:=\|w\|^{2}_{\text{adv}}+\|w\|^{2}_{\text{diff}}$
(12)

with

$\displaystyle\|w\|^{2}_{\text{adv}}$
$\displaystyle:=\displaystyle\|w\|^{2}_{0,\Omega}+\frac{1}{2}\|\,|\,\mathbf{b}\cdot\mathbf{n}\,|^{\frac{1}{2}}w\|^{2}_{0,\Gamma}+\displaystyle\frac{1}{2}\sum_{F\in\mathscr{S}^{0}_{h}}\left(|\,\mathbf{b}\cdot\mathbf{n}_{F}|\,[\\![w]\\!],[\\![w]\\!]\right)_{0,F}+\displaystyle\sum_{T\in\mathscr{P}_{h}}h_{T}\|\,\mathbf{b}\cdot\nabla{w}\,\|^{2}_{0,T},$
$\displaystyle\|w\|^{2}_{\text{diff}}$
$\displaystyle:=\displaystyle\|\,\kappa\nabla
w\,\|^{2}_{0,\Omega}+\sum_{F\in\mathscr{S}_{h}}\left(\gamma_{F}[\\![w]\\!],[\\![w]\\!]\right)_{0,F},$

where $\kappa$ denotes the (unique) symmetric positive definite tensor-valued
field such that $\kappa^{2}=K$ a.e. in $\Omega$. Following [4], we define:

$|w|^{2}_{V_{h,\beta}}:=\displaystyle\sum_{T\in\mathscr{P}_{h}}h_{T}\|\,\mathbf{b}\cdot\nabla{w}\,\|^{2}_{0,T},$
(13)

which represents the last component of the norm $\|w\|^{2}_{\text{adv}}$ that
controls the advective derivative error for small diffusivities. Finally, we
define the following extended norm $\|w\|_{V_{h},\\#}$:

$\|w\|_{V_{h},\\#}:=\displaystyle\|w\|_{V_{h}}+\left(\sum_{T\in\mathscr{P}_{h}}\|w\|^{2}_{0,\partial
T}\right)^{\frac{1}{2}}+\displaystyle\left(\sum_{T\in\mathscr{P}_{h}}h_{T}\|\kappa\nabla
w\|^{2}_{0,\partial T}\right)^{\frac{1}{2}}.$ (14)

In the remainder, the symbol $\lesssim$ indicates an inequality involving a
positive constant $C$ independent of the mesh and diffusivity. Considering the
norms above defined, the following theorem holds true (see [4, §3 & §4]):

######  Theorem 1 (Well-posedness and a priori error estimate of the dG
formulation).

The following propositions hold true:

  1. (a)

Inf-sup stability: There exists a constant
${C_{\emph{sta}}=C\Delta_{K}^{-1}}$, with ${C>0}$, uniform with respect to the
mesh size, such that:

$\displaystyle\sup_{v_{h}\in
V_{h}\backslash\\{0\\}}\frac{b_{h}(w_{h},v_{h})}{\|v_{h}\|_{V_{h}}}\geq
C_{\emph{sta}}\|w_{h}\|_{V_{h}},\quad\forall w_{h}\in V_{h},$

where $\forall\,T\in\mathscr{P}_{h}$,
${\Delta_{K}=\max_{T\in\mathscr{P}_{h}}\Delta_{K,T}}$, and

$\displaystyle\quad\Delta_{K,T}=\left\\{\begin{array}[]{l l}\vskip 3.0pt plus
1.0pt minus 1.0pt1&\text{if}\ \
\|\mathbf{b}\|_{[L^{\infty}(T)]^{d}}\gtrsim\dfrac{\lambda_{M,T}}{h_{T}},\\\
\dfrac{\lambda_{M,T}}{\lambda_{m,T}}&\text{otherwise},\end{array}\right.$

with ${\lambda_{M,T}}$ and ${\lambda_{m,T}}$ as the maximum and minimum
eigenvalues of $K|_{T}$, respectively.

  2. (b)

Boundedness: There exists a mesh-independent constant $C_{\emph{bnd}}<\infty$,
such that:

$\displaystyle b_{h}(z,v_{h})\leq
C_{\emph{bnd}}\|z\|_{V_{h},\\#}\|v_{h}\|_{V_{h}},\quad\forall(z,v_{h})\in
V_{h,\\#}\times V_{h}.$

  3. (c)

Consistency: Let $u$ be the solution of (2). If ${u\in H^{1}_{0}(\Omega)\cap
H^{2}(\mathscr{P}_{h})}$, then

$\displaystyle b_{h}(u,v_{h})=\ell_{h}(v_{h}),\quad\forall v_{h}\in V_{h}.$

Henceforth, $u^{\emph{dG}}$, solution of problem (6), is unique (cf. [4, 38]).
Moreover, the following a priori error estimate is satisfied

$\displaystyle\inf_{y_{h}\in
V_{h}}\|u-y_{h}\|_{V_{h}}\leq\|u-u^{\emph{dG}}_{h}\|_{V_{h}}\leq\left(1+\frac{C_{\emph{bnd}}}{C_{\emph{sta}}}\right)\inf_{v_{h}\in
V_{h}}\|u-v_{h}\|_{V_{h},\\#}.$ (15)

We relate the convergence rates of both left and right-hand sides of the error
estimate (15) through the following definition (cf. [38], §1.4.4):

######  Definition 1 (Optimality, quasi-optimality and suboptimality of the
error estimate).

The error estimate (15) is

  1. 1.

_optimal_ if $\|\cdot\|_{V_{h}}\simeq\|\cdot\|_{V_{h},\\#}$,

  2. 2.

_quasi-optimal_ if the two norms are different, but the lower and upper bounds
in (15) converge, for smooth $u$, at the same convergence rate as
$h\rightarrow 0$,

  3. 3.

_suboptimal_ if the upper bound converges at a slower rate than the other
bound.

##  3 Residual minimization based on discontinuous Galerkin formulations

###  3.1 Method description

We now describe the stabilized finite element formulation via residual
minimization on dual discontinuous Galerkin (dG) norms devised in [1] for the
advection-diffusion-reaction problem context. In (6), $u^{\text{dG}}_{h}$
discretely approximates the solution of (2) belonging to a discontinuous
discrete space. In [1], the formulation seeks for an approximation of $u$ in a
discrete space that may possess additional properties, for instance,
continuity and, possibly, higher smoothness. Rather than solving problem (6),
the dG-based residual minimization method implies:

  1. (a)

Consider the subspace $U_{h}=V_{h}\cap U\subset V_{h}$ (i.e., standard finite
element functions).

  2. (b)

Obtain ${u_{h}\in U_{h}}$ by minimizing the residual:

$\left\\{\begin{array}[]{l}\text{Find }u_{h}\in U_{h}\subset V_{h},\text{ such
that:}\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
u_{h}=\displaystyle\operatorname*{arg\,min}_{z_{h}\in
U_{h}}\dfrac{1}{2}\|\ell_{h}-B_{h}\,z_{h}\|^{2}_{V_{h}^{\ast}}=\displaystyle\operatorname*{arg\,min}_{z_{h}\in
U_{h}}\dfrac{1}{2}\|R^{-1}_{V_{h}}(\ell_{h}-B_{h}z_{h})\|^{2}_{V_{h}},\end{array}\right.$
(16)

where $B_{h}:V_{h},\\#\rightarrow V_{h}^{\ast}$ is:

$\displaystyle\langle B_{h}w_{h},v_{h}\rangle_{V_{h}^{\ast}\times
V_{h}}:=b_{h}(w_{h},v_{h}),$ (17)

and $R^{-1}_{V_{h}}$ denotes the inverse of the Riesz map:

$\displaystyle R_{V_{h}}:\ $ $\displaystyle V_{h}\rightarrow
V_{h}^{\ast}\vskip 3.0pt plus 1.0pt minus 1.0pt$ (18) $\displaystyle
v_{h}\rightarrow\langle R_{V_{h}}y_{h},v_{h}\rangle_{V_{h}^{\ast}\times
V_{h}}:=(y_{h},v_{h})_{V_{h}},$ (19)

with $(\cdot,\cdot)_{V_{h}}$ denoting the inner product that induces the
discrete norm ${\|\cdot\|_{V_{h}}}$ (i.e.,
${\|\cdot\|_{V_{h}}=(\cdot,\cdot)^{1/2}_{V_{h}}}$).

The second equality in (16) follows as the Riesz map (18) is an isometric
isomorphism, therefore $\|\cdot\|_{V_{h}^{\ast}}$ is equivalent to
$\|R^{-1}_{V_{h}}(\cdot)\|_{V_{h}}$. Thus, problem (16) is equivalent to the
following saddle-point problem (see [46]):

$\left\\{\begin{array}[]{l}\text{Find }(\varepsilon_{h},u_{h})\in V_{h}\times
U_{h},\text{ such that:}\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
\begin{array}[]{rcll}(\varepsilon_{h},v_{h})_{V_{h}}+b_{h}(u_{h},v_{h})&=&\ell_{h}(v_{h}),&\quad\forall
v_{h}\in V_{h},\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
b_{h}(z_{h},\varepsilon_{h})&=&0,&\quad\forall z_{h}\in
U_{h},\end{array}\end{array}\right.$ (20)

where $\varepsilon_{h}$ is a residual representative in $V_{h}$. Indeed, the
first identity in (20) implies that:

$\displaystyle\varepsilon_{h}=R^{-1}_{V_{h}}(\ell_{h}-B_{h}u_{h})=R^{-1}_{V_{h}}B_{h}(u_{h}^{\textrm{dG}}-u_{h}),$
(21)

where the second identity in (21) comes from (6).

The saddle-point formulation has several desirable properties for numerical
approximations. Firstly, the matrix associated with the inner product
$(\varepsilon_{h},v_{h})_{V_{h}}$ in (20) is always symmetric and positive-
definite, independently of the nature of the chosen dG formulation; thus,
several well-known iterative solvers are effective on the resulting saddle-
point problem. Moreover, the discrete approximation $u_{h}\in U_{h}$ inherits
the discrete stability of the dG formulation. Finally, $\varepsilon_{h}\in
V_{h}$ is a robust error representation, which under an adequate saturation
assumption becomes a reliable representative. Below, we summarize the last
properties (see [1] for details):

######  Theorem 2 (Well posedness and a priori error bound estimates for the
saddle-point problem).

The solution $(\varepsilon_{h},u_{h})\in V_{h}\times U_{h}$ of the saddle-
point problem (20) is unique and the following a priori bound applies:

$\displaystyle\|\varepsilon_{h}\|\leq\|\ell_{h}\|_{V_{h}^{*}}\quad\text{and}\quad\|u_{h}\|_{V_{h}}\leq\frac{1}{C_{\emph{sta}}}\|\ell_{h}\|_{V_{h}^{*}},$
(22)

and the following error estimate holds true:

$\|u-u_{h}\|_{V_{h}}\leq\left(1+\dfrac{C_{\emph{bnd}}}{C_{\emph{sta}}}\right)\Delta_{K}\inf_{v_{h}\in
U_{h}}\|u-v_{h}\|_{V_{h},\\#},$ (23)

where ${u\in V_{\\#}}$ represents the exact solution to the continuous problem
(1).

######  Proposition 1 (Efficiency of the residual representative).

Under the same hypotheses of Theorem 2, the following holds:

$\|\varepsilon_{h}\|_{V_{h}}\leq C_{\emph{bnd}}\,\|u-u_{h}\|_{V_{h},\\#}.$
(24)

######  Assumption 1 (Saturation).

Let ${u_{h}\in U_{h}}$ be the second component of the pair
${(\varepsilon_{h},u_{h})\in V_{h}\times U_{h}}$ solving the saddle-point
problem (20). Let ${u^{\emph{dG}}_{h}\in V_{h}}$ be the unique solution to
(6). There exists a real number $\delta\in[0,1)$, uniform with respect to the
mesh size, such that
${\|u-u^{\emph{dG}}_{h}\|_{V_{h}}\leq\delta\|u-u_{h}\|_{V_{h}}}$.

######  Proposition 2 (Reliability of the residual representative).

Let ${u_{h}\in U_{h}}$ be the second component of ${(\varepsilon_{h},u_{h})\in
V_{h}\times U_{h}}$ solving the saddle-point problem (20). Let
${u^{\text{dG}}_{h}\in V_{h}}$ be the unique solution to (6). Then the
following holds true:

$\|u_{h}-u^{\emph{dG}}_{h}\|_{V_{h}}\leq\dfrac{1}{C_{\emph{sta}}}\|\varepsilon_{h}\|_{V_{h}}.$
(25)

Moreover, if the saturation Assumption 1 is satisfied, then the following a
posteriori error estimate holds:

$\|u-u_{h}\|_{V_{h}}\leq\dfrac{1}{(1-\delta)C_{\emph{sta}}}\|\varepsilon_{h}\|_{V_{h}}.$
(26)

We can now state the requirements for having an efficient and reliable
residual representative. From (24) and (26), we have:

$\|u-u_{h}\|_{V_{h}}\lesssim\|\varepsilon_{h}\|_{V_{h}}\lesssim\|u-u_{h}\|_{V_{h},\\#}.$
(27)

Thus, to ensure the usefulness of the residual representative, we need at
least quasi-optimality, which means that the left-hand side should decay at
the same rate as the right-hand side of (27) (see Definition 1).

###  3.2 Convergence rates

In the context of the above framework, we can recover some insightful results
related to the convergence rates. Similarly as done in [1, Appendix B] (see
also [47, 48, 49]), we can prove:

$\inf_{v_{h}\in U_{h}}\|u-v_{h}\|_{V_{h},\\#}\,\lesssim\,\inf_{v_{h}\in
V_{h}}\|u-v_{h}\|_{V_{h},\\#}.$ (28)

This proves, in particular, that the residual minimization method delivers a
discrete solution with the same quality as the one delivered by the dG
formulation, as consequence of (15) and (23). Thus, if the solution is regular
enough, it follows:

$\inf_{v_{h}\in V_{h}}\|u-v_{h}\|_{V_{h},\\#}\lesssim h^{p}.$ (29)

Additionally, from (29) we can deduce a bound for the advective component
${|\cdot|_{V_{h},\beta}}$, defined in (13). Indeed, we can easily infer that:

$\|\,\mathbf{b}\cdot\nabla(u-u_{h})\,\|^{2}_{0,T}\,\lesssim\,\|\,\kappa\nabla(u-u_{h})\,\|^{2}_{0,\Omega}\,\lesssim\,\|u-u_{h}\|^{2}_{V_{h},\\#},$
(30)

where $u$ and $u_{h}$ are defined as in Theorem 2. Combining (28), (30) and
(29), multiplying both sides by the mesh size, and taking the square root, we
can recover the optimal convergence rate for ${|u-u_{h}|_{V_{h,\beta}}}$,
which reads:

$|u-u_{h}|_{V_{h,\beta}}\,=\left(\displaystyle\sum_{T\in\mathscr{P}_{h}}h_{T}\|\,\mathbf{b}\cdot\nabla(u-u_{h})\,\|_{0,T}\right)^{1/2}\,\lesssim\,h^{1/2}\|u-u_{h}\|_{V_{h},\\#}\,\lesssim\,h^{p+1/2}.$

Similarly, we can deduce a bound for the $L^{2}$-norm error in the following
way:

$\|u-u_{h}\|_{0,\Omega}\,\lesssim\,\|u-u_{h}\|_{V_{h},\\#}\,\lesssim\,\inf_{v_{h}\in
V_{h}}\|u-v_{h}\|_{V_{h},\\#}\lesssim h^{p},$ (31)

which is suboptimal. Indeed, we show through numerical examples that the
method could deliver solutions with suboptimal $L^{2}$ convergence rates.
Nevertheless, we show that the optimal convergence can be recovered by
performing the energy-norm based adaptive mesh refinement strategy.

##  4 Numerical examples

In this section, we discuss some implementation aspects and describe the test
cases that demonstrate the performance of the method under a wide range of
challenging scenarios.

###  4.1 Implementation aspects

We use FEniCS [50] to perform all the numerical simulations. We show
convergence plots of the error measured in the $L^{2}$ and $V_{h}$ norms
versus the number of degrees of freedom (DOFs) (i.e.,
$\dim(U_{h})+\dim(V_{h})$ for (20)). As Section 2.3 describes, we use the
Symmetric Weighted Interior Penalty (SWIP) method for the diffusive part of
the bilinear form, which extends the classical Symmetric Interior Penalty
(SIP) formulation. In our experience, extensive numerical testing shows that
other dG formulations have similar computational cost and convergence rates;
thus, for the sake of brevity, we detail the performance of the method
uniquely using the SWIP formulation.

####  4.1.1 Adaptive mesh refinement

We use $\varepsilon_{h}\in V_{h}$ to estimate the error and to drive the
adaptive mesh refinement process [1]. We follow a standard adaptive procedure,
which considers an iterative loop where each level of refinement we perform
the following four analysis modules:

$\text{ SOLVE }\rightarrow\text{ ESTIMATE }\rightarrow\text{ MARK
}\rightarrow\text{ REFINE. }$

That is, given a mesh partition, we first solve the saddle-point problem (20).
Later, we use a localized version of the inner product (12) (evaluated in each
mesh cell $T$) as error indicator $E_{T}$

$E_{T}^{2}:=\|\varepsilon_{h}\|^{2}_{loc,T}+\frac{1}{2}|\varepsilon_{h}|^{2}_{loc,S}\,,$
(32)

with

$\begin{array}[]{rcl}\|\varepsilon_{h}\|^{2}_{loc,T}&:=&\,\displaystyle\|\varepsilon_{h}\|_{0,T}^{2}+\displaystyle\|\kappa\nabla\varepsilon_{h}\|_{T}^{2}+\displaystyle
h_{0,T}\|\
\mathbf{b}\cdot\nabla\varepsilon_{h}\|^{2}_{0,T}\displaystyle+\sum_{F\in\mathscr{S}^{\partial}_{h}}\left({\frac{1}{2}}\,|\mathbf{b}\cdot\mathbf{n}|+\gamma_{F}\right)\left(\varepsilon_{h},\varepsilon_{h}\right)_{0,F},\\\
|\varepsilon_{h}|^{2}_{loc,S}&:=&\displaystyle\sum_{F\in\mathscr{S}^{0}_{h}}\left({\frac{1}{2}}\,|\mathbf{b}\cdot\mathbf{n}|+\gamma_{F}\right)\left([\\![\varepsilon_{h}]\\!],[\\![\varepsilon_{h}]\\!]\right)_{0,F}.\end{array}$
(33)

We then mark elements for refinement using the Dörfler bulk-chasing criterion
[51] (i.e., select all elements for which the cumulative sum of the local
values $E_{T}$ in a decreasing order remains below a user-defined fraction of
the total estimated error $\|\varepsilon_{h}\|_{V_{h}}$). Herein, we set this
fraction to be 0.5 for $d=2$ and 0.25 for $d=3$. Finally, we bisect each
marked element [52] to obtain the refined mesh to use in the next step.

####  4.1.2 Iterative solver

The matrix system that problem (20) induces has the following form

$\left[\begin{array}[]{cc}G&B\\\
B^{T}&0\end{array}\right]\left[\begin{array}[]{c}{\boldsymbol{\varepsilon}}\\\
{\bf u}\end{array}\right]=\left[\begin{array}[]{c}\bf{L}\\\ {\bf
0}\end{array}\right].$ (34)

where the superindex $T$ denotes transpose. Following [1, §5], we apply the
iterative algorithm proposed in [53]. Denoting by $\widehat{G}$ a
preconditioner for the Gram matrix $G$, and by $\widehat{S}$ a preconditioner
for the reduced Schur complement $B^{T}\widehat{G}^{-1}B$, the iterative
scheme becomes

$\left[\begin{array}[]{c}{\boldsymbol{\varepsilon}}_{i+1}\\\ {\bf
u}_{i+1}\end{array}\right]=\left[\begin{array}[]{c}{\boldsymbol{\varepsilon}}_{i}\\\
{\bf u}_{i}\end{array}\right]+\left[\begin{array}[]{cc}\widehat{G}&B\\\
B^{T}&\widehat{C}\end{array}\right]^{-1}\left\\{\left[\begin{array}[]{c}\bf{L}\\\
{\bf 0}\end{array}\right]-\left[\begin{array}[]{cc}G&B\\\
B^{T}&0\end{array}\right]\left[\begin{array}[]{c}{\boldsymbol{\varepsilon}}_{i}\\\
{\bf u}_{i}\end{array}\right]\right\\},$ (35)

with $\widehat{C}=B^{T}\widehat{G}^{-1}B-\widehat{S}$. Let ${{\bf r}_{i}={\bf
L}-G\,{\boldsymbol{\varepsilon}}_{i}-B\,{\bf u}_{i}}$ and ${{\bf
s}_{i}=-B^{T}{\boldsymbol{\varepsilon}}_{i}}$. Also, let
${\boldsymbol{\varepsilon}}$ and ${\bf u}$ be the residuals at the outer
iteration $i$. Then, the scheme requires the resolution of two interior
problems for these increments:

$\boldsymbol{\eta}_{i+1}:={\bf u}_{i+1}-{\bf
u}_{i}=\widehat{S}^{-1}\left(B^{T}\left(\widehat{G}^{-1}{\bf
r}_{i}\right)-{\bf s}_{i}\right),$ (36)

and

$\boldsymbol{\delta}_{i+1}:={\boldsymbol{\varepsilon}}_{i+1}-{\boldsymbol{\varepsilon}}_{i}=\widehat{G}^{-1}\left({\bf
r}_{i}-B\,{\boldsymbol{\eta}}_{i+1}\right).$ (37)

We construct an accurate approximation of $G^{-1}$, as in practice, this
delivers the best good computational performance; that is, low-quality
approximations lead poor conditioning of the reduced Schur complement in (36).
Bank et al. [53] proposed a relaxed approximation for the matrix $G$;
nevertheless, for this stiffer problem, we use a scheme with one outer
iteration. That is, we approximate $G^{-1}$ with a sparse Cholesky
factorization (e.g., using the module “sksparse.cholmod”, see [54]). Moreover,
we choose precondition (36) with an approximate Schur complement built as
${\widehat{S}=B^{T}(\text{diag}(G))^{-1}B}$, where ${\text{diag}(G)}$ is the
main diagonal of $G$, and the inverse is approximated through the same
procedure used for $G$. We use the LGMRES algorithm (e.g., from Scipy sparse
linear algebra package) to solve. On the coarsest mesh, our initial guess is
zero, whereas, on the next refinements, our guess is the solution of the
previous level of refinement.

(a) Problem sketch

(b) 3D view of the adaptive solution for $p=1$

(c) Convergence plot in $L^{2}$ error norm. $\alpha=2/3$

(d) Convergence plot in $V_{h}$ error norm. $\alpha=2/3$

Figure 1: Convergence rates, Poisson problem in L-shape domain. $\alpha=2/3$

###  4.2 Pure diffusion on L-shape domain

As a first example, we consider the L-shape domain
${\Omega:=(-1,1)^{2}\setminus(-1,0]^{2}}$, and the following Poisson problem:

$\displaystyle-\Delta u$ $\displaystyle=0,\quad\textrm{in }\Omega,$ (38)
$\displaystyle u$ $\displaystyle=g_{D},\quad\textrm{on }\partial\Omega,$ (39)

where $g_{D}$ corresponds to the Dirichlet trace of the analytical function in
polar coordinates ${u(r,\theta)=r^{\alpha}\sin(\alpha\,\theta)}$, with
${\theta=3\pi/2}$ for our case. This particular problem is known as reentrant
corner problem [55], where the solution has a singularity at the corner, and
its solution belongs to ${H^{1+\alpha-\epsilon}}$, ${\forall\epsilon>0}$ [56].
The dG variational formulation we use in problem (38) is the formulation (6)
with $K=1$, ${\mathbf{b}=(0,0)^{T}}$, ${\sigma=0}$ and ${f=0}$. Figure 1 shows
the convergence plots for ${\alpha=2/3}$. A uniform refinement strategy cannot
deal with the corner singularity achieving similar convergence rates for
increasing polynomial orders. However, the adaptive stabilized methodology
overcomes this limitation, recovering optimal convergence rates on both
$L^{2}$ and $V_{h}$ error norms.

(a) Problem sketch

(b) Convergence plot in $L^{2}$ error norm

(c) Convergence plot in $V_{h}$ error norm

(d) Convergence plot in $V_{h,\beta}$ error norm

Figure 2: Heterogeneous diffusion problem.

###  4.3 Heterogeneous diffusion

We solve an advection-diffusion problem with discontinuous diffusion
coefficients, based on a test from [57]; this problem shows the method
performance with heterogeneous diffusion. We split the domain in two
subdomains, $\Omega_{1}=[0,\frac{1}{2}]\times[0,1]$ and
$\Omega_{2}=[\frac{1}{2},1]\times[0,1]$ and use a constant diffusivity tensor
in each subdomain

$\displaystyle K_{i}(x,y)=\begin{pmatrix}\epsilon_{i}(x)&0\\\
0&1.0\end{pmatrix}$

where $\epsilon_{i}(x)$ represent discontinuous values across the interface
$x=\frac{1}{2}$. We set $\epsilon_{1}(x)=1\times 10^{-2}$ and
$\epsilon_{2}(x)=1.0$ with $\mathbf{b}=(1,0)^{T}$, ${\sigma=0}$ and ${f=0}$.
For this parameter choice, the exact solution is an exponential with respect
to the $x$-coordinate (i.e., independent of the $y$-coordinate). At the
interface, the solution satisfies the following conditions:

$\displaystyle\displaystyle\lim_{x\rightarrow\frac{1}{2}^{-}}u(x,y)$
$\displaystyle=\displaystyle\lim_{x\rightarrow\frac{1}{2}^{+}}u(x,y),$
$\displaystyle\displaystyle\lim_{x\rightarrow\frac{1}{2}^{-}}-\epsilon(x)\partial_{x}u(x,y)$
$\displaystyle=\displaystyle\lim_{x\rightarrow\frac{1}{2}^{+}}-\epsilon(x)\partial_{x}u(x,y).$

We set $u(0,y)=0$, $u(1,y)=1$, and by consequence of the matching conditions,
we obtain

$\displaystyle
u\left(\frac{1}{2},y\right)=\left[\frac{u(0,y)\exp\left(\frac{1}{2\epsilon_{1}}\right)}{1-\exp\left(\frac{1}{2\epsilon_{1}}\right)}+\frac{u(1,y)}{1-\exp\left(\frac{1}{2\epsilon_{2}}\right)}\right]\left[\frac{\exp\left(\frac{1}{2\epsilon_{1}}\right)}{1-\exp\left(\frac{1}{2\epsilon_{1}}\right)}+\frac{1}{1-\exp\left(\frac{1}{2\epsilon_{2}}\right)}\right]^{-1}.$

Thus, the exact solution in each subdomain becomes:

$\displaystyle u_{1}(x,y)=$
$\displaystyle\frac{u\left(\frac{1}{2},y\right)-\exp\left(\frac{1}{2\epsilon_{1}}\right)u(0,y)+\left[u(0,y)-u\left(\frac{1}{2},y\right)\right]\exp\left(\frac{x}{\epsilon_{1}}\right)}{1-\exp\left(\frac{1}{2\epsilon_{1}}\right)}$
$\displaystyle u_{1}(x,y)=$
$\displaystyle\frac{u\left(1,y\right)-\exp\left(\frac{1}{2\epsilon_{2}}\right)u(\frac{1}{2},y)+\left[u(\frac{1}{2},y)-u\left(1,y\right)\right]\exp\left(\frac{x-\frac{1}{2}}{\epsilon_{2}}\right)}{1-\exp\left(\frac{1}{2\epsilon_{2}}\right)}$

Figure 2 shows the convergence rates that using both trial and test space
functions of the same polynomial degree deliver, for $p=1,2,3$. As expected,
we recover the same convergence rates as the original dG scheme: $h^{p}$ in
the $V_{h}$ error norm and $h^{p+1/2}$ in the $V_{h,\beta}$ error norm.

Figure 3: Anisotropic diffusion problem sketch. Counterclockwise advection
field.

###  4.4 Anisotropic with high-contrast diffusion

In this example, we consider an anisotropic problem, following [4]. We
consider the unit square $\Omega=[0,1]\times[0,1]$ split into four subdomains:
$\Omega_{1}=\left[0,\frac{2}{3}\right]\times\left[0,\frac{2}{3}\right]$,
$\Omega_{2}=\left[\frac{2}{3},1\right]\times\left[0,\frac{2}{3}\right]$,
$\Omega_{3}=\left[\frac{2}{3},1\right]\times\left[\frac{2}{3},1\right]$ and
$\Omega_{4}=\left[0,\frac{2}{3}\right]\times\left[\frac{2}{3},1\right]$. The
diffusivity tensor takes different values in each subdomain , see Figure 3:

$\displaystyle K_{i}(x,y)$ $\displaystyle=\begin{pmatrix}10^{-6}&0\\\
0&1.0\end{pmatrix}$ $\displaystyle\text{for }i=1,3\ \
\forall(x,y)\in\Omega_{1},\Omega_{3},$ $\displaystyle K_{i}(x,y)$
$\displaystyle=\begin{pmatrix}1.0&0\\\ 0&10^{-6}\end{pmatrix}$
$\displaystyle\text{for }i=2,4\ \ \forall(x,y)\in\Omega_{2},\Omega_{4}.$

(a) Counterclockwise advection field

(b) Clockwise advection field

Figure 4: Anisotropic diffusion problem. Uniform mesh (25.3 K DOFs).

(a) Computed solution and refined mesh, counterclockwise advection field
(Level 12: 23.7 K DOFs)

(b) Computed solution and refined mesh, clockwise advection field (Level 10:
22.3 K DOFs)

Figure 5: Anisotropic diffusion problem. Adaptive refinement.

The advection field is solenoidal ${\mathbf{b}=(\beta_{x},\beta_{y})^{T}}$,
with ${\beta_{x}=40x(2y-1)(x-1)}$ and ${\beta_{y}=-40y(2x-1)(y-1)}$ for the
counterclockwise case, and ${\mathbf{b}=-(\beta_{x},\beta_{y})^{T}}$ for the
clockwise case. Unlike the previous example, the advective field is neither
constant nor orthogonal to the discontinuities in $K$. However, its
orientation is still along the direction of increasing diffusivity; for that
reason, internal layers develop in the solution. The forcing term is
${f(x,y)=10^{-2}\exp(-(r-0.35)^{2}/0.005)}$ with
${r^{2}=(x-0.5)^{2}+(y-0.5)^{2}}$, corresponding to a Gaussian hill with
center at $r=0.35$. Finally, we set $\sigma=1$ for the reaction term, and
$g_{D}=0$ on $\partial\Omega$ for the boundary condition. We consider two
subcases, the first on a quasi-uniform mesh with ${h=0.024}$, conforming to
the discontinuities of $K$, and the second through an adaptive scheme starting
from a uniform triangular mesh (${h=0.177}$). We solve both cases using the
same polynomial degree ${p=1}$ for trial and test spaces. Figure 4 shows
results obtained for a uniform mesh, while Figure 5 shows results for adaptive
mesh refinement with fewer degrees of freedom than in the uniform mesh case.
We use the SWIP formulation such that the jump penalty considers the principal
directions of the diffusivity. Thus, SWIP avoids overshoots and undershoots
near the material interfaces, which is not the case for standard interior
penalty schemes [4]. Finally, the adaptive strategy concentrates refinement at
the inner layer but without losing the approximation quality.

(a) Level 0. Inner contour for $u=0.5$.

(b) Level 5. Inner contour for $u=1$.

(c) Level 9. Inner contour for $u=1$.

(d) Level 13. Inner contour for $u=1$.

Figure 6: 3D advection-dominated diffusion problem. Adaptive mesh evolution.

###  4.5 3D advection-dominated diffusion

As a last example, we consider a 3D advection-diffusion problem in the unit
cube ${\Omega=(0\,,\,1)^{3}\subset\mathbb{R}^{3}}$. We set the source term
${f=0}$, the diffusion ${K=10^{-3}}$, the spiral-type advective field
${\mathbf{b}=(\beta_{x},\beta_{y},\beta_{z})^{T}=(-0.15\,\sin(4\pi
z)\,,\,0.15\,\cos(4\pi z)\,,\,1)^{T}}$, and the inflow boundary datum $g$ as

$g=\left\\{\begin{array}[]{ll}1+\tanh\left[M\left(0.15^{2}-\left(x-0.6\right)^{2}-\left(y-0.5\right)^{2}\right)\right]&\textrm{on
}z=0,\\\ 0&\textrm{elsewhere on}\ \Gamma,\\\ \end{array}\right.$

These parameters produce a solution that presents an inner layer in the bottom
of the unit cube (${z=0}$) for values of $M\gg 1$ in the inflow boundary
datum. On the other hand, at the top of the unit cube (${z=1}$), the solution
exhibits a boundary layer due to the advection-dominant regime and the
boundary condition imposed on that surface. Figure 6 shows the evolution of
the 3D mesh as the refinement strategy progresses. Figure 6(a) displays the
initial progress of the spiral shape of the inner layer inside the unit cube.
Figures 6(b), 6(c), and 6(d) reflect the refinement process to capture the
sharp internal layer that the advective field induces as well as the boundary
layer that appears at the outlet due to the small diffusion. Figure 7, shoes
the discrete solution for a refined mesh displaying the interior mesh as well
as a cut of the solution.

Figure 7: 3D advection-dominated diffusion problem. Slice over $y=0.5$. Level
13: 4’858,125 DOFs.

##  5 Conclusions

In this paper, we describe an adaptive stabilized conforming finite element
method that minimizes the residual on dual discontinuous Galerkin (dG) norms
for advection-diffusion-reaction problems. We demonstrate the performance and
robustness of the method by analyzing in detail several challenging
diffusivity distributions. In particular, we tackle problems with highly
heterogeneous and anisotropic coefficient distributions as well as advection-
dominated cases. The method computes smooth solutions that minimize a residual
on a dual dG norm. This residual minimization problem leads to an equivalent
saddle-point problem, which we solve explicitly for a discretely continuous
approximation and an efficient and reliable error estimate that can guide
adaptive mesh refinement. In the context of advection-diffusion-reaction
problems, this method recovers the optimal convergence rates for $h$-adaptive
schemes in the dG norm. Besides, the method captures sharp boundary and
internal layers as well as overcomes the classical overshooting and
undershooting problems. In summary, we use several challenging numerical
examples to show the robustness of the method under a wide range of scenarios.

##  6 Acknowledgments

The CSIRO Deep Earth Imaging Future Science Platform Postgraduate Top-Up
Scholarship supports the work of RJC and this is gratefully acknowledged. This
publication was made possible in part by the CSIRO Professorial Chair in
Computational Geoscience at Curtin University and the Deep Earth Imaging
Enterprise Future Science Platforms of the Commonwealth Scientific Industrial
Research Organisation, CSIRO, of Australia. This project has also received
funding from the European Union’s Horizon 2020 research and innovation
programme under the Marie Sklodowska-Curie grant agreement No 777778
(MATHROCKS). At Curtin University, The Curtin Corrosion Centre, the Curtin
Institute for Computation, and The Institute for Geoscience Research (TIGeR)
kindly provide continuing support.

## References

  * Calo et al. [2020] V. M. Calo, A. Ern, I. Muga, S. Rojas,  An adaptive stabilized conforming finite element method via residual minimization on dual discontinuous Galerkin norms,  Computer Methods in Applied Mechanics and Engineering 363 (2020) 112891\.
  * Ewing and Wang [2001] R. E. Ewing, H. Wang,  A summary of numerical methods for time-dependent advection-dominated partial differential equations,  Journal of Computational and Applied Mathematics 128 (2001) 423–45.
  * Calo et al. [2014] V. M. Calo, Y. Efendiev, J. Galvis, M. Ghommem,  Multiscale empirical interpolation for solving nonlinear PDEs,  Journal of Computational Physics 278 (2014) 204 –20.
  * Ern et al. [2009] A. Ern, A. F. Stephansen, P. Zunino,  A discontinuous Galerkin method with weighted averages for advection–diffusion equations with locally small and anisotropic diffusivity,  IMA Journal of Numerical Analysis 29 (2009) 235–56.
  * Hossain et al. [2012] S. S. Hossain, S. F. A. Hossainy, Y. Bazilevs, V. M. Calo, T. J. R. Hughes,  Mathematical modeling of coupled drug and drug-encapsulated nanoparticle transport in patient-specific coronary artery walls,  Computational Mechanics 49 (2012) 213–42.
  * Calo et al. [2008] V. Calo, N. Brasher, Y. Bazilevs, T. Hughes,  Multiphysics model for blood flow and drug transport with application to patient-specific coronary artery flow,  Computational Mechanics 43 (2008) 161–77.
  * Bazilevs et al. [2007] Y. Bazilevs, V. M. Calo, T. E. Tezduyar, T. J. R. Hughes,  YZ$\beta$ discontinuity capturing for advection-dominated processes with application to arterial drug delivery,  International Journal for Numerical Methods in Fluids 54 (2007) 593–608.
  * Calo et al. [2016] V. M. Calo, Y. Efendiev, J. Galvis, G. Li,  Randomized oversampling for generalized multiscale finite element methods,  Multiscale Modeling & Simulation 14 (2016) 482–501.
  * Galvis et al. [2018] J. Galvis, E. Chung, Y. Efendiev, W. Leung,  On overlapping domain decomposition methods for high-contrast multiscale problems,  in: B. P. et al. (Ed.), Domain Decomposition Methods in Science and Engineering XXIV. DD 2017, volume 125 of Lecture Notes in Computational Science and Engineering, Springer, 2018\.
  * Calo et al. [2011] V. Calo, Y. Efendiev, J. Galvis,  A note on variational multiscale methods for high-contrast heterogeneous porous media flows with rough source terms,  Advances in Water Resources 34 (2011) 1177 –85. New Computational Methods and Software Tools.
  * Brooks and Hughes [1982] A. N. Brooks, T. J. 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 (1982) 199–259.
  * Hughes et al. [1989] T. J. Hughes, L. P. Franca, G. M. Hulbert,  A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations,  Computer Methods in Applied Mechanics and Engineering 73 (1989) 173–89.
  * Hughes and Sangalli [2007] T. J. Hughes, G. Sangalli,  Variational multiscale analysis: the fine-scale Green’s function, projection, optimization, localization, and stabilized methods,  SIAM Journal on Numerical Analysis 45 (2007) 539–57.
  * Brezzi et al. [1998] F. Brezzi, D. Marini, A. Russo,  Applications of the pseudo residual-free bubbles to the stabilization of convection-diffusion problems,  Computer Methods in Applied Mechanics and Engineering 166 (1998) 51–63.
  * Reed and Hill [1973] W. H. Reed, T. Hill, Triangular mesh methods for the neutron transport equation, Technical Report, Los Alamos Scientific Lab., N. Mex.(USA), 1973.
  * Lesaint and Raviart [1974] P. Lesaint, P.-A. Raviart,  On a finite element method for solving the neutron transport equation,  Publications mathématiques et informatique de Rennes (1974) 1–40.
  * Johnson and Pitkäranta [1986] C. Johnson, J. Pitkäranta,  An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation,  Mathematics of Computation 46 (1986) 1–26.
  * Cockburn et al. [2012] B. Cockburn, G. E. Karniadakis, C.-W. Shu, Discontinuous Galerkin methods: theory, computation and applications, volume 11, Springer Science & Business Media, 2012\.
  * Brezzi et al. [2004] F. Brezzi, L. D. Marini, E. Süli,  Discontinuous Galerkin methods for first-order hyperbolic problems,  Mathematical Models and Methods in Applied Sciences 14 (2004) 1893–903.
  * Ern and Guermond [2006] A. Ern, J.-L. Guermond,  Discontinuous Galerkin methods for Friedrichs’ systems. I. General theory,  SIAM Journal on Numerical Analysis 44 (2006) 753–78.
  * Bochev and Gunzburger [2009] P. B. Bochev, M. D. Gunzburger, Least-squares finite element methods, volume 166, Springer Science & Business Media, 2009.
  * Demkowicz and Gopalakrishnan [2010] L. Demkowicz, J. Gopalakrishnan,  A class of discontinuous Petrov–Galerkin methods. Part I: The transport equation,  Computer Methods in Applied Mechanics and Engineering 199 (2010) 1558–72.
  * Demkowicz and Gopalakrishnan [2011] L. Demkowicz, J. Gopalakrishnan,  A class of discontinuous Petrov–Galerkin methods. II. Optimal test functions,  Numerical Methods for Partial Differential Equations 27 (2011) 70–105.
  * Demkowicz and Gopalakrishnan [2014] L. Demkowicz, J. Gopalakrishnan,  An overview of the discontinuous Petrov Galerkin method,  in: X. Feng, O. Karakashian, Y. Xing (Eds.), Recent Developments in Discontinuous Galerkin Finite Element Methods for Partial Differential Equations: 2012 John H Barrett Memorial Lectures, volume 157 of The IMA Volumes in Mathematics and its Applications, Springer, Cham, 2014, pp. 149–80.
  * Calo et al. [2014] V. M. Calo, N. O. Collier, A. H. Niemi,  Analysis of the discontinuous Petrov–Galerkin method with optimal test functions for the Reissner–Mindlin plate bending model,  Computers & Mathematics with Applications 66 (2014) 2570 –86.
  * Demkowicz et al. [2012] L. Demkowicz, J. Gopalakrishnan, A. H. Niemi,  A class of discontinuous Petrov–Galerkin methods. Part III: Adaptivity,  Applied Numerical Mathematics 62 (2012) 396 – 427. Third Chilean Workshop on Numerical Analysis of Partial Differential Equations (WONAPDE 2010).
  * Niemi et al. [2011] A. Niemi, J. Bramwell, L. Demkowicz,  Discontinuous Petrov-Galerkin method with optimal test functions for thin-body problems in solid mechanics,  Computer Methods in Applied Mechanics and Engineering 200 (2011) 1291–300.
  * Niemi et al. [2013] A. H. Niemi, N. O. Collier, V. M. Calo,  Automatically stable discontinuous Petrov–Galerkin methods for stationary transport problems: Quasi-optimal test space norm,  Computers & Mathematics with Applications 66 (2013) 2096 –113. ICNC-FSKD 2012.
  * Niemi et al. [2011] A. H. Niemi, N. O. Collier, V. M. Calo,  Discontinuous Petrov-Galerkin method based on the optimal test space norm for one-dimensional transport problems,  Procedia Computer Science 4 (2011) 1862 –9. Proceedings of the International Conference on Computational Science, ICCS 2011.
  * Cier et al. [2020a] R. J. Cier, T. Poulet, S. Rojas, V. M. Calo, M. Veveakis,  Adaptive stabilized finite elements: Continuation analysis of compaction banding in geomaterials,  arXiv preprint arXiv:2008.01396 (2020a).
  * Cier et al. [2020b] R. J. Cier, S. Rojas, V. M. Calo,  A nonlinear weak constraint enforcement method for advection-dominated diffusion problems,  Mechanics Research Communications (2020b) 103602.
  * Łoś et al. [2020] M. Łoś, S. Rojas, M. Paszyński, I. Muga, V. M. Calo,  A stable discontinuous Galerkin based isogeometric residual minimization for the Stokes problem,  in: International Conference on Computational Science, Springer, 2020, pp. 197–211.
  * Kyburg et al. [2020] F. Kyburg, S. Rojas, V. M. Calo,  Incompressible flow modeling using an adaptive stabilized finite element method based on residual minimization,  arXiv preprint arXiv:2011.09182 (2020).
  * Rojas et al. [2020] S. Rojas, D. Pardo, P. Behnoudfar, V. M. Calo,  Residual minimization for goal-oriented adaptivity,  arXiv preprint arXiv:2007.08824 (2020).
  * Codina [1998] R. Codina,  Comparison of some finite element methods for solving the diffusion-convection-reaction equation,  Computer methods in applied mechanics and engineering 156 (1998) 185–210.
  * Hughes et al. [2018] T. J. Hughes, G. Scovazzi, L. P. Franca,  Multiscale and stabilized methods,  Encyclopedia of Computational Mechanics Second Edition (2018) 1–64.
  * Di Pietro et al. [2008] D. A. Di Pietro, A. Ern, J.-L. Guermond,  Discontinuous Galerkin methods for anisotropic semidefinite diffusion with advection,  SIAM Journal on Numerical Analysis 46 (2008) 805–31.
  * Di Pietro and Ern [2012] D. A. Di Pietro, A. Ern, Mathematical aspects of discontinuous Galerkin methods, volume 69, Springer Science, 2012\.
  * Shahbazi [2005] K. Shahbazi,  An explicit expression for the penalty parameter of the interior penalty method,  Journal of Computational Physics 205 (2005) 401–7.
  * Epshteyn and Rivière [2007] Y. Epshteyn, B. Rivière,  Estimation of penalty parameters for symmetric interior penalty Galerkin methods,  Journal of Computational and Applied Mathematics 206 (2007) 843–72.
  * Hartmann and Houston [2008] R. Hartmann, P. Houston,  An optimal order interior penalty discontinuous Galerkin discretization of the compressible Navier–Stokes equations,  Journal of Computational Physics 227 (2008) 9670–85.
  * Bastian et al. [2012] P. Bastian, M. Blatt, R. Scheichl,  Algebraic multigrid for discontinuous galerkin discretizations of heterogeneous elliptic problems,  Numerical Linear Algebra with Applications 19 (2012) 367–88.
  * Wheeler [1978] M. F. Wheeler,  An elliptic collocation-finite element method with interior penalties,  SIAM Journal on Numerical Analysis 15 (1978) 152–61.
  * Arnold [1982] D. N. Arnold,  An interior penalty finite element method with discontinuous elements,  SIAM Journal on Numerical Analysis 19 (1982) 742–60.
  * Arnold et al. [2002] D. N. Arnold, F. Brezzi, B. Cockburn, L. D. Marini,  Unified analysis of discontinuous Galerkin methods for elliptic problems,  SIAM Journal on Numerical Analysis 39 (2002) 1749–79.
  * Cohen et al. [2012] A. Cohen, W. Dahmen, G. Welper,  Adaptivity and variational stabilization for convection-diffusion equations,  ESAIM: Mathematical Modelling and Numerical Analysis 46 (2012) 1247–73.
  * Karakashian and Pascal [2003] O. A. Karakashian, F. Pascal,  A posteriori error estimates for a discontinuous galerkin approximation of second-order elliptic problems,  SIAM Journal on Numerical Analysis 41 (2003) 2374–99.
  * Burman and Ern [2007] E. Burman, A. Ern,  Continuous interior penalty hp-finite element methods for advection and advection-diffusion equations,  Mathematics of computation 76 (2007) 1119–40.
  * Ern and Guermond [2017] A. Ern, J.-L. Guermond,  Finite element quasi-interpolation and best approximation,  ESAIM: Mathematical Modelling and Numerical Analysis 51 (2017) 1367–85.
  * Alnæs et al. [2015] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, G. N. Wells,  The FEniCS project version 1.5,  Archive of Numerical Software 3 (2015) 9–23.
  * Dörfler [1996] W. Dörfler,  A convergent adaptive algorithm for Poisson’s equation,  SIAM Journal on Numerical Analysis 33 (1996) 1106–24.
  * Bank et al. [1983] R. E. Bank, A. H. Sherman, A. Weiser,  Some refinement algorithms and data structures for regular local mesh refinement,  Scientific Computing, Applications of Mathematics and Computing to the Physical Sciences 1 (1983) 3–17.
  * Bank et al. [1989] R. E. Bank, B. D. Welfert, H. Yserentant,  A class of iterative methods for solving saddle point problems,  Numerische Mathematik 56 (1989) 645–66.
  * Chen et al. [2008] Y. Chen, T. A. Davis, W. W. Hager, S. Rajamanickam,  Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate,  ACM Transactions on Mathematical Software (TOMS) 35 (2008) 22.
  * Mitchell [2013] W. F. Mitchell,  A collection of 2D elliptic problems for testing adaptive grid refinement algorithms,  Applied Mathematics and Computation 220 (2013) 350–64.
  * Oden and Patra [1995] J. T. Oden, A. Patra,  A parallel adaptive strategy for hp finite element computations,  Computer Methods in Applied Mechanics and Engineering 121 (1995) 449–70.
  * Burman and Zunino [2006] E. Burman, P. Zunino,  A domain decomposition method based on weighted interior penalties for advection-diffusion-reaction problems,  SIAM Journal on Numerical Analysis 44 (2006) 1612–38.