Fully computable a posteriori error bounds for hybridizable discontinuous Galerkin finite element approximations

TL;DR

This paper develops fully computable a posteriori error bounds for hybridizable discontinuous Galerkin methods applied to second-order elliptic problems, enabling reliable error estimation in finite element approximations.

Contribution

It introduces constant-free, fully computable a posteriori error estimators for primal and mixed HDG methods, including local lower bounds and applicability to various formulations.

Findings

Error bounds are fully computable and constant-free.

Estimators provide local lower bounds up to a constant.

Numerical examples confirm theoretical accuracy.

Abstract

We derive a posteriori error estimates for the hybridizable discontinuous Galerkin (HDG) methods, including both the primal and mixed formulations, for the approximation of a linear second-order elliptic problem on conforming simplicial meshes in two and three dimensions. We obtain fully computable, constant free, a posteriori error bounds on the broken energy seminorm and the HDG energy (semi)norm of the error. The estimators are also shown to provide local lower bounds for the HDG energy (semi)norm of the error up to a constant and a higher-order data oscillation term. For the primal HDG methods and mixed HDG methods with an appropriate choice of stabilization parameter, the estimators are also shown to provide a lower bound for the broken energy seminorm of the error up to a constant and a higher-order data oscillation term. Numerical examples are given illustrating the theoretical results.

Figures

Equations

$$\displaystyle\left.\begin{tabular}[]{r l}${\boldsymbol{\sigma}}-a\,{\nabla}u=$&\; \!\!\!\!$0$\\ $-{\nabla\cdot}{\boldsymbol{\sigma}}=$&\; \!\!\!\!$f\in L^{2}(\Omega)$\\ \end{tabular}\right\}\text{ in }\Omega$$

$$\tilde{F}=\{K^{\prime }\in {\mathcal{T}}_h:F\in \mathcal{E}(K^{\prime })\},$$

$$\tilde{K}=\{K^{\prime }\in {\mathcal{T}}_h:\mathcal{E}(K)\cap \mathcal{E}(K^{\prime })\text{ is nonempty}\}.$$

$${\Sigma }_{h,k}:=$$

$$V_{h,k}:=$$

$$M_{h,k}:=$$

$$M_{h,k}^0:=\{\hat{v}\in M_{h,k}:\hat{v}\text{raise}-0.86108pt{|}_F=0\forall F\in {\mathcal{E}}_h^{\partial }\}.$$

$${\mathbb{V}}_{h,k,\delta }:=$$

$${\mathbb{X}}_{h,k}:=$$

$$M_{h,m}^{\mathrm{dc}}:={\Pi }_{K\in {\mathcal{T}}_h}{\mathbb{P}}_m(\partial K),$$

$$\{q\}=\frac{1}{2}(q\text{raise}-0.86108pt{|}_{K^{+}}+q\text{raise}-0.86108pt{|}_{K^{-}}),[[q]]=q\text{raise}-0.86108pt{|}_{K^{+}}{\mathbf{n}}^{+}+q\text{raise}-0.86108pt{|}_{K^{-}}{\mathbf{n}}^{-}\text{ on }F\in {\mathcal{E}}_h^o.$$

$$\{\mathbf{\varphi }\}=\frac{1}{2}(\mathbf{\varphi }\text{raise}-0.86108pt{|}_{K^{+}}+\mathbf{\varphi }\text{raise}-0.86108pt{|}_{K^{-}}),[[\mathbf{\varphi }]]=\mathbf{\varphi }\text{raise}-0.86108pt{|}_{K^{+}}\cdot {\mathbf{n}}^{+}+\mathbf{\varphi }\text{raise}-0.86108pt{|}_{K^{-}}\cdot {\mathbf{n}}^{-}\text{ on }F\in {\mathcal{E}}_h^o.$$

$$\{q\}=q,[[q]]=q\mathbf{n}\text{ on }F\in {\mathcal{E}}_h^{\partial }.$$

$$\displaystyle\mathcal{B}_{h}^{pr}\big{(}(u,\widehat{u}),(v,\widehat{v})\big{)}=$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!-\langle a\,{\nabla}v\cdot{\boldsymbol{n}}\,,\,u-\widehat{u}\rangle_{\partial{K}}+\langle\alpha_{h}(P_{M}u-\widehat{u})\,,\,P_{M}v-\widehat{v}\,\rangle_{\partial{K}}\Big{\}},$$

$$\displaystyle\mathcal{L}_{h}^{pr}\big{(}(v,\widehat{v})\big{)}=$$

$$\displaystyle\mathcal{B}_{h}^{pr}\big{(}(u_{h},\widehat{u}_{h}),(v_{h},\widehat{v}_{h})\big{)}=\mathcal{L}_{h}^{pr}\big{(}(v_{h},\widehat{v}_{h})\big{)}\quad\forall(v_{h},\widehat{v}_{h})\in\mathbb{V}_{h,k,\delta}.$$

$${\alpha }_h\text{raise}-0.86108pt{|}_F=\frac{a\text{raise}-0.86108pt{|}_K\gamma }{|F|}\sum _{F^{\prime }\in \mathcal{E}(K)}\frac{|F^{\prime }{|}^2}{|K|}\text{ for all }F\in \mathcal{E}(K),\text{ for all }K\in {\mathcal{T}}_h,$$

$$\gamma >\frac{k(k+d-1)}{d}.$$

$${\left|\left|\left|(v,\hat{v})\right|\right|\right|}_{pr}={\left(\sum _{K\in {\mathcal{T}}_h}(a\nabla v,\nabla v{)}_K\right)}^{1/2}.$$

$${\left|\left|\left|(v,\hat{v})\right|\right|\right|}_{HDG,pr}={\left({\left|\left|\left|(v,\hat{v})\right|\right|\right|}_{pr}^2+\sum _{K\in {\mathcal{T}}_h}⟨{\alpha }_h{P}_M(v-\hat{v}),P_M(v-\hat{v}){⟩}_{\partial K}\right)}^{1/2}.$$

$$\sum _{K\in {\mathcal{T}}_h}⟨{\alpha }_h{P}_M(e_u-{\hat{e}}_u),P_M(e_u-{\hat{e}}_u){⟩}_{\partial K}$$

$${\left|\left|\left|(e_u,{\hat{e}}_u)\right|\right|\right|}_{pr}\le {\left|\left|\left|(e_u,{\hat{e}}_u)\right|\right|\right|}_{HDG,pr},$$

$$c{\left|\left|\left|(e_u,{\hat{e}}_u)\right|\right|\right|}_{HDG,pr}^2\le {\left|\left|\left|(e_u,{\hat{e}}_u)\right|\right|\right|}_{pr}^2+\sum _{K\in {\mathcal{T}}_h}os{c}_{k-1}^2(f,K).$$

$$os{c}_m(f,K)=a^{-1/2}{h}_K\underset{p\in {\mathbb{P}}_m}{\inf }\Vert f-p{\Vert }_K.$$

$${\hat{\mathbf{\sigma }}}_{h,pr}:=a\nabla u_h-{\alpha }_h(P_M{u}_h-{\hat{u}}_h)\mathbf{n}\in [M_{h,k-\delta }^{\mathrm{dc}}{]}^d$$

$$(f,1{)}_K+⟨{\hat{\mathbf{\sigma }}}_{h,pr}\cdot \mathbf{n},1{⟩}_{\partial K}=0\text{ for each element K∈Th},$$

$$[[{\hat{\mathbf{\sigma }}}_{h,pr}]]\text{raise}-0.86108pt{|}_F=0\text{ on each interior facet }F\in {\mathcal{E}}_h^o.$$

$$(\nabla \cdot {\mathbf{\sigma }}_h^{\ast },v{)}_K=$$

$$⟨{\mathbf{\sigma }}_h^{\ast }\cdot \mathbf{n},\hat{v}{⟩}_{\partial K}=$$

$$({\mathbf{\sigma }}_h^{\ast },\mathbf{\tau }{)}_K=$$

$${\Sigma }_{k,sbb}(K):=\{\mathbf{\tau }\in {\mathbb{P}}_k(K{)}^d:\nabla \cdot \mathbf{\tau }=0,\mathbf{\tau }\cdot \mathbf{n}\text{raise}-0.86108pt{|}_F=0\text{ for all }F\in \mathcal{E}(K)\}$$

$$\nabla \cdot {\mathbf{\sigma }}_h^{\ast }=-{\Pi }_{k-1}f\text{ on Th},\text{ and }{\mathbf{\sigma }}_h^{\ast }\cdot \mathbf{n}={\hat{\mathbf{\sigma }}}_{h,pr}\cdot \mathbf{n}\text{ on ∂Th}.$$

$$\displaystyle S_{k}(u_{h})({\boldsymbol{x}}_{m})=\left\{\begin{tabular}[]{l l}$0$&if $m\in\mathcal{N}_{\gamma,k}^{\partial}$,\vspace{0.12cm}\\ $u_{h}({\boldsymbol{x}}_{m})$&if $m\in\mathcal{N}_{K,k}^{o}$,\vspace{0.12cm}\\ $\frac{1}{\#\Omega_{m}}\sum_{K^{\prime}\in\Omega_{m}}u_{h}\raise-0.86108pt\hbox{$|$}_{K^{\prime}}({\boldsymbol{x}}_{m})$&if $m\in\mathcal{N}_{\gamma,k}\backslash\mathcal{N}_{\gamma,k}^{\partial}$,\end{tabular}\right.$$

Full text

We obtain fully computable, constant free, a posteriori error bounds on the broken energy seminorm and the HDG energy (semi)norm of the error. The estimators are also shown to provide local lower bounds for the HDG energy (semi)norm of the error up to a constant and a higher-order data oscillation term. For the primal HDG methods and mixed HDG methods with an appropriate choice of stabilization parameter, the estimators are also shown to provide a lower bound for the broken energy seminorm of the error up to a constant and a higher-order data oscillation term. Numerical examples are given illustrating the theoretical results.

First author gratefully acknowledges the partial support of this work under AFOSR contract FA9550-12-1-0399.

1. Introduction

Recent years have seen the developments of fully computable, guaranteed error bounds for the conforming [40, 37, 6, 7, 30, 42, 43, 14, 51, 18, 15], nonconforming [29, 1, 38, 39, 33, 19, 8], discontinuous Galerkin [2, 38, 39, 21, 9, 10], and mixed finite element methods [3, 38, 5, 27]; see also unified frameworks in [4, 16, 34].

In comparison, there are relatively few works on a posteriori error estimates for the hybridizable discontinuous Galerkin (HDG) methods [22]. The a posteriori estimates for HDG methods that are currently available in the literature [17, 25, 26, 31, 20, 35] are all of residual type, in which reliability is shown up to a generic (unknown) constant. This means that, while the associated estimation may be suitable as local refinement indicators, they cannot provide a quantitative stopping criterion. Moreover, if only a single fixed mesh is used (as is often the case in practice) then the value of an a posteriori bound containing unknown constants is somewhat questionable. Here we present, for the first time, fully computable a posteriori error bounds for HDG methods, for both the primal and mixed formulations, in the setting of a linear second-order elliptic problem on conforming simplicial meshes in two and three space dimensions. The key ingredient of our analysis is the local conservation property of the HDG methods, with which cheap element-wise equilibrated fluxes can be constructed.

The remainder of this article is organized as follows. Section 2 presents the model problem and prepares the notation used throughout the article. In Section 3, we introduce the primal HDG schemes and the corresponding computation error bounds. While in Section 4, we introduce the mixed HDG schemes and the corresponding computation error bounds. Numerical results are then presented in Section 5. The proofs of the main results in Section 3 and Section 4 are presented in Section 6.

2. Preliminaries

2.1. Model Problem

Consider the following model problem:

[TABLE]

subject to $u=0$ on $\partial\Omega$, where $\Omega\in\mathbb{R}^{d}$, $d\in\{2,3\}$, is a simply connected polygonal/polyhedral domain. The datum $a\in L^{\infty}(\Omega)$ is assumed to be strictly positive and, for simplicity, is assumed piecewise constant on subdomains of $\Omega$.

2.2. Notation and finite elements

We consider a family of partitions $\mathcal{T}_{h}=\{K\}$ of the domain $\Omega$ into the union of nonoverlapping, shape-regular, simplicial elements such that the nonempty intersection of a distinct pair of elements is a single common node, single common edge or single common face (in three dimensions). The family of partitions is assumed to be locally quasi-uniform in the sense that the ratio of the diameters of any pair of neighboring elements is uniformly bounded above and below over the whole family. Furthermore, it is assumed that the partitioning is compatible with the data so that $a$ is constant on each element.

The set of all facets (edges in two dimensions and faces in three dimensions) of the elements is denoted by $\mathcal{E}_{h}$, which we partition into subsets $\mathcal{E}_{h}^{\partial}$ and $\mathcal{E}_{h}^{o}$ consisting of facets lying on the boundary $\partial\Omega$, and the remaining interior facets, respectively. Likewise, the corresponding quantities relative to an individual element $K$ are denoted by $\mathcal{E}(K),\mathcal{E}^{\partial}(K),$ and $\mathcal{E}^{o}(K)$, respectively. For each facet $F\in\mathcal{E}_{h}$, the set $\widetilde{F}$ consists of those elements for which $F$ is a facet,

[TABLE]

while, for each element $K\in\mathcal{T}_{h}$, the set $\widetilde{K}$ consists of those elements having a facet in common with $K$,

[TABLE]

Let $h_{D}$ denote the diameter of a domain $D$, $|K|$ denote the measure of an element $K$, $|F|$ denote the measure of a facet $F$, and $|\partial K|=\sum_{F\in\mathcal{E}(K)}|F|$ denote the measure of the boundary of an element $K$.

Let $\Sigma_{h,k},V_{h,k}$, and $M_{h,k}$ denote the finite dimensional spaces

[TABLE]

To simplify notation, we introduce the compound finite-dimensional space

[TABLE]

which is used for the primal HDG scheme, while the compound finite-dimensional space

[TABLE]

is used for the mixed HDG scheme.

The stabilization parameters for the HDG schemes will be taken from the space $M_{h,0}^{\mathrm{dc}}$, where, for any nonnegative integer $m$,

[TABLE]

with \mathbb{P}_{m}({{\partial K}}):=\{\mu\in L^{2}({\partial K}):\;\;\mu\raise-0.86108pt\hbox{|}_{F}\in\mathbb{P}_{m}(F)\;\;\forall F\in\mathcal{E}(K)\}.

We use the standard notation for jumps and averages [12] of functions in $M_{h,m}^{\mathrm{dc}}$ and $[M_{h,m}^{\mathrm{dc}}]^{d}$ on the mesh skeleton $\mathcal{E}_{h}$: Let $F$ be a facet shared by elements $K^{+}$ and $K^{-}$ with unit normal vectors ${\boldsymbol{n}}^{+}$ and ${\boldsymbol{n}}^{-}$ on $F$ pointing exterior to $K^{+}$ and $K^{-}$ respectively, then for $q\in M_{h,m}^{\mathrm{dc}}$

[TABLE]

We use the notation $(\cdot,\cdot)_{\omega}$ to denote the integral inner product over a region ${\omega}\subset\mathbb{R}^{d}$, and $\|\cdot\|_{\omega}$ to denote the corresponding $L^{2}$-norm. We omit the subscript in the case when $\omega$ is the physical domain $\Omega$. Finally, for each element $K\in\mathcal{T}_{h}$, we use the notation $\langle\cdot\,,\,\cdot\rangle_{\partial{K}}$ to denote the integral inner product over the element boundary ${\partial K}$.

3. The primal HDG methods and computable error bounds

3.1. Primal HDG formulation

Let $\alpha_{h}\in M_{h,0}^{\mathrm{dc}}$ be a positive stabilization parameter to be specified later, and define the bilinear form $\mathcal{B}_{h}^{pr}:\mathbb{V}_{h,k,\delta}\times\mathbb{V}_{h,k,\delta}\rightarrow\mathbb{R}$ by

[TABLE]

where $P_{M}$ denotes the $L^{2}$ projection onto the space $M_{h,k-\delta}^{\mathrm{dc}}$. We also define the linear form $\mathcal{L}_{h}^{pr}:\mathbb{V}_{h,k,\delta}\rightarrow\mathbb{R}$ by

[TABLE]

Let $u$ be the solution of (3). An approximation of (u,u\raise-0.86108pt\hbox{|}_{\mathcal{E}_{h}}) is obtained by seeking $(u_{h},\widehat{u}_{h})\in\mathbb{V}_{h,k,\delta}$ such that

[TABLE]

This scheme is known as the hybridized, symmetric, interior penalty discontinuous Galerkin method [41].

3.2. The choice of the stabilization

parameter $\alpha_{h}$.

It is well-known [41] that (13) is well-posed provided the stabilization parameter is “sufficiently large”. The following result quantifies exactly how large $\alpha_{h}$ must be; similar results for interior penalty discontinuous Galerkin methods can be found in [48, 2, 32, 9, 10, 11].

Lemma 1**.**

Suppose the stabilization parameter $\alpha_{h}$ is given by

[TABLE]

where $\gamma$ is a constant satisfying

[TABLE]

Then (13) has a unique solution $(u_{h},\widehat{u}_{h})\in\mathbb{V}_{h,k,\delta}$ for $k\geq 1$ and $\delta\in\{0,1\}$.

The proof of this and other results in this section is postponed to Section 6.

Remark 1**.**

*Observe that the stabilization parameter (14) on a facet $F$ is proportional to $h_{F}^{-1}$, which gives an optimal order a priori convergence rate $\mathcal{O}(h^{k})$ in the energy norm [41, 44]. *

3.3. The broken energy seminorm and the HDG energy norm

For a given function $(v,\widehat{v})\in\mathbb{V}_{h,k,\delta}+(H^{1}(\Omega)\times L^{2}(\mathcal{E}_{h}))$, let the broken energy seminorm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ be denoted by

[TABLE]

Our objective is to derive a fully computable estimator for the error in the HDG finite-element approximation (e_{u},\widehat{e}_{u})=(u-u_{h},u\raise-0.86108pt\hbox{|}_{\mathcal{E}_{h}}-\widehat{u}_{h}), where $u$ is the solution to (3) and $(u_{h},\widehat{u}_{h})$ is the solution to (13). Let the HDG energy norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{HDG,pr}$ be denoted by

[TABLE]

Observe that, since $P_{M}(e_{u}-\widehat{e}_{u})=-P_{M}(u_{h}-\widehat{u}_{h})$, the quantity

[TABLE]

is directly computable in terms of the HDG approximation $(u_{h},\widehat{u}_{h})$. Hence, given a constant free estimator for ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(e_{u},\widehat{e}_{u})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{pr}$, we automatically have a constant free estimator for the HDG energy norm of the error as well. The next result shows that, by analogy with the standard interior penalty methods [2, 10], these norms are equivalent in the following sense.

Lemma 2**.**

Let the stabilization parameter be given by (14) with the global constant $\gamma$ satisfying (15), then the HDG energy norm and the broken energy seminorm of the error (e_{u},\widehat{e}_{u})=(u-u_{h},u\raise-0.86108pt\hbox{|}_{\mathcal{E}_{h}}-\widehat{u}_{h}) are equivalent. That is to say,

[TABLE]

and there exists a positive constant $c$, depending only on the shape-regularity of the mesh, the polynomial degree $k$, and the local permeability ratio between neighboring elements, such that

[TABLE]

where the data oscillation on an element $K\in\mathcal{T}_{h}$ is defined to be

[TABLE]

The proof of this result is postponed to Section 6.

3.4. Local conservation

The numerical flux is defined by

[TABLE]

and satisfies the local conservation property

[TABLE]

along with

[TABLE]

These results are a straightforward consequence of the definition and (13).

3.5. The computable error bounds

We obtain computable error bounds for the discrete energy norm of the error by bounding the conforming and non-conforming errors separately [2, 10, 4]. To this end, two types of post-processing scheme will be needed.

3.5.1. Local (equilibrated) flux post-processing

We define a local flux post-processing [2, 9] as follows: Let ${\boldsymbol{\sigma}}_{h}^{*}\in\Sigma_{h,k}$ be such that, on each element $K$, there holds

[TABLE]

where

[TABLE]

denotes the divergence-free “bubble” space. The unique solvability of (24) can be established using arguments similar to those used to study the closely related Brezzi-Douglas-Marini (BDM) projection [13]. By the local conservation properties (22) and (23), we conclude that ${\boldsymbol{\sigma}}_{h}^{*}\in\Sigma_{h,k}\cap H(\mathrm{div};\Omega)$ satisfies

[TABLE]

The quantity ${\boldsymbol{\sigma}}_{h}^{*}$ is usually referred to as an equilibrated flux [7].

3.5.2. Local potential post-processing

We obtain a globally continuous potential by a simple averaging [2, 10] of the discontinuous potential $u_{h}$ as follows: Let $\mathcal{N}_{K,k}$ index a set of points $\{{\boldsymbol{x}}_{m}\}_{m\in\mathcal{N}_{K,k}}$ on $\overline{K}$ associated with a Lagrange basis for the conforming finite-element space of order $k$ on $\mathcal{T}_{h}$ and let $\mathcal{N}_{K,k}^{o}$ denote the restriction of the set $\mathcal{N}_{K,k}$ to the points that do not lie on the boundary of element $K$, with $\mathcal{N}_{\gamma,k}$ being its complementary set. Let $N_{\gamma,k}^{\partial}$ denote the restriction of the set $\mathcal{N}_{\gamma,k}$ to the points that lie on the closure of the boundary $\partial\Omega$. For $m\in\mathcal{N}_{\gamma,k}$, let $\Omega_{m}$ denote the set of elements in $\mathcal{T}_{h}$ whose closure contains the point ${\boldsymbol{x}}_{m}$.

The post-processed potential $u_{h}^{*}\in V_{h,k}\cap H^{1}(\Omega)$ is obtained through a simple averaging of the degrees of freedom for $u_{h}$: for all elements $K\in\mathcal{T}_{h}$, u_{h}^{*}\raise-0.86108pt\hbox{|}_{K}=S_{k}(u_{h})\raise-0.86108pt\hbox{|}_{K}\in\mathbb{P}_{k}(K), where the nodal values are given by

[TABLE]

and $\#\Omega_{m}$ denotes the number of elements of $\mathcal{T}_{h}$ contained in the patch $\Omega_{m}$.

3.5.3. Computable error bounds

The foregoing developments show that each of the quantities

[TABLE]

can be computed directly from the primal HDG approximation using purely local computations. The next result shows that together these quantities provide a computable, constant-free, upper bound on the broken energy seminorm of the error:

Theorem 1**.**

*Let $\eta_{CF,K}$ and $\eta_{NC,K}$ be defined as in (31), and let the stabilization parameter $\alpha_{h}$ be of the form (14) where $\gamma$ satisfies (15). Then *

[TABLE]

Moreover, there exists a positive constant $c$, depending only on the shape-regularity of the mesh, the polynomial degree $k$, and the local permeability ratio between neighboring elements, such that

[TABLE]

Furthermore,

[TABLE]

and

[TABLE]

where $c$ is as above.

The proof of this result is similar to [2, Theorem 2] and an outline of the main steps is given in Section 6 below. Numerical examples illustrating the bounds in practice are given in Section 5.

4. The mixed HDG methods and computable error bounds

4.1. The mixed HDG formulation

Whereas the primal HDG method gives an approximation for (u,u\raise-0.86108pt\hbox{|}_{\mathcal{E}_{h}}), the mixed HDG method seeks, in addition, to approximate the flux (a{\nabla}u,u,u\raise-0.86108pt\hbox{|}_{\mathcal{E}_{h}}).

Let $\alpha_{h}\in M_{h,0}^{d.c.}$ be a positive stabilization parameter to be specified later, and define the bilinear form $\mathcal{B}_{h}^{mx}:\mathbb{X}_{h,k}\times\mathbb{X}_{h,k}\rightarrow\mathbb{R}$ by

[TABLE]

along with the linear form $\mathcal{L}_{h}^{mx}:\mathbb{X}_{h,k}\rightarrow\mathbb{R}$ by

[TABLE]

An approximation of the true solution (a{\nabla}u,u,u\raise-0.86108pt\hbox{|}_{\mathcal{E}_{h}}) is obtained by seeking $({\boldsymbol{\sigma}}_{h},u_{h},\widehat{u}_{h})\in\mathbb{X}_{h,k}$ such that

[TABLE]

This scheme was originally termed the local discontinuous Galerkin-hybridizable method (LDG-H) [22] but is referred to here as the mixed HDG approximation.

4.2. The choice of the stabilization

parameter $\alpha_{h}$.

The mixed HDG scheme enjoys greater stability properties than the primal HDG scheme as the stabilization parameter $\alpha_{h}$ need only be (partially) positive in order for the scheme to be well-posed as the following result [22, Proposition 3.2] shows.

Lemma 3**.**

If the nonnegative stabilization parameter $\alpha_{h}$ is chosen such that $\alpha_{h}>0$ on at least one facet $F\in\mathcal{E}(K)$ for every element $K$, then there exists a unique solution $({\boldsymbol{\sigma}}_{h},u_{h},\widehat{u}_{h})\in\mathbb{X}_{h,k}$ for $k\geq 0$.

Remark 2** (Stabilization parameter).**

The two most common choices of stabilization parameter used in practice are:

• •

uniform stabilization

[TABLE]

• •

single-facet stabilization

[TABLE]

where $F^{*}_{K}$ is an arbitrarily chosen but fixed facet of $K$.

Each of the above choices of stabilization parameters results in an optimal a priori convergence rate $\mathcal{O}(h^{k+1})$ of the error in the energy norm [23].

4.2.1. The broken energy seminorm and the HDG energy seminorm

For a given function $({\boldsymbol{\tau}},v,\widehat{v})\in\mathbb{X}_{h,k,\delta}+(H(\mathrm{div};\Omega)\times H^{1}(\Omega)\times L^{2}(\mathcal{E}_{h}))$, let the broken energy seminorm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{mx}$ be denoted by

[TABLE]

and let the mixed HDG energy seminorm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{HDG,mx}$ be denoted by

[TABLE]

The error in the mixed HDG finite-element approximation is denoted by ({\boldsymbol{e}}_{\sigma},e_{u},\widehat{e}_{u})=({\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h},u-u_{h},u\raise-0.86108pt\hbox{|}_{\mathcal{E}_{h}}-\widehat{u}_{h}), where $({\boldsymbol{\sigma}},u)$ is the solution to (3) and $({\boldsymbol{\sigma}}_{h},u_{h},\widehat{u}_{h})$ is the solution to (36). Similarly to the primal HDG case, we have $e_{u}-\widehat{e}_{u}=-(u_{h}-\widehat{u}_{h})$ and hence the quantity

[TABLE]

can be evaluated directly given the mixed HDG approximation. Consequently, given a constant free estimator for the broken energy seminorm of the error, we automatically have a constant free estimator for the HDG energy seminorm of the error as well. The next result shows that, by analogy with the primal HDG case, the HDG energy seminorm of the error is equivalent to the broken energy seminorm when the single-facet stabilization (37d) is used.

However, the equivalence fails to hold if the uniform stabilization (37a) is employed. For instance, in the case of lowest order ($k=0$) approximation, the discrete energy norm (plus the data oscillation) cannot control the jump term $\sum_{K\in\mathcal{T}_{h}}h_{K}\langle\alpha_{h}(u_{h}-\widehat{u}_{h})\,,\,(u_{h}-\widehat{u}_{h})\rangle_{\partial{K}},$ as shown by the counterexample constructed in [24, Section 2.4.1]. The situation for general $k\in\mathbb{N}$ remains an open problem at this time.

Lemma 4**.**

Let the stabilization parameter $\alpha_{h}$ given by (37d), then the HDG energy seminorm and the broken energy seminorm of the error ({\boldsymbol{e}}_{\sigma},e_{u},\widehat{e}_{u})=({\boldsymbol{\sigma}}-{\boldsymbol{\sigma}}_{h},u-u_{h},u\raise-0.86108pt\hbox{|}_{\mathcal{E}_{h}}-\widehat{u}_{h}) are equivalent in the sense that

[TABLE]

*and there exists a positive constant $c$, depending only on the shape-regularity of the mesh and the polynomial degree $k$, such that *

[TABLE]

Proof.

The result follows at once from Lemma 7 in Section 6 below. ∎

4.3. Local conservation

Similarly to the primal HDG scheme (13), the mixed HDG scheme (36) is locally conservative but this time in the sense that the numerical flux

[TABLE]

satisfies

[TABLE]

along with

[TABLE]

4.4. The computable error bounds

We are now in a position to present computable error bounds for the discrete energy error. While the basic approach is motivated by the technique used in [3, 5, 27] for the mixed methods, the post-processing technique needed for the mixed HDG case is quite different.

4.4.1. Local (equilibrated) flux post-processing

Let ${\boldsymbol{\sigma}}_{h}^{*}\in\Sigma_{h,k+1}$ satisfy the following conditions, on each element $K$:

[TABLE]

Thanks to the local conservation properties (22) and (23), we conclude that ${\boldsymbol{\sigma}}_{h}^{*}\in\Sigma_{h,k+1}\cap H(\mathrm{div};\Omega)$ satisfies

[TABLE]

4.4.2. Local potential post-processing

A global continuous potential is constructed by averaging a higher order discontinuous approximation to the potential. However, the averaging scheme is more involved than the one used in the primal case: Firstly, we find $u_{h}^{*,\mathrm{dc}}\in V_{h,k+1}$ so that, on each element $K$, there holds

[TABLE]

The continuous potential post-processing $u_{h}^{*}=S_{k+1}(u_{h}^{*,\mathrm{dc}})\in V_{h,k+1}\cap H^{1}(\Omega)$ is then a simple averaging of the degrees of freedom for $u_{h}^{*,\mathrm{dc}}$, where $S_{k+1}(\cdot)$ is defined as in (30).

4.4.3. Computable error bounds

Each of the quantities

[TABLE]

can be computed directly from the mixed HDG approximation using only local computations. These quantities provide computable, constant-free, upper bounds on the the broken energy seminorm of the error $({\boldsymbol{e}}_{\sigma},e_{u},\widehat{e}_{u})$:

Theorem 2**.**

*Let $\eta_{CF,K}$ and $\eta_{NC,K}$ be defined as in (31), with the stabilization parameter $\alpha_{h}$ chosen to be either (37a) or (37d). Then *

[TABLE]

*Moreover, there exists a positive constant $c$, depending only on the shape-regularity of the mesh, the polynomial degree $k$, and the local permeability ratio between neighboring elements, such that *

[TABLE]

where $F_{K}^{*}$ is an arbitrary but fixed facet of $K$. Furthermore,

[TABLE]

with

[TABLE]

The proof of Theorem 2 is presented in Section 6.

Remark 3** (Single-facet stabilization).**

If the stabilization parameter is chosen as in (37d), then we can take $F_{K}^{*}\in\mathcal{E}(K)$ in (48) to be the unique facet on which $\alpha_{h}$ is non-zero which implies that

[TABLE]

Consequently, for the choice (37d), the estimator $\eta^{2}$ also gives a lower bound for the broken energy seminorm.

5. Numerical examples

In order to illustrate the results in Theorem 1–2, we consider Poisson problems in two and three dimensions approximated using the primal and mixed HDG schemes (13) and (36). The implementation is performed using the Python interface of the NGSolve software [46, 47]. Conveniently, NGSolve provides a set of basis functions for the divergence-free bubble space $\Sigma_{k,sbb}$ (25) [45] which makes the implementation of the equilibrated flux reconstructions (24) and (44) relatively straightforward.

We choose the stabilization parameter for the primal HDG schemes (13) to be

[TABLE]

which, thanks to Lemma 1, ensures well-posedness on shape-regular meshes. We adopt a shorthand notation and denote the primal HDG scheme (13) used in conjunction with the approximation space $\mathbb{V}_{h,k,0}$ as pr-Pk, while if with the approximation space $\mathbb{V}_{h,k,1}$ used we use the notation: pr-P${}_{k}^{\mathrm{red}}$. Likewise, we denote the mixed HDG scheme (36) used in conjunction with uniform stabilization (37a) as mx-Pk-U, while if the single-facet stabilization (37d) is used, we write mx-Pk-S. In all cases we use static condensation whereby the local, cell-wise, degrees of freedom are eliminated leaving only those degrees of freedom which are located on the mesh skeleton.

We take polynomial degree $k\in\{1,2,3,4\}$ for the first example, and $k\in\{1,2,3\}$ for the second example.

5.1. Example 1: Two-dimensional L-shaped problem

Here we consider the Laplace problem on a planar L-shaped domain $\Omega_{2D}=(-1,1)\times(0,1)\cup(-1,0)\times(-1,0]$ with Dirichlet boundary conditions. The initial mesh is shown in Fig. 1(A). The true solution is given by $u(r,\theta)=r^{2/3}\sin(2\theta/3)$ in polar coordinates.

The sequence of meshes was constructed by selecting for refinement the smallest number of elements whose combined contribution toward the estimator of the broken energy seminorm of the error exceeds half of the total estimated error. A sample of the meshes for the pr-Pk  scheme with $k=1$ and $k=4$ is shown in Fig. 2.

In Fig. 3, we plot the error for the primal scheme in the broken energy seminorm (16) against the total number of degrees of freedom $\dim\mathbb{V}_{h,k,\delta}$ and the number of global degrees of freedom, $\dim M_{h,k-\delta}^{0}$, remaining after local variables have been eliminated. The error in the broken energy seminorm (38) for the mixed HDG schemes is shown in Fig. 4.

In all cases, the effectivity indices are found to lie in the range $1.0$–$3.0$ as shown in Fig 5.

5.2. Example 2: Three-dimensional L-shaped problem

Here we consider the Laplace problem on a three-dimensional L-shaped domain $\Omega_{3D}=\Omega_{2D}\times(0,0.5)$, where $\Omega_{2D}$ is the two-dimensional L-shaped domain in the previous example. The initial mesh is shown in Fig. 1(B). The true solution is independent of $z$, and reduces to the same two-dimensional solution as Example 1. However, the fact that the true solution is independent of $z$ is not used in the finite element analysis and the meshes are unstructured through the thickness.

A sample of the meshes obtained for adaptive solution using the pr-Pk  scheme with $k=1$ and $k=3$ is shown in Fig. 6. The errors for the primal and mixed HDG schemes are plotted in Fig. 7 and Fig. 8, respectively. As before, the effectivity indices are found to vary in the range $1.0$–$3.0$ as shown in Fig. 9.

6. Proofs

We now turn to the proofs of the results.

As remarked earlier, the jump term in the HDG energy (semi)norm is directly computable, and, as such, we need only concern ourselves with obtaining estimates for the broken energy seminorm of the error. To this end, recall the following Helmholtz decomposition [36, 28]:

Lemma 5**.**

Let $\Omega$ be a simply connected polygon/polyhedron. Then, any ${\boldsymbol{\tau}}\in L^{2}(\Omega)^{d}$, $d\in\{2,3\}$, can be written in the form

[TABLE]

where $\phi\in H^{1}_{0}(\Omega)$ satisfies

[TABLE]

and $\psi\in H^{1}(\Omega)^{2d-3}$ satisfies

[TABLE]

Moreover, the decomposition is orthogonal

[TABLE]

We shall use the decomposition (50) in conjunction with ${\boldsymbol{\tau}}$ defined elementwise by $a{\nabla}u-a{\nabla}u_{h}$ for the primal HDG scheme (13), or by $a{\nabla}u-{\boldsymbol{\sigma}}_{h}$ for the mixed HDG scheme (36). In both cases, Lemma 5 gives an orthogonal decomposition of the broken energy seminorm error into the sum of a conforming part $\|a^{1/2}{\nabla}\phi\|^{2}$ and a nonconforming part $\|a^{-1/2}{{\nabla\times}}\psi\|^{2}$. It then suffices to obtain an a posteriori error bound for each part separately and sum to obtain an estimator for the total error.

6.1. Proof of Theorem 1

The proof of Theorem 1 follows from [10] for the symmetric interior penalty discontinuous Galerkin methods almost verbatim. Specifically, the upper bound (32) in Theorem 1 follows from (54a) and (54b) below, and the lower bound (33) follows from Lemma 2, (54) and (54d).

The following estimates follows from results in [10, Lemma 6.2-6.5] using the proof in [9, Section 6]. Let $\phi$ and $\psi$ be taken as in (50) in the case where ${\boldsymbol{\tau}}=a{\nabla}u-a{\nabla}u_{h}$, and let $\eta_{CF,K}$ and $\eta_{NC,K}$ be given by (31). Then,

[TABLE]

where $\widetilde{\mathcal{E}}_{K}=\{F\in\partial\mathcal{T}_{h}:\overline{F}\cap\overline{K}\text{ is nonempty}\,\}$.

6.2. Proof of Theorem 2

Theorem 2 is a consequence of the following three lemmas:

Lemma 6**.**

Let $\phi$ and $\psi$ be given by in the decomposition (50), where ${\boldsymbol{\tau}}$ is chosen to be ${\boldsymbol{\tau}}=a{\nabla}u-{\boldsymbol{\sigma}}_{h}$, and let $\eta_{CF,K}$ and $\eta_{NC,K}$ be given by (47). Then,

[TABLE]

where $\widetilde{\mathcal{E}}_{K}=\{F\in\partial\mathcal{T}_{h}:\overline{F}\cap\overline{K}\text{ is nonempty}\,\}.$

Proof.

Direct computation gives

[TABLE]

where the last inequality follows from the Cauchy-Schwarz and the Poincaré inequalities. This completes the proof of (55a). Turning to (55b), since $({\nabla}v,{{\nabla\times}}\psi)=0$ for all $v\in H^{1}_{0}(\Omega)$, we have

[TABLE]

Let ${\boldsymbol{\rho}}_{h}={\boldsymbol{\sigma}}_{h}^{*}-{\boldsymbol{\sigma}}_{h}$, then for all $K\in\mathcal{T}_{h}$, by (44), we have ${\boldsymbol{\rho}}_{h}\in\mathbb{P}_{k+1}(K)^{d}$ satisfies

[TABLE]

Hence, there exists a constant $c$, depending only on the polynomial degree $k$ and shape-regularity of the element $K$, such that

[TABLE]

while, using a standard bubble function technique [7, 50], we have

[TABLE]

The choice of stabilization parameter (37) means that

[TABLE]

and the proof of (55) follows from these estimates.

Finally, we have

[TABLE]

Combining the above estimates completes the proof of (55d). ∎

Lemma 7**.**

There exists a positive constant $c$, depending only on the shape-regularity of the mesh and the polynomial degree $k$, such that for any facet $F\in\mathcal{E}(K)$,

[TABLE]

Proof.

The mixed HDG scheme (36) satisfies, for every $K\in\mathcal{T}_{h}$,

[TABLE]

Let the function $z\in\mathbb{P}_{k}(K)$ satisfy

[TABLE]

where $F^{*}$ is a fixed facet of $K$. We have, by a standard scaling argument,

[TABLE]

Taking $v=z$ in (59) and rearranging terms, we obtain

[TABLE]

The proof is completed by invoking estimate (57) for the cell-wise residual term $f+{\nabla\cdot}{\boldsymbol{\sigma}}_{h}$. ∎

Lemma 8**.**

*There exists a positive constant $c$, depending only on the shape-regularity of the mesh and the polynomial degree $k$, such that *

[TABLE]

Proof.

We first prove the estimate (60). We denote ${\boldsymbol{\rho}}_{h}={\nabla}u_{h}^{*,\mathrm{dc}}-a^{-1}{\boldsymbol{\sigma}}_{h}$ and let the function ${\boldsymbol{\rho}}_{h}^{*}\in\mathbb{P}_{k+1}(K)^{d}\oplus{\boldsymbol{x}}\widetilde{\mathbb{P}}_{k+1}(K)$ be defined as follows:

[TABLE]

then, we have

[TABLE]

Moreover, by equations (46a) and (62), we have

[TABLE]

This implies that ${\nabla\cdot}{\boldsymbol{\rho}}_{h}^{*}=0$ because ${\nabla\cdot}{\boldsymbol{\rho}}_{h}^{*}\in\mathbb{P}_{k+1}(K)$. Since ${\boldsymbol{\rho}}_{h}^{*}$ has vanishing normal trace on ${\partial K}$ by equations (62a), we obtain

[TABLE]

Hence,

[TABLE]

The estimate (60) now follows from the Cauchy-Schwarz inequality.

Let $P_{M_{0}}$ be the $L^{2}$-projection onto the space $\mathbb{P}_{0}(F)$. Applying the results in [27, Lemma 3.4-3.5], we obtain

[TABLE]

The estimate (61) now immediately follows from the triangle inequality and (60).

∎

6.3. Proof of Lemma 1

Since $\mathbb{V}_{h,k,\delta}$ is a finite-dimensional space, it suffices to show that $(u_{h},\widehat{u}_{h})=(0,0)$ is the only solution to the homogeneous problem. Let

[TABLE]

then for $(v_{h},\widehat{v}_{h})\in\mathbb{V}_{h,k,\delta}$, we have

[TABLE]

Since a{\nabla}v_{h}\cdot{\boldsymbol{n}}\raise-0.86108pt\hbox{|}_{F}\in\mathbb{P}_{k-1}(F), there holds

[TABLE]

By the Cauchy-Schwarz and Young’s inequalities, for any $\epsilon>0$,

[TABLE]

where, the final inequality holds thanks to the inverse-trace inequality [52] and a{\nabla}u\raise-0.86108pt\hbox{|}_{K}\in\mathbb{P}_{k-1}(K)^{d}. Summing the above inequality over $F\in\mathcal{E}(K)$ gives

[TABLE]

Hence, we have

[TABLE]

Finally, if $\gamma$ satisfies (15) then there exists an $\epsilon$ such that

[TABLE]

Consequently, when the right hand side of (13) vanishes, $\|{\nabla}v_{h}\|_{K}=0$ for all $K\in\mathcal{T}_{h}$ and $\|P_{M}v_{h}-\widehat{v}_{h}\|_{F}=0$ for all $F\in\mathcal{E}_{h}$. Hence, the only solution to the homogeneous problem is $(v_{h},\widehat{v}_{h})=(0,0)$, and therefore there exists a unique solution $(u_{h},\widehat{u}_{h})\in\mathbb{V}_{h,k,\delta}$ to (13).

6.4. Proof of Lemma 2

The proof follows that of [2, Theorem 3].

By the conservation property (23), we have

[TABLE]

Hence,

[TABLE]

Inserting the above expression into the jump term in the HDG energy norm and regrouping gives

[TABLE]

Here, the gradient jump term can be controlled by the standard bubble function technique [50, 7]

[TABLE]

On the other hand, there holds

[TABLE]

where $\overline{\,[\![P_{M}u_{h}]\!]}$ denote the average value of $\,[\![P_{M}u_{h}]\!]$ on the facet $F$. Thanks to the trace and Poincaré inequalites, we have

[TABLE]

Hence, to show norm equivalence, it remains to show that the term

[TABLE]

can be controlled by the discrete energy seminorm plus the data oscillation. Replacing $\widehat{u}_{h}$ in the primal HDG scheme (13) with the expression (66), we obtain

[TABLE]

for all $v_{h}\in V_{h,k}$ and $K\in\mathcal{T}_{h}$.

Moreover, there holds [49]

[TABLE]

where ${\boldsymbol{S}}_{K}$ is the element stiffness matrix, i.e. ${\boldsymbol{S}}_{ij}=({\nabla}\lambda_{i},{\nabla}\lambda_{j})_{K}$ with $\{\lambda_{\ell}\}_{\ell=1}^{d+1}$ being the barycentric coordinates for the element $K$, and $\rho({\boldsymbol{S}}_{K})$ is its spectral radius. Hence the stabilization parameter $\alpha_{h}$ in (14), with $\gamma$ satisfying (15), satisfies

[TABLE]

The proof is then concluded following [2, Theorem 3] by taking special linear test functions in the equation (6.4) and using the above estimate for the stabilization parameter.