A posteriori error estimators for discontinuous Galerkin method for diffusion problems, based on the hypercircle method

In this work, we apply the hypercircle method to Discontinuous Galerkin (DG) approximations of second order diffusion problems featuring inhomogeneous Dirichlet and Neumann boundary conditions. We focus on the interior penalty discontinuous Galerkin (IPDG) approximations of diffusion problems in primal variational formulation and produce a Prager–Synge theorem for such DG methods. Using the hypercircle method, we derive an a posteriori error estimator in terms of an equilibrated flux. The estimator is proven to be reliable and efficient. Numerical results are presented which illustrate the estimator’s performance.


Introduction
The hypercircle method, also known as the two-energy principle, was introduced in 1947 by Prager and Synge [25,26] as a method to estimate the error arising from the approximation of boundary value problems for diffusive partial differential equations. The method is based on splitting the original problem into two easyto-solve problems. These problems are solved separately, a common solution is found which represents the solution of the original problem in a suitably-defined function space. In this case, the solution to the original problem is located on or inside a hypercircle in that function space. If the bilinear form of the variational formulation of the original problem is positive definite, then the hypercircle is bounded [25] and the location of the solution on a hypercircle provides an upper and lower bound for the solution on any neighbourhood of any prescribed point in the domain. If we choose the center of the hypercircle as an approximation of the original solution, then its error in a mean-square sense, is easily known.
However, the hypercircle method was forgotten for many years. It resurfaced again when there was a need for efficient and reliable a posteriori error estimators for the finite element approximation of diffusion-boundary value problems (cf. e.g. [2,6,8,27]). The equilibrated a posteriori error estimator based on the hypercircle method has a big advantage over the standard residual-type a posteriori error estimator that does not contain any generic constants except for possible data oscillations, which leads to better effectivity indices.
The hypercircle method was designed for conforming finite elements approximations, but it can be used as well as for nonconforming approaches such as discontinuous Galerkin (DG) methods. Indeed, in [7], a Prager-Synge theorem has been formulated and used to derive an equilibrated a posteriori error estimator for DG approximations of the Poisson equation. This is done by constructing equilibrated flux and post-processing the DG approximation via a related mixed formulation from [4], involving suitably defined numerical fluxes across the interfaces of the finite element meshes. Related approaches can be found in [1,12,16], which however are not directly based on a Prager-Synge-like theorem.
DG methods have been used since the 1970s and have been so popular since then, [13,14,22] and the references therein, for their flexibility in local approximation of the solution of partial differential equations (PDE's), flexibility in adaptive local mesh refinement, and for their good stability properties. When it comes to adaptive methods based on equilibrated a posteriori error estimators, DG methods have a big advantage over C 0 -conforming finite elements. In C 0 -conforming finite elements, the equilibration is done on patches of elements (union of elements sharing a common nodal point), as in [8], whereas for interior penalty DG (IPDG) methods equilibration is performed element-wise which is less expensive to compute [7,12].
In this work, we extend the work done in [7] to cover DG approximations of diffusion problems with non-zero boundary conditions and high contrast coeffiecients. Further, [7] does not address data oscillations due to approximations, this is another point that we consider here. For such problems, we apply the hypercircle method to IPDG methods of any polynomial degree, we formulate a Prager-Synge theorem and derive an a posteriori error estimator in terms of an equilibrated flux vector. The reliability of the error estimator follows directly from the Prager-Synge theorem, whereas its efficiency can be established by showing that the estimator is bounded from above by a residual-type a posteriori error estimator for IPDG methods whose efficiency has been shown in the literature, (cf., e.g., [5,17,19]). This paper is organised as follows, a brief introduction to the DG approximation is given followed by focus on the IPDG methods in Sect. 2. We introduce a posteriori error estimate based on the Prager and Synge theorem in Sect. 3.1 and we show the post processing of the DG methods to conferring FE we used in Sect. 3.2. We introduce a discussion of possible Data oscillations. In Sect. 3.3 and in Sect. 3.4 we construct an equilibrated flux. In Sect. 3.5 we show the efficiency of our estimator. Sect. 4.1 consists of an introduction to the adaptive method we are using and some numerical results are given in Sect. 4.2 to show the performance of the estimator. our conclusions were discussed in section 5.

DG methods for diffusion problems
We will use standard notations from Lebesgue and Sobolev space theory [6,9]. In particular, for a bounded domain Ω ⊂ R d for d ≥ 2 and D ⊂ Ω we denote by (·, ·) 0,D and · 0,D the L 2 -inner product and the associated L 2 -norm. We further refer to H 1 (Ω) as the Sobolev space with norm · 1,Ω and seminorm | · | 1,Ω and to H 1 2 (Γ ), Γ ⊆ ∂Ω, as the associated trace space, whereas H 1 0,Γ (Ω) stands for the subspace of H 1 -functions with vanishing trace on Γ . Moreover H (div, Ω) denotes the Hilbert space of vector fields q ∈ L 2 (Ω) 2 such that ∇ · q ∈ L 2 (Ω) equipped with the graph norm.
We start by assuming Ω ⊂ R 2 to be a bounded polygonal domain consisting of pairwise disjoint polygonal subdomains Ω k , 1 k m. We refer to Γ = Γ D ∪ Γ N , Γ D ∩ Γ N = ∅, as the boundary of Ω, consisting of a Dirichlet part Γ D and a Neumann part Γ N . We consider the diffusion problem Here, we suppose a to be a piecewise constant function with a| Ω k > 0, 1 k m, and f ∈ L 2 (Ω), . We define then, the weak formulation of (1a)-(1c) as : we seek u ∈ V such that for all For the IPDG approximation of the boundary value problem (1a)-(1c) we consider a locally quasi-uniform simplicial triangulation T h of a computational domain which aligns with the subdivision of Ω. We denote by the set of interior edges in T h , Dirichlet boundary edges, and Neumann edges respectively. We further denote by h K , K ∈ T h the diameter of K , and by h E , E ∈ E h the length of edge E. Further, for l ∈ N, we refer to P l (K ) as the set of polynomials of degree at most l on K . Due to the local quasi-uniformity of the triangulation, there exist constants 0 < c < C such that and define the average and jump across E as usual, according to Moreover, we refer to n E , E ∈ E h , E = K + ∩ K − , as the unit normal vector pointing from K + to K − and to as the exterior unit normal vector. We introduce now the DG spaces DG methods can be derived from a mixed formulation of (1a)-(1c), [4] and amount to the computation of whereũ ∂ K andp ∂ K are suitably defined numerical fluxes.

The IPDG method
the IPDG method is obtained by the following specification of the numerical fluxesũ where α > 0 is a penalty parameter. This choice of the numerical fluxes allows us to eliminate the dual variable p h from (7a),(7b) by choosing p h = a∇u h and q h = a∇v h in (7a). We are thus led to the primal variational formulation of the IPDG method: we seek u h ∈ V h such that, where the bilinear form a I P h (·, ·) : and the functional l(·) : V h → R reads as follows: Remark 2.1 Another choice of the numerical fluxesũ E andp E , E ∈ E h , gives rise to a different DG approximation. We refer to [4] for an overview. In the remaining of this paper, we will focus on the IPDG approximation, but note that both the derivation of the equilibrated a posteriori error estimator and its analysis in terms of reliability and efficiency can be done the same way.
The IPDG bilinear form (10) and the functional in (11) . This can be cured by using a lifting operator L : [18], defined according to As has been shown e.g. in [18,23], the lifting operator is continuous, so there exists a constant C L ≥ 0 depending only on the shape regularity of the triangulation such that Using the lifting operator, we define the bilinear formã I P h (·, ·) : It is easy to show that the following holds

Coercivity and continuity of the IPDG bilinear form
It is not difficult to show, cf., [19,20], that for sufficiently large α > 0 the IPDG norm and the mesh-dependent energy norm are equivalent, i.e., there exists a positive constant γ such that On the other hand, there exists a constant Γ > γ such that for any α > 0 In particular, it follows from (16a) and (16b) that for a sufficiently large α > 0, the bilinear form a I P h (·, ·) is continuous and coercive, thus the IPDG problem (9) admits a unique solution u h ∈ V h . It has been shown before [24] that to ensure the coercivity of the bilinear form and the continuity of the IPDG approximation, the penalty parameter has to be chosen according to α = a max O((k + 1) 2 ).

A posteriori error estimate based on the Prager-Synge theorem
Theorem 3.1 Prager-Synge theorem for diffusion Problems, two-energy principle Let p ∈ H (div, Ω) and v ∈ V . Further, let u ∈ V be the solution of (2) and suppose that p satisfies the equilibrium condition and the boundary condition Then it holds Proof Observing (2), (17), (18), and the fact that v − u ∈ H 1 0,Γ D (Ω), partial integration yields This orthogonality relation implies which proves (19).

A posteriori error estimator
The Prager-Synge theorem (3.1) can be applied to the IPDG approximation involving an equilibrated flux vector p eq h , if we replace f, u D , and u N by f h , u D,h , and u N ,h , respectively. In particular, we denote byû h ∈ V h the IPDG approximation of the weak solutionû ∈ H 1 (Ω) of the diffusion problem Note that in this case, the associated mixed formulation reads where the numerical fluxesû ∂ K andp ∂ K are given bŷ 3.2 Postprocessing to conforming finite element spaces h be a set of Lagrange nodal points for the elements in V h and let κ i be the number of elements that share the nodal point The following estimate has been provided by Theorem (2.2) in [19]: A flux vector p eq h ∈ V h is called an equilibrated flux vector, if p eq h ∈ H (div, Ω) and if it satisfies the equilibrium equation and the boundary condition As we shall see below, the Prager-Synge theorem (3.1) gives rise to an equilibrated a posteriori error estimator where the element-related terms η eq K ,i , 1 i 2, and edge-related terms η eq E,i , 1 i 2, are given as follows:

Data oscillations
We provide the following auxiliary result concerning the data oscillation due to the approximation of f, u D , and u N by f h , u D,h , and u N ,h respectively.

Lemma 3.2 Let z h ∈ V h be the IPDG approximation of the weak solution z
Then it holds Proof The Poincaré-Friedrichs inequality asserts Moreover, the trace inequalities for functions v K ∈ P k (K ) (cf.,e.g., [28]) In view of (30), (31a), (31b), and using (3) we deduce from (29) The assertion (28) then follows from (32) by a suitable application of Young's inequality.
In view of the previous result, the data oscillations will be denoted by where osc K ( f ), osc E (u D ), and osc E (u N ) are given by  (22). Finally, let η eq h be the equilibrated a posteriori error estimator given by (25) and let osc h be the data oscillations (33). Then there exists a constant C > 0 independent of h such that Proof Let u h ∈ V h be the solution of (9). Then we have Since z h := u h −û h is the IPDG approximation of (27a)-(27c), the second term on the right-hand side in (35) can be estimated from above by (28) and thus gives rise to the data oscillation term in (34). For the first term on the right-hand side in (35) we obtain ⎛ Again, the last term on the right-hand side can be estimated from above by (28). The Prager-Synge theorem 3.1 with v =û conf h and p = p eq h gives Using (36) and (37) in (35) allows to conclude (34).

Remark 3.4
We emphasize that the constant C in front the data oscillations is the only generic constant in the reliability estimate (34).

Construction of an equilibrated flux
We consider now the construction of an equilibrated flux p eq h ∈ BDM k (Ω, T h ) where and BDM k (K ) stands for the classical Brezzi-Douglas-Marini element on K , [10]. Denoting by λ K i , 1 ≤ i ≤ 3, the barycentric coordinates of K ∈ T h , we refer to b K := λ K 1 λ K 2 λ K 3 as the element bubble function which vanishes on ∂ K . For q h | K ∈ BDM k (K ), we define its degrees of freedom(DOF) according to Each vector field q h ∈ BDM k (K ) is uniquely determined by (38a)-(38c), [11]. An immediate consequence is the following result: Proof The proof of (39) follows by a standard scaling argument. Proof It follows from Gauss's theorem and (38b) that for p k−1 ∈ P k−1 (K ) it holds Since ∇ ·p eq h and f h belong to P k−1 (K ), from (41) we deduce (24a). For E ∈ E N h , we choose p k−1 ∈ P k−1 (K ) with supp( p k−1 ) = E. By the definition of the numerical fluxp ∂ K and the fact that both n E · p eq h | E and u N ,h | E belong to P k−1 (E), (24b) also follows from (41).

Efficiency of the equilibrated error estimator
For the symmetric IPDG approximation of second-order diffusion problems on simplicial triangulations, residual-type error estimators have been derived and analysed by many authors, [5,17,19,20]. These type of estimators are of the form, where the element residuals η res K are given by The efficiency estimate from [20] shows The efficiency of the equilibrated a posteriori error estimator η eq h follows from (43) and the following result.

Lemma 3.7 Let η
eq h be given by (25) and let η res h be the residual-type a posteriori error estimator (42). Then there holds Proof To estimate η eq K ,1 , K ∈ T h , from above, we apply Lemma3.5 to p eq h − a∇û h . From (21b) and (40a) it follows that On the other hand, due to (24a) we have Finally, in view of (40c) for each K ∈ T h it holds Using (45)-(47), from Lemma3.5 we obtain Moreover, using (23) we find

Adaptive refinement
The a posteriori error estimation's main purpose is to indicate the elements where the error is large. Then we use this information to locally refine those elements, and repeat the finite element computation. The adaptive approach consists of successive cycles of the steps In the step SOLVE, we compute the solutionû h ∈ V h of the IPDG approximation. In the second step ESTI-MATE, we compute the local components η eq K ,i and η eq E,i , 1 ≤ i ≤ 2, of the error estimator η eq h (cf. (26a)-(26d)). To come up with a purely element-based marking strategy, we set and use standard Dörfler marking [15] in step MARK: Given some constant 0 < θ ≤ 1, we choose a set In the final step REFINE, marked elements K ∈ M are subdivided by newest vertex bisection. The performance of the equilibrated a posteriori error estimator will be illustrated in the sequel by numerical examples.

Numerical examples
We present now numerical examples which validate the optimality and robustness of the equilibrated a posteriori error estimator based on the hypercircle method. We begin with verifying the optimal behaviour of the a posteriori estimator by computing experimental order of convergence by means of manufactured solutions. The efficiency of the a posteriori estimator is evaluated for the classical L-shaped domain problem and in the case of highly discontinuous coefficients. The numerical results were obtained using the computational framework FEniCS, [3,21] and the references within.

Validation, optimality of the estimator
We first validate the optimality of the a posteriori estimator by computing convergence rates in the following case. We consider (1a)-(1c) in Ω = (−8, 8) 2 such that the boundary consists of the Dirichlet boundary Γ D = ∂Ω. For simplicity, we take constant coefficient a(x) ≡ 1, ∀x ∈ Ω and we select f, u D , and u N such that is the exact solution in Ω. The optimality of the estimator and the effectivity index are also shown For the numerical approximation we use cubic elements for u 1 and we start the realizations with an initial grid and then perform a series of uniform refinements, calculating the error in the L 2 and IPDG norms. The experimental order of convergence(EOC) is computed as where E 1 , E 2 and h 1 , h 2 are the errors and maximum cell sizes, respectively, for two different realizations. The numerical results for u 1 are shown in Table 1. For u 1 the expected convergence rates of 4 and 3 are observed for the L 2 -norm and I P DG norm respectively. The optimality of the a posteriori estimator η h is also depicted The error in IPDG norm and the corresponding effectivity index are shown   Fig. 1 the approximate solution is depicted for the three levels of uniform mesh refinement.

Adaptivity and estimator performance
Next we investigate the performance of the a posteriori estimator in terms of its ability to detect regions of the solutions where high resolution is required to maintain the overall accuracy of the method. We consider the aforementioned example u 1 . We performed adaptive refinement of the mesh with θ = 0.5 in the Dörfler marking, and calculated the error in the IPDG norm and the value of the estimator. The results are shown in Table 2. As in the uniform refinement, the effectivity index stabilizes around the same value E I ∼ 5. In Fig. 2 the approximate solution is depicted for the three levels of mesh refinement.

Discontinuous coefficients and adaptivity
Next, we study the robustness of the estimator in the case where the coefficient a in (1a)-(1c) is discontinuous in the domain. Let Ω = (0, 1) 2 be such that The discontinuity of a is along the two lines x = 0.5 and y = 0.5, that is we define: a 1 = a| Ω 1 = 100, a 2 = a| Ω 2 = 1, a 3 = a| Ω 3 = 100, a 4 = a| Ω 4 = 1. The domain Ω with its subdivision and the coefficient a are shown in Fig. 3. The boundary consists of the Dirichlet boundary Γ D = ∂Ω. We also choose f, u D , and u N such that u 2 is the exact solution of (1a)-(1c).
We use linear elements and we performed several levels of adaptive refinement of the mesh with θ = 0.5 in the Dörfler marking, and calculated the error in the IPDG norm, the value of the estimator and the effectivity index. The results are shown in Table 3. The effectivity index reaches a stable value of E I ∼ 19. In Fig. 4 the approximate solution is depicted for the three levels of mesh refinement.
In the second example with discontinuous coefficients we consider again Ω = (0, 1) 2 such that Ω = Fig. 3 Subdivision of Ω and distribution of a In this case the discontinuity of a is along the two main diagonals of the domain, cf. Fig. 5. We choose a 1 = a| Ω 1 = 1, a 2 = a| Ω 2 = 100, a 3 = a| Ω 3 = 1, a 4 = a| Ω 4 = 100. The boundary Γ = ∂Ω consists of the Dirichlet boundary Γ D = ∂Ω. We also choose f, u D , and u N such that is the exact solution of (1a)-(1c).
We use linear elements and we performed several levels of adaptive refinement of the mesh with θ = 0.5 in the Dörfler marking, and calculated the error in the IPDG norm, the value of the estimator and the effectivity index. The results are shown in Table 4. The effectivity index reaches a stable value of E I ∼ 18. In Fig. 6 the approximate solution is depicted for the three levels of mesh refinement.
The last two examples show how the estimator helps in marking and refining the cells around the interfaces between the domains where the value of a varies, cf. Figs. 3 and 5. The effectivity indices shows clearly the efficiency and reliability of the a posteriori error estimator. The refinement improves with the increase in the polynomial degree k as shown below in Figs. 7 and 8.

Conclusions
Based on the hypercircle method also known as the Prager-Synge theorem, we have developed and analyzed an equilibrated a posteriori error estimator for IPDG approximations of second order diffusion problems. The error estimator relies on the construction of an equilibrated flux and has been shown to be both efficient and reliable. Numerical results for some selected examples confirm the theoretical findings.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding No funding was received to assist with the preparation of this manuscript.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Data availability statement Not applicable.