Local projection stabilized finite element methods for advection–reaction problems

An a priori analysis for generalized local projection stabilized finite element approximations for the solution of an advection–reaction equation is presented in this article. The stability and a priori error estimates are derived for the conforming, and nonconforming (Crouzeix–Raviart) approximations in the local projection streamline derivative norm. Finally, the validation of the proposed stabilization scheme and verification of the derived estimates are presented with appropriate numerical experiments.


Introduction
Advection-reaction equations arise in many engineering and industrial applications.The numerical solution of these equations has been of interest for several decades.It is well-known that applying the standard Galerkin finite element method (FEM) to the advection-reaction equations induces spurious oscillations in the numerical solution.Nevertheless, the Galerkin approximation's stability and accuracy can be enhanced by employing a stabilization method.Over the years, several stabilization methods such as the Streamline Upwind Petrov-Galerkin methods (SUPG), Least-Squares methods, Residual-free Bubbles, Continuous Interior Penalty (CIP), Subgrid B Deepika Garg d.garg@ucl.ac.uk; deepika.lpu.pbi@gmail.comSashikumaar Ganesan sashi@iisc.ac.inViscosity, and Local Projection Stabilization (LPS) have been proposed.The relation between the different approaches is also well-understood in most cases.The basic idea of stabilization is to stabilize the Galerkin variational formulation so that the discrete approximation is stable, the method is (weakly) consistent and, consequently, convergent.
The essential idea in SUPG is to add a weighted residual to the Galerkin variational formulation to make it globally stable and consistent.The SUPG method has been well-established for conforming and nonconforming FEM [9, 21, 24-26, 28, 32, 33].In the early 1970s, the Least-Square method has become popular within the numerical analysis community [6,7].However, it has already been published in Russian literature [19].The Least-Square method is inspired by the minimal residual, a technique from linear algebra [7,30].The Residual-free Bubble stabilization method is based on the Galerkin FEM with an enriched basis on each element [8].In a particular case, the Galerkin variational formulation with an enriched basis can be shown equivalent to the SUPG method with piecewise linear finite elements [1].Continuous Interior Penalty (CIP) is another efficient and well-studied stabilization technique.The basic idea in CIP stabilization (also known as edge stabilization in the literature) is to penalize the jump of the gradient across the cell interfaces [10,11,14].CIP has also been studied for the hp-finite elements [13] and the Friedrichs' systems [12].
This article focuses on stabilization by local projection for the advection-reaction equations.The local projection stabilization method has been introduced by Becker and Braack [3] and Braack and Burman [4].The stabilization term in the local projection method is based on a projection of the finite element space that approximates the unknown into a course (discontinuous) space, see [3,4].This technique has been initially studied for fluid flow problems with Stokes-like models.The macro grid approach approximates pressure and velocity using the same finite element spaces [3,22,31].Later, the LPS method on a single mesh with enriched finite element spaces was proposed and extended to various incompressible flow problems [4,23,29,34].Moreover, in a particular case, the SUPG method can be recovered from the LPS method with piecewise linear functions enriched polynomial Bubble space on triangles and with an appropriate SUPG-parameter, see [23].LPS method adds symmetric stabilization terms and contains fewer stabilization terms than residualbased stabilization methods.
furthermore, a generalization of the local projection stabilization allows defining local projection spaces on overlapping grids.Neither macro grid nor enrichment of spaces is needed in generalized local projection stabilization (GLPS).This approach has been introduced and studied for the convection-diffusion problem in [27] with conforming finite element space, recently in [18] with conforming, and nonconforming finite element spaces, and for the Oseen problem in [29].This paper presents an a priori analysis for the generalized local projection stabilization scheme with conforming and nonconforming finite element spaces for an advection-reaction equation.In the absence of the diffusion term in the advection-reaction equation, a different approach is needed in the numerical analysis compared to the earlier studies.In particular, a relatively stronger norm, the local projection streamline diffusion (LPSD) norm, is used in our analysis.Since the LPDS norm contributes to the streamline derivative term [27], it provides control with respect to streamline derivatives.Moreover, it has been shown that the LPSD norm is equivalent to the SUPG norm for an appropriate choice of mesh-dependent parameter [14].
GLPS approach is further extended to the nonconforming discretization for the advection-reaction equation.The standard nonconforming formulation together with the GLP stabilization is not sufficient to have a coercive discrete formulation.To achieve coercivity, additional weighted edge integrals with the jumps and the averages of the discrete solution at the interfaces need to be added to the GLPS bilinear form in the nonconforming finite element formulation as in SUPG [25,26].Though the nonconforming GLPS is challenging compared to the conforming scheme, it is preferred in parallel computing.Since the nonconforming shape functions have local support at most in two cells, the sparse matrix stencil will be smaller.Therefore, communication across MPI processes in parallel computing is minimal, resulting in better scalability.Moreover, numerical comparisons of the conforming and nonconforming approaches demonstrate that the nonconforming approach is marginally more robust with respect to the stabilization parameter than the conforming approach (see Fig. 4).
The article's outline is as follows: Sect. 2 introduces the model problem and GLPS formulation.In Sect.3, we derive a stability estimate of the conforming GLPS scheme and establish an optimal a priori error estimates.Section 4 focuses on the nonconforming GLPS, in which we derive the stability results and obtain an optimal a priori error estimates.Section 5 presents a set of numerical experiments to support our theoretical estimates.

The model problem
Let Ω ⊂ R 2 be a bounded polygonal domain with boundary ∂Ω.Consider the following advection-reaction equation with a boundary condition: is the source term and g ∈ L 2 (∂Ω − ) is a boundary data and ∂Ω − denotes the inflow part of the boundary of Ω namely Further, n is the unit outward normal to the boundary.We assume that there exist α > 0 such that

Variational formulation
Let L 2 (Ω) and H m (Ω), m > 0 be the standard Sobolev spaces and Note that the functions in V have traces in L 2 (∂Ω; |b • n|) [11,17], which can be defined by We now derive a variational form of the model problem in an usual way.Multiplying the model problem with a test function v ∈ V and after integrating over Ω, the variational form of the model problem (2.1) reads: where, Here, (•, •) is the L 2 (Ω) inner product, u := 1 2 (|u| − u) and u ⊕ := 1 2 (|u| + u), where |u| is the modulus function of u.The Banach-Ne čas-Babuška theorem [20, pp. 83] guarantee that the model problem is well-posed with the space V , for more details; see [20, pp. 230].

Remark 1
The analysis is presented for a two-dimensional case for simplicity.Nevertheless, the study is independent of the dimension, whereas faces instead of edges need to be used to extend to three-dimensions.

Finite element space
Let T h be a collection of non-overlapping quasi-uniform triangles obtained by a decomposition of Ω.Let h K = diam(K) for all K ∈ T h and the mesh-size h = max K ∈T h h K .Assume that the family of meshes T h is shape regular, i.e., there exist c K such that for all h > 0 and for all where ρ K denotes the radius of the largest inscribed ball in K .Let E h = E I h ∪ E B h be the set of all edges in T h , where E I h and E B h are the sets of all interior and boundary Fig. 1 The edge E = ab is shared by two neighboring triangles K + and K − and n is the unit outward normal to K + , (left), edge patch M E and node patch M a (right) edges respectively.Let h E be the length of any edge.Further, for each edge E in E h , a unit normal vector n is associated, this is taken to be the unit outward normal to Ω for all E ∈ E h .Suppose K + (E) and K − (E) are neighbours of the interior edge E ∈ E I h , then the normal vector n is oriented from K + (E) and K − (E), (Fig. 1).The average and jump of a function v on the edge E can be defined as h be the set of all vertices in V h , where V I h and V B h are the sets of all interior and boundary vertices respectively.For any a ∈ V h , M a (patch of a) denotes the union of all cells that share the vertex a. Further, define h a = diam(M a ) for all a ∈ V h .Moreover, for any E ∈ E h , M E (patch of E) denotes the union of all cells that share the edge E, (Fig. 1).
The following norm is used in the analysis.Let the piecewise constant function h T is defined by h T | K = h K and s ∈ R and m ≥ 0 for all u ∈ H m (T h ).
Suppose I (a) denotes the index set for all K l elements, so that K l ⊂ M a .Then, the local mesh-size associated to M a is defined as where card(I (a)) denotes the number of elements in M a .
We next define a piecewise polynomial space as where P k (K ), k ≥ 0, is the space of polynomials of degree at most k over the element K .Further, define a conforming finite element space of piecewise linear as and a nonconforming Crouzeix-Raviart finite element space of piecewise linear is defined as Next, the following technical results of finite element analysis are recalled.

Lemma 1
Trace inequality [17, pp. 27]: Suppose E denotes an edge of K ∈ T h .For ) where diam(M ) and |M | denote the diameter and measure of domain M .In particular, for every vertex a ∈ V h and every function v ∈ H 1 (M a ), it holds [35, pp. 8, 91] where h a is the diameter of M a .Here, the constant C equals to 1/π for convex domain M a .

Lemma 4 [5, Lemma 5] Let (H, • H ) be a Hilbert space where the norm of y ∈ H is defined by means of two semi-norms
and κ H ≤ y H with constants c 0 , c 2 > 0 and c 1 ≥ 0.Then, the bilinear form fulfills the inf-sup condition where σ > 0 is arbitrary.
Note that throughout this paper, C (sometimes subscript) denotes a generic positive constant, that may depend on the shape-regularity of the triangulation but is independent of the mesh size.Moreover, the L 2 (Ω) and L ∞ (Ω) norms are respectively denoted by u and u ∞ .

Discrete formulation
The conforming discrete solution of (2.3) is a function For any a ∈ V h , define a fluctuation operator κ a : where |M a | denotes the measure of M a .For each a ∈ V h , let β a := βh a for some +C where the constant C has the same unit as the velocity field, that is m/s.We now define a conforming local projection stabilization term as Using this stabilization, the conforming generalized local projection stabilized discrete form of (2.3) reads as: Find where, Further, we introduce a local projection norm for v h ∈ V c h as: and the local projection streamline derivative norm for v h ∈ V c h is defined as: Further, the L 2 -orthogonal projection J c h : L 2 (Ω) → V c h satisfies the following approximation properties [2].
Further, the L 2 -orthogonal projection operator exhibits the following H 1 and L 2stability properties [2,Theorem 4.1,4.2]: Moreover, the main result of this subsection is the following theorem, which ensures that the discrete bilinear form is well-posed [20].+C .Then, the discrete bilinear form (3.15) satisfies the following inf-sup condition for some positive constant γ , independent of h: Proof Using Lemma 4, the inf-sup condition can be proved in three steps: First, consider the bilinear form in (3.15) with v h = u h .Applying an integration by parts to the first term of the bilinear form and an application of (2.2) lead to, Further, the control of h Let us now estimate these four terms.Using the canonical representation of the basis function φ a at the node a ∈ V h for the mesh Using the orthogonality property of L 2 -projection (3.20) with the test function Choosing the constant C a = 1 |M a | M a b • ∇u h dx, and applying the Cauchy-Schwarz inequality, (3.21) and Young's inequality: the constant C in the above estimate depends on b ∞ .The second term is estimated by applying the Cauchy-Schwarz inequality followed by (3.21) and an inverse inequality, The constant C in (3.24) depends on b ∞ and α −1 .The third term is handled by applying the Cauchy-Schwarz inequality, trace inequality (2.6), (3.21) and Young's inequality, The constant C depends on b ∞ .Next, applying the Cauchy-Schwarz inequality to the fourth term to obtain, The second term of (3.25) is estimated by using the boundedness of local projection operator, an inverse inequality (2.7), stability of the projection estimates (3.21), and .
The constant C in the above estimates depends on b ∞ .Then, it follows that, The constant C in the above estimates depends on b ∞ .Put together, (3.23) leads to, The constant C in the estimates (3.27) depends on b ∞ .Note that (3.28) Now, estimate the four terms of (3.28) individually.Using the stability of the projection operator (3.21) and an inverse inequality, The second term is estimated by using trace inequality and (3.21) , The last two terms are handled by using boundedness of the local projection operator, an inverse inequality (2.7) and projection estimates (3.21), that is, Finally, put together, we get, The constant C in (3.29) depends on b ∞ .The selection of z h is where J c h is as defined in Lemma 5. Finally, the result follows by Lemma 4.

A priori error estimates
Lemma 6 Suppose u ∈ H 2 (Ω) and β a = βh a for some β| Proof Consider the terms in the LPSD norm defined in (3.17), Let us now estimate the terms on the right-hand side of (3.30) individually.The first and second terms are estimated by using the approximation estimates (3.18), The constant C depends on b ∞ .The third term of (3.30) is handled by using the trace inequality (3.19) over each edge, .
Note that the constant C in the above estimates depends on the b ∞ .The last term is estimated by using boundedness of the local projection operator and .
The constant C depends on b ∞ .The combination of the above estimates concludes the proof.
Lemma 7 Suppose u ∈ H 2 (Ω) and β a = βh a for some β| Proof Applying an integration by parts to the first term of the discrete bilinear form in (3.15) to obtain, The first term is estimated by using the Cauchy-Schwarz inequality and L 2 -projection property (3.18), and, The constant C in (3.32) depends on α −1/2 .The third term is handled by applying the Cauchy-Schwarz inequality, boundedness of the local projection, approximation estimates (3.18) and Applying the Cauchy-Schwarz inequality, trace inequality (2.5) and approximation estimates (3.18), Combining the above estimates leads to (3.31) and it concludes the proof.
Theorem 2 Let u ∈ H 2 (Ω) be the solution of (2.3) and u h ∈ V c h be the discrete solution of (3.14).Let β a = βh a for some β| M a = Proof By adding and subtracting the interpolation operator J c h u, we decompose the error as follows: In the second term of (3.33), using the estimate of Theorem 1, The weak formulation (2.4) and (3.15) imply that Moreover, the Cauchy-Schwarz inequality implies +C .Using the Poincaré inequality (2.9) for every vertex a ∈ V h , It follows that,

Nonconforming finite element discretization
The nonconforming discrete solution of (2.3) is a function where, Here, ∇ h denotes the piecewise (element-wise) gradient operator.For each E ∈ E h , define the fluctuation operator where the constant C has the same unit as the velocity field, that is m/s.We now define a nonconforming local projection stabilization term (4.38) using this term, the nonconforming generalized local projection stabilized discrete form of (2.3) reads as: Find u h ∈ V nc h such that where, Further, we define a nonconforming local projection norm by and nonconforming local projection streamline derivative norm by for all v h ∈ V nc h .Moreover, the L 2 -projection J nc h : L 2 (Ω) → V nc h satisfies the approximation properties stated in (3.18)-(3.21)for shape-regular triangulation, see [16,Theorem 4], [18,Lemma 7.1].+C .Then, the discrete bilinear form (4.39) satisfies the following inf-sup condition for a positive constant ν, independent of h:

.43)
Proof Using Lemma 4, the inf-sup condition can be proved in three steps: The key steps to derive the above estimates are as follows: Choosing first v h = u h as a test function in (4.40), Further, the control h 4.40), we have by adding and subtracting h Most of the estimates of (4.44) can be derived in a similar way as shown in (3.23).
The constant C in the above estimates depends on b ∞ .Now, it is sufficient to estimate the last two terms of (4.44).Using the Cauchy-Schwarz inequality, .
At the edge E, the jump term contributes to both the triangles sharing that edge; using the trace inequality (2.6) and (3.21), .
We then get, Similarly, the next term is estimated as: Combining all these estimates and (4.44) lead to, constant C in the above estimates on b ∞ .In particular, the inequality holds for where J nc h is the projection operator and the constant C in the above equation depends on the reaction coefficient, the constant in the inf-sup condition, L ∞ -norm of b.Rest of the proof can be derived in a similar way as in the proof of (3.28)-(3.29).

.46)
Proof Most of the estimates of the term (4.42) follow from Lemma 6; hence, we need to handle the last term of (4.41), The constant C depends on b ∞ .At the edge E, the jump term contributes to both the triangles sharing that edge; using the trace inequality (2.5), Squaring and summing up all the inner edges and using (3.18), we have here, the constant C depends on b ∞ .The result follows by combining all the above estimates.

.47)
Proof Using integration parts in the first term of (4.40), The first four terms of the bilinear form A nc h (u − J nc h u, v h ) can be estimated similarly as in Lemma 7.Moreover, the last two terms are handled by applying the Cauchy-Schwarz inequality, . Since ), and at the edge E, the jump term contributes to both the triangles sharing that edge; using the trace inequality (2.5), Squaring and summing up all the inner edges and using (3.18), Similarly, The constant C in the above estimates depends on b ∞ .Combining all these estimates leads to (4.47) and it concludes the proof.

Theorem 4 Let u ∈ H (Ω) be the solution of continuous problem (2.3) and u h ∈ V nc h be the solution of discrete problem (4.39). Further, let β
Proof By adding and subtracting the interpolation operator J nc h u, we decompose the error as follows: Consider the second term of (4.49) and using the Theorem 3, The weak formulation (2.4) and (4.39) imply that Moreover, the Cauchy-Schwarz inequality implies +C .Using the Poincaré inequality (2.9) for every vertex a ∈ V h , It follows that,

Numerical results
This section presents an array of numerical results to support the analysis presented in the previous sections.Numerical solutions of all examples are computed on a hierarchy of uniformly refined triangular meshes having 256, 1024, 4096, 16384, and 65536 elements, respectively.Triangulation used for computations in section 5 are shown in Fig. 2.

Example 1 (Smooth solution)
Consider the model problem (2.1) with Ω = (0, 1) 2 , coefficients b = (3, 2), μ = 2 and homogeneous Dirichlet boundary condition.The source term f is chosen such that the solution satisfies the model problem.Further, the stabilization parameters for conforming and nonconforming FEMs are chosen as β a = 0.1h a and β E = 0.1h E , respectively.
The right hand side of Fig. 2 depicts the nonconforming stabilized finite element solution computed on a mesh with h = 2 −6 .The Tables 1 and 2 present the errors of GLPS conforming and nonconforming finite element approximations respectively, in L 2 −norm, H 1 −seminorm and the local projection streamline-derivative norm defined in (3.17) and (4.42).A second-order convergence can be observed in L 2 -norm and first-order convergence in H 1 -seminorm.Moreover, the convergence order 1.5 was obtained in |||•||| L P SD norm.Also, the log-log plots of the errors in Fig. 3 show the convergence behaviour of errors in the conforming and the nonconforming approximation, and confirming our theoretical estimates.

Sensitivity with respect to problem and stabilization parameter
Even though the convergence rates in Theorems 2 and 4 hold in principle for any finite β > 0, optimizing the stabilizing parameters can have some impact on the quality of the obtained numerical solution.The left plot of Fig. 4 shows that μ has almost no influence on the error.We will set μ = 1 from now on.In centre Fig. 4, we observe that the convection field significantly influences the error bound in both cases.The convection term introduces significant numerical diffusion for b ∞ ≥ 10 (see (3.13), (4.38)), which could affect the accuracy of the numerical solution and the error shoots up for both the conforming and nonconforming FE approximations.Further, it can be noticed from the right plot of Fig. 4 that the error in terms of β behaves like a well with an approximation minimum at β ∈ [0.1, 1] for the conforming and a broader minimum at β ∈ [0.001, 0.1] for the nonconforming finite element method.These numerical experiments show that the nonconforming finite element method is slightly more robust with respect to the stabilization parameter.Both FE approximations require stabilization, so the error blows up as the β approaches zero.Further, the stabilization parameters for conforming and nonconforming finite element approximations are chosen as β a = 0.1h a and β E = 0.1h E , respectively.The Figs. 5 (a) and (b) show the nonconforming Galerkin and GLPS finite element solutions.We can observe that the spurious oscillation in the Galerkin solution is suppressed in GLPS approximation.Further, the Tables 3 and 4 present the errors and convergence behaviour of the conforming and nonconforming stabilized finite element solutions respectively.Moreover, Fig. 6 depicts the obtained optimal order of convergence in both the conforming and the nonconforming approximations.This solution possesses a circular internal layer on the circumference of the circle, centered at (0.5,0.5) and radius 0.25, in the unit square domain.The conforming and the nonconforming approximations are obtained with the stabilization parameters β a = 0.06h a and β E = 0.05h E , respectively.The Fig. 7 depicts the GLPS conforming stabilized finite element solution on a mesh with h = 2 −7 .We can observe that the conforming stabilized scheme approximates the solution well and retains its inner circular layer.A similar result is obtained with the nonconforming GLPS finite element approximation.Next, the Table 5 displays the errors and convergence behaviour in |||•||| norm as defined in (3.17) and (4.42) for the GLPS conforming and nonconforming finite approximations.Figure 7 presents the convergence plots in the conforming and nonconforming approximations and supports the theoretical estimates.Even though discontinuous boundary data [15,Example 2] is not considered in our numerical analysis, this example is considered to examine the robustness of the proposed scheme.The stabilization parameters for the conforming and the nonconforming approximations are chosen as β a = 0.7h a and β E = 0.7h E , respectively.The Fig. 8 depicts the conforming and the nonconforming stabilized finite element solutions on a mesh with h = 2 −6 .Since the boundary conditions are imposed weakly in the current scheme, boundary layers are not resolved in our computations.Nevertheless, the interior layer is captured well with the proposed GLPS method.Though small overshoots and undershoots are observed near the interior layer, there are no oscillations in the solution away from the layer.It shows the robustness of the proposed scheme.

Summary
The stability and convergence estimates for generalized local projection stabilized finite element scheme for the advection-reaction problem with conforming and nonconforming spaces was derived in this paper.The GLPS allows projection spaces on overlapping sets and avoids needing a two-level mesh or an enrichment of the finite element space.In particular, optimal a priori error estimates were established for conforming and nonconforming approximations in the local projection streamline derivative norm.the numerical examination of the sensitivity of the stabilization parameter indicated that the nonconforming approach was slightly more robust with respect to the stabilization parameter.The accuracy and robustness of the proposed scheme were shown numerically with suitable examples.Moreover, an extension of this work to the flow problem is planned.

Fig. 4
Fig. 4 Comparison of the conforming (dashed lines) and the non-conforming (solid lines) methods regarding the dependency of the error bounds on the problem parameters and the sensitivity to the stabilization parameter of the Example 2

Fig. 5 a
Fig. 5 a Nonconforming Galerkin and b nonconforming stabilized finite element β E = 0.1h E solutions of the Example (2)

Fig. 6
Fig. 6 Convergence behaviour of the stabilized finite solution of the Example 2

Fig. 7 Fig. 8
Fig.7 The errors of finite element approximations and conforming stabilized finite element solution of the example(3)

Table
Errors and convergence to the conforming FE of Example

Table 2
Errors and convergence orders to the nonconforming solution of

Table 3
Errors and orders of convergence of conforming FE solution of Example

Table 4
Errors and orders of convergence of nonconforming FE solution of Example

Table 5
Errors and convergence orders to the GLPS finite element approxomations of the Example 3