c
° TU Delft, The Netherlands, 2006
RESIDUAL BASED A POSTERIORI ERROR ANALYSIS
FOR THE ADVECTION-DIFFUSION EQUATION WITH
ANISOTROPIC AND DISCONTINUOUS DIFFUSION
TENSOR
Alexandre Ern∗ and Annette F. Stephansen†
ENPC, CERMICS
6 et 8 avenue Blaise Pascal, 77455 Marne la Vall´ee Cedex 2, France
∗e-mail: ern@cermics.enpc.fr †e-mail: stephansen@cermics.enpc.fr
Key words: Discontinuous Galerkin, a posteriori error analysis, mesh adaptation Abstract. We present a residual a posteriori error estimate for the anisotropic advection
diffusion equation with continuous or discontinuous diffusion tensor. Our numerical re-sults show that the adapted mesh based on the residual is more refined in the region where anisotropy and heterogeneity effects create difficulties for the numerical solution.
1 INTRODUCTION
The Discontinuous Galerkin (DG) scheme has triggered considerable interest in the last few years. Its popularity stems in part from the fact that it can be applied to non– matching grids and its code has proven efficient with parallelization. The scheme is also particularly useful when the solution presents discontinuities, as is often the case for non– linear hyperbolic problems and for linear problems with heterogeneous data. For elliptic problems a unified analysis for DG schemes has been proposed by Arnold et al.1
A posteriori errors based on residual estimates can be tracked back to the 1970s, when it was considered by Babuˇska and Rheinboldt2,3. The residual estimates are easy and
economic to calculate directly from the computed solution, and control the error by means of the residues. When used for mesh adaptation one tries to balance the residual on each element. We will show that the error is controlled by two terms, namely the residue itself and a space conforming error.
been studied by Ern and Guermond4,5and Houston et al6. For the problem with
discon-tinuous diffusion tensor we have proved the well-posedness of the problem separately.
2 THE ADVECTION-DIFFUSION EQUATION
We consider the following advection-diffusion equation with homogeneous Dirichlet boundary conditions:
(
−∇·(K∇u) + β·∇u + µu = f in Ω
u = 0 on ∂Ω (1)
Here µ ∈ L∞(Ω), β ∈ [L∞(Ω)]d, ∇·β ∈ L∞(Ω), the diffusion tensor K is a symmetric,
uniformly positive definite field in [L∞(Ω)]d,d and f ∈ L2(Ω).
We formulate the advection-diffusion equation as a a system of partial differential equations of first-order, with the following weak formulation:
Seek (σ, u) ∈ V such that (
(σ, τ )0,Ω− (κ∇u, τ )0,Ω = 0 ∀τ ∈ [L2(Ω)]d
−(∇·(κσ), v)0,Ω+ (β·∇u, v)0,Ω+ (µu, v)0,Ω = (f, v)0,Ω ∀v ∈ L2(Ω)
(2) The tensor κ is defined such that κκ = K, and the subscript 0 refers to the L2 scalar
product.
In the case of a continuous diffusion coefficient can be seen as a Friedrichs’s system by setting ˜β = β − ∇·κ, and considering
(
(σ, τ )0,Ω− (κ∇u, τ )0,Ω = 0
−(κ:∇σ), v)0,Ω+ ( ˜β∇u, v)0,Ω+ (µu, v)0,Ω = (f, v)0,Ω
(3) The solution is sought in the graph space V = Hσ(Ω) × H01(Ω) where Hσ(Ω) = {σ ∈
[L2(Ω)]d,d, κ:∇σ ∈ L2(Ω)}, and the problem is well-posed with respect to the norm k(σ, u)kV = k(σ, u)k0,Ω+ k∇uk0,Ω+ kκ:∇σk0,Ω (4)
When the diffusion coefficient K is discontinuous, the solution is sought in the space
V = H(κ div; Ω) × H1
0(Ω), where H(κ div; Ω) = {σ ∈ [L2(Ω)]d, ∇·(κσ) ∈ L2(Ω)}. The
problem is well-posed with respect to the norm
k(σ, u)k2
V = k(σ, u)k0,Ω+ kβ·∇u − ∇·(κσ)k20,Ω+ kκ∇uk20,Ω (5)
3 RESIDUAL A POSTERIORI ERROR ANALYSIS
In the following the symbol . indicates an inequality involving a constant C indepen-dent of triangle-size. We define the broken norm as:
kvkV,h = Ã X T ∈Th kvk2 V,T !1 2 (6)
3.1 Global upper bound
The following a posteriori error estimate can be obtained:
k(σ − σh, u − uh)kV,h. Ã X T ∈Th R2T !1/2 + inf (τ,v)∈V Ã X T ∈Th PT2 !1/2 (7) For the Friedrichs’s system the local error indicators are
R2T = kf + κ:∇σh− ˜β∇uh− µuhk20,T + kκ∇uh− σhk20,T (8)
PT2 = k(σh− τ, uh− v)k2V,T (9)
while for a discontinuous diffusion tensor we have the local error indicators
R2T = kf + ∇·(κtσh) − β·∇uh− µuhk20,T + kκ∇uh− σhk20,T (10)
PT2 = k(σh− τ, uh− v)k2V,T (11)
We see that the local error indicator can be divided into two parts; the first being the residue itself, while the other measures the distance between the discontinuous approxi-mation space and the continuous solution space, and is referred to as a conforming error. To calculate the conforming error, we make use of the Oswald interpolator, which maps any function in Vh to a continuous function:
IOs : Vh(Th) → P (Th) ≡ {wh ∈ V ∩ C0(Ω); wh|T ∈ [Pp]m}
It is defined at each node as the mean of the function taken over the elements that contain that node. On each T ∈ Th we introduce the elemental Lagrange interpolation nodes s.
At interior nodes, the Oswald interpolation operator is thus defined as
IOsvh(s) = 1 card(Ts) X T ∈Ts vh|T(s) (12)
3.2 Local lower bound
The efficiency of an a posteriori error estimate relies on being able to obtain a lower bound of the estimate locally in terms of the error measure. The error measure and the a posteriori error estimate are then equivalent norms. The residual part of our error indicator can be bound locally by the error in the V, T norm, but the conforming part presents more of a problem. The conforming error can be bound locally by the H1(T )
seminorm of the error (which is not in the V, T norm), and thus the estimate is not optimal.
3.3 An alternative residual estimate
An alternative residual estimate can be obtained by considering the error of (IOsσh, uh).
In this case the conforming error with respect to IOsσh is zero, and we have the following
estimate: k(σ − IOsσh, u − uh)kV,h . Ã X T ∈Th b R2 T !1/2 + inf (τ,v)∈V Ã X T ∈Th b P2 T !1/2 (13) The residue bR2
T is calculated as in (8) and (10) substituting σh with IOsσh, while the
space conforming error is b
P2
T = kuh− vk20,T + kβ·∇(uh− v)k20,T + kκ∇(uh− v)k20,T (14)
4 NUMERICAL RESULTS
We consider the domain, the unit square, split into two regions, where we indicate the bottom left-hand square with region 1 and the remaining part region 2. We consider
µ = 1, β = (−1, 0)t, while the diffusion tensor takes on different values in the two regions: K1 = µ 1 0 0 1e−3 ¶ K2 = µ 1e−3 0 0 1 ¶ (15) The forcing term considered is
µu + β·∇u − ∇·(K∇u) (16)
with
u = x6sin(πx) sin(πy) (17)
Due to the discontinuity in the diffusion tensor, u is not the exact solution. We will consider the alternative residual estimate which estimates the error of (IOsσh, uh) instead
of the error of (σh, uh). We solve the advection-diffusion equation with P1 discontinuous
Figure 1: Solution calculated on the original mesh IsoValue -0.0075304 0.097344 0.202218 0.307093 0.411967 0.516842 0.621716 0.72659 0.831465 0.936339 1.04121 1.14609 1.25096 1.35584 1.46071 1.56559 1.67046 1.77533 1.88021 1.98508 2.08996
calculated at the adapted mesh does not exhibit this instability, as can be seen in figure 2. The mesh has here been adapted three times. It is interesting to note that the mesh has been refined only along the inner right hand side boundary, which is the cause of the aforementioned instability. A mesh adaptation based on the gradients of the solution would have refined the boundary layer close to the upper inner boundary also, where however the method manages well.
Figure 2: Solution calculated on the adapted mesh
IsoValue -0.0165278 0.0940189 0.204566 0.315112 0.425659 0.536206 0.646753 0.757299 0.867846 0.978393 1.08894 1.19949 1.31003 1.42058 1.53113 1.64167 1.75222 1.86277 1.97331 2.08386 2.19441 REFERENCES
[1] D. Arnold, F. Brezzi, B. Cockburn, and L.D. Marini. Unified Analysis of discontinu-ous Galerkin methods for elliptic problems. SIAM J. Numer. Anal., 39, 1749–1779, 2002.
[3] I. Babuˇska and W. Rheinbolt. A posteriori error estimates for the finite element method. Int. J. Numer. Methods Engrg., 12, 1597–1615, (1978).
[4] A. Ern and J-L. Guermond. Discontinuous Galerkin methods for Friedrichs’ sym-metric systems. I. General theory. SIAM J. Numer. Anal.,, to appear.
