A posteriori error estimates for the mortar staggered DG method

TL;DR

This paper develops and proves the reliability and efficiency of two residual-type a posteriori error estimators for mortar staggered DG methods applied to second order elliptic equations, validated by numerical experiments.

Contribution

It introduces two novel residual-type error estimators for mortar staggered DG methods, avoiding saturation assumptions and validated through theoretical proofs and numerical tests.

Findings

Both error estimators are reliable and efficient.

The estimators do not require saturation assumptions.

Numerical experiments confirm the theoretical results.

Abstract

Two residual-type error estimators for the mortar staggered discontinuous Galerkin discretizations of second order elliptic equations are developed. Both error estimators are proved to be reliable and efficient. Key to the derivation of the error estimator in potential $L^2$ error is the duality argument. On the other hand, an auxiliary function is defined, making it capable of decomposing the energy error into conforming part and nonconforming part, which can be combined with the well-known Scott-Zhang local quasi-interpolation operator and the mortar discrete formulation yields an error estimator in energy error. Importantly, our analysis for both error estimators does not require any saturation assumptions which are often needed in the literature. Several numerical experiments are presented to confirm our proposed theories.

Equations

$$\begin{array}{rl}-\nabla \cdot (\rho \nabla u)& =f\text{mbox}in\Omega ,\\ u& =0\text{mbox}on\partial \Omega ,\end{array}$$

$$Q_i=L^2({\Omega }_i{)}^2,Q=\prod _{i=1}^N{Q}_i,M=L^2(\Gamma )$$

$$V_i=\{v\in H^1({\Omega }_i),v=0\text{mbox}on\partial {\Omega }_i\cap \partial \Omega \},V=\prod _{i=1}^N{V}_i.$$

$$\begin{array}{rl}\mathbf{z}& =\rho \nabla u\text{mbox}in\Omega ,\\ -\nabla \cdot \mathbf{z}& =f\text{mbox}in\Omega ,\\ u& =0\text{mbox}on\partial \Omega ,\end{array}$$

$$\begin{array}{rl}{\mathbf{z}}_i& ={\rho }_i\nabla u_i\text{mbox}in{\Omega }_i,\\ -\nabla \cdot {\mathbf{z}}_i& =f\text{mbox}in{\Omega }_i,\\ {\mathbf{z}}_i\cdot {\mathbf{n}}_i& ={\lambda }_i\text{mbox}on\partial {\Omega }_i\backslash \partial \Omega ,\\ u_i& =0\text{mbox}on\partial {\Omega }_i\cap \partial \Omega ,\\ {\lambda }_i+{\lambda }_j& =0\text{mbox}forall(i,j)\in S\end{array}$$

$$\lambda {\mid }_{{\Gamma }_{ij}}={\lambda }_i{\mathbf{n}}_i\cdot {\mathbf{n}}_{ij}={\lambda }_j{\mathbf{n}}_j\cdot {\mathbf{n}}_{ij}\forall {\Gamma }_{ij}\subset \Gamma .$$

$$\begin{array}{rl}{\rho }_i^{-1}(\mathbf{z},\mathbf{q}{)}_{{\Omega }_i}& =(\nabla u,\mathbf{q}{)}_{{\Omega }_i}\forall \mathbf{q}\in Q_i,\\ (\mathbf{z},\nabla v{)}_{{\Omega }_i}-(\lambda {\mathbf{n}}_{ij}\cdot {\mathbf{n}}_i,v{)}_{{\Gamma }_i}& =(f,v{)}_{{\Omega }_i}\forall v\in V_i,\\ \sum _{i=1}^N([[u]],\mu {)}_{{\Gamma }_i}& =0\forall \mu \in M,\end{array}$$

$$A(\mathbf{z},u,\lambda ;\mathbf{q},v,\mu )=(f,v)\forall (\mathbf{q},v,\mu )\in Q\times V\times M,$$

$$\begin{split}A(\bm{z},u,\lambda;\bm{q},v,\mu)=\sum_{i=1}^{N}\Big{(}\rho_{i}^{-1}(\bm{z},\bm{q})_{\Omega_{i}}-(\nabla u,\bm{q})_{\Omega_{i}}+(\bm{z},\nabla v)_{\Omega_{i}}-(\lambda,v\,\bm{n}_{i}\cdot\bm{n}_{ij})_{\Gamma_{i}}+(u\bm{n}_{i}\cdot\bm{n}_{ij},\mu)_{\Gamma_{i}}\Big{)}.\end{split}$$

$$Q_{h_i}=\{\mathbf{q}:\mathbf{q}{\mid }_{\tau }\in P^k(\tau {)}^2,\tau \in {\mathcal{T}}_i\text{mbox}and[[\mathbf{q}\cdot \mathbf{n}]]{\mid }_e=0,\forall e\in {\mathcal{F}}_{p,i}\}$$

$$V_{h_i}=\{v:v{\mid }_{\tau }\in P^k(\tau ),\tau \in {\mathcal{T}}_i\text{mbox}and[[v]]{\mid }_e=0,\forall e\in {\mathcal{F}}_{u,i}\},$$

$$[[\mathbf{q}\cdot \mathbf{n}]]{\mid }_e:=\mathbf{q}{\mid }_{{\tau }_1}\cdot \mathbf{n}-\mathbf{q}{\mid }_{{\tau }_2}\cdot \mathbf{n}\text{mbox}and[[v]]{\mid }_e:=v{\mid }_{{\tau }_1}-v{\mid }_{{\tau }_2}.$$

$$V_{h_i}^0=\{v\in V_{h_i}:v=0\text{mbox}on\partial \Omega \cap \partial {\Omega }_i\}.$$

$${\hat{V}}_h=\{v=(v_1,\dots ,v_N)\in V_h:{\int }_{{\Gamma }_{ij}}(v_i-v_j)\psi ds=0\forall \psi \in M_{ij},\forall {\Gamma }_{ij}\subset \Gamma \}.$$

$$c_i(v,\mu )=\sum _{{\Gamma }_{ij}\subset \partial {\Omega }_i}{\int }_{{\Gamma }_{ij}}v\mu {\mathbf{n}}_i\cdot {\mathbf{n}}_{ij}ds\forall (v,\mu )\in V_{h_i}^0\times \prod _{(i,j)\in S}{M}_{ij}.$$

$$b_i(\mathbf{z},v)$$

$$b_i^{\ast }(v,\mathbf{q})$$

$$\begin{array}{rl}{\rho }_i^{-1}({\mathbf{z}}_h,\mathbf{q}{)}_{{\Omega }_i}& =b_i^{\ast }(u_h,\mathbf{q})\forall \mathbf{q}\in Q_{h_i},i=1,\dots ,N,\\ b_i({\mathbf{z}}_h,v)-c_i(v,{\lambda }_h)& =(f,v{)}_{{\Omega }_i}\forall v\in V_{h_i}^0,i=1,\dots ,N,\\ \sum _{i=1}^N{c}_i(u_h,\mu )& =0\forall \mu \in M_h,\end{array}$$

$$A_h({\mathbf{z}}_h,u_h,{\lambda }_h;\mathbf{q},v,\mu )=(f,v)\forall (v,\mathbf{q},\mu )\in V_h\times Q_h\times M_h,$$

$$\displaystyle A_{h}(\bm{z}_{h},u_{h},\lambda_{h};\bm{q},v,\mu)=\sum_{i=1}^{N}\Big{(}\rho_{i}^{-1}(\bm{z},\bm{q})_{\Omega_{i}}-b_{i}^{*}(u_{h},\bm{q})+b_{i}(\bm{z}_{h},v)-c_{i}(v,\lambda_{h})+c_{i}(u_{h},\mu)\Big{)}.$$

$$A_h(\mathbf{z}-{\mathbf{z}}_h,u-u_h,\lambda -{\lambda }_h;\mathbf{q},v,\mu )=0\forall (v,\mathbf{q},\mu )\in V_h\times Q_h\times M_h.$$

$$\sum _{i=1}^N{\rho }_i\Vert \nabla (u-u_h){\Vert }_{0,{\Omega }_i}^2$$

$$\sum _{i=1}^N{\rho }_i\Vert u-u_h{\Vert }_{0,{\Omega }_i}^2$$

$$\sum _{i=1}^N{\rho }_i^{-1}\Vert \mathbf{z}-{\mathbf{z}}_h{\Vert }_{0,{\Omega }_i}^2$$

$$\Vert v{\Vert }_{0,e}\le C(h_{\tau }^{-\frac{1}{2}}\Vert v{\Vert }_{0,\tau }+h_{\tau }^{\frac{1}{2}}\Vert \nabla v{\Vert }_{0,\tau })\forall v\in H^1(\tau ),\forall \tau \in \mathcal{T},e\in \partial \tau .$$

$$\Vert v{\Vert }_{1/2,e}\le C\Vert v{\Vert }_{1,\tau }\forall v\in H^1(\tau ),\tau \in \mathcal{T},e\in \partial \tau .$$

$$(I_{i}v,\varphi {)}_{\tau }$$

$$(I_{i}v,\varphi {)}_e$$

$$(J_i\mathbf{\tau },\mathbf{q}{)}_{\tau }$$

$$(J_i\mathbf{\tau }\cdot {\mathbf{n}}_e,\varphi {)}_e$$

$$\begin{array}{rl}\Vert \mu -{\pi }_e\mu {\Vert }_{0,e}& \le C{h}_e^{1/2}\Vert \mu {\Vert }_{1/2,e},\\ \Vert v-I_{i}v{\Vert }_{0,\tau }& \le C{h}_{\tau }^{l+1}\Vert v{\Vert }_{l+1,\tau },l=0,1,\\ \Vert v-I_{i}v{\Vert }_{0,e}& \le C{h}_e^{3/2}\Vert v{\Vert }_{2,\tau },\\ \Vert \mathbf{q}-J_i\mathbf{q}{\Vert }_{0,\tau }& \le C{h}_{\tau }\Vert \mathbf{q}{\Vert }_{1,\tau },\\ \Vert \mathbf{q}-J_i\mathbf{q}{\Vert }_{0,e}& \le C{h}_{\tau }^{1/2}\Vert \mathbf{q}{\Vert }_{1,\tau }.\end{array}$$

$${\eta }_{\tau }^2$$

$$+\sum _{e\in \partial \tau \cap \Gamma }{h}_e^3\Vert {\lambda }_h{\mathbf{n}}_i\cdot {\mathbf{n}}_{ij}-{\mathbf{z}}_h\cdot {\mathbf{n}}_i{\Vert }_{0,e}^2+\sum _{e\in \partial \tau \cap ({\mathcal{T}}^{\Gamma ,h}\cup {\mathcal{F}}_p)}{h}_e\Vert [[u_h]]{\Vert }_{0,e}^2.$$

Full text

A posteriori error estimates for the mortar

staggered DG method

Lina Zhao111Department of Mathematics, The Chinese University of Hong Kong, Hong Kong Special Administrative Region. (lzhao@math.cuhk.edu.hk)   Eric Chung222Department of Mathematics, The Chinese University of Hong Kong, Hong Kong Special Administrative Region. (tschung@math.cuhk.edu.hk)

Abstract: Two residual-type error estimators for the mortar staggered discontinuous Galerkin discretizations of second order elliptic equations are developed. Both error estimators are proved to be reliable and efficient. Key to the derivation of the error estimator in potential $L^{2}$ error is the duality argument. On the other hand, an auxiliary function is defined, making it capable of decomposing the energy error into conforming part and nonconforming part, which can be combined with the well-known Scott-Zhang local quasi-interpolation operator and the mortar discrete formulation yields an error estimator in energy error. Importantly, our analysis for both error estimators does not require any saturation assumptions which are often needed in the literature. Several numerical experiments are presented to confirm our proposed theories.

Keywords: Staggered grids, Discontinuous Galerkin method, Nonmatching grids, A posteriori error estimates, Adaptive mesh refinement

1 Introduction

The mortar element method is a domain decomposition method with non-overlapping subdomains [5, 6]. One distinctive feature of mortar finite element method is that the meshes on adjacent subdomains are not required to be matching with each other, which makes the method well suited for problems with complicated geometries. Local features of the solution such as corner singularities or large gradients can be resolved by finer grids in the local region. Furthermore, large scale features such as geological faults and layers in subsurface flow can be modeled with nonmatching grids. Staggered discontinuous Galerkin (SDG) methods pioneered by Chung and Engquist [13, 14] earn many desirable properties such as mass conservation, superconvergence and the flexibility to deal with general quadrilateral and polygonal meshes, and have been applied to numerous partial differential equations arising from the practical applications (cf. [15, 28, 17, 20, 16, 12, 38, 40]). To further advance the applications of SDG method, a mortar formulation is developed for SDG method [24], where different triangulations in different regions of the computational domain are exploited. In the framework proposed therein, SDG discretization is employed in each subdomain and the continuity of the solution across the subdomain interfaces is imposed through the introduction of the Lagrange multipliers. The analysis developed therein shows that optimal convergence rates in both $L^{2}$ and discrete energy norms are achieved. In addition, the numerical results there illustrate that if the exact solution earns local singularities, one can only obtain optimal convergence rates in regularity, not in rate. To efficiently capture the singularities and achieve optimal approximation with minimum degrees of freedom, adaptive finite element method based on a posteriori error estimators can be utilized. Due to the nonmatching meshes across the subdomain interfaces, mortar finite element methods are favored for adaptive mesh refinement. Indeed, nonmatching grids can be used on the different subdomains of a partition, this can highly reduce the number of degrees of freedom since no further nodes must be added to avoid the nonconforming meshes on the subdomain interfaces.

A posteriori error estimators have been actively studied for mixed finite element methods and discontinuous Galerkin methods on conforming grids [7, 31, 8, 22, 1, 23, 27, 25, 33, 26, 11, 21, 9, 19, 39] since the pioneering work of Babuška and Rheinboldt [2, 3]. However, a posteriori error analysis for the discretization problems on nonmatching grids is still a largely undeveloped area. Wohlmuth introduces residual type and hierarchical type a posteriori error estimators in [36, 37] for mortar finite element methods. Wheeler and Yotov [35] propose two types of a posteriori error estimators for the mortar mixed finite element method. All these estimators are constructed with some saturation assumptions. To avoid this assumption, Bernardi and Hecht present some residual type error estimators without the presence of saturation assumptions [4]. But the mesh nodes are required to be coinciding on the interface. To extend the current framework of a posteriori error analysis developed on matching grids to nonmatching grids in a more general sense, it is important to exclude the mesh restrictions on the mortar and non-mortar sides, and avoid the saturation assumptions. Recently, some residual type a posteriori error estimators are developed based on the posteriori version of the well known Strang Lemma [34], where the aforementioned restrictions are remitted.

The purpose of this paper is to design and analyze two residual type error estimators for mortar SDG method without any saturation assumptions. We first derive a reliable and efficient error estimator for mortar SDG method in potential $L^{2}$ error, where the key ingredient is the duality argument. In contrast to the general error estimators developed for conforming grids, the jump of solution across the subdomain interfaces and a mortar flux difference term are also involved. Then we propose an error estimator in energy error, where some difficulties arise due to the following aspects: First, the interface grids are nonmatching across the adjacent subdomains, defining related conforming finite element spaces on the non-matching meshes as the methodology exploited in [10] is impossible; Second, mortar SDG method is a mixed type method, applying the Galerkin orthogonality directly like those proposed in [34] is infeasible. To overcome the aforementioned difficulties, we employ the well-known Scott-Zhang local quasi-interpolation operator defined in [30], which is conforming in each subdomain avoiding constructing a conforming operator over the whole domain that is usually cumbersome. Then we define an auxiliary function $s\in H^{1}_{0}(\Omega)$, which enables us to decompose the energy error into conforming part and nonconforming part. Combing the above primary ingredients, the error estimator in energy error can be derived. Again, in addition to element residual terms, the jump of solution across the subdomain interfaces and the mortar flux difference terms are also involved, where the presence of the additional terms is due to the mortar matching condition. We emphasize that our analysis for both error estimators does not need any mesh restrictions on the mortar and non-mortar sides, and saturation assumptions are also avoided.

The rest of the paper is organized as follows. In the next section, we briefly introduce the mortar formulation of SDG method. In Section 3, two residual type error estimators are proposed, and the reliability of the proposed error estimators are proved. Then, the efficiency of the proposed error estimators are established in Section 4. Several numerical experiments are carried out in Section 5, where the performances of the two error estimators are displayed. Numerical results demonstrate that singularities can be well captured and optimal convergence rates can be achieved under the adaptive mesh refinement. Finally, some conclusions are given at the end of this paper.

2 Mortar formulation of SDG method

In this section, we briefly describe the mortar formulation of SDG method by following the framework developed in [24]. The primary ingredient is to impose the continuity of the solution across subdomain interfaces by a mortar matching condition. To begin, we consider the following second order elliptic problem in two dimensions:

[TABLE]

where $\Omega$ is the computational domain and $f(x)$ is a given source function. We divide the domain $\Omega$ into a set of $N$ non-overlapping subdomains, $\bar{\Omega}=\cup_{i=1}^{N}\bar{\Omega}_{i}$. We assume, for simplicity, that $\{\Omega_{i}\}_{i=1}^{N}$ is a geometrically conforming partition of $\Omega$. We further assume that $\rho$ is a piecewise constant function, which equals $\rho_{i}$ in $\Omega_{i}$. Every subdomain $\Omega_{i}$ is equipped with a quasi-uniform triangulation $\mathcal{T}_{h_{i}}$ with mesh size $h_{i}>0$. The triangulations $\{\mathcal{T}_{h_{i}}\}_{i=1}^{N}$ can be non-matching across the subdomain interface $\Gamma=\cup_{i,j=1}^{N}\Gamma_{ij}$, where $\Gamma_{ij}(=\partial\Omega_{i}\cap\partial\Omega_{j})$ is the interface shared by the two subdomains $\Omega_{i}$ and $\Omega_{j}$. In addition, $S=\{(i,j):\Omega_{i}\,\mbox{and}\,\Omega_{j}\,\mbox{have non-empty intersection}\}$. In addition, we define $\Gamma_{i}=\partial\Omega_{i}\cap\Gamma=\partial\Omega_{i}\backslash\partial\Omega$.

Let $D\subset\mathbb{R}^{d},d=1,2$, we adopt the standard notations for the Sobolev spaces $H^{s}(D)$ and their associated norms $\|\cdot\|_{s,D}$, and semi-norms $|\cdot|_{s,D}$ for $s\geq 0$. The space $H^{0}(D)$ coincides with $L^{2}(D)$, for which the norm and inner products are denoted as $\|\cdot\|_{D}$ and $(\cdot,\cdot)_{D}$, respectively. If $D=\Omega$, the subscript $\Omega$ will be dropped unless otherwise mentioned. In the sequel, we use $C$ to denote a generic positive constant independent of the meshsize which can have different values at different occurrences.

Next, we define some spaces which will be utilized later

[TABLE]

and

[TABLE]

We rewrite (2.1) into a first order system by introducing an additional unknown $\bm{z}$

[TABLE]

which can be recast into the equivalent subdomain problem

[TABLE]

with an additional condition that $u_{i}$ are continuous across the subdomain interface. Let $\bm{n}_{ij}$ be the fixed unit normal direction on $\Gamma_{ij}$ common to the two subdomains $\Omega_{i}$ and $\Omega_{j}$. We define $\lambda$ on $\Gamma$ as

[TABLE]

Multiplying the equations in (2.2) by the corresponding test functions and integration by parts, we obtain the weak formulation: find $(u,\bm{z},\lambda)\in V\times Q\times M$ such that

[TABLE]

where $\lambda\,\bm{n}_{ij}\cdot\bm{n}_{i}=\bm{z}\mid_{\Omega_{i}}\cdot\bm{n}_{i}$.

We can rewrite (2.3) as: find $(\bm{z},u,\lambda)\in Q\times V\times M$ such that

[TABLE]

where

[TABLE]

We then present the construction of the SDG spaces for each $\Omega_{i}$, and the construction follows the framework given in [24]. To this end, we first introduce some notations that will be employed later. We let $\mathcal{F}_{u,i}$ be the set of edges in the initial triangulation $\mathcal{T}_{h_{i}}$ excluding the edges on the interface and $\mathcal{F}_{u,i}^{0}\subset\mathcal{F}_{u,i}$ be the set of interior edges. For each triangle $\tau\in\mathcal{T}_{h_{i}}$, we divide it into three subtriangles by connecting an interior point to the three vertices. We note that the interior point can be chosen as the centroid of the triangle to get a good regularity of the subdivided triangulation.

We denote by $\mathcal{T}_{i}$ the resulting finer triangulation and by $\mathcal{F}_{p,i}$ the set of edges generated by the subdivision process. In addition, we let $\mathcal{T}=\cup_{i=1}^{N}\mathcal{T}_{i}$, $\mathcal{T}_{h}=\cup_{i=1}^{N}\mathcal{T}_{h_{i}}$, $\mathcal{F}_{p}=\cup_{i=1}^{N}\mathcal{F}_{p,i}$, $\mathcal{F}_{u}=\cup_{i=1}^{N}\mathcal{F}_{u,i}$, $\mathcal{F}_{u}^{0}=\cup_{i=1}^{N}\mathcal{F}_{u,i}^{0}$, $\mathcal{F}^{0}=\mathcal{F}_{u}^{0}\cup\mathcal{F}_{p}$, and $\mathcal{F}=\mathcal{F}_{u}\cup\mathcal{F}_{p}$. We use $h_{\tau}$ to denote the diameter of $\tau\in\mathcal{T}$, $h_{e}$ to denote the length of edge $e$, and $h=\max_{\tau\in\mathcal{T}}h_{\tau}$.

For each edge $e$, we define a unit normal vector $\bm{n}_{e}$ as follows: If $e\in\mathcal{F}\setminus\mathcal{F}^{0}$, then $\bm{n}_{e}$ is the unit normal vector of $e$ pointing towards the outside of $\Omega$. If $e\in\mathcal{F}^{0}$, an interior edge, we then fix $\bm{n}_{e}$ as one of the two possible unit normal vectors on $e$. When there is no ambiguity, we use $\bm{n}$ instead of $\bm{n}_{e}$ to simplify the notation.

Let $k\geq 0$ be the order of polynomial used for the approximation and $P^{k}(\tau)$ be the set of polynomials with degree less than or equal to $k$ defined on $\tau$. We define the following spaces

[TABLE]

and

[TABLE]

where the jumps $\llbracket\bm{q}\cdot\bm{n}\rrbracket\mid_{e}$ and $\llbracket v\rrbracket\mid_{e}$ are defined in the standard way

[TABLE]

In the above, $\tau_{1}$ and $\tau_{2}$ are the two triangles with the common edge $e$. In the above definition, we assume $\bm{n}$ is pointing from $\tau_{1}$ to $\tau_{2}$. In addition, we define

[TABLE]

On the whole computational domain, we define $V_{h}:=\prod_{i=1}^{N}V_{h_{i}}^{0}$ and $Q_{h}:=\prod_{i=1}^{N}Q_{h_{i}}$.

We recall that $\Gamma_{ij}$ is the interface between $\Omega_{i}$ and $\Omega_{j}$, see Fig. 2. On $\Gamma_{ij}$, we introduce two different meshes called $T_{ij,i}$ and $T_{ij,j}$, which are respectively defined as the restrictions of $\mathcal{T}_{h_{i}}$ and $\mathcal{T}_{h_{j}}$ on $\Gamma_{ij}$. Among these two meshes, we select one as non-mortar mesh and the other as mortar mesh. For the non-mortar mesh, say $T_{ij,i}$, we introduce the space of Lagrange multipliers $M_{ij}=V_{h_{i}}^{0}\mid_{\Gamma_{ij}}$, which consists of piecewise polynomials of degree up to $k$ defined on $\Gamma_{ij}$ with respect to the mesh $T_{ij,i}$. Also, we denote the union of all the non-mortar mesh as $\mathcal{T}^{\Gamma,h}=\cup_{i,j=1}^{N}T_{ij,i}$. In addition, we define $M_{h}:=\Pi_{(i,j)\in S}M_{ij}$. The space $M_{h}$ is used to enforce continuity of functions in $V_{h}$. Specially, we define the following mortar SDG space for the approximation of $u$

[TABLE]

Furthermore, we use $S(\nu)$ to denote the triangle in the initial triangulation $\mathcal{T}_{h_{i}}$ with $\nu$ denoting the interior point chosen in the above subdivision process. Thus, $S(\nu)$ is the union of the three triangles in $\mathcal{T}_{i}$ having the interior point $\nu$ as a common vertex. For an edge $e\in\mathcal{F}_{u,i}^{0}$, we let $D(e)$ be the union of the two triangles in $\mathcal{T}_{i}$ sharing the edge $e$, and for an edge $e\in\mathcal{F}_{u,i}\cap\partial\Omega_{i}$, we let $D(e)$ be the triangle in $\mathcal{T}_{i}$ having the edge $e$, see Fig. 1 for an illustration. In addition for $e\in\mathcal{T}^{\Gamma,h}$, we use $D(e)$ to denote the union of the simplicial submeshes on both sides sharing the edge $e$ or part of $e$.

To derive the discrete version for (2.2), we introduce $\lambda_{h}\in M_{h}$ to approximate $\lambda$. We note that $\lambda_{h}\mid_{\Gamma_{ij}}$ is considered as an approximation of the flux $\bm{z}\cdot\bm{n}_{ij}$ on $\Gamma_{ij}$. Following [24], we define

[TABLE]

We also define the following bilinear forms

[TABLE]

where the gradient and divergence operators are elementwise operators. Integration by parts reveals that the above bilinear forms are adjoint to each other, namely, $b_{i}(\bm{q},v)=b_{i}^{*}(v,\bm{q})$  $\forall(v,\bm{q})\in V_{h_{i}}^{0}\times Q_{h_{i}}$.

With the aforementioned ingredients, the mortar SDG discretization for (2.2) reads: find $(u_{h},\bm{z}_{h},\lambda_{h})\in V_{h}\times Q_{h}\times M_{h}$ such that

[TABLE]

which can be rewritten as: find $(u_{h},\bm{z}_{h},\lambda_{h})\in V_{h}\times Q_{h}\times M_{h}$ such that

[TABLE]

where

[TABLE]

We infer from integration by parts

[TABLE]

We recall some a priori error estimates from [24] which are needed later to illustrate the efficiency of the proposed error estimators.

Lemma 2.1**.**

Let $u|_{\Omega_{i}}\in H^{1+\sigma_{i}}(\Omega_{i})$ with $\sigma_{i}>1/2$. Let $(\bm{z}_{h},u_{h},\lambda_{h})\in Q_{h}\times V_{h}\times M_{h}$ be the solution of (2.4). Then the following estimates hold

[TABLE]

3 Reliability

In this section, two residual-type error estimators are proposed. First, we develop an error estimator in potential $L^{2}$ error, which mainly relies on the duality argument. Next, we propose an energy error estimator based on an auxiliary function and the well-known Scott-Zhang local quasi-interpolation operator.

3.1 Potential error estimator

To begin, we recall the following trace inequality

[TABLE]

and

[TABLE]

We then define two interpolation operators $I_{i}:H^{1}(\Omega_{i})\rightarrow V_{h_{i}}$ and $J_{i}:H^{\epsilon}(\Omega_{i})^{2}\rightarrow Q_{h_{i}},\epsilon>1/2$ by

[TABLE]

and

[TABLE]

In addition, we let $\pi_{e}$ be the $L^{2}$ projection operator onto $M_{h}$. Then, the following inequalities hold true for smooth enough functions $\bm{q},v$ and $\mu$ (cf. [18, 24])

[TABLE]

On each element $\tau\in\mathcal{T}$, we define the local error estimator as

[TABLE]

Then the global error estimator in potential $L^{2}$ error can be defined as

[TABLE]

The main result of this section can be stated in the next theorem.

Theorem 3.1**.**

Let $(\bm{z},u)$ be the weak solution of (2.3) and $(\bm{z}_{h},u_{h},\lambda_{h})\in Q_{h}\times V_{h}\times M_{h}$ be the solution of (2.4). Let $\eta_{1}$ be defined in (3.4), then there exists a positive constant $C$ such that

[TABLE]

Proof.

Assume the auxiliary problem

[TABLE]

satisfies the elliptic regularity estimate

[TABLE]

Let $\bm{U}=-\rho\nabla w$ and $\bm{U}_{i}\cdot\bm{n}_{i}=\mu_{i}$ on $\partial\Omega_{i}\setminus\partial\Omega$, then (3.5) can be recast into the following first order system

[TABLE]

Notice that the above problem is equivalent to the following subdomain problems

[TABLE]

with the additional condition that $w_{i}$ are continuous across the subdomain interfaces.

Multiplying the first equation of (3.7) by $\bm{z}-\bm{z}_{h}$, the second equation by $u-u_{h}$ and integrating over $\Omega$ to get

[TABLE]

Integration by parts, employing the definition of $A_{h}$ (cf. (2.5)) and using $\bm{U}_{i}\cdot\bm{n}_{i}=\mu_{i}$ on $\partial\Omega_{i}\backslash\partial\Omega$ yield

[TABLE]

Employing (2.5) and (2.6), we deduce that

[TABLE]

The Cauchy-Schwarz inequality and the approximation properties (3.3) imply

[TABLE]

This together with (3.2) and the elliptic regularity estimate (3.6) completes the proof.

∎

3.2 Energy error estimator

This section is devoted to the construction of the error estimator in energy error, the primary ingredient is to define an auxiliary function which enables us to decompose the error into conforming part and nonconforming part. To this end, we first introduce the following lemma, which provides the upper bound for the nonconforming error.

Following Lemma 3.6 of [34] and Theorem 2.2 of [22], we have

Lemma 3.1**.**

There exists a positive constant $C$ independent of the mesh size such that

[TABLE]

Let $\Pi_{h}:H^{1}(\Omega_{i})\rightarrow X_{h}(\Omega_{i})$ be the well-known Scott-Zhang local quasi-interpolation operator, where $X_{h}(\Omega_{i})$ is $P^{1}$ conforming element space in each subdomain $\Omega_{i}$. In addition, $\Pi_{h}$ satisfies the following interpolation error estimates (cf. [30]).

Lemma 3.2**.**

We have the following interpolation error estimates for $v\in H^{1}_{0}(\Omega)$

[TABLE]

where $\mathcal{E}_{h}(\bar{\Omega}_{i})$ denotes all the edges of $\mathcal{T}_{i}$ in $\bar{\Omega}_{i}$, and $\omega_{\tau},\omega_{e}$ denotes the union of all the elements in $\mathcal{T}_{i}$ sharing at least a node with $\tau$ and $e$, respectively.

We define the local error estimator on each element $\tau\in\mathcal{T}$ as

[TABLE]

Then, the global error estimator in energy error can be defined as

[TABLE]

Theorem 3.2**.**

There exists a positive constant $C$ such that the following estimate holds

[TABLE]

Proof.

We first define a function $s\in H^{1}_{0}(\Omega)$ such that

[TABLE]

where the existence and uniqueness of $s$ follow from Riesz representation theorem.

By taking $v=u-s$ in (3.10), we can get

[TABLE]

We can first bound the second term by employing Lemma 3.1 yielding

[TABLE]

On the other hand, we have from the definition

[TABLE]

Integration by parts implies

[TABLE]

The penultimate term of (3.13) can be estimated by integration by parts

[TABLE]

The last term of (3.13) can be recast into the following form by exploiting integration by parts and the second equation of (2.4)

[TABLE]

Finally, we can obtain by combing the above equations

[TABLE]

which, coupling with (3.11), (3.12) and Lemma 3.2 yields the desired estimate.

∎

4 Efficiency

This section is devoted to establishing the lower bounds on the errors. To this end, we set the element bubble function for each element $\tau$ as $\psi_{\tau}$ and the edge bubble function for each edge $e$ as $\psi_{e}$. The properties of the bubble functions are given in the next lemma.

Lemma 4.1**.**

The following inequalities hold for all functions $v\in P^{k}(\tau)$.

[TABLE]

Moreover, there exists an extension operator $P_{e}$ that extends any function defined on $e\in\mathcal{F}$ to the element $\tau$ and satisfies

[TABLE]

Then, the lower bounds on the errors can be stated in the next theorem.

Theorem 4.1**.**

Let $\mathcal{T}_{h}$ be shape regular and let $f_{h}$ be the piecewise linear polynomial approximation of $f$. Then, there exists a positive constant $C$ independent of $h$ such that

[TABLE]

and

[TABLE]

In addition, the following local bounds hold for any $\tau\in\mathcal{T}_{h}$,$e\in\partial\tau$ and $\hat{e}\in\mathcal{F}_{p}\cup\mathcal{T}^{\Gamma,h}$

[TABLE]

and

[TABLE]

Proof.

Let $R_{\tau}(f)=f+\nabla\cdot\bm{z}_{h}$. Green’s theorem and the Cauchy-Schwarz inequality imply

[TABLE]

where in the last inequality, we use (4.1) and (4.2).

Combining the above inequality with (4.1), we can achieve

[TABLE]

which yields

[TABLE]

Next, fix an edge $e\in\mathcal{F}_{u}^{0}$, for any $w=\psi_{e}P_{e}\llbracket\bm{z}_{h}\cdot\bm{n}\rrbracket\in H^{1}_{0}(D(e))$, we have from integration by parts

[TABLE]

which, coupling with (4.3), inverse inequality and (4.5) yields

[TABLE]

Thus

[TABLE]

Triangle inequality implies

[TABLE]

Finally, the triangle inequality and trace inequality (3.1) yield

[TABLE]

and

[TABLE]

The preceding arguments complete the assertion. ∎

Remark 4.1**.**

It follows from Lemma 2.1 that the orders of convergence for all the terms present in the right hand side of (4.4) are comparable to $\|u-u_{h}\|_{0}$. Thus, this bound, combined with Theorem 3.1, implies that $\eta_{1}$ is an efficient and reliable estimator for the potential $L^{2}$ error. Similarly, $\eta_{2}$ is also an efficient and reliable estimator for the energy error.**

5 Numerical experiments

In this section we present several numerical experiments to demonstrate the performance of the proposed error estimators. The adaptive mesh pattern and convergence history are reported for each example. In all of our simulations, we use piecewise linear elements, i.e., $k=1$. Since the new triangulation $\mathcal{T}$ is only formed to define the method and it is not a refinement. Therefore, in our refinement algorithm, we will carry out the refinement on $\mathcal{T}_{h}$ by an estimator defined on each $\rho\in\mathcal{T}_{h}$. We define the error estimator as

[TABLE]

Moreover, for any subset $\mathcal{M}\subset\mathcal{T}_{h}$, we define

[TABLE]

Similar definitions can be applied for $\bar{\eta}_{\tau}$ defined in (3.8).

Our adaptive refinement can be implemented by the following iteration:

Start with an initial mesh $\mathcal{T}_{h}^{0}$. 2. 2.

Solve the discrete problem (2.4) for $(\bm{u}_{h}^{\ell},\bm{z}_{h}^{\ell},p_{h}^{\ell})$ with respect to $\mathcal{T}_{h}^{\ell}$. 3. 3.

Compute $\xi_{\rho},\forall\rho\in\mathcal{T}_{h}^{\ell}$. 4. 4.

Mark the minimal set $\mathcal{M}\subseteq\mathcal{T}_{h}$ satisfying $\theta\sum_{\rho\in\mathcal{T}_{h}}\xi_{\rho}^{2}\leq\sum_{\rho\in\mathcal{M}}\xi_{\rho}^{2}$ for some fixed parameter $\theta\in(0,1)$. 5. 5.

Refine marked triangles and compute $\mathcal{T}^{\ell+1}_{h}$ by red and green refinement for adaptive mesh. Update $\ell$ and go to step 2.

Example 5.1**.**

In this example, $\Omega=(0,1)^{2}$ and $\rho=1$, we consider the exact solution given by

[TABLE]

The global domain $\Omega$ is decomposed into four square subdomains and the initial grid in each subdomain is $2\times 2$.**

The contour plot of the exact solution and the adaptive mesh pattern arising from the energy error estimator are reported in Fig. 3. The adaptive mesh pattern for the error estimator in potential $L^{2}$ error is similar and is omitted for simplicity. We note that the grids are appropriately refined along the boundary layers.

The convergence history for $\|u-u_{h}\|_{0}$ and $\eta_{1}$ as well as $\|\nabla(u-u_{h})\|_{0}$ and $\eta_{2}$ are displayed in Fig. 4. We observe that the adaptive solution needs much fewer elements to provide the same accuracy.

Example 5.2**.**

In this example, we consider a non-smooth solution in $H^{3/2-\delta}(\Omega)$ with $\delta>0$ defined by $u(r,\theta)=r^{1/2}\cos(2\theta)$ with polar coordinates $(r,\theta)$ centered at $(0.5,0.5)$. In addition, we let $\rho=1$. We assume that the computational domain is decomposed into four square subdomains.**

The initial mesh and adaptive mesh pattern using the energy error estimator are shown in Fig. 6. We observe that the singularly can be well captured by the proposed error estimators.

The convergence history for $\|u-u_{h}\|_{0}$ and $\eta_{1}$ as well as $\|\nabla(u-u_{h})\|_{0}$ and $\eta_{2}$ under uniform refinement and adaptive refinement, respectively, are shown in Fig. 7. It is clear that the order of convergence for $\|u-u_{h}\|_{0}$ and $\eta_{1}$ under uniform refinement is approximately 1.5, while the order of convergence for $\|u-u_{h}\|_{0}$ and $\eta_{1}$ under adaptive refinement is approximately 2. On the other hand, the order of convergence for $\|\nabla(u-u_{h})\|_{0}$ and $\eta_{2}$ under uniform refinement is approximately 0.5, while the order of convergence for $\|\nabla(u-u_{h})\|_{0}$ and $\eta_{2}$ under adaptive refinement is approximately 1. This demonstrates that under uniform refinement we can achieve the reduced convergence rate reflecting singularity, and optimal convergence rates can be recovered by employing adaptive mesh refinement. This example highlights that adaptive mesh refinement outperforms uniform mesh refinement and can lead to optimal convergence rate even with low solution regularity.

Example 5.3**.**

Our third example is an interface problem which exhibits an interface singularity, cf. [29, 32]. Consider $\Omega=(-1,1)\times(-1,1)$, divided into four subdomains $\Omega_{i}$ along the Cartesian axes (the subregion $\{x>0,y>0\}\cap\Omega$ is denoted as $\Omega_{1}$ and the subsequent numbering is done counterclockwise). The exact solution is given by

[TABLE]

in each $\Omega_{i}$. The solution is continuous across the interfaces and the normal component of its flux $\bm{z}$ is continuous; it exhibits a singularity at the origin and it only belongs to $H^{1+\alpha}(\Omega)$. We take the piecewise constant coefficient as $\rho_{1}=\rho_{3}=5,\rho_{2}=\rho_{4}=1$ and $\alpha=0.53544095$. The values of $K_{i},S_{i}$ can be found in, e.g., [32]. **

The initial mesh is the same as Fig. 7 and the convergence history and adaptive mesh pattern are reported in Fig. 8 and Fig. 9, respectively. Again, reduced convergence rate can be achieved for uniform refinement due to low solution regularity, while optimal convergence rates can be recovered by employing adaptive mesh refinement. In addition, the singularity is well captured. This example once again illustrates that the proposed error estimators can guide adaptive mesh refinement.

6 Conclusion

In this paper, we have proposed two residual-type error estimators in potential $L^{2}$ error and energy error, respectively. The proposed error estimators are proved to be reliable and efficient. The key idea is to exploit the duality argument for potential $L^{2}$ error. To derive an error estimator in energy error, we decompose the energy error into conforming part and nonconforming part via the introduction of an auxiliary function. The numerical results demonstrate that the singularities can be well captured by the proposed error estimators, in addition, the superiority of adaptive mesh refinement over uniform mesh refinement is clearly visible in the improved convergence rate for solutions of limited regularity.

Acknowledgements

The research of Eric Chung is partially supported by the Hong Kong RGC General Research Fund (Project numbers 14304217 and 14302018), CUHK Faculty of Science Direct Grant 2018-19 and NSFC/RGC Joint Research Scheme (Project number HKUST620/15).