# On subgrid multiscale stabilized finite element method for advection-
diffusion-reaction equation with variable coefficients

Manisha Chowdhury , B.V. Rathish Kumar  Corresponding author,
Email addresses<EMAIL_ADDRESS>(M. Chowdhury) and
<EMAIL_ADDRESS>(B.V.R. Kumar)

(Indian Institute of Technology Kanpur
Kanpur, Uttar Pradesh, India)

###### Abstract

In this study a stabilized finite element method for solving advection-
diffusion-reaction equation with spatially variable coefficients has been
carried out. Here subgrid scale approach along with algebraic approximation to
the sub-scales has been chosen as stabilized method among various other
methods. Both a priori and a posteriori finite element error estimates in
$L_{2}$ norms have been derived after introducing the stabilized formulation
of the variational form. An expression of the stabilization parameter for this
2D problem has also been derived here. At last numerical experiments are
presented to verify numerical performance of the stabilized method and check
the credibility of the theoretically derived expression of the stabilization
parameter.

Keywords Advection-diffusion-reaction equation $\cdot$ Galerkin finite element
method $\cdot$ Subgrid scale methods $\cdot$ A priori error estimation $\cdot$
A posteriori error estimation

##  1 Introduction

Due to wide range of applications of the transport equation in the fields of
bio-medical engineering, chemical engineering, environmental sciences etc.,
solving more complex form of this equation has always been in great demand.
One of the appropriate mathematical model equations expressing transport
phenomena is the advection-diffusion-reaction (ADR) equation and its different
forms. Here particularly we dealing with ADR equation with spatially variable
diffusion coefficients and velocity components, due to its importance in
modelling several physical phenomena in more effective way. Since finiding
analytical solution of these type of equation is not an easy job, numerical
methods have been preferred in solving it. The article [5] has presented a
priori and a posteriori error bounds for the numerical scheme, Galerkin finite
element method. It is well-known that the Galerkin formulation suffers from
lack of stability when the diffusion coefficients are small compared to either
advection velocity components or reaction term. The purpose of this paper is
to present a study of a stabilized finite element method for this particular
type of ADR equation. Here we have considered to study the stabilization
method named by algebraic subgrid scale (ASGS) method i.e. subgrid scale (SGS)
method with an algebraic approximation to the subscales. [1] gives an essence
of stabilized multiscale subgrid method for Helmholtz equation, whereas [2]
presents SGS formulation of ADR equation with constant coefficients. We have
followed the same considerations to have SGS formulation with spatially
variable coefficients in section 2 of this paper. Before that we introduce
Galerkin finite element formulation for the equation in the same section. In
sections 3 and 4, a priori and a posteriori errors have been estimated for
stabilized form respectively. We have taken the conception of estimating
errors from [7] which provides a tutorial on the derivation of a priori and
explicit residual-based a posteriori error estimates for Galerkin finite
element discretization of general linear elliptic operators and [6] which
presents residual-based a posteriori bound for Galerkin least-squares (GLS)
formulation for Helmholtz equation. A subsection of the section 4 presents an
important highlight of this paper. It elaborately establishes the condition on
stabilization parameter for 2D ADR equation because study of stabilization
parameter only for 1D ADR equation is available in literatures [2],[3]. In the
last section numerical experiments have been presented to verify the
theoretically established results of error estimations as well as to check the
accuracy of the expression of the stabilization parameter for a test problem
whose conception is taken from literature in Hydrology [4]. This section also
presents a comparison of stability between the Galerkin method and the SGS
method for diffusion dominated and advection dominated flows and then
numerically establishes the fact that the Galerkin method suffers from lack of
stability.

##  2 Statement of the problem

Introducing the notations,
c: concentration of dispersing mass of the solute
$(D_{1},D_{2})$: diffusion coefficients along x-axis and y-axis
$(u_{1},u_{2})$: velocity components along x-axis and y-axis
$\mu$: coefficient of reaction
$q$: source of the solute mass in the flow domain

The two-dimensional ADR equation in general form with spatially variable
coefficient functions representing transportation of a solute in a bounded
domain $\Omega$ in $R^{2}$ along with homogeneous Dirichlet boundary condition
can be written as follows:
Find c : $\Omega\rightarrow$ R such that

$\frac{\partial}{\partial x}[D_{1}\frac{\partial c}{\partial
x}-u_{1}c]+\frac{\partial}{\partial y}[D_{2}\frac{\partial c}{\partial
y}-u_{2}c]-\mu c+q=0\hskip 5.69054ptin\hskip 2.84526pt\Omega$ (1)

c = 0 on $\partial\Omega$
Let us introduce two notations, $\widetilde{\bigtriangledown}$ :=
$(D_{1}\frac{\partial}{\partial x},D_{2}\frac{\partial}{\partial y})$ and
$\overline{u}=(u_{1},u_{2})$
Considering the fluid as incompressible i.e. the divergence of advection
velocity is zero, the above equation can be rewritten with the ADR operator
$\mathcal{L}$,

$\mathcal{L}=-\bigtriangledown\cdot\widetilde{\bigtriangledown}+\overline{u}\cdot\bigtriangledown+\mu\\\
$ (2)

whose adjoint operator $\mathcal{L^{*}}$ is

$\mathcal{L^{*}}=-\bigtriangledown\cdot\widetilde{\bigtriangledown}-\overline{u}\cdot\bigtriangledown+\mu$
(3)

The weak form of the boundary value problem is as follows: Find c $\in$ V such
that

$a(c,d)=l(d)\quad\forall d\in V$ (4)

where

a(c,d)=$\int_{\Omega}D_{1}\frac{\partial c}{\partial x}\frac{\partial
d}{\partial x}+\int_{\Omega}D_{2}\frac{\partial c}{\partial y}\frac{\partial
d}{\partial y}+\int_{\Omega}d$$\bar{u}$$\cdot\nabla c+\int_{\Omega}\mu cd$

l(d)=$\int_{\Omega}qd=(d,q)$
and let V represent both the usual spaces of trial solutions and weight
functions. In addition $(\cdot,\cdot):V\times V\rightarrow R$ represents
$L_{2}(\Omega)$ inner product and $\|\cdot\|$ denotes $L_{2}(\Omega)$ norm,
introduced by this inner product, i.e. $\|d\|^{2}=(d,d)=\int_{\Omega}d^{2}$.
The subscripts on the inner product and norm notations (whenever it will be
used) denote the domain of integration other than $\Omega$.

######  Remark 1.

$a(c,d)=(d,\mathcal{L}c)=(\mathcal{L^{*}}d,c)$

######  Remark 2.

Furthermore a(d,c)=$\int_{\Omega}D_{1}\frac{\partial d}{\partial
x}\frac{\partial c}{\partial x}+\int_{\Omega}D_{2}\frac{\partial d}{\partial
y}\frac{\partial c}{\partial y}-\int_{\Omega}c$$\bar{u}$$\cdot\nabla
d+\int_{\Omega}\mu cd$
Hence the bilinear form a(c,d) is not symmetric i.e. a(c,d) $\neq$ a(d,c).

###  2.1 Galerkin finite element formulation

Let the bounded domain $\Omega$ be discretized into finite number of element
subdomains $\Omega_{e}$, for e=1,2,…,$n_{el}$, where $n_{el}$ is the number of
elements. Let $h_{e}$ be the diameter of each element $\Omega_{e}$ and h
=$\underset{e=1,2,...,n_{el}}{max}$ $h_{e}$
Let $V_{h}$ $\subset$ V be a finite dimensional space containing continuous
piecewise polynomials of order p. Now the Galerkin formulation is : Find
$c_{h}\in V_{h}$ such that

$a(c_{h},d_{h})=l(d_{h})\quad\forall d_{h}\in V_{h}$ (5)

where

$a(c_{h},d_{h})$=$\int_{\Omega}D_{1}\frac{\partial c_{h}}{\partial
x}\frac{\partial d_{h}}{\partial x}+\int_{\Omega}D_{2}\frac{\partial
c_{h}}{\partial y}\frac{\partial d_{h}}{\partial
y}+\int_{\Omega}d_{h}$$\bar{u}$$\cdot\nabla c_{h}+\int_{\Omega}\mu c_{h}d_{h}$

$l(d_{h})$=$\int_{\Omega}qd_{h}$
Let $n_{pt}$ be the total number of nodes after discretizing the domain and
$S_{i}$ be the standard shape(basis) function at node each node a, for
a=1,2,…,$n_{pt}$ . Therefore the functions in $V_{h}$ can be interpolated from
these shape functions as,
$c_{h}$= $\sum_{a=1}^{n_{pt}}c_{a}S_{a}$ $d_{h}$ =
$\sum_{b=1}^{n_{pt}}d_{b}S_{b}$
Therefore

$\begin{split}a(c_{h},d_{h})&=\sum_{a,b=1}^{n_{pt}}d_{b}\hskip
2.84526pta(S_{a},S_{b})\hskip 2.84526ptc_{a}=D^{t}AC\\\
l(d_{h})&=\sum_{b=1}^{n_{pt}}d_{b}\hskip 2.84526ptl(S_{b})=D^{t}L\\\
\end{split}$ (6)

where A is a matrix of order $n_{pt}\times n_{pt}$ with components
$a(S_{a},S_{b})$ and L is a column vector with components $l(S_{b})$ for
a,b=1,2,…,$n_{pt}$
Hence the discrete problem is equivalent to solve as the following linear
system

$AC=L$ (7)

###  2.2 Stabilized integral form

Let $\Omega^{\prime}$ be the union of element interiors, that is
$\Omega^{\prime}=\cup_{e=1}^{n_{el}}\Omega_{k}$. The stabilized form of
equation (1) after applying SGS method,following [3], on it, is stated as
follows: Find $c_{h}\in V_{h}$ such that

$a_{SGS}(c_{h},d_{h})=l_{SGS}(d_{h})\quad\forall d_{h}\in V_{h}$ (8)

where

$a_{SGS}(c_{h},d_{h})$=
$a(c_{h},d_{h})+\int_{\Omega^{\prime}}(-\mathcal{L^{*}}d_{h})\tau\mathcal{L}c_{h}$

$l_{SGS}(d_{h})$= $l(d_{h})+\int_{\Omega^{\prime}}(-\mathcal{L^{*}}d_{h})\tau
q$
where $\tau$ is known as the stabilization parameter, which is measured later
in this paper.

##  3 A priori error estimation in the $L_{2}$ norm

A priori estimate provides a glimpse of convergence of the approximation
method though it cannot be computable as it depends upon the exact solution.
It also gives order of convergence before the solution is known. Before going
to estimate the error let us first introduce its splitting which plays an
important role in estimation. Let $c_{h}$ be the computed finite element
solution using stabilized SGS method and c be the exact solution. Due to
discretization of the domain we have an opportunity to introduce nodal
interpolants. Let $c_{h}^{\prime}$ be the nodal interpolant of c. Let e denote
the finite element error.
Now let us introduce the error splitting as follows:

$\begin{split}e&=c_{h}-c\\\ &=(c_{h}-c_{h}^{\prime})+(c_{h}^{\prime}-c)\\\
&=e_{h}+\lambda\end{split}$ (9)

where the first part $e_{h}=(c_{h}-c_{h}^{\prime})$ of the error belongs to
$V_{h}$ and the other part $\lambda=(c_{h}^{\prime}-c)$ is the interpolation
error belonging to V.

###  3.1 Nitsche trick

Generally Nitsche trick starts with considering an auxiliary problem. Here
this auxiliary problem involves the variational form of the ADR equation
expressed by adjoint operator $\mathcal{L^{*}}$, considering error e as the
source term. This is also known as dual problem. The reason behind replacing
the source term q by the error e is to introduce $L_{2}$ norm of e , which
will be seen in further proceedings. Now let us introduce the dual problem as
follows: Find $\eta\in V$ such that

$a(d,\eta)=(e,d)\quad\forall d\in V$ (10)

Now replacing d by e in the above equation, we will have

$a(e,\eta)=(e,e)=\|e\|^{2}$ (11)

The exact solution will also satisfy both the Galerkin finite element
formulation (5) and SGS formulation (8) and both will give the following
result:

$a(c,d_{h})=l(d_{h})\quad\forall d_{h}\in V_{h}$ (12)

Now subtracting (10) from (5), we will have

$a(e,d_{h})=0\quad\forall d_{h}\in V_{h}$ (13)

This is also known as Galerkin orthogonality.
Now we are going to bring SGS formulation into estimation by subtracting (10)
from (8) as follows:

$\begin{split}a(c_{h}-c,d_{h})+(-\mathcal{L^{*}}d_{h},\tau\mathcal{L}c_{h})&=(-\mathcal{L^{*}}d_{h},\tau
q)\quad\forall d_{h}\in V_{h}\\\
a(e,d_{h})+(-\mathcal{L^{*}}d_{h},\tau\mathcal{L}c_{h})&=(-\mathcal{L^{*}}d_{h},\tau\mathcal{L}c)\quad\forall
d_{h}\in V_{h}\\\
a(e,d_{h})+(-\mathcal{L^{*}}d_{h},\tau\mathcal{L}e)&=0\quad\forall d_{h}\in
V_{h}\\\ a_{SGS}(e,d_{h})&=0\quad\forall d_{h}\in V_{h}\\\ \end{split}$ (14)

Let $\eta_{h}$ be the interpolant of $\eta$ in $V_{h}$ and $\lambda_{\eta}$ =
$(\eta_{h}-\eta)$ be the interpolation error.
Since $\eta_{h}$ is in $V_{h}$, we can rewrite the last equation as follows,

$a(e,\eta_{h})+(-\mathcal{L^{*}}\eta_{h},\tau\mathcal{L}e)=0$ (15)

To estimate error in $L_{2}$ norm, now we are bringing it into the context by
subtracting (15) from (11). Hence,

$\begin{split}\|e\|^{2}&=a(e,\eta-\eta_{h})-(-\mathcal{L^{*}}\eta_{h},\tau\mathcal{L}e)\\\
&=(\mathcal{L^{*}}\eta_{h},\tau\mathcal{L}e)-a(e,\lambda_{\eta})\\\
&=(\mathcal{L^{*}}\eta_{h},\tau\mathcal{L}e_{h})+(\mathcal{L^{*}}\eta_{h},\tau\mathcal{L}\lambda)-a(e,\lambda_{\eta})\\\
\end{split}$ (16)

We have to find bounds for each term of the above equation so that one $\|e\|$
will be cancelled out from both sides. To have such bounds we need few
important results, which are given as follows:
Standard interpolation estimate [7] is,

$\|\lambda_{c}\|_{s}\leq\bar{C}(p,\Omega)h^{p+1-s}\|c\|_{p+1}$ (17)

The exact solution c has been assumed of regularity (p+1). In addition
$\|\cdot\|_{r}$ denotes standard $H^{r}(\Omega)$ norm (for r $\geq$ 0 and
integer), where $H^{r}(\Omega)$ is the standard Sobolev space of order r,2 on
$\Omega$. For r=0 this is $L_{2}$ norm. Similar to our previous consideration
the superscript on $\|\cdot\|_{r}$ denotes the domain of integration other
than $\Omega$.
The above result can be applied for the another interpolation error
$\lambda_{\eta}$ similarly, i.e.

$\|\lambda_{\eta}\|_{s}\leq\bar{C}(p,\Omega)h^{p+1-s}\|\eta\|_{p+1}$ (18)

Another results are strong stability conditions [5],

$\|c\|_{2}\leq C_{s}\|q\|$ (19)

and similarly for the dual problem

$\|\eta\|_{2}\leq C_{s}^{\prime}\|e\|$ (20)

Now looking at the RHS of the equation (16), we see that $a(e,\lambda_{\eta})$
is the Galerkin part, which can be bounded by the result obtained in [5]. But
it is not possible to find the desired bound for the first term due to
presence of $e_{h}$ . Therefore to get rid of $e_{h}$ we are doing the
following steps.
Replacing $c_{h}$ by $c_{h}^{\prime}+e_{h}$ in SGS formulation (6) we have

$\begin{split}a_{SGS}(c_{h}^{\prime}+e_{h},d_{h})&=l_{SGS}(d_{h})\quad\forall
d_{h}\in V_{h}\\\
a_{SGS}(c_{h}^{\prime},d_{h})+a_{SGS}(e_{h},d_{h})&=l_{SGS}(d_{h})\quad\forall
d_{h}\in V_{h}\\\
a_{SGS}(c_{h}^{\prime},d_{h})+a_{SGS}(e,d_{h})-a_{SGS}(\lambda,d_{h})&=l_{SGS}(d_{h})\quad\forall
d_{h}\in V_{h}\\\
a_{SGS}(c_{h}^{\prime},d_{h})-a_{SGS}(\lambda,d_{h})&=l_{SGS}(d_{h})\quad\forall
d_{h}\in V_{h}\\\ \end{split}$ (21)

Subtracting (21) from (8), we will have

$\begin{split}a_{SGS}(e_{h},d_{h})+a_{SGS}(\lambda,d_{h})&=0\quad\forall
d_{h}\in V_{h}\\\
a(e_{h},d_{h})+(-\mathcal{L^{*}}d_{h},\tau\mathcal{L}e_{h})+a(\lambda,d_{h})+(-\mathcal{L^{*}}d_{h},\tau\mathcal{L}\lambda)&=0\quad\forall
d_{h}\in V_{h}\\\
(\mathcal{L^{*}}d_{h},\tau\mathcal{L}e_{h})+(\mathcal{L^{*}}d_{h},\tau\mathcal{L}\lambda)&=0\quad\forall
d_{h}\in V_{h}\\\ \end{split}$ (22)

Since $\eta_{h}$ is in $V_{h}$ therefore the above equation can be rewritten
as:

$(\mathcal{L^{*}}\eta_{h},\tau\mathcal{L}e_{h})+(\mathcal{L^{*}}\eta_{h},\tau\mathcal{L}\lambda)=0\\\
$ (23)

Applying the result obtained above into equation (16), we will remain with the
Galerkin part only. As mentioned above, we now use the result from [5]
obtained to bound the Galerkin part using the standard interpolation and
stability estimates. Therefore we finally have,

$\begin{split}\|e\|^{2}&\leq\mid a(e,\lambda_{\eta})\mid\\\ &\leq
Ch^{p+1}\|e\|\|c\|_{p+1}\end{split}$ (24)

Hence,

$\boxed{\|e\|\leq Ch^{p+1}\|c\|_{p+1}}$ (25)

######  Remark 3.

If we choose the finite element space as P1, the space containing continuous
piecewise polynomial of order 1, a priori estimate using the stability
estimate (19) will be as follows:

$\begin{split}\|e\|&\leq Ch^{2}\|c\|_{2}\\\ &\leq CC_{s}h^{2}\|q\|\end{split}$
(26)

##  4 A posteriori error estimation in $L_{2}$ norm

This estimation is carried out after the solution is known and it depends upon
that computed solution, therefore it is computable. It ensures the order of
convergence obtained in a priori estimate. Here we are going to find residual
based a posteriori estimate. Residual denoted by $r_{h}$ is the term
$(\mathcal{L}c_{h}-q)$. We will start with the same Nitsche trick i.e.
considering the dual problem to proceed further. The steps will be slightly
different from what have been done in previous section because we have to get
rid of the exact solution c and bring residual into the picture. In this
section again we have to deal with the Galerkin part, bound for which will be
taken from [5].
Now let us again introduce the dual problem as follows: Find $\eta\in V$ such
that

$a(d,\eta)=(d,e)\quad\forall d\in V$ (27)

According to the definition of $L_{2}$ norm,

$\begin{split}\|e\|^{2}&=(e,e)\\\ &=a(e,\eta)\\\ &=a(c_{h}-c,\eta)\\\
&=a(c_{h},\eta)-a(c,\eta)\\\ &=a(c_{h},\eta)-l(\eta)\\\ \end{split}$ (28)

Using weak form (4) in the last line by the fact that $\eta\in V$ we get rid
of the exact solution.
The SGS stabilized form (8) can be rewritten as follows: Find $c_{h}\in V_{h}$
such that

$\begin{split}a(c_{h},d_{h})+(-\mathcal{L}^{*}d_{h},\tau((\mathcal{L}c_{h}-q))_{\Omega^{\prime}}&=l(d_{h})\quad\forall
d_{h}\in V_{h}\\\ a(c_{h},d_{h})+(-\mathcal{L}^{*}d_{h},\tau
r_{h})_{\Omega^{\prime}}&=l(d_{h})\quad\forall d_{h}\in V_{h}\\\ \end{split}$
(29)

Let $\eta_{h}$ be the interpolant of $\eta$ in $V_{h}$. Hence the above
equation can be rewritten as:

$a(c_{h},\eta_{h})+(-\mathcal{L}^{*}\eta_{h},\tau
r_{h})_{\Omega^{\prime}}-l(\eta_{h})=0$ (30)

Here we follow the same notation $\lambda_{\eta}$ for the interpolation error
$(\eta_{h}-\eta)$.
Before proceeding further let us introduce the following appropriate
interpolation estimate [8] for each element $\Omega_{k}$, k=1,2,…,$n_{el}$

$\begin{split}\|\lambda_{\eta}\|_{s,\Omega_{k}}&\leq\bar{C}(\Omega_{k})\hskip
2.84526pth_{k}^{2-s}\hskip 2.84526pt\|\eta\|_{2}\\\
&\leq\bar{C}(\Omega_{k})\hskip 2.84526ptC_{s}^{\prime}\hskip
2.84526pth_{k}^{2-s}\hskip 2.84526pt\|e\|\end{split}$ (31)

where s is an integer and the second step comes from using the strong
stability condition (20).
Now we are coming back in deducing the estimation by subtracting (30) from
(28),

$\begin{split}\|e\|^{2}&=a(c_{h},\eta-\eta_{h})-(-\mathcal{L}^{*}\eta_{h},\tau
r_{h})_{\Omega^{\prime}}-l(\eta-\eta_{h})\\\ &=(\mathcal{L}^{*}\eta_{h},\tau
r_{h})_{\Omega^{\prime}}+l(\eta_{h}-\eta)-a(c_{h},\eta_{h}-\eta)\\\
&=(\mathcal{L}^{*}\eta_{h},\tau
r_{h})_{\Omega^{\prime}}+l(\lambda_{\eta})-a(c_{h},\lambda_{\eta})\\\
&\leq\mid(\mathcal{L}^{*}\eta_{h},\tau r_{h})_{\Omega^{\prime}}\mid+\mid
a(c_{h},\lambda_{\eta})-l(\lambda_{\eta})\mid\\\ \end{split}$ (32)

Here the second term is the Galerkin part. [5] gives its bound using the above
interpolation estimate as follows:

$\mid
a(c_{h},\lambda_{\eta})-l(\lambda_{\eta})\mid\leq(\sum_{k=1}^{n_{el}}\bar{C}(\Omega_{k})\hskip
2.84526ptC_{s}^{\prime}\hskip 2.84526pth_{k}^{2}\hskip
2.84526pt\|r_{h}\|_{\Omega_{k}})\hskip 2.84526pt\|e\|$ (33)

Now our job is to find bound for the first term.

$\begin{split}\mid(\mathcal{L}^{*}\eta_{h},\tau
r_{h})_{\Omega^{\prime}}\mid&=\mid(\mathcal{L}^{*}\lambda_{\eta},\tau
r_{h})_{\Omega^{\prime}}+(\mathcal{L}^{*}\eta,\tau
r_{h})_{\Omega^{\prime}}\mid\\\ &\leq\mid(\mathcal{L}^{*}\lambda_{\eta},\tau
r_{h})_{\Omega^{\prime}}\mid+\mid(e,\tau r_{h})_{\Omega^{\prime}}\mid\\\
&=\mid\tau\mid\hskip
2.84526pt\mid\sum_{k=1}^{n_{el}}(\mathcal{L}^{*}\lambda_{\eta},r_{h})_{\Omega_{k}}\mid+\mid\tau\mid\hskip
2.84526pt\mid\sum_{k=1}^{n_{el}}(e,r_{h})_{\Omega_{k}}\mid\\\
&\leq\mid\tau\mid(\sum_{k=1}^{n_{el}}\|\mathcal{L}^{*}\lambda_{\eta}\|_{\Omega_{k}}\|r_{h}\|_{\Omega_{k}}+\sum_{k=1}^{n_{el}}\|e\|_{\Omega_{k}}\|r_{h}\|_{\Omega_{k}})\\\
&\leq\mid\tau\mid(\sum_{k=1}^{n_{el}}\|\mathcal{L}^{*}\lambda_{\eta}\|_{\Omega_{k}}\|r_{h}\|_{\Omega_{k}}+\sum_{k=1}^{n_{el}}\|r_{h}\|_{\Omega_{k}}\|e\|)\end{split}$
(34)

Now we are only remained with estimating
$\|\mathcal{L}^{*}\lambda_{\eta}\|_{\Omega_{k}}$ as follows:

$\begin{split}\|\mathcal{L}^{*}\lambda_{\eta}\|_{\Omega_{k}}&=\|-\bigtriangledown\cdot\widetilde{\bigtriangledown}\lambda_{\eta}-\overline{u}\cdot\bigtriangledown\lambda_{\eta}+\mu\lambda_{\eta}\|_{\Omega_{k}}\\\
&=\|-D_{1}\frac{\partial^{2}\lambda_{\eta}}{\partial
x^{2}}-D_{2}\frac{\partial^{2}\lambda_{\eta}}{\partial
y^{2}}+(u_{1}-\frac{\partial D_{x}}{\partial
x})\frac{\partial\lambda_{\eta}}{\partial x}+(u_{2}-\frac{\partial
D_{y}}{\partial y})\frac{\partial\lambda_{\eta}}{\partial
y}+\mu\lambda_{\eta}\|_{\Omega_{k}}\\\ &\leq D\hskip
2.84526pt(\|\frac{\partial^{2}\lambda_{\eta}}{\partial
x^{2}}\|_{\Omega_{k}}+\|\frac{\partial^{2}\lambda_{\eta}}{\partial
y^{2}}\|_{\Omega_{k}}+\|\frac{\partial\lambda_{\eta}}{\partial
x}\|_{\Omega_{k}}+\|\frac{\partial\lambda_{\eta}}{\partial
y}\|_{\Omega_{k}}+\|\lambda_{\eta}\|_{\Omega_{k}})\\\ &=D\hskip
2.84526pt\|\lambda_{\eta}\|_{2,{\Omega_{k}}}\quad(by\hskip
2.84526ptstandard\hskip 2.84526ptdefinition\hskip 2.84526ptof\hskip
2.84526ptH^{2}\hskip 2.84526ptnorm)\\\ &\leq D\hskip
2.84526pt\bar{C}(\Omega_{k})\hskip 2.84526ptC_{s}^{\prime}\hskip
2.84526pt\|e\|\end{split}$ (35)

where D is modulus of the maximum taken over all the coefficients present in
that expression and we have obtained the last line using interpolation
estimate (31).
Putting the result obtained in (35) into (34), we will have,

$\begin{split}\mid(\mathcal{L}^{*}\eta_{h},\tau
r_{h})_{\Omega^{\prime}}\mid&\leq(\sum_{k=1}^{n_{el}}(D\hskip
2.84526pt\bar{C}(\Omega_{k})\hskip 2.84526ptC_{s}^{\prime}+1)\hskip
2.84526pt\|r_{h}\|_{\Omega_{k}})\hskip 2.84526pt\mid\tau\mid\|e\|\\\
&=\sum_{k=1}^{n_{el}}\bar{D}_{k}\mid\tau\mid\|r_{h}\|_{\Omega_{k}}\hskip
2.84526pt\|e\|\\\ \end{split}$ (36)

where $\bar{D}_{k}=(D\hskip 2.84526pt\bar{C}(\Omega_{k})\hskip
2.84526ptC_{s}^{\prime}+1)$
Combining both the results obtained in (33) and (36) and substituting them
into equation (32), we will remain with the following:

$\|e\|^{2}\leq(\sum_{k=1}^{n_{el}}(\bar{C}(\Omega_{k})\hskip
2.84526ptC_{s}^{\prime}\hskip 2.84526pth_{k}^{2}\hskip
2.84526pt+\bar{D}_{k}\mid\tau\mid)\hskip
2.84526pt\|r_{h}\|_{\Omega_{k}})\hskip 2.84526pt\|e\|$ (37)

Hence,

$\boxed{\|e\|\leq(\sum_{k=1}^{n_{el}}(\bar{C}(\Omega_{k})\hskip
2.84526ptC_{s}^{\prime}\hskip 2.84526pth_{k}^{2}\hskip
2.84526pt+\bar{D}_{k}\mid\tau\mid)\hskip 2.84526pt\|r_{h}\|_{\Omega_{k}})}$
(38)

######  Remark 4.

The above a posteriori estimate is computable with an appropriate estimate of
the constants $\bar{C}(\Omega_{k}),C_{s}^{\prime},\bar{D}_{k}$ and the unknown
stabilization parameter $\tau$. An estimation of $\tau$ is derived in the next
section.

###  4.1 Estimating $\tau$

In [2] and [3], Codina gives an approach based on discrete maximum principle
(DMP) to find conditions on $\tau$ for scalar 1D case. We are going to apply
the same idea to reach to a condition on $\tau$ for general 2D case and then
use that for 2D ADR equation with spatially variable coefficients. The idea
has been taken up from the following implications:
The finite element matrix A in non-negative $\Longrightarrow$ the DMP holds
$\Longrightarrow$ stability and convergence of the finite element solution.
In both the references mentioned in the above paragraph, a condition on $\tau$
has been obtained by showing the matrix of the final algebraic system to be
non-negative. Before introducing the definition of non-negativity of matrix A,
let us mention some useful notions. We already have $n_{pt}$ as the total
number of nodes of finite element mesh. Let $n_{free}$ be the interior nodes.
The finite element discretization of the problem leads to (7), where the
values $c_{a}$ for a= $n_{free}$+1,…,$n_{pt}$ i.e. for the boundary nodes, are
known from the Dirichlet boundary conditions. Now the matrix A have dimensions
$n_{free}$ $\times$ $n_{pt}$ and the column vector on the RHS will have
$n_{free}$ components.

######  Definition 1.

Matrix A= $(a_{ij})_{n_{free}\times n_{pt}}$ is of non-negative type if the
following conditions hold:
(i) $a_{ij}\leq 0$ $\forall$ $i=1,2,...,n_{free}$ , $j=1,2,...,n_{pt}$ ,$i\neq
j$
(j) $\sum_{j=1}^{n_{pt}}a_{ij}\geq 0$ for $i=1,2,...,n_{free}$

To find a condition on $\tau$ requiring the finite element matrix to be non-
negative we first take a sup norm of the spatially variable coefficient
functions because if the condition holds for the supremum of the coefficients,
it will hold for every possible values of the coefficients as well as we can
deal with the variational form of constant coefficients. The numerical
experiments will explain this better later.
Let D and U be the sup norm values taken over the diffusion coefficients and
velocity components respectively on $\Omega$.
Hence the variational form (4) will look like as follows: Find c $\in$ V such
that

$B(c,d)=F(d)\forall d\in V$ (39)

where

B(c,d)=$D\int_{\Omega}(\frac{\partial c}{\partial x}\frac{\partial d}{\partial
x}+\frac{\partial c}{\partial y}\frac{\partial d}{\partial
y})+U\int_{\Omega}d(\frac{\partial c}{\partial x}+\frac{\partial c}{\partial
y})+\mu\int_{\Omega}cd$

l(d)=$\int_{\Omega}qd=(d,q)$
Since the assembly operator is summation of the element contributions,
therefore it is sufficient to check the element matrices to be non-negative.
Hence we will now try to find the condition for one element $\Omega_{e}$.
Let us discretize the domain into triangular elements and $\Lambda_{e}$ be the
area of each element $\Omega_{e}$ for e=1,2,…,$n_{el}$. We have taken
continuous piecewise linear polynomials as the shape functions for each
element. Letting $\\{N_{1}^{e},N_{2}^{e},N_{3}^{e}\\}$ be the set of basis
functions for each element $\Omega_{e}$, the element matrices will be as
follows:

$\begin{bmatrix}B^{e}(N_{1}^{e},N_{1}^{e})&B^{e}(N_{1}^{e},N_{2}^{e})&B^{e}(N_{1}^{e},N_{3}^{e})\\\
B^{e}(N_{2}^{e},N_{1}^{e})&B^{e}(N_{2}^{e},N_{2}^{e})&B^{e}(N_{2}^{e},N_{3}^{e})\\\
B^{e}(N_{3}^{e},N_{1}^{e})&B^{e}(N_{3}^{e},N_{2}^{e})&B^{e}(N_{3}^{e},N_{3}^{e})\\\
\end{bmatrix}$

where $B^{e}(\cdot,\cdot)$ is bilinear form defined on each element
$\Omega_{e}$
$N_{1}^{e}(x,y)=\frac{1}{2\Lambda_{e}}\\{(x_{2}^{e}y_{3}^{e}-x_{3}^{e}y_{2}^{e})+(y_{2}^{e}-y_{3}^{e})x+(x_{3}^{e}-x_{2}^{e})y\\}$
$N_{2}^{e}(x,y)=\frac{1}{2\Lambda_{e}}\\{(x_{3}^{e}y_{1}^{e}-x_{1}^{e}y_{3}^{e})+(y_{3}^{e}-y_{1}^{e})x+(x_{1}^{e}-x_{3}^{e})y\\}$
$N_{3}^{e}(x,y)=\frac{1}{2\Lambda_{e}}\\{(x_{1}^{e}y_{2}^{e}-x_{2}^{e}y_{1}^{e})+(y_{1}^{e}-y_{2}^{e})x+(x_{2}^{e}-x_{1}^{e})y\\}$
where $(x_{1}^{e},y_{1}^{e}),(x_{2}^{e},y_{2}^{e})$ and
$(x_{3}^{e},y_{3}^{e})$ are three vertices of element $\Omega_{e}$.

Introducing the notations $\beta_{i}^{e}=(y_{j}^{e}-y_{k}^{e})$ and
$\gamma^{e}=-(x_{j}^{e}-x_{k}^{e})$, the orientation of vertices can be
considered such a way in anti-clock wise direction that the following form
will hold for each element.

$B^{e}(N_{i}^{e},N_{j}^{e})=-\frac{D}{4\Lambda_{e}}K_{ij}^{e}\pm\frac{U}{6}E_{i}^{e}+\frac{\mu\Lambda_{e}}{12}$
for i,j=1,2,3.
here $K_{ij}^{e}=(\beta_{i}^{e}\beta_{j}^{e}+\gamma_{i}^{e}\gamma_{j}^{e})$
and $E_{i}^{e}=(\beta_{i}^{e}+\gamma_{i}^{e})$ both are always positive
constants.
According to the definition (1), non-negativity of each element matrix is
equivalent to

$-\frac{D}{4\Lambda_{e}}K_{ij}^{e}\pm\frac{U}{6}E_{i}^{e}+\frac{\mu\Lambda_{e}}{12}\leq
0\quad\forall i,j=1,2,3\quad and\quad i\neq j$ (40)

When SGS method is applied, the finite element matrices will have the same
form (40), but with modified coefficients such as:
$\bar{D}=D+\tau U^{2}\quad\bar{U}=U-2\mu\tau U\quad\bar{\mu}=\mu-\tau\mu^{2}$
Hence the above condition (40) will be modified as follows:

$(-\frac{D}{4\Lambda_{e}}K_{ij}^{e}\pm\frac{U}{6}E_{i}^{e}+\frac{\mu\Lambda_{e}}{12})-\tau(-\frac{U^{2}}{4\Lambda_{e}}K_{ij}^{e}\pm\frac{2\mu
U}{6}E_{i}^{e}+\frac{\mu^{2}\Lambda_{e}}{12})\leq 0\quad\forall i,j=1,2,3\quad
and\quad i\neq j.$ (41)

We will reach at a condition on $\tau$ by considering two simplified cases
following [3], as follows:
Case 1
Let us consider the reaction term $\mu$ = 0. Then (41) becomes,

$\begin{split}\tau&\geq(-\frac{D}{U^{2}}\pm\frac{2\Lambda_{e}E_{i}^{e}}{3UK_{ij}^{e}})\\\
&\geq(\frac{2\Lambda_{e}E_{i}^{e}}{3UK_{ij}^{e}}-\frac{D}{U^{2}})\\\
&\geq(\frac{2h_{e}}{3U}-\frac{D}{U^{2}})\\\
&=\frac{2h}{3U}(1-\frac{1}{Q})\end{split}$ (42)

where we take Q = $\frac{2Uh_{e}}{3D}$ and h satisfies the condition
$\frac{\Lambda_{e}E_{i}^{e}}{K_{ij}^{e}}\geq h_{e}$
Case 2
This time let us consider the advection velocity U to be 0. Then the condition
(41) reduces to,

$\begin{split}\tau&\geq(\frac{1}{\mu}-\frac{3DK_{ij}^{e}}{\mu^{2}\Lambda_{e}^{2}})\\\
&\geq(\frac{1}{\mu}-\frac{3DK_{ij}^{e}}{\mu^{2}h_{e}^{2}})\\\
&=\frac{1}{\mu}(1-\frac{3}{R})\end{split}$ (43)

where $\frac{\mu h^{2}}{D}$ is taken as R along with a similar consideration
as above, $\frac{\Lambda_{e}^{2}}{K_{ij}^{e}}\geq h_{e}^{2}$
Since,
$\frac{2h_{e}}{3U}(1-\frac{1}{Q})\leq\frac{2h_{e}}{3U}(\frac{Q}{1+Q})=(\frac{9D}{4h_{e}^{2}}+\frac{3U}{2h_{e}})^{-1}\leq(\frac{9D}{4h^{2}}+\frac{3U}{2h})^{-1}$
And,
$\frac{1}{\mu}(1-\frac{3}{R})\leq\frac{1}{\mu}(\frac{4R}{4R+9})=(\frac{9D}{4h_{e}^{2}}+\mu)^{-1}\leq(\frac{9D}{4h^{2}}+\mu)^{-1}$
Now following the approach in [3], $\tau$ can be chosen as,

$\boxed{\tau=(\frac{9D}{4h^{2}}+\frac{3U}{2h}+\mu)^{-1}}$ (44)

######  Remark 5.

From the above expression of $\tau$ we can conclude that $\tau$ is of order of
$h^{2}$.

##  5 Numerical Experiment

As mentioned earlier, Galerkin finite element method fails when diffusion
coefficients are small compared to either advection or reaction coefficients
or both. In this section we will prove this fact as well as verify the
theoretically established result through numerical data obtained by solving a
test problem with spatially variable coefficients of hydrological
significance. Here we have chosen two cases. Whereas the first case will
analyse the comparison between the Galerkin and the stabilized SGS method for
diffusion dominated flow, the second one will do the same for advection-
reaction dominated flow.
As a test problem we have considered a source of pollutant dispersed into an
incompressible fluid in a simple bounded square domain
$\Omega=(0,1)\times(0,1)$. This phenomena respects ADR equation with spatially
variable coefficients along with homogeneous Dirichlet boundary condition. The
concepts of choosing the coefficients are taken from [4], which is motivated
by the fact of real-life application from hydrology. In that study, velocity
components are taken spatial linear non-homogeneous functions, whereas the
diffusion coefficients are considered proportional to the square of respective
velocity components.

###  5.1 First case: Diffusion dominated flow

In this case we will see that the Galerkin finite element method and the
stabilized SGS method perform equally well in establishing convergence. Let us
introduce the coefficients we have worked with as follows:
$u_{1}=0.5(1+0.02x)$, $u_{2}=-0.5(1+0.02y)$, $D_{1}=(1+0.02x)^{2}$,
$D_{2}=0.1(1+0.02y)^{2}$ and $\mu=0.01$
Then D=1.0404 and U= 0.51.
The exact solution is c= sin(xy(x-1)(y-1)) on $\Omega$ along with c vanishing
on the boundary $\partial\Omega$. Working with the same expression for $\tau$
as deduced in the previous section, we have obtained results using freefem++.

######  Remark 6.

The error and order of convergence obtained for the Galerkin and the SGS
method for diffusion dominated flow have been shown in table 1 and table 2
respectively.

######  Remark 7.

The results of the table 2 shows that the order of convergence for the
stabilized SGS method is 2 which also justifies the deduced expression the
stabilization parameter $\tau$ numerically.

Mesh size | Error in $L_{2}$ norm | Order of convergence
---|---|---
10 | 0.000254151 |
20 | 6.28651 $e^{-5}$ | 2.01536
40 | 1.26455 $e^{-5}$ | 2.31363
80 | 3.01483 $e^{-6}$ | 2.06848
160 | 8.30888 $e^{-7}$ | 1.85935
320 | 1.72972 $e^{-7}$ | 2.26412
Table 1: Error and Order of convergence obtained in $L_{2}$ norm by Galerkin method for diffusion dominated flow Mesh size | Error in $L_{2}$ norm | Order of convergence
---|---|---
10 | 0.000405122 |
20 | 0.000105007 | 1.94787
40 | 2.40648 $e^{-5}$ | 2.12549
80 | 5.89785 $e^{-6}$ | 2.02866
160 | 1.55175 $e^{-6}$ | 1.92629
320 | 3.49761 $e^{-7}$ | 2.14945
Table 2: Error and Order of convergence obtained in $L_{2}$ norm by the SGS method for diffusion dominated flow Mesh size | Error in $L_{2}$ norm | Order of convergence
---|---|---
10 | 0.000541066 |
20 | 0.000375131 | 0.528409
40 | 9.51933 $e^{-5}$ | 1.97846
80 | 3.399958 $e^{-5}$ | 1.4855
160 | 1.25028 $e^{-5}$ | 1.4431
320 | 4.00531 $e^{-6}$ | 1.64227
Table 3: Error and Order of convergence obtained in $L_{2}$ norm by Galerkin method for flow having small diffusion coefficients Mesh size | Error in $L_{2}$ norm | Order of convergence
---|---|---
10 | 0.000280761 |
20 | 6.90662 $e^{-5}$ | 2.02329
40 | 1.73534 $e^{-5}$ | 1.99276
80 | 4.64904 $e^{-6}$ | 1.90022
160 | 1.37837 $e^{-6}$ | 1.85397
320 | 3.50102 $e^{-7}$ | 1.97711
Table 4: Error and Order of convergence obtained in $L_{2}$ norm by the SGS
method for flow having small diffusion coefficients

###  5.2 Second case: Advection-reaction dominated flow

We have chosen diffusion very small compared to advection and reaction
coefficients. The coefficients following the same pattern as above are,
$u_{1}=(1+0.02x)$, $u_{2}=-(1+0.02y)$, $D_{1}=10^{-7}(1+0.02x)^{2}$,
$D_{2}=10^{-8}(1+0.02y)^{2}$ and $\mu=1$
Then D= 1.0404 $\times 10^{-7}$ and U=1.02
The exact solution will be the same.

######  Remark 8.

Table 3 and table 4 represent the error and order of convergence obtained for
Galerkin and the SGS method respectively for very small diffusion
coefficients. Table 3 clearly shows oscillations in order of convergence for
Galerkin method whereas The order of convergence under the SGS method is
converging to 2.

## References

  * [1] T.J.R. Hughes, Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1995) 387-401.
  * [2] R. Codina, Comparison of some finite element methods for solving the diffusion-convection-reaction equation, Comput. Methods Appl. Mech. Engrg. 156 (1998) 185-210.
  * [3] R. Codina, On stabilized finite element methods for linear systems of convection-diffusion-reaction equations, Comput. Methods Appl. Mech. Engrg. 188 (2000) 61-82.
  * [4] A.Sanskrityayn, V.P.Singh, V.K.Bharati and N.Kumar, Analytical solution of two-dimensional advection-dispersion equation with spatio-temporal coefficients for point sources in an infinite medium using Green’s function method, Environ. Fluid Mech. 18(2018) 739-757.
  * [5] B.V.R. Kumar, M. Chowdhury, A priori and a posteriori error estimation for finite element approximation of advection-diffusion-reaction equation with spatially variable coefficients.(pre-prints), arXiv: 1811.05283v1 [math.AP] 13 Nov. 2018.
  * [6] J.R. Stewart, T.J.R. Hughes, An a posteriori error estimator and hp-adaptive strategy for finite element discretizations of the Helmholtz equation in exterior domains, Finite Elements in Analysis and Design 25 (1997) 1-26.
  * [7] J.R.Stewart and T.J.R. Hughes, A tutorial in elementary finite element error analysis: A systematic presentation of a priori and a posteriori error estimates, Comput. Methods Appl. Mech. Engrg. 158(1998), 1-22.
  * [8] C.Johnson and P.Hansbo, Adaptive finite element methods in computational mechanics, Comput. Methods Appl. Mech. Engrg., 101(1992), 143-181.