Analysis of output-based error estimation for finite element methods

Applied Numerical Mathematics 118 (2017) 182–202

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Analysis of output-based error estimation for finite element methods Hugh A. Carson ∗ , David L. Darmofal, Marshall C. Galbraith, Steven R. Allmaras Massachusetts Institute of Technology, 77 Massachusetts Ave, Cambridge, MA 02139, United States

a r t i c l e

i n f o

Article history: Received 17 October 2016 Received in revised form 19 January 2017 Accepted 8 March 2017 Available online 16 March 2017 Keywords: Output functional estimation Finite element method Discontinuous Galerkin Hybridizable Discontinuous Galerkin

a b s t r a c t In this paper, we develop a priori estimates for the convergence of outputs, output error estimates, and localizations of output error estimates for Galerkin finite element methods. Output error estimates for order p finite element solutions are constructed using the DualWeighted Residual (DWR) method with a higher-order p  > p dual solution. Specifically, we analyze these DWR estimates for Continuous Galerkin (CG), Discontinuous Galerkin (DG), and Hybridized DG (HDG) methods applied to the Poisson problem. For all discretizations, as h → 0, we prove that the output and output error estimate converge at order 2p and 2p  (assuming sufficient smoothness), while localizations of the output and output error estimate converge at 2p + d and p + p  + d. For DG, the results use a new post processing for the error associated with the lifting operator. For HDG, these rates improve an additional order when the stabilization is based upon an O (1) length scale. © 2017 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction A posteriori estimation aims to provide computable estimates of numerical error in the discrete approximate solutions of Partial Differential Equations (PDEs). The computed a posteriori estimate can then be used to correct the output and/or drive mesh adaptation. Energy-based a posteriori error estimates have been used in the Finite Element Method (FEM) since the 1980s [9,11], particularly within the structural community given the physical relevance of the energy norm. Ainsworth and Oden provided a comprehensive summary of energy norm estimates (as well as other norms) wherein they draw attention to the use of duality arguments for analyzing the estimates in an a posteriori context [2]. Methods based upon localized residuals were developed for mesh adaptation by Babuška and Rheinboldt [7]. Estimates based on the solution of local Dirichlet problems and local Neumann problems were also developed by Bank and Weiser [11] and Babuška and Rheinboldt [8]. Verfürth analyzed these for elliptic partial differential showing their equivalence [35]. In many equations, the convergence rate of the energy norm is less than that of an output functional. In a series of papers, Babuška and Miller showed that integrated output functionals from FEM could exhibit ‘super-convergence’, converging faster than the solution error (e.g. the L 2 norm of the error) [4–6]. Barrett and Elliott analyzed the linear convection– diffusion equation and reached similar conclusions [12]. The Dual-Weighted Residual (DWR) method developed for FEM by Becker and Rannacher relates the primal residual error to the computed output functional error by weighting it with an adjoint, or dual [14,15]. Adjoint based error estimation has

*

Corresponding author. E-mail addresses: [email protected] (H.A. Carson), [email protected] (D.L. Darmofal), [email protected] (M.C. Galbraith), [email protected] (S.R. Allmaras). http://dx.doi.org/10.1016/j.apnum.2017.03.004 0168-9274/© 2017 IMACS. Published by Elsevier B.V. All rights reserved.

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

183

also been developed for discretizations other than FEM, for example Finite Volume (FV) and Finite Difference (FD) methods [28–30,33]. For the purposes of mesh adaptation, error estimates must be localized to facilitate decisions about where mesh refinement (and coarsening) should be performed. A typical element-wise localization for the Continuous Galerkin (CG) FEM uses a strong form residual, which introduces an element boundary contribution due to gradient jumps [7,14,34]. For DG, localization by elemental restriction of the test function has been analyzed by Yano [36] for the second DG discretization of Bassi and Rebay (BR2) [13]. For a more complete overview of output-based error estimation and earlier mesh adaptation schemes, for fluid dynamic applications in particular, the review papers of Hartmann and Houston [31], focused on DG methods, and Fidkowski and Darmofal [27], for general discretizations, provide a comprehensive literature review. The focus of this paper is a priori analysis of the convergence rates for DWR error estimates for the Continuous Galerkin (CG), Discontinuous Galerkin (DG) and Hybridizable Discontinuous Galerkin (HDG) schemes. The framework we utilize for this analysis is new, as are the convergence results for DG and HDG estimators. For DG, a new treatment for the lifting operator that appears in many DG discretizations is shown to improve the convergence rates of the global error estimate. The convergence of the CG estimator has been studied previously [35], and so we include it here in Appendix A for completeness though the result is not new. The assumptions used in this paper do not change the key results, but avoid complexities whilst achieving the desired goal of proving convergence rates for the discretizations. We assume sufficient smoothness in line with the literature [1–3, 10,21,23,25,32,33]. The relaxation of this would give the well-known results that limit the convergence rates achievable [1,25]. The Poisson problem with homogeneous Dirichlet boundary conditions is chosen for a similar reason, in that the inclusion of Neumann or Robin boundary conditions or the addition of a convection or reaction term would not change the need for a coercivity assumption, and with such an assumption in place the convergence is dominated by the diffusive term, again this is a standard model PDE used throughout the literature [2,3,9,14,32,33]. The most nuanced assumption of the paper when considering an adaptive meshing context, is that of quasi-uniformity, though the focus of this paper is not the application of these estimates and their localization. The goal of this paper is to demonstrate the estimates meet the minimum requirement of converging as the mesh is refined in some fashion, for which we have chosen quasi-uniform refinement. For non-adaptive convergence analysis, quasi-uniform (or locally quasi-uniform [2]) refinement is a standard assumption [1,2,8–10,21,23,32]. Even in the adaptive context the assumption of quasi-uniformity may be applicable, especially in regular regions. Specifically, for a sufficiently regular solution, the optimal mesh distribution is a quasi-uniform refinement in the limit of infinite resolution. An example of this can be most clearly seen in [14] where the authors give a selection of feedback algorithms which will all ultimately achieve quasi-uniform meshes by equilibrating the error. Finally we note that the results hold with more general refinement criteria than quasi-uniform, specifically any refinement in which the maximum element size over the mesh tends to zero, i.e. h = maxκ ∈Th hκ → 0, in which case several of the local estimates could be written to track both a local elemental size (hκ ) and a global mesh size (h) dependence. The outline of the paper is as follows, in Section 2 the model Poisson primal and adjoint problems are outlined, alongside definitions that are necessary for the results in the following sections. Then, in Sections 3 and 4, a priori asymptotic convergence rates for the error, and the error in the estimate, are derived for the DG and HDG discretizations respectively. In Section 5 numerical results are presented for a two dimensional test case demonstrating the derived a priori asymptotic convergence rates for all of the discretizations and estimates considered. 2. The Poisson problem Consider the primal Poisson problem:

−u = f in 

(2.1)

u = 0 on ∂, where  ⊂ Rd is a bounded domain of dimension d. The solution u is assumed to be sufficiently smooth, which places restrictions on the domain shape and f (see for example Brenner and Scott [17]). We denote the L 2 inner product by (·, ·) and the corresponding norm over  by  ·  L 2 () . In addition inner products on a trace, such as ∂ uv, will be denoted by ·, · . Where additional clarity is useful, subscripts will be used in order to denote the domain of the integral. The problem statement is written in the variational form,

u ∈ V : R ( u , w ) ≡ l ( w ) − a ( u , w ) = 0,

∀w ∈ V

a( v , w ) = (∇ v , ∇ w ) l( w ) = ( f , w ). The boundary conditions are enforced essentially by V where, for this problem, V ≡ H 01 (). In approximating this problem discretely using a Galerkin method, we utilize a finite-dimensional space Vh, p and a modified problem statement

u h, p ∈ Vh, p : R h (u h, p , w h, p ) ≡ l( w h, p ) − ah (u h, p , w h, p ) = 0,

∀ w h, p ∈ Vh, p .

184

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

The subscript h is a notional measure of the local grid scale defined by the triangulation. While any given mesh can have a non-uniform spacing, we assume a quasi-uniform and shape regular refinement as h → 0, thus the elements can not become degenerate. For local interpolation errors, hκ ≡ diam(κ ) for an element κ , however given the assumption of a quasi-uniform grid, all hκ are within a constant factor of each other and thus will be replaced with h. The h subscript when applied to R h (·, ·) and ah (·, ·) denote that the residual and bilinear term are discretization specific and therefore may differ from R (·, ·) and a(·, ·). Those forms which display polynomial order dependence, for instance the f

lifting operator for DG BR2 will be denoted rh, p (·) as seen in Section 3. l( w ) is not modified by the methods considered in this paper, but it is possible with other methods. Consider a volume integral output functional given by J (u ) ≡ ( g , u ). The corresponding adjoint problem is,

ψ ∈ V : R ψ ( w , ψ) ≡ J ( w ) − a( w , ψ) = 0,

∀w ∈ V.

(2.2)

For the Poisson problem and volume output functional, J (u ), considered here, performing integration by parts gives the adjoint PDE

−ψ = g in 

(2.3)

ψ = 0 on ∂. As with the solution, u, the adjoint solution, ψ , is assumed smooth, which places restrictions on g and the domain shape. As with the primal problem, the adjoint problem has a corresponding discrete weak form ψ

ψh, p ∈ Vh, p : R h ( w h, p , ψh, p ) ≡ J ( w h, p ) − ah ( w h, p , ψh, p ) = 0,

∀ w h, p ∈ Vh, p .

2.1. A posteriori output error estimation In the DWR framework the output error, E , can be expressed by the Functional Error Representation Formula

E = J (u ) − J (u h , p )

= R h (uh, p , ψ − v h, p ),

∀ v h, p ∈ Vh, p

when the discretized residuals satisfy an extended form of primal and adjoint consistency, the definitions of which are described within a DG context in Yano’s thesis [36]. Definition 2.1. Extended Primal Global Consistency A discretization displays Extended Primal Global Consistency if the discretized primal residual satisfies

R h ( u , w ) = 0,

∀ w ∈ ⊕κ V (κ )

where V (κ ) is the restriction of V to an element

κ , and ⊕κ is the directed sum over all κ in Th .

Extended consistency is stronger than normal consistency [3] which only requires R h (u , w ) = 0, ∀ w ∈ Vh, p . Definition 2.2. Extended Dual Global Consistency A discretization displays Extended Dual Global Consistency if the discrete adjoint residual satisfies ψ

R h ( w , ψ) = 0,

∀ w ∈ ⊕κ V (κ ).

The Functional Error Representation Formula is obtained using Definitions 2.1 and 2.2. Lemma 2.1. Functional Error Representation Formula For a discretization that satisfies both Extended Primal and Dual Global Consistency, the error in the output functional E ≡ J (u ) − J (uh, p ) is given by

E = R h (uh, p , ψ − v h, p ),

∀ v h, p ∈ Vh, p .

Proof. E ≡ J (u ) − J (u h, p ) = J (u − u h, p ), and thus the proof is a consequence of Definitions 2.1 and 2.2 and orthogonality. 2 Another useful property of a discretization that satisfies Extended Primal and Dual Global Consistency is Global Residual Error Mapping.

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

185

Lemma 2.2. Global Residual-Error Mapping For all p1, p2 ∈ N, the global dual-weighted residuals are related by

R h ( v h, p1 , ψ − w h, p2 ) = ah (u − v h, p1 , ψ − w h, p2 ) ψ

= R h (u − v h, p1 , w h, p2 ), Proof. The proof is analogous to that of Lemma 2.1.

∀( v h, p1 , w h, p2 ) ∈ (Vh, p1 × Wh, p2 ). 2

2.2. Localization The localization, rκ ( v , w ), of the discrete primal residual is defined as

rκ ( v , w ) = lκ ( w ) − aκ ( v , w ), where rκ ( v , w ), aκ ( v , w ) and lκ ( w ) = ( f , w )κ are the discretization specific localized primal residual, primal bilinear and linear functional expressions respectively. Further rκ ( v , w ) is assumed to satisfy

Rh (v , w ) =



rκ ( v , w ),

∀( v , w ) ∈ (⊕κ V (κ ) × ⊕κ V (κ )).

κ ∈Th

Definition 2.3. Extended Primal Local Consistency A discretization displays Extended Primal Local Consistency if the localized residual satisfies

r κ ( u , w ) = 0,

∀ w ∈ V (κ ),

where V (κ ) enforces the boundary condition for those elements adjacent to the boundary. The localization of the adjoint residual is defined as ψ

ψ

rκ ( w , v ) = Jκ ( w ) − aκ ( w , v ), ψ

ψ

where rκ ( w , v ), aκ ( w , v ) and Jκ ( w ) = ( g , w )κ are the discretization specific localized adjoint residual, adjoint bilinear and output functional expressions respectively. Definition 2.4. Extended Dual Local Consistency A discretization displays Extended Dual Local Consistency if the localized adjoint residual satisfies ψ

rκ ( w , ψ) = 0,

∀ w ∈ V (κ ).

Provided a localization of a residual satisfies Extended Primal Local Consistency it also satisfies Local Residual-Error Mapping. Lemma 2.3. Local Residual-Error Mapping For all p1, p2 ∈ N, given a localization that satisfies Extended Local Primal Consistency, the local dual-weighted primal residual can be represented by

rκ ( v h, p1 , ψ − w h, p2 ) = aκ (u − v h, p1 , ψ − w h, p2 ), Proof. The proof is analogous to Lemma 2.1.

∀( v h, p1 , w h, p2 ) ∈ (Vh, p1 × Wh, p2 ).

2

In addition to Local Residual Error Mapping, the proofs of Extended Primal and Dual Global Consistency follow from Extended Primal and Dual Local Consistency. Corollary 2.1. Given a localization, rκ ( v , w ), of a residual, R h ( v , w ), that satisfies

Rh (v , w ) =



rκ ( v , w )

κ ∈Th

and Definitions 2.3 and 2.4, R h ( v , w ) automatically satisfies Definitions 2.1 and 2.2 because V ⊂ (⊕κ V (κ )).

186

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

Table 1 A priori convergence rates for global and local error and estimate error.

CG DG HDG(h) HDG(L)

E

E − E˜

h2p h2p h2p h2p +1

h2p  h2p 2p  h  h2p +1



ηκ

ηκ − η˜ κ

h2p +d h2p +d h2p +d h2p +1+d

h p + p +d  h p + p +d p + p  +d h  h p + p +1+d



2.3. Error estimates Given ψ ∈ V , where V is an infinite dimensional space, evaluating the true error via Lemma 2.1 is infeasible, since the analytic adjoint is rarely known. The approach considered here is to solve the adjoint equation, eq. (2.2), on the same mesh in a p  > p space. There are alternative methods for approximating the p  adjoint such as local higher-order interpolation [14–16,19,34], whereby the adjoint is interpolated using a patch-wise approach. The localization of some of these methods are compared by Richter and Wick [34], alongside the introduction of a partition of unity based approach. For linear problems, and the Poisson equation in particular, these local methods can be comparably accurate while being significantly cheaper [14]. Though solving in the p  space requires the solution of a linear system that is larger than that of the primal, for nonlinear primal problems the computational effort to solve the primal problem often dominates the total computational effort of the solution and error estimation process. For example, for a typical transonic Reynolds-Averaged Navier–Stokes (RANS) solution, the linear adjoint solve required approximately 10% of the computational time required for the nonlinear primal solve [37]. As a result higher-order adjoint solves can be used in practice [22,37]. In this work the following error estimate is considered,

E ≈ E˜ ≡ R h (uh, p , ψh, p  − v h, p ), where p 

∀ v h, p ∈ Vh, p ,

(2.4)

= p + p inc : p inc > 0. The true error has a localization

E=



ηκ ,

ηκ = rκ (uh, p , ψ − ψh, p )

κ ∈Th

which has a corresponding estimate

ηκ ≈ η˜ κ ≡ rκ (uh, p , ψh, p − ψh, p ).

(2.5)

Unlike for the global residual, orthogonality does not always hold for the localized residual, namely it is possible that rκ (u h, p , v h, p ) = 0, ∀ v h, p ∈ Vh, p . When the orthogonality does hold locally, then the error estimate can be calculated as, η˜ κ = rκ (uh, p , ψh, p ), avoiding the need to calculate ψh, p . Table 1 summarizes the key results of Sections 3 and 4 and Appendix A. All of the estimates and their localizations can be written using the abstraction given in eqs. (2.4) and (2.5); the specifics of the residual expressions and their localizations are given their respective sections. The specific new results are those for the DG and HDG estimates, the derivation of which are presented in Sections 3 and 4 respectively. The CG results are presented in Appendix A for completeness. 2.4. Additional notation On interior traces of the triangulation, the jump operator, J·K, and average operator, {·}, are defined for scalar, x, and vector, y , quantities as

JxK = (xnˆ )− + (xnˆ )+ 1

{x} = (x− + x+ ) 2

J y K = ( y · nˆ )− + ( y · nˆ )+ 1

{ y } = ( y − + y + ). 2

For boundary traces the jump and average operators take only the internal value, thus

JxK = (xnˆ )+

J y K = ( y · nˆ )+

{x} = x+

{ y } = y + ,

where nˆ + is the outward facing normal on the trace of the element. Some additional notation for inner products are defined as

( v , w )Th =



κ ∈Th

( v , w )κ ,

 v , w ∂ Th =



κ ∈Th

 v , w ∂ κ ,

 v , w Fh =



f ∈Fh

v , w f .

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

187

κ ∈ Th . Fhinternal ≡ {∂ κ + ∩ ∂ κ − } : ∀κ + , κ − ∈ Th denotes ≡ {∂ κ ∩ ∂} : κ ∈ Th denotes the set of boundary traces associated

Here ∂ Th ≡ {∂ κ : κ ∈ Th } denotes the union of the boundaries of all the set of internal traces associated with Th , boundary

boundary Fh

with Th and Fh ≡ Fhinternal ∪ Fh denotes the set of all traces of Th . The latter two outer products are different in that the first of the two visits internal traces twice, once from either side, whereas the second visits each trace once which implies that

1 2

u , v ∂ Th \∂ ≡ u , v Fh \∂ .

Finally, to simplify notation when taking limits, A (h)  B (h) is defined as A (h) ≤ c B (h) in the limit of h → 0 for 0 < c < ∞. 3. Discontinuous Galerkin In this paper Bassi and Rebay’s second (BR2) Discontinuous Galerkin (DG) discretization is considered [13]. The convergence of global and local DWR estimates for BR2 has already been analyzed by Yano [36], who showed that the p-dependency of the lifting operator can limit the convergence rate of the global error estimate. To resolve this a mixed formulation for the error estimate is proposed. 3.1. Discretization For the DG discretization, the approximation space Vh, p is made up of piecewise discontinuous polynomials of order p,

Vh, p ≡ { v ∈ L 2 () : v |κ ∈ P p (κ ),

∀κ ∈ Th },

and the bilinear term is given by

ah, p ( v , w ) = (∇ v , ∇ w )Th − {∇ v }, J w K Fh − J v K, {∇ w } Fh f

− ξ { rh, p (J v K)}, J w K Fh . The p subscript denotes that the bilinear term is p dependent, this means the residual operator R h, p ( v , w ) is also p dependent. For stability ξ must be set to greater than the number of traces of the element where the residual is evaluated, i.e. 3 for 2D simplices [18]. The lifting operator is defined on trace f ∈ Fh by the expression [3],

rhf, p ( v ) ∈ [Vh, p ]d : ( rhf, p ( v ), g )Th = − v , { g } f ,

∀ g ∈ [Vh, p ]d

where v ∈ [ L 1 ( f )]d . The p-dependence of the DG residual operator is contained in the expression for the lifting operator

rhf, p ( v ) [36], namely

rhf,q ( v h,q ) = rhf, p ( v h,q ),

v h,q ∈ Vh,q ,

∀q, p ∈ N : q < p .

The p dependence suggests two similar estimates, one where the residual is evaluated in p, E˜ 1 , and one where the residual is evaluated in p  , E˜ 2 [36].

E˜1 = R h, p (uh, p , ψh, p  ),

E˜2 = R h, p  (uh, p , ψh, p  ).

The sub-optimal performance is a result of not accounting for the error in the computation of r f [20]. An alternative strategy is to append the lifting operator definition in the weak form to the residual such that the bilinear term for the BR2 discretization is written as

ah ( v , r F ; w , s F ) = (∇ v , ∇ w )Th

− {∇ v }, J w K Fh − J v K, {∇ w } Fh − ξ { r f }, J w K Fh    − ( r f , s f )Th + J v K, { s f } f , f ∈Fh



d where r F and s F ∈ ShF, p and ShF, p ≡ r f is the element in r F corresponding to face f . The infinite f ∈Fh [Vh, p ] , and  F d

F ) and W = ( w , s F ) tuples dimensional equivalent is denoted S ≡ r F ),  = (ψ, σ f ∈Fh [V ] . For succinctness U = (u ,

F is the adjoint lifting operator, and W h, p ≡ (Vh, p × ShF, p ). Thus R h (Uh, p , W) ≡ are specified as arguments, where σ F F R h (u h, p , rh, p ; w , s ) and the discretization may be written as Uh, p ∈ W h, p : R h (Uh, p , Wh, p ) = 0,

∀Wh, p ∈ W h, p .

188

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

F can be post processed. Given this residual is used only for the estimation process, the explicit calculation of r F and σ The exact r f is defined by substituting in the strong form solution u ∈ C 1 (), thus Ju K = 0. Then setting w = 0 shows that

r f = 0 for the exact solution. The adjoint residual for the DG formulation can be written using the tuples as ψ

R h (W, V) = J ( w ) − ah (W, V).

Substituting in ψ ∈ C 1 (), JψK = 0, then setting w = 0 shows that σ

f

= 0.

3.2. Localization The DG residual is simply localized,

rκ (V, W) = R h (V, W|κ ), where W|κ is defined to be 0 outside of elementally restricted.

κ . Unlike CG, DG displays extended local consistency when the test functions are

Lemma 3.1. DG Extended Local Consistency

rκ (U, W) = 0,

∀W ∈ W(κ )

ψ

∀W ∈ W(κ ),

rκ (W, ) = 0,

where W(κ ) is the elemental restrictions of W . Proof. For the exact solution U, Ju K = 0 and r f = 0 leaving

rκ (U, W) = ( f , w )κ − (∇ u , ∇ w )κ + ∇ u · nˆ , w ∂ κ

= ( f + u , w )κ = 0, The proof for the dual is analogous.

∀ w ∈ V (κ ).

2

Lemma 3.2. DG Extended Global Consistency

R h (U, W) = 0,

∀W ∈ ⊕κ W(κ )

ψ

∀W ∈ ⊕κ W(κ ).

R h (W, ) = 0,

Proof. The proof follows directly from Lemma 3.1.

2

3.3. Bilinear Error Bounds Assuming sufficient smoothness of the solution and coercivity of the PDE, error bounds [3,18] can be combined with a d

smooth compact function, φ , satisfying φ = 1, ∀x ∈ κ and φ L 2 ()  h 2 . The convergence of the weighted error |φ( v − v h, p )| H m () can then be used to establish rates of convergence for the localized error, and for a given v ∈ C 0 , bounds on interface jumps, J v h, p K L 2 ( f ) can be used d

| v − v h, p | H m (κ )  hn+1−m+ 2 | v | H n+1 () , | v − v h, p | H m ( f )  h J v h, p KL 2 ( f )  h

n+ 12 −m+ d2

− 12

| v | H n+1 () ,

m ∈ (0, 1)

(3.1)

m ∈ (0, 1)

 v − v h , p  L 2 (κ f ) ,

where κ f is the set of κ adjacent to f . For the case of sufficient regularity, n is equal to p, but were v to be insufficiently smooth then the maximum n < p, this would reduce the maximum convergence rate, thus were v ∈ / H 3 (), the maximum f

f

h , p . n + 1 would be 2, because | v | H 3 () ≮ ∞ [1]. In addition, error bounds are required for rh, p and σ Lemma 3.3. DG Auxiliary Variable Bounds 1

 rh, p ( v ) L 2 (κ f )  h− 2 J v h, p KL 2 ( f ) . f

Proof. The proof can be found in Lemma 2 of [18].

2

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

189

Lemma 3.4. DG Local Bilinear Error Bound

aκ (U − Uh, p1 ,  −  h, p2 ) ≡ ah (U − Uh, p1 , ( −  h, p2 )|κ )

 h p1+ p2+d |u | H p1+1 () |ψ| H p2+1 () , where (Uh, p1 ,  h, p2 ) are finite element solutions of polynomial order p1, p2 ∈ N. Proof.

aκ (U − Uh, p1 ,  −  h, p2 )

= (∇(u − uh, p1 ), ∇(ψ − ψh, p2 ))κ + {∇(u − uh, p1 )} · nˆ , ψ − ψh, p2 ∂ κ  

 

I

+ Ju − uh, p1 K, α ∇(ψ − ψh, p2 ) ∂ κ −  

III

+



f

f

( rh, p1 , σ h, p2 )κ −

f ∈∂ κ











II f ξ { rh, p1 } · nˆ , ψ

f ∈∂ κ



− ψh, p2 f





IV f

Ju − uh, p1 K, α σ h, p2 f ,

f ∈∂ κ





V



VI

where α = 12 ∀ f ∈ Fhinternal and individually.

α = 1 ∀ f ∈ Fhboundary . In order to bound the whole, it suffices to bound the components

I : (∇(u − u h, p1 ), ∇(ψ − ψh, p2 ))κ

≤ ∇(u − uh, p1 ) L 2 (κ ) ∇(ψ − ψh, p2 ) L 2 (κ ) II : {∇(u − u h, p1 )}, ψ − ψh, p2 ∂ κ

≤ {∇(u − uh, p1 )} · nˆ L 2 (∂ κ ) ψ − ψh, p2 L 2 (∂ κ ) III : Ju − u h, p1 K, α ∇(ψ − ψh, p2 ) ∂ κ

 Ju − uh, p1 KL 2 (∂ κ ) ∇(ψ − ψh, p2 ) L 2 (∂ κ )  f IV : ξ { rh, p1 } · nˆ , ψ − ψh, p2 f f ∈∂ κ





f

{ rh, p1 }L 2 ( f ) ψ − ψh, p2 L 2 ( f )

f ∈∂ κ





h−2 u − u h, p1  L 2 (κ f ) ψ − ψh, p2  L 2 (κ )

f ∈∂ κ

V:



f

f

( rh, p1 , σ h, p2 )κ

f ∈∂ κ

≤



f

f

 rh, p1 L 2 (κ ) σ h, p2 L 2 (κ )

f ∈∂ κ

VI :



f

Ju − uh, p1 K, α σ h, p2 f

f ∈∂ κ





f

Ju − uh, p1 KL 2 ( f ) σ h, p2 L 2 ( f )

f ∈∂ κ





h, p2 L 2 (κ ) h−1 u − u h, p1  L 2 (κ f ) σ f

f ∈∂ κ

All of these are bound by h p1+ p2+d |u | H p1+1 () |ψ| H p2+1 () thus,

aκ (U − Uh, p1 ,  −  h, p2 )  h p1+ p2+d |u | H p1+1 () |ψ| H p2+1 () .

2

Using Corollary 2.1 the Global Bilinear Error Bound is given by Lemma 3.5.

190

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

Lemma 3.5. DG Global Bilinear Form Error Bound

ah (U − Uh, p1 ,  −  h, p2 )  h p1+ p2 |u | H p1+1 () |ψ| H p2+1 () . Proof.

ah (U − Uh, p1 ,  −  h, p2 ) =



aκ (U − Uh, p1 ,  −  h, p2 )

κ ∈Th





h p1+ p2+d |u | H p1+1 () |ψ| H p2+1 ()

κ ∈Th

 h p1+ p2 |u | H p1+1 () |ψ| H p2+1 () .

2

3.4. Global Error Bound In this section Lemma 3.5 is used to derive asymptotic convergence rates for E and E − E˜ for the DG discretization. Theorem 3.1. DG Global Error Bound

|E |  h2p |u | H p+1 () |ψ| H p+1 () . Proof.

|E | = | R h (Uh, p ,  − Wh, p )|, ∀Wh, p ∈ W h, p = | R h (Uh, p ,  −  h, p )| = |ah (U − Uh, p ,  −  h, p )| via Lemma 2.2  h2p |u | H p+1 () |ψ| H p+1 () via Lemma 3.5.

2

Theorem 3.2. DG Global Estimate Error Bound  |E − E˜|  h2p |u | H p +1 () |ψ| H p +1 () .

Proof.

|E − E˜| = | R h (Uh, p ,  − Wh, p ) − R h (Uh, p ,  h, p  − Wh, p )|, ∀Wh, p ∈ W h, p = | R h (Uh, p ,  −  h, p  )| ψ

= | R h (U − Uh, p ,  h, p  )| via Lemma 2.2 ψ

= | R h (U − Wh, p  ,  h, p  )|, ∀Wh, p  ∈ W h, p 

via orthogonality

= |ah (U − Wh, p  ,  −  h, p  )|, ∀Wh, p  ∈ W h, p  h

2p 

|u | H p +1 () |ψ| H p +1 () via Lemma 3.5.

via Lemma 2.2

2

3.5. Local Error Bound In this section Lemma 3.4 is used to derive asymptotic convergence rates for ηκ and ηκ − η˜ κ for the DG discretization. In common with the global results, all that is necessary to bound the relevant formulae are Lemmas 2.3 and 3.4. Unlike for CG (see Theorem A.3), the DG localized residual is the same as the global residual, aside from the local restriction of test functions, as such rκ (Uh, p , Wh, p ) = 0, ∀Wh, p ∈ W h, p . Theorem 3.3. DG Local Error Bound

|ηκ |  h2p +d |u | H p+1 () |ψ| H p+1 () .

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

191

Proof.

|ηκ | = |rκ (Uh, p ,  − Wh, p )|, ∀Wh, p ∈ W h, p (κ ) = |rκ (Uh, p ,  −  h, p )| = |aκ (U − Uh, p ,  −  h, p )| via Lemma 2.3  h2p +d |u | H p+1 () |ψ| H p+1 () via Lemma 3.4.

2

Theorem 3.4. DG Local Estimate Error Bound 

|ηκ − η˜ κ |  h p + p +d |u | H p+1 () |ψh, p  | H p+1 () . Proof.

|ηκ − η˜ κ | = |rκ (Uh, p ,  − Wh, p ) − rκ (Uh, p ,  h, p  − Wh, p )|, ∀Wh, p ∈ W h, p (κ ) = |rκ (Uh, p ,  −  h, p  )| = |aκ (U − Uh, p ,  −  h, p  )| via Lemma 2.3 

 h p + p +d |u | H p+1 () |ψ| H p +1 () via Lemma 3.4.

2

4. Hybridizable Discontinuous Galerkin The Hybridizable Discontinuous Galerkin, or HDG, method is a mixed method, in which the first derivative on the scalar variable is introduced as an auxiliary variable, q ,

q − ∇ u = 0 in 

−∇ · q = f

in 

u = 0 on ∂. The particular HDG method here is derived from the Local Discontinuous Galerkin (LDG) method and is sometimes referred to as the LDG-H method [24,25]. 4.1. Discretization In addition to replacing the solution gradient with an additional variable HDG replaces the trace of the scalar on the ˆ Thus uˆ = u and ψˆ = ψ on the element traces, for the exact solutions u element boundaries with an additional variable u. ˆ and W ≡ ( w , s, w

, uˆ ),  ≡ (ψ, ϕ , ψ) ˆ ), thus R h (Uh, p , W) ≡ R h (uh, p , q h, p , uˆ h, p ; w , s, w ˆ ). and ψ . For ease of notation U ≡ (u , q The discretization is defined by

Vh, p ≡ { v ∈ L 2 () : v |κ ∈ P p (κ ),

∀κ ∈ Th }

Qh, p ≡ { q ∈ [ L ()] : q |κ ∈ [P (κ )]d , 2

d

p

∀κ ∈ Th }

Mh, p ≡ {m ∈ [ L 2 (Fh )]d−1 : m| f ∈ [P p ( f )]d−1 ,

∀ f ∈ Fh },

with bilinear term

ˆ ) = ( t , ∇ w )Th +  F v · nˆ , w ∂ Th ah ( v , t , vˆ ; w , s, w + ( t , s)Th + ( v , ∇ · s)Th −  vˆ , s · nˆ ∂ Th ˆ Fh − J F v K, w F v = − t + (τ nˆ )( v − vˆ ).

τ is a stabilization parameter that must be greater than 0 for stability [26]. In this paper two methods of choosing the length scale are analyzed. The first HDG(h) uses the local length scale of the mesh, and the second HDG(L) uses a globally specified length scale L.

192

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

The adjoint equation can also be written in mixed form

ϕ − ∇ψ = 0 in  −∇ · ϕ = g in  ψ = 0 on ∂. The spaces Vh, p , Qh, p and Mh, p are discrete counterparts to the infinite dimensional spaces V , Q and M,

V ≡ L 2 ()

Q ≡ [ L 2 ()]d M ≡ {m ∈ [ L 2 ( f )]d−1 ,

∀ f ∈ Fh }.

The collection of elemental spaces is denoted

W ≡ (⊕κ (V (κ ) × Q(κ )) × M) with corresponding elemental restriction

W(κ ) ≡ (V (κ ) × Q(κ ) × M(∂ κ )). 4.2. Localization Consider the following localization

rκ (U, W) = lκ ( w ) − aκ (U, W) aκ (U, W) = ( q, ∇ w )κ +  F u · nˆ , w ∂ κ

+ ( q, s)κ + (u , ∇ · s)κ − uˆ , s · nˆ ∂ κ 1

ˆ ∂ κ , − J F u K, w 2

which closely follows the global form except for apportioning half the jump term to the two elements adjacent to each trace. Using this apportionment, the sum of the local residuals can be readily shown to equal the global residual. Lemma 4.1. HDG Extended Local Consistency

rκ (U, W) = 0,

∀W ∈ W(κ )

ψ

∀W ∈ W(κ ).

rκ (W, ) = 0,

Proof. F u ≡ − q + (τ nˆ )(u − uˆ ) occurs only in trace integrals thus F u ≡ − q,

ˆ ) = ( f , w )κ − ( q, ∇ w )κ −  F u · nˆ , w ∂ κ rκ (u , q , uˆ ; w , s, w − ( q, s)κ − (u , ∇ · s)κ + uˆ , s · nˆ ∂ κ 1

ˆ ∂ κ + J F u K, w 2

= ( f + ∇ · q , w )κ − ( F u + q ) · nˆ , w ∂ κ − ( q − ∇ u , s)κ − u − uˆ , s · nˆ ∂ κ 1

ˆ ∂ κ = 0, + J F u K, w 2

∀W ∈ W(κ ).

The proof for the dual proceeds analogously. Extended Global Consistency follows naturally from Lemma 4.1. Lemma 4.2. HDG Extended Global Consistency

R h (U, W) = 0,

∀W ∈ W

ψ

∀W ∈ W.

R h (W, ) = 0,

Proof. The proof is a direct consequence of Lemma 4.1

2

2

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

193

p

4.3. The h projector p

p

p

To prove the convergence rates for the output error, the h ( q, u ) ≡ (s q ,  w u ) projector as introduced by Cockburn et al. [25] is used. The projector is defined for use with the primal problem, but it can also be applied to the adjoint. The projector gives bounds on the local and the global solution errors of the HDG discretization in terms of semi-norms on the exact quantities, whilst accounting for τ dependencies. p

Definition 4.1. The h Projector

∀ s ∈ [P p −1 (κ )]d

(s t , s)κ = ( t , s)κ , p

∀ w ∈ P p −1 (κ )

p

( w v , w )κ = ( v , w )κ , ˆ ∂ κ =  t · nˆ + τ v , w ˆ ∂ κ , s t · nˆ + τ  w v , w p

p

ˆ ∈ P p (∂ κ ). ∀w

p

The domain of h is Qh, p × Vh, p . p

p

The h projector can be related to the P M projector, defined as the L 2 projection in p onto M( f )

s t · nˆ + τ ( w v − P v ). P M ( t · nˆ ) =  M p

p

p

p

p

Theorem 4.1. h Projector Error Convergence rates [25, Theorem 2.1] p1 s t − t  H m (κ )  h p1+1−m | t | H p1+1 (κ ) + h p1+1−m τ | v | H p1+1 (κ )

 w v − v  H m (κ )  h p1+1−m | v | H p1+1 (κ ) + p1

h p1+1−m

τ

|∇ · t | H p1 (κ ) .

Given t ≡ ∇ v for a mixed formulation, combined with the regularity and smoothness assumptions from Section 2 and a smooth compact extension (see Section 3), d

p1 s t − t  H m (κ )  h p1+1−m+ 2 (| v | H p1+2 () + τ | v | H p1+1 () ) d

 w v − v  H m (κ )  h p1+1−m+ 2 (| v | H p1+1 () + p1

1

τ

| v | H p1+2 () ).

Using these expressions, norms for the local finite element error can be derived, assuming that hτ  1 p p  t − th, p  H m (κ ) ≤  t − s t  H m (κ ) + s t − th, p  H m (κ )

  t − s t  H m (κ ) p

d

 h p +1−m+ 2 (| v | H p+2 () + τ | v | H p+1 () ), p

p

 v − v h, p  H m (κ ) ≤  v −  w u  H m (κ ) +  w v − v h, p  H m (κ ) p p   v −  w v  H m (κ ) + h t − s t  H m (κ ) d

 h p +1−m+ 2 (| v | H p+1 () +

1

τ

| v | H p+2 () ).

The final piece is to bound the errors associated with the element traces, p p J t − F v KL 2 ( f ) ≤ J t − P M t KL 2 ( f ) + J P M t − F v KL 2 ( f ) 1

p p  h− 2 ( t − P M t L 2 (κ f ) + s t − t L 2 (κ f ) ) 1

 h− 2  t − th, p L 2 (κ f ) , p

p

 v − vˆ h, p L 2 (∂ κ ) ≤  v − P M v L 2 (∂ κ ) +  P M v − vˆ h, p L 2 (∂ κ ) 1

  v − P M v L 2 (∂ κ ) + h− 2  P M v − vˆ h, p h(κ ) p

p

1

p p   v − P M v L 2 (∂ κ ) + h 2 s t − t L 2 (κ )

 h p+ h

1+d 2

d p + 1+ 2

(| v | H p+1 () + h(| v | H p+2 () + τ | v | H p+1 () )) | v | H p+1 () .

194

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

4.4. Bilinear Error Bounds Lemma 4.3. HDG Local Bilinear Error Bound

aκ (U − Uh, p1 ,  −  h, p2 )

 1  h p1+ p2+1+d |u | H p1+1 () + |u | H p1+2 () )(|ψ| H p2+2 () + τ |ψ| H p2+1 () )

τ

+ (|u | H p1+2 () + τ |u | H p1+1 () )(|ψ| H p2+1 () + 

1

τ



1

τ

 |ψ| H p2+2 () )

+ τ h p1+ p2+1+d ,

where (Uh, p1 ,  h, p2 ) are finite element solutions of order p1 and p2 and τ is the stabilization parameter, assumed to satisfy hτ  1. Proof.

aκ (U − Uh, p1 ,  −  h, p2 )

= ( q − q h, p1 , ∇(ψ − ψh, p2 ))κ − ( q − F u ) · nˆ , ψ − ψh, p2 ∂ κ  

 

I

II

+ ( q − q h, p1 , ϕ − ϕ h, p2 )κ + (u − uh, p1 , ∇ · (ϕ − ϕ h, p2 ))κ  

 

III

IV

1

− u − uˆ h, p1 , (ϕ − ϕ h, p2 ) · nˆ ∂ κ − J q − F u K, ψ − ψˆ h, p2 ∂ κ \∂ .  

2 

V

VI

To bound the whole it suffices to bound each of the components individually using the triangle inequality

I : ( q − q h, p1 , ∇(ψ − ψh, p2 ))κ

≤  q − q h, p1 L 2 (κ ) ψ − ψh, p2  H 1 (κ )  h p1+ p2+1+d (|u | H p1+2 () + τ |u | H p1+1 () )(|ψ| H p2+1 () +

1

τ

|ψ| H p2+2 () )

II : ( q − F u ) · nˆ , ψ − ψh, p2 ∂ κ

≤ ( q − F u ) · nˆ L 2 (∂ κ ) ψ − ψh, p2 L 2 (∂ κ )  h p1+ p2+1+d (|u | H p1+2 () + τ |u | H p1+1 () )(|ψ| H p2+1 () +

1

τ

|ψ| H p2+2 () )

III : ( q − q h, p1 , ϕ − ϕ h, p2 )κ

≤  q − q h, p1 L 2 (κ ) ϕ − ϕ h, p2 L 2 (κ )  h p1+ p2+2+d (|u | H p1+2 () + τ |u | H p1+1 () )(|ψ| H p2+2 () + τ |ψ| H p2+1 () ) IV : (u − u h, p1 , ∇ · (ϕ − ϕ h, p2 ))κ

≤ u − uh, p1 L 2 (κ ) ϕ − ϕ h, p2  H 1 (κ )  h p1+ p2+1+d (|u | H p1+1 () +

1

τ

|u | H p1+2 () )(|ψ| H p2+2 () + τ |ψ| H p2+1 () )

V : u − uˆ h, p1 , (ϕ − ϕ h, p2 ) · nˆ ∂ κ

≤ u − uˆ h, p1 L 2 (∂ κ ) (ϕ − ϕ h, p2 ) · nˆ L 2 (∂ κ )    h p1+ p2+1+d |u | H p1+1 () |ψ| H p2+2 () + τ |ψ| H p2+1 () VI : J q − F u K, ψ − ψˆ h, p2 ∂ κ \∂

≤ J q − F u KL 2 (∂ κ ) ψ − ψˆ h, p2 L 2 (∂ κ )  h p1+ p2+1+d (|u | H p1+2 () | + τ |u | H p1+1 () )|ψ| H p2+1 () .

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

195

Combining these together,

aκ (U − Uh, p1 ,  −  h, p2 )

 1  h p1+ p2+1+d (|u | H p1+1 () + |u | H p1+2 () )(|ψ| H p2+2 () + τ |ψ| H p2+1 () ) τ  1 + (|u | H p1+2 () + τ |u | H p1+1 () )(|ψ| H p2+1 () + |ψ| H p2+2 () ) τ 1  p1+ p2+1+d  +τ h . 2

τ

Using Lemma 4.3 the Global Bilinear Error Bound can be proven. Lemma 4.4. HDG Global Bilinear Error Bound

ah (U − Uh, p1 ,  −  h, p2 )

 1  h p1+ p2+1 (|u | H p1+1 () + |u | H p1+2 () )(|ψ| H p2+2 () + τ |ψ| H p2+1 () ) τ  1 + (|u | H p1+2 () + τ |u | H p1+1 () )(|ψ| H p2+1 () + |ψ| H p2+2 () ) τ 1   + τ h p1+ p2+1 .

τ

Proof. The proof follows directly from the local result.

2

4.5. Global Error Bound In this section Lemma 4.4 is used to derive asymptotic convergence rates for E and E − E˜ for the HDG discretization. Theorem 4.2. HDG Global Error Bound

 1 |E |  h2p +1 |u | H p+1 () + |u | H p+2 () )(|ψ| H p+2 () + τ |ψ| H p+1 () ) τ  1 + (|u | H p+2 () + τ |u | H p+1 () )(|ψ| H p+1 () + |ψ| H p+2 () ) τ 1  2p +1  +τ h ,

τ

where τ is the stabilization parameter, assumed to satisfy hτ  1. Proof.

|E | = | R h (Uh, p ,  − Wh, p )|, ∀Wh, p ∈ W h, p = | R h (Uh, p ,  −  h, p )| = |ah (U − Uh, p ,  −  h, p )| via Lemma 2.2  1  h2p +1 |u | H p+1 () + |u | H p+2 () )(|ψ| H p+2 () + τ |ψ| H p+1 () ) τ  1 + (|u | H p+2 () + τ |u | H p+1 () )(|ψ| H p+1 () + |ψ| H p+2 () ) τ 1   + τ h2p +1 via Lemma 4.4. 2

τ

Theorem 4.3. HDG Global Estimate Error Bound

1

 |E − E˜|  h2p +1 (|u | H p +1 () + |u | H p +2 () )(|ψ| H p +2 () + τ |ψ| H p +1 () ) τ 1    h2p +1 +τ ,

τ

where τ is the stabilization parameter, assumed to satisfy hτ  1.

196

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

Proof.

|E − E˜| = | R h (Uh, p ,  − Wh, p ) − R (Uh, p ,  h, p  − Wh, p )|, ∀Wh, p ∈ W h, p = | R h (Uh, p ,  −  h, p  )| ψ

= | R h (U − Uh, p ,  h, p  )| via Lemma 2.2 ψ

= | R h (U − Wh, p  ,  h, p  )|, ∀Wh, p  ∈ W h, p 

via orthogonality

= |ah (U − Wh, p  ,  −  h, p  )|, ∀Wh, p  ∈ W h, p  via Lemma 2.2  1   h2p +1 |u | H p +1 () + |u | H p +2 () )(|ψ| H p +2 () + τ |ψ| H p +1 () )

τ

+ (|u | H p +2 () + τ |u | H p +1 () )(|ψ| H p +1 () + 

1

τ



+τ h

2p  +1

via Lemma 4.4.

1

τ

 |ψ| H p +2 () )

2

4.6. Local Error Bound In this section Lemma 4.3 is used to derive asymptotic convergence rates for ηκ and ηκ − η˜ κ for the HDG discretization. In common with the global results, all that is necessary to bound the relevant formulae are Lemmas 2.3 and 3.4. The HDG localized residual is the same as the global residual, aside from the local restriction of test functions and the factor of 12 on the jump condition, as such rκ (Uh, p , Wh, p ) = 0, ∀Wh, p ∈ W h, p (κ ). Theorem 4.4. HDG Local Error Bound

 1 |ηκ |  h2p +1+d |u | H p+1 () + |u | H p+2 () )(|ψ| H p+2 () + τ |ψ| H p+1 () )

τ

+ (|u | H p+2 () + τ |u | H p+1 () )(|ψ| H p+1 () + 

1

τ



+τ h

2p +1+d

1

τ

 |ψ| H p+2 () )

,

where τ is the stabilization parameter, assumed to satisfy hτ  1. Proof.

|ηκ | = |rκ (Uh, p ,  − Wh, p )|, ∀Wh, p ∈ W h, p (κ ) = |rκ (Uh, p ,  −  h, p )| = |aκ (U − Uh, p ,  −  h, p )| via Lemma 2.3  1  h2p +1+d |u | H p+1 () + |u | H p+2 () )(|ψ| H p+2 () + τ |ψ| H p+1 () )

τ

+ (|u | H p+2 () + τ |u | H p+1 () )(|ψ| H p+1 () + 

1

τ



+ τ h2p +1+d via Lemma 4.3.

1

τ

 |ψ| H p+2 () )

2

Theorem 4.5. HDG Local Estimate Error Bound

 1  |ηκ − η˜ κ |  h p + p +1+d |u | H p+1 () + |u | H p+2 () )(|ψ| H p +2 () + τ |ψ| H p +1 () )

τ

+ (|u | H p+2 () + τ |u | H p+1 () )(|ψ| H p +1 () + 

1

τ



+τ h

p + p  +1+d

,

where τ is the stabilization parameter, assumed to satisfy hτ  1.

1

τ

 |ψ| H p +2 () )

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

197

Fig. 1. 2D Data.

Proof.

|ηκ − η˜ κ | = |rκ (Uh, p ,  − Wh, p ) − rκ (Uh, p ,  h, p  − Wh, p )|, ∀Wh, p ∈ W h, p (κ ) = |rκ (Uh, p ,  −  h, p  )| = |aκ (U − Uh, p ,  −  h, p  )| via Lemma 2.3  1   h p + p +1+d |u | H p+1 () + |u | H p+2 () )(|ψ| H p +2 () + τ |ψ| H p +1 () ) τ  1 + (|u | H p+2 () + τ |u | H p+1 () )(|ψ| H p +1 () + |ψ| H p +2 () ) τ 1    + τ h p + p +1+d via Lemma 4.3. 2

τ

5. Numerical results To demonstrate the theoretical convergence rates, exact primal and adjoint solutions with sufficient smoothness are required. To do this, u , ψ ∈ C ∞ () are chosen and the required forcing functions f and g are then determined. The primal and adjoint equations are

u = (1 − x)(1 − y )(1 − e 3xy ),

in (0, 1)2

ψ = sin(π x) sin(π y ),

in (0, 1)2

The forcing functions, f and g, required to give u and ψ are shown in Figs. 1(a) and 1(b). Averages of absolute values of the local quantities are used to measure the convergence of the local error estimates, thus here



|ηκ¯ | =

κ ∈Th |ηκ |

N



,

|ηκ¯ − η˜ κ¯ | =

κ ∈Th |ηκ − η˜ κ |

N

.

Table 1 shows a summary of the expected rates of convergence for all of the error estimates and schemes analyzed in Appendix A, Sections 3 and 4, and is thus the reference against which the numerical results are compared. p inc = 2 is used here to demonstrate that the results hold for all p inc given sufficient smoothness of u and ψ , in practice a p inc = 1 would be used for practical estimation. Global Estimation — Figs. 2(a) and 2(b) show the global error and error in the error estimate for CG, DG, HDG(h) and 1 HDG(L). The global length scale was taken to be 10 for HDG(L). All of the estimates considered are converging exactly as expected from Table 1. HDG(h) has approximately an order of magnitude less error and error in the estimate than CG and DG for p = 1, 2. CG and DG are approximately equal, with DG having slightly less error than CG for a given mesh. HDG(L) has approximately half the error of HDG(h) and the ratio increases as the grid is refined with HDG(L) performing better relatively because of the additional O (h) convergence of HDG(L) compared to HDG(h). Local Estimation — Figs. 3(a) and 3(b) show the local error and error in the error estimate, with expected rates of convergence for |ηκ¯ | and |ηκ¯ − η˜ κ¯ | given in Table 1. |ηκ¯ − η˜ κ¯ | are all converging an average p inc /2 faster than expected. HDG(h) has approximately an order of magnitude less error than CG and DG in |ηκ |, with HDG(L) having about half an order of magnitude less error than HDG(h). HDG(L) has increasingly better error than HDG(h) with mesh refinement

198

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

Fig. 2. Global error, |E |, and estimate error, |E − E˜ | | Solid: p = 1, Dashed: p = 2 | : CG, : DG, : HDG(h), ♦: HDG(L).

Fig. 3. [CG, DG and HDG Local Error and Estimates – 2D] Local error, |ηκ¯ |, and estimate error, |ηκ¯ − η˜ κ¯ | | Solid: p = 1, Dashed: p = 2 | : CG, : DG, : HDG(h), ♦: HDG(L).

because of the additional O (h) convergence rate. For |ηκ¯ − η˜ κ¯ |, HDG(L) also has about half an order of magnitude less error than HDG(h), with HDG(L) becoming increasingly smaller error than HDG(h) due to the additional O (h) convergence rate. The choice of global length scale, L, for HDG(L) is relatively simple for a linear one dimensional case, but for multidimensional cases or a PDE with anisotropic features (for example a boundary layer), it is difficult to know a priori what length scale to choose. Choosing a length scale too small ensures the asymptotic convergence rate, but at the cost of the initial error being larger, and vice versa, thus the correct choice of L is difficult a priori. Fig. 3(b) shows that instead choosing the reference length scale of the order of the grid spacing is a viable alternative that would require no a priori knowledge of the important underlying features. In summary, all of the schemes converge as predicted by the a priori analysis, or slightly better. All of the data used for generating these plots are presented in [20]. Acknowledgements This research was supported through: a Research Agreement with Saudi Aramco, a Founding Member of the MIT Energy Initiative with technical monitors Dr. Ali Dogru and Dr. Nick Burgess; and a Research Agreement with The Boeing Company with technical monitor Dr. Mori Mani. Appendix A. Continuous Galerkin The DWR framework was originally developed for the CG method [14], and combining DWR with classical theory for finite element methods [17], the asymptotic convergence rate of the output functional and its estimates can be bound as shown in this section.

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

199

A.1. Discretization For the CG discretization, the approximation space Vh, p is made up of piecewise continuous polynomials of order p denoted P p (κ ),

Vh, p ≡ { v ∈ C 0 () ∩ L 2 () : v |κ ∈ P p (κ ), ∀κ ∈ Th : v |∂ = 0}. The bilinear form is

ah ( v , w ) = (∇ v , ∇ w )Th . A.2. Localization The weighted strong form of the residual is used for the localization, following the method established by Babuška and Rheinboldt in their paper on residual based methods [7] and used by Becker and Rannacher [14].

R h ( v , w ) = ( f , w )Th − (∇ v , ∇ w )Th

= ( f , w )Th + ( v , w )Th − ∇ v · nˆ , w ∂ Th 1

= ( f +  v , w )Th − J∇ v K, w ∂ Th .

(A.1)

2

From eq. (A.1), aκ ( v , w ) and lκ ( w ), and thus rκ ( v , w ), are defined as

lκ ( w ) = ( f , w )κ aκ ( v , w ) = −( v , w )κ +

1 2

J∇ v K, w ∂ κ .

Lemma A.1. CG Extended Local Consistency

rκ (u , w ) = ( f + u , w )κ − ψ

rκ ( w , ψ) = ( g + ψ, w )κ −

1 2 1 2

J∇ u K, w ∂ κ = 0, J∇ψK, w ∂ κ = 0,

∀ w ∈ V (κ ) ∀ w ∈ V (κ ),

where V (κ ) is the restriction of V to an element κ . Proof. The volume integrals are zero since the strong form solutions satisfy eqs. (2.1) and (2.3). The trace integrals are zero because of the assumption that u and ψ are in at least C 1 (). 2 Given Extended Local Consistency is satisfied, Extended Global Consistency is also satisfied via Corollary 2.1. A.3. Bilinear Error Bounds For the following bounds, we assume shape regularity and also sufficient smoothness in u and ψ such that bounds of the form in eq. (3.1) may be used. Another useful result [35] for the error analysis is

(u − uh, p ) L 2 (κ )  h−1 |u − uh, p | H 1 (κ ) . Lemma A.2. CG Local Bilinear Error Bound

aκ (u − u h, p1 , ψ − ψh, p2 )  h p1+ p2+d |u | H p1+1 () |ψ| H p2+1 () , where (u h, p1 , ψh, p2 ) ∈ (Vh, p1 × Vh, p2 ) are finite element solutions of polynomial order p1, p2 ∈ N. Proof.

aκ (u − u h, p1 , ψ − ψh, p2 ) = −((u − u h, p1 ), ψ − ψh, p2 )κ +

1 2

J∇(u − uh, p1 )K, ψ − ψh, p2 ∂ κ 1

≤ (u − uh, p1 ) L 2 (κ ) ψ − ψh, p2 L 2 (κ ) + J∇(u − uh, p1 )K L 2 (∂ κ ) ψ − ψh, p2 L 2 (∂ κ ) 2

 h p1+ p2+d |u | H p1+1 () |ψ| H p2+1 () .

2

200

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

Lemma A.3. CG Global Bilinear Error Bound

ah (u − u h, p1 , ψ − ψh, p2 )  h p1+ p2 |u | H p1+1 () |ψ| H p2+1 () . Proof.

ah (u − u h, p1 , ψ − ψh, p2 ) =



aκ (u − u h, p1 , ψ − ψh, p2 )

κ ∈Th





h p1+ p2+d |u | H p1+1 () |ψ| H p2+1 ()

κ ∈Th

 h p1+ p2 |u | H p1+1 () |ψ| H p2+1 () .

2

A.4. Global Error Bounds In this section Lemma A.3 is used to derive asymptotic convergence rates for E and E − E˜ for the CG discretization. Theorem A.1. CG Global Error Bound Let u h, p ∈ Vh, p be the finite element solution of order p. Then the global error satisfies

|E |  h2p |u | H p+1 () |ψ| H p+1 () . Proof.

|E | = | R h (uh, p , ψ − v h, p )|,

∀ v h, p ∈ Vh, p

= | R h (uh, p , ψ − ψh, p )| = |ah (u − uh, p , ψ − ψh, p )| via Lemma 2.2  h2p |u | H p+1 () |ψ| H p+1 () via Lemma A.3.

2

Theorem A.2. CG Global Estimate Error Bound Let ψh, p  ∈ Vh, p  be the adjoint finite element solutions of order p  . Then the error in the global estimate satisfies  |E − E˜|  h2p |u | H p +1 () |ψ| H p +1 () .

Proof.

|E − E˜| = | R h (uh, p , ψ − v h, p ) − R h (uh, p , ψh, p  − v h, p )|,

∀ v h, p ∈ Vh, p

= | R h (uh, p , ψ − ψh, p  )| ψ

= | R h (u − uh, p , ψh, p  )| via Lemma 2.2 ψ

= | R h (u − v h, p  , ψh, p  )|,

∀ v h, p  ∈ Vh, p 

= |ah (u − v h, p  , ψ − ψh, p  )|, h

2p 

via orthogonality

∀ v h, p  ∈ Vh, p 

|u | H p +1 () |ψ| H p +1 () via Lemma A.3.

via Lemma 2.2

2

A.5. Local Error Bound In this section Lemma A.2 is used to derive asymptotic convergence rates for ηκ and ηκ − η˜ κ for the CG discretization. In common with the global results, all that is necessary to bound the relevant formulae are Lemmas A.2 and 2.3. The integration by parts of the global residual to arrive at the localized residual means that rκ (u h, p , v h, p ) = 0, ∀ v h, p ∈ Vh, p thus the adjoint solution in p must be explicitly subtracted to achieve the optimal convergence. Theorem A.3. CG Local Error Bound

|ηκ |  h2p +d |u | H p+1 () |ψ| H p+1 () .

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

201

Proof.

|ηκ | = |rκ (uh, p , ψ − ψh, p )| = |aκ (u − uh, p , ψ − ψh, p )| via Lemma 2.3  h2p +d |u | H p+1 () |ψ| H p+1 () via Lemma A.2.

2

Theorem A.4. CG Local Estimate Error Bound 

|ηκ − η˜ κ |  h p + p +d |u | H p+1 () |ψ| H p+1 () . Proof.

|ηκ − η˜ κ | = |rκ (uh, p , ψ − ψh, p ) − rκ (uh, p , ψh, p  − ψh, p )| = |rκ (uh, p , ψ − ψh, p  )| = |aκ (u − uh, p , ψ − ψh, p  )| via Lemma 2.3 

 h p + p +d |u | H p+1 () |ψ| H p +1 () via Lemma A.2.

2

References [1] Mark Ainsworth, Bill Senior, An adaptive refinement strategy for hp-finite element computations, Appl. Numer. Math. 26 (1) (1998) 165–178. [2] Mark Ainsworth, J. Tinsley Oden, A posteriori error estimation in finite element analysis, Comput. Methods Appl. Mech. Eng. 142 (1) (1997) 1–88. [3] N. Douglas Arnold, Franco Brezzi, Bernardo Cockburn, Donatella L. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (5) (2002) 1749–1779. [4] I. Babuška, A. Miller, The post-processing approach in the finite element method – part 1: calculation of displacements, stresses and other higher derivatives of the displacements, Int. J. Numer. Methods Eng. 20 (6) (1984) 1085–1109. [5] I. Babuška, A. Miller, The post-processing approach in the finite element method – part 2: the calculation of stress intensity factors, Int. J. Numer. Methods Eng. 20 (6) (1984) 1111–1129. [6] I. Babuška, A. Miller, The post-processing approach in the finite element method – part 3: a posteriori error estimates and adaptive mesh selection, Int. J. Numer. Methods Eng. 20 (12) (1984) 2311–2324. [7] Ivo Babuška, Werner C. Rheinboldt, A-posteriori error estimates for the finite element method, Int. J. Numer. Methods Eng. 12 (10) (1978) 1597–1615. [8] Ivo Babuška, Werner C. Rheinboldt, Error estimates for adaptive finite element computations, SIAM J. Numer. Anal. 15 (4) (1978) 736–754. [9] Ivo Babuška, Werner C. Rheinboldt, A posteriori error analysis of finite element solutions for one-dimensional problems, SIAM J. Numer. Anal. 18 (3) (1981) 565–589. [10] Ivo Babuška, Barna A. Szabo, Norman I. Katz, The p-version of the finite element method, SIAM J. Numer. Anal. 18 (3) (1981) 515–545. [11] Randolph E. Bank, Alan Weiser, Some a posteriori error estimators for elliptic partial differential equations, Math. Comput. 44 (170) (1985) 283–301. [12] John W. Barrett, Charles M. Elliott, Total flux estimates for a finite-element approximation of elliptic equations, IMA J. Numer. Anal. 7 (2) (1987) 129–148. [13] F. Bassi, S. Rebay, GMRES discontinuous Galerkin solution of the compressible Navier–Stokes equations, in: Discontinuous Galerkin Methods: Theory, Computation and Applications, Springer, 2000, pp. 197–208. [14] Roland Becker, Rolf Rannacher, A feed-back approach to error control in finite element methods: basic analysis and examples, East-West J. Numer. Math. 4 (1996) 237–264. [15] Roland Becker, Rolf Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numer. 10 (2001) 1–102. [16] Rodolfo Bermejo, Jaime Carpio, A space-time adaptive finite element algorithm based on dual weighted residual methodology for parabolic equations, SIAM J. Sci. Comput. 31 (5) (2009) 3324–3355. [17] Susanne Brenner, Ridgway Scott, The Mathematical Theory of Finite Element Methods, Vol. 15, Springer Science & Business Media, 2007. [18] F. Brezzi, G. Manzini, D. Marini, P. Pietra, A. Russo, Discontinuous finite elements for diffusion problems, in: Atti Convegno in Onore di F. Brioschi (Milano 1997), Istituto Lombardo, Accademia di Scienze e Lettere, 1999, pp. 197–217. [19] Jaime Carpio, Juan Luis Prieto, Rodolfo Bermejo, Anisotropic “goal-oriented” mesh adaptivity for elliptic problems, SIAM J. Sci. Comput. 35 (2) (2013) A861–A885. [20] Hugh A. Carson, A Priori Analysis of Global and Local Output Error Estimates for CG, DG and HDG Finite Element Discretizations, Master’s thesis, Mass. Inst. of Tech., Department of Aeronautics and Astronautics, 2016. [21] Paul Castillo, Bernardo Cockburn, Ilaria Perugia, Dominik Schötzau, An a priori error analysis of the local discontinuous Galerkin method for elliptic problems, SIAM J. Numer. Anal. 38 (5) (2000) 1676–1706. [22] Marco A. Ceze, Krzysztof J. Fidkowski, High-order output-based adaptive simulations of turbulent flow in two dimensions, AIAA J. 54 (2) (2016) 2611–2625, American Institute of Aeronautics and Astronautics. [23] Bernardo Cockburn, Bo Dong, Johnny Guzmán, Marco Restelli, Riccardo Sacco, A hybridizable discontinuous Galerkin method for steady-state convection–diffusion-reaction problems, SIAM J. Sci. Comput. 31 (5) (2009) 3827–3846. [24] Bernardo Cockburn, Jayadeep Gopalakrishnan, Raytcho Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2) (2009) 1319–1365. [25] Bernardo Cockburn, Jayadeep Gopalakrishnan, Francisco-Javier Sayas, A projection-based error analysis of HDG methods, Math. Comput. 79 (271) (2010) 1351–1367. [26] Ngoc Cuong Nguyen, Jaime Peraire, Bernardo Cockburn, An implicit high-order hybridizable discontinuous Galerkin method for linear convection– diffusion equations, J. Comput. Phys. 228 (9) (2009) 3232–3254. [27] Krzysztof J. Fidkowski, David L. Darmofal, Review of output-based error estimation and mesh adaptation in computational fluid dynamics, AIAA J. 49 (4) (2011) 673–694. [28] Michael B. Giles, Niles A. Pierce, Adjoint error correction for integral outputs, in: Error Estimation and Adaptive Discretization Methods in Computational Fluid Dynamics, Springer, 2003, pp. 47–95.

202

H.A. Carson et al. / Applied Numerical Mathematics 118 (2017) 182–202

[29] Michael B. Giles, Niles Pierce, Endre Süli, Progress in adjoint error correction for integral functionals, Comput. Vis. Sci. 6 (2–3) (2004) 113–121. [30] Michael B. Giles, Endre Süli, Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality, Acta Numer. 11 (2002) 145–236. [31] Ralf Hartmann, Joachim Held, Tobias Leicht, Florian Prill, Error estimation and adaptive mesh refinement for aerodynamic flows, in: ADIGMA-a European Initiative on the Development of Adaptive Higher-Order Variational Methods for Aerospace Applications, Springer, 2010, pp. 339–353. [32] Todd A. Oliver, David L. Darmofal, Analysis of dual consistency for discontinuous Galerkin discretizations of source terms, SIAM J. Numer. Anal. 47 (5) (2009) 3507–3525. [33] Niles A. Pierce, Michael B. Giles, Adjoint recovery of superconvergent functionals from PDE approximations, SIAM Rev. 42 (2) (2000) 247–264. [34] Thomas Richter, Thomas Wick, Variational localizations of the dual weighted residual estimator, J. Comput. Appl. Math. 279 (2015) 192–208. [35] Rüdiger Verfürth, A posteriori error estimation and adaptive mesh-refinement techniques, J. Comput. Appl. Math. 50 (1) (1994) 67–83. [36] Masayuki Yano, An Optimization Framework for Adaptive Higher-Order Discretizations of Partial Differential Equations on Anisotropic Simplex Meshes, PhD thesis, Massachusetts Institute of Technology, Department of Aeronautics and Astronautics, June 2012. [37] Masayuki Yano, David Darmofal, An optimization framework for anisotropic simplex mesh adaptation: application to aerodynamic flows, (2012), no. 2012–0079.