Igor Mozolevski

Serge Prudhomme

A robust goal-oriented estimator based on the construction of equilibrated fluxes for discontinuous Galerkin finite element approximations of convection-diffusion problems

Keywords: Convection-diffusion equation, Discontinuous Galerkin finite element method, A posteriori error estimation, Adaptive refinement, Asymptotically exact error estimate, Quantity of interest 2000 MSC: 65N15, 65N30, 65N50

We propose an a posteriori error estimator with respect to quantities of interest for discontinuous Galerkin approximations of convection-diffusion boundary-value problems. The error estimator is based on the construction of equilibrated fluxes in Raviart-Thomas finite element spaces and on the solution of the dual problem. We show that it is asymptotically exact in both the elliptic and hyperbolic regimes if the dual problem is approximated by a discontinuous Galerkin method of order one greater than that of the primal problem. We show in this case that the effectivity index behaves as (1+Pe 1/2 )o(h), where Pe is the Péclet number and h the mesh diameter. It follows that the quality of the effectivity index may deteriorate for large values of Pe, but we put in evidence that it suffices to increase the approximation order of the dual problem to keep the effectivity index close to unity even on coarse meshes. Two-dimensional numerical examples demonstrate the robustness of the error estimator in both the diffusion and advection regimes.

Introduction

In most convection-diffusion problems of practical interest, solutions exhibit thin boundary or internal layers in which gradients may assume extremely large values [START_REF] Roos | Robust numerical methods for singularly perturbed differential equations[END_REF]. Therefore, accurate numerical approximations of the boundary layer solution by the finite element method call for efficient adaptive strategies of the mesh size and/or the approximation order. A fair amount of literature has been indeed devoted to a posteriori error estimation for convection-diffusion-reaction problems in recent years (see e.g. [START_REF] Verfürth | A posteriori error estimation techniques for finite element methods[END_REF] and references therein).

Discontinuous Galerkin (dG) finite element methods are widely used for the solution of convection-diffusion problems in the advection-dominated regime owing to their local conservative properties and their high degree of localisation. Moreover, dG methods possess the attractive property that they spare one from using excessive numerical stabilization, hence reducing numerical dissipation over classical FE methods [START_REF] Houston | Discontinuous hp-finite element methods for advection-diffusion problems[END_REF]. A posteriori error estimates with respect to energy norms for dG approximations of convection-reaction-diffusion problems have been presented in [START_REF] Ern | A posteriori energy-norm error estimates for advectiondiffusion equations approximated by weighted interior penalty methods[END_REF][START_REF] Schötzau | A robust a-posteriori error estimator for discontinuous Galerkin methods for convection-diffusion equations[END_REF][START_REF] Ern | Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF][START_REF] Creusé | A posteriori error estimator based on gradient recovery by averaging for convection-diffusion-reaction problems approximated by discontinuous Galerkin methods[END_REF][START_REF] Baccouch | Asymptotically exact a posteriori LDG error estimates for onedimensional transient convection-diffusion problems[END_REF]. However, the literature on error estimates with respect to quantities of interest is rather limited for these problems [START_REF] Harriman | hp-version discontinuous Galerkin methods with interior penalty for partial differential equations with nonnegative characteristic form[END_REF], while there exist several works on the subject for the case of conforming finite element methods (see e.g. [START_REF] Paraschivoiu | A hierarchical duality approach to bounds for the outputs of partial differential equations[END_REF][START_REF] Sauer-Budge | Computing bounds for linear functionals of exact weak solutions to the advection-diffusion-reaction equation[END_REF][START_REF] Parés | Exact bounds for linear outputs of the advectiondiffusion-reaction equation using flux-free error estimates[END_REF][START_REF] Kuzmin | Goal-oriented a posteriori error estimates for transport problems[END_REF][START_REF] Parés | Computable exact bounds for linear outputs from stabilized solutions of the advection-diffusion-reaction equation[END_REF] and references therein) following early contributions on goal-oriented error estimation [START_REF] Eriksson | Introduction to adaptive methods for differential equations[END_REF][START_REF] Paraschivoiu | A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations[END_REF][START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF][START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF][START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF].

In practice, one is interested in robust a posteriori goal-oriented error estimates that are fully computable from the problem data (without introducing unknown theoretical constants) and can be reliably used in adaptive methods independently of the Péclet number Pe for the convection-diffusion problem under consideration. We have recently proposed in [START_REF] Mozolevski | Goal oriented estimates for elliptic equations based on equilibrated fluxes reconstruction[END_REF] an approach dealing with goal-oriented error estimation for finite element approximations of elliptic problems that combines the dual-weighted residual method and equilibrated-flux reconstruction. However, one may expect for convection-diffusion problems that the sharpness of the goal-oriented error bound may significantly degrade with increasing Péclet number, even if the bounding property is retained. The objective of the paper is thus to extend the technique developed in [START_REF] Mozolevski | Goal oriented estimates for elliptic equations based on equilibrated fluxes reconstruction[END_REF] to dG approximations of convection-diffusion problems in order to propose error estimates that are robust with respect to the Péclet number and are asymptotically exact, i.e. the effectivity index remains close to unity whenever the mesh size approaches zero and/or the approximation order increases. In particular, we show that the proposed error estimator is asymptotically exact and that the effectivity index converges to unity at a rate proportional to the square root of the Péclet number (when Pe becomes large).

The paper is organized as follows. In Section 2, we present the model problem and the corresponding dual problem in the continuous setting. In Section 3, we briefly introduce the discontinuous Galerkin method for convection-diffusion problems and discuss the reconstruction of the total (diffusive and advective) flux within the Raviart-Thomas finite element space. An error representation for the quantity of interest is introduced in Section 4 using the global equilibrated fluxes obtained from the primal dG solution and the discrete diffusion fluxes derived from the dual dG solution. Section 5 describes numerical examples chosen to demonstrate the robustness of the estimator with respect to the Péclet number and its asymptotical exactness for smooth solutions involving boundary layers. We also illustrate the potential usage of the error estimator for mesh adaptation. Concluding remarks are included in Section 6.

Model problem

Let Ω be an open bounded polyhedral domain in R d , d = 2 or 3, with boundary Γ. We consider the boundary-value problem composed of the convection-diffusion equation and homogeneous Dirichlet conditions:

∇ • (-ε∇u + βu) = f in Ω, u = 0 on Γ. ( 1
)
Here diffusion coefficient ε is a strictly positive constant, vector field β ∈ [Lip(Ω)] d is Lipschitzian and divergence-free, ∇ • β = 0, in Ω, and source term f is assumed in L 2 (Ω).

Moreover, we suppose that the diffusion coefficient satisfies ε ≪ 1 and that the diameter L of domain Ω and β ∞ = β L ∞ (Ω) are of order unity, so that the Péclet number of the problem, Pe = L β ∞ ε -1 , is of order ε -1 . Note that we choose here homogeneous Dirichlet boundary conditions only for the sake of simplicity in the exposition; other types of boundary conditions, such as inhomogeneous Dirichlet or mixed Dirichlet-Neumann boundary conditions, could be considered as well.

A weak formulation of Problem (1) reads:

Find u ∈ H 1 0 (Ω) such that B(u, v) = F (v), ∀v ∈ H 1 0 (Ω) (2)
where bilinear form B(•, •) defined on [H 1 0 (Ω)] 2 and linear form F (•) defined on H 1 0 (Ω) are given by:

B(u, v) = Ω (ε∇u -βu) • ∇v, (3)
F (v) = Ω f v. (4)
Problem (2) admits a unique weak solution in H 1 0 (Ω) owing to the Lax-Milgram theorem, since bilinear form B is continuous and coercive on [H 1 0 (Ω)] 2 and linear functional is continuous on H 1 0 (Ω). Moreover, the following weak stability estimate holds:

√ ε∇u L 2 (Ω) ≤ C f L 2 (Ω) .
Given a linear continuous functional Q on L 2 (Ω), that is, the quantity of interest, the dual problem in weak form reads:

Find p ∈ H 1 0 (Ω) such that B(v, p) = Q(v), ∀v ∈ H 1 0 (Ω) (5)
Above problem can be recast in strong form, thanks to the regularity of the dual solution, as:

∇ • (-ε∇p -βp) = q in Ω, p = 0 on Γ, ( 6
)
where q is the Riesz representer of Q in L 2 (Ω), i.e.

Q(v) = Ω qv, ∀v ∈ L 2 (Ω). ( 7
)
Similarly to the primal problem, the dual problem is well posed in H 1 0 (Ω) and the weak stability estimate holds:

√ ε∇p L 2 (Ω) ≤ C q L 2 (Ω) .

Discontinuous Galerkin method for convection-diffusion equations

Let T h , h > 0, be a shape regular family of simplicial meshes, with no hanging nodes, on Ω (see e.g. [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]), where h = max T∈T h h T denotes the mesh size and h T is the diameter of mesh element T . The set of all mesh faces is denoted by F h , which is decomposed into the set of all interfaces between adjacent elements, F i,h and the set of all boundary faces, F Γ,h . We define on F h the piecewise constant vector field n : F h → R d , such that n| F corresponds to one arbitrarily chosen unit normal vector to face F if F ∈ F i,h and coincides with the external normal to Γ if F ∈ F Γ,h . We also introduce the piecewise constant function of face diameters h F h :

F h → R, where h F h | F is the diameter of F ∈ F h .
For a given partition T h , we define the broken Sobolev space of order l ∈ N as:

H l (T h ) = u ∈ L 2 (Ω) : u| T ∈ H l (T ), ∀T ∈ T h . ( 8
)
Since a function v in H 1 (T h ) admits (possibly different) traces at faces F ∈ F h , we define the average and jump of v at a face F ∈ F i,h by

{{v}} = 1 2 v| T -+ v| T + | F , [[v]] = v| T --v| T + | F , (9)
where the adjacent elements T -and T + sharing face F are chosen in such a way that n| F points from T -to T + . We extend the definition to boundary edges

F ∈ F Γ,h , F ⊂ ∂T , by setting {{v}} = [[v]] = v| T | F . (10)
We now introduce the space V k h of piecewise polynomial functions of order k ∈ N on Ω:

V k h := {v h ∈ L 2 (Ω) : v h | T ∈ P k (T ), ∀T ∈ T h }, (11)
where P k (T ) denotes the vector space of polynomial functions with total degree less than or equal to k on T . Finally, for any ω ⊆ T h and any γ ⊆ F h , we denote:

ω f = T ∈ω T f, ∀f ∈ L 2 (ω), (12)
γ g = F ∈γ F g, ∀g ∈ L 2 (γ). ( 13
)
The symmetric dG formulation of order k ∈ N for the primal problem (1) reads (see e.g. [START_REF] Houston | Discontinuous hp-finite element methods for advection-diffusion problems[END_REF]):

Find u h ∈ V k h such that B h (u h , v h ) = F (v h ), ∀v h ∈ V k h , (14)
where the bilinear form B h is defined as

B h (u h , v h ) = T h ε∇ h u h • ∇ h v h - T h u h (β • ∇ h v h ) - F h {{(ε∇ h u h ) • n}}[[v h ]] + [[u h ]]{{(ε∇ h v h ) • n}} + F h δ ε [[u h ]][[v h ]] + F i,h (β • n){{u h }} + δ β [[u h ]] [[v h ]] + F Γ,h (β • n) + u h v h , ∀u h , v h ∈ V k h . (15)
The discrete gradient ∇ h is defined locally on

V k h as (∇ h u h )| T = ∇(u h | T ).
The stabilization functions δ ε , δ β : F h → R for the elliptic and convective terms of the operator, respectively, are chosen piecewise constant as

δ ε | F = δεk 2 h -1 F , F ∈ F h , and δ β | F = |β • n|/2, F ∈ F i,h
, where the penalty parameter δ > 0 should be taken sufficiently large to ensure that the formulation is stable (see e.g. [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]). Note that the penalization with respect to the advective term corresponds to the usual up-winding stabilization for hyperbolic operators. Moreover, the Dirichlet boundary condition is weakly imposed so that it degenerates to an in-flow boundary condition in the hyperbolic limit ε = 0. In the sequel we will use the same discontinuous Galerkin method (possibly of different order) for the discretization of the dual problem [START_REF] Eriksson | Introduction to adaptive methods for differential equations[END_REF] with the respective change of sign of β:

Find p h ∈ V k h such that B h (v h , p h ) = Q(v h ), ∀v h ∈ V k h . (16)
For v ∈ H l (T h ) ∩ V k h , we introduce the mesh-dependent norm

|v | 2 = T h ε|∇ h v| 2 + β ∞ L v 2 + F h ε h F h + |β • n| [[v]] 2 + h F h ε{{∇ h }} 2 . (17)
In the error analysis presented below, we will use the following assumptions on the problem data and on the solution to the primal and dual problems.

Assumption 1. The solutions u and p to the primal and dual boundary-value problems, respectively, are assumed sufficiently smooth, i.e. u, p ∈ H 2 (T h ), u, p, n • ∇u, and n • ∇p are continuous across each face F ∈ F i,h .

Assumption 2. The term β • ∇ h v h is assumed to be in V k h for all v h ∈ V k h .
Under these assumptions, the following theorem provides an a priori error estimate for high-order discontinuous Galerkin approximations. We refer the reader to [START_REF] Houston | Discontinuous hp-finite element methods for advection-diffusion problems[END_REF] for the proof.

Theorem 1. Let Assumptions 1 and 2 hold, let u ∈ H k+1 (Ω), k ≥ 1, be the solution to (1), and let u h be the discontinuous Galerkin solution of order k to the discrete problem [START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF]. Then there exists a constant C, which depends on Ω, the shape regularity of T h , and the approximation order k, but is independent of ε, β and h, such that

|u -u h | ≤ C ε 1 2 + β 1 2 ∞ h 1 2 h k u H k+1 . ( 18
)
Note 1. From now on, we will denote by C a generic constant that may depend on Ω, the shape regularity of T h , and the approximation order k, as in the statement of the theorem.

In practical applications, one is generally interested in quantifying the approximation error u -u h in terms of some given goal functional [START_REF] Ern | A posteriori energy-norm error estimates for advectiondiffusion equations approximated by weighted interior penalty methods[END_REF] rather than the dG norm [START_REF] Parés | Exact bounds for linear outputs of the advectiondiffusion-reaction equation using flux-free error estimates[END_REF]. We now consider a representation of the error in the goal functional

E Q (u h ) = Q(u) -Q(u h ) = Q(u -u h )
for given discrete solution u h . An a priori error estimate for E Q can be obtained using the adjoint consistency and the Galerkin orthogonality of the discontinuous Galerkin method, as follows. We first introduce the mesh-dependent dual problem:

Find p ∈ H 2 (T h ) such that B h (v, p) = Q(v), ∀v ∈ H 2 (T h ). ( 19
)
Integrating by parts [START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF], one can show that the solution to this problem is in fact the solution to the dual problem (6), i.e. p = p, see [START_REF] Harriman | hp-version discontinuous Galerkin methods with interior penalty for partial differential equations with nonnegative characteristic form[END_REF]. In such a case, namely when the dual problem obtained by transposing the arguments in the bilinear form B h results in the boundary-value problem involving the formal adjoint of the partial differential operator, the approach is thus referred to as being adjoint-consistent, see e.g. [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF]. Under Assumption 1, the dG method ( 14) is consistent and the following Galerkin orthogonality properties hold:

B h (u -u h , v h ) = 0, ∀v h ∈ V k h , B h (w h , p -p h ) = 0, ∀w h ∈ V k h .
where u h and p h are solutions to ( 14) and ( 16), respectively. Using the Galerkin orthogonality, we obtain the standard error representation with respect to the quantity of interest:

E Q (u h ) = Q(u) -Q(u h ) = B h (u -u h , p) = B h (u -u h , p -p h )
from which an a priori error estimate follows (see [START_REF] Harriman | hp-version discontinuous Galerkin methods with interior penalty for partial differential equations with nonnegative characteristic form[END_REF] for a proof and additional details).

Theorem 2. Let Assumptions 1 and 2 hold and let u, p ∈ H k+1 (T h ), k ≥ 1 be the solutions to the primal (1) and dual problem (6), respectively. Let u h be the discontinuous Galerkin solution of order k to the discrete problem [START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF]. Then, there exists a constant C, that depends on Ω, the shape regularity of T h , and the approximation order k, but is independent of ε, β and h, such that

|Q(u) -Q(u h )| ≤ C ε 1 2 + β 1 2 ∞ h 1 2 2 h 2k u H k+1 p H k+1 . ( 20
)
Remark 1. Theorems 1 and 2 provide robust error estimates, optimal in both the elliptic and hyperbolic regimes; in particular,

|Q(u) -Q(u h )| = O(h 2k ) if the problem is dominated by diffusion and |Q(u) -Q(u h )| = O(h 2k+1 ) if dominated by advection.
The a priori error estimate [START_REF] Roos | Robust numerical methods for singularly perturbed differential equations[END_REF] describes the asymptotical behavior of the error in the quantity of interest as the mesh size goes to zero. Unfortunately, it also involves the unknown exact solutions to the primal and dual problems. We develop below a computable a posteriori error estimator (i.e. exempt of unknown functions or constants), which is robust with respect to the small parameter ε and is asymptotically exact. The goal-oriented error estimator can be constructed in different ways depending on how the dual solution is computed: for instance, using residual or weighted residual-type error estimators [START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF][START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF], gradient-averaging techniques [START_REF] Kuzmin | Goal-oriented a posteriori error estimates for transport problems[END_REF], computable exact bounds for goal functionals [START_REF] Parés | Computable exact bounds for linear outputs from stabilized solutions of the advection-diffusion-reaction equation[END_REF], or functional a posteriori techniques [START_REF] Nicaise | Functional a posteriori error estimates for the reactionconvection-diffusion problem[END_REF]. We propose here to extend the ideas of our previous work [START_REF] Mozolevski | Goal oriented estimates for elliptic equations based on equilibrated fluxes reconstruction[END_REF], in which we presented an error estimator for elliptic problems based on a H(div)conforming approximation of the discrete flux using the Raviart-Thomas finite elements space. The error estimator is composed of the oscillations in the problem data, the error term between the discrete gradient of the primal solution and its conforming reconstruction in the Raviart-Thomas space, and the interface jumps in the primal solution, all weighted by the dual solution.

Asymptotically exact goal-oriented error estimator

Let us denote by H(div, Ω) the space of vector-valued functions whose divergence is well defined in Ω and by RT m (T h ) the Raviart-Thomas finite element space of order m ∈ N 0 , that is:

H(div, Ω) = {t ∈ [L 2 (Ω)] 2 : ∇ • t ∈ L 2 (Ω)}, ( 21
)
RT m (T h ) = {t ∈ H(div, Ω) : t| T ∈ [P m (T )] 2 + x P m (T ), ∀T ∈ T h }. ( 22
)
We derive a new error representation in the goal functional following the method introduced in [START_REF] Mozolevski | Goal oriented estimates for elliptic equations based on equilibrated fluxes reconstruction[END_REF] for elliptic boundary-value problems. Since

E Q (u h ) = B(u, p) - Ω u h q = Ω f p - Ω (∇ • (-ε∇p) -β • ∇p)u h ,
we obtain after integration by parts:

E Q (u h ) = Ω f p - T ∈T h T (ε∇p • ∇ h u h ) - ∂T (ε∇p • n)u h + Ω (β • ∇p)u h = Ω f p + T h σ h (u h ) • ∇p + F h (ε∇p • n)[[u h ]].
Note that for any t ∈ H(div, Ω), we have

E Q (u h ) = T h (f -∇ • t)p + T h (σ h (u h ) -t) • ∇p + F h (n • ε∇p)[[u h ]], (23)
where σ h (u h ) = -ε∇ h u h + βu h denotes the total discrete flux for the dG solution u h . A vector σh (u h ) ∈ RT k-1 (T h ) is said to be a reconstructed equilibrated flux for the discontinuous Galerkin solution

u h ∈ V k h to Problem (14) if ∇ • σh (u h ) = π k-1 h (f ), (24)
where

π k-1 h : L 2 (Ω) → V k-1 h
denotes the L 2 -orthogonal projection operator. So choosing t = σh (u h ) in [START_REF] Verfürth | A posteriori error estimation techniques for finite element methods[END_REF] and using a higher-order approximation ph ∈ V m h , m > k of the exact solution p to the dual problem, we obtain the error representation

E Q (u h ) = T h (f -π k-1 h (f ))p h + T h (σ h (u h ) -σh (u h )) • ∇ h ph + F h {{n • ε∇p h }}[[u h ]] + ∆(u, p; u h , ph ), (25)
where the remainder term ∆(u, p; u h , ph ) is given by

∆(u, p; u h , ph ) = T h (f -π k-1 h (f ))(p -ph ) - T h σ h (u h ) -σh (u h ) ∇p -∇ h ph + F h n • ε∇p -{{n • ε∇ h ph }} [[u h ]]. (26)
Note that the exact diffusive flux has been replaced along the interior faces F ∈ F i,h by the mean flux computed from adjacent elements. We will show below that the remainder term (26) is indeed of higher order when the mesh size tends to zero. We therefore identify the following computable error estimator from the error representation (25):

η(u h , ph ) = T ∈T h T (f -π k-1 h (f ))p h + T ∈T h T (σ h (u h ) -σh (u h )) • ∇p h + F h {{n • ε∇p h }}[[u h ]]. (27)
It remains to propose a method for the construction of the equilibrated fluxes σh (u h ) ∈ RT k-1 with respect to the dG solution u h ∈ V k h . The idea is to compute the degrees of freedom for σh (u h ) from the numerical fluxes, as follows (see e.g. [START_REF] Ern | Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF][START_REF] Ern | A posteriori energy-norm error estimates for advectiondiffusion equations approximated by weighted interior penalty methods[END_REF])

For k ≥ 1, ∀F ∈ F h , ∀q h ∈ P k-1 (F ), F σh (u h ) • n q h = F ({{σ h (u h )}} + (δ ε + δ β )[[u h ]]) q h ; For k ≥ 1, ∀F ∈ F Γ,h , ∀q h ∈ P k-1 (F ), F σh (u h ) • n q h = F -(ε∇ h u h ) • n + (δ ε + (β • n) + )u h q h ; For k ≥ 2, ∀T ∈ T h , ∀r h ∈ [P k-2 (T )] d , T σh (u h ) • r h = T σ h (u h ) • r h + ∂T χ F h (ε r h • n)[[u h ]]. (28)
Here χ F h represents the face indicator function on F h , that is χ

F h | F = 1/2, if F ∈ F i,h , and χ F h | F = 1, if F ∈ F Γ,h . Lemma 1. Let u h ∈ V k h , k ∈ N,
be the solution to the discontinuous Galerkin problem (14) and let σh (u h ) ∈ RT k-1 be defined by (28). Then σh (u h ) is the reconstructed equilibrated flux of u h .

Proof: It is sufficient to check that, owing to ( 15) and (28),

T (∇ • σh (u h ))w = - T ( σh (u h )) • ∇w + ∂T n • ( σh (u h ))w = B h (u h , χ T w) = T f w ∀w ∈ P k-1 , which implies that σh (u h ) = π k-1 h (f ).
In order to study the asymptotic behavior of the remainder term (26), it is necessary to estimate the error between the equilibrated flux and the discrete gradient of the dG solution. We thus require the following auxiliary result. Lemma 2. For any t ∈ RT k (T ), T ∈ T h , k ∈ N 0 , it holds:

t [L 2 (T )] d ≤ C sup r ∈ [P k-1 (T )] d r [L 2 (T )] d =1 T t • r + h 1 2 T F ∈∂T sup q ∈ P k (F ) q L 2 (F ) = 1 F (t • n)q (29)
where constant C only depends on k, d, and the shape regularity of T h .

The proof is standard and is based on the equivalence of norms on finite-dimensional spaces, the Piola transformation, and scaling arguments.

Lemma 3. Let Assumptions 1 and 2 hold and let u h ∈ V k h . Then there exists a constant C > 0, that depends on Ω, the shape regularity of T h , and k, but is independent of ε, β, and h, such that

h T π k h (f ) + ε∆u h -β • ∇u h L 2 (T ) ≤ C h T f -π k h (f ) L 2 (T ) + ε∇(u -u h ) L 2 (T ) + β L ∞ (T ) u -u h L 2 (T ) , ∀T ∈ T h .
Lemma 4. Let Assumptions 1 and 2 hold and let u h ∈ V k h . Then there exists a constant C > 0, that depends on Ω, the shape regularity of T h , and k, but is independent of ε, β, and h, such that

h 1 2 F (-ε[[∇u h ]] + β[[u h ]]) • n L 2 (F ) ≤ C h T f -π k h (f ) L 2 (ω F ) + ε∇(u -u h ) L 2 (ω f ) + β L ∞ (ω F ) u -u h L 2 (ω F ) , ∀F ∈ F h ,
where

ω F = ∪ F ∈∂T ′ T ′ .
The proofs of these lemmas make use of appropriate bubble functions and Green's formula and are similar to the proofs for standard residual error estimates [START_REF] Verfürth | A posteriori error estimation techniques for finite element methods[END_REF].

For an arbitrary T ∈ T h , we introduce the norm of v on ω T = ∪ ∂T ∩∂T ′ =∅ T ′ as:

|v | 2 ω T = ω T ε|∇ h v| 2 + β L ∞ (ω T ) L v 2 + ∂T ε h F h + |β • n| [[v]] 2 + h T ε({{∇ h }}) 2 . ( 30
)
Lemma 5. Let Assumptions 1 and 2 hold and let u h ∈ V k h be the dG solution to [START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF]. Then there exists a constant C > 0, that depends on Ω, the shape regularity of T h , and the approximation order k, but is independent of ε, β, and h, such that

σ h (u h ) -σh (u h ) [L 2 (T )] d ≤ C h T f -π k h (f ) L 2 (ω F ) + ε 1 2 + L 1 2 β 1 2 L ∞ (ω T ) |u -u h | ω T .
Proof: Owing to Lemma 2, it suffices to estimate the two terms in (29). From (28), one obtains for a given F ∈ F i,h :

(σ h (u h ) -σh (u h )) • n| F = (-ε[[∇u h ]] + β[[u h ]]) • n| F -(δ ε + δ β )[[u h ]]| F .
Then, invoking Lemma 4 and using Equation (30), we have for any q ∈ P k (F ), q L 2 (F ) = 1, that

h 1 2 T F ((σ h (u h ) -σh (u h )) • n)q ≤ h 1 2 T (-ε[[∇u h ]] + β[[u h ]]) • n L 2 (F ) + h 1 2 T (δ ε + δ β )[[u h ]] L 2 (F ) ≤ C h T f -π k h (f ) L 2 (ω F ) + ε∇(u -u h ) L 2 (ω f ) + β L ∞ (ω F ) u -u h L 2 (ω F ) + h 1 2 T δεh -1 f [[u h ]] L 2 (F ) + 1 2 h 1 2 T |β • n|[[u h ]] L 2 (F ) ≤ C h T f -π k h (f ) L 2 (ω F ) + ε 1 2 |u -u h | ω T + L 1 2 β 1 2 L ∞ (ω T ) |u -u h | ω T + ε 1 2 |u -u h | ω T + h 1 2 T β 1 2 L ∞ (ω T ) |u -u h | ω T ≤ C h T f -π k h (f ) L 2 (ω F ) + ε 1 2 + L 1 2 β 1 2 L ∞ (ω T ) |u -u h | ω T .
Again, from (28), we have for

F ∈ F Γ,h (σ h (u h ) -σh (u h )) • n| F = -(δ ε + δ β )u h | F .
Since the latter expression can be estimated in a similar manner as before, one gets the desired estimate for the first term.

For the second term, we simply have:

T (σ h (u h ) -σh (u h )) • r = ∂T χ F h (ε r h • n)[[u -u h ]] ≤ C ε(u -u h ) L 2 (∂T ) r • n L 2 (∂T ) ≤ Ch -1 2 T ε(u -u h ) L 2 (∂T ) ≤ Cε 1 2 |u -u h | ω T ,
where we have made use of the inverse inequality:

r • n L 2 (∂T ) ≤ Ch -1 2 T r [L 2 (T )] d .
This completes the proof.

Definition 1. Let the effectivity index of error estimator η(u h , ph ) be defined as:

I η = I η (u h , ph ) = η(u h , ph ) E Q (u h ) .
The error estimator η = η(u h , ph ) is said to be asymptotically exact if

lim h→0 + I η (u h , ph ) = 1.
According to above definition, it follows from (25) that the error estimator η(u h , ph ) is asymptotically exact if and only if lim

h→0 + ∆(u, p; u h , ph ) E Q (u h ) = 0.
In order to simplify the analysis of the asymptotical behavior of the error estimator (27), we suppose that f ∈ V k-1 h .

Theorem 3. Let Assumptions 1 and 2 hold and assume that f ∈ V k-1 h . Let u ∈ H k+1 (T h ), k ≥ 1, and p ∈ H m+1 (T h ), m > k be the solutions to the primal problem (1) and the dual problem (6), respectively. Let u h ∈ V k h and ph ∈ V m h be the discontinuous Galerkin solutions to the discrete primal and dual problems. Moreover, let the goal functional be properly of order 2k, i.e.

|Q(u) -

Q(u h )| ≥ c ε 1 2 + β 1 2 ∞ h 1 2 2 h 2k u H k+1 p H k+1 (31)
where constant c > 0 is independent of ε, β, and h. Then there exists a constant C that depends on domain Ω, the shape regularity of T h , and the approximation order k, but is independent of ε, β, and h, such that

∆(u, p; u h , ph ) Q(u) -Q(u h ) ≤ C 1 + Pe 1 2 h m-k p H m+1 p H k+1 . ( 32
) Proof: For f ∈ V k-1 h , Equation (26) reduces to ∆(u, p; u h , ph ) = T h σ h (u h ) -σh (u h ) ∇p -∇ h ph + F h εn • (∇p -{{∇ h ph }} [[u h ]] := ∆ 1 + ∆ 2 .
Then, from Lemma 5 and the a priori error estimate for the discrete solution to the dual problem (see Theorem 1), we obtain

|∆ 1 | = T h ε -1 2 (σ h (u h ) -σh (u h )) ε 1 2 ∇p -∇ h ph ≤ ε -1 2 |σ h (u h ) -σh (u h )| L 2 (T h ) ε 1 2 |∇p -∇ h ph | L 2 (T h ) ≤ Cε -1 2 ε 1 2 + L 1 2 β 1 2 ∞ |u -u h | |p -ph | ≤ C 1 + Pe 1 2 ε 1 2 + β 1 2 ∞ h 1 2 2 h k+m u H k+1 p H m+1 .
Similarly using the a priori error estimates obtain

|∆ 2 | = F h (εh) 1 2 n • (∇p -{{∇ h ph }} ε h 1 2 [[u -u h ]] ≤ (εh) 1 2 n • (∇p -{{∇ h ph }} L 2 (F h ) ε h 1 2 [[u -u h ]] L 2 (F h ) ≤ |u -u h | |p -ph | ≤ C ε 1 2 + β 1 2 ∞ h 1 2 2 h k+m u H k+1 p H m+1 .
Combining these estimates with (31) completes the proof.

Corollary 1. Under the assumptions of Theorem 3, the error estimator (27) is asymptotically exact if the order of the discontinuous Galerkin approximation of the dual problem is greater than the order of the discontinuous Galerkin approximation of the primal problem.

Remark 2. In the elliptic regime, i.e. Pe ∼ 1, the estimate (32) coincides with that obtained in [START_REF] Mozolevski | Goal oriented estimates for elliptic equations based on equilibrated fluxes reconstruction[END_REF]. The property that the estimator is asymptotically exact is therefore compatible with the results of [START_REF] Mozolevski | Goal oriented estimates for elliptic equations based on equilibrated fluxes reconstruction[END_REF].

Remark 3. In the advection dominated regime, i.e. P e ≫ 1, the estimate (32) shows that the effectivity index may deteriorate for coarse meshes. Such an issue can be resolved, for example, by increasing the approximation order m of the dual problem m = k + 1 + l where h l ∼ Pe -1/2 , see numerical examples below.

Numerical Experiments

Test case 1: smooth primal and dual solutions

In the first test case, we examine the robustness of the error estimator (25) with respect to the Péclet number when the primal and the dual problems possess smooth solutions. We consider, as in [START_REF] Ern | Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF], Problem [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF] in Ω = [0, 1] 2 with right-hand side term f that gives the exact solution u = x(1 -x)y(1 -y) tanh(10 -20x). We suppose here that we are interested in the goal functional Q(u) = Ω uq, where q = ∇ • (-ε∇pβp). In this case, p corresponds to the solution to the dual problem [START_REF] Eriksson | Introduction to adaptive methods for differential equations[END_REF] and is chosen in this example as p = x(1 -x)y(1 -y). For any T ∈ T h , we introduce the local error estimator η(T ) and define the three contributions as

η(T ) = T (f -π k-1 h (f ))p h + T [σ h (u h ) -σh (u h )] • ∇p h + ∂T χ F h {{n • ε∇p h }}[[u h ]] := η O (T ) + η ∇ (T ) + η H (T ), (33)
so that the estimator η (27) can be rewritten as: Here, contribution η O (T ) measures the data oscillation, contribution η ∇ (T ) estimates the deviation of the total discrete gradient from the reconstructed gradient in element T , and contribution η H (T ) measures the jump in the approximate solution along the boundary of element T . We compute the error estimator η, the contributions η O , η ∇ , η H , and the effectivity index I η for several values of the Péclet number on a sequence of uniformly refined meshes. We set β = (1, 0) and choose the following values of the viscosity ǫ = 1, 10 -2 , 10 -4 , and 10 -6 . We use the discontinuous Galerkin finite element method of order one to approximate the primal problem and of order two for the approximation of the dual problem. We also estimate the rate of convergence O(Q) in the goal functional Q.

η = T ∈T h η(T ) = T ∈T h η O (T ) + T ∈T h η ∇ (T ) + T ∈T h η H (T ) = η O + η ∇ + η H .
#T h Q(u) -Q(u h ) I η (u h , ph ) O(Q(u) -Q(u h )) 128 1.962e
As one can observe in Tables 12345678, the convergence rate in the goal functional is of order two in the elliptic regime and of order three in the hyperbolic regime, as predicted by Theorem 3. The error estimator (27) is deemed robust as the effectivity index approaches unity (e.g. the estimator is asymptotically exact) for all considered Péclet numbers in both the diffusion and advective-dominated regimes, which confirms the theoretical results of Corollary 1. Note that all contributions to the error estimator has the same order in the elliptic regime, whereas in the hyperbolic regime, the jumps contribution is of lesser order, possibly due to the up-winding stabilization.

Test Case 2: boundary layer

The main challenge in mesh adaptation strategies for convection-diffusion problems is to be able to capture thin boundary layers. Some techniques based on the theory of singularly Table 3: Test case 1. Error estimator and contributions, for ε = 10 -2 , calculated using the dG method of order k = 1 for the primal problem and m = 2 for the dual problem. 4: Test case 1. Effectivity index and rate of convergence, for ε = 10 -2 , calculated using the dG method of order k = 1 for the primal problem and m = 2 for the dual problem. perturbed convection-diffusion equations have been proposed in the past and are known as layer-adapted mesh refinement, see e.g. [START_REF] Roos | Robust numerical methods for singularly perturbed differential equations[END_REF]. However, these techniques require a priori knowledge of the solution within the boundary layer with respect to the small parameter, which may be an issue in practical applications. Moreover, in the case of goal-oriented adaptivity, the problem becomes more complex as a boundary layer may also appear in the dual solution, which implies that one has to deal with layer-adapted meshes simultaneously for both the primal and dual problems.

#T h η O (u h , ph ) η ∇ (u h , ph ) η H (u h , ph ) η(u h ,
#T h Q(u) -Q(u h ) I η (u h , ph ) O(Q(u) -Q(u h )) 128 3.
In this section, we present some numerical experiments regarding the performance of the goal-oriented error estimator in the presence of boundary layers. Note that the error estimator (27) does not require any prior knowledge about the position and/or the singularity of the solution within the boundary layer. Instead, it measures the respective error components on each element, uses this information to detect boundary layers in the domain and to refine the mesh according to the regularity of the solution in the boundary layers.

We consider Problem [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF] in Ω = (0, 1) 2 with β = (0, 1) and f chosen such that the exact solution to the problem is (see [START_REF] Zarin | Interior penalty discontinuous approximations of convectiondiffusion problems with parabolic layers[END_REF])

u(x, y) = x(1 -e -(1-x)/ √ ε )(1 -e -y/ √ ε )(1 -e -(1-y)/ √ ε ).
The goal functional is defined here as the average value of u in Ω, i.e. Q(u) = Ω u, so that the representer of the quantity of interest is given by q = 1. The solutions to the primal and dual problems both exhibit boundary layers: near the outflow boundary, Γ o = {(x, y) ∈ Γ; x = 1, 0 < y < 1} and near the lateral boundaries, Γ w = {(x, y) ∈ Γ; y = 0 or y = Effectivity index and rate of convergence, for ε = 10 -4 , calculated using the dG method of order k = 1 for the primal problem and m = 2 for the dual problem. We apply the adaptive procedure starting from an initial mesh of 128 elements and using the maximal marking strategy to select the elements to be refined and the longest edge bisection method to refine the triangular elements. In the maximal marking approach, an element T ∈ T h such that |η(T )| ≥ θ max

#T

h Q(u) -Q(u h ) I η (u h , ph ) O(Q(u) -Q(u h )) 128
T ∈T h |η(T )|
is marked as a candidate for refinement. Elements are further checked and possibly marked for refinement in order to eliminate hanging nodes. Figures 1,2, and 3 show the convergence history of the quantity of interest, as well as the adapted mesh when the relative error has reached 10 -2 , for ε = 10 -2 , 10 -3 , and 10 -4 , respectively. The primal problem was approximated using a first-order dG method and the reconstructed fluxes were obtained in the lowest order Raviart-Thomas space. The dual problem was approximated by a dG finite element method of order 2, 3, and 4, respectively, in order to ensure that the estimator remains accurate in all three simulations. We observe in all cases that the adaptive algorithm based on the proposed error estimator clearly detects the boundary layers and efficiently refines the mesh according to their thickness. We therefore conclude that the approximations of the goal functional converge to the exact value in a robust manner. Effectivity index and rate of convergence, for ε = 10 -6 , calculated using the dG method of order k = 1 for the primal problem and m = 2 for the dual problem.

#T

h Q(u) -Q(u h ) I η (u h , ph ) O(Q(u) -Q(u h )) 128

Test Case 3

To demonstrate the potential of the method, we consider the slightly more complex problem presented in [START_REF] Sauer-Budge | Computing bounds for linear functionals of exact weak solutions to the advection-diffusion-reaction equation[END_REF]. The boundary-value problem is given by (1) in the unit square Ω = (0, 1) 2 with ǫ = 2 • 10 -3 , β = (y -1/2, 1/2 -x), f = 2χ [0.750,0.875] 2 , and homogeneous Dirichlet boundary conditions. We also consider two quantities of interest defined as follows: 1) the average value of the solution over the square region ω = [0.125, 0.250] 2 , that is

Q 1 (u) = |ω| -1
ω u; 2) the average value of the solution over the whole domain Ω, that is

Q 2 (u) = Ω u.
The adaptive process starts here from a uniform mesh with 128 elements and utilizes the error estimator (25) calculated from the first-order dG method and reconstructed fluxes based on minimal order Raviart-Thomas finite elements for the primal problem and the second-order dG method for the dual problem. Mesh refinement relies on the bisection method and the maximal marking strategy with θ = 0.75.

Figure 4 shows the convergence history of the goal functional Q 1 and the mesh of 1,687 elements obtained with a 1% relative error. The exact value of the output functional was calculated using an overkill uniform mesh of 131,072 elements. Figure 5 displays the solution to the primal and dual problems on the adaptive mesh that corresponds to the relative error of 1%. We report in Figures 6 and7 the same results for goal functional Q 2 . The adaptive mesh obtained with a relative error of 1% has 949 elements.

We observe that the approximation of the goal functional in both cases converges to the exact value. Convergence of Q 1 is however slower than that of Q 2 since Q 1 has a smaller support, which implies that the corresponding dual solution is less regular than for Q 2 . It also results in a more refined final mesh than in the case of Q 2 . Nonetheless, the approximation of the quantities of interest is achieved on relatively coarse meshes and the adaptive process allows one to avoid solving for globally accurate primal and dual solutions, see Figures 5 and7.

Conclusions

We have developed a goal-oriented error estimator for high-order discontinuous Galerkin approximations of convection-diffusion boundary-value problems. We have proven that the estimator is asymptotically exact in both the elliptic and hyperbolic regimes and robust with respect to the Péclet number of the problem. More precisely, we have shown that if the approximation order m of the dual problem is greater than the approximation order k of the primal problem, i.e. m ≥ k +1, the effectivity index behaves as (1+Pe 1/2 )o(h m-k ). This implies that the effectivity index may deteriorate on coarse meshes for large Péclet numbers. In order to circumvent this issue, we have suggested to increase the approximation order of the dual problem and have confirmed the robustness of such a strategy on several numerical examples. The proposed goal-oriented error estimator is based on the construction of equilibrated fluxes. It is expected that it can be extended to heterogeneous and non-linear convection-diffusion problems important for practical applications (e.g. multiphase flow and transport modeling in heterogeneous porous media) and this will be the subject of future work.

Figure 1 :

 1 Figure 1: Test case 2. ǫ = 10 -2 . Convergence history and adaptively refined mesh corresponding to 1% relative error for k = 1, kRT = 0, and m = 2.

Figure 2 :

 2 Figure 2: Test case 2. ǫ = 10 -3 . Convergence history and adaptively refined mesh corresponding to 1% relative error for k = 1, kRT = 0, and m = 3.

Figure 3 :

 3 Figure 3: Test case 2. ǫ = 10 -4 . Convergence history and adaptively refined mesh corresponding to 1% relative error for k = 1, kRT = 0, and m = 4.

Figure 4 :

 4 Figure 4: Test case 3. Convergence history of the goal functional Q1 (left) and mesh corresponding to 1% relative error in the approximation of the goal functional (right).

Table 1 :

 1 Test case 1. Error estimator and contributions, for ε = 1, calculated using the dG method of order k = 1 for the primal problem and m = 2 for the dual problem.

	#T h 128	η O (u h , ph ) η ∇ (u h , ph ) η H (u h , ph ) η(u h , ph ) 1.129e-02 2.517e-02 -1.355e-02 2.291e-02
	512	2.763e-03	6.882e-03 -3.587e-03 6.058e-03
	2048	6.909e-04	1.771e-03 -9.616e-04 1.501e-03
	8192	1.728e-04	4.482e-04 -2.514e-04 3.697e-04
	32768	4.322e-05	1.127e-04 -6.353e-05 9.236e-05
	131072 1.081e-05	2.824e-05 -1.593e-05 2.312e-05

Table 2 :

 2 Test case 1. Effectivity index and rate of convergence, for ε = 1, calculated using the dG method of order k = 1 for the primal problem and m = 2 for the dual problem.

Table

		ph )
	128	3.789e-02 -3.698e-02 2.152e-03 3.058e-03
	512	8.588e-03 -8.563e-03 7.349e-04 7.599e-04
	2048	2.133e-03 -2.202e-03 2.370e-04 1.679e-04
	8192	5.335e-04 -5.601e-04 6.649e-05 3.987e-05
	32768	1.334e-04 -1.411e-04 1.756e-05 9.866e-06
	131072 3.335e-05 -3.539e-05 4.506e-06 2.466e-06

Table 5 :

 5 Test case 1. Error estimator and contributions, for ε = 10 -4 , calculated using the dG method of order k = 1 for the primal problem and m = 2 for the dual problem.

	#T h 128	η O (u h , ph ) η ∇ (u h , ph ) η H (u h , ph ) η(u h , ph ) 4.104e-02 -4.008e-02 4.764e-05 1.010e-03
	512	9.281e-03 -9.073e-03 2.233e-05 2.301e-04
	2048	2.303e-03 -2.285e-03 1.177e-05 3.010e-05
	8192	5.762e-04 -5.785e-04 6.323e-06 3.996e-06
	32768	1.441e-04 -1.467e-04 3.179e-06 5.685e-07
	131072 3.602e-05 -3.744e-05 1.509e-06 8.715e-08

Table 6 :

 6 Test case 1.

Table 7 :

 7 Test case 1. Error estimator and contributions, for ε = 10 -6 , calculated using the dG method of order k = 1 for the primal problem and m = 2 for the dual problem.

	#T h 128	η O (u h , ph ) η ∇ (u h , ph ) η H (u h , ph ) η(u h , ph ) 4.108e-02 -4.009e-02 5.207e-07 9.913e-04
	512	9.289e-03 -9.064e-03 2.649e-07 2.246e-04
	2048	2.305e-03 -2.277e-03 1.463e-07 2.845e-05
	8192	5.767e-04 -5.732e-04 7.573e-08 3.505e-06
	32768	1.442e-04 -1.438e-04 3.623e-08 4.352e-07
	131072 3.605e-05 -3.601e-05 1.677e-08 5.431e-08

Table 8 :

 8 Test case 1.

Acknowledgements. SP is grateful for the support by a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada. He is also a participant of the KAUST SRI Center for Uncertainty Quantification in Computational Science and Engineer-