Local projection stabilized finite element methods for advection–reaction problems

Abstract

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.

1 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 Viscosity, 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,25,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 residual-based 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.

2 Finite elements for the advection–reaction equation

2.1 The model problem

Let $\varOmega \subset {{\mathbb {R}}}^2$ be a bounded polygonal domain with boundary $\partial {\varOmega }$. Consider the following advection–reaction equation with a boundary condition:

$$\begin{aligned} \begin{array} {rcl} \mu {u} +\textbf{b} \cdot \nabla {u}&{} =&{} f \ \; \text{ in } \; \varOmega , \\ u &{}=&{} g \ \; \text{ on } \; \ \partial {\varOmega }^{-}. \end{array} \end{aligned}$$

(2.1)

Here, u is an unknown scalar function, ${\textbf {b}} \in [W^1_\infty (\varOmega )]^2$ is the advective velocity, $\mu \in L^\infty (\varOmega )$ the reaction coefficient, $f\in L^2(\varOmega )$ is the source term and $g \in L^2(\partial \varOmega ^{-})$ is a boundary data and $\partial {\varOmega }^{-}$ denotes the inflow part of the boundary of $\varOmega$ namely

$$\begin{aligned} \partial {\varOmega }^{-}:= \{x \in \partial {\varOmega } \ | \ {\textbf {b}}(x)\cdot {\textbf {n}}(x) <0 \ \}. \end{aligned}$$

Further, ${\textbf {n}}$ is the unit outward normal to the boundary. We assume that there exist $\alpha > 0$ such that

$$\begin{aligned} \mu _{0} := \left( \mu - \frac{1}{2} \text{ div } {{\textbf {b}}}\right) \ge \alpha >0 \ \ a.e. \text { in } \ \varOmega . \end{aligned}$$

(2.2)

2.2 Variational formulation

Let $L^{2}(\varOmega )$ and $H^{\textrm{m}}(\varOmega ),~{\textrm{m}}>0$ be the standard Sobolev spaces and

$$\begin{aligned} V= \{v \in L^{2}(\varOmega ) \ | {\textbf {b}} \cdot \nabla {v} \in L^{2}(\varOmega ) \ \}. \end{aligned}$$

Note that the functions in V have traces in $L^{2}(\partial {\varOmega }; |{\textbf {b}} \cdot {\textbf {n}}|)$ [11, 17], which can be defined by

$$\begin{aligned} L^{2}(\partial {\varOmega }; |{\textbf {b}} \cdot {\textbf {n}}|):=\left\{ v \ \text {is measuarable function on} \ \partial {\varOmega } \ | \ \int _{\partial {\varOmega }} |\textbf{b}\cdot \textbf{n}|v^{2} \;\textit{ds}< \infty \ \right\} . \end{aligned}$$

We now derive a variational form of the model problem in an usual way. Multiplying the model problem with a test function $v\in V$ and after integrating over $\varOmega$, the variational form of the model problem (2.1) reads:

Find ${u} \in V$ such that

$$\begin{aligned} a(u,v) = l(v)\ \ \text {for all} \ v \in V, \end{aligned}$$

(2.3)

where,

$$\begin{aligned} a(u,v)&:= ({\textbf {b}} \cdot \nabla {u}, v) + (\mu {u}, v) +\int _{\partial {\varOmega }}({{\textbf {b}} \cdot {\textbf {n}}})^{\circleddash } uv \,ds,\nonumber \\ l(v)&:= (f,v) + \int _{\partial \varOmega ^-} (\textbf{b} \cdot \textbf{n})^{\circleddash } gv \,ds. \end{aligned}$$

(2.4)

Here, $(\cdot ,\cdot )$ is the $L^{2}(\varOmega )$ inner product, $u^\ominus :=\frac{1}{2}(|u|-u)$ and $u^\oplus :=\frac{1}{2}(|u|+u)$, where $\left| u\right|$ is the modulus function of u. The Banach–Ne$\check{c}$as–Babu$\check{s}$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.

2.3 Finite element space

Let $\mathscr {T}_{h}$ be a collection of non-overlapping quasi-uniform triangles obtained by a decomposition of $\varOmega$. Let $h_K=\textrm{diam}(K)$ for all $K \in \mathscr {T}_{h}$ and the mesh-size $h = \text{ max}_{K\in \mathscr {T}_h} h_K$. Assume that the family of meshes $\mathscr {T}_{h}$ is shape regular, i.e., there exist $c_K$ such that for all $h>0$ and for all $K \in \mathscr {T}_h$

$$\begin{aligned} \frac{h_K}{\rho _K}\le c_K, \end{aligned}$$

where $\rho _K$ denotes the radius of the largest inscribed ball in K. Let $\mathscr {E}_{h}=\mathscr {E}_{h}^{I}\cup \mathscr {E}_{h}^{B}$ be the set of all edges in $\mathscr {T}_{h}$, where $\mathscr {E}_{h}^{I}$ and $\mathscr {E}_{h}^{B}$ are the sets of all interior and boundary edges respectively. Let $h_E$ be the length of any edge. Further, for each edge E in $\mathscr {E}_{h}$, a unit normal vector ${\textbf {n}}$ is associated, this is taken to be the unit outward normal to $\varOmega$ for all $E\in \mathscr {E}_{h}$. Suppose $K^{+}(E)$ and $K^{-}(E)$ are neighbours of the interior edge $E\in \mathscr {E}_{h}^{I}$, then the normal vector ${\textbf {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

$$\begin{aligned} \{\hspace{-1.0pt}v \hspace{-1.0pt}\}=\frac{1}{2} \left( v^+|_E+v^-|_E\right) ,[\hspace{-0.5pt}v \hspace{-0.5pt}]:= v^+|_E-v^-|_E, \qquad \end{aligned}$$

where $v^{\pm }:=v|_{K_{\pm }}$. Let $\mathscr {V}_h:=\mathscr {V}_{h}^I\cup \mathscr {V}_{h}^B$ be the set of all vertices in $\mathscr {V}_{h}$, where $\mathscr {V}_{h}^{I}$ and $\mathscr {V}_{h}^{B}$ are the sets of all interior and boundary vertices respectively. For any $a\in \mathscr {V}_h$, ${\mathscr {M}}_a$ (patch of a) denotes the union of all cells that share the vertex a. Further, define $h_a=\textrm{diam}(\mathscr {M}_{a})$ for all $a\in \mathscr {V}_h$. Moreover, for any $E\in \mathscr {E}_{h}$, ${\mathscr {M}}_E$ (patch of E) denotes the union of all cells that share the edge E, (Fig. 1).

Fig. 1

The edge $E=ab$ is shared by two neighboring triangles $K^{+}$ and ${K^{-}}$ and $\textbf{n}$ is the unit outward normal to $K^+,$ (left), edge patch ${\mathscr {M}}_E$ and node patch ${\mathscr {M}}_a$ (right)

The following norm is used in the analysis. Let the piecewise constant function $h_{\mathscr {T}}$ is defined by $h_{\mathscr {T}} |_K = h_K$ and $s \in \mathbb {R}$ and ${\textrm{m}} \ge 0$

$$\begin{aligned} \left\| h_{\mathscr {T}}^{s} u \right\| _{\textrm{m}} = \left( \sum _{K \in \mathscr {T}_h } h_{K}^{2s} \left\| u \right\| ^{2}_{H^{\textrm{m}}(K)} \right) ^{\frac{1}{2}} \text { for all} \ u \in \ H^{\textrm{m}}({\mathscr {T}}_h). \end{aligned}$$

Suppose I(a) denotes the index set for all $K_l$ elements, so that $K_l \subset {\mathscr {M}}_a$. Then, the local mesh-size associated to ${\mathscr {M}}_a$ is defined as

$$\begin{aligned} {\hat{h}}_a := \frac{1}{\text{ card }(I(a))} \sum _{l\in I(a)} h_l, \quad \text { for each} \ a\in \mathscr {V}_h, \end{aligned}$$

where $\text {card}(I(a))$ denotes the number of elements in ${\mathscr {M}}_a$.

We next define a piecewise polynomial space as

$$\begin{aligned} \mathbb {P}_{k}(\mathscr {T}_h):=\left\{ v\in L^{2}(\varOmega ): v|_K\in \mathbb {P}_k(K)\quad \forall K\in \mathscr {T}_h\right\} , \end{aligned}$$

where $\mathbb {P}_k(K)$, $k\ge 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

$$\begin{aligned} V^{c}_{h} := \left\{ v \in H^{1}(\varOmega ) \ : \ v|_K \in \mathbb {P}_1(K) ~~ \forall ~ K \in \mathscr {T}_{h} \right\} , \end{aligned}$$

and a nonconforming Crouzeix–Raviart finite element space of piecewise linear is defined as

$$\begin{aligned} V_h^{nc}&:=\left\{ v\in L^2(\varOmega ): v|_K\in \mathbb {P}_{1}(K), \quad \int _{E} [\hspace{-0.5pt}v \hspace{-0.5pt}] \;\textit{ds}=0, \quad \text {for all} \quad E\in \mathscr {E}^{I}_h \right\} . \end{aligned}$$

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 \in \mathscr {T}_{h}$. For $v|_{K} \in H^{1}(K)$ and $v_h\in \mathbb {P}_{k}(\mathscr {T}_{h})$, there holds,

$$\begin{aligned} \Vert v\Vert _{L^2(E)}&\le C \left( h_K^{-1/2} \Vert v\Vert _{L^2(K)} + h_K^{1/2} \Vert \nabla v\Vert _{L^2(K)}\right) , \end{aligned}$$

(2.5)

$$\begin{aligned} \Vert v_h\Vert _{L^2(E)}&\le C h_K^{-1/2} \Vert v_h\Vert _{L^2(K)} . \end{aligned}$$

(2.6)

Lemma 2

Inverse inequality [17, pp. 26]: Let $v \in \mathbb {P}_{k}(\mathscr {T}_{h})$ for all $k \ge 0$. Then,

$$\begin{aligned} \left\| \nabla {v} \right\| _{L^{2}(K)} \le C h^{-1}_{K} \left\| v \right\| _{L^{2}(K)}. \end{aligned}$$

(2.7)

Lemma 3

Poincaré inequality [35, pp. 91]: For a bounded and connected polygonal domain ${\mathscr {M}}$ and for any $v\in H^1({\mathscr {M}})$,

$$\begin{aligned} \left\| v-\frac{1}{|{\mathscr {M}}|}\int _{{\mathscr {M}}} v\;\textit{dx} \right\| _{L^2({\varOmega })}\le C\ \text {diam}({\mathscr {M}})\left\| \nabla v \right\| _{L^2({\mathscr {M}})}, \end{aligned}$$

(2.8)

where $\text {diam}({\mathscr {M}})$ and $|{\mathscr {M}}|$ denote the diameter and measure of domain ${\mathscr {M}}$.

In particular, for every vertex $a\in \mathscr {V}_h$ and every function $v \in H^1({\mathscr {M}_a})$, it holds [35, pp. 8, 91]

$$\begin{aligned} \left\| v-\frac{1}{|{\mathscr {M}}_a|}\int _{{\mathscr {M}}_a} v\;\textit{dx} \right\| _{L^{2}({\mathscr {M}}_a)}\le Ch_a\left\| \nabla v \right\| _{L^{2}({\mathscr {M}}_a)}, \end{aligned}$$

(2.9)

where $h_a$ is the diameter of ${\mathscr {M}}_a$. Here, the constant C equals to $1/ \pi$ for convex domain ${\mathscr {M}}_a$.

Lemma 4

[5, Lemma 5] Let (H,$\left\| \cdot \right\| _H$) be a Hilbert space where the norm of $y \in H$ is defined by means of two semi-norms $|\cdot |_a$ and $|\cdot |_b$ as $\left\| y \right\| ^2_H=|\cdot |^2_a+|\cdot |^2_b$. For a given bilinear form $B:H\times H \rightarrow \mathbb {R}$, we assume

$$\begin{aligned}&\forall ~~ y \in H:\quad \quad&B(y,y) \ge c_0 |y|^2_a\end{aligned}$$

(2.10)

$$\begin{aligned}&\forall ~~ y \in H \exists \varkappa \in H:\quad \quad&B(y,\varkappa ) \ge c_2 |y|^2_b-c_1 |y|^2_a \end{aligned}$$

(2.11)

and $\left\| \varkappa \right\| _H \le \left\| y \right\| _H$ with constants $c_0, c_2 >0$ and $c_1 \ge 0$. Then, the bilinear form fulfills the inf-sup condition

$$\begin{aligned}\forall ~~ y \in H \ \exists z \in H \setminus \{0\}:\quad \quad B(y,y) \ge \gamma |y|_H |z|_H,\end{aligned}$$

with the constant

$$\begin{aligned}\gamma =min\Bigl \{\frac{c_2}{1+\sigma }, \frac{c_2}{1+(c_1+c_2)/c_0}\Bigr \} >0,\end{aligned}$$

where $\sigma >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(\varOmega )$ and $L^{\infty }(\varOmega )$ norms are respectively denoted by $\left\| u \right\|$ and $\left\| u \right\| _{\infty }$.

3 Conforming finite element discretization

3.1 Discrete formulation

The conforming discrete solution of (2.3) is a function $u_h \in V^{c}_{h}$ such that

$$\begin{aligned} a(u_h,v_h) = l(v_h) \ \hbox { for all} \ v_h \in V_{h}^{c}. \end{aligned}$$

(3.12)

For any $a \in \mathscr {V}_{h}$, define a fluctuation operator $\kappa _{a}: V\rightarrow L^{2}({\mathscr {M}}_{a})$ such that

$$\begin{aligned} \kappa _{a}(u):= \textbf{b}\cdot \nabla {u} \ - \frac{1}{\vert {{\mathscr {M}}_{a}\vert }} \int _{{\mathscr {M}}_{a}} {\textbf{b}\cdot \nabla {u}} \;\textit{dx}, \end{aligned}$$

where $\vert {{\mathscr {M}}_{a}}\vert$ denotes the measure of ${{\mathscr {M}}_{a}}$. For each $a \in \mathscr {V}_h$, let $\beta _{a}:=\beta h_{a}$ for some $\beta = \frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_a)}+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

$$\begin{aligned} S_{h}^c(u_h, v_h):= \sum _{a \in \mathscr {V}_{h}} \beta _{a} \big ({\kappa _{a}({u_h})}, {\kappa _{a}({v_h})\big )}_{L^{2}(\mathscr {M}_{a})}. \end{aligned}$$

(3.13)

Using this stabilization, the conforming generalized local projection stabilized discrete form of (2.3) reads as:

Find $u_{h} \in V^{c}_{h}$ such that

$$\begin{aligned} A^{c}_{h}(u_{h},v_{h}) = l(v_{h}) \quad \hbox {for all} \ v_{h} \in V^{c}_{h}, \end{aligned}$$

(3.14)

where,

$$\begin{aligned} A^{c}_{h}(u_h,v_h) = a(u_h,v_h)+S_{h}^c(u_h,v_h). \end{aligned}$$

(3.15)

Further, we introduce a local projection norm for $v_{h} \in V^{c}_{h}$ as:

$$\begin{aligned} \left| \! \left| \! \left| v_{h}\right| \! \right| \! \right| ^{2}_{LP} = \alpha \left\| v_{h} \right\| ^{2} + \int _{\partial \varOmega } \frac{|\textbf{b} \cdot \textbf{n}|}{2} v^{2}_h \;\textit{ds}+S_{h}^{c}(v_h,v_h), \end{aligned}$$

(3.16)

and the local projection streamline derivative norm for $v_{h} \in V^{c}_{h}$ is defined as:

$$\begin{aligned} \left| \! \left| \! \left| v_{h}\right| \! \right| \! \right| _{LPSD}^{2} = \left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b} \cdot \nabla {v}_h) \right\| ^{2} + \left| \! \left| \! \left| v_{h}\right| \! \right| \! \right| ^{2}_{LP}. \end{aligned}$$

(3.17)

Further, the $L^2$-orthogonal projection ${J}^{c}_h: L^2 (\varOmega ) \rightarrow \text {V}^{c}_{h}$ satisfies the following approximation properties [2].

Lemma 5

$L^2$-orthogonal projections: [18, Lemma 7.1] The $L^2$-projection $J^{c}_h:L^2(\varOmega )\rightarrow \text {V}^{c}_{h}$ satisfies

$$\begin{aligned} \left\| h_{\mathscr {T}}^{-1} ({v-J^{c}_hv}) \right\| +\left\| \nabla (v-J^{c}_hv) \right\|&\le C\left\| h_{\mathscr {T}} v \right\| _{2}\qquad \forall ~v\in H^2(\varOmega ),\end{aligned}$$

(3.18)

$$\begin{aligned} \left( \sum _{E\in \mathscr {E}_h}\left\| {v}-{J}^{c}_h{v} \right\| ^2_{L^2(E)} \right) ^{1/2}&\le C\left\| h_{\mathscr {T}}^{3/2}{v} \right\| _{2}\qquad \forall ~{v}\in H^2(\varOmega ), \end{aligned}$$

(3.19)

$$\begin{aligned} ({v}-{J}^{c}_h {v},{v}_h)&=0 \qquad \forall ~{v}_h\in {V}^{c}_h . \end{aligned}$$

(3.20)

Further, the $L^{2}$-orthogonal projection operator exhibits the following $H^{1}$ and $L^{2}$-stability properties [2, Theorem 4.1,4.2]:

$$\begin{aligned} \left\| {J}^{c}_h{v} \right\| \le \left\| v \right\| , \quad \left\| h_{\mathscr {T}}^{-1}{J}^{c}_h{v} \right\| \le C \left\| h_{\mathscr {T}}^{-1}v \right\| , \quad \left\| \nabla {J}^{c}_h{v} \right\| \le C\left\| \nabla {v} \right\| . \end{aligned}$$

(3.21)

Moreover, the main result of this subsection is the following theorem, which ensures that the discrete bilinear form is well–posed [20].

Theorem 1

(Stability) Assume $\beta _a=\beta h_a$ for some $\beta |_{{\mathscr {M}}_a} =\frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_a)}+C}$. Then, the discrete bilinear form (3.15) satisfies the following inf-sup condition for some positive constant $\gamma$, independent of h:

$$\begin{aligned} \inf _{u_{h} \in V^{c}_{h} } \sup _{z_{h} \in V^{c}_{h} \setminus \{0\} }\frac{A^{c}_{h}(u_{h},z_{h})}{\left| \! \left| \! \left| u_{h}\right| \! \right| \! \right| _{LPSD}\left| \! \left| \! \left| z_{h}\right| \! \right| \! \right| _{LPSD}} \ge \gamma >0. \end{aligned}$$

Proof

Using Lemma 4, the inf-sup condition can be proved in three steps:

• Step 1: $A^c_{h}({u}_h,{u}_h) \ge c_0 \left| \! \left| \! \left| {u}_h\right| \! \right| \! \right| ^2_{LP},$

• Step 2: $\forall {u}_h \in {V}^c_h, \exists {v}_h \in {V}^c_h$: $A^c_{h}({u}_h,{v}_h) \ge c_0 \left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b} \cdot \nabla {u}_h) \right\| ^{2}- c_1 \left| \! \left| \! \left| {u}_h\right| \! \right| \! \right| ^2_{LP},$

• Step 3: $\left| \! \left| \! \left| {v}_h\right| \! \right| \! \right| _{LPSD} \le \left| \! \left| \! \left| {u}_h\right| \! \right| \! \right| _{LPSD}.$

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,

$$\begin{aligned} A^{c}_{h}(u_{h},u_{h})&\ge \alpha \left\| u_h \right\| ^{2} + \int _{\partial \varOmega } \frac{|\textbf{b} \cdot \textbf{n}|}{2} u^{2}_h \;\textit{ds}+ S_{h}^c(u_{h},u_{h}) = \left| \! \left| \! \left| u_h\right| \! \right| \! \right| ^{2}_{LP}. \end{aligned}$$

(3.22)

Further, the control of $\left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b} \cdot \nabla {v}_h) \right\| ^{2}$ can be obtained by choosing $v_{h} = J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}))$ in (3.15), that is,

$$\begin{aligned}&A^{c}_{h}(u_{h},J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}))) \nonumber \\&\quad =\left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b}\cdot \nabla {u_h}) \right\| ^{2} + \left( \textbf{b}\cdot \nabla {u_{h}}, J^{c}_{h}\big (h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}})\big ) - h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}})\right) \nonumber \\&\qquad + \left( \mu {u_{h}}, J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}))\right) \nonumber \\&\qquad + \int _{\partial \varOmega } ({\textbf{b} \cdot \textbf{n}})^{\circleddash }u_h J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})) \;\textit{ds}+S_{h}^c(u_{h},J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}))) \nonumber \\&\quad = \left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b}\cdot \nabla {u_h}) \right\| ^{2} +(a)+(b)+(c) +(d). \end{aligned}$$

(3.23)

Let us now estimate these four terms. Using the canonical representation of the basis function $\phi _a$ at the node $a \in \mathscr {V}_h$ for the mesh $\mathscr {T}_h$ i.e. $\sum _{a \in \mathscr {V}_{h}}{\phi _{a}} = 1$,

$$\begin{aligned} (a)&=\sum _{K \in \mathscr {T}_{h}}\int _{K}\left( \sum _{a\in \mathscr {V}_h} \phi _a\right) (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}})) - h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}})) (\textbf{b}\cdot \nabla {u_{h}})\;\textit{dx}\\&=\sum _{a\in \mathscr {V}_h} \int _{\mathscr {M}_a} (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}})) - h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}})) (\textbf{b}\cdot \nabla {u_{h}}) \phi _a\;\textit{dx}. \end{aligned}$$

Using the orthogonality property of $L^2$-projection (3.20) with the test function $C_a \phi _a \in V^c_h$, where $C_a$ is a constant and $\left\| \phi _{a} \right\| _{\infty } \le 1$,

$$\begin{aligned} (a)&\le \sum _{a\in \mathscr {V}_h} \left\| J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}))- h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}) \right\| _{L^2(\mathscr {M}_a)} \left\| \textbf{b}\cdot \nabla {u_{h}} -C_a \right\| _{L^2(\mathscr {M}_a)}. \end{aligned}$$

Choosing the constant $C_a = \frac{1}{\vert {\mathscr {M}_{a}\vert }} \int _{\mathscr {M}_{a}} {\textbf{b}\cdot \nabla {u_h}} \;\textit{dx}$, and applying the Cauchy–Schwarz inequality, (3.21) and Young’s inequality:

$$\begin{aligned} (a) \le&\left( \sum _{a\in \mathscr {V}_h} \beta ^{-1}_{a} \left\| J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}))- h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}) \right\| ^2_{L^2(\mathscr {M}_a)} \right) ^{1/2}\\&\left( \sum _{a\in \mathscr {V}_h} \beta _{a} \int _{\mathscr {M}_a}\kappa ^{2}_a ({u_{h}})\;\textit{dx}\right) ^{1/2}\\ \le&\left\| \textbf{b} \right\| ^{\frac{1}{2}}_{{\infty }} \left( \sum _{a\in \mathscr {V}_h} \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}) \right\| ^2_{L^2(\mathscr {M}_a)} \right) ^{1/2} [S_{h}^{c}(u_h, u_h)]^{\frac{1}{2}} \\ \le&\left\| \textbf{b} \right\| ^{\frac{1}{2}}_{{\infty }} \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}) \right\| [S_{h}^{c}(u_h, u_h)]^{\frac{1}{2}} \\ \le&C S_{h}^c(u_h, u_h) + \frac{1}{6}\left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}}) \right\| ^2, \end{aligned}$$

the constant C in the above estimate depends on $\left\| \textbf{b} \right\| ^{\frac{1}{2}}_{{\infty }}$. The second term is estimated by applying the Cauchy–Schwarz inequality followed by (3.21) and an inverse inequality,

$$\begin{aligned} (b) \le \left\| \mu {u_{h}} \right\| \left\| J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}) \right\| \le C \alpha \left\| u_h \right\| ^{2}. \end{aligned}$$

(3.24)

The constant C in (3.24) depends on $\left\| \textbf{b} \right\| _{\infty }$ and $\alpha ^{-1}$. The third term is handled by applying the Cauchy–Schwarz inequality, trace inequality (2.6), (3.21) and Young’s inequality,

$$\begin{aligned} (c)&\le \sum _{E \in \mathscr {E}^{B}_h} \left\| ({\textbf{b} \cdot \textbf{n}})^{\circleddash }u_h \right\| _{L^{2}(E)} \left\| J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}})) \right\| _{L^{2}(E)} \\ {}&\le C \sum _{E \in \mathscr {E}^{B}_h} \int _{E} \frac{|\textbf{b} \cdot \textbf{n}|}{2} u^{2}_h \;\textit{ds}+ \frac{1}{6}\left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}}) \right\| ^2. \end{aligned}$$

The constant C depends on $\left\| \textbf{b} \right\| _{\infty }$. Next, applying the Cauchy–Schwarz inequality to the fourth term to obtain,

$$\begin{aligned} (d) \le [S_{h}^{c}(u_{h},u_{h})]^{\frac{1}{2}}[S_{h}^{c}\big (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}})),J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}))\big )]^{\frac{1}{2}}. \end{aligned}$$

(3.25)

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 $\beta _{a}=\beta h_{a}$ with $\beta = 1/( \left\| \textbf{b} \right\| _{W^{1,\infty }({\mathscr {M}}_a)}+C)$,

$$\begin{aligned} S_{h}^c&(J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}})),J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}))) \\&\le \sum _{a\in \mathscr {V}_h} \beta _{a} \left\| {\textbf{b} \cdot \nabla (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}})))}-\frac{1}{|\mathscr {M}_a|} \int _{\mathscr {M}_a}{\textbf{b} \cdot \nabla (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}))} \;\textit{dx} \right\| _{L^{2}(\mathscr {M}_a)}^{2}\\&\le C \sum _{a\in \mathscr {V}_h} \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}) \right\| _{L^{2}(\mathscr {M}_a)}^2. \end{aligned}$$

The constant C in the above estimates depends on $\left\| \textbf{b} \right\| _{{\infty }}$. Then, it follows that,

$$\begin{aligned} S_{h}^c(u_{h},J^{c}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}))) \le C S_{h}^c(u_h, u_h) + \frac{1}{6} \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla {u_{h}}) \right\| ^2. \end{aligned}$$

(3.26)

The constant C in the above estimates depends on $\left\| \textbf{b} \right\| _{{\infty }}$. Put together, (3.23) leads to,

$$\begin{aligned} A^{c}_{h}(u_{h},J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})))&\ge \frac{1}{2}\left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b} \cdot \nabla {u}_h) \right\| ^{2} - {C} \left| \! \left| \! \left| u_h\right| \! \right| \! \right| _{LP}^{2}. \end{aligned}$$

(3.27)

The constant C in the estimates (3.27) depends on $\left\| \textbf{b} \right\| _{{\infty }}$. Note that

$$\begin{aligned} \left| \! \left| \! \left| J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}}))\right| \! \right| \! \right| _{LPSD}&= \alpha \left\| J^{c}_{h}( h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})) \right\| ^{2} \nonumber \\&\quad +\sum _{E \in \mathscr {E}_{h}^{B} } \int _{E} \frac{|\textbf{b} \cdot \textbf{n}|}{2} (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})))^{2} \;\textit{ds}\nonumber \\ {}&\quad +S_{h}^c (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})),J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}}))) \nonumber \\&\quad + \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b} \cdot \nabla (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})))) \right\| ^2. \end{aligned}$$

(3.28)

Now, estimate the four terms of (3.28) individually. Using the stability of the projection operator (3.21) and an inverse inequality,

$$\begin{aligned} \alpha \left\| J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})) \right\| ^{2}&\le \alpha C\left\| \textbf{b} \right\| ^2_{{\infty }}\left\| u_h \right\| ^{2} \le C\left\| \textbf{b} \right\| ^2_{{\infty }} \left| \! \left| \! \left| u_h\right| \! \right| \! \right| ^{2}_{LPSD}. \end{aligned}$$

The second term is estimated by using trace inequality and (3.21) ,

$$\begin{aligned}&\sum _{E \in \mathscr {E}_{h}^{B} } \int _{E} \frac{|\textbf{b} \cdot \textbf{n}|}{2} (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})))^{2} \;\textit{ds}\\&\quad \le C\left\| \textbf{b} \right\| _{{\infty }} \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}}) \right\| ^{2} \le C \left\| \textbf{b} \right\| _{{\infty }} \left| \! \left| \! \left| u_h\right| \! \right| \! \right| ^{2}_{LPSD}. \end{aligned}$$

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,

$$\begin{aligned} S_{h}^c (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})),J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})))&+ \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b} \cdot \nabla (J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}})))) \right\| ^2 \\&\le C \left\| \textbf{b} \right\| _{{\infty }} \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}}) \right\| ^{2} \\&\le \left\| \textbf{b} \right\| _{{\infty }} \left| \! \left| \! \left| u_h\right| \! \right| \! \right| ^{2}_{LPSD}. \end{aligned}$$

Finally, put together, we get,

$$\begin{aligned} \left| \! \left| \! \left| J^{c}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla {u_{h}}))\right| \! \right| \! \right| _{LPSD} \le C \left| \! \left| \! \left| {u_{h}}\right| \! \right| \! \right| _{LPSD}. \end{aligned}$$

(3.29)

The constant C in (3.29) depends on $\left\| \textbf{b} \right\| _{\infty }$. The selection of $z_h$ is

$$\begin{aligned} z_h = u_h+ \frac{1}{C+1} J^{c}_{h}(h_{\mathscr {T}} (\textbf{b} \cdot \nabla {u_{h}})), \end{aligned}$$

where $J^{c}_h$ is as defined in Lemma 5. Finally, the result follows by Lemma 4. $\square$

3.2 A priori error estimates

Lemma 6

Suppose $u \in H^{2}(\varOmega )$ and $\beta _a=\beta h_a$ for some $\beta |_{{\mathscr {M}}_a} =\frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_a)}+C}$. Then,

$$\begin{aligned} \left| \! \left| \! \left| u-J^{c}_h u\right| \! \right| \! \right| _{LPSD}\le C\left\| h_{\mathscr {T}}^{3/2}u \right\| _2. \end{aligned}$$

Proof

Consider the terms in the LPSD norm defined in (3.17),

$$\begin{aligned}&\left| \! \left| \! \left| u-J^{c}_h u\right| \! \right| \! \right| ^{2}_{LPSD}=\left\| u-J^{c}_h{u} \right\| ^{2}+ \left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b}\cdot \nabla ({u-J^{c}_{h}{u}})) \right\| ^{2}\nonumber \\&\quad + \int _{\partial \varOmega } \frac{|\textbf{b} \cdot \textbf{n}|}{2} (u-J^{c}_{h}{u})^{2}\;\textit{ds}+S_{h}^c(u-J^{c}_h{u}, u-J^{c}_h{u}). \end{aligned}$$

(3.30)

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

$\Vert u-J^{c}_h{u} \Vert \le \left\| h_{\mathscr {T}}^{2} u \right\| _{2} \text {and} \ \left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b}\cdot \nabla ({u-J^{c}_{h}{u}})) \right\| \le C \left\| h_{\mathscr {T}}^{\frac{3}{2}} u \right\| _{2}. \$

The constant C depends on $\left\| \textbf{b} \right\| _{\infty }$. The third term of (3.30) is handled by using the trace inequality (3.19) over each edge,

$$\begin{aligned} \int _{\partial \varOmega } \frac{|\textbf{b} \cdot \textbf{n}|}{2} (u-J^{c}_{h}{u})^{2}\;\textit{ds}\le C\left\| h_{\mathscr {T}}^{\frac{3}{2}}u \right\| ^{2}_{2}. \end{aligned}$$

Note that the constant C in the above estimates depends on the $\left\| \textbf{b} \right\| _{\infty }$. The last term is estimated by using boundedness of the local projection operator and $\beta _{a} = \beta h_{a}$ with $\beta = 1/({\left\| \textbf{b} \right\| _{W^{1,\infty }(\mathscr {M}_{a})}+C})$,

$$\begin{aligned} S_{h}^c(u-J^{c}_h{u}, u-J^{c}_h{u})&:= \sum _{a \in \mathscr {V}_{h}} \beta _{a} \left\| {\textbf{b}\cdot \nabla ({u-J^{c}_{h}{u}})}- \frac{1}{|\mathscr {M}_a|}\int _{\mathscr {M}_a}{\textbf{b}\cdot \nabla ({u-J^{c}_{h}{u}})\;\textit{dx}} \right\| ^{2}_{L^{2}(\mathscr {M}_a)} \\&\le \sum _{a \in \mathscr {V}_{h}} \beta h_{a} \left\| {\textbf{b}\cdot \nabla ({u-J^{c}_{h}{u}})} \right\| _{L^{2}(\mathscr {M}_a)}^{2} \le C \left\| h_{\mathscr {T}}^{\frac{1}{2}}{\nabla ({u-J^{c}_{h}{u}})} \right\| ^{2} \\ {}&\le C \left\| h_{\mathscr {T}}^{3/2}u \right\| ^{2}_{2}. \end{aligned}$$

The constant C depends on $\left\| \textbf{b} \right\| _{\infty }$. The combination of the above estimates concludes the proof. $\square$

Lemma 7

Suppose $u \in H^{2}(\varOmega )$ and $\beta _a = \beta {h_a}$ for some $\beta |_{{\mathscr {M}}_a} =\frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_a)}+C}$. Then,

$$\begin{aligned} A^{c}_{h}(u-J^{c}_h{u},{v}_h) \le Ch^{\frac{3}{2}} \left\| u \right\| _{2} \left| \! \left| \! \left| v_h\right| \! \right| \! \right| _{LPSD} \quad \forall ~v_h \in V^{c}_{h}. \end{aligned}$$

(3.31)

Proof

Applying an integration by parts to the first term of the discrete bilinear form in (3.15) to obtain,

$$\begin{aligned} A^{c}_{h}({u-J^{c}_{h}{u}},v_{h}) =&-(u-J^{c}_{h}{u}, \textbf{b}\cdot \nabla {v_h}) +((\mu - \text{ div } \textbf{b})(u-J^{c}_h{u}),{v}_h)+S_{h}^c({u-J^{c}_{h}{u}},v_{h}) \\ {}&+\int _{\partial \varOmega } (\textbf{b} \cdot \textbf{n}) (u-J^{c}_{h}{u})v_h \;\textit{ds}+ \int _{\partial \varOmega } (\textbf{b} \cdot \textbf{n})^{\circleddash } (u-J^{c}_{h}{u})v_h\;\textit{ds}\\ =&-(u-J^{c}_{h}{u}, \textbf{b}\cdot \nabla {v_h}) +((\mu - \text{ div } \textbf{b})(u-J^{c}_h{u}),{v}_h)+S_{h}^c({u-J^{c}_{h}{u}},v_{h}) \\ \ {}&+ \int _{\partial \varOmega } (\textbf{b} \cdot \textbf{n})^{\oplus } (u-J^{c}_{h}{u})v_h \;\textit{ds}=(a)+(b)+(c)+(d). \end{aligned}$$

The first term is estimated by using the Cauchy–Schwarz inequality and $L^2$-projection property (3.18),

$$\begin{aligned} (a) \le \left\| u-J^{c}_{h}{u} \right\| \left\| \textbf{b}\cdot \nabla {v_h} \right\|&\le \left\| h_{\mathscr {T}}^{2}u \right\| _{2} \left\| (\textbf{b}\cdot \nabla {v_h}) \right\| \le \left\| h_{\mathscr {T}}^{\frac{3}{2}}u \right\| _{2} \left\| h_{\mathscr {T}}^{\frac{1}{2}}(\textbf{b}\cdot \nabla {v_h}) \right\| \\&\le C \left\| h_{\mathscr {T}}^{3/2}u \right\| _2\left| \! \left| \! \left| v_h\right| \! \right| \! \right| _{LPSD}. \end{aligned}$$

and,

$$\begin{aligned} (b)&\le \frac{\left\| \mu -\text{ div } \textbf{b} \right\| _{\infty }}{\surd {\alpha }} \left\| u-J^{c}_h{u} \right\| \surd {\alpha } \left\| v_h \right\| \le C \left\| h_{\mathscr {T}}^{2}u \right\| _{2} \left| \! \left| \! \left| v_h\right| \! \right| \! \right| _{LPSD}. \end{aligned}$$

(3.32)

The constant C in (3.32) depends on $\alpha ^{-1/2}$. The third term is handled by applying the Cauchy–Schwarz inequality, boundedness of the local projection, approximation estimates (3.18) and $\beta _{a} = \beta h_{a}$ with $\beta = 1/({\left\| \textbf{b} \right\| _{W^{1,\infty }(\mathscr {M}_{a})}+C})$,

$$\begin{aligned} (c)&=\sum _{a\in \mathscr {V}_h} \beta _a\big ( \kappa _a(u-J^{c}_hu), \kappa _a( v_h)\big )_{L^{2}(\mathscr {M}_a)} \\&\quad \le \left( \sum _{a\in \mathscr {V}_h} \beta _a\left\| \kappa _a(u-J^{c}_hu) \right\| _{L^2(\mathscr {M}_a)}^2 \right) ^{1/2} \left| \! \left| \! \left| v_h\right| \! \right| \! \right| _{LPSD} \\&\quad \le \left( \sum _{a\in \mathscr {V}_h} \beta _a \left\| \textbf{b}\cdot \nabla (u-J^{c}_hu) \right\| ^2_{L^2(\mathscr {M}_a)} \right) ^{1/2} \left| \! \left| \! \left| v_h\right| \! \right| \! \right| _{LPSD} \\ {}&\le C\left\| h_{\mathscr {T}}^{3/2}u \right\| _2\left| \! \left| \! \left| v_h\right| \! \right| \! \right| _{LPSD}. \end{aligned}$$

Applying the Cauchy–Schwarz inequality, trace inequality (2.5) and approximation estimates (3.18),

$$\begin{aligned} (d) \le C\left( \sum _{E \in \mathscr {E}^{B}_{h}} \left\| u-J^{c}_{h}{u} \right\| ^{2}_{L^{2}(E)} \right) ^{\frac{1}{2}} \left( \int _{\partial \varOmega }{\frac{|\textbf{b} \cdot \textbf{n}|}{2}} v^{2}_h \;\textit{ds}\right) ^{\frac{1}{2}} \le C \left\| h_{\mathscr {T}}^{3/2}u \right\| _2\left| \! \left| \! \left| v_h\right| \! \right| \! \right| _{LPSD}. \end{aligned}$$

Combining the above estimates leads to (3.31) and it concludes the proof. $\square$

Theorem 2

Let $u\in H^2(\varOmega )$ be the solution of (2.3) and $u_h\in V^c_h$ be the discrete solution of (3.14). Let $\beta _a=\beta h_a$ for some $\beta |_{{\mathscr {M}}_a} =\frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_a)}+C}$. Then,

$$\begin{aligned} \left| \! \left| \! \left| u-u_h\right| \! \right| \! \right| _{LPSD}\le C \left\| h_{\mathscr {T}}^{3/2}u \right\| _2. \end{aligned}$$

Proof

By adding and subtracting the interpolation operator $J^{c}_h u$, we decompose the error as follows:

$$\begin{aligned} \left| \! \left| \! \left| u-u_h\right| \! \right| \! \right| _{LPSD}\le \left| \! \left| \! \left| u-J^{c}_h u\right| \! \right| \! \right| _{LPSD}+\left| \! \left| \! \left| J^{c}_h u-u_h\right| \! \right| \! \right| _{LPSD}. \end{aligned}$$

(3.33)

In the second term of (3.33), using the estimate of Theorem 1,

$$\begin{aligned}&c\left| \! \left| \! \left| u_h-J^{c}_h u\right| \! \right| \! \right| _{LPSD} \le \text{ sup}_{w_h \in V^{c}_{h}}{\frac{A^{c}_h( u_h-J^{c}_h u,w_h)}{\left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{LPSD}}}\nonumber \\&\quad =\text{ sup}_{w_h \in V^{c}_{h}}{\frac{A^{c}_h( u_h-u,w_h)+ A^{c}_h(u-J^{c}_h u,w_h)}{\left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{LPSD}}}. \end{aligned}$$

(3.34)

The weak formulation (2.4) and (3.15) imply that

$$\begin{aligned} A^{c}_h( u_h-u,w_h)&= -S_{h}^c(u,w_h). \end{aligned}$$

Moreover, the Cauchy–Schwarz inequality implies

$$\begin{aligned} S_{h}^c(u,w_h)&=\sum _{a\in \mathscr {V}_h} \beta _a \left( \kappa _a( u), \kappa _a( w_h)\right) _{L^{2}(\mathscr {M}_a)} \\&\le \left( \sum _{a\in \mathscr {V}_h} \beta _a\left\| \kappa _a( u) \right\| _{L^2(\mathscr {M}_a)}^2 \right) ^{1/2} \left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{LPSD}. \end{aligned}$$

Note that $\beta _{a} = \beta h_{a}$ with $\beta |_{{\mathscr {M}}_a} =\frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_a)}+C}$. Using the Poincaré inequality (2.9) for every vertex $a\in \mathscr {V}_h$,

$$\begin{aligned}&S_{h}^c(u,w_h) \le C \Big (\sum _{a\in \mathscr {V}_h} \beta _a h_a^2\left\| \nabla (\textbf{b}\cdot \nabla u) \right\| ^2_{L^2(\mathscr {M}_a)} \Big )^{1/2} \left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{LPSD} \\&\quad \le C \left( \sum _{a\in \mathscr {V}_h} \beta \left( \left\| \textbf{b} \right\| ^{2}_{W^{1,{\infty }}({\mathscr {M}}_a)}\left\| h_{a}^{3/2}u \right\| ^{2}_{H^{2}(\mathscr {M}_a)}+\left\| \nabla \textbf{b} \right\| ^{2}_{W^{{1,\infty }}({\mathscr {M}}_a)} \left\| h_{a}^{3/2}u \right\| ^{2}_{H^{1}(\mathscr {M}_a)}\right) \right) ^{\frac{1}{2}} \left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{LPSD} \\&\quad \le \left( \sum _{a\in \mathscr {V}_h} \beta \left( \left\| \textbf{b} \right\| ^{2}_{W^{1,{\infty }}({\mathscr {M}}_a)}+\left\| \nabla \textbf{b} \right\| ^{2}_{W^{{1,\infty }}({\mathscr {M}}_a)} \right) \left\| h_{a}^{3/2}u \right\| ^{2}_{H^{2}(\mathscr {M}_a)} \right) ^{\frac{1}{2}} \left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{LPSD}\\&\quad \le C\left\| h_{\mathscr {T}}^{3/2}u \right\| _2\left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{LPSD}. \end{aligned}$$

It follows that,

$$\begin{aligned} A^{c}_h( u_h-u,w_h)\le C\left\| h_{\mathscr {T}}^{3/2}u \right\| _2\left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{LPSD}. \end{aligned}$$

(3.35)

Using the estimate (3.35) and Lemma 7 in (3.34), we obtain

$$\begin{aligned} \left| \! \left| \! \left| u_h-J^{c}_h u\right| \! \right| \! \right| _{LPSD}\le \left\| h_{\mathscr {T}}^{3/2}u \right\| _2. \end{aligned}$$

(3.36)

Finally, Lemma 6 and (3.36) lead (3.33) to an a priori estimate. $\square$

4 Nonconforming finite element discretization

The nonconforming discrete solution of (2.3) is a function $u_h \in V^{nc}_h$ such that

$$\begin{aligned} a^{nc}_{h}(u_h,v)=l(v) \qquad \forall ~ v \in \ V^{nc}_h, \end{aligned}$$

(4.37)

where,

$$\begin{aligned} a^{nc}_{h}(u_h,v) :&= (\textbf{b}\cdot \nabla _h{u_h}, v) + (\mu {u_h}, v) +\int _{\partial \varOmega } (\textbf{b} \cdot \textbf{n})^{\circleddash } uv \;\textit{ds}\\&\quad -\sum _{E\in \mathscr {E}_h^I}\int _E (\textbf{b}\cdot \textbf{n})[\hspace{-0.5pt}u_h \hspace{-0.5pt}]\{\hspace{-1.0pt}v \hspace{-1.0pt}\} \;\textit{ds}+ \sum _{E\in \mathscr {E}_h^I}\int _E \frac{|\textbf{b}\cdot \textbf{n}|}{2}[\hspace{-0.5pt}u_h \hspace{-0.5pt}][\hspace{-0.5pt}v \hspace{-0.5pt}] \;\textit{ds}. \end{aligned}$$

Here, $\nabla _h$ denotes the piecewise (element-wise) gradient operator. For each $E \in \mathscr {E}_{h}$, define the fluctuation operator $\kappa _{E}: V +V^{nc}_h\rightarrow L^{2}(\mathscr {M}_{E})$ such that

$$\begin{aligned} \kappa _{E}(u_h):= \textbf{b}\cdot \nabla _h{u_h} \ - \frac{1}{\vert {\mathscr {M}_{E}\vert }} \int _{\mathscr {M}_{E}} {\textbf{b}\cdot \nabla _h{u_h}} \;\textit{dx}, \end{aligned}$$

where, $\vert {\mathscr {M}_{E}}\vert$ denotes the measure of $\mathscr {M}_{E}$. For each $E \in \mathscr {E}_h$, let $\beta _{E}:=\beta h_{E}$ for stabilization parameter $\beta = 1/ (\left\| \textbf{b} \right\| _{W^{1,\infty }({\mathscr {M}}_E)}+C)$, 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

$$\begin{aligned} S_{h}^{nc}(u_h, v_h):= \sum _{E \in \mathscr {E}_{h}} \beta _{E} \big ({\kappa _{E}({u_h})}, {\kappa _{E}({v_h})\big )}_{L^{2}(\mathscr {M}_{E})}, \end{aligned}$$

(4.38)

using this term, the nonconforming generalized local projection stabilized discrete form of (2.3) reads as:

Find $u_{h} \in V^{nc}_{h}$ such that

$$\begin{aligned} A^{nc}_{h}(u_{h},v_{h}) = l(v_{h}) \qquad \forall ~ v_{h} \in V^{nc}_{h}, \end{aligned}$$

(4.39)

where,

$$\begin{aligned} \begin{array} {rcl} A^{nc}_{h}(u_h,v_h)= & {} a^{nc}_{h}(u_h,v_h)+S_{h}^{nc}(u_h,v_h). \end{array} \end{aligned}$$

(4.40)

Further, we define a nonconforming local projection norm by

$$\begin{aligned} \left| \! \left| \! \left| v_{h}\right| \! \right| \! \right| ^{2}_{NLP} = \alpha \left\| v_{h} \right\| ^{2} + \int _{\partial \varOmega } \frac{| \textbf{b} \cdot \textbf{n}|}{2} v^{2}_h \;\textit{ds}+S_{h}^{nc}(v_h,v_h) + \sum _{E\in \mathscr {E}_h^I}\int _E \frac{|\textbf{b}\cdot \textbf{n}|}{2}[\hspace{-0.5pt}v_h \hspace{-0.5pt}]^2 \;\textit{ds}, \nonumber \\ \end{aligned}$$

(4.41)

and nonconforming local projection streamline derivative norm by

$$\begin{aligned} \left| \! \left| \! \left| v_{h}\right| \! \right| \! \right| _{NLPSD}^{2} = \left| \! \left| \! \left| v_{h}\right| \! \right| \! \right| ^{2}_{NLP} +\left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b} \cdot \nabla _h{v_h}) \right\| ^{2}, \end{aligned}$$

(4.42)

for all $v_{h} \in V^{nc}_{h}$. Moreover, the $L^2$-projection ${J}^{nc}_h: L^2 (\varOmega ) \rightarrow \text {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].

Theorem 3

(Stability) Assume $\beta _E=\beta h_E$ for some $\beta |_{{\mathscr {M}}_E} =\frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_E)}+C}$. Then, the discrete bilinear form (4.39) satisfies the following inf-sup condition for a positive constant $\nu$, independent of h:

$$\begin{aligned} \inf _{u_{h} \in V^{nc}_{h} } \sup _{z_{h} \in V^{nc}_{h} \setminus \{0\} }\frac{A^{nc}_{h}(u_{h},z_{h})}{\left| \! \left| \! \left| u_{h}\right| \! \right| \! \right| _{NLPSD}\left| \! \left| \! \left| z_{h}\right| \! \right| \! \right| _{NLPSD}} \ge \nu . \end{aligned}$$

(4.43)

Proof

Using Lemma 4, the inf-sup condition can be proved in three steps:

• Step 1: $A^{nc}_{h}({u}_h,{u}_h) \ge c_0 \left| \! \left| \! \left| {u}_h\right| \! \right| \! \right| ^2_{NLP},$

• Step 2: $\forall {u}_h \in {V}^{nc}_h, \exists {v}_h \in {V}^{nc}_h$: $A^{nc}_{h}({u}_h,{v}_h) \ge c_0 \left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b} \cdot \nabla _h{u_h}) \right\| ^{2}- c_1 \left| \! \left| \! \left| {u}_h\right| \! \right| \! \right| ^2_{NLP},$

• Step 3: $\left| \! \left| \! \left| {v}_h\right| \! \right| \! \right| _{NLPSD} \le \left| \! \left| \! \left| {u}_h\right| \! \right| \! \right| _{NLPSD}.$

The key steps to derive the above estimates are as follows: Choosing first $v_h=u_h$ as a test function in (4.40),

$$\begin{aligned} A^{nc}_{h}(u_{h},&u_{h}) \ge \alpha \left\| u_h \right\| ^{2} + \int _{\partial \varOmega } \frac{|\textbf{b} \cdot \textbf{n}|}{2} u^{2}_h \;\textit{ds}+\sum _{E\in \mathscr {E}_h^I}\int _E (\textbf{b}\cdot \textbf{n}) [\hspace{-0.5pt}u_h \hspace{-0.5pt}]^{2} \;\textit{ds}+ S_{h}^{nc}(u_{h},u_{h}). \end{aligned}$$

Further, the control of $\left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b} \cdot \nabla _h{v_h}) \right\| ^{2}$ is obtained by choosing $v_{h} = J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}))$ in (4.40), we have by adding and subtracting $\left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}) \right\| ^{2}$,

$$\begin{aligned} A^{nc}_{h}\big (u_{h},J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})\big )&= \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}) \right\| ^{2} \nonumber \\&\quad + (\textbf{b}\cdot \nabla _h{u_{h}}, J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) - h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \nonumber \\ {}&\quad + \big (\mu {u_{h}}, J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}))\big )\nonumber \\&\quad + \int _{\partial \varOmega } ({\textbf{b} \cdot \textbf{n}})^{\circleddash }u_h J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}))\;\textit{ds}\nonumber \\&\quad +S_{h}^{nc}(u_{h},J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \nonumber \\&\quad +\sum _{E\in \mathscr {E}_h^I}\int _E \frac{|\textbf{b}\cdot \textbf{n}|}{2} [\hspace{-0.5pt}u_h \hspace{-0.5pt}][\hspace{-0.5pt}J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \hspace{-0.5pt}]\;\textit{ds}\nonumber \\&\quad -\sum _{E\in \mathscr {E}_h^I}\int _E (\textbf{b}\cdot \textbf{n}) [\hspace{-0.5pt}u_h \hspace{-0.5pt}]\{\hspace{-1.0pt}J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \hspace{-1.0pt}\} \;\textit{ds}. \end{aligned}$$

(4.44)

Most of the estimates of (4.44) can be derived in a similar way as shown in (3.23).

$$\begin{aligned}&(\textbf{b}\cdot \nabla _h{u_{h}},J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}))-h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \\&\quad \le C S_{h}^{nc}(u_h, u_h) + \frac{1}{10}\left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b} \cdot \nabla _h{u_{h}}) \right\| , \\&\quad (\mu {u_{h}}, J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}))) \le C \alpha \left\| u_h \right\| ^{2}, \\&\quad S_{h}^{nc}(u_{h},J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \\ {}&\quad \le C S_{h}^{nc}(u_h, u_h) + \frac{1}{10} \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}) \right\| ^2, \\&\quad \int _{\partial \varOmega } ({\textbf{b} \cdot \textbf{n}})^{\circleddash }u_h J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}))\;\textit{ds}\\ {}&\quad \le C \int _{\partial \varOmega } \frac{|\textbf{b} \cdot \textbf{n}|}{2} u^{2}_h \;\textit{ds}+ \frac{1}{10}\left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b} \cdot \nabla _h{u_{h}}) \right\| ^2. \end{aligned}$$

The constant C in the above estimates depends on $\left\| \textbf{b} \right\| _{{\infty }}$. Now, it is sufficient to estimate the last two terms of (4.44). Using the Cauchy–Schwarz inequality,

$$\begin{aligned}&\sum _{E\in \mathscr {E}_h^I}\int _E \frac{|\textbf{b}\cdot \textbf{n}|}{2} [\hspace{-0.5pt}u_h \hspace{-0.5pt}][\hspace{-0.5pt}J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \hspace{-0.5pt}]\;\textit{ds}\\&\quad \le C\left( \sum _{E \in \mathscr {E}^{I}_h} \int _{E} \frac{|\textbf{b} \cdot \textbf{n}|}{2} [\hspace{-0.5pt}u_h \hspace{-0.5pt}]^{2} \;\textit{ds}\right) ^{\frac{1}{2}} \left( \sum _{E \in \mathscr {E}^{I}_h} \left\| [\hspace{-0.5pt}J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \hspace{-0.5pt}] \right\| _{L^{2}(E)}^{2} \right) ^{\frac{1}{2}}. \end{aligned}$$

At the edge E, the jump term contributes to both the triangles sharing that edge; using the trace inequality (2.6) and (3.21),

$$\begin{aligned} \left\| [\hspace{-0.5pt}J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \hspace{-0.5pt}] \right\| _{L^2(E)}&\le C \left\| h_{\mathscr {T}}^{-1/2}J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \right\| _{L^2(\mathscr {M}_{E})} \\&\le C \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}) \right\| ^2_{L^2(\mathscr {M}_{E})}. \end{aligned}$$

We then get,

$$\begin{aligned}&\sum _{E\in \mathscr {E}_h^I}\int _E \frac{|\textbf{b}\cdot \textbf{n}|}{2} [\hspace{-0.5pt}u_h \hspace{-0.5pt}][\hspace{-0.5pt}J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \hspace{-0.5pt}]\;\textit{ds}\\&\quad \le C\left( \sum _{E \in \mathscr {E}^{I}_h} \int _{E} \frac{|\textbf{b} \cdot \textbf{n}|}{2} [\hspace{-0.5pt}u_h \hspace{-0.5pt}]^{2} \;\textit{ds}\right) ^{\frac{1}{2}} \\&\quad \left( \sum _{E \in {\mathscr {E}}_h} \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}) \right\| _{L^{2}(\mathscr {M}_{E})}^{2} \right) ^{\frac{1}{2}}\\ {}&\quad \le C \sum _{E \in \mathscr {E}^{I}_h} \int _{E} \frac{|\textbf{b} \cdot \textbf{n}|}{2} [\hspace{-0.5pt}u_h \hspace{-0.5pt}]^{2} \;\textit{ds}+\frac{1}{10} \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}) \right\| ^2. \end{aligned}$$

Similarly, the next term is estimated as:

$$\begin{aligned}&\sum _{E\in \mathscr {E}_h^I}\int _E (\textbf{b}\cdot \textbf{n}) [\hspace{-0.5pt}u_h \hspace{-0.5pt}]\{\hspace{-1.0pt}J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})) \hspace{-1.0pt}\}\;\textit{ds}\\&\quad \le C \sum _{E \in \mathscr {E}^{I}_h} \int _{E} \frac{|\textbf{b} \cdot \textbf{n}|}{2} [\hspace{-0.5pt}u_h \hspace{-0.5pt}]^{2} \;\textit{ds}+\frac{1}{10} \left\| h^{\frac{1}{2}}_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}}) \right\| ^2. \end{aligned}$$

Combining all these estimates and (4.44) lead to,

$$\begin{aligned} A^{nc}_{h}(u_{h},J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b}\cdot \nabla _h{u_{h}})))&\ge \frac{1}{2}\left\| h^{\frac{1}{2}}_{\mathscr {T}} (\textbf{b} \cdot \nabla _h{u}_h) \right\| ^{2} - {C} \left| \! \left| \! \left| u_h\right| \! \right| \! \right| _{NLP}^{2}. \end{aligned}$$

(4.45)

The constant C in the above estimates depends on $\left\| \textbf{b} \right\| _{{\infty }}$. In particular, the inequality holds for

$$\begin{aligned} z_h = u_h+\frac{1}{C+1} J^{nc}_{h}(h_{\mathscr {T}}(\textbf{b} \cdot \nabla _h{u_{h}})), \end{aligned}$$

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^{\infty }$-norm of $\textbf{b}$. Rest of the proof can be derived in a similar way as in the proof of (3.28)–(3.29). $\square$

4.1 A priori error estimates

Lemma 8

Suppose $u \in H^{2}(\varOmega )$ and $\beta _E=\beta h_E$ for some $\beta |_{{\mathscr {M}}_a} =\frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_a)}+C}$. Then,

$$\begin{aligned} \left| \! \left| \! \left| u-J^{nc}_h u\right| \! \right| \! \right| _{NLPSD}\le C \left\| h_{\mathscr {T}}^{3/2}u \right\| _2. \end{aligned}$$

(4.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),

$$\begin{aligned} \sum _{E\in \mathscr {E}_h^I}\int _E \frac{|\textbf{b}\cdot \textbf{n}|}{2}[\hspace{-0.5pt}u-J^{nc}_h u \hspace{-0.5pt}]^2\;\textit{ds}\le C \sum _{E\in \mathscr {E}_h^I} \left\| [\hspace{-0.5pt}u-J^{nc}_h u \hspace{-0.5pt}] \right\| ^2_{L^{2}(E)}. \end{aligned}$$

The constant C depends on $\left\| \textbf{b} \right\| _{\infty }$. At the edge E, the jump term contributes to both the triangles sharing that edge; using the trace inequality (2.5),

$$\begin{aligned} \left\| [\hspace{-0.5pt}u-J^{nc}_h u \hspace{-0.5pt}] \right\| _{L^2(E)}&\le C\left( h_{K}^{-1/2}\left\| u-J^{nc}_h u \right\| _{L^2(\mathscr {M}_{E})}+ h_{K}^{1/2}\left\| \nabla _h(u-J^{nc}_h u) \right\| _{L^2(\mathscr {M}_{E})}\right) . \end{aligned}$$

Squaring and summing up all the inner edges and using (3.18), we have

$$\begin{aligned} \sum _{E\in \mathscr {E}_h^I}\int _E \frac{|\textbf{b}\cdot \textbf{n}|}{2}[\hspace{-0.5pt}u-J^{nc}_h u \hspace{-0.5pt}]^2\;\textit{ds}\le C \left\| h_{\mathscr {T}}^{3/2}u \right\| ^{2}_{2}, \end{aligned}$$

here, the constant C depends on $\left\| \textbf{b} \right\| _{\infty }$. The result follows by combining all the above estimates. $\square$

Lemma 9

Suppose $u \in H^{2}(\varOmega )$ and $\beta _E = \beta {h_E}$ for some $\beta |_{{\mathscr {M}}_a} =\frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_a)}+C}$. Then,

$$\begin{aligned} A^{nc}_{h}(u-J^{nc}_h{u},{v}_h) \le C \left\| h_{\mathscr {T}}^{\frac{3}{2}} u \right\| _{2} \left| \! \left| \! \left| v_h\right| \! \right| \! \right| _{NLPSD} \quad \forall ~v_h \in V^{nc}_{h}. \end{aligned}$$

(4.47)

Proof

Using integration by parts in the first term of (4.40),

$$\begin{aligned} A^{nc}_{h}(u-J^{nc}_h{u},v_h)&= -({u-J^{nc}_h{u}}, \textbf{b}\cdot \nabla _h v_h) + ((\mu -\text{ div}_h \textbf{b}){(u-J^{nc}_h{u})}, v_h) \\ {}&\quad +S_{h}^{nc}(u-J^{nc}_h{u},v_h) +\int _{\partial \varOmega } (\textbf{b} \cdot \textbf{n})^{\oplus } (u-J^{nc}_h{u})v_h \;\textit{ds}\\ {}&\quad -\sum _{E\in \mathscr {E}_h^I}\int _E (\textbf{b}\cdot \textbf{n})\{\hspace{-1.0pt}u-J^{nc}_h{u} \hspace{-1.0pt}\}[\hspace{-0.5pt}v_h \hspace{-0.5pt}]\;\textit{ds}\\ {}&\quad + \sum _{E\in \mathscr {E}_h^I}\int _E \frac{|\textbf{b}\cdot \textbf{n}|}{2}[\hspace{-0.5pt}u-J^{nc}_h{u} \hspace{-0.5pt}][\hspace{-0.5pt}v_h \hspace{-0.5pt}]\;\textit{ds}. \end{aligned}$$

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,

$$\begin{aligned}&\sum _{E\in \mathscr {E}_h^I}\int _E (\textbf{b}\cdot \textbf{n}) \{\hspace{-1.0pt}u-J^{nc}_hu \hspace{-1.0pt}\}[\hspace{-0.5pt}v_h \hspace{-0.5pt}] \;\textit{ds}\le C\left( \sum _{E\in \mathscr {E}_h^I} \left\| \{\hspace{-1.0pt}u-J^{nc}_hu \hspace{-1.0pt}\} \right\| ^2_{L^2(E)}\right) ^{1/2}\\ {}&\quad \times \left( \sum _{E\in \mathscr {E}_h^I}\int _E \frac{|\textbf{b}\cdot \textbf{n}|}{2} [\hspace{-0.5pt}v_h \hspace{-0.5pt}]^2 \;\textit{ds}\right) ^{1/2}. \end{aligned}$$

Since $\beta _{E}=\beta h_{E}$ with $\beta = 1/({\left\| \textbf{b} \right\| _{W^{1,\infty }(\mathscr {M}_{E})}+C})$, and at the edge E, the jump term contributes to both the triangles sharing that edge; using the trace inequality (2.5),

$$\begin{aligned} \left\| \{\hspace{-1.0pt}u-J^{nc}_hu \hspace{-1.0pt}\} \right\| _{L^2(E)},&\left\| [\hspace{-0.5pt}u-J^{nc}_hu \hspace{-0.5pt}] \right\| _{L^2(E)} \\&\le C\left( h_{K}^{-1/2}\left\| u-J^{nc}_hu \right\| _{L^2(\mathscr {M}_{E})}+ h_{K}^{1/2}\left\| \nabla _h(u-J^{nc}_hu) \right\| _{L^2(\mathscr {M}_{E})}\right) . \end{aligned}$$

Squaring and summing up all the inner edges and using (3.18),

$$\begin{aligned} \sum _{E\in \mathscr {E}_h^I}\int _E (\textbf{b}\cdot \textbf{n}) \{\hspace{-1.0pt}u-J^{nc}_hu \hspace{-1.0pt}\}[\hspace{-0.5pt}v_h \hspace{-0.5pt}] \;\textit{ds}&\le C \left\| h_{\mathscr {T}}^{3/2}u \right\| _2\left| \! \left| \! \left| v_h\right| \! \right| \! \right| _{NLPSD}. \end{aligned}$$

Similarly,

$$\begin{aligned} \sum _{E\in \mathscr {E}_h^I}\int _E \frac{|\textbf{b}\cdot \textbf{n}|}{2} [\hspace{-0.5pt}u-J^{nc}_hu \hspace{-0.5pt}][\hspace{-0.5pt}v_h \hspace{-0.5pt}] \;\textit{ds}\le C \left\| h_{\mathscr {T}}^{3/2}u \right\| _2\left| \! \left| \! \left| v_h\right| \! \right| \! \right| _{NLPSD}. \end{aligned}$$

The constant C in the above estimates depends on $\left\| \textbf{b} \right\| _{\infty }$. Combining all these estimates leads to (4.47) and it concludes the proof. $\square$

Theorem 4

Let $u \in H^{2}(\varOmega )$ be the solution of continuous problem (2.3) and $u_h \in V^{nc}_{h}$ be the solution of discrete problem (4.39). Further, let $\beta _E=\beta h_E$ for some $\beta |_{{\mathscr {M}}_a} =\frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_a)}+C}$. Then,

$$\begin{aligned} \left| \! \left| \! \left| u-u_h\right| \! \right| \! \right| _{NLPSD}\le C \left\| h_{\mathscr {T}}^{3/2}u \right\| _2. \end{aligned}$$

(4.48)

Proof

By adding and subtracting the interpolation operator $J^{nc}_h u$, we decompose the error as follows:

$$\begin{aligned} \left| \! \left| \! \left| u-u_h\right| \! \right| \! \right| _{NLPSD}\le \left| \! \left| \! \left| u-J^{nc}_h u\right| \! \right| \! \right| _{NLPSD}+\left| \! \left| \! \left| J^{nc}_h u-u_h\right| \! \right| \! \right| _{NLPSD}. \end{aligned}$$

(4.49)

Let $w_h:= J^{nc}_h u-u_h \in V^{nc}_{h}$. Consider the second term of (4.49) and using the Theorem 3,

$$\begin{aligned}&c\left| \! \left| \! \left| u_h-J^{nc}_h u\right| \! \right| \! \right| _{NLPSD} \le \text{ sup}_{w_h \in V^{nc}_{h}}{\frac{A^{nc}_h( u_h-J^{nc}_h u,w_h)}{\left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{NLPSD}}}\nonumber \\&=\text{ sup}_{w_h \in V^{nc}_{h}}{\frac{A^{nc}_h( u_h-u,w_h)+ A^{nc}_h(u-J^{nc}_h u,w_h)}{\left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{NLPSD}}}. \end{aligned}$$

(4.50)

The weak formulation (2.4) and (4.39) imply that

$$\begin{aligned} A^{nc}_h( u_h-u,w_h)&= -S_{h}^{nc}(u,w_h). \end{aligned}$$

Moreover, the Cauchy–Schwarz inequality implies

$$\begin{aligned} S_{h}^{nc}(u,w_h)&=\sum _{E\in \mathscr {E}_h} \beta _E \left( \kappa _E( u), \kappa _E( w_h)\right) _{L^{2}(\mathscr {M}_E)} \\ {}&\le \left( \sum _{E\in \mathscr {E}_h} \beta _E\left\| \kappa _E( u) \right\| _{L^2(\mathscr {M}_E)}^2 \right) ^{1/2} \left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{NLPSD}. \end{aligned}$$

Note that $\beta _{E} = \beta h_{E}$ with $\beta |_{{\mathscr {M}}_E} =\frac{1}{\left\| \textbf{b} \right\| _{W^{{1,\infty }}({\mathscr {M}}_a)}+C}$. Using the Poincaré inequality (2.9) for every vertex $a\in \mathscr {V}_h$,

$$\begin{aligned} S_{h}^{nc}(u,w_h)&\le C \Big (\sum _{E\in \mathscr {E}_h} \beta _E h_E^2\left\| \nabla (\textbf{b}\cdot \nabla u) \right\| ^2_{L^2(\mathscr {M}_E)} \Big )^{1/2} \left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{NLPSD} \\&\le C\left\| h_{\mathscr {T}}^{3/2}u \right\| _2\left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{NLPSD}. \end{aligned}$$

It follows that,

$$\begin{aligned} A^{nc}_h( u_h-u,w_h)\le C\left\| h_{\mathscr {T}}^{3/2}u \right\| _2\left| \! \left| \! \left| w_h\right| \! \right| \! \right| _{NLPSD}. \end{aligned}$$

(4.51)

Using the estimate (4.51) and Lemma 9 in (4.50), we obtain

$$\begin{aligned} \left| \! \left| \! \left| u_h-J^{c}_h u\right| \! \right| \! \right| _{NLP}\le \left\| h_{\mathscr {T}}^{3/2}u \right\| _2. \end{aligned}$$

(4.52)

Finally, Lemma 8 and (4.52) lead (4.48) to an a priori estimate. $\square$

5 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 $\varOmega =(0, 1)^2$, coefficients $\textbf{b}=(3,2)$, $\mu =2$ and homogeneous Dirichlet boundary condition. The source term f is chosen such that the solution

$$\begin{aligned} u(x,y)=100x^2(1-x)^2y(1-y)(1-2y), \end{aligned}$$

satisfies the model problem. Further, the stabilization parameters for conforming and nonconforming FEMs are chosen as $\beta _a = 0.1 h_a$ and $\beta _E = 0.1 h_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 $\left| \! \left| \! \left| \cdot \right| \! \right| \! \right| _{LPSD}$ 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.

Table 1 Errors and convergence orders to the conforming FE solution of Example

Table 2 Errors and convergence orders to the nonconforming FE solution of Example

Fig. 2

Left: triangulation used for computations in Examples 1-4, right: nonconforming stabilized finite element solution of the Example 1

Fig. 3

Convergence history of GLPS method; left conforming approximations; right: nonconforming approximations of the Example (1)

Example 2

(Advection problem) Consider the model problem (2.1) with $\varOmega =(0, 1)^2$, coefficients $\textbf{b}=(0,1)$, $\mu =1$ and inflow boundary condition $g(x)=0$. The source term f is chosen such that the solution

$$\begin{aligned} u(x,y)=\frac{1}{2}\left( \tanh \left( \frac{y-.5}{0.04}\right) +1\right) \end{aligned}$$

satisfies the model problem.

5.1 Sensitivity with respect to problem and stabilization parameter

Even though the convergence rates in Theorems 2 and 4 hold in principle for any finite $\beta >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 $\mu$ has almost no influence on the error. We will set $\mu =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 $\left\| \textbf{b} \right\| _{\infty } \ge 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 $\beta$ behaves like a well with an approximation minimum at $\beta \in [0.1,1]$ for the conforming and a broader minimum at $\beta \in [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 $\beta$ approaches zero.

Fig. 4

Comparison of the GLPS 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

Further, the stabilization parameters for conforming and nonconforming finite element approximations are chosen as $\beta _a = 0.1 h_a$ and $\beta _E = 0.1 h_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.

Fig. 5

a Nonconforming Galerkin and b nonconforming stabilized finite element with $\beta _E = 0.1 h_E$ solutions of the Example (2)

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

Fig. 6

Convergence behaviour of the stabilized finite element solution of the Example 2

Example 3

(Circular internal layer) Consider the model problem (2.1) with $\varOmega =(0, 1)^2$, coefficients $\textbf{b}=(2,3)$ and $\mu =2$. The source term f and the inflow boundary condition are chosen such that the solution

$$\begin{aligned} u(x,y)= 16x(1-x)y(1-y)\left( \frac{1}{2}+ \frac{\tan ^{-1}\left( 200 \left( (0.25)^2-(x-.5)^2-(y-.5)^2)\right) \right) }{\pi } \right) \end{aligned}$$

satisfies the model equation.

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 $\beta _a = 0.06 h_a$ and $\beta _E = 0.05 h_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 $\left| \! \left| \! \left| \cdot \right| \! \right| \! \right|$ 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.

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

Fig. 7

The errors of GLPS finite element approximations and conforming stabilized finite element solution of the example (3)

Example 4

(Non-smooth solution)

Consider the model problem (2.1) with $\varOmega =(-1, 1)^2$, coefficients $\textbf{b}=(1,0)$, $\mu =0$, $f=0$, the inflow boundary condition

$$\begin{aligned} g(x) = \left\{ \begin{array}{ll} 1 &{} \quad y > 0 \\ 0 &{} \quad y < 0, \end{array} \right. \end{aligned}$$

and the exact solution

$$\begin{aligned} u(x,y) = \left\{ \begin{array}{ll} 1 &{} \quad y > 0 \\ 0 &{} \quad y < 0. \end{array} \right. \end{aligned}$$

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 $\beta _a = 0.7 h_a$ and $\beta _E = 0.7 h_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.

Fig. 8

Left side conforming stabilized solution and right side nonconforming stabilized solution of the example 4 with $h=2^{-6}$

6 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. Furthermore, 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.