Robust a Posteriori Error Estimates for HDG method for Convection-Diffusion Equations

# Robust a Posteriori Error Estimates for HDG method for Convection-Diffusion Equations

Huangxin Chen School of Mathematical Sciences and Fujian Provincial Key Laboratory on Mathematical Modeling & High Performance Scientific Computing, Xiamen University, Fujian, 361005, P.R. China Jingzhi Li Faculty of Science, South University of Science and Technology of China,Shenzhen, 518055, China  and  Weifeng Qiu Department of Mathematics, City University of Hong Kong, 83 Tat Chee Avenue, Kowloon, Hong Kong, China
###### Abstract.

We propose a robust a posteriori error estimator for the hybridizable discontinuous Galerkin (HDG) method for convection-diffusion equations with dominant convection. The reliability and efficiency of the estimator are established for the error measured in an energy norm. The energy norm is uniformly bounded even when the diffusion coefficient tends to zero. The estimators are robust in the sense that the upper and lower bounds of error are uniformly bounded with respect to the diffusion coefficient. A weighted test function technique and the Oswald interpolation are key ingredients in the analysis. Numerical results verify the robustness of the proposed a posteriori error estimator.

###### Key words and phrases:
hybridizable discontinuous Galerkin method, a posteriori error estimates, convection-diffusion equations
###### 2000 Mathematics Subject Classification:
65N30, 65L12
The authors would also like to thank the associate editor and all the referees for constructive criticism leading to a better presentation of the material in this paper. The first author would like to thank the support from the City University of Hong Kong where this work was carried out during his visit, and he also thanks the supports from the NSF of China (Grant No. 11201394) and the NSF of Fujian Province (Grant No. 2013J05016). The work of the second author was supported by the NSF of China (Grant No. 11201453). The work of the third author was supported by the GRF of Hong Kong (Grant No. 9041980).

## 1. Introduction

Given a bounded, polyhedral domain , we consider the convection-diffusion equations

 −ϵΔu+β⋅∇u+cu= f in Ω, (1.1a) u= g on ∂Ω. (1.1b)

The data and the right-hand sides in (1.1) satisfy the following assumptions:

1. .

2. , , and .

3. .

4. There is a function and a positive constant such that .

Assumption (A1) includes the case of the convection-dominated regime. According to [4], Assumption (A4) is satisfied if has no closed curves and .

It is well known that solutions of (1.1) may develop layers (cf. [24, 29]). In particular, the solutions may have singular interior layer of width or outflow layer of width . Standard numerical methods, e.g., standard finite element method or central finite difference method, are not robust when the quantity is small compared to the mesh size. In order to stabilize the numerical method, several remedies are proposed for addressing the issue, for instance, streamline diffusion method [6], residual free bubble methods [7, 8, 10], local projection schemes [37], subgrid scale method [16, 5], continuous interior penalty (CIP) methods [11, 12], discontinuous Galerkin methods [4, 32, 33], and recently discontinuous Petrov-Galerkin (DPG) methods [9, 13, 23], HDG method [31] and the first order least squares method [15]. One can refer to [40, 44] for more other stabilization techniques. But in order to capture the potential interior or outflow layer of the solutions to the problem (1.1), the local Péclet number near the layers should be small enough, where is the mesh size. Hence it would be quite expensive for the stabilized numerical methods used on the quasi-uniform mesh to capture the layers when is small. If the mesh in the vicinity of the layers can be locally refined, the cost of numerical computations could be reduced. Therefore the adaptive finite element method is a natural choice for the efficient solution of convection-diffusion equations with dominant convection.

The adaptive finite element method based on a posteriori error estimates have been well established for second-order elliptic problems (cf. [2, 47]). In recent years the a posteriori error estimates are also extended to convection-diffusion equation. An early attempt was proposed by Eriksson and Johnson in [26], using regularization and duality techniques. Verfürth [48] proposed semi-robust estimators in the energy norm for the standard Galerkin approximation and the streamline upwind Petrov-Galerkin (SUPG) discretization. In [49] Verfürth improved his results by giving the estimates which are robust with dominant convection in a norm incorporating the standard energy norm and a dual norm of the convective derivative. Very recently, Tobiska and Verfürth [45] derived the same robust a posteriori error estimators for a wide range of stabilized finite element methods such as streamline diffusion methods, local projection schemes, subgrid scale technique and CIP method. However, the energy norm of error used in Verfürth’s estimates is defined through a dual norm which is not easy to compute. Sangalli [41] proposed different norms for the a posteriori error estimates that allow for robust estimators, but the analysis is only valid in the one dimensional case. As to other approaches for the robust error estimations, one can refer to [50] for mixed finite element methods, [51] for cell-centered finite volume scheme, [3] for nonconforming finite element method, [27, 28, 42, 52] for interior penalty discontinuous Galerkin method.

Discontinuous Galerkin (DG) methods have several attractive features compared with conforming finite element methods. For example, DG methods have elementwise conservation of mass, and they work well on arbitrary meshes. However, the dimension of the approximation DG space is much larger than the dimension of the corresponding conforming space. The HDG method [18, 19, 36] was recently introduced to address this issue. HDG methods retain the advantages of standard DG methods and result in significant reduced degrees of freedom. New variables on all interfaces of mesh are introduced such that the numerical solution inside each element can be computed in terms of them, and the resulting algebraic system is only due to the unknowns on the skeleton of the mesh. In [31], the HDG method was proposed and analyzed for the problem (1.1) on shape-regular mesh. The stabilization parameter of the HDG method in [31] can be determined clearly, meanwhile the penalty parameters of the DG schemes [4] need to be chosen empirically. Moreover, the condition number of stiffness matrix of a new way of implementing HDG method in [31] was proven to be bounded by and independent of the diffusion coefficient . These properties are important for the efficient solution of the problem (1.1) and encourage us to consider the corresponding HDG method on adaptive meshes.

The a posteriori error analysis for the HDG method for second order elliptic problems has been presented in [20, 21], where the error incorporates only the flux and a postprocessed solution used in the estimators. To our best knowledge, no a posteriori error estimates for the HDG discretizations of convection-diffusion problems have been studied in the literature so far. In this work, our objective is to show that the HDG scheme proposed in [31] gives rise to robust a posteriori error estimates for the problem (1.1). In comparison with the postprocessing technique utilized in [50, 20, 21], we establish the estimators without any postprocessed solution since the solution of (1.1) is always not smooth and there is no superconvergence result for the HDG method when .

We notice that the a posteriori error estimators for nonconforming finite element method [3] and interior penalty DG method [27, 28] are only semi-robust in the sense that they yield lower and upper bounds of the error which differ by a factor equal at most to . In [42, 52], the a posteriori error estimator is robust in the sense that the ratio of the constants in the upper and lower bounds of error is independent of the diffusion coefficient. However, the energy norm of error in [42, 52] contains the jump term on each interior interface of meshes. In contrast, the a posteriori error estimator in this paper is robust based on the energy norm in (2.9), which contains a jump term instead. One can refer to (2.7) for the definition of paparemeter . So, our a posteriori error estimator will not enlarge the error estimate too much as the error estimator in [42, 52], when the mesh size is not small compared with the diffusion coefficient.

To derive the reliability and efficiency of the estimators for the error measured in an energy norm which incorporates a scaling flux and the scalar solution of the HDG discretization, two techniques are utilized. The first one is to use the Oswald interpolation operator to approximate a discontinuous polynomial by a continuous and piecewise polynomial function and to control the approximation by the jumps (cf. [34, 35]). For most of a posteriori error estimates mentioned above (e.g. [49, 50, 45, 27, 28, 42, 52]), the analysis only gives the estimates for the energy error without -error of the scalar solution when . The second one is to address this issue and to employ a weighted function to derive the estimates for the error which contains -error of the scalar solution. This idea goes back to Nävert’s work [46] for convection-diffusion problems and was used to obtain the -stability of the original DG method for pure hyperbolic equation [39] and extended to convection-diffusion equations using the IP-DG method [30, 4], the HDG method [31] and the first order least squares method [15].

In the numerical experiments, the convection-diffusion problems with interior or outflow layers are tested based on the proposed a posteriori error estimator. The robustness of the a posteriori error estimator based on the HDG method is observed for the problems with different diffusion coefficient. We also find that the convergence of the adaptive HDG method is almost optimal, i.e., the convergence rate is almost , where is the number of elements, depends on the polynomial order .

The outline of the paper is as follows: We introduce some notations, the HDG method, a posteriori error estimator and main results in the next section. In section 3, we collect some auxiliary results for analysis. Section 4 and section 5 are devoted to the proofs of reliability and efficiency, respectively. In the final section, we give some numerical results to confirm our theoretical analysis.

## 2. Notation, HDG method, error estimator, and main results

In this section, we begin with some basic notation and hypotheses of meshes. Secondly, we introduce the HDG method for (1.1) in [31]. Then, we define the corresponding a posteriori error estimator. Finally, we give the main results of reliability and efficiency.

### 2.1. Notation and the mesh

Let be a conforming, shape-regular simplicial triangulation of . For any element , denotes the set of its edges in the two dimensional case and of its faces in the three dimensional case. Elements of will be generally referred to as faces, regardless of dimension, and denoted by . We define . We denote by the set of all faces in the triangulation (the skeleton), while the set of all interior (boundary) faces of the triangulation will be denoted (). Correspondingly, we refer to the set of vertices and to the set of interior vertices. For any , let be the diameter of element . Similarly, for any , we define . Throughout this paper, we use the standard notations and definitions for Sobolev spaces (see, e.g., Adams[1]). We also use the notation and to denote the -norm on the elements and faces , respectively.

### 2.2. The HDG method

The HDG method is based on a first order formulation of the convection-diffusion equation (1.1), which can be rewritten in a mixed form as finding such that

 ϵ−1q+∇u =0in Ω, (2.1a) divq+β⋅∇u+cu =fin Ω, (2.1b) u =gon ∂Ω. (2.1c)

For any element and any face , we define

 V(T):=(Pp(T))d,W(T):=Pp(T),M(F):=Pp(F),

where is the space of polynomials of total degree not larger than on . The finite element spaces are given by

 Vh: ={v∈L2(Ω):v|T∈V(T) for all T∈Th}, Wh: ={w∈L2(Ω):w|T∈W(T) for all T∈Th}, Mh: ={μ∈L2(Eh):μ|F∈M(F) for all F∈Eh}, Mh(g): ={μ∈Mh:∫∂Ω(μ−g)ξds=0 for all ξ∈Mh},

where and .

The HDG method seeks finite element approximations satisfying

 (ϵ−1qh,r)Th−(uh,divr)Th+⟨ˆuh,r⋅n⟩∂Th=0, (2.2a) −(qh+βuh,∇w)Th+((c−divβ)uh,w)Th+⟨(ˆqh+ˆβuh)⋅n,w⟩∂Th (2.2b) =(f,w)Th, ⟨ˆuh,μ⟩∂Ω=⟨g,μ⟩∂Ω, (2.2c) ⟨(ˆqh+ˆβuh)⋅n,μ⟩∂Th∖∂Ω=0, (2.2d)

for all , where the normal component of numerical flux is given by

 (ˆqh+ˆβuh)⋅n=qh⋅n+(β⋅n)ˆuh+τ(uh−ˆuh)on ∂Th, (2.3)

and the stabilization function is a piecewise, nonnegative constant defined on . Here, we define and . One of the advantages of the HDG method is the elimination of both and from the system (2.2) to obtain a formulation in terms of numerical trace only, one can refer to [17, 31, 36, 38] for the implementation.

The stabilization function in (2.3) is chosen as

 τ|F=max(supx∈Fβ(x)⋅n,0)+min(ρ0ϵhT,1)∀F∈∂T, T∈Th. (2.4)

Here, . We emphasize that the choice of in (2.4) is the second type of stabilization function in [31]. According to [14], the HDG method (2.2) with stabilization function (2.4) has a unique solution. Compared with a recent work [15], the intrinsic idea of choosing in the above HDG method is similar to the strategy of setting the ultra-weakly imposed boundary condition in the first order least squares method for (1.1).

### 2.3. A posteriori error estimator

We define the elementwise residual function as

 Rh|T=f−divqh−β⋅∇uh−cuh,∀T∈Th. (2.5)

We define

 αS=min{hSϵ−12,1}, (2.6)

where can be any element or any face . In addition, for any , we introduce

 γF=min{ϵhF+(hFϵ+ϵ−12αF)∥β∥L∞(F)+hF,ϵ+∥β∥L∞(F)hF+hF}. (2.7)

Now, we are ready to introduce the a posteriori error estimator in the following.

###### Definition 2.1.

(A posteriori error estimator) The a posteriori error estimator is defined as

 η =(ΣT∈Thη2T+ΣF∈E0h(η0F)2+ΣF∈E∂h(η∂F)2)12 where (2.8a) ηT =(α2T∥Rh∥20,T+ϵ−1∥qh+ϵ∇uh∥20,T)12,∀T∈Th, (2.8b) η0F =(ϵ−12αF∥\llbracketqh⋅n\rrbracket∥20,F+γF∥\llbracketuh\rrbracket∥20,F)12,∀F∈E0h, (2.8c) η∂F =γ12F∥uh−g∥0,F,∀F∈E∂h. (2.8d)

Here, for any interior face in , we define the jump of scalar function and the jump of normal component of vector field by

 \llbracketϕ\rrbracket=ϕ+−ϕ−,\llbracketσ⋅n\rrbracket=σ+⋅n++σ−⋅n−, % respectively.

### 2.4. Reliability and efficiency of a posteriori error estimator

From now on, we use to denote generic constants, which are independent of the diffusion coefficient and the mesh size.

For any , we define the energy norm

 \interleave(p,w)\interleave2h (2.9) = ∑T∈Th(ϵ−1∥p∥20,T+∥w∥20,T+ϵ∥∇w∥20,T+α2T∥divp+β⋅∇w∥20,T) +∑F∈E0h(ϵ−12αF∥\llbracketp⋅n\rrbracket∥20,F+γF∥\llbracketw\rrbracket∥20,F)+∑F∈E∂hγF∥w∥20,F.

We outline main results by showing the reliability and efficiency of the a posteriori error estimator, respectively, in the following theorems. Using the techniques developed in this paper, we would like to emphasize that all the a posteriori error estimates in this paper will also hold for the mixed hybrid method in [25].

###### Theorem 2.2.

(Reliability) If the Dirichlet data is contained in , then

 \interleave(q−qh,u−uh)\interleaveh≤C0η. (2.10)
###### Remark 2.1.

In the diffusion dominated case and , the parameter in (2.7) would be , and the proposed a posterior error estimators coincide with the ones in [21]. In this particular case, the total energy norm can be defined as .

###### Theorem 2.3.

(Efficiency) We have

 η0F ≤C1(ϵ−12αF∥\llbracket(q−qh)⋅n\rrbracket∥20,F+γF∥\llbracketu−uh\rrbracket∥20,F)12,∀F∈E0h, (2.11a) η∂F ≤C1(γF∥g−uh∥20,F)12,∀F∈E∂h, (2.11b) ηT ≤C1(ϵ−1∥q−qh∥20,T+∥u−uh∥20,T+ϵ∥∇(u−uh)∥20,T (2.11c) +α2T∥div(q−qh)+β⋅∇(u−uh)∥20,T+osc2h(Rh,T))12,∀T∈Th.

Here, the data oscillation term , and is the orthogonal projection onto .

In addition, for any with or any with , we have the following efficiency results with explicit dependence on .

###### Theorem 2.4.

(Efficiency on refined element) For any , if , then

 (η0F)2 ≤C2∑T∈ωF(ϵ−1∥q−qh∥20,T+ϵ∥∇(u−uh)∥20,T (2.12) +∥u−uh∥20,T+osc2h(Rh,T)).

For any , if , then

 η2T ≤C2(ϵ−1∥q−qh∥20,T+ϵ∥∇(u−uh)∥20,T (2.13) +∥u−uh∥20,T+osc2h(Rh,T)).

Here, is the union of elements sharing the common face , and is the data oscillation term introduced in Theorem 2.3.

###### Remark 2.2.

For the diffusion dominated case , and , the above efficiency estimates (2.12) and (2.13) also coincide with the efficiency results in [21].

## 3. Auxiliary results

We collect some auxiliary results in this section for the proof of reliability and efficiency.

For every element , we denote by the union of all elements that share at least one point with . For any face , the set is defined analogously, meanwhile, is defined to be the union of elements sharing the common face . The following Clément-type interpolation is crucial for the proof of reliability.

###### Definition 3.1.

(cf. [47]) One can define a linear mapping via

 πhv:=∑z∈N0h(1|Ωz|∫Ωzvdx)ϕz,

where is nodal bases function for every vertex , is the support of a nodal bases function which consists of all elements that share the vertex , and is the corresponding conforming finite element space defined by

The interpolation in Definition 3.1 has the following approximation properties.

###### Lemma 3.2.

(cf. [48, 49]) For any and , the following estimates hold for any function :

 ∥(I−πh)v∥0,T ≤C1αT(∥v∥20,ΩT+ϵ∥∇v∥20,ΩT)12, (3.1a) ∥(I−πh)v∥0,F ≤C2ϵ−14α12F(∥v∥20,ΩT+ϵ∥∇v∥20,ΩT)12. (3.1b)
###### Remark 3.1.

The proof of Lemma 3.2 can be obtained by the Lemmas 3.1 and 3.2 in [48]. For the particular case , the constant is set to be for any or in [49], and the associated energy error excludes the -error . In this paper, a weighted function technique used in [4] shall be employed, such that we can also obtain the estimates for the energy error including the -error even when . Hence, is always set as (2.6).

For the approximation of function in and the Dirichlet boundary data by continuous finite element space, we need to introduce Oswald interpolation. If is contained in , then the continuous finite element space is not empty. Hence, we can introduce the Oswald interpolation . Given a function , the operator is prescribed at the Lagrangian nodes in the interior of by the average of the values of at this node. For the nodes at the boundary , is prescribed at the Lagrangian nodes on by the value of at this node. The following estimate has been analyzed for nonconforming mesh and conforming mesh in [34, 35] and was extended to variable polynomial degree in [53].

###### Lemma 3.3.

(cf. [34, 35, 53, 21]) If the Dirichlet data is contained in , then for any and any multi-index with , the following estimate holds:

 ∑T∈Th∥Dα(vh−Ioshvh)∥20,T (3.2) ≤ C3⎛⎜⎝∑F∈E0hh1−2|α|F∥\llbracketvh\rrbracket∥20,F+∑F∈E∂hh1−2|α|F∥g−vh∥20,F⎞⎟⎠.

In order to prove the local efficiency of the estimators, proper element and face bubble functions are useful. We set the element bubble function in the element as , where denotes the linear nodal basis function at th vertex in . Besides the element bubble function, as mention in Lemma 3.3 in [48] and Lemma 3.6 in [49], there also exists proper face bubble function with such that the following lemma holds.

###### Lemma 3.4.

(cf. [48, 49]) For any element , a polynomial and any face , a polynomial , the following estimates hold:

 ∥ϕ∥20,T ≤C4(ϕ,BTϕ)T, ∥BTϕ∥0,T ≤C5∥ϕ∥0,T, ∥ψ∥20,F ≤C6⟨ψ,BFψ⟩F, ∥BFψ∥0,ωF ≤C7ϵ14α12F∥ψ∥0,F, ∥BFψ∥0,ωF+ϵ12∥∇BFψ∥0,ωF ≤C8ϵ14α−12F∥ψ∥0,F.

## 4. Proof of reliability

In this section, we give the proof of Theorem 2.2, which shows reliability of the a posteriori error estimator in Definition 2.1.

In view of the assumption (A4), we define a weighted function

 φ:=e−ψ+χ, (4.1)

where is a positive constant to be determined later. Let . We have the following Lemma 4.1.

###### Lemma 4.1.

Let be given in (4.1) with . Then the following estimate holds:

 C(ϵ−1∥eq∥2Th+∥eu∥2Th) (4.2) ≤ ϵ−1(eq,φeq)Th−(eu,∇φ⋅eq)Th −12(β⋅∇φeu,eu)Th+((c−12divβ)eu,φeu)Th.
###### Proof.

According to the assumptions (A3)-(A4), we have

 ϵ−1(eq,φeq)Th+(eu,e−ψ∇ψ⋅eq)Th +12(β⋅∇ψe−ψeu,eu)Th+((c−12divβ)eu,φeu)Th ≥ ϵ−1χ(eq,eq)Th+(eu,e−ψ∇ψ⋅eq)Th+b02(e−ψeu,eu)Th.

By the Cauchy-Schwarz and Young’s inequalities, for any , we have

 (eu,e−ψ∇ψ⋅eq)Th ≤ 12(δ−1∥∇ψ∥2L∞(Ω)∥e−ψ∥L∞(Ω)(eq,eq)Th+δ∥e−ψ∥L∞(Ω)(eu,eu)Th).

Then, by taking and , we can conclude that the proof is complete. ∎

We are now ready to state a key result of the upper bound estimate of .

###### Lemma 4.2.

If the Dirichlet data is contained in , then

 (ϵ−1∥eq∥2Th+∥eu∥2Th)≤C(∑T∈Thη2T+∑F∈E0h(η0F)2+∑F∈E∂h(η∂F)2). (4.3)
###### Proof.

According to Lemma 4.1, we have

 C(ϵ−1∥eq∥2Th+∥eu∥2Th) (4.4) ≤ ϵ−1(eq,φeq)Th−(eu,∇φ⋅eq)Th −12(β⋅∇φeu,eu)Th+((c−12divβ)eu,φeu)Th.

Let . By the definition of , clearly, we have . Adding and subtracting into , we have

 ϵ−1(eq,φeq)Th (4.5) = −(q−qh,φ(∇u−∇uIh))Th−ϵ−1(q−qh,φ(ϵ∇uIh+qh))Th = (∇φ⋅(q−qh),u−uIh)Th+(φdiv(q−qh),u−uIh)Th −⟨(q−qh)⋅n,φ(u−uIh)⟩∂Th−ϵ−1(q−qh,φ(ϵ∇uIh+qh))Th = (∇φ⋅(q−qh),u−uIh)Th+(Rh,φ(u−uIh))Th +⟨qh⋅n,φ(u−uIh)⟩∂Th−(β⋅∇(u−uh)+c(u−uh),φ(u−uIh))Th −ϵ−1(q−qh,φ(ϵ∇uIh+qh))Th.

In the last step of (4.5), we have used the equation (2.1b) and the fact that (since jumps of vanish on all interior faces and on ). Inserting (4.5) into (4.4) and subtracting and adding into in the first term of the right-hand side of (4.5), we have

 C(ϵ−1∥eq∥2Th+∥eu∥2Th) (4.6) ≤ (∇φ⋅(q−qh),uh−uIh)Th+(Rh,φ(u−uIh))Th +⟨qh⋅n,φ(u−uIh)⟩∂Th−(β⋅∇(u−uh)+c(u−uh),φ(u−uIh))Th −ϵ−1(q−qh,φ(ϵ∇uIh+qh))Th −12(β⋅∇φeu,eu)Th+((c−12divβ)eu,φeu)Th.

For any , the equation (2.2b) in the HDG method (2.2) can be rewritten as follows after integration by parts:

 (f,w)Th =−(qh,∇w)Th+(div(βuh),w)Th−⟨βuh⋅n,w⟩∂Th +(cuh−divβuh,w)Th+⟨(ˆqh+ˆβuh)⋅n,w⟩∂Th =(divqh,w)Th−⟨qh⋅n,w⟩∂Th+(β⋅∇uh+cuh,w)Th −⟨β⋅nuh,w⟩∂Th+⟨(ˆqh+ˆβuh)⋅n,w⟩∂Th.

Note that the equation (2.2d) indicates . Hence, we have

 −⟨(qh+βuh)⋅n,w⟩∂Th=(Rh,w)Th. (4.7)

Combining (4.6) and (4.7), we obtain

 C(ϵ−1∥eq∥