A posteriori analysis of a positive streamwise invariant discretization of a convection-diffusion equation

A posteriori analysis of a positive streamwise invariant discretization of a convection-diffusion equation by Christine Bernardi1, Tom´as Chac´on Rebollo2, and Marco Restelli2 Abstract: We consider the finite element discretization of a convection-diffusion equation, where the convection term is handled via a fluctuation splitting algorithm. We prove a posteriori error estimates which allow us to perform mesh adaptivity in order to optimize the discretization of these equations. Numerical results confirm the interest of such an approach. R´esum´e: Nous consid´erons une discr´etisation par ´el´ements finis d’une ´equation de convection–diffusion, o`u un algorithme de d´ecentrage est utilis´e pour traiter le terme de convection. Nous prouvons des estimations d’erreur a posteriori qui permettent d’adapter le maillage pour optimiser la discr´etisation de ces ´equations. Des r´esultats num´eriques confirment l’int´erˆet d’une telle approche. 1Laboratoire Jacques-Louis Lions, C.N.R.S. & Universit´e Pierre et Marie Curie - Paris 6, boˆıte 187, 4 place Jussieu, 75252 Paris Cedex 05, France. b[email protected] 2Departamento de Ecuaciones Diferenciales y An´alisis Numerico, Universidad de Sevilla, Tarfia s/n, 41012 Sevilla, Spain. [email protected], [email protected] 1. Introduction. We are interested in the discretization of the following convection-diffusion equation, set in a bounded connected domain Ω in R2or R3, −α∆p+ (u·∇)p=fin Ω,(1.1) and provided with appropriate boundary conditions, in both cases where the diffusion term or the convection term is the leading one. The discretization that we propose relies on Lagrange finite element methods of low order. However, as standard for convection dominated problems, a further upwind scheme must be inserted into the discrete problem in order to have good convergence properties. After the pioneering papers [15] and [10] where the method of characteristics is introduced, a large amount of work has been performed concerning the construction of new upwind schemes and their analysis. We quote among others the discontinuous Galerkin methods [7], several kinds of finite volume discretizations [9][19], the Nscheme [14] and the positive streamwise invariant (PSI) algorithm introduced in [8]. We use here the abstract framework proposed in [5] for handling the more general fluctuation splitting methods, also called residual distribution methods [1], which includes the PSI scheme and other secondorder non-linear residual-distribution schemes (see also [5] for a more general bibliography on these upwind schemes). On the opposite, it seems that very few works deal with the a posteriori analysis of such methods (see [2] for the method of characteristics and [19] for finite volume methods). In this paper, we perform the a posteriori analysis of the PSI method, in the framework proposed in [5]. We observe that the techniques of proof are very similar to their analogue in the case of a simple Galerkin method and we establish optimal upper bounds of the error as a function of residual type error indicators (see [18] for the description of these indicators). Conversely, these indicators are bounded from above by the local error. So our estimates are fully optimal and the error indicators constitute a very efficient tool for mesh adaptivity. We next describe our adaptivity strategy. The numerical experiments that we present are aimed to prove that mesh adaptivity is a very efficient tool for computing some solutions in realistic cases, for instance when they present boundary layers. They confirm the interest of using the PSI algorithm in strongly convection dominated regimes. An outline of the paper is as follows: •In Section 2, we present the abstract framework for the fluctuation splitting methods. Next, we recall from [5, §6 & 7] the main results concerning the well-posedness of the discrete problem and the convergence of its solution. •Section 3 is devoted to the a posteriori analysis of the discrete problem. •Numerical experiments are presented in Section 4. 1 2. The continuous and discrete problems. Let Ω be a connected bounded open set in Rd,d= 2 or 3, with a Lipschitz–continuous boundary ∂Ω. The generic point in Ω is denoted by x, while the generic point on ∂Ω is denoted by τ. As usual, nstands for the unit outward normal vector to Ω on ∂Ω. With standard notation for the Sobolev spaces, we introduce a vector function uin W1,r(Ω)3, r > d 2, and we make the further assumption that div u= 0 in Ω.(2.1) We are thus in a position to introduce a partition of ∂Ω without overlap into Γ−and Γ+, with Γ−=τ∈∂Ω; (u·n)(τ)<0,Γ+=τ∈∂Ω; (u·n)(τ)≥0, and we assume from now on that Γ−has a positive measure in ∂Ω. We are interested in the following equation          −α∆p+ (u·∇)p=fin Ω, p= 0 on Γ−, α ∂np= 0 on Γ+, (2.2) where the diffusion coefficient αis a positive constant. The unknown is the scalar function pand the datum fis assumed to be regular enough. Note that we have chosen to work with zero boundary conditions only for simplicity. Let H1 ∗(Ω) denote the space H1 ∗(Ω) = nq∈H1(Ω); q= 0 on Γ−o.(2.3) It is readily checked that, when the boundaries of Γ−and Γ+are smooth enough, system (2.2) admits the equivalent variational formulation Find pin H1 ∗(Ω) such that ∀q∈H1 ∗(Ω), a(p, q) = hf, qi,(2.4) where the bilinear form a(·,·) is defined by a(p, q) = αZΩ (grad p)(x)·(grad q)(x) dx+ZΩ(u·∇)p(x)q(x) dx,(2.5) while h·,·i stands for the duality pairing between H1 ∗(Ω) and its dual space H1 ∗(Ω)0. Owing to the assumption on u, we deduce from standard Sobolev imbeddings that the form a(·,·) is continuous on H1(Ω) ×H1(Ω). Moreover, we have for any function qin H1 ∗(Ω) a(q, q) = α|q|2 H1(Ω) +1 2ZΓ+ (u·n)(τ)q2(τ)dτ≥α|q|2 H1(Ω). Thus, using the Poincar´e–Friedrichs inequality (indeed, Γ−has a positive measure), we obtain the following ellipticity property: There exists a positive constant csuch that ∀q∈H1 ∗(Ω), a(q, q)≥c α kqk2 H1(Ω).(2.6) 2 All this together with the Lax–Milgram lemma yields that, for any datum fin H1 ∗(Ω)0, problem (2.4) has a unique solution. We assume from now on that Ω is a polygon (d= 2) or a polyhedron (d= 3). Let (Th)hbe a regular family of triangulations of Ω (by triangles or tetrahedra), in the sense that: •For each h, Ω is the union of all elements of Th; •The intersection of two different elements of Th, if not empty, is a vertex or a whole edge or a whole face of both of them; •The ratio of the diameter hKof any element Kof Thto the diameter of its inscribed circle or sphere is smaller than a constant σindependent of h. As usual, hstands for the maximum of the diameters hK,K∈ Th. We make the further and non restrictive assumption that both Γ−and Γ+are the union of whole edges (d= 2) or whole faces (d= 3) of elements of Th. In what follows, c,c0,. . . are generic constants that may vary from line to line but are always independent of h. In view of the finite element discretization, we introduce the space Xh=nqh∈H1(Ω); ∀K∈ Th, qh|K∈ P1(K)o, where P1(K) is the space of restrictions to Kof affine functions with dvariables. We also need its subspace X∗ h=Xh∩H1 ∗(Ω). In view of the fluctuation splitting method, we introduce another finite-dimensional space Y∗ hmade of piecewise constant functions. We assume that the dimension of Y∗ his the same as the dimension of X∗ hand moreover we introduce a one-to-one linear mapping Πhfrom X∗ honto Y∗ h. Next, assuming that the datum fbelongs to L2(Ω), we consider the discrete problem Find phin X∗ hsuch that ∀qh∈X∗ h, ah(ph, qh) = ZΩ f(x)(Πhqh)(x) dx,(2.7) where the bilinear form ah(·,·) is defined by ah(ph, qh) = αZΩ (grad ph)(x)·(grad qh)(x) dx+ZΩ(u·∇)ph(x) (Πhqh)(x) dx.(2.8) To go further, we now describe the links between the spaces X∗ hand Y∗ hin the framework of the fluctuation splitting method. Let bj, 1 ≤j≤J, denote all the vertices of the elements of Th. With these vertices, we associate the standard Lagrange functions ϕj: For 1≤j≤J, the function ϕjbelong to Xhand satisfies ϕj(bj) = 1 and ϕj(bi)=0,1≤i≤J, i 6=j. (2.9) The functions ϕj, 1 ≤j≤J, form a basis of Xhand without restriction we assume that the ϕj, 1 ≤j≤J∗, form a basis of X∗ h(with J∗< J). The fluctuation splitting method relies on the choice of linearly independent functions λj, 1 ≤j≤J∗, which are piecewise constant and have a small support in a neighbourhood 3 of bj(“small” means contained in the union of triangles or tetrehedra in Ththat contain bj). Then, the operator Πhis defined by Πhqh= J∗ X j=1 qh(bj)λj,(2.10) and it is readily checked that it is one-to-one from X∗ honto the space Y∗ hspanned by the functions λj, 1 ≤j≤J∗. Of course, since our aim is to construct upwind schemes, the operator Πh, hence the functions λj, depend on the convection vector u. They can also depend on the fluxes RKu(x)·(∇ph)(x) dx, and in this case we use the notation Πh(ph) and λj(ph) when needed (note that in this case the discrete problem is nonlinear). For each Kin Th, we denote by bjKthe vertices of K, where jKis the global index corresponding to the local index j, 1 ≤j≤d+1, and by λK jthe restriction of the function λjKto K. From now on, we assume that the λK jare zero when bjis not a vertex of K. Example 2.1. The Nscheme [14] For each Kin Th, we introduce the normal vector to the edge opposite to bjK nK j=dmeas(K)grad ϕjK. Denoting by uKthe mean value of uon K, we set aK j=1 duK·nK j,(aK j)+= max{aK j,0},(aK j)−= min{aK j,0}, bK= d+1 X j=1 (aK j)−. Finally, for any function ρhin X∗ h, we define the functions λK j(ρh) by λK j(ρh) = 1 RKu(x)·(∇ρh)(x) dx d+1 X i=1 (aK j)+bK(aK i)−ρh(bK j)−ρh(bK i).(2.11) Despite their complex form, the functions λK jare easy to compute, for instance when uis constant on K. Note however that the local flux RKu(x)·(∇ρh)(x) dxcan vanish and, in this case, the functions λK jare not defined. Example 2.2. The PSI method With the previous notation, the functions λK jare now defined by λK j(ρh) = 1−µK j RKu(x)·(∇ρh)(x) dx d+1 X i=1 (aK j)+bK(aK i)−ρh(bK j)−ρh(bK i),(2.12) where the constants µK jare chosen such that 0≤µK j≤1,1≤j≤d+ 1,and d+1 X j=1 λK j(ρh)=1, and also such that all λK j,K∈ Th, are bounded independently of ρh. It can be noted that the only difference between the Nscheme and the PSI method relies on the choice of the 4 constants µK jwhich are added to stabilize the discretization. Examples of constants µK j satisfying the previous properties are given in [1] for instance. To conclude, we recall from [5, §6 & 7] the main results concerning the well-posedness of problem (2.7) and the convergence of the discrete solution phto p. We are led to make a few assumptions for that. Assumption 2.3. The following property holds for each Kin Th λK j≥0,1≤j≤d+ 1,and d+1 X j=1 λK j= 1.(2.13) We now introduce the matrices linked to our discrete problem: Dstands for the matrix with coefficients RΩ(grad ϕi)(x)·(grad ϕj)(x) dx, 1 ≤i, j ≤J∗, and, for any ρhin X∗ h,C(ρh) stands for the matrix with coefficients RΩ(u·∇)ϕi(x)λi(ρh)(x) dx, 1≤i, j ≤J∗. We need a further definition. Definition 2.4. A matrix Awith coefficients Aij , 1 ≤i, j ≤N, is a M-matrix if Aii >0,1≤i≤N, and Aij ≤0,1≤i, j ≤N, i 6=j, Aii ≥X j6=i|Aij|. Assumption 2.5. For any ρhin X∗ hand any real number α > 0, the matrix α D +C(ρh) is a M-matrix. Assumption 2.6. The mapping: ρh7→ C(ρh) is continuous from X∗ hinto the space of square real matrices with dimension J∗×J∗. Indeed, the following results are established in [5, Thms 6.1 & 7.1]. Proposition 2.7. If Assumptions 2.3,2.5and 2.6are satisfied, for any datum fin L2(Ω), problem (2.7) has a solution phin X∗ hfor hsmall enough. Moreover, this solution satisfies kphkH1(Ω) ≤ckfkL2(Ω).(2.14) Proposition 2.8. If Assumptions 2.3,2.5and 2.6are satisfied, any sequence of solutions (ph)hof problem (2.7) converges to the solution pof problem (2.4) strongly in H1(Ω) when htends to zero. We refer to [5, Thm 7.3] for an explicit a priori estimate of the error between pand phunder the previous assumptions: The scheme is first-order acurate in the H1(Ω)-norm. We also observe that the Nscheme and the PSI method both satisfy Assumption 2.3. On the other hand, they satisfy Assumption 2.5 only when the angles between edges (d= 2) or the dihedral angles between faces (d= 3) of all elements of Thare ≤π 2. And the PSI method always satisfies Assumption 2.6, but the N-scheme does not. To conclude, let us remark that a relevant feature of the schemes verifying Assumptions 2.3, 2.5 and 2.6 is the verification of the discrete maximum principle, although this only occurs with piecewise affine discretizations. 5 3. A posteriori analysis of the discrete problem. We need some further notation. With each Kin Th, we associate (i) the set E0 Kof all edges (d= 2) or faces (d= 3) of Kwhich are not contained in ∂Ω; (ii) the sets E+ Kof all edges (d= 2) or faces (d= 3) of Kwhich are contained in Γ+; (iii) the domain ωKequal to the union of all elements of Ththat share at least an edge (d= 2) or a face (d= 3) with K. For each ein E0 K, we denote by [·]ethe jump through e(the introduction of a vector normal to eis necessary to make precise the sign of this jump, however we do not need it in what follows). For each ein E0 Kand in E+ K,hestands for the length (d= 2) or diameter (d= 3) of e. We introduce a piecewise constant approximation fhof the datum f: For each Kin Th, we set fh|K=1 meas(K)ZK f(x) dx. We also need an approximation uhof u. Several choices are possible for that: •Assuming that uis continuous on Ω, we can define uhas the Lagrange interpolate of u with values in Xh; •A more consistent choice consists in taking uhin the Raviart–Thomas finite element space [16] Zh=nvh∈H(div,Ω); ∀K∈ Th,vh|K∈ PRT (K)o, where H(div,Ω) is the domain of the divergence operator in L2(Ω)dand PRT (K) stands for the space of restrictions to Kof polynomials of the form µ+νx,µ∈Rd,ν∈R. Indeed, using the corresponding interpolation operator (see [16]) leads to a function uh which is exactly divergence-free. We are now in a position to define the error indicators: For each Kin Th, ηK=hKkfh−(uh·∇)phkL2(K)+X e∈E0 K αh 1 2 ek[∂nph]ekL2(e)+X e∈E+ K αh 1 2 ek∂nphkL2(e).(3.1) It can be noted that all the functions which appear in the norms are polynomials with low degree, so that these indicators are very easy to compute once the discrete solution is known. We successively prove an upper and a lower bound of the error as a function of these indicators. We set: εK=hKk(u−uh)·∇phkL2(K).(3.2) Proposition 3.1. If Assumption 2.3is satisfied, the following a posteriori estimate holds for the error between the solution pof problem (2.4) and any solution phof problem (2.7) αkp−phkH1(Ω) ≤cX K∈Thη2 K+ε2 K+h2 Kkf−fhk2 L2(K)1 2.(3.3) Proof: Using (2.4) yields the residual equation, for any qun H1 ∗(Ω), a(p−ph, q) = ZΩ f(x)q(x) dx−a(ph, q). 6 Then, it follows from (2.7) that, for any qhin X∗ h, a(p−ph, q) = ZΩ f(x)(q−Πhqh)(x) dx −αZΩ (grad ph)(x)·grad (q−qh)(x) dx−ZΩ(u·∇)ph(x) (q−Πhqh)(x) dx. By integrating by parts on each K, we observe that, since phis affine on K, αZK (grad ph)(x)·grad (q−qh)(x) dx=X e∈E0 K∪E+ K αZe (∂nph)(τ)(q−qh)(τ) dτ. By also inserting the approximations fhand uhin the previous equation, we obtain, with the appropriate signs for the jumps [∂nph]e, a(p−ph, q) =X K∈Th ZK (f−fh)(x)(q−Πhqh)(x) dx−ZK((u−uh)·∇)ph(x) (q−Πhqh)(x) dx +ZKfh−(uh·∇)ph(x) (q−Πhqh)(x) dx −1 2X e∈E0 K αZe [∂nph]e(τ)(q−qh)(τ) dτ−X e∈E+ K αZe (∂nph)(τ)(q−qh)(τ) dτ. Next, we take qhequal to the image of qby a Cl´ement type regularization operator, see [6]. We recall from [3, Chap. IX, Thms 3.7, 3.8 & Cor. 3.9] that this operator satisfies, for all Kin Thand ein E0 Kor E+ K kq−qhkL2(K)+hK|qh|H1(K)≤c hK|q|H1(∆K),kq−qhkL2(e)≤c h 1 2 e|q|H1(∆K), where ∆Kis the union of elements of Ththat intersect K. On the other hand, it is proven in [5, Lemma 5.2] that, owing to Assumption 2.3, for any qhin X∗ hand Kin Th, kqh−ΠhqhkL2(K)≤c hK|qh|H1(K). Combining all this with Cauchy–Schwarz inequalities yields a(p−ph, q)≤cX K∈Th (ηK+εK+hKkf−fhkL2(K))|q|H1(∆K), whence (note that each Kin This only included in a finite number of ∆K0, where “finite” means bounded as a function of the regularity parameter σof the family of triangulations) a(p−ph, q)≤cX K∈Th (η2 K+ε2 K+h2 Kkf−fhk2 L2(K))1 2|q|H1(Ω). To conclude, we take qequal to p−phand use the ellipticity property (2.6) of the form a(·,·). 7 Figure 7: The reference solution −1− 0 . 5 0 0 . 5 1 0 0.2 0.4 0.6 0.8 1 −1− 0 . 5 0 0 . 5 1 0 0.2 0.4 0.6 0.8 1 −1− 0 . 5 0 0 . 5 1 0 0.2 0.4 0.6 0.8 1 −1− 0 . 5 0 0 . 5 1 0 0.2 0.4 0.6 0.8 1 −1− 0 . 5 0 0 . 5 1 0 0.2 0.4 0.6 0.8 1 −1− 0 . 5 0 0 . 5 1 0 0.2 0.4 0.6 0.8 1 Figure 8: The adapted meshes We now compute the triangulations Thby iterating the adaptation step with the constant ρequal to 0.9. Figure 8 presents from left to right and top to bottom the initial mesh next the resulting adapted meshes, after 5, 10, 15, 20 and 25 iterations of the adaptation step. It can be noted that the final meshes are refined exactly at the location 14 of the sharp front, which proves the efficiency of our strategy. Finally, Figure 9 presents the Hilbertian sum of the indicators ηh(circles), and the errors in the norms of H1(Ω) (stars) and L2(Ω) (crosses) as a function of h(defined as previously as equal to J−1 2) in bilogarithmic scales, all these quantities being divided by their mean value. The expected order for uniform refinement is 1 in the H1(Ω) norm, so it seems that clustering the points around the sharp gradient region yields a higher convergence rate. 0.1 1 h 0.01 0.1 Figure 9: The Hilbertian sum of the error indicators as a function of h. For the last test case, the domain is [0,1]2and we still use the coefficient α= 10−4 and data f= 0 and r+= 0. We now consider the veloicity u= (−∂yψ, ∂xψ)T, where is the stream function ψis given by ψ(x, y) = −y−1 3x2(1 −2y) + 1 5πxcos(kπx2) sin(pπy),(4.6) with k= 16 and p= 4. 0 0 .2 0 .4 0 . 6 0 . 8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 10: The stream lines of the velocity u 15 Such a velocity field is chosen in order to mimic a turbulent flow, and the corresponding stream lines are plotted in Figure 10. The Dirichlet boundary conditions on Γ−are p−(x, y) = y5(1 −y)5,if x= 0, 0,otherwise. (4.7) Due to the choice of the velocity field, Γ−and Γ+are not connected (the inflow boundary Γ−is indicated by a grey thick line on Figure 10). Hence, to simplify the treatment of the boundary conditions, the Dirichlet boundary conditions are enforced weakly employing the method by Nitsche, for which we refer to [13] and [4], which allows the use of computational meshes such that the intersection Γ−∩Γ+does not necessarily coincide with a subset of the mesh nodes. In fact, the initial mesh considered in the mesh adaptation process uses 16 uniform subdivisions for each side of the computational domain, so that the small scale inflow/outflow pattern is not resolved. The computed solution after 20 iterations of the adaption step is depicted in Figure 11, and the corresponding adapted grid is shown in Figure 12 (left part), together with a zoom on the right upper part of the domain (right part). Figure 11: The computed solution 16 0 0 .2 0 .4 0 . 6 0 . 8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 . 8 0 . 85 0 . 9 0 . 95 1 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1 Figure 12: The adapted mesh Even in these rather hard situations, the PSI method seems very efficient to solve the advection-dominated problems and is fully compatible with mesh adaptivity. 17 References [1] R. Abgrall, M. Mezine — Construction of second-order accurate monotone and stable residual distribution schemes for steady problems, J. Comput. Physics 195 (2004), 474–507. [2] C. Bernardi, F. Hecht, R. Verf¨urth — A posteriori error analysis of the method of characteristics, to appear in Math. Models and Methods in Applied Sciences. [3] C. Bernardi, Y. Maday, F. Rapetti — Discr´etisations variationnelles de probl`emes aux limites elliptiques, Collection “Math´ematiques et Applications” 45, Springer-Verlag (2004). [4] E. Burman — A unified analysis for conforming and nonconforming stabilized finite element methods using interior penalty, SIAM J. Numer. Anal. 43 (2005), 2012–2033. [5] T. Chac´on Rebollo, M. G´omez M´armol, G. Narbona Reina — Numerical analysis of the PSI solution of advection–diffusion problems through a Petrov–Galerkin formulation, Math. Models and Methods in Applied Sciences 17 (2004), 1905–1936. [6] P. Cl´ement — Approximation by finite element functions using local regularization, R.A.I.R.O. Anal. Num´er. 9R2 (1975), 77–84. [7] B. Cockburn — An introduction to the discontinuous Galerkin method for advection-dominated problems, Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, Lecture Notes in Math. 1697, Springer (1998), pp. 151–268. [8] H. Deconinck, R. Struijs, G. Bourgeois. P.L. Roe — Compact advection schemes on unstructued meshes, Comput. Fluid Dynamics, VKI Lecture Series 1993-04 (1993). [9] B. Despr´es — Lax theorem and finite volume schemes, Math. Comput. 73 (2004), 1203–1234. [10] J. Douglas, Jr., T.F. Russell — Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures, SIAM J. Numer. Anal. 19 (1982), 871–885. [11] P.J. Frey, P.-L. George — Maillages, applications aux ´el´ements finis, Herm`es (1999). [12] F. Hecht, O. Pironneau — FreeFem++, see www.freefem.org. [13] J. Nitsche — ¨ Uber ein Variationsprinzip zur L¨osung von Dirichlet-Problemen bei Verwendung von Teilr¨aumen, die keinen Randbedingungen unterworfen sind, Abh. Math. Sem. Univ. Hamburg 36 (1971), 9–15. [14] B. Perthame — Convergence of N-schemes for linear advection equations, Trends in Applications of Mathematics to Mechanics, Pitman Monogr. Surveys Pure Appl. Math. 77 (1985), pp. 323–333. [15] O. Pironneau — On the transport-diffusion algorithm and its applications to the Navier–Stokes equations, Numer. Math. 38 (1982), 309–332. [16] P.-A. Raviart, J.-M. Thomas — A mixed finite element method for second order elliptic problems, Mathematical Aspects of Finite Element Methods, Lecture Notes in Mathematics 606, Springer (1977), pp. 292–315. [17] R.M. Smith, A. G. Hutton — The numerical treatment of advection: A performance comparison of current methods, Numer. Heat Transfer, Part A: Applications 5(1982), 439–461. [18] R. Verf¨urth — A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques, Wiley & Teubner (1996). [19] M. Vohral´ık — Residual flux-based a posteriori error estimates for finite volume and related locally conservative methods, Numer. Math. 111 (2008), 121–158. [20] M. Vohral´ık — Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients, submitted for publication (2008). 18