RESIDUAL-BASED A POSTERIORI ERROR ESTIMATION FOR MULTIPOINT FLUX MIXED FINITE ELEMENT METHODSThis work was supported in part by the Education Science Foundation of Chongqing (KJ120420), National Natural Science Foundation of China (11171239), The Project-sponsored by Scientific Research Foundation for the Returned Overseas Chinese Scholars and Open Fund of Key Laboratory of Mountain Hazards and Earth Surface Processes, CAS.

# Residual-based a posteriori error estimation for multipoint flux mixed finite element methods

Abstract. A novel residual-type a posteriori error analysis technique is developed for multipoint flux mixed finite element methods for flow in porous media in two or three space dimensions. The derived a posteriori error estimator for the velocity and pressure error in norm consists of discretization and quadrature indicators, and is shown to be reliable and efficient. The main tools of analysis are a locally postprocessed approximation to the pressure solution of an auxiliary problem and a quadrature error estimate. Numerical experiments are presented to illustrate the competitive behavior of the estimator.

Key words. multipoint flux mixed finite element method, postprocessed approximation, a posteriori error estimate

AMS subject classifications. 65N06, 65N12, 65N15, 65N30, 76S05,

## 1 Introduction

Let be a bounded polygonal () or polyhedral () domain with a Lipschitz continuous boundary . We consider the following first-order system of diffusion-type partial differential equations:

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩u=−K∇pin  Ω,∇⋅u=fin Ω,p=gon ΓD,u⋅n=0on ΓN. (1)

Here are partitions of the boundary corresponding to the Dirichlet and Neumann conditions, respectively, with , and , is the outward unit normal vector on , and is a symmetric and uniformly positive definite tensor with

 k0ξTξ≤ξTK(x)ξ≤k1ξTξ,  ∀ x∈Ω, ∀ ξ∈Rd (2)

for . This system has been widely used in physics to model diffusion processes such as heat or mass transfer and flow in porous media. In flow in porous media, denotes the pressure, is the Darcy velocity, and represents the permeability divided by the viscosity.

The main goal of this paper is to derive residual-based a posteriori error estimation for multipoint flux mixed finite element (MFMFE) methods for the model (1). The MFMFE approach was developed for single phase flow in porous media in [30, 39, 40]. It is motivated by the multipoint flux approximation (MPFA) approach [2, 1, 26, 32, 33], which is a control volume method developed by the oil industry as a reliable discretization for single-phase Darcy flow. One main advantage of this method lies in that, by introducing sub-edge (or sub-face) fluxes, it provides a local explicit flux with respect to the flow pressure, and allows for local flux elimination around grid vertices and reduction to a cell-centered pressure scheme. The MFMFE method is based on the lowest order Brezzi-Douglas-Marini (BDM1) [17] or Brezzi-Douglas-Duran-Fortin (BDDF1) [16] finite element space. By using special quadrature rules, local velocity elimination is also attained which leads to a symmetric and positive definite cell-centered system for the pressure on quadrilateral, simplicial and hexahedral meshes. In [41], a coupling discretization of MFMFE method and continuous Galerkin finite element method was applied to the poroelasticity system that describes fluid flow in deformable porous media.

It is well-known that adaptive algorithms for the numerical solution of partial differential equations are nowadays standard tools in science and engineering. A posteriori error estimation, as an essential ingredient of adaptivity, provides adaptive mesh refinement strategy and quantitative estimates of the numerical solution obtained. For second-order elliptic problems, the theory of a posteriori error estimation has reaches a degree of maturity for finite element of conforming, nonconforming and mixed types (see [3, 4, 5, 6, 12, 13, 7, 14, 15, 20, 21, 11, 22, 24, 31, 34, 37] and the references therein). To the authors’ knowledge, no a posteriori estimation for the MFMFE method has been proposed in the literature so far.

In this paper, we develop a novel technique to derive residual-based a posteriori error estimation for the MFMFE method for the porous media model in two or three-dimensional case. Since the MFMFE method employs a special quadrature rule, its a posteriori error estimator should include a term to control the error of quadrature. This is different from the standard analytical technique based on the discrete -inner product. Moreover, we can not directly utilize the analytical technique developed by Carstensen in [21] for nonconforming finite elements to estimate

 infβ∈H1(Ω)||∇β−K−1uh||,

because the BDM1 finite element for the velocity approximation, , does not have the same continuity of mean of trace across the interior sides as the nonconforming finite elements do. To overcome this difficulty, we shall construct a locally postprocessed approximation to the pressure solution, obtained by the MFMFE scheme, of a special auxiliary problem, and use a derived estimate of quadrature error. We note that the idea of postprocessing in this contribute follows from the works [34, 38].

The rest of this paper is organized as follows. In section 2, we introduce some notations and the continuous problem. Section 3 shows the MFMFE method. Section 4 includes main results. Sections 5-6 are respectively devoted to the a posteriori error estimation and the analysis of efficiency. Finally, we illustrate the performance of the obtained estimation in section 7 by numerical experiments.

## 2 Notations and continuous problem

Let be a shape regular triangulation of in the sense of [23] which satisfies the angle condition, namely there exists a constant such that for all

 C−10hdT≤|T|≤C0hdT,

where . Let be a piecewise constant function with .

We denote by the set of element sides (or faces) in , by the set of sides (or faces) of element , by and respectively the sets of the interior and Dirichlet boundary sides (or faces) of all elements in , by the union of all elements in sharing side (or face) , and by the set of nodes in .

For a domain , let be the inner product on , and the dual pair between and . Let be the usual Sobolev space consisting of functions defined on with all derivatives of order up to belonging to , with norm . When , and , especially for . We omit the subscript if . For a tensor-valued function , let for any norm . Introduce

 H(div;A):={v∈L2(A)d:∇⋅v∈L2(A)},

and define the ”broken Sobolev space”

 H1(∪Th):={φ∈L2(Ω):φ|T∈H1(T),∀T∈Th}.

We denote by the jump of over an interior side with diameter , shared by the two neighboring (closed) elements . Especially, if .

Since we consider two and three-dimensional cases () simultaneously, the Curl of a function with if and if is defined by

 Curlψ:=(−∂2ψ,∂1ψ)  if d=2  and  Curlψ:=∇×ψ  if d=3,

where denotes the usual vector product of two vectors in . Given a unit normal vector along the side , we define the tangential component of a vector with respect to by

 γtE(v):={v⋅(−n2,n1)if  d=2,v×nEif  d=3.

Throughout the paper, denotes the local version of differential operator defined by for all . We also use the notation to represent where is a generic, positive constant independent of the mesh size of . Moreover, abbreviates .

Denote

 V:={v∈H(div;Ω) :v⋅n=0  on ΓN},  W:=L2(Ω),

then the weak formulation of the model (1) is as follows: Find , such that

 (K−1u,v)=(p,∇⋅v)−ΓD,   ∀ v∈ V, (3)
 (∇⋅u,w)=(f,w),   ∀ w∈ W. (4)

It is well-known that this problem admits a unique solution [18].

## 3 Multipoint flux mixed finite element method

We follow the notations and definitions employed in [39, 30] to describe the MFMFE method. Let be the reference element which is a unit triangle in two-dimensional case or unit tetrahedron in three-dimensional case, and be the set of polynomials of degree . The lowest order mixed finite element spaces on are defined as

 ^V(^T)=P1(^T)d,   ^W(^T)=P0(^T).

Since for any and any edge (or face) of , the degrees of freedom for can be chosen to be the values of at any two points on each edge of if is the unit triangle, or any three points on each face of if is the unit tetrahedron [18, 17]. In the MFMFE method, these points are chosen to be the vertices of for the requirement of accuracy and certain orthogonality for the trapezoidal quadrature rules. Such a choice allows for local velocity elimination and leads to a cell-centered stencil for the pressure [39, 30].

The lowest order spaces on are given by

 Vh:={v∈V :   v|T=1JTDFT^v∘F−1T,  ^v∈^V(^T)   ∀ T∈Th},Wh:={w∈W :   w|T=^w∘F−1T,  ^w∈^W(^T)   ∀ T∈Th},

where is the inverse mapping of the bijection , is the Jacobian matrix with respect to on the element with . Note that the vector transformation is is known as the Piola transformation.

For , it holds

 ∫TK−1q⋅vdx=∫^T^K−11JTDFT^q⋅1JTDFT^vJTd^x=∫^T1JT(DFT)T^K−1DFT^q⋅^vd^x=∫^TK−1^q⋅^vd^x

with The quadrature formula on an element is then defined as [39, 30]

 (K−1q,v)Q,T:=(K−1^q,^v)^Q,^T:=|^T|ss∑i=1K−1(^ri)^q(^ri)⋅^v(^ri), (5)

where () are the corresponding vertices of with for the unit triangle and for the unit tetrahedron.

Define the global quadrature formula as

 (K−1q,v)Q=∑T∈Th(K−1q,v)Q,T, (6)

then the MFMFE method is formulated as follows: Find and such that

 (K−1uh,vh)Q=(ph,∇⋅vh)−ΓD,   ∀ vh∈ Vh, (7)
 (∇⋅uh,wh)=(f,wh),   ∀ wh∈ Wh. (8)

The existence and uniqueness of the solution to the scheme (7)-(8) follow from [39, 30]. As shown in [39, 30], the algebraic system that arises from (7)-(8) is of the form

 (  A    BT−B   0)( U P)=( G F), (9)

where with and , and , are respectively the bases of and . The matrix is block-diagonal with symmetric and positive definite blocks, and the local elimination of leads to a system for with a symmetric and positive definite matrix . For the details, we refer to [39, 30].

## 4 Main results

Let be the discretization indicator defined by

 η2h:=||h(f−∇⋅uh)||2+∑T∈Th∑E∈εThEJ2tE, (10)

where

 J2tE:=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩||[γtE(K−1uh)]||2Eif  E∈ε0h∩∂T,||γtE(K−1uh)−∂g/∂s||2E+h2E||∂2g∂s2||2Eif  E∈∂T∩εD,0if  E∈∂T∩ΓN, (11)

and and denote respectively the first and second order tangential derivatives of function along side . Introduce the quadrature indicator

 η2Q:=∑T∈Thh2T||uh||21,T. (12)

We note this indicator is owing to the use of the special quadrature formula (5) in the MFMFE method.

We now state in Theorems 4-4 a posteriori error estimates for the errors of velocity and pressure in norm, respectively.

{theorem}

Let be the weak solution of the continuous problem (3)-(4), and be the solution of the MFMFE method (7)-(8). Assume . Then it holds

 ||K−1/2(u−uh)||≲(η2h+η2Q)1/2. (13)
{theorem}

Assume . Under the assumptions of Theorem 4, it holds

 ||Qhp−ph||≲hmax(ηh+ηQ)+||h(f−∇⋅uh)||, (14)
 ||p−ph||≲hmax(ηh+ηQ)+||hK−1uh||+||h(f−∇⋅uh)||. (15)

Here , and denotes the projection operator onto .

###### Remark 4.1

We note that the two terms and in in the estimator are of high order with respect to the lowest order scheme, which are usually omitted in computation. In fact, from (8) it follows , and turns out to be an oscillation term of high order.

###### Remark 4.2

The above estimates (13)-(15) also apply to the original mixed finite element discretization where the special quadrature rule (5) is not used in the scheme (7)-(8). In this case, the estimator is not involved, and then in the estimates (13)-(15). In this sense, our work can be regarded as a generalization of Carstensen’s [20] to the three-dimensional case. We note that our estimator is a bit different from that in [20] due to no occurrence of the term ( denotes the piecewise operator acting on element by element in ). Here we also consider more general boundary conditions.

We finally state in Theorem 4 the efficiency of the a posteriori error estimators. Note that the efficiency of a reliable a posteriori error estimator means that its converse estimate holds up to high order terms and different multiplicative constants. For the sake of simplicity, we assume that is a matrix of piecewise polynomial functions.

{theorem}

Under the assumptions of Theorems 4-4, it holds

 ηh+ηQ+h−1max||hK−1uh||≲||K−1/2(u−uh)||+||h−1(p−ph)||+h.o.t..

where denotes some high-order term depending on given data.

## 5 A posteriori error analysis

This section is devoted to the proofs of Theorems 4-4.

Introduce the global quadrature error and the element quadrature error as follows:

 σ(K−1uh,vh)|T=σT(K−1uh,vh):=(K−1uh,vh)T−(K−1uh,vh)Q,T,,  for all  T∈Th. (16)

Let denote the lowest order element space on .

We state two estimates on the quadrature error derived in [39, 30] as follows. If for all element , then it holds

 |σ(K−1qh,vh)|≲∑T∈ThhT||qh||1,T||vh||T (17)

for all , . Moreover, if for all element , then it holds

 |σ(K−1qh,vh)|≲∑T∈Thh2T||qh||1,T||vh||1,T (18)

for all .

Denote respectively by and the standard projection operators from onto and for some (cf. [20, 39]). It holds the following estimates:

 ||h−1(q−Π0q)||≲||q||H1(∪Th)   for all  q∈(H1(∪Th))d∩H(div;Ω), (19)
 ||Π0v||1,T≲||v||1,T,||Πv||1,T≲||v||1,T   for all  v∈(H1(T))d,  ∀T∈Th. (20)

Note that bound (19) can be found in [20], and bounds (20) are the direct results of Lemma 3.1 in [39].

To derive a reliable a posteriori error estimate for the velocity error, we need to introduce an auxiliary problem as following:

 ⎧⎪⎨⎪⎩∇⋅(K∇ϑ)=∇⋅uhin  Ω,ϑ=−gon  ΓD,K∇ϑ⋅n=0on  ΓN. (21)

Since is a symmetric and uniformly positive definite tensor, by the Lax-Milgram theorem there exists a unique solution to this problem, provided that . As is divergence-free, a decomposition of two or three-dimensional vector fields (see Theorem 3.4 and Remark 3.10 in [28]) implies that there exists a stream function such that

 K∇ϑ−uh=Curl ψ.

Since and vanish on , we easily know on .

Introduce , then and it holds

 u−uh=−K∇p−K∇ϑ+Curl ψ=K∇z+Curl ψ. (22)

 ||K−1/2(u−uh)||2=∫ΩK−1(u−uh)⋅(u−uh)=∫Ω(∇z+K−1Curl ψ)⋅(K∇z+Curl ψ)=∫ΩK∇z⋅∇z+2∫Ω∇z⋅Curl ψ+∫ΩK−1Curl ψ⋅Curl ψ. (23)

Using integration by parts and noticing on and on , we have

 ∫Ω∇z⋅Curl ψ=−∫Ω∇⋅(Curl ψ)z+∫ΓD∪ΓNCurl ψ⋅nz=0. (24)

Notice that , on and on . The relation (24) and integration by parts yield

 ∫ΩK∇z⋅∇z=∫Ω(u−uh)⋅∇z=−∫Ω∇⋅(u−uh)z. (25)

Let denote the projection of onto . From (4) and (8) it follows

 (∇⋅(u−uh),Qhz)=0. (26)

In view of , the above two relations, (25) and (26), imply

 ∫ΩK∇z⋅∇z=−∫Ω∇⋅(u−uh)(z−Qhz)=∑T∈Th∫T(−f+∇⋅uh)(z−Qhz)≲∑T∈ThhT||f−∇⋅uh||T||∇z||T≲||h(f−∇⋅uh)|| ||K1/2∇z||,

which results in

 ||K1/2∇z||≲||h(f−∇⋅uh)||. (27)

By (22) and (24) we have

 ||K−1/2(u−uh)||2=||K1/2∇z||2+||K−1/2Curl ψ||2. (28)

Recalling for all , in light of (22) we have, for any ,

 ∫ΩK−1Curl ψ⋅Curl ψ=∫Ω(K−1(u−uh)−∇z)⋅Curl ψ =∫ΩK−1(u−uh−K∇v)⋅Curl ψ =∫ΩK−1(u−K∇v−K∇β)⋅Curl ψ+∫ΩK−1(K∇β−uh)⋅Curl ψ ≤(||K−1(u−K∇v−K∇β)||+||∇β−K−1uh||)||Curl ψ||,

which implies

 ||K−1/2Curl ψ||≲infv∈H1D(Ω)||K−1(u−K∇v−K∇β)||+infβ∈H1(Ω)||∇β−K−1uh||. (29)

Finally, from (27)-(29) it follows

 ||K−1/2(u−uh)||≲{infv∈H1D(Ω)||K−1(u−K∇v−K∇β)||+infβ∈H1(Ω)||∇β−K−1uh||+||h(f−∇⋅uh)||}. (30)

In what follows, we shall follow the routines of [21] to estimate the first and second terms on the right-hand side of (30). To this end, we assume that and for all and denote by the nodal piecewise linear interpolation of on which satisfies for all . Let be the nodal basis of the lowest order finite element space associated to , i.e., for all , for , and . Denote by . We then introduce a subspace of , , as follows (see [21]):

 ~S:==⎧⎪⎨⎪⎩∑z∈Nφzvz:∀ z∈N,vz∈C(ωz), vz|ωz is a piecewisepolynomial, and vz=−gh,D on ΓD∩ωz.⎫⎪⎬⎪⎭
{lemma}

For , it holds

 infv∈H1D(Ω)||K−1(u−K∇v−K∇β)||≲{∑E⊂ΓDh3E||∂2g/∂s2||2E}1/2. (31)
{proof}

The definition of shows on . Noticing , we have

 infv∈H1D(Ω)||K−1(u−K∇v−K∇β)||=infw∈H1(Ω),w|ΓD=g−gh,D||∇w||.

The desired result (31) immediately follows from an estimate in the proof of Lemma 3.4 in [21].

On the other hand, it holds

 infβ∈H1(Ω)||∇β−K−1uh||≤infvh∈~S||∇vh−K−1uh||. (32)

It is sophisticated to give a computational upper bound for the right-hand side term of (32) with the help of and given data. To this end, let denote the piecewise mean value of on , i.e. for all Then is symmetric and has the following ellipticity:

 k−11ξTξ≤ξT¯¯¯¯¯¯¯¯¯¯K−1ξ≤k−10ξTξ   for all x∈Ω,  ξ∈Rd.

Recall that is the lowest order element space on . and is the piecewise constant space.Introduce the following auxiliary problem: Find such that

 (¯¯¯¯¯¯¯¯¯¯K−1~uh,vh)=(~ph,∇⋅vh)−ΓD,  ∀ vh∈V0h, (33)
 (∇⋅~uh,wh)=(f,wh),   ∀ wh∈Wh. (34)

It is well-known that this problem admits a unique solution (see [18]). {lemma} Let