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

TL;DR

This paper develops a stabilized finite element method using subgrid scales for solving advection-diffusion-reaction equations with variable coefficients, providing error estimates and a stabilization parameter expression verified through numerical experiments.

Contribution

Introduces a subgrid multiscale stabilized finite element method with algebraic sub-scale approximation for variable coefficient advection-diffusion-reaction equations, including error analysis and stabilization parameter derivation.

Findings

Error estimates in $L_2$ norms confirm method accuracy.

Derived stabilization parameter improves numerical stability.

Numerical experiments validate theoretical results.

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.

Tables

Table 1. Table 1: Error and Order of convergence obtained in L 2 subscript 𝐿 2 L_{2} norm by Galerkin method for diffusion dominated flow

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 2. Table 2: Error and Order of convergence obtained in L 2 subscript 𝐿 2 L_{2} norm by the SGS 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 3. Table 3: Error and Order of convergence obtained in L 2 subscript 𝐿 2 L_{2} norm by Galerkin method for flow having small diffusion coefficients

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 4. Table 4: Error and Order of convergence obtained in L 2 subscript 𝐿 2 L_{2} norm by the SGS 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

Equations

$$\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=0in\Omega$$

$$\mathcal{L}=-▽\cdot \widetilde{▽}+\overline{u}\cdot ▽+\mu$$

$${\mathcal{L}}^{\ast }=-▽\cdot \widetilde{▽}-\overline{u}\cdot ▽+\mu$$

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

$$a(c_h,d_h)=l(d_h)\forall d_h\in V_h$$

$$\begin{array}{rl}a(c_h,d_h)& =\sum _{a,b=1}^{n_{pt}}{d}_{b}a(S_a,S_b)c_a=D^{t}AC\\ l(d_h)& =\sum _{b=1}^{n_{pt}}{d}_{b}l(S_b)=D^{t}L\end{array}$$

$$AC=L$$

$$a_{SGS}(c_h,d_h)=l_{SGS}(d_h)\forall d_h\in V_h$$

$$\begin{array}{rl}e& =c_h-c\\ & =(c_h-c_h^{\prime })+(c_h^{\prime }-c)\\ & =e_h+\lambda \end{array}$$

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

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

$$a(c,d_h)=l(d_h)\forall d_h\in V_h$$

$$a(e,d_h)=0\forall d_h\in V_h$$

$$\begin{array}{rl}a(c_h-c,d_h)+(-{\mathcal{L}}^{\ast }{d}_h,\tau \mathcal{L}{c}_h)& =(-{\mathcal{L}}^{\ast }{d}_h,\tau q)\forall d_h\in V_h\\ a(e,d_h)+(-{\mathcal{L}}^{\ast }{d}_h,\tau \mathcal{L}{c}_h)& =(-{\mathcal{L}}^{\ast }{d}_h,\tau \mathcal{L}c)\forall d_h\in V_h\\ a(e,d_h)+(-{\mathcal{L}}^{\ast }{d}_h,\tau \mathcal{L}e)& =0\forall d_h\in V_h\\ a_{SGS}(e,d_h)& =0\forall d_h\in V_h\end{array}$$

$$a(e,{\eta }_h)+(-{\mathcal{L}}^{\ast }{\eta }_h,\tau \mathcal{L}e)=0$$

$$\begin{array}{rl}\Vert e{\Vert }^2& =a(e,\eta -{\eta }_h)-(-{\mathcal{L}}^{\ast }{\eta }_h,\tau \mathcal{L}e)\\ & =({\mathcal{L}}^{\ast }{\eta }_h,\tau \mathcal{L}e)-a(e,{\lambda }_{\eta })\\ & =({\mathcal{L}}^{\ast }{\eta }_h,\tau \mathcal{L}{e}_h)+({\mathcal{L}}^{\ast }{\eta }_h,\tau \mathcal{L}\lambda )-a(e,{\lambda }_{\eta })\end{array}$$

$$\Vert {\lambda }_c{\Vert }_s\le \bar{C}(p,\Omega )h^{p+1-s}\Vert c{\Vert }_{p+1}$$

$$\Vert {\lambda }_{\eta }{\Vert }_s\le \bar{C}(p,\Omega )h^{p+1-s}\Vert \eta {\Vert }_{p+1}$$

$$\Vert c{\Vert }_2\le C_s\Vert q\Vert$$

$$\Vert \eta {\Vert }_2\le C_s^{\prime }\Vert e\Vert$$

$$\begin{array}{rl}{a}_{SGS}(c_h^{\prime }+e_h,d_h)& =l_{SGS}(d_h)\forall d_h\in V_h\\ a_{SGS}(c_h^{\prime },d_h)+a_{SGS}(e_h,d_h)& =l_{SGS}(d_h)\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)\forall d_h\in V_h\\ a_{SGS}(c_h^{\prime },d_h)-a_{SGS}(\lambda ,d_h)& =l_{SGS}(d_h)\forall d_h\in V_h\end{array}$$

$$\begin{array}{rl}{a}_{SGS}(e_h,d_h)+a_{SGS}(\lambda ,d_h)& =0\forall d_h\in V_h\\ a(e_h,d_h)+(-{\mathcal{L}}^{\ast }{d}_h,\tau \mathcal{L}{e}_h)+a(\lambda ,d_h)+(-{\mathcal{L}}^{\ast }{d}_h,\tau \mathcal{L}\lambda )& =0\forall d_h\in V_h\\ ({\mathcal{L}}^{\ast }{d}_h,\tau \mathcal{L}{e}_h)+({\mathcal{L}}^{\ast }{d}_h,\tau \mathcal{L}\lambda )& =0\forall d_h\in V_h\end{array}$$

$$({\mathcal{L}}^{\ast }{\eta }_h,\tau \mathcal{L}{e}_h)+({\mathcal{L}}^{\ast }{\eta }_h,\tau \mathcal{L}\lambda )=0$$

$$\begin{array}{rl}\Vert e{\Vert }^2& \le \mid a(e,{\lambda }_{\eta })\mid \\ & \le C{h}^{p+1}\Vert e\Vert \Vert c{\Vert }_{p+1}\end{array}$$

$$\boxed{\Vert e\Vert \le C{h}^{p+1}\Vert c{\Vert }_{p+1}}$$

$$\begin{array}{rl}\Vert e\Vert & \le C{h}^2\Vert c{\Vert }_2\\ & \le C{C}_s{h}^2\Vert q\Vert \end{array}$$

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

$$\begin{array}{rl}\Vert e{\Vert }^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{array}$$

$$\begin{array}{rl}a(c_h,d_h)+(-{\mathcal{L}}^{\ast }{d}_h,\tau ((\mathcal{L}{c}_h-q){)}_{{\Omega }^{\prime }}& =l(d_h)\forall d_h\in V_h\\ a(c_h,d_h)+(-{\mathcal{L}}^{\ast }{d}_h,\tau r_h{)}_{{\Omega }^{\prime }}& =l(d_h)\forall d_h\in V_h\end{array}$$

$$a(c_h,{\eta }_h)+(-{\mathcal{L}}^{\ast }{\eta }_h,\tau r_h{)}_{{\Omega }^{\prime }}-l({\eta }_h)=0$$

Full text

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

[TABLE]

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}$,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

This is also known as Galerkin orthogonality.

Now we are going to bring SGS formulation into estimation by subtracting (10) from (8) as follows:

[TABLE]

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,

[TABLE]

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

[TABLE]

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,

[TABLE]

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.

[TABLE]

Another results are strong stability conditions [5],

[TABLE]

and similarly for the dual problem

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

Hence,

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}$

[TABLE]

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),

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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,

[TABLE]

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:

[TABLE]

Hence,

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

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,

[TABLE]

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,

[TABLE]

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,

[TABLE]

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.

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.