A Simple Approach to Reliable and Robust A Posteriori Error Estimation for Singularly Perturbed Problems

TL;DR

This paper introduces a simple flux reconstruction method for finite element solutions that provides reliable, fully computable a posteriori error bounds for singularly perturbed reaction-diffusion problems, ensuring robustness and efficiency.

Contribution

The paper proposes a straightforward flux reconstruction technique that yields guaranteed error bounds and efficient local error indicators, even in the presence of singular perturbations.

Findings

Provides fully computable upper bounds on error in energy norm.

Ensures local efficiency and robustness of error indicators.

Works even when Galerkin orthogonality is violated.

Abstract

A simple flux reconstruction for finite element solutions of reaction-diffusion problems is shown to yield fully computable upper bounds on the energy norm of error in an approximation of singularly perturbed reaction-diffusion problem. The flux reconstruction is based on simple, independent post-processing operations over patches of elements in conjunction with standard Raviart--Thomas vector fields and gives upper bounds even in cases where Galerkin orthogonality might be violated. If Galerkin orthogonality holds, we prove that the corresponding local error indicators are locally efficient and robust with respect to any mesh size and any size of the reaction coefficient, including the singularly perturbed limit.

Figures

Equations

$$-\Delta u+{\kappa }^{2}u=f\text{in }\Omega ;u=0\text{on }{\Gamma }_{\mathrm{D}};\partial u/\partial \mathbf{\nu }=g_{\mathrm{N}}\text{on }{\Gamma }_{\mathrm{N}},$$

$$\underset{{\mathcal{T}}_h\in \mathcal{G}}{\sup }\underset{K\in {\mathcal{T}}_h}{\max }\frac{h_K}{{\rho }_K}\le C.$$

$$\tilde{K}=\mathrm{int}\bigcup \left\{K^{\prime }\in {\mathcal{T}}_h:K^{\prime }\cap K\ne \varnothing \right\}.$$

$$\text{if}{h}_K{\kappa }_K>1\text{then}{\kappa }_K\le C{\kappa }_{K^{\prime }}\text{for all }{K}^{\prime }\subset \tilde{K}.$$

$$C^{-1}{\kappa }_{K^{\prime }}\le {\kappa }_K\le C{\kappa }_{K^{\prime }}$$

$$C^{-1}\min \{h_{K^{\prime }},{\kappa }_{K^{\prime }}^{-1}\}\le \min \{h_K,{\kappa }_K^{-1}\}\le C\min \{h_{K^{\prime }},{\kappa }_{K^{\prime }}^{-1}\}.$$

$$\min \{h_K,{\kappa }_K^{-1}\}=h_K\text{if }{\kappa }_K=0.$$

$$\mathcal{B}(u,v)=\mathcal{F}(v)\forall v\in V,$$

$$\mathcal{B}(u,v)={\int }_{\Omega }(\nabla u\cdot \nabla v+{\kappa }^{2}uv)\mathrm{d}\mathbf{x};\mathcal{F}(v)={\int }_{\Omega }fv\mathrm{d}\mathbf{x}+{\int }_{{\Gamma }_{\mathrm{N}}}{g}_{\mathrm{N}}v\mathrm{d}\mathbf{s}.$$

$${\mathcal{B}}_K(u,v)={\int }_K(\nabla u\cdot \nabla v+{\kappa }_K^{2}uv)\mathrm{d}\mathbf{x};{\mathcal{F}}_K(v)={\int }_{K}fv\mathrm{d}\mathbf{x}+{\int }_{{\Gamma }_{\mathrm{N}}\cap \partial K}{g}_{\mathrm{N}}v\mathrm{d}\mathbf{s}.$$

$$\mathcal{B}(u_h,v_h)=\mathcal{F}(v_h)\forall v_h\in V_h.$$

$$-\mathrm{div}\mathbf{\tau }$$

$$\mathbf{\tau }\cdot \mathbf{\nu }$$

$${\eta }_K(\mathbf{\tau })={\left({\Vert \mathbf{\epsilon }\Vert }_K^2+{\kappa }_K^{-2}{\Vert r\Vert }_K^2\right)}^{1/2}+\sum _{\gamma \subset {\Gamma }_{\mathrm{N}}\cap \partial K}{C}_{\mathrm{T}}^{K,\gamma }{\Vert R_{\mathrm{N}}\Vert }_{\gamma }.$$

$${\eta }^2(\mathbf{\tau })=\sum _{K\in {\mathcal{T}}_h}[{\eta }_K(\mathbf{\tau })+{\mathrm{osc}}_K(f,g_{\mathrm{N}}){]}^2$$

$${\mathrm{osc}}_K(f,g_{\mathrm{N}})=\min \left\{\frac{h_K}{\pi },\frac{1}{{\kappa }_K}\right\}{\Vert f-{\Pi }_{K}f\Vert }_K+\sum _{\gamma \subset {\Gamma }_{\mathrm{N}}\cap \partial K}\min \{C_{\mathrm{T}}^{K,\gamma },{\overline{C}}_{\mathrm{T}}^{K,\gamma }\}{\Vert g_{\mathrm{N}}-{\Pi }_{\gamma }{g}_{\mathrm{N}}\Vert }_{\gamma }$$

$${\left(C_{\mathrm{T}}^{K,\gamma }\right)}^2$$

$${\left({\overline{C}}_{\mathrm{T}}^{K,\gamma }\right)}^2$$

$$|||u-u_h|||\le \eta (\mathbf{\tau }).$$

$$\begin{gathered} \mathcal{B}(u-u_h,v)=\sum _{K\in {\mathcal{T}}_h}[(\mathbf{\epsilon },\nabla v{)}_K+(r,v{)}_K+\sum _{\gamma \subset {\Gamma }_{\mathrm{N}}\cap \partial K}(R_{\mathrm{N}},v{)}_{\gamma } \\ +(f-{\Pi }_{K}f,v{)}_K+\sum _{\gamma \subset {\Gamma }_{\mathrm{N}}\cap \partial K}(g_{\mathrm{N}}-{\Pi }_{\gamma }{g}_{\mathrm{N}},v{)}_{\gamma }] \end{gathered}$$

$$(f-{\Pi }_{K}f,v{)}_K+\sum _{\gamma \subset {\Gamma }_{\mathrm{N}}\cap \partial K}(g_{\mathrm{N}}-{\Pi }_{\gamma }{g}_{\mathrm{N}},v{)}_{\gamma }\le {\mathrm{osc}}_K(f,g_{\mathrm{N}})|||v||{|}_K.$$

$$(\mathbf{\epsilon },\nabla v{)}_K+(r,v{)}_K\le {\left({\Vert \mathbf{\epsilon }\Vert }_K^2+{\kappa }_K^{-2}{\Vert r\Vert }_K^2\right)}^{1/2}|||v||{|}_K,$$

$$(R_{\mathrm{N}},v{)}_{\gamma }\le C_{\mathrm{T}}^{K,\gamma }{\Vert R_{\mathrm{N}}\Vert }_{\gamma }|||v||{|}_K,$$

$$\begin{gathered} \mathcal{B}(u-u_h,v)\le \sum _{K\in {\mathcal{T}}_h^{+}}\left[{\left({\Vert \mathbf{\epsilon }\Vert }_K^2+{\kappa }_K^{-2}{\Vert r\Vert }_K^2\right)}^{1/2}+\sum _{\gamma \subset {\Gamma }_{\mathrm{N}}\cap \partial K}{C}_{\mathrm{T}}^{K,\gamma }{\Vert R_{\mathrm{N}}\Vert }_{\gamma }\right]|||v||{|}_K \\ +\sum _{K\in {\mathcal{T}}_h^0}{\Vert \mathbf{\epsilon }\Vert }_K|||v||{|}_K+\sum _{K\in {\mathcal{T}}_h}{\mathrm{osc}}_K(f,g_{\mathrm{N}})|||v||{|}_K=\sum _{K\in {\mathcal{T}}_h}[{\eta }_K(\mathbf{\tau })+{\mathrm{osc}}_K(f,g_{\mathrm{N}})]|||v||{|}_K \end{gathered}$$

$$\mathbf{W}({\omega }_{\mathbf{n}})=\{\mathbf{\tau }\in \mathbf{H}(\mathrm{div},{\omega }_{\mathbf{n}}):\mathbf{\tau }{|}_K\in {\mathbf{RT}}_1(K),\mathbf{\tau }\cdot {\mathbf{\nu }}_{\mathbf{n}}=0\text{ on }\gamma \in {\mathcal{E}}_{\mathbf{n}}^{\mathrm{E}}\},$$

$$\begin{gathered} E_{\mathbf{n}}({\mathbf{\tau }}_{\mathbf{n}})={\Vert {\mathbf{\tau }}_{\mathbf{n}}-{\theta }_{\mathbf{n}}\nabla u_h\Vert }_{{\omega }_{\mathbf{n}}}^2 \\ +{\Vert {\kappa }^{-1}\left[\Pi ({\theta }_{\mathbf{n}}(\Pi f-{\kappa }^2{u}_h))-\nabla {\theta }_{\mathbf{n}}\cdot \nabla u_h+\mathrm{div}{\mathbf{\tau }}_{\mathbf{n}}\right]\Vert }_{{\omega }_{\mathbf{n}}^{+}}^2 \\ +{\Vert C_{\mathrm{T}}^{\mathrm{N}}\left[{\Pi }_{\mathrm{N}}({\theta }_{\mathbf{n}}{\Pi }_{\mathrm{N}}{g}_{\mathrm{N}})-{\mathbf{\tau }}_{\mathbf{n}}\cdot \mathbf{\nu }\right]\Vert }_{{\Gamma }_{\mathbf{n}}^{\mathrm{N}+}}^2, \end{gathered}$$

$$\Pi ({\theta }_{\mathbf{n}}(\Pi f-{\kappa }^2{u}_h))-\nabla {\theta }_{\mathbf{n}}\cdot \nabla u_h+\mathrm{div}{\mathbf{\tau }}_{\mathbf{n}}$$

$${\Pi }_{\mathrm{N}}({\theta }_{\mathbf{n}}{\Pi }_{\mathrm{N}}{g}_{\mathrm{N}})-{\mathbf{\tau }}_{\mathbf{n}}\cdot \mathbf{\nu }$$

$$\begin{gathered} ({\mathbf{\tau }}_{\mathbf{n}},{\mathbf{w}}_h{)}_{{\omega }_{\mathbf{n}}}+({\kappa }^{-2}\mathrm{div}{\mathbf{\tau }}_{\mathbf{n}},\mathrm{div}{\mathbf{w}}_h{)}_{{\omega }_{\mathbf{n}}^{+}}+((C_{\mathrm{T}}^{\mathrm{N}}{)}^2{\mathbf{\tau }}_{\mathbf{n}}\cdot \mathbf{\nu },{\mathbf{w}}_h\cdot \mathbf{\nu }{)}_{{\Gamma }_{\mathbf{n}}^{\mathrm{N}+}} \\ +(q_h,\mathrm{div}{\mathbf{w}}_h{)}_{{\omega }_{\mathbf{n}}^0}+(d_h,{\mathbf{w}}_h\cdot \mathbf{\nu }{)}_{{\Gamma }_{\mathbf{n}}^{\mathrm{N}0}}=({\theta }_{\mathbf{n}}\nabla u_h,{\mathbf{w}}_h{)}_{{\omega }_{\mathbf{n}}} \\ -{\left({\kappa }^{-2}\left[{\theta }_{\mathbf{n}}(\Pi f-{\kappa }^2{u}_h)-\nabla {\theta }_{\mathbf{n}}\cdot \nabla u_h\right],\mathrm{div}{\mathbf{w}}_h\right)}_{{\omega }_{\mathbf{n}}^{+}}+((C_{\mathrm{T}}^{\mathrm{N}}{)}^2{\theta }_{\mathbf{n}}{\Pi }_{\mathrm{N}}{g}_{\mathrm{N}},{\mathbf{w}}_h\cdot \mathbf{\nu }{)}_{{\Gamma }_{\mathbf{n}}^{\mathrm{N}+}} \end{gathered}$$

$$(\mathrm{div}{\mathbf{\tau }}_{\mathbf{n}},{\phi }_h{)}_{{\omega }_{\mathbf{n}}^0}$$

$$({\mathbf{\tau }}_{\mathbf{n}}\cdot \mathbf{\nu },{\psi }_h{)}_{{\Gamma }_{\mathbf{n}}^{\mathrm{N}0}}$$

$$\mathbf{\tau }=\sum _{\mathbf{n}\in {\mathcal{N}}_h}{\mathbf{\tau }}_{\mathbf{n}}.$$

$$0=\sum _{\mathbf{n}\in {\mathcal{N}}_K}({\theta }_{\mathbf{n}}({\Pi }_{K}f-{\kappa }_K^2{u}_h)-\nabla {\theta }_{\mathbf{n}}\cdot \nabla u_h+\mathrm{div}{\mathbf{\tau }}_{\mathbf{n}},{\phi }_h{)}_K=({\Pi }_{K}f-{\kappa }_K^2{u}_h+\mathrm{div}\mathbf{\tau },{\phi }_h{)}_K$$

Full text

Keywords: finite element analysis, robust a posteriori error estimate, singularly perturbed problems, flux reconstruction

MSC: 65N15, 65N30, 65J15

1 Introduction

This paper is concerned with developing fully computable bounds for the error in the finite element approximation of the following linear reaction-diffusion problem

[TABLE]

where the domain $\Omega\subset\mathbb{R}^{d}$, $d\geq 2$, is a polytope and ${\boldsymbol{\nu}}$ denotes the unit outward normal vector on the boundary $\partial\Omega$. Here, the portions ${\Gamma_{\mathrm{D}}}$ and ${\Gamma_{\mathrm{N}}}$ of the boundary $\partial\Omega$ are open, disjoint and satisfy $\overline{\Gamma}_{\mathrm{D}}\cup\overline{\Gamma}_{\mathrm{N}}=\partial\Omega$. For simplicity, we assume that the data $f\in L^{2}(\Omega)$ and $g_{\mathrm{N}}\in L^{2}({\Gamma_{\mathrm{N}}})$ and that the reaction coefficient $\kappa\geq 0$ is piecewise constant. The problem has a unique solution provided that either ${\Gamma_{\mathrm{D}}}$ has a positive measure or $\kappa$ is not identically zero.

Given a conforming approximation $u_{h}$ of the true solution $u$ of problem (1), we present a novel a posteriori error estimator for the energy norm of the error $|\!|\!|u-u_{h}|\!|\!|$. The estimator is rather easy to evaluate (via a fast element by element algorithm) and provides a guaranteed upper bound on the true error measured in the energy norm. In the case where $u_{h}$ is the Galerkin finite element approximation of $u$, we prove that the estimator is locally efficient and provides an upper bound which does not degenerate in the singularly perturbed limit (i.e. when $\kappa\to\infty$).

This error estimator is evaluated using a reconstructed flux. In contrast to our previous work [1, 3, 4], the reconstructed flux is obtained by solving small local problems on patches of elements by Raviart–Thomas finite elements. This approach is technically simpler and yields a more accurate error estimator.

The idea of the flux reconstruction by solving small problems on patches comes from [6]. It can be seen as approximate minimization of the error bound using an overlapped domain decomposition method with subdomains chosen as patches of elements. Local minimization problems on these patches have equilibration constraints and are solved by by mixed finite elements. Our result can be seen as a robust generalization of this idea to reaction-diffusion problems.

The general idea of flux reconstructions, however, dates back to the method of hypercircle [26, 30] and later to [5, 15, 17, 20, 35]. In the last two decades it was vastly developed, see e.g. [1, 8, 14, 21, 25, 27] and references there in. Interestingly, error estimates based on flux reconstructions can be utilized to estimate various components of the error such as the discretization, iteration, and algebraic errors [10, 12, 16, 23]. This enables us to adaptively equilibrate all components of the error and develop algorithms that do not perform excessive iterations of linear and nonlinear solvers in cases when the iteration and algebraic errors are already on the level of the discretization error.

The first robust, reliable, and locally efficient a posteriori error estimate for problem (1) was derived by Verfürth in [33, 32]. Locally efficient and robust guaranteed error bounds for the vertex-centred finite volume discretization of (1) were proposed in [9]. Robust reliability estimate for singularly perturbed problem on anisotropic meshes is proved in [18]. A similar result for a guaranteed and fully computable error bound is provided in [19]. An interesting alternative idea for guaranteed upper bounds for reaction-diffusion problems was recently published in [24]. Preprint [29] proofs robustness of a simple a posteriori error estimator, however their approach considerably differs from the one presented below due to equilibration of fluxes even if the reaction term dominates and due to the presence of weights in the estimator.

The rest of the paper is organized as follows. Section 2 briefly introduces the finite element approximation of problem (1) and the corresponding notation. Section 3 defines the a posteriori error estimator and proves that it is the guaranteed upper bound on the error. Section 4 introduces flux reconstruction based on local minimization problems and Sections 5 and 6 define two auxiliary flux reconstructions that are used in Section 7 to prove the robust local efficiency of the proposed error indicators. Section 8 proposes alternative flux reconstruction that does not require equilibration condition. Section 9 provides a couple of numerical examples and Section 10 draws the conclusions.

2 Model Problem and Its Discretization

2.1 Partitions

Let $\mathcal{G}=\{\mathcal{T}_{h}\}$ be a family of partitionings of the domain $\Omega$ into simplicial elements. The intersection of each distinct pair of elements in a given partition $\mathcal{T}_{h}\in\mathcal{G}$ is assumed to consist of a single common vertex or a single common facet of both elements. The diameter and inradius of an element $K$ are denoted by $h_{K}$ and $\rho_{K}$, respectively. The family $\mathcal{G}$ is assumed to be regular in the sense that there exists a constant $C>0$ such that

[TABLE]

This assumption permits meshes in which the elements are locally refined such as might arise from an adaptive refinement algorithm. The patch consisting of an element $K\in\mathcal{T}_{h}$ and those elements in $\mathcal{T}_{h}$ sharing at least one common point with $K$ is defined by

[TABLE]

The regularity condition (2) means that the number of elements in any patch is uniformly bounded over the family $\mathcal{G}$, as is the number of patches containing a particular element. Further, condition (2) implies the following local quasi-uniformity and shape regularity properties: there exist constants $c>0$ and $C>0$ such that for all elements $K^{\prime}\subset\widetilde{K}$, all $K\in\mathcal{T}_{h}$, and all $\mathcal{T}_{h}\in\mathcal{G}$ estimates $ch_{K}\leq h_{K^{\prime}}\leq Ch_{K}$ and $c\rho_{K}\leq\rho_{K^{\prime}}\leq C\rho_{K}$ hold.

Here, and throughout, we adopt the convention whereby the symbol $C$ is used to denote a generic constant throughout the paper, whose actual numerical value can differ in different occurrences, but it is always independent of $\kappa$ and any mesh-size.

The notation $(\cdot,\cdot)_{\omega}$ and $\left\|\cdot\right\|_{\omega}$ is used to denote the $L^{2}(\omega)$ scalar product and norm over a subset $\omega\subset\Omega$, and we omit the subscript in the case when $\omega=\Omega$. The $L^{2}(K)$-orthogonal projector onto the space of affine functions $\mathbb{P}_{1}(K)$ over element $K\in\mathcal{T}_{h}$ is denoted by $\Pi_{K}:L^{2}(K)\rightarrow\mathbb{P}_{1}(K)$, whist $\Pi$ is used to denote the concatenation of the elementwise projections $\Pi_{K}$, i.e. $(\Pi f)|_{K}=\Pi_{K}f$ for all $K\in\mathcal{T}_{h}$. Similarly, for a facet $\gamma\subset{\Gamma_{\mathrm{N}}}\cap\partial K$, $\Pi_{\gamma}:L^{2}(\gamma)\rightarrow\mathbb{P}_{1}(\gamma)$ denotes the $L^{2}(\gamma)$-orthogonal projector, and $\Pi_{\mathrm{N}}$ denotes the concatenation of the facetwise projections $\Pi_{\gamma}$: $(\Pi_{\mathrm{N}}g_{\mathrm{N}})|_{\gamma}=\Pi_{\gamma}g_{\mathrm{N}}$ for all facets $\gamma\subset{\Gamma_{\mathrm{N}}}$.

2.2 Assumptions on the Reaction Coefficient $\kappa$

For simplicity, we shall assume that the reaction coefficient is constant on every element over the entire set of partitions in $\mathcal{G}$ and we denote by $\kappa_{K}=\kappa|_{K}$ its constant value in $K\in\mathcal{T}_{h}$. Moreover, we shall assume that the reaction coefficient $\kappa$ varies slowly between neighbouring elements in the sense that that there exists a constant $C>0$ such that the following condition holds for all triangulations $\mathcal{T}_{h}\in\mathcal{G}$ and all elements $K\in\mathcal{T}_{h}$:

[TABLE]

We state, without proof, some elementary consequences of the above assumption:

Lemma 1**.**

Suppose that condition (4) holds. Then

if $\kappa_{K}=0$, then $\kappa_{K^{\prime}}<1/h_{K^{\prime}}$ on the patch $\widetilde{K}$; 2. 2.

if $h_{K}\kappa_{K}>1$, then $\kappa_{K^{\prime}}>0$ on the patch $\widetilde{K}$; 3. 3.

there exists a constant $C>0$ such that for all $\mathcal{T}_{h}\in\mathcal{G}$, all $K\in\mathcal{T}_{h}$, if $h_{K}\kappa_{K}>1$, then

[TABLE]

for all $K^{\prime}\in\widetilde{K}$; 4. 4.

there exists a constant $C>0$ such that for all $\mathcal{T}_{h}\in\mathcal{G}$, all $K\in\mathcal{T}_{h}$, and all elements $K^{\prime}\subset\widetilde{K}$,

[TABLE]

The quantity $\min\{h_{K},\kappa_{K}^{-1}\}$ appears extensively throughout the paper and we shall adopt the convention whereby

[TABLE]

2.3 Finite Element Discretization

The weak formulation of problem (1) reads: find $u\in V=\{v\in H^{1}(\Omega):v=0\text{ on }{\Gamma_{\mathrm{D}}}\}$ such that

[TABLE]

where $\mathcal{B}:V\times V\to\mathbb{R}$ and $\mathcal{F}:V\to\mathbb{R}$ are defined by

[TABLE]

It will be useful to introduce local counterparts of these forms

[TABLE]

The associated global and local energy norms $|\!|\!|\cdot|\!|\!|$ and $|\!|\!|\cdot|\!|\!|_{K}$ are defined by $|\!|\!|v|\!|\!|^{2}=\mathcal{B}(v,v)$ and $|\!|\!|v|\!|\!|_{K}^{2}=\mathcal{B}_{K}(v,v)$, respectively.

Let $V_{h}=\{v_{h}\in V:v_{h}|_{K}\in\mathbb{P}_{1}(K)\ \forall K\in\mathcal{T}_{h}\}$, where $\mathbb{P}_{1}(K)$ is the space of affine functions on $K$, then the finite element approximation $u_{h}\in V_{h}$ of (1) is defined by

[TABLE]

3 A Posteriori Error Estimator

Every partition $\mathcal{T}_{h}\in\mathcal{G}$ can be split into disjoint subsets $\mathcal{T}_{h}^{+}=\{K\in\mathcal{T}_{h}:\kappa_{K}>0\}$ and $\mathcal{T}_{h}^{0}=\{K\in\mathcal{T}_{h}:\kappa_{K}=0\}$. Let $\boldsymbol{\tau}\in\boldsymbol{H}(\operatorname{div},\Omega)$ be any vector field satisfying the conditions

[TABLE]

Let $\boldsymbol{\varepsilon}=\boldsymbol{\tau}-\boldsymbol{\nabla}u_{h}$ in $\Omega$, $r=\Pi f-\kappa^{2}u_{h}+\operatorname{div}\boldsymbol{\tau}$ in $\Omega$, and $R_{\mathrm{N}}=\Pi_{\mathrm{N}}g_{\mathrm{N}}-\boldsymbol{\tau}\cdot{\boldsymbol{\nu}}$ on ${\Gamma_{\mathrm{N}}}$, then the local error indicator over an element $K\in\mathcal{T}_{h}$ is defined to be

[TABLE]

Observe that $r$ and $R_{\mathrm{N}}$ vanish if $K\in\mathcal{T}_{h}^{0}$ and that the second and the third term is taken to be zero on such elements. The error estimator is then defined by

[TABLE]

where the oscillation term is given by

[TABLE]

and the constants

[TABLE]

arose in the corrigendum of [4, Lemma 1].

The following result, based on [4, Lemma 2], shows that the estimator provides an upper bound on the error:

Theorem 2**.**

Let $u_{h}\in V$ be arbitrary. If $\boldsymbol{\tau}\in\boldsymbol{H}(\operatorname{div},\Omega)$ satisfies equilibration conditions (8)–(9) then

[TABLE]

Proof.

The weak formulation (6) and the divergence theorem yield identity

[TABLE]

for all $v\in V$. The last two terms are estimated in the same way as in the proof of [4, Lemma 2]:

[TABLE]

For elements $K\in\mathcal{T}_{h}^{+}$ we bound

[TABLE]

where the trace inequality [4, Lemma 1] is employed.

Due to equilibration conditions (8)–(9), we arrive at

[TABLE]

Cauchy–Schwarz inequality, notation (11), and choice $v=u-u_{h}$ finish the proof. ∎

It will not have escaped the reader’s notice that nothing in the above argument relies on $u_{h}$ being a finite element approximation. Consequently, the upper bound presented in Theorem 2 holds true for arbitrary conforming approximation $u_{h}\in V$. However, the local efficiency and robustness results proved in Theorem 11 will require $u_{h}$ to be a Galerkin finite element approximation exactly satisfying the condition (7).

4 Flux Reconstruction by Patchwise Minimization

Let $\mathcal{N}_{h}$ denote the nodes in the partition $\mathcal{T}_{h}$. In particular, given a node ${\boldsymbol{n}}\in\mathcal{N}_{h}$, the subset $\mathcal{T}_{\boldsymbol{n}}=\{K\in\mathcal{T}_{h}:{\boldsymbol{n}}\in K\}$ consists of elements that touch the node, while $\mathcal{E}^{\mathrm{N}}_{\boldsymbol{n}}=\{\gamma\subset{\Gamma_{\mathrm{N}}}:{\boldsymbol{n}}\in\gamma\}$ consists of facets on the Neumann boundary ${\Gamma_{\mathrm{N}}}$ which touch the node ${\boldsymbol{n}}$.

The flux reconstructions used in the current work are constructed over patch $\omega_{\boldsymbol{n}}=\operatorname{int}\bigcup\mathcal{T}_{\boldsymbol{n}}$. Specifically, let

[TABLE]

where $\mathcal{E}^{\mathrm{E}}_{\boldsymbol{n}}=\{\gamma\subset\partial\omega_{\boldsymbol{n}}:{\boldsymbol{n}}\not\in\gamma\}$, ${\boldsymbol{\nu}}_{\boldsymbol{n}}$ denotes the unit outward facing normal vector on the boundary of the patch $\omega_{\boldsymbol{n}}$, and $\mathbf{RT}_{1}(K)=[\mathbb{P}_{1}(K)]^{d}\oplus\boldsymbol{x}\mathbb{P}_{1}(K)$ is the standard Raviart–Thomas space.

Let $\boldsymbol{\tau}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ denote the minimizer of the quadratic functional

[TABLE]

over $\boldsymbol{\tau}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ satisfying constraints

[TABLE]

where $\theta_{\boldsymbol{n}}$ is the usual piecewise affine and continuous hat function satisfying $\theta_{\boldsymbol{n}}({\boldsymbol{n}}^{\prime})=\delta_{{\boldsymbol{n}}{\boldsymbol{n}}^{\prime}}$ for all nodes ${\boldsymbol{n}}^{\prime}\in\mathcal{N}_{h}$, $\omega_{\boldsymbol{n}}^{+}=\operatorname{int}\bigcup(\mathcal{T}_{\boldsymbol{n}}\cap\mathcal{T}_{h}^{+})$, $\omega_{\boldsymbol{n}}^{0}=\operatorname{int}\bigcup(\mathcal{T}_{\boldsymbol{n}}\cap\mathcal{T}_{h}^{0})$, ${\Gamma_{\boldsymbol{n}}^{\mathrm{N}+}}=\bigcup\{\gamma\in\mathcal{E}^{\mathrm{N}}_{\boldsymbol{n}}:\kappa_{K_{\gamma}}>0\}$, $K_{\gamma}$ is the element adjacent to the facet $\gamma$, ${\Gamma_{\boldsymbol{n}}^{\mathrm{N}0}}=\bigcup\{\gamma\in\mathcal{E}^{\mathrm{N}}_{\boldsymbol{n}}:\kappa_{K_{\gamma}}=0\}$, and $C_{\mathrm{T}}^{\mathrm{N}}$ stands for piecewise constant function over facets defined as $C_{\mathrm{T}}^{\mathrm{N}}|_{\gamma}=C_{\mathrm{T}}^{K_{\gamma},\gamma}$ for all facets $\gamma\subset{\Gamma_{\mathrm{N}}}$ such that $\kappa_{K_{\gamma}}>0$.

The minimizer of (18) satisfying constraints (19)–(20) could, equally well, be characterised as the unique solution of the following problem: Find $\boldsymbol{\tau}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ and Lagrange multipliers $q_{h}\in\mathbb{P}_{1}^{*}(\omega_{\boldsymbol{n}}^{0})$ and $d_{h}\in\mathbb{P}_{1}^{*}({\Gamma_{\boldsymbol{n}}^{\mathrm{N}0}})$ satisfying

[TABLE]

for all $\boldsymbol{w}_{h}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ and

[TABLE]

where $\mathbb{P}_{1}^{*}(\omega_{\boldsymbol{n}}^{0})$ and $\mathbb{P}_{1}^{*}({\Gamma_{\boldsymbol{n}}^{\mathrm{N}0}})$ are spaces of discontinuous and piecewise affine functions over $\omega_{\boldsymbol{n}}^{0}$ and ${\Gamma_{\boldsymbol{n}}^{\mathrm{N}0}}$, respectively.

The condition in definition (17) imposed on the facets $\mathcal{E}^{\mathrm{N}}_{\boldsymbol{n}}$ means that $\boldsymbol{\tau}_{\boldsymbol{n}}$ can be extended by zero onto $\Omega$ thereby obtaining a vector field in $\boldsymbol{H}(\operatorname{div},\Omega)$, which we again denote by $\boldsymbol{\tau}_{\boldsymbol{n}}$. With this convention in place, the reconstructed flux $\boldsymbol{\tau}\in\boldsymbol{H}(\operatorname{div},\Omega)$ is taken to be the sum

[TABLE]

The resulting globally defined vector field $\boldsymbol{\tau}$ can be used in (12) to obtain an upper bound on the energy norm of the error, because it satisfies the equilibration conditions as stated in the following lemma.

Lemma 3**.**

*Reconstructed flux $\boldsymbol{\tau}\in\boldsymbol{H}(\operatorname{div},\Omega)$ given by (24) satisfies equilibration conditions (8)–(9). *

Proof.

Let $K\in\mathcal{T}_{h}^{0}$. Equality (22), definition (24), and partition of unity $\theta_{\boldsymbol{n}}$ yield

[TABLE]

for all $\varphi_{h}\in\mathbb{P}_{1}(K)$. Since $\Pi_{K}f-\kappa_{K}^{2}u_{h}|_{K}+\operatorname{div}\boldsymbol{\tau}|_{K}\in\mathbb{P}_{1}(K)$, equilibration condition (8) follows.

Similarly, given a facet $\gamma\subset{\Gamma_{\mathrm{N}}}\cap\partial K$, the equality (23) implies

[TABLE]

for all $\psi_{h}\in\mathbb{P}_{1}(\gamma)$. Equilibration condition (9) then follows, because $\Pi_{\gamma}g_{\mathrm{N}}-\boldsymbol{\tau}\cdot{\boldsymbol{\nu}}|_{\gamma}\in\mathbb{P}_{1}(\gamma)$. ∎

The next two sections are concerned with showing that reconstructed flux $\boldsymbol{\tau}$ defined in (24) yields locally efficient and robust error indicators. The main idea used in the proof is based on comparing $\boldsymbol{\tau}$ with two judiciously chosen flux reconstructions $\boldsymbol{\sigma}^{(1)}_{\widetilde{K}}$ and $\boldsymbol{\sigma}^{(2)}_{\widetilde{K}}$. Each of these reconstructions is defined in the neighbourhood $\widetilde{K}$ of the element $K$, see (3). While the flux reconstruction $\boldsymbol{\sigma}^{(1)}_{\widetilde{K}}$ is based on equilibrated interface fluxes $g_{K}$ introduced in [1] and analysed in [4], the second reconstruction $\boldsymbol{\sigma}^{(2)}_{\widetilde{K}}$ is new.

5 The First Auxiliary Flux Reconstruction

In this section we introduce the auxiliary flux reconstruction $\boldsymbol{\sigma}^{(1)}_{\widetilde{K}}$ and prove its properties. This flux reconstruction is based on equilibrated interface fluxes. If $u_{h}\in V_{h}$ is the Galerkin solution given by (7) then results of [2] guarantee the existence of interface fluxes $g_{K}$ satisfying for all elements $K\in\mathcal{T}_{h}$ the following properties

[TABLE]

and equilibration condition

[TABLE]

where $\mathcal{N}_{K}$ stands for the set of $d+1$ vertices of $K$. These fluxes do not yield robust a posteriori error estimators for large values of the reaction coefficient $\kappa$, as it was shown in [1]. However, we will use them only in elements where $h_{K}\kappa_{K}\leq 1$.

Below, we will utilize two estimates from [4]. First, quantity $R=g_{K}-\boldsymbol{\nabla}u_{h}\cdot{\boldsymbol{\nu}}_{K}$ defined on $\partial K$ for all elements $K\in\mathcal{T}_{h}$ satisfies

[TABLE]

see [4, estimate (31)]. Second, residual $r_{h}=\Pi f-\kappa^{2}u_{h}+\Delta u_{h}$ is bounded as

[TABLE]

see [4, estimate (29)] and also [1, Lemma 5].

The definition of the first auxiliary flux reconstruction proceeds as follows. If an element $K\in\mathcal{T}_{h}$ is such that $\kappa_{K}h_{K}\leq 1$, then we define

[TABLE]

where $\boldsymbol{\sigma}_{\boldsymbol{n}}^{(1)}$ is given piecewise as

[TABLE]

and vector fields $\boldsymbol{\sigma}^{(1)}_{{\boldsymbol{n}},K}$ are determined by the following lemma.

Lemma 4**.**

Let $u_{h}\in V_{h}$ satisfies (7). Let $K\in\mathcal{T}_{h}$ be an element and let ${\boldsymbol{n}}\in\mathcal{N}_{K}$ be its vertex. Then there exists $\boldsymbol{\sigma}^{(1)}_{{\boldsymbol{n}},K}\in\mathbf{RT}_{1}(K)$ such that

[TABLE]

and

[TABLE]

Proof.

To establish the existence and uniqueness of $\boldsymbol{\sigma}^{(1)}_{{\boldsymbol{n}},K}$ we recall [7] that $\mathbf{RT}_{1}(K)$ is unisolvent with respect to the degrees of freedom defined by

[TABLE]

and

[TABLE]

Observing that $\nabla:\mathbb{P}_{1}(K)/\mathbb{R}\to\mathbb{P}_{0}^{d}(K)$ is surjective, we may rewrite the former set of degrees of freedom in the equivalent form

[TABLE]

which, on integrating by parts and using the second set of degrees of freedom, shows that $\mathbf{RT}_{1}(K)$ is unisolvent with respect to the degrees of freedom defined by (34) augmented with the following

[TABLE]

Since the data in conditions (32) belong to $\mathbb{P}_{1}(\gamma)$ and since the equilibration condition (27) implies the following compatibility condition

[TABLE]

we deduce that $\boldsymbol{\sigma}^{(1)}_{{\boldsymbol{n}},K}$ exists and is unique. ∎

The following lemma shows that vector fields $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}$ lie in $\boldsymbol{W}(\omega_{\boldsymbol{n}})$.

Lemma 5**.**

Let $u_{h}\in V_{h}$ satisfies (7). Let ${\boldsymbol{n}}\in\mathcal{N}_{h}$ be a vertex and let $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}$ be defined by (31). Then $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ and $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}\cdot{\boldsymbol{\nu}}_{\boldsymbol{n}}=\Pi_{\gamma}(\theta_{\boldsymbol{n}}\Pi_{\gamma}g_{N})$ on facets $\gamma\in\mathcal{E}^{\mathrm{N}}_{\boldsymbol{n}}$.

Proof.

The fact that $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}\in\boldsymbol{H}(\operatorname{div},\omega_{\boldsymbol{n}})$ follows from the continuity of its normal components over element interfaces. Indeed, if $\gamma\subset\partial K\cap\partial K^{\prime}$ for elements $K,K^{\prime}\in\mathcal{T}_{\boldsymbol{n}}$ is an interior facet then

[TABLE]

by (26). Similarly, we verify the boundary conditions on $\partial\omega_{\boldsymbol{n}}$ required in (17). Clearly, $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}\cdot{\boldsymbol{\nu}}_{\boldsymbol{n}}=\Pi_{\gamma}(\theta_{\boldsymbol{n}}g_{K})=0$ on facets $\gamma\in\mathcal{E}^{\mathrm{E}}_{\boldsymbol{n}}$, because $\theta_{\boldsymbol{n}}=0$ on $\gamma$, and $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}\cdot{\boldsymbol{\nu}}_{\boldsymbol{n}}=\Pi_{\gamma}(\theta_{\boldsymbol{n}}\Pi_{\gamma}g_{N})$ on facets $\gamma\in\mathcal{E}^{\mathrm{N}}_{\boldsymbol{n}}$ by (25). ∎

An interesting consequence of Lemma 5 and identity (33) is that

[TABLE]

The vanishing normal components of $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}$ on exterior facets $\gamma\in\mathcal{E}^{\mathrm{E}}_{\boldsymbol{n}}$ guarantee that $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}$ can be extended by zero and as such belongs to $\boldsymbol{H}(\operatorname{div},\Omega)$. Consequently, $\boldsymbol{\sigma}^{(1)}_{\widetilde{K}}$ defined in (30) lies in $\boldsymbol{H}(\operatorname{div},\Omega)$ as well. The following lemma shows that $\boldsymbol{\sigma}^{(1)}_{\widetilde{K}}$ satisfies equilibration conditions (8)–(9).

Lemma 6**.**

Let $u_{h}\in V_{h}$ satisfies (7). Let $K\in\mathcal{T}_{h}$ be a fixed element and let $\boldsymbol{\sigma}^{(1)}_{\widetilde{K}}\in\boldsymbol{H}(\operatorname{div},\Omega)$ be defined by (30). Then

[TABLE]

Proof.

To prove (36), we use definitions (33) and (31) to find that

[TABLE]

holds in $K$. This identity together with the partition of unity $\sum_{{\boldsymbol{n}}\in\mathcal{N}_{K}}\theta_{\boldsymbol{n}}=1$, definition (30), and properties of the projection $\Pi_{K}$ yields

[TABLE]

in $K$.

To prove (37), we consider $K\in\mathcal{T}_{h}$ to be an element adjacent to the Neumann boundary ${\Gamma_{\mathrm{N}}}$ and $\gamma\subset{\Gamma_{\mathrm{N}}}\cap\partial K$ to be its facet. On this $\gamma$ we clearly have

[TABLE]

by (30), (31), (32), (25), and the fact that $\sum_{{\boldsymbol{n}}\in\mathcal{N}_{\gamma}}\theta_{\boldsymbol{n}}=1$. Here, $\mathcal{N}_{\gamma}$ stands for the set of $d$ vertices of the facet $\gamma$. ∎

Now, we formulate and prove the main result of this section. For an element $K\in\mathcal{T}_{h}$, we introduce neighbourhood $\widetilde{\widetilde{K}}=\bigcup\{K^{\prime}\in\mathcal{T}_{h}:K^{\prime}\cap\widetilde{K}\neq 0\}$, where $\widetilde{K}$ is given by (3).

Theorem 7**.**

Let $K\in\mathcal{T}_{h}$ be an element where $\kappa_{K}h_{K}\leq 1$. Let ${\boldsymbol{n}}\in\mathcal{N}_{K}$ be its vertex and let $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}$ be defined by (31). Then

[TABLE]

Proof.

Let $\boldsymbol{\varepsilon}_{\boldsymbol{n}}=\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}-\theta_{\boldsymbol{n}}\boldsymbol{\nabla}u_{h}$ and $K^{\prime}\in\mathcal{T}_{\boldsymbol{n}}$ be an element. By (33), the quantity $\boldsymbol{\varepsilon}_{\boldsymbol{n}}$ satisfies

[TABLE]

where $r_{h}=\Pi_{K^{\prime}}f-\kappa_{K^{\prime}}^{2}u_{h}+\Delta u_{h}$ was introduced above (29). Consequently,

[TABLE]

Now, boundary conditions (32) imply

[TABLE]

and

[TABLE]

where $\gamma_{\boldsymbol{n}}$ stands for the facet of $K^{\prime}$ opposite to the vertex ${\boldsymbol{n}}$. Since quantity $(\left\|\operatorname{div}\boldsymbol{w}\right\|_{K^{\prime}}^{2}+\left\|\boldsymbol{w}\cdot{\boldsymbol{\nu}}_{K^{\prime}}\right\|_{\partial K^{\prime}}^{2})^{1/2}$ is a norm in the finite dimensional space $\mathbf{RT}_{1}(K^{\prime})$, we can use the scaling argument

[TABLE]

Thus, using $\boldsymbol{w}=\boldsymbol{\varepsilon}_{\boldsymbol{n}}|_{K^{\prime}}$, inequalities (41) and (43), we obtain

[TABLE]

Hence, estimates (29) and (28) applied in (44) yield

[TABLE]

Finally, the bound (39) follows by using the local quasi-uniformity of the mesh. ∎

6 The Second Auxiliary Flux Reconstruction

For elements $K\in\mathcal{T}_{h}$, where $h_{K}\kappa_{K}>1$, we define the second auxiliary flux reconstruction as

[TABLE]

where

[TABLE]

and $f_{\boldsymbol{n}}$ is the $L^{2}(\omega_{\boldsymbol{n}})$-orthogonal projection of $f$ onto the space $\mathbb{P}_{1}(\omega_{\boldsymbol{n}})$ of affine functions on $\omega_{\boldsymbol{n}}$. Note that $\boldsymbol{\sigma}_{\boldsymbol{n}}^{(2)}$ is supported in $\omega_{\boldsymbol{n}}$ and that it is continuous. To simplify the notation we introduce piecewise linear and discontinuous function $f_{\boldsymbol{n}}^{\kappa}$ in the patch $\omega_{\boldsymbol{n}}$ by the rule

[TABLE]

For completeness, we define $f_{\boldsymbol{n}}^{\kappa}=\Pi f$ in those patches $\omega_{\boldsymbol{n}}$, where $h_{K}\kappa_{K}\leq 1$ for all $K\in\mathcal{T}_{\boldsymbol{n}}$.

It is clear that $\boldsymbol{\sigma}^{(2)}_{\boldsymbol{n}}|_{K}\in\mathbf{RT}_{1}(K)$ for all $K\in\mathcal{T}_{\boldsymbol{n}}$ and that $\boldsymbol{\sigma}^{(2)}_{\boldsymbol{n}}\cdot{\boldsymbol{\nu}}_{K}=0$ on facets $\gamma\in\mathcal{E}^{\mathrm{E}}_{\boldsymbol{n}}$. Thus $\boldsymbol{\sigma}^{(2)}_{\boldsymbol{n}}$ can be extended by zero such that $\boldsymbol{\sigma}^{(2)}_{\boldsymbol{n}}\in\boldsymbol{H}(\operatorname{div},\Omega)$ and consequently $\boldsymbol{\sigma}^{(2)}_{\widetilde{K}}\in\boldsymbol{H}(\operatorname{div},\Omega)$.

Lemma 8**.**

Let $u_{h}\in V_{h}$ be arbitrary. Let $K\in\mathcal{T}_{h}$ be an element such that $h_{K}\kappa_{K}>1$. Let ${\boldsymbol{n}}\in\mathcal{N}_{K}$ be its vertex and let $\boldsymbol{\sigma}^{(2)}_{\boldsymbol{n}}$ be defined by (46). Then

[TABLE]

holds for all $K^{\prime}\in\mathcal{T}_{\boldsymbol{n}}$.

Proof.

Let $K^{\prime}\in\mathcal{T}_{\boldsymbol{n}}$ be fixed. Recall that Lemma 1 implies $\kappa_{K^{\prime}}>0$, $\kappa_{K}\leq C\kappa_{K^{\prime}}$ and $\kappa_{K^{\prime}}^{-1}h_{K^{\prime}}^{-1}\leq C$. The inverse inequality yields

[TABLE]

Similarly,

[TABLE]

where estimate (48) and the fact that $\Delta u_{h}|_{K^{\prime}}=0$ were used. Consequently, using bound $\kappa_{K^{\prime}}^{-1}h_{K^{\prime}}^{-1}\leq C$ and triangle inequality, we obtain

[TABLE]

where we employ (29). ∎

Lemma 9**.**

Let $u_{h}\in V_{h}$ be arbitrary. Let $K\in\mathcal{T}_{h}$ be an element such that $\kappa_{K}h_{K}>1$. Let ${\boldsymbol{n}}\in\mathcal{N}_{K}$ be its vertex and let $\gamma\in\mathcal{E}^{\mathrm{N}}_{\boldsymbol{n}}$ be a facet on ${\Gamma_{\mathrm{N}}}$ adjacent to an element $K_{\gamma}$. Then

[TABLE]

Proof.

Using the definition (46) of $\boldsymbol{\sigma}^{(2)}_{\boldsymbol{n}}$, properties of projection $\Pi_{\gamma}$ and hat function $\theta_{\boldsymbol{n}}$, we obtain

[TABLE]

To bound this norm, we consider a special test function. We define function $v$ on the boundary $\partial{K_{\gamma}}$ as $v=0$ on $\partial{K_{\gamma}}\setminus\gamma$ and $v=\beta\cdot(\Pi_{\gamma}g_{\mathrm{N}}-\kappa_{K}^{-2}\boldsymbol{\nabla}f_{\boldsymbol{n}}\cdot{\boldsymbol{\nu}})$ on $\gamma$, where $\beta=\prod_{{\boldsymbol{n}}\in\mathcal{N}_{\gamma}}\theta_{\boldsymbol{n}}$ is a bubble function defined on the facet $\gamma$. Then we introduce minimum energy extension $\mathcal{E}v$ to the interior of ${K_{\gamma}}$ satisfying $\mathcal{E}v\in H^{1}({K_{\gamma}})$, $\mathcal{E}v=v$ on $\partial{K_{\gamma}}$, and $\mathcal{B}_{K_{\gamma}}(\mathcal{E}v,w)=0$ for all $w\in H^{1}_{0}({K_{\gamma}})$, see [1, Section 3.1]. Extending $\mathcal{E}v$ further by zero to the rest of the domain $\Omega$, we have $\mathcal{E}v\in V$.

Using $\mathcal{E}v$ as a test function in (6) together with identity $\int_{\gamma}\boldsymbol{\nabla}f_{\boldsymbol{n}}\cdot{\boldsymbol{\nu}}v\,\mathrm{d}\boldsymbol{s}=\int_{K_{\gamma}}\boldsymbol{\nabla}f_{\boldsymbol{n}}\cdot\boldsymbol{\nabla}\mathcal{E}v\,\mathrm{d}\boldsymbol{x}$, we derive

[TABLE]

Since $\Pi_{\gamma}g_{\mathrm{N}}-\kappa_{K}^{-2}\boldsymbol{\nabla}f_{\boldsymbol{n}}\cdot{\boldsymbol{\nu}}\in\mathbb{P}_{1}(\gamma)$ and $\left(\int_{\gamma}\beta\varphi^{2}\,\mathrm{d}\boldsymbol{x}\right)^{1/2}$ is a norm in $\mathbb{P}_{1}(\gamma)$, we use the equivalence of norms in the finite dimensional space $\mathbb{P}_{1}(\gamma)$ to get

[TABLE]

This estimate together with identity (51), Cauchy–Schwarz inequality, inverse inequality, and bound $h_{K_{\gamma}}^{-1}\kappa_{K_{\gamma}}^{-1}\leq C$ provided by assumption (4) yields

[TABLE]

Here, the energy norm $|\!|\!|\mathcal{E}v|\!|\!|_{K_{\gamma}}$ is bounded by [1, Lemma 4] as

[TABLE]

and norm $\left\|v\right\|_{\gamma}$ by using the equivalence of norms in $\mathbb{P}_{1}(\gamma)$ as

[TABLE]

Finally, using these bounds, inequalities $h_{K_{\gamma}}^{-1}\kappa_{K_{\gamma}}^{-1}\leq C$ and $C_{\mathrm{T}}^{K,\gamma}\leq C\kappa_{K_{\gamma}}^{-1/2}$, and estimate (29) in (52), we obtain

[TABLE]

This estimate and (50) finish the proof. ∎

7 Efficiency of Patchwise Minimizations

We first formulate a lemma stating that error indicators can be bounded by the value of the quadratic functional $E_{\boldsymbol{n}}$.

Lemma 10**.**

Let $u_{h}\in V$ be arbitrary. Let $d$ stand for the dimension. Let error indicators $\eta_{K}$ be defined by (10). Let reconstructed flux $\boldsymbol{\tau}\in\boldsymbol{H}(\operatorname{div},\Omega)$ given by (24) satisfy equilibration conditions (8)–(9) and let its local components $\boldsymbol{\tau}_{\boldsymbol{n}}$ be in $\boldsymbol{W}(\omega_{\boldsymbol{n}})$. Then

[TABLE]

Proof.

Let simplex $K\in\mathcal{T}_{h}^{+}$ be fixed. Since there is $d+1$ facets on the boundary of $K$, we can bound $\eta_{K}(\boldsymbol{\tau})$ by Cauchy–Schwarz inequality as

[TABLE]

Using the partition of unity $\theta_{\boldsymbol{n}}$, definition of $\boldsymbol{\varepsilon}$, definition (24) of $\boldsymbol{\tau}$, and Cauchy–Schwarz inequality, we obtain

[TABLE]

Similarly, we estimate

[TABLE]

and for facets $\gamma\subset{\Gamma_{\mathrm{N}}}\cap\partial K$

[TABLE]

Statement (53) now follows by a combination of these estimates.

If $K\in\mathcal{T}_{h}^{0}$ then bound (53) is easy to verify, because of identity $\eta_{K}(\boldsymbol{\tau})=\left\|\boldsymbol{\varepsilon}\right\|_{K}$ and estimate (54). ∎

Note that this lemma holds true for any flux $\boldsymbol{\tau}\in\boldsymbol{H}(\operatorname{div},\Omega)$ given by (24). Its local components $\boldsymbol{\tau}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ are not required to minimize the quadratic functional $E_{\boldsymbol{n}}$ defined in (18).

The following theorem presents the main result of this paper. It states the efficiency and robustness of error indicators computed by (10) from the patchwise flux reconstruction $\boldsymbol{\tau}\in\boldsymbol{H}(\operatorname{div},\Omega)$ defined in (24). To formulate it, we introduce the union $\widetilde{\Gamma}_{\mathrm{N}}^{K}$ of facets in the triangulation $\mathcal{T}_{h}$ that lie on ${\Gamma_{\mathrm{N}}}$ and have at least one common point with $K$, i.e., $\widetilde{\Gamma}_{\mathrm{N}}^{K}=\bigcup\{\gamma\subset{\Gamma_{\mathrm{N}}}:\gamma\cap K\neq\emptyset\}$.

Theorem 11**.**

Let $u\in V$ be the weak solution (6) and let $u_{h}\in V_{h}$ be its Galerkin approximation satisfying (7). Let flux reconstruction $\boldsymbol{\tau}\in\boldsymbol{H}(\operatorname{div},\Omega)$ be given by (24) and let its local components $\boldsymbol{\tau}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ solve local problems (21)–(23). Then there exists a constant $C>0$ independent of reaction coefficient $\kappa$ and any mesh size such that the local efficiency estimate

[TABLE]

holds true for all elements $K\in\mathcal{T}_{h}$.

Proof.

We consider two cases. First, let $K\in\mathcal{T}_{h}$ be such that $\kappa_{K}h_{K}\leq 1$. Since $\boldsymbol{\tau}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ minimizes the functional $E_{\boldsymbol{n}}$, both $\boldsymbol{\tau}_{\boldsymbol{n}}$ and $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}$ satisfy constraints (19)–(20), and fluxes $\boldsymbol{\sigma}^{(1)}_{\boldsymbol{n}}$ satisfy (35), we obtain from (53) the estimate

[TABLE]

Applying Theorem 7 in this inequality immediately yields

[TABLE]

Second, let $K\in\mathcal{T}_{h}$ be such that $\kappa_{K}h_{K}>1$. By assumption (4), the reaction coefficient satisfies $\kappa_{K^{\prime}}>0$ for all elements $K^{\prime}\subset\widetilde{K}$. Since $\boldsymbol{\sigma}^{(2)}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ and $\boldsymbol{\tau}_{\boldsymbol{n}}$ is the minimizer of $E_{\boldsymbol{n}}$, the inequality (53) yields

[TABLE]

Estimates (47) and (49) then give

[TABLE]

Assumption (4) and a combination of (56) and (57) finishes the proof. ∎

The oscillation term $\min\{h_{K},\kappa_{K}^{-1}\}\left\|f_{\boldsymbol{n}}^{\kappa}-\Pi f\right\|_{\omega_{\boldsymbol{n}}}$ is not standard, however, as well as the other oscillation terms in (55) it is of higher order than the error $|\!|\!|u-u_{h}|\!|\!|_{\widetilde{\widetilde{K}}}$ and does not spoil the robust local efficiency result.

8 Avoiding Equilibration

Theorem 2 requires the flux $\boldsymbol{\tau}$ to satisfy equilibration conditions (8)–(9). Therefore, the local minimization problems (18) is constrained by (19)–(20). However, in practical computations constraints (19)–(20) are often not satisfied exactly due to round-off errors. Consequently, assumptions of Theorem 2 are not valid and error estimator $\eta(\boldsymbol{\tau})$ is not guaranteed to provide the upper bound on the error. This problem can be avoided by introducing Friedrichs–Poincaré and trace inequalities and two additional parameters in the definition of the estimator.

Given arbitrary flux $\boldsymbol{\tau}\in\boldsymbol{H}(\operatorname{div},\Omega)$, we define modified local error indicators

[TABLE]

where small parameter $\kappa_{0}>0$ replaces the zero value of $\kappa$ in a sense and small parameter $\zeta_{0}>0$ has a similar meaning. Notice that $\tilde{\eta}_{K}(\boldsymbol{\tau})$ differs from $\eta_{K}(\boldsymbol{\tau})$ only for elements $K\in\mathcal{T}_{h}^{0}$.

The modified error estimator is then defined as

[TABLE]

where $C_{\mathrm{FP}}>0$ and $C_{\mathrm{T}}>0$ are constants from Friedrichs–Poincaré and trace inequalities

[TABLE]

The following theorem presents a modification of Theorem 2 that avoids the equilibration conditions (8)–(9).

Theorem 12**.**

Let $u_{h}\in V$ and $\boldsymbol{\tau}\in\boldsymbol{H}(\operatorname{div},\Omega)$ be arbitrary. Then

[TABLE]

for all $\kappa_{0}>0$ and $\zeta_{0}>0$.

Proof.

Using identity (13) and estimates (14)–(16), we arrive at

[TABLE]

Using Cauchy–Schwarz inequality

[TABLE]

we obtain

[TABLE]

Friedrichs–Poincaré and trace inequalities (60), notation (59), and choice $v=u-u_{h}$ finish the proof. ∎

Constants $\kappa_{0}$ and $\zeta_{0}$ should be small. Ideally so small that $1+\kappa_{0}^{2}C_{\mathrm{FP}}^{2}+\zeta_{0}^{2}C_{\mathrm{T}}^{2}\approx 1$. In this case the influence of Friedrichs–Poincaré and trace constants $C_{\mathrm{FP}}$ and $C_{\mathrm{T}}$ on the value of the error bound (61) is negligible. In the case of pure Dirichlet boundary conditions, i.e., ${\Gamma_{\mathrm{N}}}=\emptyset$, the parameter $\zeta_{0}$ is not needed and estimate (61) holds with $C_{\mathrm{T}}=0$. In case $\kappa>0$ everywhere in $\Omega$ the set $\mathcal{T}_{h}^{0}$ is empty and parameters $\kappa_{0}$, $\zeta_{0}$ and constants $C_{\mathrm{FP}}$, $C_{\mathrm{T}}$ are not needed. Estimate (61) then holds with $C_{\mathrm{FP}}=C_{\mathrm{T}}=0$ and local error indicators (10) and (58) coincide. However, if $\kappa$ vanishes at some parts of $\Omega$ and constants $C_{\mathrm{FP}}$ and $C_{\mathrm{T}}$ are needed, then they can be computed analytically in some special cases and numerically, in general. Even their guaranteed numerical bounds are available, see e.g. [28, 31].

Since Theorem 12 does not require any equilibration condition, the patchwise flux reconstruction procedure simplifies. Modified fluxes $\tilde{\boldsymbol{\tau}}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ minimize the quadratic functional

[TABLE]

over the space $\boldsymbol{W}(\omega_{\boldsymbol{n}})$, where ${\Gamma_{\boldsymbol{n}}^{\mathrm{N}}}$ is the union of the facets belonging to $\mathcal{E}^{\mathrm{N}}_{\boldsymbol{n}}$ and piecewise constant parameters $\tilde{\kappa}$ and $\tilde{\zeta}$ are given by

[TABLE]

for all elements $K\in\mathcal{T}_{h}$ and all facets $\gamma\subset{\Gamma_{\mathrm{N}}}$, where we recall that $K_{\gamma}$ denotes the element adjacent to the facet $\gamma$.

The minimizer $\tilde{\boldsymbol{\tau}}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ of (62) could, equally well, be characterised as the unique solution of the following problem:

[TABLE]

for all $\boldsymbol{w}_{h}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$. Note that the large values $\kappa_{0}^{-1}$ and $\zeta_{0}^{-1}$ play here the role of penalty parameters to impose constraints (19)–(20) and (22)–(23) in a weak sense. Consequently, the difference between $\boldsymbol{\tau}$ and $\tilde{\boldsymbol{\tau}}$ is small in practical computations.

Summing up fluxes $\tilde{\boldsymbol{\tau}}_{\boldsymbol{n}}$ as in Section 4 results in a modified reconstructed flux

[TABLE]

that can be directly used in Theorem 12 to obtain a guaranteed upper bound on the error.

The modified reconstructed flux is locally efficient and robust as it is stated in the following corollary.

Corollary 13**.**

Let $u\in V$ be the weak solution (6) and let $u_{h}\in V_{h}$ be its Galerkin approximation satisfying (7). Let flux reconstruction $\tilde{\boldsymbol{\tau}}\in\boldsymbol{H}(\operatorname{div},\Omega)$ be given by (65) and let its local components $\tilde{\boldsymbol{\tau}}_{\boldsymbol{n}}\in\boldsymbol{W}(\omega_{\boldsymbol{n}})$ solve local problems (64). Then there exists a constant $C>0$ independent of reaction coefficient $\kappa$ and any mesh size such that the local efficiency estimate

[TABLE]

holds true for all elements $K\in\mathcal{T}_{h}$.

Proof.

The proof follows the same lines as the proof of Theorem 11. In particular, we use the fact that

[TABLE]

∎

9 Numerical Examples

Example 1

In this example, we consider problem (1) in a domain with reentrant corner: $\Omega=\{(\varrho,\phi):0\leq\varrho<1\text{ and }\phi\in(\pi/2,2\pi)\}$, where $\varrho$ and $\phi$ are standard polar coordinates, see Figure 1 (left). The boundary conditions are homogeneous Dirichlet only, i.e., ${\Gamma_{\mathrm{D}}}=\partial\Omega$ and ${\Gamma_{\mathrm{N}}}=\emptyset$. The reaction coefficient $\kappa$ is assumed positive, constant in $\Omega$, and its specific values are provided below. Choosing the right-hand side as $f=\kappa^{2}\varrho^{2/3}\sin(2\phi-\pi)/3$, the exact solution is explicitly given by

[TABLE]

where $I_{\alpha}$ stands for the modified Bessel function of the first kind. This solution exhibits singularity at the origin and a boundary layer at $\varrho=1$ for large values of $\kappa$.

We first compute the finite element solution (7) using the mesh shown in Figure 1 (right) for $\kappa=10^{-3}$, $10^{-2},\dots,10^{6}$. For each value of $\kappa$ we compute flux reconstruction (24) by solving local problems (21)–(23) and evaluate the error estimator $\eta(\boldsymbol{\tau})$ given by (11). Note that since ${\Gamma_{\mathrm{N}}}=\emptyset$ and $\kappa>0$, the procedure considerably simplifies. The set $\mathcal{T}_{h}^{0}$ is empty, equilibration conditions (8)–(9) do not apply as well as constraints (22)–(23). Reconstructed fluxes $\boldsymbol{\tau}$ and $\tilde{\boldsymbol{\tau}}$ given by (24) and (65), respectively, are identical and $\eta_{K}(\boldsymbol{\tau})=\tilde{\eta}_{K}(\tilde{\boldsymbol{\tau}})$ for all $K\in\mathcal{T}_{h}$. In particular constants $\kappa_{0}$, $\zeta_{0}$, $C_{\mathrm{FP}}$, and $C_{\mathrm{T}}$ are not needed.

Figure 2 (left) presents the index of effectivity

[TABLE]

for the chosen values of $\kappa$. All values of $I_{\text{eff}}$ are above 1 confirming that $\eta(\boldsymbol{\tau})$ is the guaranteed upper bound on the error. On the other hand they are not far from 1 in the whole range of values of $\kappa$ showing the robust efficiency. All these indices of effectivity are below 1.12, which illustrates high accuracy of computed error estimators. For comparison, we also present indices of effectivity for the error estimator proposed in our previous work [4], see the dashed lines in Figure 2. Its accuracy legs behind the current approach.

To illustrate the robustness with respect to the mesh size, we also solve this problem on a sequence of uniformly refined meshes for a fixed value $\kappa=100$ and plot the resulting indices of effectivity in Figure 2 (right). In this case we observe robust efficiency and high accuracy as well.

Error indicators $\eta_{K}(\boldsymbol{\tau})$ given in (10) can be utilized for adaptive mesh refinement and error estimator $\eta(\boldsymbol{\tau})$ for a guaranteed stopping criterion. We use the standard adaptive algorithm: SOLVE – ESTIMATE – STOP – MARK – REFINE. Given an initial mesh, the SOLVE step computes the finite element solution by (7), the ESTIMATE step evaluates the flux reconstruction $\boldsymbol{\tau}$ defined by (24) and error indicators $\eta_{K}(\boldsymbol{\tau})$ introduced in (10). In the STOP step, the error estimator $\eta(\boldsymbol{\tau})$ given by (11) is computed and the algorithm is stopped if $\eta(\boldsymbol{\tau})$ (and consequently the error $|\!|\!|u-u_{h}|\!|\!|$) is below the required tolerance. In the MARK step, the Dörfler strategy [11] is used to mark elements, where $\eta_{K}(\boldsymbol{\tau})$ indicate large error. Finally, the longest edge bisection algorithm [22, 34] is applied in the REFINE step to refine the marked elements and create a new mesh.

Several examples of adaptively refined meshes are provided in Figure 3. The optimal speed of convergence of both the error $|\!|\!|u-u_{h}|\!|\!|$ and error estimator $\eta(\boldsymbol{\tau})$ during the adaptive algorithm is presented in Figure 4 (left). Figure 4 (right) shows corresponding indices of effectivity. They are all above and quite close to 1, confirming the robust efficiency of the error estimator even on highly graded meshes.

Example 2

This example illustrates the behaviour of the proposed error estimator for a discontinuous right-hand side $f$, piecewise constant reaction coefficient $\kappa$, and homogeneous Neumann boundary conditions. We consider problem (1) in a square $\Omega=(-1,1)^{2}$ with ${\Gamma_{\mathrm{N}}}=\partial\Omega$ and $g_{\mathrm{N}}=0$ on ${\Gamma_{\mathrm{N}}}$. Right-hand side $f$ equals to $\kappa^{2}$ in the disc $B_{1/2}=\{\varrho\leq 1/2\}$, where $\varrho$ is the distance from the origin, and it vanishes elsewhere. The exact solution of this problem is not known, but for large $\kappa$ it is supposed to be close to the characteristic function of the disc $B_{1/2}$ with a steep interior layer close to the boundary of $B_{1/2}$. Since the exact solution is not known, we approximate the true error by $u_{h}^{\mathrm{ref}}-u_{h}$, where the reference solution $u_{h}^{\mathrm{ref}}$ is computed by finite elements of order 5 on the same mesh as $u_{h}$.

We choose $\kappa=100$ and use the modified flux reconstruction (65) and the modified error estimator (59). Since $\Omega$ is a square, we can compute the Friedrichs–Poincaré and trace constants analytically. We use $C_{\mathrm{FP}}^{2}=2/\pi^{2}$ and $C_{\mathrm{T}}^{2}=\sqrt{2}\coth(\sqrt{2}/2)$. Parameters $\kappa_{0}$ and $\zeta_{0}$ are chosen as square roots of the machine epsilon: $\kappa_{0}=\zeta_{0}\approx 10^{-8}$.

We solve this problem by the adaptive algorithm described above starting with a mesh with two triangles. This setting does not satisfy assumptions listed in Subsection 2.2, because discontinuities in $\kappa$ are not compatible with the mesh. Therefore, for the purpose of computation, we use the value of $\kappa$ in the centroid of each element as the constant value in the element. In this way we construct certain approximate solution and the corresponding error estimator, which is guaranteed by Theorem 12 to be above the true error. The obtained indices of effectivity show robust and efficient performance of the estimator even in this case.

Figure 5 (left) shows the energy norm of the approximate error $|\!|\!|u_{h}^{\mathrm{ref}}-u_{h}|\!|\!|$, the computed error bound $\tilde{\eta}(\tilde{\boldsymbol{\tau}})$, and the oscillation term $\operatorname{osc}^{2}(f,g_{\mathrm{N}})=\sum_{K\in\mathcal{T}_{h}}\operatorname{osc}_{K}^{2}(f,g_{\mathrm{N}})$ during the adaptive process. Figure 5 (right) presents the corresponding indices of effectivity $I_{\text{eff}}=\tilde{\eta}(\tilde{\boldsymbol{\tau}})/|\!|\!|u_{h}^{\mathrm{ref}}-u_{h}|\!|\!|$. We may observe that the error bound is really above the error and that the error estimator estimates it robustly on all meshes. The oscillation term is of comparable size as the error at the beginning of the adaptive process, which leads to higher values of the index of effectivity. However, starting from meshes with around $10^{3}$ degrees of freedom the interior layer is well resolved, the oscillation term decreases faster than the error, and the index of effectivity decreases towards one. For illustration we present three adaptively refined meshes in Figure 6.

10 Conclusions

In this paper we present an a posteriori error estimator that is fully computable and provides a locally efficient upper bound on the energy norm of the error. This error estimator can be computed by a fast and easily parallelizable algorithm by solving small and independent problems on patches of elements. We proved its robustness both with respect to the mesh size and the reaction coefficient $\kappa$. We demonstrated by numerical examples that the corresponding local error indicators can be successfully used in the standard adaptive algorithm to guide the mesh adaptation and that the error estimator provides sharp results on rough, fine, and adaptively refined meshes as well as in the singularly perturbed case when $\kappa$ is large.

Further research questions about this error estimator may include its robustness for higher order finite element approximations [29] and its possible modifications to guarantee robustness on anisotropically refined meshes.

The proposed flux reconstruction can be used not only for the presented reaction-diffusion problems, but also for related eigenvalue problems. It was recently shown [31] that any flux reconstruction for boundary value problems can be directly used in the Lehmann–Goerisch method for guaranteed bounds on eigenvalues.