A posteriori error estimation of residual type for anisotropic ... - CiteSeerX
A posteriori error estimation of residual type for anisotropic diffusion–convection–reaction problems Thomas Apel∗
Serge Nicaise†
May 22, 2009 Abstract: This paper presents an a posteriori residual error estimator for diffusion– convection–reaction problems with anisotropic diffusion, approximated by a SUPG finite element method on isotropic or anisotropic meshes in Rd , d = 2 or 3. The equivalence between the energy norm of the error and the residual error estimator is proved. Numerical tests confirm the theoretical results. Key words: anisotropic diffusion, SUPG, a posteriori error estimate. AMS subject classification: 65N30, 65N15
1
Introduction
This paper is devoted to the singularly perturbed diffusion–convection–reaction problem with special focus on anisotropic diffusion: for f ∈ L2 (Ω) and g ∈ L2 (ΓN ), let u be the solution of −div (A∇u) + b · ∇u + cu = f in Ω, u = 0 on ΓD , (1) A∇u · n = g on ΓN , where the matrix A and the functions b and c satisfy assumptions (A1) to (A6) below, and Ω ⊂ Rd , d = 2 or 3, is a bounded domain with a polygonal (d = 2) or polyhedral (d = 3) boundary Γ. This boundary is divided into two parts ΓD and ΓN , where Dirichlet and Neumann boundary conditions are imposed, respectively. ∗
Institut f¨ ur Mathematik und Bauinformatik, Universit¨at der Bundeswehr M¨ unchen, D–85577 Neubiberg, Germany, [email protected] † Universit´e de Valenciennes et du Hainaut Cambr´esis, LAMAV, ISTV, F–59313 - Valenciennes Cedex 9, France, [email protected]
1
We are particularly interested in the case when A tion, for instance the cases µ ¶ ε ε 0 0 A= (d = 2), or A = 0 1 0
becomes small in some direc 0 0 1 0 0 1
(d = 3),
ε > 0. In the case when ε is small with respect to b and c, the problem is singularly perturbed and the solution may generate sharp boundary or interior layers, where the solution of the limit problem (corresponding to ε = 0) is not smooth or does not satisfy the boundary condition. Let us quote [17, 18, 19] for the a priori error analysis in two dimensions. It is shown that anisotropic finite elements must be used in order to achieve convergence uniform in the perturbation parameter ε. There is a vast amount of literature on a posteriori error estimation. For singularly perturbed problems with convection we cite [2, 8, 11, 16, 21, 25, 26], where anisotropic finite element meshes were considered in [8, 16, 21] only. An anisotropic diffusion tensor is considered only in [7]. In this paper we combine all those ingredients and derive a residual type error estimator. We prove the reliability and efficiency of this error estimator where the dependence on ε is traced. The lower bound mainly depends on the local mesh Peclet number PeT := hmin,A,T kA−1/2 bk∞,T , therefore the efficiency is achieved if PeT ≤ c which is always satisfied in the absence of convection. The reliability is based on the introduction of an alignment measure as it was done by Kunert [12, 13]. The quantity is of the order one if the mesh is well adapted to the problem, see the discussion in Subsection 3.3. Let us mention that, to our knowledge, no approach is known that leads to twosides estimates on anisotropic meshes without any assumption on the mesh. The classical results as summarized in [1, 24] are obtained for isotropic meshes only. The dual weighted residual method, see [5] for an overview, is applied in [8, 9] on anisotropic meshes, but there is no estimate from below. The more recent approach in [20] is not yet analyzed for anisotropic meshes and two different error estimators are used for the upper and lower bounds. Let us finally mention the approach by Picasso [21] who considers anisotropic meshes and proves reliability for an estimator that depends on ∇(u − uh ) where ∇u is replaced in practice by a recovered gradient ∇R u. We note that we can control in the same way the alignment measure. In this paper we develop an estimator of residual type for problems with convection, reaction and anisotropic diffusion. For the discretization we use the h-version of the streamline upwind Petrov–Galerkin method (SUPG). Without the stabilization term, the method reduces to a standard Galerkin method and produces non-physical oscillations. We note that our error estimator works as well in this case. 2
In comparison with the paper [7], where a posteriori error estimation is investigated for an isotropic discretization of a problem with anisotropic diffusion but without convection, our residual error estimator allows to prove an optimal lower bound. The factor ε−1/2 in the upper bound in [7] is retained in our analysis, since the alignment measure is of the same order in the isotropic case. Our experiments show, however, that the effectivity index is bounded uniformly in ε on adequately refined anisotropic meshes. In this sense, our analysis is sharper. The outline of the paper is as follows. In Sections 2 and 3 we introduce the discretization, notation, and estimates for bubble functions and the interpolation. The a posteriori error estimator is introduced in Section 4 where also the the upper and lower bounds are proved. The paper is completed with a numerical test in Section 5 and with conclusions. As usual, we denote by L2 (.) the Lebesgue spaces and by H s (.), s ≥ 0, the standard Sobolev spaces. The usual norm and seminorm of H s (D) are denoted by k · ks,D and | · |s,D . For the sake of brevity the L2 (D)-norm will be denoted by k · kD and in the case D = Ω, we will drop the index Ω. The space HΓ1D (Ω) is defined, as usual, by HΓ1D (Ω) := {v ∈ H 1 (Ω) : v = 0 on Γ}. In the sequel the symbol | · | will denote either the Euclidean norm in Rd , d = 2 or 3, or the length of a line segment, or the measure of a domain of Rd . Finally the notation a . b means here and below that there exists a positive constant C independent of a and b (of the mesh size of the triangulation, as well as the diffusion matrix A) such that a ≤ C b. The notation a ∼ b means that a . b and b . a hold simultaneously.
2
Discretization of the diffusion–convection–reaction equation
Let Ω be a bounded open subset of Rd , d = 2 or 3, with a polygonal (d = 2) or polyhedral (d = 3) boundary Γ. This boundary is divided into two parts ΓD and ΓN , where Dirichlet and Neumann boundary conditions are imposed, respectively. We consider the standard elliptic problem: for f ∈ L2 (Ω) and g ∈ L2 (ΓN ), let u be the solution of (1) where A, b and c satisfy the following assumptions: (A1) (A2) (A3) (A4) (A5) (A6)
b ∈ W 1,∞ (Ω)d , c ∈ L∞ (Ω), ∃c0 ≥ 0 : c − 21 div b ≥ c0 and if c0 = 0 then c ≡ 0, b · n ≥ 0 on ΓN , A ∈ Rd×d is symmetric, ∃α0 > 0 : Aξ · ξ ≥ α0 , ∀ξ ∈ Rd . ΓD 6= ∅ if c0 = 0. 3
Note that the assumption ”if c0 = 0 then c ≡ 0” is not necessary for our proofs but simplifies the presentation. Now we define the weighted H 1 (semi-)norm Z 2 (2) |||u|||ω := (|A∇u|2 + c0 |u|2 ). ω
on a subdomain ω of Ω. Let us further introduce the space HΓ1D (Ω) = {v ∈ H 1 (Ω) : u = 0 on ΓD }. and the forms Z B(u, v) = ZΩ F (v) =
(A∇u · ∇v + b · ∇uv + cuv) dx, Z f v dx + gv dΓ(x).
Ω
ΓN
For further purposes, we denote by Bω the restriction of B on a subset ω of Ω, namely Z Bω (u, v) = (A∇u · ∇v + b · ∇uv + cuv) dx. ω
With this notation, the variational formulation of problem (1) reads: Find u ∈ solution of (3) B(u, v) = F (v) ∀v ∈ HΓ1D (Ω). HΓ1D (Ω)
The assumptions (A1) to (A6) guarantee that B is continuous and coercive, i.e., B satisfies (4) (5)
B(v, v) ≥ |||v|||2
∀v ∈ HΓ1D (Ω),
¡ ¢ −1/2 |Bω (u, v)| ≤ |||u|||ω max{1, c−1 bk∞,ω kvkω , 0 kck∞,ω }|||v|||ω + kA
for all u, v ∈ H 1 (ω). If c0 = 0, then the term c−1 0 kck∞,ω disappears and the maximum in (5) is equal to one. By the Lax-Milgram lemma, problem (3) has a unique solution u ∈ HΓ1D (Ω). Note that the case div b = 0, c = 0, i. e. c0 = 0, is admitted. It is excluded in several other publications. To approximate problem (3) by a finite element scheme we fix a family {Th }h>0 of meshes of Ω that satisfies the usual conformity conditions, cf. [6, Chapter 2]. In 2D we assume that all elements of Th are triangles, while in 3D the mesh is made up of tetrahedra. For T ∈ Th we denote by hT the diameter of T , and h = maxT ∈Th hT . 4
Let Vh be the subspace of HΓ1D (Ω) defined by Vh = {vh ∈ HΓ1D (Ω) : vh|T ∈ Pk (T ) ∀T ∈ Th }, where k is a positive integer. Problem (3) is now approximated by a Streamline Upwind Petrov Galerkin scheme (SUPG): Find uh ∈ Vh such that (6)
Bh (uh , vh ) = Fh (vh ) ∀vh ∈ Vh ,
where X
Bh (uh , vh ) = B(uh , vh ) +
δT (−div (A∇uh ) + b · ∇uh + cuh , b · ∇vh )T ,
T ∈Th
Fh (vh ) = F (vh ) +
X
δT (f, b · ∇vh )T .
T ∈Th
The parameters δT ≥ 0 should satisfy similar assumptions as in [16] where the case of isotropic diffusion was investigated: hmin,A,T , kA−1/2 bk∞,T ≤ 2(1 − α)µ−2 h2min,T kAk−1 2→2 , µ ¶−1 ≤ 2(1 − α)c0 max c(x)2
(7)
δT ≤
(8)
δT
(9)
δT
x∈T
if c 6≡ 0,
for all T ∈ Th and some α ∈ (0, 1). The element quantities hmin,T and hmin,A,T are introduced below, and µ is the constant in the inverse inequality k∇ · ∇vh kT ≤ µh−1 min,T k∇vh kT . Note further that (8) and (9) guarantee the coercivity Bh (vh , vh ) ≥ α|||vh |||2 with the above α ∈ (0, 1), compare [23, Lemma 3.25] for the case of isotropic meshes. The optimal choice of δT was discussed for the slightly different Galerkin– Least-Squares method and for the case of isotropic diffusion in [4]. This choice satisfies the conditions (7)–(9). We note also that the choice δT = 0 (pure Galerkin method) satisfies these conditions. Meanwhile it is well known that this choice is suited within boundary layers if adequately refined anisotropic meshes are used there [23, p. 391 ff.]. Outside the layers, the choice δT = 0 leads in general to non-physical oscillations. Therefore this choice is not advisable, but the error estimator still works and estimates the large error.
5
3
Analytical tools
Let us define Eh as the set of edges (d = 2) or faces (d = 3) of the triangulation and let Ehint = {E ∈ Eh /E ⊂ Ω} be the set of interior edges/faces of Th , while Ehext = Eh \ Ehint is the set of boundary edges/faces of Th . For an edge/face E of a 2D/3D element T we denote by nT,E the unit outward normal vector to T along E. Furthermore we fix one of the two normal vectors of E and denote it by nE . The jump of some function v across an edge/face E at a point y ∈ E is defined as ½ limα→+0 v(y + αnE ) − v(y − αnE ) ∀E ∈ Ehint , [[v(y)]]E := v(y) ∀E ∈ Ehext . Finally we will need local subdomains, also called patches. For any T ∈ Th , let, as usual, ωT be the union of all elements having a common vertex with T . Similarly let ωE be the union of the elements having E as edge/face.
3.1
Some anisotropic quantities
As explained in the introduction, anisotropic discretizations can be very advantageous or, in certain situations, even mandatory. More information and arguments concerning anisotropy can be found in [3, 12]. Let us shortly recall some useful anisotropic quantities from Kunert [12], see also [15, 16]. We start with an arbitrary (anisotropic) tetrahedron T . We enumerate its vertices so that P0 P1 is the longest edge, meas2 (4P0 P1 P2 ) ≥ meas2 (4P0 P1 P3 ), and meas1 (P1 P2 ) ≥ meas1 (P0 P2 ). Further, we introduce three orthogonal vectors pi,T of length hi,T := |pi,T |, as described in Figure 1. P3
p3
P2 p2
p1
P0
P1 Figure 1: Notation of tetrahedron T .
6
The minimal element size is particularly important; thus we define hmin,T := h3,T . The three main anisotropic directions pi,T play an important role in several proofs. They span the matrix CT := (p1,T , p2,T , p3,T ) ∈ R3×3 . This matrix may be considered as a transformation matrix which defines implicitly a reference element TˆT via TˆT := CT−1 (T − P~0 ), cf. Figure 2. Note in particular that the reference element TˆT is of size O(1). Pˆ3
Pˆ2
Pˆ0
Pˆ1
Figure 2: Reference tetrahedron TˆT . In 2D the notation is similar. For a triangle T the enumeration is as in the bottom triangle P0 P1 P2 of Figure 1. We set hmin,T := h2,T , and CT becomes a 2 × 2 matrix. The new idea is now to transform any T ∈ Th by the matrix A−1/2 . More precisely, we transform T into TA by the affine transformation (10)
FA,T : T → TA : x → A−1/2 (x − gT ) + gT ,
where gT is the center of gravity of T . This element TA is a triangle in 2D or a tetrahedron in 3D that can be isotropic or not. Therefore we use its anisotropic quantities hi,TA , hmin,TA , CTA as introduced before.
7
Remark 3.1 In 2D, take
µ
ε 0 0 1
A=
¶ ,
c = 1, b = 0 and Ω = (0, 1)2 , ΓN = ∅. If f is smooth then u has boundary layers near x1 = 0 and x1 = 1. The transformation ˜ : x → A−1/2 x, Ω→Ω replaces the problem (1) into ˜ = (0, √1 ) × (0, 1), −∆˜ u + u˜ = f˜ in Ω ε ˜ u˜ = 0 on ∂ Ω. ˜ Then the triangle Take for simplicity a uniform triangular triangulation of Ω. ˜ with nodes (0, 0), (h, 0), and (0, h) becomes the triangle T ∈ Ω with vertices T˜ ∈ Ω √ (0, 0), ( εh, 0), and (0, h) by the inverse transformation. This element T is a good one to capture adequately the boundary layer near x1 = 0. Moreover by using FA,T , the triangle T is transformed into an isotropic element TA which is a translation of T˜, and therefore CT˜ = h−1 Id. For further use, we denote (11)
hmin,TA0 . hmin,A,T = min 0 T ⊂ωT
Note that the composition of the transformation TˆTA → TA : xˆ → CTA xˆ + P0 , −1 with FA,T , see (10), yields the following transformation from TˆTA to T :
(12)
TˆTA → T : xˆ → CA,T xˆ + bT
with (13)
CA,T = A1/2 CTA .
Note that TˆTA depends on T and A but is of unit size in the sense of Figure 2. Finally we introduce a scaling factor αT that will be used quite often: (14)
−1/2
αT = min{c0
8
, hmin,A,T }.
−1/2
Here and below, we use the convention that c0 = +∞ if c0 = 0. For an edge/face E of an element T introduce the height hE,T = We finally require, as usual, that (15) (16)
|T | . |E|
|T | ∼ |T 0 | if T ∩ T 0 6= ∅, the number of elements containing a vertex x is bounded uniformly.
Remark 3.2 If A = εId (the case treated by Kunert in [15, 16]), then TA,T is simply a homothetic transformation of T with a factor ε−1/2 . Therefore the matrix CεId,T defined by (13) is equal to CT . Moreover we have hmin,εId,T = ε−1/2 hmin,T . This last property implies that −1/2
αT = min{c0
, ε−1/2 hmin,T },
which is exactly the scaling factor introduced in [16].
3.2
Bubble functions, extension operator, and inverse inequalities
For our further analysis we require standard bubble functions and extension operators that satisfy certain properties recalled here for the sake of completeness. We need two types of bubble functions, namely bT and bE associated with an element T and an edge E, respectively. For a triangle or a tetrahedron T , denoting by λaTi , i = 1, · · · , d + 1, the barycentric coordinates of T and by aE i , i = 1, · · · , d, the vertices of the edge/face E ⊂ ∂T we recall that bT =
d+1 Y
λaTi and bE =
d Y
λaEi .
i=1
i=1
We note that bT = 0 on ∂T,
bE = 0 on ∂ωE ,
kbT k∞,T = kbE k∞,ωE ∼ 1.
In 2D, denote by T¯ the standard reference element with vertices (0, 0), (1, 0), and (0, 1). For an edge E¯ of T¯ included into the x¯1 axis, the extension Fext (vE¯ ) of ¯ to T¯ is defined by Fext (vE¯ )(¯ vE¯ ∈ C(E) x1 , x¯2 ) = vE¯ (¯ x1 ). The extension operator Fext (vE ) of vE ∈ C(E) to T for an edge E ⊂ ∂T is obtained using the affine transformation mapping T to T¯ and E to E¯ and the extension operator defined above. We proceed similarly in 3D. Now we may state the so-called inverse inequalities that are proved using classical scaling techniques, cf. [24] for the isotropic case and [12] for the anisotropic case. 9
Lemma 3.3 (Inverse inequalities) Let vT ∈ Pk0 (T ) and vE ∈ Pk1 (E), for some nonnegative integers k0 and k1 . Then the following inequalities hold, with the constants in the inequality depending on the polynomial degrees k0 or k1 but not on T , E or vT , vE . 1/2
(17)
kvT bT kT ∼ kvT kT ,
(18) (19)
kvE bE kE ∼ kvE kE , |||vT bT |||T . αT−1 kvT kT .
1/2
Proof: The equivalence (17) and (18) are proved in [12], see also Lemma 1 of [16]. For the last estimate, we write |||vT bT |||2T = c0 kvT bT k2T + kA1/2 ∇(vT bT )k2T ≤ c0 kvT k2T + kCT−T k2 kCTTA A1/2 ∇(vT bT )k2T A Now recalling that kCT−T k ∼ h−1 min,TA and using the affine transformation (12), we A obtain ˆ vT ˆbT )k2ˆ . |||vT bT |||2T . c0 kvT k2T + h−2 min,A,T |T |k∇(ˆ T Since Tˆ is regular in Ciarlet’s sense we can use the inverse inequality with hTˆ ∼ 1 to deduce that |||vT bT |||2T . c0 kvT k2T + h−2 vT k2Tˆ . min,A,T |T |kˆ Going back to T again using the affine transformation (12), we obtain (19). As usual for singularly perturbed problems, we need to use squeezed edge/face bubble functions bE,γ . Here according to our previous point of view, they are defined through the transformation FA,T from (10). Namely for a fixed edge/face E of T , the mapping (10) transforms E into an edge/face EA of TA . Now for a parameter γ ∈ (0, 1], we define the squeezed element TEA ,γ of TA as in [16]. The squeezed element TE,γ of T is simply the element obtained by the inverse transformation −1 TE,γ = FA,T TEA ,γ .
Note that TE,γ is the usual squeezed element on T with the parameter γ depending on A. For the sake of simplicity we do not write this dependence. The squeezed edge/face bubble function bE,γ is defined on the two elements T1,E,γ and T2,E,γ sharing γ, as the usual edge/face bubble function on these elements and extended by zero outside T1,E,γ ∪ T2,E,γ 10
Lemma 3.4 (Further Inverse inequalities) Under the assumptions of Lemma 3.3, we have 1/2
(20)
kbE,γ Fext (vE )kT . γ 1/2 hE,T kvE kE ,
(21)
kA1/2 ∇(bE,γ Fext (vE ))kT . γ −1/2 hE,T h−1 min,A,T kvE kE .
1/2
Proof: Scaling arguments yield ˜ext (˜ kbE,γ Fext (vE )kT = |T |1/2 |TA |−1/2 k˜bE,γ F vE )kTA , where we write v˜(˜ x) = v(x). Now using Lemma 2 of [16] in TA , we have ˜ext (˜ |TA |−1/2 k˜bE,γ F vE )kTA . γ 1/2 |EA |−1/2 k˜ vE kEA . Again scaling arguments lead to |EA |−1/2 k˜ vE kEA . |E|−1/2 kvE kE .
(22)
The three above estimates imply (20). For the second estimate, scaling arguments yield ˜ ˜bE,γ F ˜ext (˜ kA1/2 ∇(bE,γ Fext (vE ))kT = |T |1/2 |TA |−1/2 k∇( vE ))kTA . Again Lemma 2 of [16] applied in TA leads to kA1/2 ∇(bE,γ Fext (vE ))kT . |T |1/2 γ 1/2 |EA |−1/2 k˜ vE kEA min{γhEA ,TA , hmin,TA }−1 . Using the estimate (22), we arrive at 1/2
kA1/2 ∇(bE,γ Fext (vE ))kT . γ 1/2 hE,T kvE kE min{γhEA ,TA , hmin,TA }−1 . The estimate (21) will be proved if we can show that min{γhEA ,TA , hmin,TA }−1 . γ −1 h−1 min,A,T , or equivalently (23)
min{γhEA ,TA , hmin,TA } & γhmin,A,T .
But it was proved in Lemma 3.1 of [15] that hEA ,TA & hmin,TA . Since γ ∈ (0, 1] we then have γhEA ,TA & γhmin,TA ≥ γhmin,A,T and hmin,TA ≥ γhmin,A,T . This leads to (23). 11
3.3
Anisotropic interpolation error estimates
In order to obtain an accurate discrete solution uh , it is obviously helpful to align the elements of the mesh according to the anisotropy of the solution. It turns out that this intuitive alignment is also necessary to prove sharp upper error bounds. In particular the proof employs specific interpolation error estimates. These interpolation estimates hold for isotropic meshes, but do not hold for general anisotropic meshes; instead the mesh has to have the aforementioned anisotropic alignment with the function to be interpolated. In order to quantify this alignment, we introduce a so-called alignment measure m1 (v, A, Th ) which was originally introduced in [13] for the identity matrix A and that we extend here to any matrix A. Definition 3.5 (Alignment measure) Let v ∈ H 1 (Ω), and T = {Th } be a family of triangulations of Ω. Define the alignment measure m1 : H 1 (Ω) × T 7→ R by (24)
m1 (v, A, Th ) :=
³X
> h−2 min,A,T kCA,T
∇vk2T
´1/2 .
kA1/2 ∇vk.
T ∈Th
One has m1 (v, A, Th ) & 1 since > kCA,T ∇vkT = kCT>A A1/2 ∇vkT & hmin,A,T kA1/2 ∇vkT .
For arbitrary isotropic meshes one obtains that m1 (v, Id, Th ) ∼ 1. The same is achieved for anisotropic meshes Th that are aligned with the anisotropic function v. Therefore the alignment measure is not an obstacle for reliable a posteriori error estimation. We refer to [13, 14] for discussions concerning this alignment measure. Now we recall the definition of the Cl´ement interpolation operator that maps a function from HΓ1D (Ω) into Vh,1 := {vh ∈ HΓ1D (Ω) : vh|T ∈ P1 , ∀T ∈ Th } ⊂ Vh . For that purpose, let the basis function ϕx ∈ Vh,1 associated with the node x be determined by the condition ϕx (y) = δx,y
∀y ∈ Nh ,
where Nh is the set of nodes of the triangulation included into Ω and ΓN . Then, the Cl´ement interpolation operator will be defined via these basis functions:
12
Definition 3.6 (Cl´ ement interpolation operator) The Cl´ement interpolation operator ICl : HΓ1D (Ω) → Vh,1 is defined by ¶ Xµ 1 Z ICl v := v ϕx , |ω x | ωx x∈N h
with ωx being the union of elements T of Th having x has vertex. Lemma 3.7 (Global interpolation error bounds) For each edge/face E, let us set (25) βE = max (hE,T h−1 min,A,T ). T ⊂ωE
Let v ∈ HΓ1D (Ω), then the following estimates hold: (26)
X
(27) (28)
X
αT−1
T ∈Th
X
|||ICl v||| . m1 (v, A, Th )|||v|||, αT−2 kv
− ICl vk2T . m1 (v, A, Th )2 |||v|||2 ,
T ∈Th
βE kv − ICl vk2E . m1 (v, A, Th )2 |||v|||2 .
E∈∂T \ΓD : βE =hE,T h−1 min,A,T
Note that in (28), every edge/face E ∈ Eh with E 6⊂ ΓD appears in the double sum at least once. Proof: Lemma 3.1 of [12] says that (29)
kICl vk . kvk, X > kCA,T kv − ICl vkT . 0 ∇vkT 0 ,
(30)
T 0 ⊂ωT > kCA,T ∇(v − ICl v)kT .
(31)
X
> kCA,T 0 ∇vkT 0 .
T 0 ⊂ωT
Multiplying the estimate (31) by h−1 min,A,T and summing the squares for all T yields X X > 2 > 2 h−2 h−2 min,A,T kCA,T ∇vkT . min,A,T kCA,T ∇(v − ICl v)kT . T ∈Th
T ∈Th
By the definition of the alignment measure, we get X 2 1/2 2 > (32) ∇vk2 . h−2 min,A,T kCA,T ∇(v − ICl v)kT . m1 (v, A, Th ) kA T ∈Th
13
In the same manner by (31) and (30) we have X > 2 2 1/2 h−2 (33) ∇vk2 , min,A,T kCA,T ∇ICl vkT . m1 (v, A, Th ) kA T ∈Th
X
(34)
2 2 1/2 h−2 ∇vk2 . min,A,T k(v − ICl v)kT . m1 (v, A, Th ) kA
T ∈Th
Now we remark that kA1/2 ∇ICl vk2 =
X
1/2
kAT ∇ICl vk2T
T ∈Th
≤
X
> kCT−1 k2 kCA,T ∇ICl vk2T A
T ∈Th
.
X
> 2 h−2 min,A,T kCA,T ∇ICl vkT .
T ∈Th
By the estimate (33), we conclude that kA1/2 ∇ICl vk2 . m1 (v, A, Th )2 kA1/2 ∇vk2 .
(35)
The estimates (29) and (35) prove (26). Let us go on with the estimate (27): X X X 2 αT−2 kv − ICl vk2T ≤ c0 kv − ICl vk2T + h−2 min,A,T kv − ICl vkT . T ∈Th
T ∈Th
T ∈Th
By (29) and (34), we obtain (27). For the last estimate, for any edge/face E of Eh , take any element T ∈ Th such that βE = hE,T h−1 min,A,T . Now we use a standard trace inequality on T to get > kv − ICl vk2E . h−1 E,T kv − ICl vkT (kv − ICl vkT + kCA,T ∇(v − ICl v)kT ).
Multiplying this estimate by βE (= hE,T h−1 min,A,T ), we get > βE kv − ICl vk2E . h−1 min,A,T kv − ICl vkT (kv − ICl vkT + kCA,T ∇(v − ICl v)kT ).
Multiplying this estimate by αT−1 and summing on E and then on T , we have obtained X X X αT−1 βE kv − ICl vk2E . αT−1 h−1 min,A,T kv − ICl vkT · T ∈Th
T ∈Th
E⊂∂T \ΓD : βE =hE,T h−1 min,A,T
> (kv − ICl vkT + kCA,T ∇(v − ICl v)kT ).
14
Now using the discrete Cauchy-Schwarz inequality, we arrive at !1/2 Ã X X X αT−1 βE kv − ICl vk2E . αT−2 kv − ICl vk2T · T ∈Th
E⊂∂T \ΓD : βE =hE,T h−1 min,A,T
T ∈Th
Ã
X
!1/2 2 > 2 h−2 min,A,T (kv − ICl vkT + kCA,T ∇(v − ICl v)kT )
.
T ∈Th
We conclude thanks to (27), (31) and (34). Remark 3.8 If A = εId, then the estimate (28) implies the estimate (21) of [16], since √ hE,T √ βE,T = ε & ε. hmin,T
4 4.1
Error estimator Definition of the error estimator
We investigate a residual error estimator. The exact element residual is defined by RT := f − Auh on T. Similarly the exact edge/face residual is [[A∇uh · nE ]]E g − A∇uh · n RE = 0
on E ∈ Ehint , on E ∈ Ehext ∩ ΓN , on E ∈ Ehext ∩ ΓD .
As usual, these exact residuals are replaced by some finite-dimensional approximation rT ∈ Pk0 (T ) and rE ∈ Pk1 (E) called approximate element residuals. Definition 4.1 (Residual error estimator) The local and global residual error estimators are defined by X X η 2 := ηT2 . ηT2 := αT2 krT k2T + αT βE−1 krE k2E , E∈∂T \ΓD : βE =hE,T h−1 min,A,T
The local and global approximation terms are defined by X X ζT2 := αT2 kRT 0 − rT 0 k2T 0 + αT βE−1 kRE − rE k2E , T 0 ⊂ωT
E∈∂T \ΓD
15
T ∈Th
ζ 2 :=
X T ∈Th
ζT2 .
4.2
Upper error bound
Theorem 4.2 Assume that δT satisfies (7). Let u be a solution of (3) and uh a solution of (6). Then the error is bounded as follows: (36)
|||u − uh ||| . m1 (u − uh , A, Th )(η + ζ).
Proof: By (4) we have (37) |||u − uh |||2 ≤ B(u − uh , u − uh ) ≤ B(u − uh , v − ICl v) + B(u − uh , ICl v), where for shortness we write v = u − uh . For the first term, elementwise integration by parts yields X X B(u − uh , v − ICl v) = (RT , v − ICl v)T + (RE , v − ICl v)E . T ∈Th
E∈Eh
By the continuous and by the discrete Cauchy-Schwarz inequality and the use of Lemma 3.7, we arrive at (38)
B(u − uh , v − ICl v) . m1 (v, A, Th )(η + ζ)|||v|||.
For the second term of the right-hand side of (37), we first estimate kA1/2 ∇ICl vkT . Indeed we first write kA1/2 ∇ICl vkT = kCT−> CT>A A1/2 ∇ICl vkT A > . h−1 min,A,T kCA,T ∇ICl vkT
. h−1 min,A,T kICl vkT , this last estimate coming from the inverse inequality on TˆTA and scaling arguments. This finally implies that −1/2
kA1/2 ∇ICl vkT . h−1 min,A,T c0
|||ICl v|||T .
On the other hand, we trivially have kA1/2 ∇ICl vkT . |||ICl v|||T = h−1 min,A,T hmin,A,T |||ICl v|||T , and, by the definition of αT , we have obtained (39)
kA1/2 ∇ICl vkT . h−1 min,A,T αT |||ICl v|||T . 16
Now using (3) and (6), we get B(u − uh , ICl v) = −
X
δT (RT , b · ∇ICl v)T ,
T ∈Th
and by the Cauchy-Schwarz inequality we obtain X B(u − uh , ICl v) ≤ δT kRT kT kA−1/2 bk∞,T kA1/2 ∇ICl vkT . T ∈Th
Using (39), we obtain B(u − uh , ICl v) .
X
δT kRT kT kA−1/2 bk∞,T h−1 min,A,T αT |||ICl v|||T ,
T ∈Th
and by the assumption on δT , we arrive at X B(u − uh , ICl v) . αT kRT kT |||ICl v|||T . T ∈Th
The discrete Cauchy-Schwarz inequality and the estimate (26) lead to B(u − uh , ICl v) . m1 (v, A, Th )(η + ζ)|||v|||. This estimate and (38) in the identity (37) lead to the conclusion.
4.3
Lower error bound
Theorem 4.3 For all elements T , the following local lower error bound holds: (40) where (41)
ηT . κωT |||u − uh |||ωT + ζT , © ª −1 −1/2 0 } + αT 0 kA 0 κω = max max{1, c kck bk . ∞,T ∞,T 0 0 T ⊂ω
Proof: As already mentioned elementwise integration by parts yields X X B(u − uh , w) = (RT , w)T + (RE , w)E , ∀w ∈ HΓ1D (Ω). (42) T ∈Th
E∈Eh
Element residual: For a fixed element T define wT = rT bT which belongs to From the definition of RT and using (42) with w = wT we have Z Z Z Z rT wT = (rT − RT )wT + RT wT = (rT − RT )wT + B(u − uh , wT ).
HΓ1D (Ω). T
T
T
T
17
Using the equivalence (17) and the estimate (5) we obtain krT k2T . krT − RT kT + |||u − uh |||T ( max{1, c−1 0 kck∞,T }|||wT |||T −1/2 + kA bk∞,T kwT kT ). By the inverse inequalities (17) and (19) we get krT kT . krT − RT kT + −1 −1/2 bk∞,T kwT kT ). |||u − uh |||T ( max{1, c−1 0 kck∞,T }αT +kA Multiplying this estimate by αT , we arrive at (43)
αT krT kT . κωT |||u − uh |||T + ζT .
Edge/face residual Fix an arbitrary edge/face E ∈ Eh \ ΓD . We apply (42) with w = wE , where wE := Fext (rE )bE,γE,T on TE,γE ⊂ T ⊂ ωE , where TE,γE,T is the squeezed element associated with T defined with the parameter γE,T ∈ (0, 1] that will be fixed later on. This yields (rE , wE )E = (rE − RE , wE )E + (RE , wE )E = (rE − RE , wE )E + B(u − uh , wE ) −
X
(RT , wE )T .
T ⊂ωE
Using the equivalence (18) and the estimate (5) we obtain X kRT kT kwE kT krE k2E . krE − RE kE kwE kE + +
X
T ⊂ωE −1/2 |||u − uh |||T ( max{1, c−1 bk∞,T kwE kT ). 0 kck∞,T }|||wE |||T + kA
T ⊂ωE
Using (20), (21) and (43) we have X 1/2 1/2 krE kE . krE − RE kE + γE,T hE,T αT−1 ζT +
X
³
T ⊂ωE 1/2
1/2
1/2
1/2
1/2
|||u − uh |||T kA−1/2 bk∞,T γE,T hE,T + c0 max{1, c−1 0 kck∞,T }γE,T hE,T
T ⊂ωE −1/2 1/2
1/2
1/2
−1 −1 + max{1, c−1 0 kck∞,T }γE,T hE,T hmin,A,T + γE,T hE,T αT κωT
18
´
1/2 −1/2
Multiplying by αT βE (44)
1/2 −1/2 αT βE krE kE
, we arrive at
. ζωE +
X
τE,T ζT + |||u− uh |||ωE
T ⊂ωE
3 X ³X
´ κT i +τE,T κωT ,
i=1
T ⊂ωE
where we have set κT 1 = kA−1/2 bk∞,T τE,T αT , 1/2
κT 2 = c0 max{1, c−1 0 kck∞,T }τE,T αT , −1/2 1/2
1/2 −1/2
−1 κT 3 = max{1, c−1 0 kck∞,T }γE,T hE,T hmin,A,T αT βE 1/2
1/2
.
−1/2 −1/2 βE
τE,T = γE,T hE,T αT
The conclusion follows from the estimates (43) and (44) if we can show that (45) (46)
κT i . κωE , i = 1, 2, 3, T ⊂ ωE , τE,T . 1,
κωE being defined in (41). In view of this definition we see that (45) holds for i = 1 if τE,T αT . αT , or equivalently if (46) holds. Since (46) is equivalent to γE,T . αT h−1 E,T βE and as we −1 have chosen βE = maxT ⊂ωE {hE,T hmin,A,T }, the estimate (46) holds if (47)
−1 −1 γE,T . αT h−1 E,T hE,T hmin,A,T = αT hmin,A,T
is satisfied. Similarly (45) holds for i = 2 if 1/2 1/2
1/2
1/2 −1/2
c0 γE,T hE,T αT βE
. 1,
i.e., −1 γE,T . αT−1 h−1 E,T c0 βE .
Again by the definition of βE , this last estimate holds if (48)
−1 γE,T . αT−1 h−1 min,A,T c0 .
Finally, (45) holds for i = 3 if 1/2 −1/2
−1/2 1/2
γE,T hE,T h−1 min,A,T αT βE 19
. 1.
Again by the definition of βE , this estimate holds if −1/2 1/2
1/2 −1/2
γE,T hE,T h−1 min,A,T αT βE,T . 1, or equivalently (49)
αT h−1 min,A,T . γE,T .
To satisfy these three conditions, (47)–(49), we take γE,T = αT h−1 min,A,T . First this right-hand side is ≤ 1 by the definition of αT . Secondly conditions (47) and (49) are automatically satisfied, and thirdly (48) reduces to αT2 c0 . 1, which also holds thanks to the definition of αT . Remark 4.4 In comparison with the case A = εId (see [10, 3, 16] for anisotropic meshes and [22, 25] for isotropic meshes), we could define a local mesh Peclet number by PeT = hmin,A,T kA−1/2 bk∞,T . Then the factor κT reduces to κT = max{1, c−1 0 kck∞,T 0 } + PeT . Furthermore in the case A = εId, our new Peclet number reduces to the standard one. Thus the lower bound is optimal if the local mesh Peclet number is sufficiently small. This is in congruence with Verf¨ urth’s paper [25].
5
Numerical results
The aim is to test the behavior of the estimated error in relationship with the true error. Therefore we use a test example with a known exact solution. We consider the problem −div (A∇u) = 0 in Ω = (0, 1)2 , u = g on Γ, with
µ A=
ε 0 0 1
¶ ,
ε = 10−k , k = 2, . . . , 10, and choose the boundary data g such that the exact solution contains a typical boundary layer of that problem, √
u(x, y) = e−πx/ 20
ε
sin πy,
1 0.2
0.8 0.3
0.4
0.6 0.5
0.6
0.7
0.4 0.8
0.9
0.2 1
0
Figure 3: Illustration of the solution ε = 10−2 N 64 256 1 024 4 096 16 384
|||e|||
η
0.25042 0.76160 0.13428 0.41209 0.068559 0.21368 0.034470 0.10851 0.017259 0.05463
ε = 10−10 η |||e||| 3.04 3.07 3.12 3.15 3.17
|||e|||
η
0.0072484 0.0047663 0.0028324 0.0015332 0.00078715
0.018535 0.011821 0.0074057 0.0042239 0.0022193
η |||e||| 2.56 2.48 2.61 2.76 2.82
Table 1: Test results see Figure 3. The√ mesh is piecewise uniform with an anisotropic part in the boundary strip Ω1 = (0, a0 ε| ln ε|) × (0, 1) where we made good experience with choosing a0 = 0.5. Both Ω1 and Ω2 = Ω \ Ω1 were subdivided into n × n congruent rectangles which were afterwards split into two triangles each (union-jack grid). In this way the aspect ratio of the elements is about ε−1/2 in Ω1 and about unity in Ω2 . In Table 1 we summarize the test results. We see that the meshes are appropriately chosen such that the error is proportional to the mesh size (when the mesh size halves then the number of elements multiplies by four). But most importantly for this paper, the effectivity index η/|||e||| is about 2.5 to 3.2 independent of N and ε, and does not blow up. This experiment illustrates the efficiency and reliability of our estimator.
21
6
Conclusions
We have proposed and rigorously analyzed a new a posteriori error estimate for the finite element approximation of anisotropic diffusion–convection–reaction equations with anisotropic finite elements. We have shown that this estimate is reliable and efficient. Numerical experiments confirm our theoretical predictions. Acknowledgement The authors thank Matthias G¨otz for performing the numerical tests.
References [1] M. Ainsworth and J. T. Oden. A Posteriori Error Estimation in Finite Element Analysis. Wiley, New York, 2000. [2] L. Angermann. Balanced a-posteriori error estimates for finite volume type discretizations of convection-dominated elliptic problems. Computing, 55(4):305– 323, 1995. [3] Th. Apel. Anisotropic finite elements: Local estimates and applications. Advances in Numerical Mathematics. Teubner, Stuttgart, 1999. [4] Th. Apel and G. Lube. Anisotropic mesh refinement in stabilized Galerkin methods. Numer. Math., 74:261–282, 1996. [5] R. Becker and R. Rannacher. An optimal control approach to a posteriori error estimation in finite element methods. Acta Numerica, 37:1–225, 2001. [6] P. G. Ciarlet. The finite element method for elliptic problems. North-Holland, Amsterdam, 1978. Reprinted by SIAM, Philadelphia, 2002. [7] F. Fierro and A. Veeser. A posteriori error estimators, gradient recovery by averaging, and superconvergence. Numer. Math., 103:267–298, 2006. [8] L. Formaggia, S. Perotto, and P. Zunino. An anisotropic a-posteriori error estimate for a convection-diffusion problem. Comput. Vis. Sci., 4(2):99–104, 2001. [9] E. H. Georgoulis, E. Hall, and P. Houston. Discontinuous Galerkin methods for advection–diffusion–reaction problems on anisotropically refined meshes. submitted, 2006. 22
[10] R. Hangleiter and G. Lube. Stabilized Galerkin methods and layer-adaptedgrids for elliptic problems. Comput. Methods Appl. Mech. Eng., 166(1–2):165–182, 1998. [11] D. Kay and D. Silvester. The reliability of local error estimators forconvection– diffusion equations. IMA J. Numer. Anal., 21(1):107–122, 2001. [12] G. Kunert. A posteriori error estimation for anisotropic tetrahedral and triangular finite element meshes. PhD thesis, TU Chemnitz, 1999. Logos, Berlin, 1999. [13] G. Kunert. An a posteriori residual error estimator for the finiteelement method on anisotropic tetrahedral meshes. Numer. Math., 86:471–490, 2000. [14] G. Kunert. A local problem error estimator for anisotropic tetrahedral finite element meshes. SIAM J. Numer. Anal., 39:668–689, 2001. [15] G. Kunert. Robust a posteriori error estimation for a singularlyperturbed reaction-diffusion equation on anisotropictetrahedral meshes. Adv. Comp. Math., 15:237–259, 2001. [16] G. Kunert. A posteriori error estimation for convection dominated problems on anisotropic meshes. Math. Methods Appl. Sci., 26:589–617, 2003. [17] J. Li. Quasioptimal uniformly convergent finite element methods for the elliptic boundary layer problem. Comput. Math. Appl., 33:11–22, 1997. [18] J. Li. Convergence and superconvergence analysis of finite element methods on highly nonuniform anisotropic meshes for singularly perturbed reactiondiffusion problems. Appl. Numer. Math., 36:129–154, 2001. [19] J. Li and M. F. Wheeler. Uniform convergence and superconvergence of mixed finite element methods on anisotropically refined grids. SIAM J. Numer. Anal., 38:770–798, 2000. [20] P. Neittaanma¨aki and S. Repin. Reliable methods for computer simulation: error control and a posteriori error estimates., volume 33 of Studies in Mathematics and its applications. Elsevier, Amsterdam, 2004. [21] M. Picasso. An anisotropic error indicator based on Zienkiewicz-Zhu error estimator:Application to elliptic and parabolic problems. SIAM J. Sci. Comput., 24(4):1328–1355, 2003. 23
[22] H.-G. Roos, M. Stynes, and L. Tobiska. Numerical methods for singularly perturbed differential equations. Convection-diffusion and flow problems. Springer, Berlin, 1996. [23] H.-G. Roos, M. Stynes, and L. Tobiska. Robust numerical methods for singularly perturbed differential equations. Convection-diffusion-reaction and flow problems. Springer, Berlin, 2008. [24] R. Verf¨ urth. A review of a posteriori error estimation andadaptive mesh– refinement techniques. Wiley and Teubner, Chichester and Stuttgart, 1996. [25] R. Verf¨ urth. A posteriori error estimators forconvection–diffusion equations. Numer. Math., 80(4):641–663, 1998. [26] R. Verf¨ urth. Robust a posteriori error estimators for the singularly perturbed reaction–diffusion equation. Numer. Math., 78:479–493, 1998.
24
Serge Nicaise†
May 22, 2009 Abstract: This paper presents an a posteriori residual error estimator for diffusion– convection–reaction problems with anisotropic diffusion, approximated by a SUPG finite element method on isotropic or anisotropic meshes in Rd , d = 2 or 3. The equivalence between the energy norm of the error and the residual error estimator is proved. Numerical tests confirm the theoretical results. Key words: anisotropic diffusion, SUPG, a posteriori error estimate. AMS subject classification: 65N30, 65N15
1
Introduction
This paper is devoted to the singularly perturbed diffusion–convection–reaction problem with special focus on anisotropic diffusion: for f ∈ L2 (Ω) and g ∈ L2 (ΓN ), let u be the solution of −div (A∇u) + b · ∇u + cu = f in Ω, u = 0 on ΓD , (1) A∇u · n = g on ΓN , where the matrix A and the functions b and c satisfy assumptions (A1) to (A6) below, and Ω ⊂ Rd , d = 2 or 3, is a bounded domain with a polygonal (d = 2) or polyhedral (d = 3) boundary Γ. This boundary is divided into two parts ΓD and ΓN , where Dirichlet and Neumann boundary conditions are imposed, respectively. ∗
Institut f¨ ur Mathematik und Bauinformatik, Universit¨at der Bundeswehr M¨ unchen, D–85577 Neubiberg, Germany, [email protected] † Universit´e de Valenciennes et du Hainaut Cambr´esis, LAMAV, ISTV, F–59313 - Valenciennes Cedex 9, France, [email protected]
1
We are particularly interested in the case when A tion, for instance the cases µ ¶ ε ε 0 0 A= (d = 2), or A = 0 1 0
becomes small in some direc 0 0 1 0 0 1
(d = 3),
ε > 0. In the case when ε is small with respect to b and c, the problem is singularly perturbed and the solution may generate sharp boundary or interior layers, where the solution of the limit problem (corresponding to ε = 0) is not smooth or does not satisfy the boundary condition. Let us quote [17, 18, 19] for the a priori error analysis in two dimensions. It is shown that anisotropic finite elements must be used in order to achieve convergence uniform in the perturbation parameter ε. There is a vast amount of literature on a posteriori error estimation. For singularly perturbed problems with convection we cite [2, 8, 11, 16, 21, 25, 26], where anisotropic finite element meshes were considered in [8, 16, 21] only. An anisotropic diffusion tensor is considered only in [7]. In this paper we combine all those ingredients and derive a residual type error estimator. We prove the reliability and efficiency of this error estimator where the dependence on ε is traced. The lower bound mainly depends on the local mesh Peclet number PeT := hmin,A,T kA−1/2 bk∞,T , therefore the efficiency is achieved if PeT ≤ c which is always satisfied in the absence of convection. The reliability is based on the introduction of an alignment measure as it was done by Kunert [12, 13]. The quantity is of the order one if the mesh is well adapted to the problem, see the discussion in Subsection 3.3. Let us mention that, to our knowledge, no approach is known that leads to twosides estimates on anisotropic meshes without any assumption on the mesh. The classical results as summarized in [1, 24] are obtained for isotropic meshes only. The dual weighted residual method, see [5] for an overview, is applied in [8, 9] on anisotropic meshes, but there is no estimate from below. The more recent approach in [20] is not yet analyzed for anisotropic meshes and two different error estimators are used for the upper and lower bounds. Let us finally mention the approach by Picasso [21] who considers anisotropic meshes and proves reliability for an estimator that depends on ∇(u − uh ) where ∇u is replaced in practice by a recovered gradient ∇R u. We note that we can control in the same way the alignment measure. In this paper we develop an estimator of residual type for problems with convection, reaction and anisotropic diffusion. For the discretization we use the h-version of the streamline upwind Petrov–Galerkin method (SUPG). Without the stabilization term, the method reduces to a standard Galerkin method and produces non-physical oscillations. We note that our error estimator works as well in this case. 2
In comparison with the paper [7], where a posteriori error estimation is investigated for an isotropic discretization of a problem with anisotropic diffusion but without convection, our residual error estimator allows to prove an optimal lower bound. The factor ε−1/2 in the upper bound in [7] is retained in our analysis, since the alignment measure is of the same order in the isotropic case. Our experiments show, however, that the effectivity index is bounded uniformly in ε on adequately refined anisotropic meshes. In this sense, our analysis is sharper. The outline of the paper is as follows. In Sections 2 and 3 we introduce the discretization, notation, and estimates for bubble functions and the interpolation. The a posteriori error estimator is introduced in Section 4 where also the the upper and lower bounds are proved. The paper is completed with a numerical test in Section 5 and with conclusions. As usual, we denote by L2 (.) the Lebesgue spaces and by H s (.), s ≥ 0, the standard Sobolev spaces. The usual norm and seminorm of H s (D) are denoted by k · ks,D and | · |s,D . For the sake of brevity the L2 (D)-norm will be denoted by k · kD and in the case D = Ω, we will drop the index Ω. The space HΓ1D (Ω) is defined, as usual, by HΓ1D (Ω) := {v ∈ H 1 (Ω) : v = 0 on Γ}. In the sequel the symbol | · | will denote either the Euclidean norm in Rd , d = 2 or 3, or the length of a line segment, or the measure of a domain of Rd . Finally the notation a . b means here and below that there exists a positive constant C independent of a and b (of the mesh size of the triangulation, as well as the diffusion matrix A) such that a ≤ C b. The notation a ∼ b means that a . b and b . a hold simultaneously.
2
Discretization of the diffusion–convection–reaction equation
Let Ω be a bounded open subset of Rd , d = 2 or 3, with a polygonal (d = 2) or polyhedral (d = 3) boundary Γ. This boundary is divided into two parts ΓD and ΓN , where Dirichlet and Neumann boundary conditions are imposed, respectively. We consider the standard elliptic problem: for f ∈ L2 (Ω) and g ∈ L2 (ΓN ), let u be the solution of (1) where A, b and c satisfy the following assumptions: (A1) (A2) (A3) (A4) (A5) (A6)
b ∈ W 1,∞ (Ω)d , c ∈ L∞ (Ω), ∃c0 ≥ 0 : c − 21 div b ≥ c0 and if c0 = 0 then c ≡ 0, b · n ≥ 0 on ΓN , A ∈ Rd×d is symmetric, ∃α0 > 0 : Aξ · ξ ≥ α0 , ∀ξ ∈ Rd . ΓD 6= ∅ if c0 = 0. 3
Note that the assumption ”if c0 = 0 then c ≡ 0” is not necessary for our proofs but simplifies the presentation. Now we define the weighted H 1 (semi-)norm Z 2 (2) |||u|||ω := (|A∇u|2 + c0 |u|2 ). ω
on a subdomain ω of Ω. Let us further introduce the space HΓ1D (Ω) = {v ∈ H 1 (Ω) : u = 0 on ΓD }. and the forms Z B(u, v) = ZΩ F (v) =
(A∇u · ∇v + b · ∇uv + cuv) dx, Z f v dx + gv dΓ(x).
Ω
ΓN
For further purposes, we denote by Bω the restriction of B on a subset ω of Ω, namely Z Bω (u, v) = (A∇u · ∇v + b · ∇uv + cuv) dx. ω
With this notation, the variational formulation of problem (1) reads: Find u ∈ solution of (3) B(u, v) = F (v) ∀v ∈ HΓ1D (Ω). HΓ1D (Ω)
The assumptions (A1) to (A6) guarantee that B is continuous and coercive, i.e., B satisfies (4) (5)
B(v, v) ≥ |||v|||2
∀v ∈ HΓ1D (Ω),
¡ ¢ −1/2 |Bω (u, v)| ≤ |||u|||ω max{1, c−1 bk∞,ω kvkω , 0 kck∞,ω }|||v|||ω + kA
for all u, v ∈ H 1 (ω). If c0 = 0, then the term c−1 0 kck∞,ω disappears and the maximum in (5) is equal to one. By the Lax-Milgram lemma, problem (3) has a unique solution u ∈ HΓ1D (Ω). Note that the case div b = 0, c = 0, i. e. c0 = 0, is admitted. It is excluded in several other publications. To approximate problem (3) by a finite element scheme we fix a family {Th }h>0 of meshes of Ω that satisfies the usual conformity conditions, cf. [6, Chapter 2]. In 2D we assume that all elements of Th are triangles, while in 3D the mesh is made up of tetrahedra. For T ∈ Th we denote by hT the diameter of T , and h = maxT ∈Th hT . 4
Let Vh be the subspace of HΓ1D (Ω) defined by Vh = {vh ∈ HΓ1D (Ω) : vh|T ∈ Pk (T ) ∀T ∈ Th }, where k is a positive integer. Problem (3) is now approximated by a Streamline Upwind Petrov Galerkin scheme (SUPG): Find uh ∈ Vh such that (6)
Bh (uh , vh ) = Fh (vh ) ∀vh ∈ Vh ,
where X
Bh (uh , vh ) = B(uh , vh ) +
δT (−div (A∇uh ) + b · ∇uh + cuh , b · ∇vh )T ,
T ∈Th
Fh (vh ) = F (vh ) +
X
δT (f, b · ∇vh )T .
T ∈Th
The parameters δT ≥ 0 should satisfy similar assumptions as in [16] where the case of isotropic diffusion was investigated: hmin,A,T , kA−1/2 bk∞,T ≤ 2(1 − α)µ−2 h2min,T kAk−1 2→2 , µ ¶−1 ≤ 2(1 − α)c0 max c(x)2
(7)
δT ≤
(8)
δT
(9)
δT
x∈T
if c 6≡ 0,
for all T ∈ Th and some α ∈ (0, 1). The element quantities hmin,T and hmin,A,T are introduced below, and µ is the constant in the inverse inequality k∇ · ∇vh kT ≤ µh−1 min,T k∇vh kT . Note further that (8) and (9) guarantee the coercivity Bh (vh , vh ) ≥ α|||vh |||2 with the above α ∈ (0, 1), compare [23, Lemma 3.25] for the case of isotropic meshes. The optimal choice of δT was discussed for the slightly different Galerkin– Least-Squares method and for the case of isotropic diffusion in [4]. This choice satisfies the conditions (7)–(9). We note also that the choice δT = 0 (pure Galerkin method) satisfies these conditions. Meanwhile it is well known that this choice is suited within boundary layers if adequately refined anisotropic meshes are used there [23, p. 391 ff.]. Outside the layers, the choice δT = 0 leads in general to non-physical oscillations. Therefore this choice is not advisable, but the error estimator still works and estimates the large error.
5
3
Analytical tools
Let us define Eh as the set of edges (d = 2) or faces (d = 3) of the triangulation and let Ehint = {E ∈ Eh /E ⊂ Ω} be the set of interior edges/faces of Th , while Ehext = Eh \ Ehint is the set of boundary edges/faces of Th . For an edge/face E of a 2D/3D element T we denote by nT,E the unit outward normal vector to T along E. Furthermore we fix one of the two normal vectors of E and denote it by nE . The jump of some function v across an edge/face E at a point y ∈ E is defined as ½ limα→+0 v(y + αnE ) − v(y − αnE ) ∀E ∈ Ehint , [[v(y)]]E := v(y) ∀E ∈ Ehext . Finally we will need local subdomains, also called patches. For any T ∈ Th , let, as usual, ωT be the union of all elements having a common vertex with T . Similarly let ωE be the union of the elements having E as edge/face.
3.1
Some anisotropic quantities
As explained in the introduction, anisotropic discretizations can be very advantageous or, in certain situations, even mandatory. More information and arguments concerning anisotropy can be found in [3, 12]. Let us shortly recall some useful anisotropic quantities from Kunert [12], see also [15, 16]. We start with an arbitrary (anisotropic) tetrahedron T . We enumerate its vertices so that P0 P1 is the longest edge, meas2 (4P0 P1 P2 ) ≥ meas2 (4P0 P1 P3 ), and meas1 (P1 P2 ) ≥ meas1 (P0 P2 ). Further, we introduce three orthogonal vectors pi,T of length hi,T := |pi,T |, as described in Figure 1. P3
p3
P2 p2
p1
P0
P1 Figure 1: Notation of tetrahedron T .
6
The minimal element size is particularly important; thus we define hmin,T := h3,T . The three main anisotropic directions pi,T play an important role in several proofs. They span the matrix CT := (p1,T , p2,T , p3,T ) ∈ R3×3 . This matrix may be considered as a transformation matrix which defines implicitly a reference element TˆT via TˆT := CT−1 (T − P~0 ), cf. Figure 2. Note in particular that the reference element TˆT is of size O(1). Pˆ3
Pˆ2
Pˆ0
Pˆ1
Figure 2: Reference tetrahedron TˆT . In 2D the notation is similar. For a triangle T the enumeration is as in the bottom triangle P0 P1 P2 of Figure 1. We set hmin,T := h2,T , and CT becomes a 2 × 2 matrix. The new idea is now to transform any T ∈ Th by the matrix A−1/2 . More precisely, we transform T into TA by the affine transformation (10)
FA,T : T → TA : x → A−1/2 (x − gT ) + gT ,
where gT is the center of gravity of T . This element TA is a triangle in 2D or a tetrahedron in 3D that can be isotropic or not. Therefore we use its anisotropic quantities hi,TA , hmin,TA , CTA as introduced before.
7
Remark 3.1 In 2D, take
µ
ε 0 0 1
A=
¶ ,
c = 1, b = 0 and Ω = (0, 1)2 , ΓN = ∅. If f is smooth then u has boundary layers near x1 = 0 and x1 = 1. The transformation ˜ : x → A−1/2 x, Ω→Ω replaces the problem (1) into ˜ = (0, √1 ) × (0, 1), −∆˜ u + u˜ = f˜ in Ω ε ˜ u˜ = 0 on ∂ Ω. ˜ Then the triangle Take for simplicity a uniform triangular triangulation of Ω. ˜ with nodes (0, 0), (h, 0), and (0, h) becomes the triangle T ∈ Ω with vertices T˜ ∈ Ω √ (0, 0), ( εh, 0), and (0, h) by the inverse transformation. This element T is a good one to capture adequately the boundary layer near x1 = 0. Moreover by using FA,T , the triangle T is transformed into an isotropic element TA which is a translation of T˜, and therefore CT˜ = h−1 Id. For further use, we denote (11)
hmin,TA0 . hmin,A,T = min 0 T ⊂ωT
Note that the composition of the transformation TˆTA → TA : xˆ → CTA xˆ + P0 , −1 with FA,T , see (10), yields the following transformation from TˆTA to T :
(12)
TˆTA → T : xˆ → CA,T xˆ + bT
with (13)
CA,T = A1/2 CTA .
Note that TˆTA depends on T and A but is of unit size in the sense of Figure 2. Finally we introduce a scaling factor αT that will be used quite often: (14)
−1/2
αT = min{c0
8
, hmin,A,T }.
−1/2
Here and below, we use the convention that c0 = +∞ if c0 = 0. For an edge/face E of an element T introduce the height hE,T = We finally require, as usual, that (15) (16)
|T | . |E|
|T | ∼ |T 0 | if T ∩ T 0 6= ∅, the number of elements containing a vertex x is bounded uniformly.
Remark 3.2 If A = εId (the case treated by Kunert in [15, 16]), then TA,T is simply a homothetic transformation of T with a factor ε−1/2 . Therefore the matrix CεId,T defined by (13) is equal to CT . Moreover we have hmin,εId,T = ε−1/2 hmin,T . This last property implies that −1/2
αT = min{c0
, ε−1/2 hmin,T },
which is exactly the scaling factor introduced in [16].
3.2
Bubble functions, extension operator, and inverse inequalities
For our further analysis we require standard bubble functions and extension operators that satisfy certain properties recalled here for the sake of completeness. We need two types of bubble functions, namely bT and bE associated with an element T and an edge E, respectively. For a triangle or a tetrahedron T , denoting by λaTi , i = 1, · · · , d + 1, the barycentric coordinates of T and by aE i , i = 1, · · · , d, the vertices of the edge/face E ⊂ ∂T we recall that bT =
d+1 Y
λaTi and bE =
d Y
λaEi .
i=1
i=1
We note that bT = 0 on ∂T,
bE = 0 on ∂ωE ,
kbT k∞,T = kbE k∞,ωE ∼ 1.
In 2D, denote by T¯ the standard reference element with vertices (0, 0), (1, 0), and (0, 1). For an edge E¯ of T¯ included into the x¯1 axis, the extension Fext (vE¯ ) of ¯ to T¯ is defined by Fext (vE¯ )(¯ vE¯ ∈ C(E) x1 , x¯2 ) = vE¯ (¯ x1 ). The extension operator Fext (vE ) of vE ∈ C(E) to T for an edge E ⊂ ∂T is obtained using the affine transformation mapping T to T¯ and E to E¯ and the extension operator defined above. We proceed similarly in 3D. Now we may state the so-called inverse inequalities that are proved using classical scaling techniques, cf. [24] for the isotropic case and [12] for the anisotropic case. 9
Lemma 3.3 (Inverse inequalities) Let vT ∈ Pk0 (T ) and vE ∈ Pk1 (E), for some nonnegative integers k0 and k1 . Then the following inequalities hold, with the constants in the inequality depending on the polynomial degrees k0 or k1 but not on T , E or vT , vE . 1/2
(17)
kvT bT kT ∼ kvT kT ,
(18) (19)
kvE bE kE ∼ kvE kE , |||vT bT |||T . αT−1 kvT kT .
1/2
Proof: The equivalence (17) and (18) are proved in [12], see also Lemma 1 of [16]. For the last estimate, we write |||vT bT |||2T = c0 kvT bT k2T + kA1/2 ∇(vT bT )k2T ≤ c0 kvT k2T + kCT−T k2 kCTTA A1/2 ∇(vT bT )k2T A Now recalling that kCT−T k ∼ h−1 min,TA and using the affine transformation (12), we A obtain ˆ vT ˆbT )k2ˆ . |||vT bT |||2T . c0 kvT k2T + h−2 min,A,T |T |k∇(ˆ T Since Tˆ is regular in Ciarlet’s sense we can use the inverse inequality with hTˆ ∼ 1 to deduce that |||vT bT |||2T . c0 kvT k2T + h−2 vT k2Tˆ . min,A,T |T |kˆ Going back to T again using the affine transformation (12), we obtain (19). As usual for singularly perturbed problems, we need to use squeezed edge/face bubble functions bE,γ . Here according to our previous point of view, they are defined through the transformation FA,T from (10). Namely for a fixed edge/face E of T , the mapping (10) transforms E into an edge/face EA of TA . Now for a parameter γ ∈ (0, 1], we define the squeezed element TEA ,γ of TA as in [16]. The squeezed element TE,γ of T is simply the element obtained by the inverse transformation −1 TE,γ = FA,T TEA ,γ .
Note that TE,γ is the usual squeezed element on T with the parameter γ depending on A. For the sake of simplicity we do not write this dependence. The squeezed edge/face bubble function bE,γ is defined on the two elements T1,E,γ and T2,E,γ sharing γ, as the usual edge/face bubble function on these elements and extended by zero outside T1,E,γ ∪ T2,E,γ 10
Lemma 3.4 (Further Inverse inequalities) Under the assumptions of Lemma 3.3, we have 1/2
(20)
kbE,γ Fext (vE )kT . γ 1/2 hE,T kvE kE ,
(21)
kA1/2 ∇(bE,γ Fext (vE ))kT . γ −1/2 hE,T h−1 min,A,T kvE kE .
1/2
Proof: Scaling arguments yield ˜ext (˜ kbE,γ Fext (vE )kT = |T |1/2 |TA |−1/2 k˜bE,γ F vE )kTA , where we write v˜(˜ x) = v(x). Now using Lemma 2 of [16] in TA , we have ˜ext (˜ |TA |−1/2 k˜bE,γ F vE )kTA . γ 1/2 |EA |−1/2 k˜ vE kEA . Again scaling arguments lead to |EA |−1/2 k˜ vE kEA . |E|−1/2 kvE kE .
(22)
The three above estimates imply (20). For the second estimate, scaling arguments yield ˜ ˜bE,γ F ˜ext (˜ kA1/2 ∇(bE,γ Fext (vE ))kT = |T |1/2 |TA |−1/2 k∇( vE ))kTA . Again Lemma 2 of [16] applied in TA leads to kA1/2 ∇(bE,γ Fext (vE ))kT . |T |1/2 γ 1/2 |EA |−1/2 k˜ vE kEA min{γhEA ,TA , hmin,TA }−1 . Using the estimate (22), we arrive at 1/2
kA1/2 ∇(bE,γ Fext (vE ))kT . γ 1/2 hE,T kvE kE min{γhEA ,TA , hmin,TA }−1 . The estimate (21) will be proved if we can show that min{γhEA ,TA , hmin,TA }−1 . γ −1 h−1 min,A,T , or equivalently (23)
min{γhEA ,TA , hmin,TA } & γhmin,A,T .
But it was proved in Lemma 3.1 of [15] that hEA ,TA & hmin,TA . Since γ ∈ (0, 1] we then have γhEA ,TA & γhmin,TA ≥ γhmin,A,T and hmin,TA ≥ γhmin,A,T . This leads to (23). 11
3.3
Anisotropic interpolation error estimates
In order to obtain an accurate discrete solution uh , it is obviously helpful to align the elements of the mesh according to the anisotropy of the solution. It turns out that this intuitive alignment is also necessary to prove sharp upper error bounds. In particular the proof employs specific interpolation error estimates. These interpolation estimates hold for isotropic meshes, but do not hold for general anisotropic meshes; instead the mesh has to have the aforementioned anisotropic alignment with the function to be interpolated. In order to quantify this alignment, we introduce a so-called alignment measure m1 (v, A, Th ) which was originally introduced in [13] for the identity matrix A and that we extend here to any matrix A. Definition 3.5 (Alignment measure) Let v ∈ H 1 (Ω), and T = {Th } be a family of triangulations of Ω. Define the alignment measure m1 : H 1 (Ω) × T 7→ R by (24)
m1 (v, A, Th ) :=
³X
> h−2 min,A,T kCA,T
∇vk2T
´1/2 .
kA1/2 ∇vk.
T ∈Th
One has m1 (v, A, Th ) & 1 since > kCA,T ∇vkT = kCT>A A1/2 ∇vkT & hmin,A,T kA1/2 ∇vkT .
For arbitrary isotropic meshes one obtains that m1 (v, Id, Th ) ∼ 1. The same is achieved for anisotropic meshes Th that are aligned with the anisotropic function v. Therefore the alignment measure is not an obstacle for reliable a posteriori error estimation. We refer to [13, 14] for discussions concerning this alignment measure. Now we recall the definition of the Cl´ement interpolation operator that maps a function from HΓ1D (Ω) into Vh,1 := {vh ∈ HΓ1D (Ω) : vh|T ∈ P1 , ∀T ∈ Th } ⊂ Vh . For that purpose, let the basis function ϕx ∈ Vh,1 associated with the node x be determined by the condition ϕx (y) = δx,y
∀y ∈ Nh ,
where Nh is the set of nodes of the triangulation included into Ω and ΓN . Then, the Cl´ement interpolation operator will be defined via these basis functions:
12
Definition 3.6 (Cl´ ement interpolation operator) The Cl´ement interpolation operator ICl : HΓ1D (Ω) → Vh,1 is defined by ¶ Xµ 1 Z ICl v := v ϕx , |ω x | ωx x∈N h
with ωx being the union of elements T of Th having x has vertex. Lemma 3.7 (Global interpolation error bounds) For each edge/face E, let us set (25) βE = max (hE,T h−1 min,A,T ). T ⊂ωE
Let v ∈ HΓ1D (Ω), then the following estimates hold: (26)
X
(27) (28)
X
αT−1
T ∈Th
X
|||ICl v||| . m1 (v, A, Th )|||v|||, αT−2 kv
− ICl vk2T . m1 (v, A, Th )2 |||v|||2 ,
T ∈Th
βE kv − ICl vk2E . m1 (v, A, Th )2 |||v|||2 .
E∈∂T \ΓD : βE =hE,T h−1 min,A,T
Note that in (28), every edge/face E ∈ Eh with E 6⊂ ΓD appears in the double sum at least once. Proof: Lemma 3.1 of [12] says that (29)
kICl vk . kvk, X > kCA,T kv − ICl vkT . 0 ∇vkT 0 ,
(30)
T 0 ⊂ωT > kCA,T ∇(v − ICl v)kT .
(31)
X
> kCA,T 0 ∇vkT 0 .
T 0 ⊂ωT
Multiplying the estimate (31) by h−1 min,A,T and summing the squares for all T yields X X > 2 > 2 h−2 h−2 min,A,T kCA,T ∇vkT . min,A,T kCA,T ∇(v − ICl v)kT . T ∈Th
T ∈Th
By the definition of the alignment measure, we get X 2 1/2 2 > (32) ∇vk2 . h−2 min,A,T kCA,T ∇(v − ICl v)kT . m1 (v, A, Th ) kA T ∈Th
13
In the same manner by (31) and (30) we have X > 2 2 1/2 h−2 (33) ∇vk2 , min,A,T kCA,T ∇ICl vkT . m1 (v, A, Th ) kA T ∈Th
X
(34)
2 2 1/2 h−2 ∇vk2 . min,A,T k(v − ICl v)kT . m1 (v, A, Th ) kA
T ∈Th
Now we remark that kA1/2 ∇ICl vk2 =
X
1/2
kAT ∇ICl vk2T
T ∈Th
≤
X
> kCT−1 k2 kCA,T ∇ICl vk2T A
T ∈Th
.
X
> 2 h−2 min,A,T kCA,T ∇ICl vkT .
T ∈Th
By the estimate (33), we conclude that kA1/2 ∇ICl vk2 . m1 (v, A, Th )2 kA1/2 ∇vk2 .
(35)
The estimates (29) and (35) prove (26). Let us go on with the estimate (27): X X X 2 αT−2 kv − ICl vk2T ≤ c0 kv − ICl vk2T + h−2 min,A,T kv − ICl vkT . T ∈Th
T ∈Th
T ∈Th
By (29) and (34), we obtain (27). For the last estimate, for any edge/face E of Eh , take any element T ∈ Th such that βE = hE,T h−1 min,A,T . Now we use a standard trace inequality on T to get > kv − ICl vk2E . h−1 E,T kv − ICl vkT (kv − ICl vkT + kCA,T ∇(v − ICl v)kT ).
Multiplying this estimate by βE (= hE,T h−1 min,A,T ), we get > βE kv − ICl vk2E . h−1 min,A,T kv − ICl vkT (kv − ICl vkT + kCA,T ∇(v − ICl v)kT ).
Multiplying this estimate by αT−1 and summing on E and then on T , we have obtained X X X αT−1 βE kv − ICl vk2E . αT−1 h−1 min,A,T kv − ICl vkT · T ∈Th
T ∈Th
E⊂∂T \ΓD : βE =hE,T h−1 min,A,T
> (kv − ICl vkT + kCA,T ∇(v − ICl v)kT ).
14
Now using the discrete Cauchy-Schwarz inequality, we arrive at !1/2 Ã X X X αT−1 βE kv − ICl vk2E . αT−2 kv − ICl vk2T · T ∈Th
E⊂∂T \ΓD : βE =hE,T h−1 min,A,T
T ∈Th
Ã
X
!1/2 2 > 2 h−2 min,A,T (kv − ICl vkT + kCA,T ∇(v − ICl v)kT )
.
T ∈Th
We conclude thanks to (27), (31) and (34). Remark 3.8 If A = εId, then the estimate (28) implies the estimate (21) of [16], since √ hE,T √ βE,T = ε & ε. hmin,T
4 4.1
Error estimator Definition of the error estimator
We investigate a residual error estimator. The exact element residual is defined by RT := f − Auh on T. Similarly the exact edge/face residual is [[A∇uh · nE ]]E g − A∇uh · n RE = 0
on E ∈ Ehint , on E ∈ Ehext ∩ ΓN , on E ∈ Ehext ∩ ΓD .
As usual, these exact residuals are replaced by some finite-dimensional approximation rT ∈ Pk0 (T ) and rE ∈ Pk1 (E) called approximate element residuals. Definition 4.1 (Residual error estimator) The local and global residual error estimators are defined by X X η 2 := ηT2 . ηT2 := αT2 krT k2T + αT βE−1 krE k2E , E∈∂T \ΓD : βE =hE,T h−1 min,A,T
The local and global approximation terms are defined by X X ζT2 := αT2 kRT 0 − rT 0 k2T 0 + αT βE−1 kRE − rE k2E , T 0 ⊂ωT
E∈∂T \ΓD
15
T ∈Th
ζ 2 :=
X T ∈Th
ζT2 .
4.2
Upper error bound
Theorem 4.2 Assume that δT satisfies (7). Let u be a solution of (3) and uh a solution of (6). Then the error is bounded as follows: (36)
|||u − uh ||| . m1 (u − uh , A, Th )(η + ζ).
Proof: By (4) we have (37) |||u − uh |||2 ≤ B(u − uh , u − uh ) ≤ B(u − uh , v − ICl v) + B(u − uh , ICl v), where for shortness we write v = u − uh . For the first term, elementwise integration by parts yields X X B(u − uh , v − ICl v) = (RT , v − ICl v)T + (RE , v − ICl v)E . T ∈Th
E∈Eh
By the continuous and by the discrete Cauchy-Schwarz inequality and the use of Lemma 3.7, we arrive at (38)
B(u − uh , v − ICl v) . m1 (v, A, Th )(η + ζ)|||v|||.
For the second term of the right-hand side of (37), we first estimate kA1/2 ∇ICl vkT . Indeed we first write kA1/2 ∇ICl vkT = kCT−> CT>A A1/2 ∇ICl vkT A > . h−1 min,A,T kCA,T ∇ICl vkT
. h−1 min,A,T kICl vkT , this last estimate coming from the inverse inequality on TˆTA and scaling arguments. This finally implies that −1/2
kA1/2 ∇ICl vkT . h−1 min,A,T c0
|||ICl v|||T .
On the other hand, we trivially have kA1/2 ∇ICl vkT . |||ICl v|||T = h−1 min,A,T hmin,A,T |||ICl v|||T , and, by the definition of αT , we have obtained (39)
kA1/2 ∇ICl vkT . h−1 min,A,T αT |||ICl v|||T . 16
Now using (3) and (6), we get B(u − uh , ICl v) = −
X
δT (RT , b · ∇ICl v)T ,
T ∈Th
and by the Cauchy-Schwarz inequality we obtain X B(u − uh , ICl v) ≤ δT kRT kT kA−1/2 bk∞,T kA1/2 ∇ICl vkT . T ∈Th
Using (39), we obtain B(u − uh , ICl v) .
X
δT kRT kT kA−1/2 bk∞,T h−1 min,A,T αT |||ICl v|||T ,
T ∈Th
and by the assumption on δT , we arrive at X B(u − uh , ICl v) . αT kRT kT |||ICl v|||T . T ∈Th
The discrete Cauchy-Schwarz inequality and the estimate (26) lead to B(u − uh , ICl v) . m1 (v, A, Th )(η + ζ)|||v|||. This estimate and (38) in the identity (37) lead to the conclusion.
4.3
Lower error bound
Theorem 4.3 For all elements T , the following local lower error bound holds: (40) where (41)
ηT . κωT |||u − uh |||ωT + ζT , © ª −1 −1/2 0 } + αT 0 kA 0 κω = max max{1, c kck bk . ∞,T ∞,T 0 0 T ⊂ω
Proof: As already mentioned elementwise integration by parts yields X X B(u − uh , w) = (RT , w)T + (RE , w)E , ∀w ∈ HΓ1D (Ω). (42) T ∈Th
E∈Eh
Element residual: For a fixed element T define wT = rT bT which belongs to From the definition of RT and using (42) with w = wT we have Z Z Z Z rT wT = (rT − RT )wT + RT wT = (rT − RT )wT + B(u − uh , wT ).
HΓ1D (Ω). T
T
T
T
17
Using the equivalence (17) and the estimate (5) we obtain krT k2T . krT − RT kT + |||u − uh |||T ( max{1, c−1 0 kck∞,T }|||wT |||T −1/2 + kA bk∞,T kwT kT ). By the inverse inequalities (17) and (19) we get krT kT . krT − RT kT + −1 −1/2 bk∞,T kwT kT ). |||u − uh |||T ( max{1, c−1 0 kck∞,T }αT +kA Multiplying this estimate by αT , we arrive at (43)
αT krT kT . κωT |||u − uh |||T + ζT .
Edge/face residual Fix an arbitrary edge/face E ∈ Eh \ ΓD . We apply (42) with w = wE , where wE := Fext (rE )bE,γE,T on TE,γE ⊂ T ⊂ ωE , where TE,γE,T is the squeezed element associated with T defined with the parameter γE,T ∈ (0, 1] that will be fixed later on. This yields (rE , wE )E = (rE − RE , wE )E + (RE , wE )E = (rE − RE , wE )E + B(u − uh , wE ) −
X
(RT , wE )T .
T ⊂ωE
Using the equivalence (18) and the estimate (5) we obtain X kRT kT kwE kT krE k2E . krE − RE kE kwE kE + +
X
T ⊂ωE −1/2 |||u − uh |||T ( max{1, c−1 bk∞,T kwE kT ). 0 kck∞,T }|||wE |||T + kA
T ⊂ωE
Using (20), (21) and (43) we have X 1/2 1/2 krE kE . krE − RE kE + γE,T hE,T αT−1 ζT +
X
³
T ⊂ωE 1/2
1/2
1/2
1/2
1/2
|||u − uh |||T kA−1/2 bk∞,T γE,T hE,T + c0 max{1, c−1 0 kck∞,T }γE,T hE,T
T ⊂ωE −1/2 1/2
1/2
1/2
−1 −1 + max{1, c−1 0 kck∞,T }γE,T hE,T hmin,A,T + γE,T hE,T αT κωT
18
´
1/2 −1/2
Multiplying by αT βE (44)
1/2 −1/2 αT βE krE kE
, we arrive at
. ζωE +
X
τE,T ζT + |||u− uh |||ωE
T ⊂ωE
3 X ³X
´ κT i +τE,T κωT ,
i=1
T ⊂ωE
where we have set κT 1 = kA−1/2 bk∞,T τE,T αT , 1/2
κT 2 = c0 max{1, c−1 0 kck∞,T }τE,T αT , −1/2 1/2
1/2 −1/2
−1 κT 3 = max{1, c−1 0 kck∞,T }γE,T hE,T hmin,A,T αT βE 1/2
1/2
.
−1/2 −1/2 βE
τE,T = γE,T hE,T αT
The conclusion follows from the estimates (43) and (44) if we can show that (45) (46)
κT i . κωE , i = 1, 2, 3, T ⊂ ωE , τE,T . 1,
κωE being defined in (41). In view of this definition we see that (45) holds for i = 1 if τE,T αT . αT , or equivalently if (46) holds. Since (46) is equivalent to γE,T . αT h−1 E,T βE and as we −1 have chosen βE = maxT ⊂ωE {hE,T hmin,A,T }, the estimate (46) holds if (47)
−1 −1 γE,T . αT h−1 E,T hE,T hmin,A,T = αT hmin,A,T
is satisfied. Similarly (45) holds for i = 2 if 1/2 1/2
1/2
1/2 −1/2
c0 γE,T hE,T αT βE
. 1,
i.e., −1 γE,T . αT−1 h−1 E,T c0 βE .
Again by the definition of βE , this last estimate holds if (48)
−1 γE,T . αT−1 h−1 min,A,T c0 .
Finally, (45) holds for i = 3 if 1/2 −1/2
−1/2 1/2
γE,T hE,T h−1 min,A,T αT βE 19
. 1.
Again by the definition of βE , this estimate holds if −1/2 1/2
1/2 −1/2
γE,T hE,T h−1 min,A,T αT βE,T . 1, or equivalently (49)
αT h−1 min,A,T . γE,T .
To satisfy these three conditions, (47)–(49), we take γE,T = αT h−1 min,A,T . First this right-hand side is ≤ 1 by the definition of αT . Secondly conditions (47) and (49) are automatically satisfied, and thirdly (48) reduces to αT2 c0 . 1, which also holds thanks to the definition of αT . Remark 4.4 In comparison with the case A = εId (see [10, 3, 16] for anisotropic meshes and [22, 25] for isotropic meshes), we could define a local mesh Peclet number by PeT = hmin,A,T kA−1/2 bk∞,T . Then the factor κT reduces to κT = max{1, c−1 0 kck∞,T 0 } + PeT . Furthermore in the case A = εId, our new Peclet number reduces to the standard one. Thus the lower bound is optimal if the local mesh Peclet number is sufficiently small. This is in congruence with Verf¨ urth’s paper [25].
5
Numerical results
The aim is to test the behavior of the estimated error in relationship with the true error. Therefore we use a test example with a known exact solution. We consider the problem −div (A∇u) = 0 in Ω = (0, 1)2 , u = g on Γ, with
µ A=
ε 0 0 1
¶ ,
ε = 10−k , k = 2, . . . , 10, and choose the boundary data g such that the exact solution contains a typical boundary layer of that problem, √
u(x, y) = e−πx/ 20
ε
sin πy,
1 0.2
0.8 0.3
0.4
0.6 0.5
0.6
0.7
0.4 0.8
0.9
0.2 1
0
Figure 3: Illustration of the solution ε = 10−2 N 64 256 1 024 4 096 16 384
|||e|||
η
0.25042 0.76160 0.13428 0.41209 0.068559 0.21368 0.034470 0.10851 0.017259 0.05463
ε = 10−10 η |||e||| 3.04 3.07 3.12 3.15 3.17
|||e|||
η
0.0072484 0.0047663 0.0028324 0.0015332 0.00078715
0.018535 0.011821 0.0074057 0.0042239 0.0022193
η |||e||| 2.56 2.48 2.61 2.76 2.82
Table 1: Test results see Figure 3. The√ mesh is piecewise uniform with an anisotropic part in the boundary strip Ω1 = (0, a0 ε| ln ε|) × (0, 1) where we made good experience with choosing a0 = 0.5. Both Ω1 and Ω2 = Ω \ Ω1 were subdivided into n × n congruent rectangles which were afterwards split into two triangles each (union-jack grid). In this way the aspect ratio of the elements is about ε−1/2 in Ω1 and about unity in Ω2 . In Table 1 we summarize the test results. We see that the meshes are appropriately chosen such that the error is proportional to the mesh size (when the mesh size halves then the number of elements multiplies by four). But most importantly for this paper, the effectivity index η/|||e||| is about 2.5 to 3.2 independent of N and ε, and does not blow up. This experiment illustrates the efficiency and reliability of our estimator.
21
6
Conclusions
We have proposed and rigorously analyzed a new a posteriori error estimate for the finite element approximation of anisotropic diffusion–convection–reaction equations with anisotropic finite elements. We have shown that this estimate is reliable and efficient. Numerical experiments confirm our theoretical predictions. Acknowledgement The authors thank Matthias G¨otz for performing the numerical tests.
References [1] M. Ainsworth and J. T. Oden. A Posteriori Error Estimation in Finite Element Analysis. Wiley, New York, 2000. [2] L. Angermann. Balanced a-posteriori error estimates for finite volume type discretizations of convection-dominated elliptic problems. Computing, 55(4):305– 323, 1995. [3] Th. Apel. Anisotropic finite elements: Local estimates and applications. Advances in Numerical Mathematics. Teubner, Stuttgart, 1999. [4] Th. Apel and G. Lube. Anisotropic mesh refinement in stabilized Galerkin methods. Numer. Math., 74:261–282, 1996. [5] R. Becker and R. Rannacher. An optimal control approach to a posteriori error estimation in finite element methods. Acta Numerica, 37:1–225, 2001. [6] P. G. Ciarlet. The finite element method for elliptic problems. North-Holland, Amsterdam, 1978. Reprinted by SIAM, Philadelphia, 2002. [7] F. Fierro and A. Veeser. A posteriori error estimators, gradient recovery by averaging, and superconvergence. Numer. Math., 103:267–298, 2006. [8] L. Formaggia, S. Perotto, and P. Zunino. An anisotropic a-posteriori error estimate for a convection-diffusion problem. Comput. Vis. Sci., 4(2):99–104, 2001. [9] E. H. Georgoulis, E. Hall, and P. Houston. Discontinuous Galerkin methods for advection–diffusion–reaction problems on anisotropically refined meshes. submitted, 2006. 22
[10] R. Hangleiter and G. Lube. Stabilized Galerkin methods and layer-adaptedgrids for elliptic problems. Comput. Methods Appl. Mech. Eng., 166(1–2):165–182, 1998. [11] D. Kay and D. Silvester. The reliability of local error estimators forconvection– diffusion equations. IMA J. Numer. Anal., 21(1):107–122, 2001. [12] G. Kunert. A posteriori error estimation for anisotropic tetrahedral and triangular finite element meshes. PhD thesis, TU Chemnitz, 1999. Logos, Berlin, 1999. [13] G. Kunert. An a posteriori residual error estimator for the finiteelement method on anisotropic tetrahedral meshes. Numer. Math., 86:471–490, 2000. [14] G. Kunert. A local problem error estimator for anisotropic tetrahedral finite element meshes. SIAM J. Numer. Anal., 39:668–689, 2001. [15] G. Kunert. Robust a posteriori error estimation for a singularlyperturbed reaction-diffusion equation on anisotropictetrahedral meshes. Adv. Comp. Math., 15:237–259, 2001. [16] G. Kunert. A posteriori error estimation for convection dominated problems on anisotropic meshes. Math. Methods Appl. Sci., 26:589–617, 2003. [17] J. Li. Quasioptimal uniformly convergent finite element methods for the elliptic boundary layer problem. Comput. Math. Appl., 33:11–22, 1997. [18] J. Li. Convergence and superconvergence analysis of finite element methods on highly nonuniform anisotropic meshes for singularly perturbed reactiondiffusion problems. Appl. Numer. Math., 36:129–154, 2001. [19] J. Li and M. F. Wheeler. Uniform convergence and superconvergence of mixed finite element methods on anisotropically refined grids. SIAM J. Numer. Anal., 38:770–798, 2000. [20] P. Neittaanma¨aki and S. Repin. Reliable methods for computer simulation: error control and a posteriori error estimates., volume 33 of Studies in Mathematics and its applications. Elsevier, Amsterdam, 2004. [21] M. Picasso. An anisotropic error indicator based on Zienkiewicz-Zhu error estimator:Application to elliptic and parabolic problems. SIAM J. Sci. Comput., 24(4):1328–1355, 2003. 23
[22] H.-G. Roos, M. Stynes, and L. Tobiska. Numerical methods for singularly perturbed differential equations. Convection-diffusion and flow problems. Springer, Berlin, 1996. [23] H.-G. Roos, M. Stynes, and L. Tobiska. Robust numerical methods for singularly perturbed differential equations. Convection-diffusion-reaction and flow problems. Springer, Berlin, 2008. [24] R. Verf¨ urth. A review of a posteriori error estimation andadaptive mesh– refinement techniques. Wiley and Teubner, Chichester and Stuttgart, 1996. [25] R. Verf¨ urth. A posteriori error estimators forconvection–diffusion equations. Numer. Math., 80(4):641–663, 1998. [26] R. Verf¨ urth. Robust a posteriori error estimators for the singularly perturbed reaction–diffusion equation. Numer. Math., 78:479–493, 1998.
24
