New residual-based a posteriori error estimators for lowest-order Raviart-Thomas element approximation to convection-diffusion-reaction equations

A new technique of residual-type a posteriori error analysis is developed for the lowest-order Raviart-Thomas mixed finite element discretizations of convection-diffusion-reaction equations in two- or three-dimension. Both centered mixed scheme and upwind-weighted mixed scheme are considered. The a posteriori error estimators, derived for the stress variable error plus scalar displacement error in $L^{2}$-norm, can be directly computed with the solutions of the mixed schemes without any additional cost, and are robust with respect to the coefficients in the equations. Local efficiency dependent on local variations in coefficients is obtained without any saturation assumption, and holds from the cases where convection or reaction is not present to convection- or reaction-dominated problems. The main tools of analysis are the postprocessed approximation of scalar displacement, abstract error estimates, and the property of modified Oswald interpolation. Numerical experiments are reported to support our theoretical results and to show the competitive behavior of the proposed posteriori error estimates.


Introduction.
Let Ω ⊂ R d be a bounded polygonal or polyhedral domain in R d , d = 2 or 3. We consider the following homogeneous Dirichlet boundary value problem for the convection-diffusion-reaction equations: −∇ · (S∇p) + ∇ · (pw) + rp = f in Ω, p = 0 on ∂Ω, (1.1) where S ∈ L ∞ (Ω; R d×d ) denotes an inhomogeneous and anisotropic diffusion-dispersion tensor, w is a (dominating) velocity field, r a reaction function, f a source term. The choice of boundary conditions is made for ease of presentation, since similar results are valid for other boundary conditions. This type of equations arise in many chemical and biological settings. For instance, in hydrology these equations govern the transport and degradation of adsorbing contaminants and microbe-nutrient systems in groundwater. Reliable and efficient a posteriori error estimators are an indispensable tool for adaptive algorithms. For second-order elliptic problems without convection term, the theory of a posteriori error estimation has reached a degree of maturity for finite elements of conforming, nonconforming and mixed types; see [1-9, 11-14, 18, 20, 22-23, 27, 31-33] and the references therein. For convection-diffusion(-reaction) problems, on the contrary, the theory is still under development.
The mathematical analysis of robustness of a-posteriori estimators for the convectiondiffusion-reaction equations was first addressed by Verfürth [35] in the singular perturbation case, namely S = εI with I the identical matrix and 0 < ε ≪ 1. The proposed estimators for the standard Galerkin approximation and the SUPG disctetization give global upper and local lower bounds on the error measured in the energy norm, and are robust when the Péclet number becomes small. In [36] Verfürth improved the results of [35] in the sense that the derived estimates are fully robust with respect to convection dominance and uniform with respect to the size of the zero-order reaction term. Sangalli [30] developed an a posteriori estimator for the residual-free bubbles methods applied to convection-diffusion problems. Later he presented a residual-based a posteriori estimator for the one-dimensional convectiondiffusion-reaction model problem [31]. In [22] Kunert carried out a posteriori error estimation for the SUPG approach to a singularly perturbed convection-diffusion problem on anisotropic meshes. One may also refer to [25,26] for a posteriori error estimation in the framework of finite volume approximations.
For the convection-diffusion-reaction model (1.1), following an idea of postprocessing in [24] Vohralík [37] established residual a posteriori error estimates for lowest-order Raviart-Thomas mixed finite element discretizations on simplicial meshes. Global upper bounds and local lower bounds for the postprocessed approximation error, p−p h , in the energy norm were derived withp h the postprocessed approximation to the finite element solution p h , and the local efficiency of the estimators was shown to depend only on local variations in the coefficients and on the local Péclet number. Moreover, the developed general framework allows for asymptotic exactness and full robustness with respect to inhomogeneities and anisotropies.
In this paper, we develop a new technique for residual-based a posteriori estimation of the lowest-order Raviart-Thomas mixed finite element schemes (centered mixed scheme and upwind-mixed scheme) over both the stress error, u − u h , and the displacement error, p − p h , of the mixed finite element solutions (u h , p h ) for the problem (1.1) with u := −S∇p. The derived reliability results are robust with respect to the coefficients. Local efficiency dependent only on local variations in the coefficients is obtained without any saturation assumption, holds for the convection or reaction dominated equations. Compared with the standard analysis to the diffusion equations, our analysis avoids, by using the postprocessed approximatioñ p h as a transition, Helmholtz decomposition of stress variables and dual arguments of displacement error in L 2 -norm, and then does not need any weak regularity assumption on the diffusion-dispersion tensor. We note that although being employed in our analysis, the postprocessed displacement approximation and its modified Oswald interpolation are not involved in our estimators.
The rest of this paper is organized as follows. In Section 2 we give notations, assumptions of data, and the weak problem. We introduce in Section 3 the mixed finite element schemes (include the centered and upwind-weighted mixed scheme) and the post-processed techniques. Section 4 includes the main results. Section 5 collects some preliminary results and remarks. Section 6 and 7 analyze respectively the a posteriori error estimates and the local efficiency. Finally, we present several numerical examples in Section 8 to test our estimators.

Notations, assumptions and weak problem.
For a domain A ⊂ R d , we denote by L 2 (A) and L 2 (A) =: (L 2 (A)) d the spaces of square-integrable functions, by (·, ·) A the L 2 (A) or L 2 (A) inner product, by || · || A the associated norm, and by |A| the Lebesgue measure of A. Let H k (A) be the usual Sobolev space consisting of functions defined on A with all derivatives of order up to k square-integrable; product on ∂A for the duality paring between H −1/2 (∂A) and H 1/2 (∂A).
Let T h be a shape regular triangulation in the sense of [15] which satisfies the angle condition, namely there exists a constant c 0 such that for all K ∈ T h with h K := diam(K), Let C Q , c Q be positive constants dependent only on a quantity Q, and c i (i = 1, 2, · · · ) positive constants determined only by the shape regularity parameter, c 0 , of T h . We denote by ε h the set of element sides in T h , by ε int h and ε ext h the sets of all interior and exterior sides of T h , respectively. For K ∈ T h , denote by ε K the set of sides of K, especially by ε ext K the set of the boundary sides of K. Furthermore, we denote by ω σ andω σ the union of all elements in T h sharing a side σ and the union of all elements sharing at least one point of σ, respectively. For an element K ∈ T h the setω K is defined analogously. We also use the "broken Sobolev space" where × denotes the usual vector product of two vectors in R 3 . Given a unit normal vector n = (n 1 , n 2 ) along the side σ, we define the tangential component of a vector v ∈ R d by We note that throughout the paper, the local versions of differential operators ∇, curl are understood in the distribution sense, namely, curl h : We need in our analysis the following inequalities, Poincaré inequality and Friedrichs inequalities [10,27]: for K ∈ T h and ϕ ∈ H 1 (K), Here ϕ K := (1, ϕ) K /|K| and ϕ σ :=< 1, ϕ > σ /|σ| denote the integrable means of ϕ over K and over σ ∈ ε K , respectively. The constant C P,d can be evaluated as d/π for a simplex by using its convexity. Following [37], we suppose that there exists an original triangulation T 0 of Ω such that data of the problem (1.1) are given in the following way. Assumptions of data : (D1) S K := S| K is a constant, symmetric, and uniformly positive definite tensor such that As pointed out in [37], all the assumptions are made for the sake of simplicity and are usually satisfied in practice. If data do not satisfy these assumptions, we may employ the interpolation or projection of data with additional occurrence of data oscillation. Here the bilinear form and T h is a refinement of T 0 . We define as following an energy (semi) norm corresponding to the bilinear form B: We note that the weak problem (2.3) admits a unique solution under the Assumptions (D1)-(D6) [37].

Mixed finite element schemes and postprocessing.
Since it is of interest in many applications, the stress variable u := −S∇p are usually approximated by using the mixed finite elements for the problem (1.1). We introduce in this section the centered and upwindweighted mixed finite element schemes, and show the postprocessed techniques presented by Vohralík in [37].
We define the lowest order Raviart-Thomas finite element and piecewise constant space respectively as following: such that q h (x) = a + bx, for all x ∈ K. , Here n is the unit outer normal vector along σ ∈ ε h , and P 0 (K) denotes the set of constant functions on each K ∈ T h . We note that ∇ · (RT 0 (T h )) ⊂ P 0 (T h ). The centered mixed finite element scheme [17,37] reads as: 2) The upwind-weighted mixed finite element scheme [16,37] reads as: where w K,σ :=< 1, w · n > σ for σ ∈ ε K , with n the unit normal vector of σ, outward to K, ϕ K = (1, ϕ h ) K /|K| = ϕ h | K for all K ∈ T h , andp σ is the weighted upwind value given bŷ when σ is an interior side sharing by elements K and L, and bŷ when σ is a boundary side included in ε K . Here p K and p L denotes respectively the restrictions of p h over K and L, ν σ ∈ [0, 1/2] denotes the coefficient of the amount of upstream weighting which may be chosen as [37] where c S,σ is the harmonic average of c S,K and c S,L if σ ∈ ∂K ∩ ∂L and c S,K otherwise. We now introduce the postprocessed technique in [37], where a postprocessed approximationp h to the displacement p is constructed which links p h and u h on each simplex in the following way: We refer to [37] for the existence ofp h . We note that the new quantityp h ∈ W 0 (T h ) but / ∈ H 1 0 (Ω) (see LEMMA 6.1 in [37]), where

Main results.
With the stress variable u = −S∇p, we define the global and local errors, E and E K , of the stress and displacement variables as Denote respectively by η D,K and η R,K the elementwise displacement and residual estimator with Here the residual weight factors Note that in (4.4), if c w,r,K = 0, α K should be understood as h K / √ c S,K . Let ν σ be given in (3.7) for each side σ ∈ ε h . We denotê when σ is an interior side sharing by elements K and L, and when σ is a boundary side included in ε K . We thus define an elementwise upwind estimator η U,K by In order to reflect the change of the maximum eigenvalue of the coefficients matrix S over the patchω σ of a side σ ∈ ε h , we introduce a quantity Similarly, the change of one variation c w,r,K of the coefficients over the patchω K of an element K ∈ T h is described by the quantity Thus we define η N C,K as the elementwise nonconforming estimator by Since the convection occurs in the equations, we need to define two numbers Λ ∇·w,K and Λ w,σ similar to Péclet numbers describing the convection-dominated. To this end, for each K ∈ T h we denote C ∇·w,K := |∇ · w| K |, Λ ∇·w,K := max and for each σ ∈ ε h we set Λ w,σ := min{λ w,σ , p w,σ } with We then define η C,K as an elementwise convection estimator by We now state a posteriori error estimates for the global error of stress and displacement. be the solution of the centered mixed scheme (3.1)- (3.2). Let E be the error of the stress and displacement in the weighted norm defined in (4.1), η D,K , η R,K , η N C,K , and η C,K are the corresponding elementwise displacement estimator, residual estimator, convection estimator, and nonconforming estimator, defined in (4.2)-(4.3) and (4.8)-(4.9), respectively. Then it holds We note that the constants c 1 in (4.10) and c 2 in (4.11) only depend on the spatial dimension and the shape regularity parameter of the triangulation T h , and are independent of the coefficients S, w, r. In this sense, the proposed estimators are robust with respect to all the coefficients.

)). In his estimators, the term min
which shows that occurrence of ||hS −1 u h || Ω is reasonable in the a posteriori error estimators η D,K defined in (4.2). In addition, we note that the postprocessing (3.8) can remove the term ||hcurl(S −1 u h )|| Ω , which is also contained in Carstensen's estimators.
The global error estimates above show that the a posteriori indicator over each element consists of a series of estimators. Thus, the local efficiency of each component ensures the local efficiency of the a posteriori indicator over an element. Here, we point out the local efficiency is in the sense that its converse estimate holds up to a different multiplicative constant.   (Local efficiency for the nonconforming and convection estimators) Let η N C,K and η C,K be the elementwise nonconforming and convection estimators defined in (4.8) and (4.9), respectively. Then it holds and Λ w,r,K , Λ ∇·w,K are the same as in (4.8)- (4.9). We finally need the following quantities for the local efficiency of the upwind estimator over an element, where ν σ is given in (3.7) for each side σ ∈ ε h .

Preliminary results and remarks.
In this section, firstly we show the abstract error estimates developed by Vohralík in [37], and then make some remarks on Vohralík's a posteriori error estimators. To this end, for any ϕ ∈ H 1 0 (Ω) we define where ϕ K is the mean of ϕ over K, s ∈ H 1 0 (Ω) is arbitrarily given,p h is the postprocessed approximation solution given by (3.8)-(3.9), andp σ is the weighted upwind value defined in ifp h is the postprocessed solution, given by (3.8)- (3.9), of the centered mixed finite element scheme (3.1)- (3.2), and holds [37], the modified Oswald interpolation,

the abstract error estimates (5.4)-(5.5) so as to obtain computable estimates of the terms.
We now state our abstract error estimates for the global error of stress and displacement in the weighted norm.
LEMMA 5.3. For any K ∈ T h and ϕ ∈ H 1 (K), it holds where ϕ K denotes the mean of ϕ over K, and α K is defined as in (4.4).

A posteriori error analysis.
We devote this section to computable estimates of T R (ϕ), T U (ϕ) and T C (ϕ, s) defined in (5.1), (5.3) and (5.2), respectively, with the help of u h and p h . Moreover, we derive an estimate of |||p h − s||| by substituting s with the modified Oswald interpolation I MO (p h ) (see [37]), and by using the postprocessing technique as a transition. Finally, we give the proof of THEOREMs 4.1-4.2. LEMMA 6.1. (Residual estimator) Let T R (ϕ) be defined as in (5.1) with |||ϕ||| Ω = 1, and η R,K be defined as in (4.3). Then it holds   .3) with |||ϕ||| Ω = 1, and η U,K be defined as in (4.7). Then it holds Proof. We denote byp σ the mean ofp h over σ ∈ ε h , i.e.,p σ :=< 1,p h > σ /|σ|. The definitions of T U (ϕ) and w K,σ , together with Assumption (D2) of the velocity field w, imply For an element K ∈ T h , it holds σ ∈ ε K ∩ ε L or σ ∈ ε ext K . For the former case, recalling p K = p h | K , p L = p h | L , from the postprocessing (3.9) we obtain For the latter case, we similarly havê For convenience, in what follows we denotê when σ ∈ ε ext K . In light of the definitions ofp σ andp σ in (3.5)-(3.6) and (4.5)-(4.6), and from (6.4)-(6.6) we have Since ϕ ∈ H 1 0 (Ω), andp σ ,p ωσ are constants over a side σ ∈ ε h , it holds From Friedrichs inequality (2.2) and the postprocessing (3.8) we have The trace inequality (see LEMMA 3.1 in [35]) and local shape regularity of elements indicate
For the first term, |||p h − s||| Ω , in the right side of the abstract error estimate (5.6) or (5.7), we follow [37] to take s := I MO (p h ) in the sequel, where I MO (p h ) is the modified Oswald interpolation ofp h . Recall an estimate on the modified Oswald interpolation [20], or 3 ) denotes the set of polynomials of degree at most d on each simplex, σ ∩ K = ∅ when σ contains a vertex of K.
By definition we have . LEMMAs 6.3-6.4 show respectively computable estimates of the two right-side terms of the above identity with the help of u h and p h . LEMMA 6.3. Let γ tσ (·) be defined as in Section 2.1, and s := I MO (p h ). Then it holds where Λ σ is given in Section 4, and t σ denotes the unit tangent vector along σ.
Proof. From the estimate (6.14) we have where σ ∩ K = ∅ when σ contains a vertex of K.
Since the mean ofp h over interior side is continuous and its mean on exterior side vanishes, i.e., A combination of (6.16)-(6.18) yields Summing (6.19) over each element K, noticing that the number of summation over a side σ ∈ ε h is bounded by a positive constant c 17 , and combining the definition of Λ σ , we obtain REMARK 6.1. The node with respect to which the quasi-monotone condition is violated is called singular node (cf. [28]). We can derive an alternative form of (6.20) as following: with C S,ωK := max K ′ ∈ωK C S,K ′ . LEMMA 6.4. Let Λ w,r,K be the same as in (4.8) and s := I MO (p h ). Then it holds A combination of (6.22), (6.23) and the postprocessing (3.8) yields that Summing (6.24) over each element K, noticing that the mesh is local quasi-uniform, and combining the definition of Λ w,r,K , we finally get REMARK 6.2. (Alternative form) For a side σ ∈ ε h , we denote Λ w,r,σ := max K,K∩σ =∅ {c w,r,K }.A combination of (6.17), (6.18) and (6.22) yields which leads to an alternative form of the estimate (6.21),

This inequality shows that the term
The following corollary is a combined result of LEMMAs 6.3-6.4. COROLLARY 6.5. Let η N C,K be defined as in (4.8) and s := I MO (p h ). Then it holds LEMMA 6.6. (Convection estimator.) Let T C (ϕ, s) be defined as in (5.2) with |||ϕ||| Ω = 1 and s := I MO (p h ), and η C,K be defined as in (4.9). Then it holds Proof. By triangle inequality and Hölder inequality we obtain (6.28) Apply (6.19) and (6.24) to the inequality (6.28), and combine the definitions of λ w,σ and Λ ∇·w,K , we then arrive at (6.29) Since the modified Oswald interpolation s = I MO (p h ) preserves the mean ofp h on the side, and w · n is constant over a side, it holds where ϕ K is the mean of ϕ over K. Write v :=p h − s, then we have (6.30) A combination of (6.30), (2.1), (6.19), (6.24) and Hölder inequality yields (6.31) This estimate, together with the definitions of p w,σ and Λ ∇·w,K , indicates T C (ϕ, s) from (6.31) The desired result (6.27) follows from (6.29) and (6.32) with c 24 = max{c 25 , c 26 }.
(6.33) Thus the a posteriori error estimate (4.10) is reduced to In addition, Remark 6.1 implies an alternative estimate Note that being an oscillation term, the first term in the right side of (6.34) or (6.35) may not be computed in practice. REMARK 6.5. (A posteriori error estimate of divergence of the stress variable.) The continuous weak formulation of (1.1) reads as: Find (u, p) ∈ H(div, Ω) × L 2 (Ω) such that Notice that (6.36) can be equivalently written as: For each K ∈ T h (∇ · u, ϕ) K − (S −1 u · w, ϕ) K + ((r + ∇ · w)p, ϕ) K = (f, ϕ) K for all ϕ ∈ L 2 (K). (6.37) Meanwhile, the centered mixed finite element scheme (3.2) can be equivalently written as: For every K ∈ T h (∇ · u h , ϕ) K − (S −1 u h · w, ϕ) K + ((r + ∇ · w)p h , ϕ) K = (f, ϕ) K for all ϕ ∈ P 0 (K). (6.38) Let R K denote the mean of the elementwise residual over K ∈ T h , and set 0 ≤ ι ≤ 1. We define the data oscillation osc h as For any ϕ ∈ L 2 (K), let ϕ K denote the mean of ϕ over K ∈ T h , then a combination of (6.37) and (6.38) yields This means that (6.39) From (6.39) and (4.10) we obtain the following a posteriori error estimate of the divergence of the stress variable for the centered mixed finite element scheme: where the constant β c := max K∈T h max{ Cw,r,K √ cw,r,K h ι K }.

Analysis of local efficiency.
Using standard arguments we easily derive lemmas 7.1-7.2. LEMMA 7.1.
Denote v := f − ∇ · u h + (S −1 u h ) · w − (r + ∇ · w)p h , and let E K be the local error for the stress and displacement defined in (4.1). Under Assumption (D5) for f , it holds Proof. For all K ∈ T h , let ψ K denote the bubble function on K with zero boundary values on K and 0 ≤ ψ K ≤ 1. The relation Integration by parts implies A combination of (7.4), (7.5) and inverse inequality imply The desired result (7.3) then follows with c 31 := √ 2c 32 . LEMMA 7.4. It holds This relation, together with trace theorem and the postprocessing (3.9), indicates which, together with the local shape regularity of elements and the postprocessing (3.8), implies If σ ∈ ε K ∩ ε ext h , from σp h | K ds = 0 and ν σ ≤ 1/2 the second assertion of the lemma follows.

Proof of THEOREM 4.4: Notice that
From the definitions of β * ,K and c ωσ in this theorem, we respectively apply LEMMAs 7.2-7.3 to the above two inequalities so as to obtain The assertion (4.13) follows from (7.12) and (7.13) by taking c 4 := 2(d + 1) max{c 36 , c 37 }.
Proof of THEOREM 4.5: The local shape regularity of elements implies Then the desired estimate (4.14) follows from (7.14), LEMMA 7.4 and the definitions of the constants λ σ , ρ σ , and E D,ωσ .

Numerical experiments.
In this section, we test our proposed posteriori error estimators on three model problems. where ρ, θ are the polar coordinates.
It is well known that this model possesses singularity at the origin. The original mesh consists of 6 right-angled triangles. We employ the centered mixed scheme described in section 3.1 to compute the approximaton solution, mark elements in terms of Dörfler marking with the marking parameter θ = 0.5, and then use the "longest edge" refinement to recover an admissible mesh. Specially, the uniform refinement means that all elements should be marked. We note that in the given case, the residual estimators η R,K vanish over all K ∈ T h . We see in the first figure of Fig 8.1 with 1635 elements that the refinement concentrates around the origin, which means the predicted error estimator captures well the singularity of the solution. The second graph of Fig 8.1 reports the estimated and actual errors of the numerical solutions on uniformly and adaptively refined meshes. It can be seen that one can substantially reduce the number of unknowns necessary to obtain the prescribed accuracy by using the a posteriori error estimates and adaptively refined meshes, and that the error of the flux in L 2 norm uniformly reduces with a fixed factor on two successive meshes, and that the adaptive mixed finite element method is a contraction with respect to the energy error. [18,29,37]. We consider the problem (1.1) in a square domain Ω = (−1, 1) × (−1, 1) with w = r = 0 and f = 0, where Ω is divided into four subdomains Ω i (i = 1, 2, 3, 4) corresponding to the axis quadrants (in the counterclockwise direction), and the diffusion-dispersion tensor S is piecewise constant matrix with S = s i I in Ω i . We suppose the exact solution of this model has the form p(r, θ) = r α (a i sin(αθ) + b i cos(αθ)) in each Ω i with Dirichlet boundary conditions. Here r, θ are the polar coordinates in Ω, a i and b i are constants depending on Ω i , and α is a parameter. We note that the stress solution u = −S∇p is not continuous across the interfaces, and only its normal component is continuous. It finally exhibits a strong singularity at the origin. We consider two sets of coefficients in the following table: The origin mesh consists of 8 right-angled triangles. We use the centered scheme compute the approximation solution, and mark elements in terms of Dörfler marking with the marking parameter θ = 0.7 in the first case and θ = 0.94 in the second case. We note that the elementwise estimators ξ K are used as the a posteriori error indicators, since the residual estimators η R,K vanish over K ∈ T h .

Model problem with inhomogeneous diffusion tensor
In Table 8.1 we show for Case 1 some results of the actual error E k , the a posteriori indicator η k , the experimental convergence rate, EOC E , of E k , and the experimental convergence rate, EOC η , of η k , where , and DOF k denotes the number of elements with respect to the k−th iteration. We can see that the convergence rates EOC E and EOC η are close to 0.5 as the iteration number k = 15, which means the optimal decay of the actual error and a posteriori error indicator η k is almost attained after 15 iterations with optimal meshes.  8.2 shows an adaptively refined mesh with 4763 elements and the estimated and actual errors against the number of elements in adaptively refined meshes for Case 1 . Fig 8.3 shows an adaptively refined mesh with 1093 elements and the actual error against the number of elements in adaptively refined meshes for Case 2.
From the first figures of Fig 8.2-8.3, we can see that the refinement again concentrates around the origin, which means the adaptive mixed finite element method detects the region of rapid variation. In the second graphs of Fig 8.1-8.3 each includes an optimal convergence line, which shows in both cases, the energy error performs a trend of descend with an optimal order convergent rate. Simultaneously, from the second graphs of Fig 8.1-8.3, we also see that the proposed estimators are efficient with respect to the strongly discontinuously coefficients.
We note that the energy error is approximated with a 7-point quadrature formula in each triangle.  with a a positive constant. This solution is, in fact, one-dimensional and possesses an internal layer of width a which we shall set, respectively, equal to 0.1, 0.05, 0.02, 0.001. We still start computations from an origin mesh which consists of 8 right-angled triangles, and refine it either uniformly (up to five refinements) or adaptively.
In Fig 8.4 with ε = 0.01, a = 0.05 and Fig 8.5 with ε = 0.001, a = 0.05, we can see that the refinement concentrates at an internal layer of width a = 0.05, and is away from the center of the shock. Both the convection-dominated regime on coarse grids and diffusiondominated regime obtain the progressive refinement. The effect is still rather good even if the approximation to displacement is piecewise constant.  We again see that the refinement focuses around layer of width a = 0.001, this indicates that the estimators actually capture interior layers and resolve them in convection-domianed regions. In addition, the postprocessing approximation to the scalar displacement obtains a satisfactory result.
In Fig 8.7 with ε = 0.1, a = 0.02 (left), the estimated and actual errors are plotted against the number of elements in uniformly and adaptively refined meshes. Again, we see that one can substantially reduce the unknowns necessary to attain the prescribed precision by using the proposed estimators and adaptively refined grids. The second graph of Fig 8.7 shows the actual error against the number of elements in adaptively refined meshes for different ε in case a = 0.1, and also concludes a line with optimal convergence −1/2. In addition, we also see that the almost same error decay occurs in cases: ε = 0.01 and ε = 0.001.