SIAM J. NUMER. ANAL. Vol. 44, No. 3, pp. 949–978
c 2006 Society for Industrial and Applied Mathematics
A POSTERIORI ERROR ESTIMATIONS OF SOME CELL CENTERED FINITE VOLUME METHODS FOR DIFFUSION-CONVECTION-REACTION PROBLEMS∗ SERGE NICAISE† Abstract. This paper presents an a posteriori residual error estimator for diffusion-convectionreaction problems approximated by some cell centered finite volume methods on isotropic or anisotropic meshes in Rd , d = 2 or 3. For that purpose we built a reconstructed approximation, which is an appropriate interpolant of the finite volume solution. The error is then the difference between the exact solution and this interpolant. The residual error estimator is based on the jump of the normal derivative of the interpolant. We then prove the equivalence between the energy norm of the error and the residual error estimator. Some numerical tests confirm our theoretical results. Key words. finite volume method, cell centered method, a posteriori error estimates AMS subject classifications. 65N30, 65N15 DOI. 10.1137/040611483
1. Introduction. The finite volume method is a well-adapted method for the discretization of various partial differential equations and is largely used by engineers [31]. The mathematical analysis of the method has started only recently. Existence and uniqueness results as well as a priori error estimates are now available for quite a large class of problems; see [10] and the references cited there. For finite element methods, a posteriori error estimates are now well understood for a large class of equations; see, for instance, [34]. On the other hand, for finite volume methods, such estimates are not well developed and up to now only a few such results had been obtained. Let us quote [14, 29, 1, 12, 13] for cell centered finite volume methods, [23, 25, 33, 30] for vertex centered methods, and [4, 5, 21, 22, 20] for finite volume element methods. Recently we obtained a posteriori error estimates of residual type for some cell centered finite volume methods for the Laplace equation in a bounded domain of Rd , d = 2 or 3 [27, 26]. This estimator is based on the use of a reconstructed approximation, namely, an interpolant of Morley type of the finite volume solution. The first goal of the present paper is to extend the previous analysis to diffusion-convectionreaction problems that eventually develop boundary or interior layers. As for the Laplace equation, this requires the introduction of an interpolant, also of Morley type, of the finite volume solution possessing the desired conservation properties. For that purpose, we introduce new finite elements with appropriate degrees of freedom. In contrast with [27, 26] our interpolant is in H 1 (Ω), and therefore the residual error estimator is naturally based on the jump of normal derivatives of the interpolant of the solution. We finally show the equivalence between the energy norm of the error and the residual error estimator. The proof of the upper error bound uses a quasiorthogonality relation based on the conservation properties of the interpolant. The ∗ Received by the editors July 12, 2004; accepted for publication (in revised form) December 14, 2005; published electronically May 5, 2006. http://www.siam.org/journals/sinum/44-3/61148.html † MACS, ISTV, Universit´ e de Valenciennes et du Hainaut-Cambr´esis, F-59313 Valenciennes Cedex 9, France (
[email protected]).
949
950
SERGE NICAISE
proof of the lower error bound is more standard and simply uses some Green formulas and inverse inequalities as for finite element methods [34, 36, 15, 19]. In certain situations the solution of the diffusion-convection-reaction problem exhibits strong directional features. For instance, if the Peclet number (see below) is large, then boundary or interior layers may occur. Other examples are edge singularities appearing along concave edges in three-dimensional domains. In these cases the use of anisotropic meshes is recommended. Such meshes consist of elements in which the aspect ratio can be very large. Our second goal is to present a residual error estimate valid for anisotropic meshes satisfying minimal assumptions (contrary to the case of a standard residual error estimate for finite element methods; see [15, 16, 18]). This estimate is such that the size of the constant appearing in the upper bound is independent of the coefficients of the operator, while the size of the constant appearing in the lower bound is explicitly given with respect to the coefficients of the operator. These facts are further confirmed numerically. In contrast with standard practice [35, 19], we do not assume a strong coerciveness assumption (see below), and therefore our residual a posteriori analysis differs from [35, 19]. The outline of the paper is as follows: In section 2 we describe the so-called “cell centered” method for the diffusion-convection-reaction model problem on a triangulation of the domain consisting of triangles, rectangles, or tetrahedra. Some inverse inequalities are recalled in section 3, where we give some further useful interpolation error estimates. Section 4 is devoted to the introduction of some new finite elements (of Morley type) used later on. In section 5 we introduce the Morley interpolant of the approximate solution and prove its main properties. The upper and lower error bounds are then deduced in section 6. The upper error bound is based on the properties of the Morley interpolant, while the lower error bound is proved in a quite standard way. Finally, section 7 is devoted to numerical tests that confirm our theoretical considerations. 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 Γ. We consider the following standard elliptic problem: For f ∈ L2 (Ω) let u be the solution of Au := div (−ε∇u + vu) + bu = f in Ω, (2.1) u = 0 on Γ, where ε is a fixed positive constant, v a fixed vector function assumed to be suffi¯ Rd ), and b ∈ L∞ (Ω, R) is a fixed function. This ciently regular, namely v ∈ C 1 (Ω, problem is a linearized diffusion-convection-reaction problem appearing in many physical models. In the case of a large Peclet number P e ≡ ε−1 v∞ and/or large number Γ ≡ ε−1 div v + b∞ , 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 Dirichlet boundary condition. In this paper we further assume that 1 div v + b ≥ 0. 2
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
951
This assumption guarantees that there exists a unique weak solution u ∈ H01 (Ω) of problem (2.1), i.e., satisfying (2.2) (ε∇u · ∇v + div (vu)v + buv) dx = f v dx ∀v ∈ H01 (Ω). Ω
Ω
Moreover this solution satisfies (2.3)
|||u||| f ,
where ||| · ||| is the natural energy norm defined by 1 |||u|||2 := ε|∇u|2 + div v + b |u|2 . 2 Ω This norm slightly differs from the one used in [35, 19] since in those papers the authors assume that the strong coercivity assumption 1 div v(x) + b(x) ≥ c0 > 0 ∀x ∈ Ω, 2 holds. This condition excludes the study of convection-diffusion equations (2.1) with b = 0 and div v = 0 that we want to consider here. Therefore their residual a posteriori analysis slightly differs from the one presented here. 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 · s,D and | · |s,D . For the sake of brevity the L2 (D)-norm will be denoted by · D , and in the case when D = Ω, we will drop the index Ω. The space H01 (Ω) is defined, as usual, by H01 (Ω) := {v ∈ H 1 (Ω) : v = 0 on Γ}. In what follows, 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 meshsize of the triangulation, as well as the diffusion constant ε) such that a ≤ C b. The notation a ∼ b means that a b and b a hold simultaneously. 2.1. Discretization of the domain. To approximate problem (2.1) by a finite volume scheme we fix a mesh Th of Ω that satisfies the usual conformity conditions; cf. [6, Chap. 2]. In two dimensions we assume that all elements of Th are triangles or rectangles, while in three dimensions the mesh is made up of tetrahedra only. For K ∈ Th we denote by hK the diameter of K, and h = maxK∈Th hK . We define Eh as the set of edges (d = 2) or faces (d = 3) of the triangulation, Ehint = {E ∈ Eh /E ⊂ Ω} as the set of interior edges/faces of Th , and Ehext = Eh \Ehint as the set of exterior edges/faces of Th . For an edge E of a two-dimensional element K we introduce nK,E = (nx , ny ), which is the unit outward normal vector to K along E. Similarly, for a face E of a tetrahedron K let nK,E = (nx , ny , nz ) be the unit outward normal vector to K on E. Furthermore, for each edge/face E we fix one of the two normal vectors 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). As usual, let ωK be the union of all elements having a common edge/face with K. Similarly, let ωE be the union of both elements having E as an edge/face.
952
SERGE NICAISE
2.2. The finite volume scheme. Integrating (2.1) on a control volume K and using the divergence formula, we obtain (2.4) (−ε∇u + vu) · nK,E ds + bu dx = f (x) dx ∀K ∈ Th , E∈EK
E
K
K
where EK is the set of edges/faces of K. The continuous diffusion flux −ε∇u · nK,E is approximated using finite differences and the principle of conservation of flux, the expression vu · nK,E is approximated by a first order upwind scheme, and the reaction term K bu by a simple quadrature formula (see [10]). These approximations lead to the following system: Find uh := (uK )K∈Th (uK being the approximation of u(xK ) for K ∈ Th with xK being the “center” of the box K) that is a solution of
D R −εFK,E (2.5) (uh ) + vK,E FEC (uh ) + βK FK (uh ) = f (x) dx ∀K ∈ Th , K
E∈EK
where vK,E = ME (v · nK,E ), βK = MK b; MK g and ME g denote the mean of g on K and E, respectively, i.e., 1 1 MK g = g(x) dx ∀K ∈ Th , ME g = g(x) dσ(x) ∀E ∈ Eh , |K| K |E| E R (uh ) are defined as while the quantities FEC (uh ) and FK
(2.6)
FEC (uh ) := |E|uE,+ ,
where for E ∈ Ehint , uE,+ = uKE,+ , with KE,+ being the upstream control volume, ¯ ∩ Γ, uE,+ = uK if vK,E ≥ 0, and uE,+ = 0 i.e., vKE,+ ,E ≥ 0, while for E ∈ K otherwise. Similarly, R (uh ) = |K|uK . FK D (uh ), but the principle For our purposes, we do not need the exact form of FK,E of conservation of flux is required: D D ¯ ∩ L. ¯ (uh ) = −FL,E (uh ) for E = K FK,E
If the mesh Th is admissible in the sense of [10, Def. 9.1], i.e., satisfies standard orthogonality conditions (see Figure 2.1), then the numerical diffusion flux is defined by ⎧ ⎨ |E| (uL − uK ) if E = K ∩ L, dE D FK,E (uh ) := (2.7) ⎩ − |E| u if E ⊂ K ∩ ∂Ω, dE K when dE = d(xK , xL ) if E = K ∩L, with K, L ∈ Th and dE = d(xK , Γ) if E = ∂K ∩Γ. D (uh ) is proposed in [7, 8] using the For general meshes, a possible choice for FK,E diamond cell method. D For a restricted admissible mesh in the sense of [10, Def. 9.4], if FK,E (uh ) is given by (2.7), then system (2.5) is well defined as proved in [9]; see also [10] in the particular case when div v ≥ 0 and b a positive constant. For a general mesh, as is the case here, we simply assume that system (2.5) has a unique solution.
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
K
953
× xK
nL,E E
nK,E × xL L Fig. 2.1. The standard orthogonality condition.
P3
p3
P2 p2 p1
P0 P1 Fig. 2.2. Notation of tetrahedron K.
2.3. 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 [2, 15]. In this subsection we introduce and describe anisotropic quantities and present their basic properties. We start with an arbitrary (anisotropic) tetrahedron K ∈ Th . We enumerate its vertices so that P0 P1 is the longest edge, meas2 ( P0 P1 P2 ) ≥ meas2 ( P0 P1 P3 ), and meas1 (P1 P2 ) ≥ meas1 (P0 P2 ). Further, we introduce three orthogonal vectors pi,K of length hi,K := |pi,K |, as described in Figure 2.2. The minimal element size is particularly important; thus define hmin,K := h3,K . The three main anisotropic directions pi,K play an important role in several proofs. They span the orthogonal matrix CK := (p1,K , p2,K , p3,K ) ∈ R3×3 . This matrix may be considered as a transformation matrix which defines implicitly ˆ via K ˆ := C −1 (K − P0 ); cf. Figure 2.3. In order to facilitate the reference element K K the understanding of this mapping, the circumscribing box of K has been drawn in Figure 2.2. This box is mapped onto the unit cube given in Figure 2.3. Note in ˆ is of size O(1). particular that the reference element K
954
SERGE NICAISE
Pˆ3
Pˆ2
Pˆ0
Pˆ1
ˆ Fig. 2.3. Reference tetrahedron K.
In two dimensions the notation is similar. For a triangle K the enumeration is as in the bottom triangle P0 P1 P2 of Figure 2.2. For a rectangle K we simply take P0 P1 as the longest edge of K, P0 P2 as the edge perpendicular to P0 P1 , and simply take p1,K = P0 P1 , p2,K = P0 P2 . In both cases, hmin,K := h2,K , and CK becomes a 2 × 2 matrix. For an edge/face E of an element K, introduce the height hE,K = |K| |E| . 3. Analytical tools. 3.1. 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, bK and bE associated with an element K and an edge E, respectively. For a triangle or a tetrahedron K, denoting by λaK , i = 1, . . . , d + 1, the barycentric coordinates of K and by aE i , i = 1, . . . , d, the i vertices of the edge/face E ⊂ ∂K we recall that bK =
d+1
λaK and bE = i
i=1
d
λaE . i
i=1
For a rectangle K, we enumerate its vertices in a clockwise sense here. Denoting , i = 1, . . . , 4, the barycentric coordinates of K, namely, λaK is the unique by λaK i i K (a ) = δ element in Q1 (K) such that λaK , we recall that i,j j i bK = λaK λaK and bE = λaK (λaK + λaK ) 1 3 1 2 3 K if the endpoints of the edge E are aK 1 and a2 . We note that
bK = 0 on ∂K,
bE = 0 on ∂ωE ,
bK ∞,K = bE ∞,ωE ∼ 1.
In two dimensions for an edge E ⊂ ∂K, using temporarily the local coordinates system (x, y) such that E is included in the x-axis, the extension Fext (vE ) of vE ∈ C(E) to K is defined by Fext (vE )(x, y) = vE (x). We proceed similarly in three dimensions.
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
955
Now we may recall the so-called inverse inequalities that are proved using classical scaling techniques (cf. [34] for the isotropic case and [15] for the anisotropic case). Lemma 3.1 (inverse inequalities). Let vK ∈ Pk0 (K) and vE ∈ Pk1 (E) for some nonnegative integers k0 and k1 . Then the following inequalities hold, with the constants in the inequalities depending on the polynomial degrees k0 or k1 but not on K, E or vK , vE : 1/2
vK bK K ∼ vK K , ∇(vK bK )K h−1 min,K vK K ,
(3.1) (3.2)
1/2
(3.3)
vE bE E ∼ vE E ,
(3.4)
Fext (vE )bE K hE,K vE E ,
1/2
∇(Fext (vE )bE )K hE,K h−1 min,K vE E . 1/2
(3.5)
3.2. 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. However, these interpolation estimates do not hold for general 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, Th ) which was originally introduced in [16]. Definition 3.2 (alignment measure). Let v ∈ H 1 (Ω), and F = {Th } be a family of triangulations of Ω. Define the alignment measure m1 : H 1 (Ω) × F → R by 1/2 −2 2 m1 (v, Th ) := (3.6) ∇v. hmin,K CK ∇vK K∈Th
By definition one has m1 (v, Th ) ≥ 1. For arbitrary isotropic meshes one obtains that m1 (v, 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 error estimation. Since the focus of our work is a posteriori error estimation, we refer to [17, 16] for discussions concerning this alignment measure. Lemma 3.3 (local interpolation error bounds). Let v ∈ H01 (Ω); then, v − MK vK CK ∇vK
(3.7) (3.8)
hE,K v −
ME v2E
CK ∇v2K
∀K ∈ Th , ∀E ∈ EK , K ∈ Th .
Proof. The first inequality (3.7) has been proved in [16, Lemma 4]. The same ˆ into L2 (E) ˆ yield the second scaling argument and the compact embedding of H 1 (K) estimate. Lemma 3.4 (global interpolation error bounds). Let v ∈ H 1 (Ω); then, 2 2 2 (3.9) h−2 min,K v − MK vK m1 (v, Th ) ∇v , (3.10)
K∈Th
int K∈Th E∈EK ∩Eh
2 2 2 hE,K h−2 min,K v − ME vE m1 (v, Th ) ∇v .
956
SERGE NICAISE
Proof. These are direct consequences of the previous lemma and the definition of the alignment measure. For further analysis of our estimator, we need specific interpolation estimates related to the diffusion-convection-reaction problem [36, 18, 19]. Namely, the error estimates have to be related to the energy norm ||| · ||| and a local quantity relying on the local meshsize and the local behavior of the functions involved in the operator. This quantity is defined by −1/2
αK := min{cK
, ε−1/2 hmin,K },
where we set cK = min x∈K
1 div v(x) + b(x) . 2
Here we use the convention that if cK = 0, then the minimum is the second term, namely, αK := ε−1/2 hmin,K . We are now able to prove the following error estimate (compare with Lemma 3.9 of [18]). Lemma 3.5 (global interpolation error bounds). Let v ∈ H01 (Ω). Then we have
(3.11) ε
(3.12)
1/2
K∈Th
−2 αK v − MK v2K m1 (v, Th )2 |||v|||2 ,
K∈Th −1 αK
v − ME v2E m1 (v, Th )2 |||v|||2 .
int E∈EK ∩Eh
Proof. For the first estimate, we split up its left-hand side as follows:
−2 αK v − MK v2K =
cK v − MK v2K
K∈Th :εh−2 min,K ≤cK
K∈Th
+
2 εh−2 min,K v − MK vK .
K∈Th :εh−2 min,K >cK
By the definition of cK and estimate (3.9), we conclude that
−2 αK v − MK v2K
K∈Th
1 div v + b 2
1/2 v2 + εm1 (v, Th )2 ∇v2
m1 (v, Th )2 |||v|||2 . For the second estimate, we first use the same splitting
int K∈Th E∈EK ∩Eh
K∈Th : −2 εhmin,K ≤cK
int E∈EK ∩Eh
−1 αK v − ME v2E =
+
−1 αK v − ME v2E
K∈Th : −2 εhmin,K >cK
int E∈EK ∩Eh
−1 αK v − ME v2E .
957
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
As before, this leads to −1 αK v − ME v2E ε1/2 ε1/2 int K∈Th E∈EK ∩Eh
+ε
1/2
cK
K∈Th : −2 εhmin,K ≤cK
h−1 min,K
v2E
int E∈EK ∩Eh
v − ME v2E .
int E∈EK ∩Eh
K∈Th : −2 εhmin,K >cK
Using the estimate hmin,K hE,K proved in Lemma 3.1 of [18] and the estimate (3.10), we obtain 1/2 −1 (3.13) ε1/2 αK v − ME v2E ε1/2 cK v2E int K∈Th E∈EK ∩Eh
K∈Th : −2 εhmin,K ≤cK
int E∈EK ∩Eh
+ εm1 (v, Th )2 ∇v2 . For the first term of this right-hand side, using the trace inequality (see, for instance, Lemma 2.4 of [15] or Lemma 3.5 of [18]) v2E h−1 E,K vK (vK + CK ∇vK ),
we may write 1/2 ε1/2 cK
int E∈EK ∩Eh
K∈Th : −2 εhmin,K ≤cK
v2E
2 ε1/2 cK h−1 E,K (vK + vK CK ∇vK ). 1/2
K∈Th : −2 εhmin,K ≤cK
(3.14) Fix for a moment an element K such that εh−2 min,K ≤ cK . Using the property hmin,K hE,K and Young’s inequality with a parameter ηK > 0, we may write 2 1/2 2 cK h−1 ε1/2 cK h−1 E,K (vK + vK CK ∇vK ) ε min,K vK 1/2
1/2
1/2 ε1/2 cK h−2 ε1/2 cK min,K ηK 2 CK + vK + ∇v2K . 2ηK 2 1/2
Since ε1/2 h−1 min,K ≤ cK , the first term has the correct factor, namely, 1/2
2 ε1/2 cK h−1 E,K (vK + vK CK ∇vK ) 1/2
1/2 ε1/2 cK h−2 ε1/2 cK min,K ηK 2 CK + vK + ∇v2K . 2ηK 2 1/2
cK v2K
For the last two terms we choose ηK =
ε1/2 1/2 cK
which yields
−2 2 2 2 ε1/2 cK h−1 E,K (vK + vK CK ∇vK ) cK vK + εhmin,K CK ∇vK . 1/2
This estimate in (3.14) leads to 1/2 2 ε1/2 cK v2E (cK v2K + εh−2 min,K CK ∇vK ) K∈Th : −2 εhmin,K ≤cK
int E∈EK ∩Eh
K∈Th
m1 (v, Th )2 |||v|||2 . Inserting this estimate in (3.13) we arrive at (3.12).
958
SERGE NICAISE
4. Some finite elements of Morley type. For our further analysis, we need a continuous interpolant p satisfying ∂p D (4.1) ds = FK,E (uh ) ∀E ∈ EK , ∂n K,E E C v · nK,E p ds = vK,E FK,E (uh ) ∀E ∈ EK , (4.2) E
R bp = βK FK (uh )
(4.3) K
for any K ∈ Th . This means that we need to use C 0 -finite elements having as degrees of freedom the three left-hand sides above. Since these left-hand sides clearly depend on the restriction of v and b on K, the finite elements depend on these restrictions. The interpolant’s construction is quite complicated and could be expansive from a computational point of view. Such elements are even not unique. We therefore build them in a generic way and also give simpler elements in the case of constant coefficients. 4.1. Convection-diffusion elements for constant coefficients. For a convection-diffusion problem with constant coefficients (i.e., for b = 0 and v ∈ Rd , v = 0) we need an interpolant p satisfying only (4.1) and (4.2). In this case we need finite elements having as degrees of freedom the mean of p and the normal derivative of p on each edge. Therefore we need a kind of Morley element [24, 6]. In two dimensions we modify the element of Nilssen, Tai, and Winther [28] and extend this new element to rectangles and to tetrahedra. 4.1.1. Triangles. Here K is a (nondegenerate) triangle with vertices aK i , i = 1, 2, Nf := 3. Now, inspired by [28, sect. 4], we take
(4.4)
PK = P2 (K) ⊕ P1 (K)bK = {q + pbK : q ∈ P2 (K), p ∈ P1 (K)}, ∂p K ΣK = {p(ai )}i=1,...,Nf ∪ p ds ∪ ds . E E ∂nK,E E∈EK E∈EK
Lemma 4.1. The above triple (K, PK , ΣK ) is a C 0 -finite element. Proof. It suffices to show that w ∈ PK satisfying (4.5)
l(w) = 0 ∀l ∈ ΣK
is equal to zero. But w is of the form w = q + pbK with q ∈ P2 (K), p ∈ P1 (K), and K since bK is identically equal to zero on the edges E ∈ EK , the conditions w(ai ) = w ds = 0 are equivalent to E K q(ai ) = q ds = 0. E
Since the restriction of q to E is of degree 2, these conditions imply that q|E ≡ 0 for all E ∈ EK and consequently q ≡ 0. The conclusion follows from Lemma 4.1 of [28] since it is proved there that w = pbK with p ∈ P1 (K) satisfying E ∂n∂w ds = 0 for all E ∈ EK is identically equal to K,E zero. The continuity of the element follows from the first part of the proof.
959
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
4.1.2. Rectangles. Here K is a (nondegenerate) rectangle with vertices aK i , i = 1, . . . , Nf := 4 with edges parallel to the axes. Now we take PK = {q + pbK : q ∈ Q2 (K), p(x, y) = αx + βy + γy 2 , α, β, γ ∈ R}. As before, the degrees of freedom are defined by (4.4). Lemma 4.2. The above triple (K, PK , ΣK ) is a C 0 -finite element. Proof. By an affine transformation, it suffices to show the result on the reference ˆ = (0, 1)2 . For the sake of brevity we will omit the signˆ. As dim PK = card element K ΣK , it suffices to show that w ∈ PK satisfying (4.5) is equal to zero. But w is of the form w = q+pbK with q ∈ Q2 (K) and p of the form p(x, y) = αx+βy+γy 2 , α, β, γ ∈ R. Since bK is identically equal to zero on the edges E ∈ EK , the conditions w(aK i ) = w ds = 0 are equivalent to E q(aK ) = q ds = 0. i E
Since the restriction of q to E is of degree 2, these conditions imply, as before, that q|E ≡ 0 for all E ∈ EK . Fix the basis {λi (x)λj (y)}0≤i,j≤2 of Q2 (K), where λi ∈ P2 (0, 1) satisfy λi (xj ) = δij , where xi = 2i , i = 0, 1, 2. Writing q in this basis, we observe that the property q|E ≡ 0 for all E ∈ EK implies that q(x, y) = δλ1 (x)λ1 (y) = δbK (x, y) for some δ ∈ R. Returning to w, we conclude that w(x, y) = (δ + αx + βy + γy 2 )bK (x, y). Now, one readily sees that the remainding condition E ∂n∂w ds = 0 for all E ∈ EK K,E is equivalent to a homogeneous 4 × 4 linear system in α, β, γ, δ whose unique solution is α = β = γ = δ = 0. Therefore w is identically equal to zero. 4.1.3. Tetrahedra. Here K is a (nondegenerate) tetrahedron with vertices aK i , i = 1, 2, 3, Nf := 4. Inspired by the above triangular example, we choose ΣK defined by (4.4) and PK = q + pbK + αE bE : p, q ∈ P1 (K), αE ∈ R . E∈EK
As before, we can prove that the above triple (K, PK , ΣK ) is a C 0 -finite element. 4.2. Diffusion-convection-reaction elements for constant coefficients. We build our finite elements by slightly modifying the elements from the previous subsection. 4.2.1. Triangles. We take PK = {q + (p + αbK )bK : q ∈ P2 (K), p ∈ P1 (K), α ∈ R}, K (4.6) ΣK = {p(ai )}i=1,...,Nf ∪ p ∪ p ds ∪ K
E
E∈EK
E
∂p ds ∂nK,E
. E∈EK
Lemma 4.3. The above triple (K, PK , ΣK ) is a C 0 -finite element. Proof. It suffices to show that w ∈ PK satisfying (4.5) is equal to zero. Here w is of the form w = q + (p + αbK )bK with q ∈ P2 (K), p ∈ P1 (K), and α ∈ R. As in Lemma 4.1, the conditions w(aK i ) = E w ds = 0 imply that q ≡ 0.
960
SERGE NICAISE
Next, we consider the condition ∂w ds = 0 ∂n K,E E which is reduced to
E
∀E ∈ EK ,
∂(pbK ) ds = 0 ∂nK,E
∀E ∈ EK .
Thanks to Lemma 4.1 of [28], p is identically equal to zero. Finally, the condition K w becomes α K b2K = 0, which leads to the conclusion. 4.2.2. Rectangles. In the setting of subsection 4.1.2, we take ΣK defined by (4.6) and PK = {q + pbK : q ∈ Q2 (K), p(x, y) = αx + βy + γx2 + δy 2 , α, β, γ, δ ∈ R}. Arguments similar to those used in the proof of Lemma 4.2 allow one to show that the above triple (K, PK , ΣK ) is a C 0 -finite element. 4.2.3. Tetrahedra. Inspired by the above triangular example and the tetrahedral example from subsection 4.1, we choose ΣK defined by (4.6) and 2 PK = q + pbK + αbK + αE bE : p, q ∈ P1 (K), αE , α ∈ R . E∈EK
As before, the above triple (K, PK , ΣK ) is a C 0 -finite element. 4.3. Reaction-diffusion elements for constant coefficients. If we consider only reaction-diffusion equations (i.e., v = 0), we may restrict ourselves to an interpolant p satisfying (4.1) and (4.3). Therefore we need to use finite elements having as degrees of freedom the left-hand side of (4.1) and the mean on K. Again we need to slightly modify the previous elements. 4.3.1. Triangles or tetrahedra. We take
(4.7)
PK = {q + (p + αbK )bK : q ∈ P1 (K), p ∈ P1 (K), α ∈ R}, ∂p K ΣK = {p(ai )}i=1,...,Nf ∪ p ∪ ds . K E ∂nK,E E∈EK
Lemma 4.4. The above triple (K, PK , ΣK ) is a C 0 -finite element. Proof. As usual, it suffices to show that w ∈ PK satisfying (4.5) is equal to zero. Here w is of the form w = q + (p + αbK )bK with q ∈ P1 (K), p ∈ P1 (K), and α ∈ R. As in Lemma 4.1, the conditions w(aK i ) = 0 imply that q ≡ 0 since q is of degree at most 1. The remainder of the proof is the same as the one of Lemma 4.3. 4.3.2. Rectangles. In the setting of subsection 4.1.2, we take ΣK defined by (4.7) and PK = {q + pbK : q ∈ Q1 (K), p(x, y) = η + αx + βy + γx2 + δy 2 , α, β, γ, δ, η ∈ R}. Arguments similar to those used in the proof of Lemma 4.2 allow one to show that the triple (K, PK , ΣK ) is a C 0 -finite element.
961
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
4.4. General case. 4.4.1. Triangles or tetrahedra. Inspired by the previous examples, we take (4.8)
v,b PK =
q0 +
αE v˜K,E bE bK
E∈EK
+
βE bE + γbK bK : q0 ∈ P1 (K), αE , βE , γ ∈ R ,
E∈EK
(4.9)
Σv,b K
=
{p(aK i )}i=1,...,Nf ∪ E
∪
∂p ds ∂nK,E
v · nK,E p ds
E
∪
E∈EK
E∈EK
bp ,
K
where v˜K,E is any extension from E to K of the function v · nK,E such that v˜K,E ≡ 0 if v · nK,E ≡ 0 on E. Here there is a slight abuse of notation in the sense that the degree of freedom E v · nK,E p ds (or K bp) disappears if v · nK,E ≡ 0 on E (or b ≡ 0). v,b 0 As before, we can prove that the triple (K, PK , Σv,b K ) is a C -finite element. v,b defined by (4.8), but with q0 in Q1 (K) 4.4.2. Rectangles. Now we take PK v,b 0 defined by (4.9). Again this triple (K, PK , Σv,b and Σv,b K K ) is a C -finite element.
5. The Morley interpolant. v,b 5.1. Definition. For any K ∈ Th , we fix a C 0 -finite element (K, PK , Σv,b K ) such that ∂p v · nK,E p ds ∪ ds ∪ bp ⊂ Σv,b K . E E ∂nK,E K E∈EK E∈EK
We refer to the previous section for its existence. We now introduce the finite element space: Vh :=
v,b vh ∈H01 (Ω) : vh|K ∈ PK ∀K ∈ Th , ∂vh|K ∂vh|K ds = ds ∀E ∈ Eh , K, L ∈ Th : E = K ∩ L . E ∂nE E ∂nE
Definition 5.1. For uh = (uK )K∈Th we define its interpolant (of Morley type) IM uh as the unique element vh in Vh satisfying (5.1) E
∂vh|K D ds = FK,E (uh ) ∂nK,E
∀E ∈ EK , K ∈ Th ,
v · nK,E vh ds = vK,E FEC (uh )
(5.2)
∀E ∈ Ehint ,
E
R bvh dx = βK FK (uh )
(5.3) K
∀K ∈ Th .
962
SERGE NICAISE
5.2. Some useful properties. The key point of our a posteriori analysis is the following basic property of the Morley interpolant. Lemma 5.2. If uh is solution of (2.5), then IM uh satisfies (5.4) (A(IM uh ) − f ) dx = 0 ∀K ∈ Th . K
Proof. By Green’s formula we have A(IM uh ) dx = (div (−ε∇IM uh + vIM uh ) + bIM uh ) dx K K ∂(IM uh ) = −ε + v · nK,E IM uh ds + bIM uh dx. ∂nK,E E K E∈EK
Using the properties (5.1), (5.2), and (5.3) satisfied by IM uh we obtain
D R A(IM uh ) dx = (uh ) + vK,E FEC (uh ) + βK FK (uh ) −εFK,E K
E∈EK
and we conclude by using (2.5). Now we prove a quasi-orthogonality relation that will be used for the upper error bound. We first define the gradient jump of IM uh in the normal direction by ∂ JE,n (uh ) = ε (IM uh ) ∀E ∈ Ehint . ∂nE E Lemma 5.3. If u is a solution of (2.2) and uh is a solution of (2.5), then for any χ ∈ H01 (Ω), setting e = u − IM uh , we have that (ε∇e · ∇χ + (div (ve) + be)χ) dx = (f − A(IM uh ))(χ − MK χ) dx (5.5) Ω
K∈Th
+
K
int E∈Eh
JE,n (uh )(χ − ME χ) ds.
E
Proof. For the sake of brevity denote the left-hand side of (5.5) by I1 (χ). By (2.2) and using Green’s formula on each element K we obtain ∂(IM uh ) I1 (χ) = f χ dx − A(IM uh )χ dx − ε χ ds. ∂nK Ω K ∂K K∈Th
K∈Th
The continuity of IM uh and of χ across the edges/faces (in the sense of trace) and the fact that χ = 0 on Γ lead to I1 (χ) = (f − A(IM uh ))χ dx + JE,n (uh )χ ds. K∈Th
K
int E∈Eh
E
Using the identity (5.4) we arrive at I1 (χ) = (f − A(IM uh ))(χ − MK χ) dx + JE,n (uh )χ ds. K∈Th
K
int E∈Eh
E
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
963
The conclusion now follows from the fact that JE,n (uh ) ds = 0 ∀E ∈ Ehint , E D D which is due to (5.1) and the principle of conservation of flux: FK,E (uh ) = −FL,E (uh ) if E = K ∩ L, K, L ∈ Th . Corollary 5.4. Let the assumptions of Lemma 5.3 be satisfied. Then, the following estimate holds: 2 αK (5.6) |I1 (χ)| f − A(IM uh )2K K∈Th
+ε
−1/2
αK
1/2 JE,n (uh )2E
m1 (χ, Th )|||χ|||.
int E∈EK ∩Eh
Proof. The identity (5.5) and the Cauchy–Schwarz inequality yield |I1 (χ)| f − A(IM uh )K χ − MK χK K∈Th
+
JE,n (uh )E χ − ME χE .
int K∈Th E∈EK ∩Eh
The discrete Cauchy–Schwarz inequality then leads to 1/2 1/2 −2 2 |I1 (χ)| αK f − A(IM uh )2K αK χ − MK χ2K K∈Th
+
K∈Th
ε−1/2 αK
K∈Th
·
−1 ε1/2 αK
K∈Th
int E∈EK ∩Eh
1/2
JE,n (uh )2E χ −
1/2 ME χ2E
.
int E∈EK ∩Eh
We conclude by applying Lemma 3.5. Remark 5.5. The fundamental properties above are based only on the definition of the scheme (2.5), on the continuity of the interpolant, and on the interpolation properties (5.1), (5.2), and (5.3). Therefore our further analysis is valid for any finite v,b , Σv,b element (K, PK K ) such that the associated interpolant satisfies these properties. But the finite element and the definition of the interpolant should be appropriately chosen in order to guarantee the convergence of IM uh to the exact solution u. This convergence analysis is yet to be performed but, in any case, it is outside the scope of this paper. 6. Error estimators. 6.1. Residual error estimators. The exact element residual is defined by RK := f − A(IM uh ) on K. As usual, this exact residual is replaced by a certain finite-dimensional approximation rK ∈ Pk (K) called an approximate element residual. Definition 6.1 (residual error estimator). The local and global residual error estimators are defined by 2 2 2 := αK rK 2K + ε−1/2 αK JE,n (uh )2E , η 2 := ηK . ηK int E∈EK ∩Eh
K∈Th
964
SERGE NICAISE
The local and global approximation terms are defined by 2 2 2 ζK := αK RK − rK 2K , ζ 2 := ζK . K ⊂ωK
K∈Th
6.2. Upper error bound. Theorem 6.2. Let u be a solution of (2.2) and uh a solution of (2.5), and denote the error by e := u − IM uh . Then the error is bounded as follows: |||e||| m1 (e, Th )(η + ζ).
(6.1)
Proof. As e belongs to H01 (Ω), Green’s formula yields 1 2 div v + b |e|2 dx. (div (ve)e + b|e| ) dx = 2 Ω Ω Therefore we have
|||e|||2 =
(ε∇e · ∇e + (div (ve) + be)e) dx, Ω
or equivalently, with the notation from Lemma 5.3, |||e|||2 = I1 (e). The conclusion immediately follows from estimate (5.6). 6.3. Lower error bound. Theorem 6.3. Assume that there exists δ1 > 2 and δ2 ∈ [3/2, 2) such that −div v ≤ δ1 b in {x ∈ Ω : div v(x) > 0}, (6.2) −div v ≤ δ2 b in {x ∈ Ω : div v(x) < 0}. Then for all elements K, the following local lower error bound holds: (6.3)
ηK (1 + P eωK + ΓωK )|||e|||ωK + ζK ,
where |||e|||2ω := ω (ε|∇e|2 + ( 12 div v + b)|e|2 ); P eωK = maxK ⊂ωK P eK is the local patch Peclet number, with the local mesh Peclet number being defined as usual by P eK = ε−1 v∞,K hmin,K ; see [11, 2, 19] for anisotropic meshes and [32, 35] for isotropic meshes. The elementwise quantity ΓωK is defined similarly; namely, ΓωK = maxK ⊂ωK ΓK , where (cf. [3, 2]) 1/2
ΓK = max{1, αK b + div v∞,K }. Proof. Since b + 12 div v ≥ 0, we readily see that (6.2) is equivalent to 1 |b + div v| ≤ γ b + div v in Ω, (6.4) 2 δ2 −1 with γ = 2 max{ δδ11 −1 −2 , 2−δ2 } > 0.
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
965
Element residual. For a fixed element K define wK = rK bK , which belongs to H01 (Ω). From the definition of RK we have rK wK = (rK − RK )wK + RK wK = (rK − RK )wK + A(u − IM uh )wK K K K K K = (rK − RK )wK + (div (ve) + be)wK − div (ε∇e)wK . K
K
K
Integrating by parts in that last term we obtain rK wK = (rK − RK )wK + (ε∇e · ∇wK + (div (ve) + be)wK ) . K
K
K
Leibniz’s rule and the Cauchy–Schwarz inequality yield rK wK = rK − RK K wK K + (ε∇wK K + v∞,K wK K )∇eK K
+ (div v + b)eK wK K . Using property (6.4), we obtain rK wK ≤ rK − RK K wK K + (ε∇wK K + v∞,K wK K )∇eK K
+γ
1/2
div v +
1/2 b∞,K
1 div v + b 2
1/2 eK wK K .
By the inverse inequalities (3.1), (3.2) we get αK rK K ξK + αK (εh−1 min,K + v∞,K )∇eK 1/2 1 1/2 + αK div v + b∞,K div v + b eK . 2 Using the property αK ≤ ε−1/2 hmin,K , we conclude that (6.5)
αK rK K (1 + P eK + ΓK )|||e|||K + ζK .
Normal jump. Fix an arbitrary edge/face E ∈ Ehint . Recall that JE,n (uh ) ∈ Pk (E) for some k ∈ N and let wE := Fext (JE,n (uh ))bE,γE,K on KE,γE,K ⊂ K ⊂ ωE , where KE,γE,K is the squeezed element associated with K (defined in [36, 15, 19] for triangles and tetrahedra and easily extended to rectangles; we do not go into detail here, though we note that the main properties of KE,γE,K are that it be included in K to have E as edge/face and be of height ∼ γE,K hE,K ) and the parameter γE,K is fixed and is equal to hmin,K ε1/2 γE,K = min 1, . , hE,K c1/2 hE,K K Here the function bE,γE,K is the edge/face bubble associated with E in KE,γE,K . The difference with the choice made in [19] relies on the factor cK .
966
SERGE NICAISE
By elementwise partial integration we get ∂e JE,n (uh )wE = −ε wE = − (ε∇e∇wE + div (ε∇e)wE ) ∂nK E K⊂ωE ∂K K⊂ωE K (−ε∇e∇wE + AewE − (div (ve) + be)wE ) . =
K⊂ωE
K
By Leibniz’s rule, the Cauchy–Schwarz inequality, and condition (6.4) we obtain RK K wE K + (ε∇wE K + v∞,K wE K )∇eK JE,n (uh )wE E
K⊂ωE
1/2
+ div v + b∞,K
1 div v + b 2
1/2
eK wE K .
The inverse inequalities from Lemma 2 of [19] (see our Lemma 3.1 and the above properties of KE,γE,K ) in the previous estimate lead to 1/2 1/2 γE,K hE,K [(ε min{γE,K hE,K , hmin,K }−1 + v∞,K )∇eK JE,n (uh )E K⊂ωE
+ div v +
1/2 b∞,K
1 div v + b 2
1/2 eK + RK K ].
Multiplying this estimate by ε−1/4 αK and using the definition of γE,K we arrive at 1/2
(6.6) 1/2 ε−1/4 αK JE,n (uh )E
(1 + P eωE )ε
1/2
∇eωE
1 + ΓωE div v + b 2
1/2 eωE + ζK .
The conclusion follows from estimates (6.5) and (6.6). Remark 6.4. Condition (6.2) is not restrictive since it is satisfied in the so-called convection-dominated case b + 12 div v ≥ c0 > 0 in Ω. It also holds if div v + b ≥ 0 and b ≥ 0 (a.e.) in Ω. 7. Numerical results. In this section we present two examples that illustrate the efficiency and reliability of our estimator and show that our estimator is appropriate for adaptivity. Additionally, for both examples we provide the order of convergence of the error |||u − IM uh |||. In both cases the order is approximately h, which confirms that IM uh is a good approximation to u. The first and second examples concern problems with solutions which exhibit boundary layers and for which the use of anisotropic meshes is recommended. The third example even treats the case when the meshes are not fully aligned with the solution. Finally, for a convection-diffusion problem we present two sequences of meshes obtained by an adaptive process. As the meshes are refined in the region of large variation of the solutions we can verify the reliability of our estimator. For both examples we investigate the main theoretical results which are the upper and lower error bounds. In order to present the inequalities (6.1) and (6.3) appropriately, we consider the ratios qup :=
|||u − IM uh ||| , η+ξ
qlow := max
K∈Th
(1 + P eωK
ηK + ΓωK )|||e|||ωK + ζK
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
967
as a function of the degrees of freedom. The first ratio qup is frequently referred to as the effectivity index. It measures the reliability of the estimator and is related to the global upper bound on the error. The second ratio is related to the local lower error bound and measures the efficiency of the estimator. From our theoretical considerations, both ratios should be bounded from above, which is confirmed below experimentally. Hence our estimator is reliable and efficient. 7.1. A reaction-diffusion problem (Test 1). As a particular problem, in (2.1) we take v = (0, 0) and b = 1 in the unit square; i.e., we consider −ε Δu + u = 0
in
Ω = (0, 1)2 .
We prescribe the exact solution to be √ u(x, y) = exp (−x/ ε), with the Dirichlet boundary datum being √fixed accordingly. This solution exhibits an exponential boundary layer of width O( ε| ln ε|) along the y-axis. To approximate the solution to this problem appropriately we use anisotropic meshes. Each mesh is the tensor product of a Shishkin-type mesh in the x-direction and a uniform mesh in the y-direction,√both with n subintervals. More precisely we fix the transition point τ := min{1/2, ε| ln ε|} and define the rectangular mesh of nodes (xi , yj ), 0 ≤ i, j ≤ n, with (n is assumed to be even) xi = i2τ h for 0 ≤ i ≤ n/2, xi = τ + (2i − n)(1 − τ )h for n/2 ≤ i ≤ n, yj = jh for 0 ≤ j ≤ n, where the meshsize h is equal to 1/n. Such a mesh is illustrated in Figure 7.1 (left) for h = 1/8 and τ ≈ 0.25. For each triangulation, we compute the finite volume approximation uh as a solution to (2.5), and then compute its interpolant IM uh using the finite element of subsection 4.3.2 (with the values of IM uh at any interior node a being fixed by IM uh (a) = 14 a∈K uK ). Figure 7.2 presents the energy norm |||u − IM uh ||| with respect to n for different values of ε. The rate of convergence is approximately h, which confirms that IM uh converges quite well to u. Now we investigate the upper and lower error bounds. According to previous results on finite element methods [15, 19], the meshes are appropriately chosen to resolve the boundary layer, and consequently the alignment measure m1 (e, Th ) should be moderated. Figure 7.3 presents the values of qup (top) and of qlow (bottom) with respect to n for different values of ε. We observe that qup is bounded from above by 0.15 and that qlow is bounded from above by 0.5. 7.2. Convection-diffusion problems (Test 2). Here we consider the convectiondiffusion problem (7.1)
−ε Δu + div (vu) = 0 in Ω = (0, 1)2 ,
where v is either (−1, 0) or (−1, −1) . We prescribe the exact solution to be, respectively, u(x, y) = exp (−x/ε) or u(x, y) = exp (−(x + y)/ε),
968
SERGE NICAISE
1
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
Fig. 7.1. Shishkin-type meshes on the unit square with n = 8. −3 ε = 10−4 ε = 10−6 ε = 10−8 ε = 10−10 ε = 10−12
−4
| ln( | u − I_M u_h |1, h ) |
−5
−6
−7
−8
−9 1 −10 1 −11
2
2.5
3
3.5
4 | ln(n) |
4.5
5
5.5
6
Fig. 7.2. Illustration of the convergence of ∇h (u − IM uh ) for Test 1.
with the Dirichlet boundary datum being fixed accordingly. In the first case the solution exhibits an exponential boundary layer of width O(ε| ln ε|) along the y-axis, while in the second case, two exponential boundary layers appear along the x and y axes. Therefore these problems are approximated using similar anisotropic meshes as before with the transition point τ := min{1/2, ε| ln ε|} (with the tensor product of two Shishkin type meshes in the second case; see Figure 7.1 (right)). Once we have computed the finite volume approximation uh , solution of (2.5), we compute its interpolant IM uh using the finite element of subsection 4.1.2 and the same values at the nodes as before. Figure 7.4 shows the energy norm |||u − IM uh ||| with respect to n for different values of ε and confirms a rate of convergence of approximately h. As before, the meshes are appropriately chosen so the alignment measure m1 (e, Th ) should be moderated. The values of qup and of qlow with respect to n for different values of ε are plotted in Figures 7.5 and 7.6. Here we observe similar values for both examples; moreover we see that qup is bounded from above by 0.1 and that qlow is approximately 2.2. 7.3. Meshes not aligned with the solutions (Test 3). Here we consider the convection-diffusion problem (7.1) with v = (− cos θ, sin θ), with θ = 0, 0.1π, and 0.2π, and its prescribed solution given by u(x, y) = exp ((− cos θx + sin θy)/ε)
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
969
ε = 10−4 −6 ε = 10 −8 ε = 10 −10 ε = 10 ε = 10−12
0.25
q
up
0.2
0.15
0.1
0.05
0
0
0.8
50
100
150 n
200
250
300
150 n
200
250
300
ε = 10−4 −6 ε = 10 −8 ε = 10 ε = 10−10 ε = 10−12
0.7
0.6
q
low
0.5
0.4
0.3
0.2
0.1
0
0
50
100
Fig. 7.3. qup (top) and qlow (bottom) with respect to n for Test 1.
and ε = 0.05. This solution has a layer along the line − cos θx + sin θy = 0 (the case θ = 0 corresponds to the case of the previous subsection but is presented in order to compare the different results). Here we take the sequence of meshes from subsection 7.1 with the transition point τ = 4ε = 0.2. This means that for θ = 0.1π and 0.2π, the meshes are not fully aligned with the solution and are less and less aligned as θ increases. Nevertheless the convergence rate is h according to Figure 7.7. Moreover the values of qup and qlow with respect to n for the different values of θ are plotted in Figures 7.8, where we see that qup varies between 0.1 and 0.14 and that qlow is bounded from above by 2.75. From this figure, we further remark that the effectivity index qup grows as θ increases, implying that the matching function m1 (e, Th ) grows as well. This phenomenon was expected since the meshes are less and less aligned as θ increases. 7.4. An adaptive algorithm (Test 4). In view of the upper bound (6.1), we use the following (standard) adaptive process: Starting from an initial mesh Th0 , we
970
SERGE NICAISE –2 01 ε = 10– ε = 10–03 05 – ε = 10
–4
–6
ln error
–8
–10
1 –12 1
–14
–16
1
1.5
2
2.5
3
3.5
4
4.5
log(n)
–2 ε = 10–01 03 ε = 10– –05 ε = 10
–4
–6
ln error
–8
–10
1 –12 1 –14
–16
1
1.5
2
2.5
3
3.5
4
4.5
log(n)
Fig. 7.4. Illustration of the convergence of ∇h (u − IM uh ) for Test 2 for v = (−1, 0) (top) and for v = (−1, −1) (bottom).
mark the elements K for which ηK , ηK > c max K
for a chosen constant c ∈ (0, 1). All marked elements are divided into four subelements (standard regular refinement rule), with the other elements being subdivided only to guarantee the conformity of the new mesh. We refine the meshes up to the requested accuracy. We test this adaptive algorithm for the convection-diffusion problem (7.1) with v = (1, 1) and the prescribed solution given by u(x, y) = exp (−100((x − 0.5)2 + (y − 0.5)2 )), with the right-hand side and the boundary datum being fixed accordingly. This exact solution is a Gaussian function whose center is the point (0.5, 0.5). For this example we use either meshes consisting of rectangles satisfying the admissibility condition or
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
971
0.13 01 ε = 10– ε = 10–03 05 – ε = 10
0.125
0.12
q
up
0.115
0.11
0.105
0.1
0.095
1
1.5
2
2.5
3
3.5
4
4.5
log(n)
2.35
ε = 10–01 03 ε = 10– ε = 10–05
2.3
2.25
2.2
q low
2.15
2.1
2.05
2
1.95
1.9
1.85
1
1.5
2
2.5
3
3.5
4
4.5
log(n)
Fig. 7.5. qup (top) and qlow (bottom) with respect to n for Test 2 for v = (−1, 0) .
meshes consisting of triangles that do not satisfy the admissibility condition (in which case we use the diamond cell technique; see subsection 2.2). The Morley interpolant is computed using the finite elements from subsection 4.1.2 or 4.1.1. The meshes obtained after four iterations are shown in Figure 7.9. In the case of rectangular meshes, anisotropic meshes appear due to the admissibility condition, but this does not affect the convergence of our adaptive process. This also underlines the incompatibility between standard refinement rules and this admissibility condition. From this figure, we can conclude that the meshes are refined in the region of large variation of the solution. Again this confirms the reliability of our estimator. Figure 7.10 shows the energy norm |||u − IM uh ||| with respect to the degrees of freedom for the adaptive algorithm in comparison with uniformly refined meshes. In both cases, we see that the adaptive algorithm gives rise to better error bounds and to convergence in h. The values of qup and qlow with respect to the degrees of freedom are plotted in Figures 7.11 and 7.12 and are compared with uniformly refined meshes. There we observe that qup is bounded from above by 0.1 and that qlow is bounded from above
972
SERGE NICAISE 0.13 01 ε = 10– ε = 10–03 05 – ε = 10
0.125
0.12
q
up
0.115
0.11
0.105
0.1
0.095
1
1.5
2
2.5
3
3.5
4
4.5
log(n)
2.5 ε = 10–01 ε = 10–03 –05 ε = 10
2.4
2.3
q
low
2.2
2.1
2
1.9
1.8
1
1.5
2
2.5
3
3.5
4
4.5
log(n)
Fig. 7.6. qup (top) and qlow (bottom) with respect to n for Test 2 for v = (−1, −1) .
by 1.42 (resp., 3.5). From these figures, we can say that the different values of qup remain similar for adaptive meshes and for uniformly refined meshes, while the values of qlow have a different behavior but stay in a relatively small range of variations. For rectangular meshes, uniform meshes lead to smallest values of qlow , probably due to the overrefinement of adaptive meshes, while the converse holds for triangular meshes because adaptive triangular meshes fit the solution well. 8. Conclusions. We have proposed and rigorously analyzed a new a posteriori error estimate for a cell centered finite volume approximation of diffusion-convectionreaction equations. We have shown that this estimate is reliable and efficient. This estimate is based on the construction of an appropriate interpolant of Morley type. Some numerical experiments confirm our theoretical predictions. Acknowledgment. I am very grateful to Karim Djadel (Cermics, ENPC, France), who made the numerical experiments presented in section 7.
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS –1 θ = 0. θ = 0.1 π θ = 0.2 π
–1.5
–2
–2.5
ln error
–3
–3.5
–4
–4.5
–5
1
–5.5 1 –6
2
2.5
3
3.5 log(n)
4
4.5
5
Fig. 7.7. Illustration of the convergence of ∇h (u − IM uh ) for Test 3. 0.16 θ = 0. θ = 0.1 π θ = 0.2 π 0.15
q
up
0.14
0.13
0.12
0.11
0.1
2
2.5
3
3.5 log(n)
4
4.5
5
5 θ = 0. θ = 0.1 π θ = 0.2 π 4.5
4
q
low
3.5
3
2.5
2
1.5
2
2.5
3
3.5 log(n)
4
4.5
5
Fig. 7.8. qup (top) and qlow (bottom) with respect to n for Test 3.
973
974
1 0.8 0.6 0.4 0.2 0 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
SERGE NICAISE
Fig. 7.9. The last meshes in the refinement sequences for Test 4. 0 adaptative uniform
–0.5
ln error
–1
–1.5
–2 1/2 –2.5 1
–3
5
6
7
8
9
10
11
12
log(DoF)
0 adaptative uniform –0.5
–1
ln error
–1.5
–2
–2.5 1/2 –3 1 –3.5
5
6
7
8
9
10
11
12
log(DoF)
Fig. 7.10. Illustration of the convergence of ∇h (u − IM uh ) for Test 4 for rectangular meshes (top) and for triangular meshes (bottom).
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
975
0.106 adaptative uniform
0.104
q up
0.102
0.1
0.098
0.096
0.094
5
6
7
8
9
10
11
12
log(DoF)
adaptative uniform
1.8
1.7
1.6
q low
1.5
1.4
1.3
1.2
1.1
1
5
6
7
8
9
10
11
12
log(DoF)
Fig. 7.11. qup (top) and qlow (bottom) with respect to n for Test 4 for rectangular meshes.
976
SERGE NICAISE
0.112 adaptative uniform 0.111
0.11
q up
0.109
0.108
0.107
0.106
0.105
0.104
5
6
7
8
9
10
11
12
log(DoF)
6 adaptative uniform 5.5
5
q up
4.5
4
3.5
3
2.5
5
6
7
8
9
10
11
12
log(DoF)
Fig. 7.12. qup (top) and qlow (bottom) with respect to n for Test 4 for triangular meshes.
A POSTERIORI ANALYSIS OF FINITE VOLUME METHODS
977
REFERENCES [1] A. Agouzal and F. Oudin, A posteriori error estimator for finite volume methods, Appl. Math. Comput., 110 (2000), pp. 239–250. [2] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., Teubner, Stuttgart, 1999. [3] T. Apel and G. Lube, Anisotropic mesh refinement in stabilized Galerkin methods, Numer. Math., 74 (1996), pp. 261–282. [4] A. Bergam and Z. Mghazli, Estimateurs a posteriori d’un sch´ ema de volumes finis pour un probl` eme non lin´ eaire, C. R. Acad. Sci. Paris S´er. I Math. 331 (2000), pp. 475–478. ¨rth, A posteriori estimators for the finite volume [5] A. Bergam, Z. Mghazli, and R. Verfu discretization of an elliptic problem, Numer. Math., 95 (2003), pp. 599–624. [6] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, North–Holland, Amsterdam, 1978. `re, J.-P. Villa, and P. Villedieu, Convergence rate of a finite volume scheme for [7] Y. Coudie a two dimensional convection-diffusion problem, M2AN Math. Model. Numer. Anal., 33 (1999), pp. 493–516. `re and P. Villedieu, Convergence rate of a finite volume scheme for the lin[8] Y. Coudie ear convection-diffusion equation on locally refined meshes, M2AN Math. Model. Numer. Anal., 34 (2000), pp. 1123–1149. [9] J. Droniou, Error estimates for the convergence of a finite volume discretization of convectiondiffusion equations, J. Numer. Math., 11 (2003), pp. 1–32. ¨t, and R. Herbin, Finite volume methods, in Handbook of Numerical [10] R. Eymard, T. Galloue Analysis, vol. 7, P. Ciarlet and J.-L. Lions, eds., North–Holland, Amsterdam, 2000, pp. 723– 1020. [11] R. Hangleiter and G. Lube, Stabilized Galerkin methods and layer-adapted grids for elliptic problems, Comput. Methods Appl. Mech. Engrg., 166 (1998), pp. 165–182. [12] R. Herbin and M. Ohlberger, A posteriori error estimate for finite volume approximations of convection diffusion problems, in Finite Volumes for Complex Applications III, R. Herbin and D. Kr¨ oner, eds., Lab. Anal. Topol. Probab. CNRS, Paris, 2002, pp. 739–746. [13] N. Jullian, An error indicator for cell-centered finite volumes for linear convection-diffusion problems, in Finite Volumes for Complex Applications III, R. Herbin and D. Kr¨ oner, eds., Lab. Anal. Topol. Probab. CNRS, Paris, 2002, pp. 763–770. ¨ ner and M. Ohlberger, A posteriori error estimates for upwind finite volume schemes [14] D. Kro for nonlinear conservation laws in multi dimensions, Math. Comp., 69 (1999), pp. 25–39. [15] G. Kunert, A Posteriori Error Estimation for Anisotropic Tetrahedral and Triangular Finite Element Meshes, Logos Verlag, Berlin, 1999; Ph.D. thesis version available online at http://archiv.tu-chemnitz.de/pub/1999/0012/index.html. [16] G. Kunert, An a posteriori residual error estimator for the finite element method on anisotropic tetrahedral meshes, Numer. Math., 86 (2000), pp. 471–490. [17] G. Kunert, A local problem error estimator for anisotropic tetrahedral finite element meshes, SIAM J. Numer. Anal., 39 (2001), pp. 668–689. [18] G. Kunert, Robust a posteriori error estimation for a singularly perturbed reaction-diffusion equation on anisotropic tetrahedral meshes, Adv. Comp. Math., 15 (2001), pp. 237–259. [19] G. Kunert, A posteriori error estimation for convection dominated problems on anisotropic meshes, Math. Methods Appl. Sci., 26 (2003), pp. 589–617. [20] G. Kunert, Z. Mghazli, and S. Nicaise, A Posteriori Error Estimation for a Finite Volume Discretization on Anisotropic Meshes, preprint SFB393/03–016, TU Chemnitz, Chemnitz, Germany, 2003. Available online at http://archiv.tu-chemnitz.de/pub/2003/0005/ index.html. [21] R. Lazarov and S. Tomov, Adaptive finite volume element method for convection-diffusionreaction problems in 3-D, in Scientific Computing and Application, Y. W. P. Minev and Y. Lin, eds., Nova Sci. Publ., Huntington, NY, 2001, pp. 91–106. [22] R. Lazarov and S. Tomov, A posteriori error estimates for finite volume approximations of convection-diffusion-reaction equations, Comput. Geosci., 6 (2002), pp. 483–503. [23] J. Mackenzie, T. Sonar, and G. Warnecke, A posteriori error estimates for the cell-vertex finite volume method, in Adaptive Methods: Algorithms, Theory and Applications, W. Hackbusch and G. Wittum, eds., Vieweg, Braunschweig, 1994, pp. 221–235. [24] L. Morley, The triangular equilibrium element in the solution of plate bending problems, Aero. Quarterly, 19 (1968), pp. 149–169.
978
SERGE NICAISE
¨li, A posteriori and a priori error analysis of finite volume methods, [25] K. W. Morton and E. Su in The Mathematics of Finite Elements and Applications, J. R. Whiteman, ed., John Wiley and Sons, New York, 1994, pp. 267–288. [26] S. Nicaise, A posteriori error estimations of some cell-centered finite volume methods, SIAM J. Numer. Anal., 43 (2005), pp. 1481–1503. [27] S. Nicaise, A posteriori residual error estimation of a cell-centered finite volume method, C. R. Acad. Sci. Paris, S´er., 338 (2004), pp. 419–424. [28] T. K. Nilssen, X.-C. Tai, and R. Winther, A robust nonconforming H 2 -element, Math. Comp., 70 (2000), pp. 489–505. [29] M. Ohlberger, A posteriori error estimate for finite volume approximations to singularly perturbed nonlinear convection-diffusion equations, Numer. Math., 87 (2001), pp. 737–761. [30] M. Ohlberger, A posteriori error estimates for vertex centered finite volume approximations of convection-diffusion-reaction equations, M2AN Math. Model. Anal. Numer., 35 (2001), pp. 355–387. [31] S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Ser. Comput. Methods Mech. Thermal Sci., McGraw–Hill, New York, 1980. [32] H.-G. Roos, M. Stynes, and L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations. Convection-Diffusion and Flow Problems, Springer, Berlin, 1996. ¨li, A dual graph-norm refinement indicator for finite volume approxima[33] T. Sonar and E. Su tions of the Euler equations, Numer. Math., 78 (1998), pp. 619–658. ¨rth, A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement [34] R. Verfu Techniques, Wiley–Teubner, Chichester, Stuttgart, 1996. ¨rth, A posteriori error estimators for convection-diffusion equations, Numer. Math., [35] R. Verfu 80 (1998), pp. 641–663. ¨rth, Robust a posteriori error estimators for the singularly perturbed Helmholtz [36] R. Verfu equation, Numer. Math., 78 (1998), pp. 479–493.