A SYMMETRIC NODAL CONSERVATIVE FINITE ELEMENT ... - LNCC

A SYMMETRIC NODAL CONSERVATIVE FINITE ELEMENT METHOD FOR THE DARCY EQUATION ´ ERIC ´ GABRIEL R. BARRENECHEA, LEOPOLDO P. FRANCA 1 , AND FRED VALENTIN

2

Abstract. This work introduces and analyzes novel stable Petrov-Galerkin Enriched Methods (PGEM) for the Darcy problem based on the simplest but unstable continuous P1 /P0 pair. Stability is recovered inside a Petrov-Galerkin framework where element-wise dependent residual functions, named multi-scale functions, enrich both velocity and pressure trial spaces. Unlike the velocity test space that is augmented with bubble-like functions, multiscale functions correct edge residuals as well. The multi-scale functions turn out to be the well-known lowest order Raviart-Thomas basis functions for the velocity and discontinuous quadratics polynomial functions for the pressure. The enrichment strategy suggests the way to recover the local mass conservation property for nodal-based interpolation spaces. We prove that the method and its symmetric version are well-posed and achieve optimal error estimates in natural norms. Numerical validations confirm claimed theoretical results.

1. Introduction The Darcy equation arising in a porous media field belongs to the family of mixed problems [13] for which numerical methods are limited by the choice of pair of approximation spaces. From classical stable elements as the Raviart-Thomas family (RTk ) [27], Brezzi-DouglasMarini elements (BDMk ) [12], high order stable elements given in [25, 5, 6] to more recent stabilized or least square finite element methods [26, 23, 14, 10, 11] the range of possibilities to tackle the Darcy equation has increased over the past years. Methods for this problem should combine stability and accuracy while preserving physical properties inherited from the continuous problem. Properties that are only fulfilled by few of them. To the best of our knowledge, symmetric stable nodal based finite element methods for the Darcy equation preserving mass locally remain an open problem (see [11] for a recent discussion). For example, least-squares finite element methods (cf. [10]) lead to a symmetric positive-definite system, but, in their original nodal version they are not locally mass conservative, and in [11] nodal unknowns for the velocity are forbidden. Furthermore, in Key words and phrases. Darcy model; Enriched space; simplest element; Petrov-Galerkin approach. 1 This author is supported by NSF/USA No. 0610039. 2

This author is supported by CNPq /Brazil Grant No. 306255/2008-1 and 304051/2006-3, FAPERJ/Brazil

Grant No. E-26/100.519/2007, and by NSF/USA No. 0610039. 1

2

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

both cases, the lowest order piecewise constant space for the pressure is not allowed. This possibility is considered in [25], but the degrees of freedom for the velocity are not nodal. Some of the previous requirements are satisfied by the so-called Petrov-Galerkin Enriched Methods (PGEM). PGEM have been developed in [19, 20, 21] and further analyzed in [18, 4]. The method is constructed by enriching polynomial functions with two types of enhancement: we add bubble functions to the test function and we add a special function to the trial function. The latter depends on the residual of the polynomial part over edges, and thus, it is no longer a bubble-like function. This gives us a Petrov-Galerkin framework. To get stabilized method forms of PGEM we use static condensation, thanks to the use of bubbles as test functions [3, 8]. Number and type of degrees of freedom stay unchanged whereas basis functions incorporate unsolved sub-scales modifying their form yet preserving the polynomial basis function support. Interested readers can find a review on the subject in [2]. When applied to the Darcy equation, the Petrov-Galerkin approach leads to different finite element methods [7]. One of them is obtained by searching the velocity solution into a subspace of the Raviart-Thomas space built with the Raviart-Thomas interpolation over the linear continuous trial functions. The space for the pressure stays untouched. The underlying PGEM appears to be stable for the simplest pair of interpolation spaces P1 /P0 while preserving the mass-conservative feature, a desirable property for porous media practitioners. Performance of PGEM over several numerical tests given in [7] attests its stability and accuracy while keeping loss of local mass negligible. Based on the previous considerations, the current work introduces a variant of the strategy proposed in [7] and leads to a final method which is symmetric, locally mass conservative and whose degrees of freedom are piecewise constants for the pressure and nodal values for the velocity. Indeed, we keep the trial space for the velocity and pressure as in [7], but the test space is built differently: it is first mapped using the Raviart-Thomas interpolation operator and then enhanced with bubble functions, an approach which allows the static condensation procedure. This new perspective opens the door to two new finite element methods, one of them fulfilling all the requirements of symmetry, nodal degrees of freedom for the velocity and locally mass conservative. Both methods prove to be well-posed and achieve optimal error estimates in natural norms. Since the starting point of our approach is a Petrov-Galerkin method, the terminology PGEM is still used in this work, even if the final methods differ from the ones presented in [7]. Finally, the approach suggests a general way of rendering some finite element methods locally mass conservative. We end this introduction by summarizing the plan and main results of this paper:

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

3

• In Section 2 we introduce the new finite element methods, namely, the non-symmetric

(13) and its symmetric counterpart (16). We then derive the methods in a PGEM framework in Section 2.1, and in Section 2.2 we explicit the local problems solved by our enrichment basis functions (cf. (35) and (37)). In particular, we recover the local basis of the RT0 space as the solution of (35)-(36). Finally, in Section 2.3 we prove

that the discrete enhanced solution is locally mass conservative; • Section 3 is devoted to the error analysis. We analyze in detail method (16) (the analysis of (13) is treated in Theorem 9, Section 3.3). Well posedeness and consis-

tency error are proved first (cf. Lemmas 3 and 4) followed by convergence results (see Theorems 7 and 8 in Section 3.2). Furthermore, we use the characterization of the RT0 interpolation operator as the solution of (35)-(36) to obtain an alternative proof for the classical RT0 error estimate (see Corollary 6); • The numerical tests are in Section 4 where two analytical solutions confirm theoretical results;

• Conclusions and future perspectives are drawn in Section 5;

• Finally, we relax the assumption on the source term g (initially assumed piecewise

constant) to propose in Appendix A an error estimate for a smooth datum g (cf. Theorem 10).

1.1. Some notations. This section introduces definitions and notations used throughout. In what follows, Ω denotes an open bounded domain in R2 with polygonal boundary ∂Ω, and x = (x1 , x2 ) is a typical point in Ω. As usual, L2 (Ω) is the space of square integrable functions over Ω, L20 (Ω) represents functions belonging to L2 (Ω) with zero average in Ω, and H div (Ω) is composed by functions that belong to L2 (Ω)2 with divergence in L2 (Ω). The space H0div (Ω) stands for the space of functions belonging to H div (Ω) which have normal component vanishing on ∂Ω. From now on we denote by {Th } a family of regular triangulations of Ω built up using

triangles K with boundary ∂K composed by edges F . The set of internal edges of the triangulation Th is denoted by Eh . The characteristic length of K and F are denoted by hK

and hF , respectively, and h := max{hK : K ∈ Th } > 0, and due to the mesh regularity there exists a positive constant C such that hF ≤ hK ≤ C hF , for all F ⊆ ∂K. Also, for each ′

F = K ∩ K ∈ Eh we choose, once and for all, an unit normal vector n which coincides with

the unit outward normal vector when F ⊆ ∂Ω. The standard outward normal vector at the edge F with respect to the element K is denoted by nK F . Moreover, for a function q, one

4

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

denotes JqK its jump, defined by (see Figure 1): JqK(x) := lim+ q(x + δn) − lim− q(x + δn) ,

(1)

δ→0

δ→0

and JqK = 0 if F ⊆ ∂Ω.

Next, we denote by H0div (K) the space whose functions belong to H0div (Ω) with support in

K and vanishing normal component on ∂K, and L20 (K) the space of functions which belong to L20 (Ω) with support and zero mean in K. Then, we can define the corresponding global spaces H0div (Th ) := ⊕

X

K∈Th

H0div (K) and L20 (Th ) := ⊕

X

L20 (K).

K∈Th

Finally, (· , · )D stands for the inner product in L2 (D) (or in L2 (D)2 , when necessary), and

k· ks,D (|· |s,D ) the norm (seminorm) in H s (D) (or H s (D)2 , if necessary), and k· kdiv,D the norm in H div (D).

F n

K’

K

Figure 1. The normal vector.

1.2. Preliminaries. In this work we consider the following Darcy problem: Find (u, p) such that σ u + ∇p = f ,

(2)

∇· u = g

in Ω,

u· n = 0 on ∂Ω, where σ =

µ κ

∈ R+ is assumed constant in Ω, with µ and κ denoting the viscosity and

permeability, respectively. Here, u is the so-called Darcy velocity, p is the pressure, f and

g are given source terms. We suppose f piecewise constant since it is usually related to the gravity force. Moreover, we assume that the given data have enough regularity and the usual compatibility condition Z

Ω

g = 0,

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

5

holds. Remark. When we consider (2) with a prescribed flux b on ∂Ω such that Z Z g= b, Ω

∂Ω

we can recover the homogeneous case since there exists a function w b belonging to H div (Ω) such that w b · n = b on ∂Ω (cf. [22]), and thus we replace the right hand side f by f − σw b

and g by g − ∇· wb . 

Remark. In the more general case, σ can always be approximated by performing projections onto the piecewise constant space. This possibility has been considered in [7]. On the other hand, despite the fact that the methods are presented in the two dimensional case, their extension to the three dimensional framework is straightforward.  The standard symmetric mixed variational formulation associated with (2) reads: Find (u, p) ∈ H0div (Ω) × L20 (Ω) such that (3)

As ((u, p), (v, q)) = Fs (v, q) ∀(v, q) ∈ H0div (Ω) × L20 (Ω),

where As ((u, p), (v, q)) := (σ u, v)Ω − (p, ∇· v)Ω − (q, ∇· u)Ω ,

Fs (v, q) := (f , v)Ω − (g, q)Ω .

The well-posedness of (3) follows from the classical Babuska-Brezzi theory for variational problems with constraints (see [13] for more details). Remark. An equivalent and still well-posed non-symmetric version of (3) arises from adding the weak form of the second equation to the first one in (2). The bilinear form and the linear form are now denoted by A(., .) and F(.), respectively, and are given by (4) A((u, p), (v, q)) := (σ u, v)Ω − (p, ∇· v)Ω + (q, ∇· u)Ω ,

F(v, q) := (f , v)Ω + (g, q)Ω .

Next, the classical discrete mixed formulation of this problem is: Find (uh , ph ) ∈ Vh × Qh

such that (5)

As ((uh , ph ), (v h , qh )) = Fs (v h , qh ) ∀(v h , qh ) ∈ Vh × Qh ,

where Vh and Qh are finite-dimensional approximations of H0div (Ω) and L20 (Ω), respectively. It is well known that the pair of interpolation spaces for pressure and velocity must satisfy the discrete Babuska-Brezzi (or inf-sup) condition [13] in order to lead to a stable discrete version of problem (5). For the Darcy model containing a zero order term, this restriction

6

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

has been proved to be unnecessary. In fact, in [24] standard continuous Lagrangian finite element spaces have been proved to be stable and convergent. The lowest order Raviart-Thomas space is one of the simplest examples of a stable mass conservative element, and is composed by the velocity space VRT0 := {v ∈ H0div (Ω) : v |K ∈ RT0 (K) ∀K ∈ Th }, where the local space RT0 (K) is defined by (6)

RT0 (K) := P0 (K)2 + x P0 (K),

and (7)

Qh := {q0 ∈ L20 (Ω) : q0 |K ∈ P0 (K) ∀ K ∈ Th } .

Hence, only the normal component of the velocity is continuous and the inf-sup condition is satisfied since ∇· VRT0 = Qh . Associated with the space RT0 (K) there exists a natural local

interpolation operator πK : [H 1 (K)]2 → RT0 (K), defined by (cf. [13, 17]) Z Z (8) πK (v)· n = v· n , F

F

for all F ∈ ∂K, or, equivalently (9)

πK (v) :=

X

F ⊆∂K

R

F

v· n ϕF , hF

where ϕF is the Raviart-Thomas’ basis function given by (10)

ϕF (x) = ±

hF (x − xF ), 2|K|

and xF denotes the node opposite to the edge F . Hence, a global interpolation operator noted π : [H 1 (Ω)]2 → VRT0 follows by defining π(v) |K = πK (v) in each K ∈ Th . Remark. The sign before the Raviart-Thomas basis function ϕF depends on whether the normal vector n on F ⊆ ∂K points inwards or outwards K.  A lifting operator from L1 (F ) to VRT0 will be needed in the sequel, it is denoted by ℓ and P is such that ℓ(q) := F ∈Eh ℓF (q) where R αF F q (11) ϕF , ℓF (q) = σ where the coefficient αF is a given positive constant which is independent of hF and σ, but can vary with F . We finally denote, for K ∈ Th , ℓK (q) := ℓ(q)|K .

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

7

2. The finite element methods We begin by introducing the standard finite element space Vh := [Vh ]2 ∩ H0div (Ω) for the

velocity variable, where

Vh := {v ∈ C 0 (Ω) : v|K ∈ P1 (K), ∀K ∈ Th },

(12)

whereas the pressure is discretized using the space Qh defined in (7). We start by presenting the methods that will be analyzed in this work. First, the Petrov-Galerkin Enriched Method reads: Find (u1 , p0 ) ∈ Vh × Qh such that (13)

B((u1 , p0 ), (v 1 , q0 )) = Fs (π(v 1 ), q0 ) ,

for all (v 1 , q0 ) ∈ Vh × Qh , where (14) B((u1 , p0 ), (v 1 , q0 )) := As ((π(u1 ), p0 ), (π(v 1 ), q0 )) + (ℓ(Jp0 K), σ π(v 1 ))Ω −

X

τF (Jp0 K, Jq0 K)F ,

F ∈Eh

and π and ℓ are the operators defined through (8) and (11), respectively, and the coefficient τF stands for τF :=

(15)

αF hF . σ

In Section 3 this problem is proved to be well-posed for an appropriate choice of αF . ˆ 1 , pˆ0 ) ∈ Alternatively, a symmetric related formulation can also be derived and reads: Find (u

Vh × Qh such that

ˆ 1 , pˆ0 ), (v1 , q0 )) = Fs (π(v 1 ), q0 ) , Bs ((u

(16)

for all (v 1 , q0 ) ∈ Vh × Qh , where (17)

Bs ((u1 , p0 ), (v 1 , q0 )) := As ((π(u1 ), p0 ), (π(v 1 ), q0 )) −

X

τF (Jp0 K, Jq0 K)F .

F ∈Eh

Remark. The latter method is based on the error analysis (see §3) which points out that we

can remove from (13) the term (ℓ(Jp0 K), σ π(v 1 ))Ω without introducing a loss of accuracy. Moreover, we recover the symmetric form of the so-called reduced PGEM method presented

in [7], by replacing the term (π(u1 ), π(v 1 ))Ω by (u1 , v 1 )Ω in (13). This reduced method turns out to be optimally convergent [14], and in §2.3 we show how to render it locally mass conservative. 

8

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

2.1. Derivation of the methods. The starting point towards our final method is the following Petrov-Galerkin method for (2): Find uh := u1 + ue ∈ Vh + H0div (Ω) and ph := p0 + pe ∈ Qh ⊕ L20 (Th ) such that (18)

As ((uh , ph ), (vh , qh )) = Fs (v h , qh ),

for all v h := π(v 1 ) + v b ∈ π(Vh ) ⊕ H0div (Th ) and for all qh := q0 + qe ∈ Qh ⊕ L20 (Th ). Here

π(Vh ) stands for the subspace of VRT0 built as the image of space Vh through the operator π. This scheme is equivalent to the following system: for all (v 1 , q0 ) ∈ Vh × Qh and for all

(v b , qe ) ∈ H0div (Th ) × L20 (Th ) (19)

As ((uh , ph ), (π(v1 ), q0 )) = Fs (π(v 1 ), q0 ),

(20)

As ((uh , ph ), (v b , qe )) = Fs (v b , qe ).

From now on, and just in order to derive the method, we will assume that g is a piecewise constant function (even if the method is analyzed and implemented for more general functions g). With this assumption in mind, starting from (20) and proceeding as in [7], the following strong problem is obtained for (ue , pe ): (21) (22)

σ ue + ∇pe = f − σu1 , ∇· ue = Jp0 K + (ΠK (∇ · u1 ) − ∇ · u1 ) in K, Z
σ ue · n = αF Jp0 K + (f − σu1 ) · n − ΠF ((f − σu1 ) · n) , F

on each F ⊆ ∂K ∩ Eh , and ue · n = 0 on F ⊆ ∂Ω. Here, ΠF and ΠK stand for the L2 R R 1 projection operators over the constant space, i.e, ΠF (v) = h1F F v and ΠK (v) = |K| v. K

Finally, the constant Jp0 K is chosen in order to make (21)-(22) compatible, and is given by 3

1 X αFi hFi Jp0 K := |K| i=1 σ

Z

Fi

Jp0 K n· nK Fi .

Remark. For higher order velocity interpolation the vanishing right hand side term ΠK (∇ · u1 ) − ∇ · u1 left in the local problem (21) needs to be taken into account. Moreover, written

in this form, it will help us to bound the consistency errors (see equation (62) below). On the other hand, since we have assumed that f is a constant (or piecewise constant) function,

the divergence equation in (21) may be rewritten as follows (23)

∇· ue = Jp0 K +

1 (∇ · (f − σu1 ) − ΠK (∇ · (f − σu1 ))) , σ

which is the form that we will consider from now on. 

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

9

Remark. The boundary condition (22) assures the continuity of the normal component of the velocity on each edge, thus keeping our approach conforming. This fact may be kept even for discontinuous coefficients, as it has been done in [7], where the mean value of σ has been included in the boundary condition.  p u Now, let MK := (MuK , MpK ) : H 1 (K)2 → H0div (K) × L20 (K) and DK := (DK , DK ) :

L2 (∂K) → H div (K)×L20 (K), defined as follows: (v e , ηe ) := (MuK (v), MpK (v)) is the solution

of

(24)

σv e + ∇ηe = v ,

σ ∇ · v e = ∇ · v − ΠK (∇ · v) in K,

σ v e · n = v · n − ΠF (v · n) on each F ⊆ ∂K, p u and (we , ξe ) := (DK (q), DK (q)) solves

(25)

Z 3 1 X αFi hFi σ w e + ∇ξe = 0 , ∇ · w e = q n· nK Fi |K| i=1 σ Fi Z σ w e · n = αF q on each F ⊆ ∂K ∩ Ω.

in K,

F

Then, using these operators and (23), we can characterize the solution (ue , pe ) = (uM e + M D uD e , pe + pe ) of (21)-(22) as follows

(26)

M (uM e , pe ) = MK (f − σu1 )

(27)

D (uD e , pe ) = DK (Jp0 K)

∀K ∈ Th ,

∀K ∈ Th .

Next, we turn back to equation (19). First, since pe ∈ L20 (K) and ∇· π(v1 ) |K ∈ R we

obtain (28)

(pe , ∇· π(v1 ))K = 0 for all K ∈ Th .

Therefore, the problem (19) becomes: Find (u1 , p0 ) ∈ Vh × Qh such that (29)

As ((u1 + ue , p0 ), (π(v1 ), q0 )) = Fs (π(v 1 ), q0 ) ∀(v 1 , q0 ) ∈ Vh × Qh ,

where ue is characterized with respect to u1 and p0 by (26)-(27). It is also convenient to rewrite the problem above in an equivalent form integrating it by parts in each K ∈ Th (30)

X

K∈Th

(q0 , ∇· ue )K =

X

F ∈Eh

τF (Jp0 K, Jq0 K)F .

10

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

Remark. The term related to f in (26) vanishes. Indeed, since f is constant in K then f · n − ΠF (f · n) = 0 on each edge F , and there also exists a polynomial function qe belonging to L20 (K) such that

∇qe = f

in each K ∈ Th ,

which leads to MuK f = 0. Therefore, no enriching contribution comes from (26) but for the one related to u1 . 

Finally, based on the previous remark, and replacing (30) and (26)-(27) in (29), we arrive at the following final form of PGEM: Find (u1 , p0 ) ∈ Vh × Qh such that X u (DK (Jp0 K), σ π(v 1 ))K As ((u1 − σMuK (u1 ), p0 ), (π(v 1 ), q0 )) + K∈Th

−

X

τF (Jp0 K, Jq0 K)F = Fs (π(v 1 ), q0 ),

F ∈Eh

for all (v 1 , q0 ) ∈ Vh × Qh , which is precisely the method (13) since, as it will be shown in

terms of the basis function in the next section, the following holds (31)

u πK (u1 ) ≡ (I − σMuK )(u1 ) and ℓK (Jp0 K) ≡ DK (Jp0 K).

Next, the symmetric method (16) follows by neglecting the non diagonal term (ℓ(Jˆ p0 K), σ π(v 1 ))Ω as this term does not undermine convergence estimates (see Section 3). Remark. Following analogous steps and just replacing the forms As (., .) and Fs (.) by A(., .) and F(.), and switching the velocity test space π(Vh ) to Vh in (18) we arrive at the following ˆ 1 , pˆ0 ) ∈ Vh × Qh such that method: Find (u (32)

ˆ 1 ), pˆ0 ), (v1 , q0 )) + (ℓ(Jˆ A((π(u p0 K), σ v 1 )Ω +

X

τF (Jˆ p0 K, Jq0 K)F = F(v 1 , q0 ) ,

F ∈Eh

for all (v 1 , q0 ) ∈ Vh × Qh , which is exactly the PGEM proposed in [7].  2.2. The local problems. This section is devoted to show the relationship between the local problems (26)-(27) and the Raviart-Thomas interpolation operator π and the lifting P P operator ℓ. First, we decompose u1 = 2k=1 3i=1 uki ψ ki , where uki are the nodal values of u1

and ψ ki , that denotes the (vector-valued) hat function. Then, we look for solutions of (26) and (27) in the form

(33)

uD e |K =

3 α X Fj j=1

R

Jp0 K Fj σ

ϕj

and pD e |K =

3 α X Fj j=1

R

Fj

σ

Jp0 K

ηj ,

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

11

and (34)

uM e |K

= −σ

2 X 3 X

uki

ϕki

and

k=1 i=1

pM e |K

= −σ

2 X 3 X

uki η ki .

k=1 i=1

Here, indexes j and i are representing edges and local nodal numeration, respectively. Next, by replacing (33) and (34) in the local problems (27) and (26) respectively, and factoring out the coefficients it follows that the multi-scale basis functions (ϕj , ηj ) and (ϕki , ηki ) must satisfy the following well-posed local Darcy problems (35) (36)

hF σ ϕj + ∇ηj = 0 , ∇· ϕj = j n· nK in K, Fj |K| ( 1 if j = i, on each Fi ⊆ ∂K, ϕj · n = 0 otherwise

and σ ϕki + ∇η ki = ψ ki ,

(37) (38)

∇· ϕki = 0 in K,

σ ϕki · n = ψ ki · n − ΠF (ψ ki · n) on each F ⊆ ∂K.

The enrichment functions emanating from the problem (35)-(36) are nothing but the well known basis functions of the space RT0 (K), i.e, the lowest order Raviart-Thomas approximation of H div (K) defined by (6), and they are given by (cf. (10)) (39)

ϕj (x) = ±

hFj (x − xFj ) for j = 1, 2, 3 . 2|K|

Consequently,  σ hFj  |x|2 − x· xFj + Cj for j = 1, 2, 3, 2 |K| 2 R where the constant Cj is set up so K ηj = 0. The solution of (37)-(38) can also be analyti(40)

ηj (x) = ∓

cally computed, providing

(41)

σϕki

=

ψ ki

−

3 X j=1

ΠF (ψ ki · n) ϕj = ψ ki − πK (ψ ki ) ,

where πK has been defined in (9). A similar local problem has been used in [15] to obtain multiscale basis functions for the Darcy problem with oscillating coefficients (see also [1] for the extension to porous media with stochastic coefficients).

12

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

2.3. The local mass conservation feature: a general strategy. Computing numerical solutions through PGEM (13) and the symmetric method (16) does not assure local mass ˆ 1 ) is considered. Regarding conservative velocity field if only linear part of solution u1 (or u the first method, the required feature is achieved by locally updating u1 with uD e given by (27). In fact, since discontinuous pressure interpolations are used it emerges from (13) that (see [7] for a related idea) Z

K

∇· (u1 +

uD e )

=

Z

g, K

for all K ∈ Th .

Remark. Summing up, we see that in order to obtain a stable pair of interpolation spaces with

discrete velocity field locally mass-conservative, it is fundamental to enrich the linear part of the discrete velocity u1 with an element of the Raviart-Thomas’ space VRT0 , namely, the multi-scale function uD e computed from (33). We stress the fact that the computation of (33) follows directly from the discrete solution, without the need of any extra local computation. Moreover, the exact velocity u is approximated in each K ∈ Th by D uh = u1 + uM e + ue

=

2 X 3 X k=1 i=1

=

uki

ψ ki

−σ

2 X 3 X k=1 i=1

uki

ϕki

Z 3 X αFl Jp0 K ϕl + σ Fl l=1

3   X ΠFl (u1 · n) + τFl ΠFl (Jp0 K) ϕl l=1

= πK (u1 ) + ℓK (Jp0 K),

so, as expected, the continuity of the normal velocity component through the internal edges is assured, but not the tangential one.  Now, it turns out that such local mass recovering is not only restricted to methods arising exactly from the enhancing approach. For example, stabilized methods based on pressure jumps as the one presented in [14], or the symmetric method (16) are elegible to recover the local mass conservation feature adding ℓK (Jp0 K) to the computed velocity field, where p0 is the constant pressure solution. We illustrate this fact for the symmetric formulation (16). Choosing v 1 = 0 in (16) and q0 = 1 in K and −|K|/|K ′ | in K ′ (where K ∩ K ′ = F ∈ Eh ),

we obtain after integration by parts and the definition of ℓK that Z Z |K| ˆ 1 + ℓK (Jˆ ˆ 1 + ℓK (Jˆ ∇· (u p0 K)) − ∇· (u p0 K)) = |K ′ | K ′ K Z Z Z Z X |K| |K| ˆ1 + ˆ1 − ∇·u g, τF (Jˆ p0 K, Jq0 K)F = g− ∇·u |K ′ | K ′ |K ′ | K ′ K K ′ F ⊆∂K∪∂K

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

and then, following closely the arguments given in [7], we obtain that the value

13

R

K

ˆ1 − ∇ · (u

ℓK (Jˆ p0 K)) − g vanishes in each K, and hence it must vanish on each element. It is important

to emphasize once more that ℓ(Jˆ p0 K) does not perturb too much the solution in the sense that the order of error estimates is still preserved (see Theorems 8 and 9 for details). 3. Error analysis

In the sequel C denotes a generic positive constant, independent of h or σ, with values that may vary in each occurrence. Before performing an error analysis of (13) and (16), we need to consider interpolation inequalities to approximate variables. 3.1. Interpolation, stability and consistency results. We start by presenting the Cl´ement interpolation operator (cf. [16, 22, 17]) Ch : H 1 (Ω) → Vh (with the obvious extension to vector-valued functions), satisfying, for all K ∈ Th and all F ∈ Eh , (42) (43)

kv − Ch (v)km,K ≤ C ht−m |v|t,ωK K t− 21

kv − Ch (v)k0,F ≤ C hF

|v|t,ωF

∀v ∈ H t (ωK ) , ∀v ∈ H t (ωF ) ,

for t = 1, 2, m = 0, 1, where ωK = {K ′ ∈ Th : K ∩ K ′ 6= ∅} and ωF = {K ∈ Th : K ∩ F 6= ∅}. Now, in order to take into account the approximation of the pressure and the consistency

error, we consider the L2 (Ω) projection onto Qh which is denoted by Πh : L2 (Ω) → Qh . This projection satisfies (cf. [17]) (44)

kq − Πh (q)km,Ω ≤ C h1−m |q|1,Ω

∀ q ∈ H 1 (Ω) ,

for m = 0, 1. Moreover, using the result above and the following local trace inequality: given K ∈ Th , F ⊆ ∂K, there exists C such that for all v ∈ H 1 (K)   1 2 2 2 (45) kvk0,K + hK |v|1,K , kvk0,F ≤ C hK we obtain (46)

hX

F ∈Eh

hF kJq −

Πh (q)Kk20,F

i1/2

≤ C h |q|1,Ω .

Moreover, we will systematically use the Raviart-Thomas interpolation operator π defined through (9) as π(v) =

X

F ∈Eh

ΠF (v · n|F ) ϕF ,

satisfying (see [17] or Corollary 6 for an alternative proof) (47)

kv − π(v)k0,Ω ≤ C h |v|1,Ω ,

14

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

for all v ∈ H 1 (Ω)2 and (48)

k∇ · v − ∇ · π(v)k0,Ω ≤ C h |∇ · v|1,Ω ,

for all v ∈ H 1 (Ω)2 such that ∇ · v ∈ H 1 (Ω).

Now, we define the following mesh dependent norm

(49)

k(v, q)kh :=

 X

σ

K∈Th

kvk20,K

+

X

F ∈Eh

τF kJqKk20,F

1/2

,

and we present an interpolation result in this norm. Lemma 1. Let us suppose that (v, q) ∈ H 1 (Ω)2 × H 1 (Ω). Then, there exists C such that (50)

k(v − π(Ch (v)), q − Πh (q))kh ≤ C h

√

 1 σ |v|1,Ω + √ |q|1,Ω . σ

Proof. From the definition of the norm, (42) and (46) there follows that k(v − π(Ch (v)), q − Πh (q))k2h = σ kv − π(Ch (v))k20,Ω +

X

F ∈Eh

τF kJq − Πh (q)Kk20,F

≤ 2 σ kv − π(v)k20,Ω + kπ(v − Ch (v))k20,Ω ) + (51)

X

F ∈Eh

τF kJq − Πh (q)Kk20,F

  h2 2 2 2 ≤ C σ h |v|1,Ω + |q|1,Ω + 2σ kπ(v − Ch (v))k20,Ω . σ

Next, from its definition it is easy to prove that the Raviart-Thomas operator π satisfies (see [9] for a related result and Lemma 5 for an alternative proof) (52)

kπ(v)k0,Ω ≤ C (kvk0,Ω + h |v|1,Ω ) ,

for all v ∈ H 1 (Ω)2 , and then the result follows applying (52) and (42) in (51).



Before heading to stability, an auxiliary result is stated next. Lemma 2. Let π be the Raviart-Thomas interpolator, then there exists a positive constant C1 such that, for all v 1 ∈ Vh and q0 ∈ Qh , it holds (53)

(ℓ(Jq0 K), σ π(v 1 ))Ω ≤ C1

where α := max{αF : F ∈ Eh }.

(

X

F ∈Eh

τF kJq0 Kk20,F

) 21

√

σα h kπ(v 1 )k0,Ω ,

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

15

Proof. In order to prove (53) we first note that ℓ(Jq0 K) ∈ VRT0 . Hence, using successively

the Cauchy-Schwarz inequality, (11) and kϕF k0,K ≤ C1 hF we get (ℓ(Jq0 K), σ π(v 1 ))Ω ≤

X

K∈Th

kℓK (Jq0 K)k0,K σ kπ(v 1 )k0,K

X X αF Z ≤ |Jq0 K| kϕF k0,K σ kπ(v1 )k0,K σ F K∈Th F ⊆∂K X X 1 τF hF2 kJq0 Kk0,F σ kπ(v 1 )k0,K ≤ C1 K∈Th F ⊆∂K

≤ C1

(

X

F ∈Eh

τF kJq0 Kk20,F

) 21

h

√

σ α kπ(v 1 )k0,Ω ,

and the result follows.



We are ready to prove PGEM are well-posed. Lemma 3. The bilinear forms Bs (., .) defined in (17) satisfies Bs ((v 1 , q0 ), (v 1 , −q0 )) = k(π(v1 ), q0 )k2h Moreover, assuming α ≤

1 C12 h2

∀(v 1 , q0 ) ∈ Vh × Qh .

where C1 is the positive constant from Lemma 2, the bilinear

form B(., .) defined in (14) satisfies, 1 k(π(v1 ), q0 )k2h 2

B((v 1 , q0 ), (v 1 , −q0 )) ≥

∀(v 1 , q0 ) ∈ Vh × Qh .

Hence, the problems (13) and (16) are well-posed. Proof. The first equality follows directly from the definition of the bilinear form Bs (., .). For the second one, we recall that (54)

B((v1 , q0 ), (v 1 , −q0 )) = σ kπ(v 1 )k2Ω + (ℓ(Jq0 K), σ π(v 1 ))Ω +

X

F ∈Eh

τF kJq0 Kk20,F .

As for the second term, we use Lemma 2 to obtain (ℓ(Jq0 K), σ π(v 1 ))Ω ≤ C1 ≤

(

X

F ∈Eh

τF kJq0 Kk20,F

) 21

h

√

σα kπ(v 1 )k0,Ω

X 1 1 C12 α h2 τF kJq0 Kk20,F + σ kπ(v1 )k20,Ω , 2 2 F ∈E h

16

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

Hence by successively applying inequality above into (54) and assuming α ≤ B((v1 , q0 ), (v 1 , −q0 )) ≥

1 C12 h2

it holds

X C 2 α h2 σ (1 − 1 kπ(v1 )k20,Ω + ) τF kJq0 Kk20,F 2 2 F ∈E h

1 ≥ k(π(v1 ), q0 )k2h . 2 The continuity of bilinear forms is straightforward and thus the well-posedness of (13) and (16) stems from the Necas Theorem [17].



Remark. We remark that the metric k(π(v 1 ), q0 )kh defines a norm in the space Vh ×Qh since

the Raviart-Thomas interpolation operator is injective when restricted to Vh . On the other hand, we remark that the hypothesis on α is not really restrictive, since it only applies for coarse meshes. For a sufficiently refined mesh the choice for α is essentially unlimited.  Neither of the methods proposed in the previous section are formally consistent as points out the next result. Lemma 4. Let (u, p) ∈ H0div (Ω) × [H 1 (Ω) ∩ L20 (Ω)] be the weak solution of (3), (u1 , p0 ) the solution of (13) and (ˆ u1 , pˆ0 ) the solution of (16), respectively. Then, X
(σ MuK (u), σ πK (v 1 ))K , B (u − u1 , p − p0 ), (v 1 , q0 ) = − K∈Th




ˆ 1 , p − pˆ0 ), (v 1 , q0 ) = − Bs (u − u

for all (v 1 , q0 ) ∈ Vh × Qh .

X

K∈Th

(σ MuK (u), σ πK (v 1 ))K ,

Proof. The result follows from the definition of B(., .) and Bs (., .), and noting that JpK = 0 a.e. across all the internal edges



3.2. Error estimates for the symmetric formulation. We begin this section by proving the following technical result concerning the operator MuK . Lemma 5. Let v ∈ H 1 (K)2 . Then, there exists a constant C such that kMuK (v)k0,K ≤ C σ −1 hK |v|1,K . Proof. Let v ∈ H 1 (K)2 , and w := MuK (v). Then, from the definition of MuK (cf. (24)), w

satisfies (55)

σ w + ∇ξ = v ,

σ ∇ · w = ∇· v − ΠK (∇· v)

σ w · n = v · n − ΠF (v · n)

in K ,

on each F ⊆ ∂K ,

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

17

where ξ ∈ L20 (K). Multiplying the first equation by w and integrating it by parts we arrive at

σ kwk20,K = (v, w)K + = (v, w)K + −

1 X 1 (ξ, v · n − ΠF (v · n))F (ξ, ∇· v − ΠK (∇· v))K − σ σ F ⊆∂K 1 (ξ − ΠK (ξ), ∇· v − ΠK (∇· v))K σ

1 X (ξ − ΠF (ξ), v · n − ΠF (v · n))F σ F ⊆∂K

≤ kvk0,K kwk0,K + +

(56)

1 kξ − ΠK (ξ)k0,K k∇· v − ΠK (∇· v)k0,K σ

1 X kξ − ΠF (ξ)k0,F kv · n − ΠF (v · n)k0,F . σ F ⊆∂K

Next, from (55) it holds k∇ξk0,K ≤ kvk0,K + σ kwk0,K and hence, using the local trace result (45), the approximation property of the projection operators ΠF and ΠK , and the inequality

above to obtain C hK |ξ|1,K |v|1,K σ C hK ≤ kvk0,K kwk0,K + (kvk0,K + σ kwk0,K ) |v|1,K σ σ ≤ C σ −1 (kvk20,K + h2K |v|21,K ) + kwk20,K , 2

σ kwk20,K ≤ kvk0,K kwk0,K +

(57)

and then we have proved that kMuK (v)k0,K ≤ C σ −1 (kvk0,K + hK |v|1,K ) .

(58)

Finally, let us denote v 0 = ΠK (v). Since v 0 is a constant in each element, there holds that MuK (v 0 ) = 0 and then, from (58) it follows that kMuK (v)k0,K = kMuK (v − v 0 )k0,K ≤ C σ −1 (kv − v 0 k0,K + hK |v|1,K ) , and the result follows using the approximation properties of the projection.



The previous lemma results in an alternative proof of the following classical interpolation error estimate: Corollary 6. There exists C such that kv − πK (v)k0,K ≤ C hK |v|1,K , for all v ∈ H 1 (K)2 .

18

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

Proof. The result follows from the previous lemma and the fact that v − πK (v) = σMuK (v).



Theorem 7. Let (u, p) ∈ H0div (Ω) ∩ H 2 (Ω)2 × H 1 (Ω) ∩ L20 (Ω) be the solution of (2), and

ˆ 1 , pˆ0 ) ∈ Vh × Qh the solution of method (16). Then, defining u ˆ h := π(u ˆ 1 ), the following (u error estimate holds


√ 1 ˆ h , p − pˆ0 )kh ≤ C h ( σ kuk2,Ω + √ |p|1,Ω . k(u − u σ

(59)

Proof. Let (v 1 , q0 ) = (Ch (u), Πh (p)). From the triangle inequality we have

(60)

ˆ 1 ), p − pˆ0 )kh ≤ k(u − π(v 1 ), p − q0 )kh + k(π(v 1 − u ˆ 1 ), q0 − pˆ0 )kh . k(u − π(u

The first term is easily estimated using Lemma 1. Next, let us estimate the second term on the right hand side. For that, we use the coercivity of Bs (., .) (cf. Lemma 3) and the consistency result (cf. Lemma 4) to obtain

ˆ 1 ), q0 − pˆ0 )k2h = Bs ((v 1 − u ˆ 1 , q0 − pˆ0 ), (v 1 − u ˆ 1 , pˆ0 − q0 )) k(π(v 1 − u X ˆ 1 ))K ˆ 1 , pˆ0 − q0 )) − σ 2 (MuK (u), πK (v 1 − u = −Bs ((u − v 1 , p − q0 ), (v 1 − u K∈Th

ˆ 1 ))Ω + (p − q0 , ∇ · (π(v 1 − u ˆ 1 )))Ω + (ˆ = −σ (π(u − v 1 ), π(v 1 − u p0 − q0 , ∇ · (π(u − v 1 )))Ω X X ˆ 1 ))K . σ 2 (MuK (u), πK (v 1 − u τF (Jp − q0 K, Jˆ p0 − q0 K)F − + F ∈Eh

K∈Th

Next, from the properties of the projection operator Πh it holds

(61)

ˆ 1 )))Ω = 0 , (p − q0 , ∇ · (π(v 1 − u

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

19

and then, integrating by parts and using the Cauchy-Schwarz’s inequality we arrive at ˆ 1 ), q0 − pˆ0 )k2h ≤ σ kπ(u − v 1 )k0,Ω kπ(v 1 − u ˆ 1 )k0,Ω + k(π(v 1 − u X

+

F ∈Eh

τF kJp − q0 Kk0,F kJˆ p0 − q0 Kk0,F + σ 2

ˆ 1 )k0,Ω + C ≤ σ kπ(u − v 1 )k0,Ω kπ(v1 − u X

+

F ∈Eh

+

F ∈Eh

≤ C

F ∈Eh

σ kπ(u − v 1 )k20,Ω + σ h2 |u|22,Ω + X

K∈Th

ˆ 1 ))K . (MuK (u), πK (v 1 − u

X 1

F ∈Eh

X

F ∈Eh

F ∈Eh

(Jˆ p0 − q0 K, (u − v 1 ) · n)F

ˆ 1 ))K (MuK (u), πK (v 1 − u 3

kJˆ p0 − q0 Kk0,F hF2 |u|2,ωF

K∈Th

X

τF kJp − q0 Kk0,F kJˆ p0 − q0 Kk0,F + σ 2

(

+ σ2

K∈Th

X

τF kJp − q0 Kk0,F kJˆ p0 − q0 Kk0,F + σ 2

ˆ 1 )k0,Ω + C ≤ σ kπ(u − v 1 )k0,Ω kπ(v1 − u X

X

X

ˆ 1 ))K (MuK (u), πK (v 1 − u 1

σ 2 τF2 kJˆ p0 − q0 Kk0,F hF |u|2,ωF

X

K∈Th

ˆ 1 ))K (MuK (u), πK (v 1 − u

τF kJp − q0 Kk20,F

) 12

ˆ 1 ), q0 − pˆ0 )kh k(π(v 1 − u

The last term on the right hand side is estimated next. Using Lemma 5 we may bound the consistency term as follows X

K∈Th

ˆ 1 ))K ≤ (MuK (u), πK (v 1 − u

(62)

X

K∈Th

ˆ 1 )k0,K kMuK (u)k0,K kπK (v 1 − u

ˆ 1 )k0,Ω . ≤ C σ −1 h |u|1,Ω kπ(v1 − u

Collecting all the above results it follows from Lemma 1 that √  1 h ˆ 1 ), q0 − pˆ0 )k2h ≤ C ˆ 1 ), pˆ0 − q0 )kh , σhkuk2,Ω + √ |p|1,Ω k(π(v1 − u k(π(v1 − u 2 σ and hence (59) follows.



Next, we analyze the error in the H div (Ω) norm for the velocity and the L2 (Ω) norm for the pressure.

20

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

ˆ 1 , pˆ0 ) be the solutions of (3) and (16), respectively, and Theorem 8. Let (u, p) and (u ˆ h := π(u ˆ 1 ). Under the hypothesis of Theorem 7, there exists C such that u
1 ˆ h kdiv,Ω ≤ C h kuk2,Ω + |p|1,Ω , (63) ku − u σ
1 ˆ h − ℓ(Jˆ ku − u p0 K)kdiv,Ω ≤ C h ( kuk2,Ω + |p|1,Ω , (64) σ
(65) kp − pˆ0 k0,Ω ≤ C h σ kuk2,Ω + |p|1,Ω . Proof. First, let v 1 := Ch (u); then

ˆ 1 )k20,Ω = (∇ · (u − u ˆ 1 ), ∇ · (u − u ˆ 1 ))Ω k∇ · (u − u ˆ 1 ), ∇ · (u − v 1 ))Ω + (∇ · (u − u ˆ 1 ), ∇ · (u ˆ 1 − v 1 ))Ω . = (∇ · (u − u

(66)

ˆ 1 − v 1 ) ∈ Qh and v 1 = 0) we get Next, from Lemma 4 (considering q0 := ∇ · (u X ˆ 1 − v 1 )K)F ˆ 1 ), ∇ · (u ˆ 1 − v 1 ))Ω = − τF (Jp − pˆ0 K, J∇· (u (∇ · (u − u F ∈Eh

≤ ≤

(67)

X

F ∈Eh

ˆ 1 − v 1 )Kk0,F τF kJp − pˆ0 Kk0,F kJ∇· (u

X τF X ˆ 1 − v 1 )Kk20,F . τF kJ∇· (u kJp − pˆ0 Kk20,F + γ γ F ∈E F ∈E h

h

Next, using the local trace result (45), (42) and the mesh regularity to obtain i X α hK h X 2 ˆ 1 − v 1 )Kk20,F ≤ C γ ˆ τF kJ∇ · (u γ h−1 k∇ · ( u − v )k 1 1 0,K K σ K∈Th F ∈Eh i Cγ X h ˆ 1 )k20,K + k∇ · (u − v 1 )k20,K k∇ · (u − u ≤ σ K∈T h i Cγ X h ˆ 1 )k20,K + Ch2K |u|22,ωK . ≤ k∇ · (u − u (68) σ K∈T h

Hence, choosing γ =

σ 4C

in (68) and using ab ≤ (a2 /4) + b2 , the mesh regularity and (42)

again, (66) and (67) become ˆ 1 )k20,Ω ≤ k∇ · (u − u

1 ˆ 1 )k20,Ω + k∇ · (u − v 1 )k20,Ω k∇ · (u − u 4 i X τF Cγ X h ˆ 1 )k20,K + Ch2K |u|22,ωK k∇ · (u − u + kJp − pˆ0 Kk20,F + γ σ K∈T F ∈E h

h

≤ C σ −1 (69)

X

F ∈Eh

τF kJp − pˆ0 Kk20,F + C h2 |u|22,Ω +

ˆ h , p − pˆ0 )k2h + C h2 |u|22,Ω + ≤ C σ −1 k(u − u

1 ˆ 1 )k20,Ω k∇ · (u − u 2

1 ˆ 1 )k20,Ω , k∇ · (u − u 2

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

21

and the result follows by applying Theorem 7 and extracting the square root. Now, we use the local mass conservation feature to prove (64). In fact, we get Z ˆ 1 − ℓK (Jˆ ∇ · (u − u p0 K)) = 0 , K

ˆ 1 + ℓK (Jˆ and then, since ∇ · (u p0 K))|K ∈ R, we end up with ˆ 1 − ℓK (Jˆ ˆ 1 − ℓK (Jˆ ˆ 1 − ℓK (Jˆ k∇ · (u − u p0 K))k20,K = (∇ · (u − u p0 K)), ∇ · (u − u p0 K)))K ˆ 1 − ℓK (Jˆ = (∇ · (u − u p0 K)), ∇ · (u − π(u)))K ˆ 1 − ℓK (Jˆ ≤ Ch |∇ · u|1,K k∇ · (u − u p0 K)k0,K .

(70)

As seen in Lemma 2 we can prove that (71)

kℓ(Jˆ p0 K)k0,Ω ≤

X

F ∈Eh

τF kϕF k0,Ω |Jˆ p0 K| ≤ C h σ

− 21

(

X

F ∈Eh

τF kJp − pˆ0 Kk20,F

) 12

,

and then using (59) we obtain ˆ h − ℓ(Jˆ ku − u p0 K)k0,Ω ≤ C h ( kuk2,Ω +

(72)

and (64) follows from (70) and (72).


1 |p|1,Ω , σ

Finally, we consider the estimate for the pressure. From the continuous inf-sup condition (see [13]), there exists a w ∈ H01 (Ω) such that ∇· w = p − pˆ0 in Ω and kwk1,Ω ≤ C kp − pˆ0 k0,Ω .

(73)

Let w 1 = Ch (w). Since ∇ · w1 = ∇ · (π(w 1 )), using Lemma 4 and recalling that πK =

I − σMuK , we obtain

kp − pˆ0 k20,Ω = (∇ · w, p − pˆ0 )Ω = (∇ · (w − w 1 ), p − pˆ0 )Ω + (∇ · w 1 , p − pˆ0 )Ω X [(w − w 1 , ∇p)K + (w − w1 , (p − pˆ0 ) I· n)∂K ] = K∈Th

+

X

K∈Th

≤C

h X

σ ((I − σMuK )(u) − πK (u1 ), πK (w1 ))K + σ 2 (MuK (u), πK (w 1 ))K

K∈Th

h X

K∈Th

h2K |p|21,K + σ k(u − π(u1 ), p − pˆ0 )k2h

2 h−2 K kw − w 1 k0,K +

X

F ∈Eh

i 12

·

2 2 h−1 F kw − w 1 k0,F + kπ(w 1 )k0,Ω

i 21

.

22

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

Now, using (42)-(43), the regularity of the mesh, (52) and (73) we obtain i1/2 h X X −1 2 2 −2 2 τF kw − w1 k0,F + kπ(w1 )k0,Ω ≤ C |w|1,Ω ≤ C kp − pˆ0 k0,Ω , hK kw − w 1 k0,K + F ∈Eh

K∈Th

hence, dividing by kp − pˆ0 k0,Ω and using Theorem 7, we arrive at h i h i 21 ˆ 1 ), p − pˆ0 )k2h + σ 2 h2 |p|21,Ω ≤ C h σ kuk2,Ω + |p|1,Ω , kp − pˆ0 k0,Ω ≤ C σ k(u − π(u

and the result follows.



Remark. Having assumed the coefficient σ to be constant in Ω (and then independent of any small scale), the H 2 (Ω)-norm of the exact solution u does not blow up as is usual in Darcy problems with highly oscillating coefficients (see [15] for further details).  3.3. An error estimate for the method (13). We end the error analysis by proving the following error result concerning the method (13). Theorem 9. Let (u, p) ∈ H0div (Ω) ∩ H 2 (Ω)2 × H 1 (Ω) ∩ L20 (Ω) be the solution of (2), and ˆ 1 , pˆ0 ) the solution of methods (13) and (16), respectively. Then, defining uh := (u1 , p0 ), (u ˆ h := π(u ˆ 1 ), the following error estimate holds π(u1 ) and u ˆ h , p − pˆ0 )kh . k(u − uh , p − p0 )kh ≤ 3 k(u − u Proof. First, from Lemmas 3 and 4 and Lemma 2 it follows 1 ˆ 1 ),p0 − pˆ0 )k2h ≤ B((u1 − u ˆ 1 , p0 − pˆ0 ), (u1 − u ˆ 1 , p0 − pˆ0 )) k(π(u1 − u 2 X ˆ 1 ))K ˆ 1 , p − pˆ0 ), (u1 − u ˆ 1 , p0 − pˆ0 )) + (σ MuK (u), σ πK (u1 − u = B((u − u K∈Th

ˆ 1 , p − pˆ0 ), (u1 − u ˆ 1 , pˆ0 − p0 )) + = Bs ((u − u +

X

K∈Th

=

X

K∈Th

≤ C1

K∈Th

ˆ 1 ))K (ℓK (Jp − pˆ0 K), σ πK (u1 − u

ˆ 1 ))K (ℓK (Jp − pˆ0 K), σ πK ((u1 − u

(

X

F ∈Eh

≤ C12 α h2 (74)

ˆ 1 ))K σ (MuK (u), σ πK (u1 − u

X

τF kJp − pˆ0 Kk20,F

X

F ∈Eh

) 21

√

τF kJp − pˆ0 Kk20,F +

ˆ h , p − pˆ0 )k2h + ≤ k(u − u

ˆ 1 )k0,Ω σα h kπ(u1 − u

1 ˆ 1 )k20,Ω σkπ(u1 − u 4

1 ˆ 1 ), p0 − pˆ0 )k20,Ω k(π(u1 − u 4

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

where we used α ≤

1 C12 h2

23

and C1 is the positive constant from Lemma 2. The result follows

applying the triangular inequality.



Remark. We end this section by remarking that the same analysis from Theorem 8 may be ˆ h kdiv,Ω and kp − pˆ0 k0,Ω as well.  carried out to prove error estimates on ku − u 4. Numerical experiments Now, we are interested in the numerical validation of the PGEM in its symmetric version (16). The method (13) behaves similarly as shown in Section 3.3. Two numerical tests with available analytical solutions are performed and the theoretical results validated. The assumed vanishing boundary condition to generate the methods is adopted by the first numerical test, a property which is no longer shared by the second case. In all the computations the value for αF has been set to one. 4.1. An analytical problem. The domain is Ω = (0, 1) × (0, 1) and we set σ = 1 and the exact pressure is given by p(x, y) = cos(2πx) cos(2πy). Next, the exact velocity is

determined from the Darcy law and the boundary condition is taken to be its normal component on the boundary, thus b = 0. Consequently, the divergence velocity field is set as g = 8 π 2 cos(2πx) cos(2πy). MESH

Figure 2. Mesh for the analytical problem. In Figures 4-6 we report the errors on velocity and pressure in a sequence of structured meshes. One observes optimal convergence of all quantities as h → 0 in their respective

24

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

natural norms, which are in accordance with the theoretical results. Moreover, in Figure 4 we also plot the error kπ(u) − π(u1 )k0,Ω which is smaller than ku − u1 k0,Ω . Here we denote hP i1/2 2 |Jp − p0 K|h := . F ∈Eh τF kJp − p0 Kk0,F

First, we adopt the structured mesh described in Figure 2 which contains 4096 triangular

elements (h = 3.125 × 10−2 ). We depict in Figure 3 the isolines (free of oscillations) of the pressure and |u1 | obtained from (16).

Furthermore, in Table 4.1 we study the local mass conservation feature of (16) whether we look at either u1 or u1 + ue . We define the quantities

(75)

Me := max

K∈Th

R

K


∇· (u1 + ue ) − g dx |K|

and M1 := max

K∈Th

R

K


∇· u1 − g dx , |K|

and we find a loss of mass, as expected, when just the linear part of the solution is used. Nevertheless, we recover the local mass conservation property updating the linear velocity field by the multi-scale velocity ue . Similar results were obtained in [7] using the nonsymmetric PGEM (32).

h

0.25

0.125

Me 2.3 × 10−14 2.3 × 10−13

6.25 × 10−2 10−12

3.125 × 10−2 10−11

1.5625 × 10−2 1.2 × 10−10

7.8125 × 10−3 10−9

M1 0.81 0.35 0.09 0.03 6.4 × 10−3 1.6 × 10−3 Table 4.1: Relative local mass conservation errors with the symmetric method (16).

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

PRESSURE

VELOCITY

Figure 3. Isolines of the pressure (left) and |u1 | (right) using the symmetric

method (16).

10

error

1

ku − u1k0,Ω kπ(u) − π(u1)k0,Ω h2

0.1 0.01 0.001 1e-04 0.001

0.01

0.1

1

h Figure 4. Convergence history of ku − u1 k0,Ω and kπ(u) − π(u1 )k0,Ω for the

symmetric method (16).

25

26

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

100

100

k∇·(u − u1)k0,Ω

h

10

h

error

error

10

k(u − u1,p − p0)kh k(u − π(u1),p − p0)kh

1

1 0.1

0.01 0.001

0.01

0.1

0.1 0.001

1

0.01

h

0.1

1

h

Figure 5. Convergence history of k∇· (u − u1 )k0,Ω and k(u − u1 , p − p0 )kh , and k(u − π(u1 ), p − p0 )kh for the symmetric method (16).

10

10

kp − p0k0,Ω

h

|Jp − p0K|h

h

1 error

error

1 0.1

0.1 0.01

0.001 0.001

0.01

0.1 h

1

0.01 0.001

0.01

0.1 h

Figure 6. Convergence history of kp−p0 k0,Ω and |Jp − p0 K|h for the symmetric method (16).

1

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

27

4.2. A second analytical problem. The problem is set up as in the first test, differing y3x x3 y − . Once again, the velocity by replacing the previous exact pressure by p(x, y) = 3 3 is computed from the Darcy law and the boundary condition b is taken to be its no vanishing normal component on the boundary. On the other hand, clearly the velocity field is divergence free (g = 0). We validate the symmetric method (16) using a sequence of structured meshes. Optimality is reached whatever the norm is considered as shown in Figures 8-10. Similar quadratic convergence is observed for ku − u1 k0,Ω and kπ(u) − π(u1 )k0,Ω . Table 4.2 highlights that the local mass conservation property is recovered as soon as u1 is updated by u1 + ue .

h

0.25

0.125

Me 3.1 × 10−15 1.5 × 10−14

6.25 × 10−2 1.9 × 10−13

3.125 × 10−2 7.1 × 10−13

1.5625 × 10−2 5.6 × 10−12

7.8125 × 10−3 7.6 × 10−11

M1 0.14 0.08 0.04 0.02 0.01 5.3 × 10−3 Table 4.2: Relative local mass conservation errors with the symmetric method (16). Next, the structured mesh of Figure 2 is once more adopted, and the solution is oscillationfree as it can be seen through the isolines of p0 and u1 in Figure 7.

PRESSURE

VELOCITY

Figure 7. Isolines of the pressure (left) and |u1 | (right) using the symmetric

method (16).

28

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

0.1

ku − u1k0,Ω kπ(u) − π(u1)k0,Ω

error

0.01

h2

0.001

1e-04

1e-05 0.001

0.01

0.1

1

h Figure 8. Convergence history of ku − u1 k0,Ω and kπ(u) − π(u1 )k0,Ω for the

symmetric method (16).

1

1

k∇·(u − u1)k0,Ω

h

h

error

0.1

error

0.1

k(u − u1,p − p0)kh k(u − π(u1),p − p0)kh

0.01

0.001 0.001

0.01

0.01

0.1

1

0.001 0.001

0.01

h

0.1 h

Figure 9. Convergence history of k∇· (u − u1 )k0,Ω and k(u − u1 , p − p0 )kh , and k(u − π(u1 ), p − p0 )kh for the symmetric method (16).

1

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

0.1

0.1

kp − p0k0,Ω

h

29

|Jp − p0K|h h

error

error

0.01 0.01

0.001

1e-04 0.001

0.01

0.1

1

0.001 0.001

0.01

h

0.1

1

h

Figure 10. Convergence history of kp − p0 k0,Ω and |Jp − p0 K|h for the symmetric method (16).

5. Conclusion New enriched finite element methods make the simplest pair of nodal based interpolation spaces stable for the Darcy model. It has been proved that such methods lead to optimal error estimates in natural norms in addition to be locally mass conservative. Such fundamental property is recovered inside a general framework which relies on updating the linear part of velocity with a particular Raviart-Thomas function. Such strategy can prevent other jump-based stabilized finite element methods from local loss of mass while keeping them stable and accurate. Alternatives to deal with higher order interpolations should include additional control on the jumps of gradient of the pressure, in the form of a new enrichment function leading to a term like (ℓ(J∇p1 · nK), ∇q1 )Ω , a feature that may be incorporated into

the current Petrov-Galerkin framework.

Appendix A. The error for general g If we do not suppose that g is a piecewise constant function, instead we admit that g ∈ H 1 (Ω), then a third enrichment function (uge , pge ) appears as the solution of (76)

σ uge + ∇pge = 0 in K, uge · n = 0

∇· uge = g − ΠK (g) in K, in ∂K .

30

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

Now, denoting (uge , pge ) = GK (g − ΠK (g)), then, the original method (13) becomes: Find

˜ 1 , p˜0 ) ∈ Vh × Qh such that for all (v 1 , q0 ) ∈ Vh × Qh it holds (u X ˜ 1 , p˜0 ), (v 1 , q0 )) = Fs (π(v 1 ), q0 ) − (GK (g − ΠK (g)), σ πK (v 1 ))K , (77) B((u K∈Th

where we recall that ΠK (g) =

1 |K|

R

K

g.

Remark. Regarding the symmetric method its right hand side must be enhanced with the same contribution as well, and then the following results also apply to such version.  Remark. Next, we see that there exists C > 0 such that (78)

kuge k0,K ≤ C hK kg − ΠK (g)k0,K

∀ K ∈ Th .

Indeed, we remark that multiplying (77) by uge and integrating it by parts we obtain that 1 kuge k20,K = − (∇pge , uge )K σ 1 g = (pe , ∇· uge )K σ 1 = (pge , g − ΠK (g))K σ 1 ≤ kpge k0,K kg − ΠK (g)k0,K σ hK g ≤C |p |1,K kg − ΠK (g)k0,K σ e = C hK kuge k0,K kg − ΠK (g)k0,K , and the result follows.  Now, applying the Strang Lemma (cf. [17]) and (78) we arrive at P − K∈Th (GK (g − ΠK (g)), σ πK (v 1 ))K ˜ 1 − u1 ), p˜0 − p0 )kh ≤ k(π(u sup k(π(v 1 ), q0 )kh (v1 ,q0 )∈Vh ×Q0h −{0} # 12 " X kGK (g − ΠK (g))k20,K ≤ C K∈Th

≤ C h kg − ΠK (g)k0,Ω ≤ C h2 |g|1,Ω , and then, as claimed, we see that the error is not affected by the fact that we projected g onto the piecewise constant space. Therefore, following the same strategy in §3 we can

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

31

estimate pressure and divergence of velocity in their respective natural norms, and so they are summarized in the following theorem: ˜ h := π(u ˜ 1 ). Then, under the Theorem 10. Let us suppose that g ∈ H 1 (Ω) and denote u

hypothesis of Theorem 9, the following error estimates hold


1 |p|1,Ω + h |g|1,Ω , σ
1 ≤ C h |u|2,Ω + |p|1,Ω + h |g|1,Ω , σ
≤ C h σ |u|2,Ω + |p|1,Ω + σ h |g|1,Ω .

˜ h kdiv,Ω ≤ C h |u|2,Ω + ku − u ˜ h − ℓ(˜ ku − u p0 )kdiv,Ω kp − p˜0 k0,Ω Acknowledgments.

A part of this work was performed during the stay of Gabriel Barrenechea at

the Applied Mathematics Dpt., LNCC, Brazil, in the framework of joint Chile(CONICYT)-Brazil(CNPq) project No. 2005-073 (Chile)-490639/2005-4 (Brazil). The authors want to thank Mark Ainsworth and Rodolfo Rodr´ıguez for helpful discussions and comments.

References [1] J. E. Aarnes and Y. Efendiev, Mixed multiscale finite element methods for stochastic porous media flows, SIAM J. Sci. Comput., 30 (2008), pp. 2319–2339. [2] R. Araya, G. R. Barrenechea, L. P. Franca, and F. Valentin, Stabilization arising from PGEM: A review and further developments, Applied Numerical Mathematics, 227 (2009), pp. 93–101. [3] R. Araya, G. R. Barrenechea, and F. Valentin, Stabilized finite element methods based on multiscale enrichment for the Stokes problem, SIAM J. Numer. Anal., 44 (2006), pp. 322–348. [4] R. Araya and F. Valentin, A multiscale a posteriori error estimate, Comput. Methods Appl. Mech. Engrg., 194 (2005), pp. 2077–2094. [5] T. Arbogast, Analysis of a two-scale locally conservative subgrid upscaling for elliptic problems, SIAM J. Numer. Anal., 42 (2004), pp. 576–598. [6] T. Arbogast and M. F. Wheeler, A family of rectangular mixed elements with a continuos flux for second order elliptic problems, SIAM J. Numer. Anal., 42 (2005), pp. 1914–1931. [7] G. R. Barrenechea, L. P. Franca, and F. Valentin, A Petrov-Galerkin enriched method: a mass conservative finite element method for the Darcy equation, Computer Methods in Applied Mechanics and Engineering, 196 (2007), pp. 2449–2464. [8] G. R. Barrenechea and F. Valentin, Relationship between multiscale enrichment and stabilized finite element methods for the generalized Stokes problem, CRAS, 341 (2005), pp. 635–640. ´ n-Rebollo, F. Hecht, and Z. Mghazli, Mortar finite element discretiza[9] C. Bernardi, T. Chaco tion of a model coupling Darcy and Stokes equations, M2AN, 42 (2008), pp. 375–410. [10] P. Bochev and M. Gunzburger, On least-squares finite element methods for the Poisson equation and their connection to the Dirichlet and Kelvin principles, SIAM J. Numer. Anal., 42 (2005), pp. 340– 362.

32

[11]

G.R. BARRENECHEA, L.P. FRANCA, AND F. VALENTIN

, A locally conservative least-squares method for Darcy flows, Communications in Numerical Methods in Engineering, 24 (2008), pp. 97–110.

[12] F. Brezzi, J. Douglas, and L. Marini, Two families of mixed finite elements for second order elliptic problems, Numer. Math., 47 (1985), pp. 217–235. [13] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, vol. 15 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, New-York, 1991. [14] E. Burman and P. Hansbo, A unified stabilized method for Stokes’ and Darcy’s equations, J. Comput. Appl. Math., 198 (2007), pp. 35–51. [15] Z. Chen and T. Hou, A mixed multiscale finite element method for elliptic problems with oscillating coefficients, Mathematics of Computation, 72 (2002), pp. 541–576. ´ment, Approximation by finite element functions using local regularization, RAIRO Anal. [16] P. Cle Num´er., (1975), pp. 77–84. [17] A. Ern and J.-L. Guermond, Theory and Practice of Finite Elements, Springer-Verlag, 2004. [18] L. P. Franca, A. L. Madureira, L. Tobiska, and F. Valentin, Convergence analysis of a multiscale finite element method for singularly perturbed problems, SIAM Multiscale Model. and Simul., 4 (2005), pp. 839–866. [19] L. P. Franca, A. L. Madureira, and F. Valentin, Towards multiscale functions: enriching finite element spaces with local but not bubble–like functions, Comput. Methods Appl. Mech. Engrg., 194 (2005), pp. 3006–3021. [20] L. P. Franca, J. V. A. Ramalho, and F. Valentin, Multiscale and Residual–Free Bubble functions for reaction–advection–diffusion Problems, International Journal for Multiscale Enginnering, 3 (2005), pp. 297–312. [21]

, Enriched finite element methods for unsteady reaction-diffusion problems, Communications in Numerical Methods in Engineering, 22 (2006), pp. 619–625.

[22] V. Girault and P. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms, vol. 5 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, New-York, 1986. [23] T. J. R. Hughes, A. Masud, and J. Wan, A stabilized mixed discontinuos Galerkin method for Darcy flow, Comput. Methods Appl. Mech. Engrg., 195 (2006), pp. 3347–3381. [24] J. Li, T. Arbogast, and Y. Huang, Mixed methods using standard conforming finite elements, Comput. Methods Appl. Mech. Engrg., 198 (2009), pp. 680–692. [25] K. Mardal, X. Tai, and R. Winther, A robust finite element method for Darcy-Stokes flow, SIAM J. Numer. Anal., 40 (2002), pp. 1605–1631. [26] A. Masud and T. J. R. Hughes, A stabilized mixed finite element method for Darcy flow, Comput. Methods Appl. Mech. Engrg., (2002), pp. 4341–4370. [27] R. Raviart and J. Thomas, A mixed finite element method for 2nd order elliptic problems, Mathematical aspect of finite element methods, no. 606 in Lecture Notes in Mathematics, Springer-Verlag, New York, 1977, pp. 292–315.

A SYMMETRIC NODAL CONSERVATIVE FEM FOR THE DARCY EQUATION

33

Department of Mathematics, University of Strathclyde, 26 Richmond Street, Glasgow G1 1XH, UK E-mail address: [email protected] Department of Mathematical Sciences, University of Colorado at Denver, P.O. Box 173364, Campus Box 170 Denver, Colorado 80217-3364, USA E-mail address: [email protected] ´tica Aplicada, Laborato ´ rio Nacional de Computac ˜ o Cient´ıfica, Departamento de Matema ¸a ´ ´ Av. Getulio Vargas, 333, 25651-070 Petropolis - RJ, Brazil E-mail address: [email protected]