A posteriori error analysis for discontinuous finite volume methods of elliptic interface problems

Journal of Computational and Applied Mathematics 255 (2014) 529–543

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

A posteriori error analysis for discontinuous finite volume methods of elliptic interface problems Lin Mu a,∗ , Rabeea Jari b a

Department of Mathematics, Michigan State University, East Lansing, MI, 48824, United States

b

Department of Applied Science, University of Arkansas at Little Rock, Little Rock, AR, 72204, United States

article

info

Article history: Received 9 March 2012 Received in revised form 3 February 2013 Keywords: A posteriori error estimator Finite volume method Elliptic interface problems

abstract This paper presents two types of a posterior error estimations for discontinuous finite volume methods to elliptic interface problems. One is the residual-based estimator and other one is the recovery-based estimator. For both a posterior error estimators, we establish reliability and efficiency bounds. Numerical experiments are conducted to validate our theoretical conclusions. © 2013 Elsevier B.V. All rights reserved.

1. Introduction The finite volume method (FVM) is a discretization technique for solving partial differential equations. Due to the property of local conservation, FVM is widely used in many fields such as computational fluid dynamics. The notion of a posteriori error estimate and adaptive mesh refinement was first introduced into finite element analysis in 1978 by Babuška and Rheinboldt [1,2]. A posteriori error estimate plays an important role in adaptive mesh refinement procedure. It uses computable quantities such as a numerical solution to estimate the error between the exact solution and the numerical solution. The computable a posteriori error estimator can be used as an indicator to guide mesh refinement to reduce the global error. There are two kinds of a posteriori error estimators: one is the residual-based error estimator (see [3–11]) and other is the recovery-based error estimator (see [12–14]). The recovery-based methods have been widely used in engineering, and studied by many researchers, for example ([15–22]), because of their many appealing features. The recent work of a posteriori error estimate for finite volume methods can be found in [23–28]. In this paper, we will study residual-based and recover-based a posteriori error estimators for the discontinuous finite volume method for elliptic interface problems. The numerical schemes of discontinuous finite volume methods have been introduced in [29,28]. We obtain reliability and efficiency bounds for both the residual-based and recovery-based estimators. Our numerical experiments show that these estimators can be used as an effective indicator to reduce the error and to catch the singularity of problems. A posteriori error analysis is investigated for discontinuous finite element methods for the elliptic interface problem in [30]. This paper is motivated by the work in [30]. The remainder of this article is organized as follows. In Section 2, the elliptic interface problem and its discontinuous finite volume formulation are introduced. In Section 3 we introduce two interpolation operators and their properties respectively. Then in Section 4, the residual-based a posteriori error estimation is analyzed. The recovery procedure and the recoverybased a posteriori error estimator are established in Section 5. Finally, in Section 6, we provide numerical results to support our conclusions.

∗

Corresponding author. Tel.: +1 3479253157. E-mail addresses: [email protected], [email protected] (L. Mu), [email protected] (R. Jari).

0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.05.020

530

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

Fig. 1. Primal and dual partition.

1.1. Preliminaries and the problems In this paper, let Ω be a bounded polygonal domain in R2 with boundary ∂ Ω . We consider the following interface problem

− ∇ · (k(x)∇ u) = f

in Ω

(1)

with Dirichlet boundary condition u=0

on ∂ Ω ,

(2)

where f is a given scalar-valued function, and k(x) is positive and piecewise constant on polygonal subdomains of Ω with possible large jumps across subdomain boundaries (interfaces): k(x) = ki > 0

in Ωi

for i = 1, . . . , n. Here, {Ωi }ni=1 is a partition of the domain Ω with Ωi being an open polygonal domain. Let kmax = maxx∈Ω k(x) and kmin = minx∈Ω k(x). This paper will use standard definitions for the Sobolev spaces H s (D) and their associated inner products (·, ·)s,D , norms ∥ · ∥s,D and seminorms | · |s,D , for s ≥ 0, where D is a bounded polygonal domain. The space H 0 (D) coincides with L2 (D), with inner products (·, ·)D . When D = Ω , we will use (·, ·) to denote the inner product. Next, we define two operators. For a vector-valued function τ = (τ1 , τ2 )t , define the divergence operator as follows

∇ · τ :=

∂τ2 ∂τ1 + . ∂ x1 ∂ x2

For a scalar-valued function v , define the operator ∇ ⊥ by

  ∂v ∂v t ∇ ⊥v = − , . ∂ x2 ∂ x1 Moreover, another Hilbert space is defined as H (div; Ω ) = {τ ∈ L2 (Ω )2 : ∇ · τ ∈ L2 (Ω )} equipped with the norm

 ∥τ∥H (div;Ω ) = ∥τ∥

2 0 ,Ω

+ ∥∇ · τ∥

2 0,Ω

 21

.

Let Th be a finite element triangulation of the domain Ω . Assume that the triangulation Th is regular; i.e. diam(K ) ≤ κρK , ∀K ∈ Th . Here diam(K ) is the diameter of the element K and ρK denotes the diameter of the largest inscribed circle of K . Moreover, we will assume the interfaces do not cut through any element K ∈ Th . For a given triangulation Th , its dual partition Th∗ is the union of smaller triangles T1 , T2 , and T3 . These smaller triangles are formed by connecting the barycenter and the three corners of the triangles (shown as Fig. 1 T1 = ∆A2 CA3 , T2 = ∆A3 CA1 , T3 = ∆A1 CA2 ). The trial function space associated with Th for the discontinuous finite volume method is defined as Vh = {v ∈ L2 (Ω ) : v|K ∈ P1 (K ), ∀K ∈ Th }. The test function space is defined by Qh = {p ∈ L2 (Ω ) : p|T ∈ P0 (T ), ∀T ∈ Th∗ }. Let Eh be the union of edges of triangles K ∈ Th , and Eh0 := Eh \ ∂ Ω be the interior edges. Similarly, for e ∈ Eh , we can define (·, ·)s,e and ∥ · ∥s,e with s ≥ 0.

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

531

1.2. Jumps and averages Let e be an interior edge shared by elements K1 and K2 in Th , k1 and k2 be the diffusion coefficients on K1 and K2 , and n1 , n2 be the unit normal vectors on e exterior to K1 and K2 respectively. Define jump [·] on edge e for scalar q and vector v by

[q]e := q|∂ K1 n1 + q|∂ K2 n2 ,

and

[v]e := v|∂ K1 · n1 + v|∂ K2 · n2 .

If e is an edge on the boundary of Ω , with e ∈ ∂ K ∩ ∂ Ω , we define [q]e = q|∂ K n, and [v]e = v|∂ K · n. Let ω1 ∈ [0, 1] and ω2 ∈ [0, 1] be the two weights defined on the edge e ∈ Eh0 satisfying

ω1 + ω2 = 1.

(3)

Define the following weighted averages {·} on edge e for scalar q and vector v by



ω1 q|∂ K1 + ω2 q|∂ K2 , q|∂ K ,

e ∈ Eh0 e ∈ ∂K ∩ ∂Ω,



ω2 q|∂ K1 + ω1 q|∂ K2 , q|∂ K ,

e ∈ Eh0 e ∈ ∂K ∩ ∂Ω,



ω1 v|∂ K1 + ω2 v|∂ K2 , v|∂ K ,

e ∈ Eh0 e ∈ ∂K ∩ ∂Ω,



ω2 v|∂ K1 + ω1 v|∂ K2 , v|∂ K ,

e ∈ Eh0 e ∈ ∂K ∩ ∂Ω.

{q}ω = ω

{q} = {v}ω = and

{v}ω =

Based on the piecewise constant k(x), there are three choices of the weights ω1 , ω2 as ω1,i and ω2,i (i = 1, 2, 3):

√ ω1,1 =

1 2

,

ω1 , 2 =

k2 k1 + k2

,

ω1,3 = √

k2

√

k1 +

k2

√

ω2,1 =

1 2

,

ω2 , 2 =

k1 k1 + k2

,

ω2,3 = √

k1 k1 +

√ .

(4)

k2

A simple calculation leads to the following identity:

 K ∈Th

qv · nds =

e∈∂ K

 e∈Eh

[q]e · {v}ω ds + e

 e∈Eh0

{q}ω [v]e ds.

(5)

e

When there is no ambiguity, the subscript e in the designation will be dropped. Let We = {k}ω be the weighted average of k on edge e ∈ Eh0 . For the above three choices of weights ω1 and ω2 , We is then the arithmetic, the harmonic, and the geometric averages: We,1 =

k1 + k2 2

,

We,2 =

2k1 k2 k1 + k2

,

and

We,3 =



k1 k2 ,

(6)

respectively. It is simple to prove that We,2 ≤ We,3 ≤ We,1 . For the boundary edge e ∈ ∂ K ∩ ∂ Ω , set

ω = 1,

and We = k|K .

It is easy to check that

(ω1,i )2 k1 We

≤ 1,

and

(ω2,i )2 k2 We

≤ 1,

for the above three choices of ω1 , ω2 and We .

(i = 1, 2, 3)

(7)

532

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

2. Finite volume approximation Multiplying both sides of (1) by a test function p ∈ Qh , and integrating by parts gives

(f , p) = −

 ∂T

T ∈Th∗

=−

 ∂K

K ∈Th

(k∇ u · n)pds (k∇ u · n)pds −

3   K ∈Th j=1

(k∇ u · n)pds. Aj+1 CAj

Using (5), we obtain

 K ∈Th

∂K



(k∇ u · n)pds =

e∈Eh0

[k∇ u]{p}ω ds + e

 e∈Eh

{k∇ u}ω · [p]ds. e

Because of the continuity of the flux [31], we have [k∇ u] = 0 on interior edge, thus we get

−

3   K ∈Th j=1

(k∇ u · n)pds − Aj+1 CAj

 e∈Eh

{k∇ u}ω · [p]ds = (f , p). e

Let V (h) = Vh + H 2 (Ω ) ∩ H01 (Ω ) and define an operator γ : V (h) → Qh :

γ v(x)|T =

1 he



v|T ds,

∀T ∈ Th∗ .

(8)

e

Through simple calculation, the mapping γ satisfies the following properties: P1. Quadrature-like and restriction assumptions for γ :



(v − γ v)dx = 0, ∀v ∈ Vh , ∀K ∈ Th K (v − γ v)ds = 0, ∀v ∈ H 2 (Th ), ∀e ∈ ∂ K , ∀K ∈ Th e

if [v] = 0,

then [γ v] = 0, ∥[γ v]∥0,e ≤ ∥[v]∥0,e , ∀e ∈ Eh .

P2. Approximation property of γ :

∥γ w − w∥0,K ≤ ChK |w|1,K ,

∀ K ∈ Th .

Define the bilinear form by ah (v, w) = Ah (v, w) −

 e∈Eh

[γ w] · {k∇v}ω ds + e

 e∈Eh

[γ v] · {k∇w}ω ds + δ e



We [γ v]e · [γ w]e

(9)

e∈Eh

with Ah (v, w) = −

3  

k(x)∇v · nγ w ds, Aj+1 CAj

K ∈Th j=1

and δ being a positive constant to be determined. Then the finite volume approximation is to find uh ∈ Vh such that ah (uh , v) = (f , γ v),

∀v ∈ Vh .

(10)

In the following context, we will use {v} to replace {v}ω when there is no ambiguity. It is easy to see that u, the solution of (1), satisfies ah (u, v) = (f , γ v),

∀v ∈ Vh .

(11)

We define a norm ||| · ||| for V (h) as follows:

|||v|||2 = ∥k1/2 ∇v∥20,h +

 e

We [γ v]2e ,

(12)

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

where ∥k1/2 ∇v∥0,h =

 

K ∈Th

∥k1/2 ∇v∥20,K

1/2

533

. It is easy to prove that ||| · ||| defines a norm. Here and thereafter, we use

C with or without subscripts in this paper to denote a generic positive constant, possibly different at different occurrences, which is independent of the mesh parameter hK and the ratio kmax /kmin but may depend on the domain Ω . Finally, the discrete gradient and divergence is defined piecewisely as follow:

(∇h v)|K := ∇(v|K ),

(∇h · τ)|K := ∇ · (τ|K )

for K ∈ Th , respectively. Lemma 2.1. For any v, w ∈ V (h), we have Ah (v, w) = (k∇h v, ∇h w) +

 ∂K

K

where (k∇h v, ∇h w) =



K ∈Th

(γ w − w)k∇v · nds +



(∇ · k∇v, w − γ w)K ,

(13)

K

(k∇v, ∇w)K .

The proof of Lemma 2.1 can be found in [29]. Furthermore, the following trace inequality can be found in [32]. For

w ∈ H 2 (K ) and for an edge e of K , the following inequality is true   1 2 2 ∥w∥20,e ≤ C h− K ∥w∥0,K + hK ∥∇w∥0,K ,

(14)

where C depends only on the minimum angle of K . Lemma 2.2. For v, w ∈ Vh , we have ah (v, w) ≤ C |||v||||||w|||.

(15)

Additionally, for v ∈ Vh , ah (v, v) ≥ |||v|||2 .

(16)

Proof. Cauchy–Schwartz inequality gives



[γ v] · {k∇w}ω ds ≤ C



e

e∈Eh



We [γ v]

2 e

1/2 

e∈Eh

e∈Eh

he We

2 0 ,e

∥{k∇w}∥

1/2 

.

For e ∈ Eh0 , let e = ∂ K1 ∩ ∂ K2 , k∇v is constant on each element, then we can obtain he We

∥{k∇w}∥20,e ≤ 2h2e



(ω1 )2 k21 (∇w|K1 )2 We

+

(ω2 )2 k22 (∇w|K2 )2



We

  (ω1 )2 k1 (ω2 )2 k2  , k1 (∇w|K1 )2 + k2 (∇w|K2 )2 ≤ 2h2e max We We   2 2 2 ≤ 2he k1 (∇w|K1 ) + k2 (∇w|K2 ) . 

Similarly, e ∈ Eh ∩ ∂ Ω , since {k∇w} = (k∇w)|K , we get he We

∥{k∇w}∥20,e = h2e kK (∇w|K )2 .

Here we will borrow the idea from [33]: for any element K ∈ Th , let {λi }3i=1 be the standard barycentric coordinates of   K , SK = SiK,j 3×3 , and SiK,j = (∇λi , ∇λj )K . Denote the maximum of the largest eigenvalue of SK over Th by ρT . It is proved h in [30] that



h2e kK (∇wK )2 ≤ 4ρT (k∇h w, ∇h w)K . h

e∈∂ K

Combining the above conclusions, we can obtain

 he e∈Eh

We

∥{k∇w}∥20,e ≤ 4



h2e kK (∇w|K )2

K ∈Th e∈∂ K

≤ 16ρTh (k∇h w, ∇h w) = 16ρTh ∥k1/2 ∇w∥20,h .

534

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

The definition of ah (v, w), Lemma 2.1 and the inequalities above imply that ah (v, w) ≤ C



1/2

∥k

∇v∥0,h ∥k

1/2

∇w∥0,h +



1/2

∥k

∇v∥0,K

1/2 



∥k1/2 ∇w∥0,K

1/2 

We [γ v]2e

1/2

+δ

e∈Eh

K ∈Th

We [γ w]

e∈Eh

K ∈Th

+

1/2 2 e



We [γ v]2e

1/2 

e∈Eh

We [γ w]2e

1/2 

e∈Eh

≤ C |||w||| |||v|||. Replacing v ∈ Vh in (9), using Lemma 2.1, we can obtain ah (v, v) = Ah (v, v) + δ



We [γ v]2e .

e∈Eh

By choosing δ ≥ 1, we get (16) directly.



3. Oswald and Clement-type interpolations For a given function v ∈ Vh , define the Oswald interpolation operator Ih : Vh → Vh ∩ H01 (Ω ). The following Lemma can be found in [30]. Lemma 3.1. For any v ∈ Vh , there exists a positive constant C independent of the ratio kmax /kmin such that kK ∥v − Ih v∥

2 0,K

≤C



he We ∥[v]∥

2 0,e



,

(17)

e∈∂ K

and that 1/2

∥kK ∇(v − Ih v)∥20,K ≤ C



We

e∈∂ K

he

∥[v]∥20,e

 (18)

for all K ∈ Th . Furthermore, we define the Clement-type interpolation operator Jh : H01 (Ω ) → Vh ∩ H01 (Ω ) and cite the following Lemma in [30]. Lemma 3.2. For any K ∈ Th and v ∈ H01 (Ω ), the following estimates: −1/2

∥v − Jh v∥0,K ≤ ChK kK

∥k1/2 ∇v∥0,∆K

(19)

and −1/2

∥∇(v − Jh v)∥0,K ≤ CkK

∥k1/2 ∇v∥0,∆K

(20)

hold, where ∆K is the union of all elements that share at least one vertex with K . For any e ∈ Eh0 and v ∈ H01 (Ω ), the following estimate:

 ∥v − Jh v∥0,e ≤ C

he We,1

∥k1/2 ∇v∥0,∆e

holds, where ∆e is the union of all elements that share at least one vertex with edge e. 4. Residual-based a posteriori error estimators In this section, we study the following residual-based a posteriori error estimator:

ηR =



ηR2,K

1/2

,

K ∈Th

where the local indicator on K ∈ Th is defined by

ηR2,K = ηR2f ,K + ηJ2σ ,K + ηJ2u,K

(21)

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

535

with

ηR2f ,K = ηJ2u ,K =

h2K ∥f ∥20,K kK

 We e∈∂ K

he

,

ηJ2σ ,K =

 he ∥[k∇ uh ]∥20,e

1 2

We

e∈∂ K ∩Eh0

∥[uh ]∥20,e .

To analyze the estimator ηR , we employ a standard technique that uses the Helmholtz decomposition: for any given vector-valued function τ ∈ L2 (Ω )2 , there exist α ∈ H01 (Ω ) and β ∈ H 1 (Ω ) such that

τ = k(x)∇α + ∇ ⊥ β.

(22)

Integrating by parts gives

(∇α, ∇ ⊥ β) = 0. Let eh = u − uh , then there exists α ∈ H01 (Ω ) and β ∈ H 1 (Ω ) such that

(23)

k∇h eh = k∇α + ∇ ⊥ β

(24)

and that

∥k1/2 ∇ eh ∥20,h = ∥k1/2 ∇α∥20,Ω + ∥k−1/2 ∇ ⊥ β∥20,Ω .

(25)

Theorem 4.1. Let u ∈ H01 (Ω ) and uh ∈ Vh be the solutions of (1) and (10), respectively. Then there exists a constant C independent of h such that the residual-based a posteriori error estimator satisfies the following global reliability bound:

|||u − uh ||| ≤ C ηR .

(26)

Proof. We have the fact ah (u, v) = (f , v) ∀v ∈ H01 (Ω ).

(27)

The difference of (27) and (11) gives ah (u − uh , v) = (f , v − γ v) ∀v ∈ H01 (Ω ) ∩ Vh .

(28)

It follows from (23), the error equation in (28) with v = Jh α , and the continuity of Jh α ∈ 1/2

∥k

∇α∥

2 0 ,Ω

H01

(Ω ) ∩ Vh that

= (k∇h eh , ∇h α) = (k∇h eh , ∇h (α − Jh α))    We [γ Jh α] · {k∇ eh } − {k∇ Jh α} · [γ eh ] − δ [γ eh ] · [γ Jh α] ds + (f , Jh α − γ Jh α) + he

e

e∈Eh

= (k∇h eh , ∇h (α − Jh α)) −

 e∈Eh

{k∇ Jh α} · [γ eh ]ds + (f , Jh α − γ Jh α).

(29)

e

For e ∈ Eh0 , the definition of weighted average, Cauchy–Schwartz inequality, (14) and (7) and (19) give



{k∇ Jh α} · [γ eh ]ds = e



ω1 (k1 ∇ Jh α|K1 )[γ eh ] + ω2 (k2 ∇ Jh α|K2 )[γ eh ]ds e

 ≤

1/2 1/2

ω1 he k1 1/2 We

1/2 1/2

1/2

∥k1 ∇ Jh α|K1 ∥0,e +

ω2 he k2 1/2 We

1/2

∥k2 ∇ Jh α|K2 ∥0,e



1/2

We

1/2 he

∥[γ uh ]∥0,e

 W 1/2  e 1/2 1/2 ≤ C ∥k1 ∇ Jh α∥0,K1 + ∥k2 ∇ Jh α∥0,K2 ∥[γ uh ]∥0,e . 1/2 he

Similarly, we have

  e∈Eh ∩∂ Ω

{k∇ Jh α}[γ eh ]ds ≤

 

e

e∈Eh ∩∂ Ω

We he

∥[γ uh ]∥20,e

1/2

∥k1/2 ∇α∥0,Ω .

(30)

Thus, by summing over all e ∈ Eh and P1, we get

 e∈Eh

{k∇ Jh α}[γ eh ]ds ≤ e

 e∈Eh

We he

∥[uh ]∥

2 0,e

1/2

∥k1/2 ∇α∥0,Ω .

(31)

536

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

Furthermore, the property of γ and (19) gives

(f , Jh α − γ Jh α) ≤ C

h2K ∥f ∥20,K



kK

K

≤C

1/2 

∥Jh α − γ Jh α∥0,K 2

kK

1/2 

hK

K

 2 1/2  hK ∥f ∥20,K K

kK

1/2

∥k

∇α∥

2 0,K

1/2 

.

(32)

K

Denote eJ = α − Jh α . By integrating by parts, Cauchy–Schwartz inequality, (19) and (21), we have

(k∇h eh , ∇h (α − Jh α)) = −



(∇ · (k∇h eh ), eJ ) +

K ∈Th



=

K ∈Th

(f , eJ ) −

≤C

hK ∥f ∥0,K 1/2



∂K

(k∇ eh · n)eJ ds

[k∇ uh ]{eJ }ds e

 1/2

∥k

kK

K ∈Th

≤C

 e∈Eh

K ∈Th





∇α∥0,∆K +

 e∈Eh0

he We

∥[k∇ uh ]∥0,e ∥k1/2 ∇α∥0,∆e



1/2 (ηR2f ,K + ηJ2σ ,K ) ∥k1/2 ∇α∥0,Ω .

K ∈Th

Thus,

∥k1/2 ∇α∥0,Ω ≤ C ηR . Since u ∈

H01

( Ω ) , I h uh ∈

H01

(33)

(Ω ) on ∂ Ω , and β ∈ H (Ω ), then integrating by parts gives 1

(∇(u − Ih uh ), ∇ β) = 0. ⊥

This, together with the Cauchy–Schwarz inequality and (18), implies

∥k−1/2 ∇ ⊥ β∥20,Ω = (∇h eh , ∇ ⊥ β) =



(∇(u − uh ), ∇ ⊥ β)K

K ∈Th

=



(∇(Ih uh − uh ), ∇ ⊥ β)K

K ∈Th

≤



∥k1/2 ∇(Ih uh − uh )∥0,K ∥k−1/2 ∇ ⊥ β∥0,K

K ∈Th

≤



1/2

∥k

∇(Ih uh − uh )∥

2 0,K

1/2

∥k−1/2 ∇ ⊥ β∥0,Ω

K ∈Th

≤ C ∥k−1/2 ∇ ⊥ β∥0,Ω



ηJ2u ,K

1/2

.

K ∈Th

Thus −1/2

∥k

∇ β∥0,Ω ≤ C ⊥



η

2 Ju ,K

1/2

.

(34)

K ∈Th

Combining (25), (33) and (34) will get the conclusion in (26).



Lemma 4.2. There exists a positive constant C independent of the mesh-size and the ratio kmax /kmin such that



1/2

ηRf ,K ≤ C ∥k

∇ eh ∥0,K

 + 1/2 ∥f − fK ∥0,K , hK

kK

∀ K ∈ Th

(35)

and



he We

∥[k∇ uh ]∥0,e ≤ C

  hK ∥k1/2 ∇ eh ∥0,K + 1/2 ∥f − fK ∥0,K , K ∈ωe

kK

∀e ∈ Eh0

where fK is the average of f over K and ωe is the union of all elements that share the common edge e.

(36)

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

537

Proof. For any element K ∈ Th , let bK be the standard cubic bubble function whose support is K . Then for any polynomial q with degree at most m, there exist positive constants cm and Cm , such that cm ∥q∥

2 0 ,K

 ≤ K

q2 bK dx ≤ ∥q∥20,K ,

1 ∥∇(qbK )∥0,K ≤ Cm h− K ∥q∥0,K .

(37)

Integrating by parts leads to

(k∇h e, ∇h (fK bK ))K = (f − fK , fK bK )K + (fK , fK bK )K . It follows from the Cauchy–Schwarz inequality, the above equality, and (37) that C ∥fK ∥20,K ≤ (fK , fK bK )K = (k∇ eh , ∇(fK bK ))K − (f − fK , fK bK )

≤ ∥k1/2 ∇ eh ∥0,K ∥k1/2 fK ∇ bK ∥0,K + ∥f − fK ∥0,K ∥fK bK ∥0,K   1/2 k ≤ C K ∥k1/2 ∇ eh ∥0,K ∥fK ∥0,K + ∥f − fK ∥0,K ∥fK ∥0,K .

(38)

hK

Hence hK 1/2

kK



1/2

∥fK ∥0,K ≤ C ∥k

∇ eh ∥0,K +



hK 1/2

kK

∥f − fK ∥0,K .

(39)

Now, (35) is a direct consequence of the triangle inequality and the above inequality. For any edge e ∈ Eh0 , let be be the standard piecewise quadratic edge bubble function corresponding to the edge e whose support is ωe . Then for any polynomial q with degree at most m, there exist positive constants dm , Dm and Em , such that dm ∥q∥

2 0,e

 ≤ e

q2 be ds ≤ ∥q∥20,e ,

∥∇(qbe )∥0,ωe ≤ Dm he−1/2 ∥q∥0,e ,

∥qbe ∥0,ωe ≤ Em he1/2 ∥q∥0,e .

(40)

Similarly, integrating by parts leads to

(k∇h eh , ∇h ([k∇ uh ]e be ))ωe = (f , [k∇ uh ]e be )ωe + ∥[k∇ uh ]be1/2 ∥20,e . It follows from the above equality, the Cauchy–Schwarz inequality, and (40) that C ∥[k∇ uh ]∥20,e ≤ ∥[k∇ uh ]b1e /2 ∥0,e

≤ ∥k∇h eh ∥0,ωe ∥[k∇ uh ]∇ be ∥0,ωe + ∥f ∥0,ωe ∥[k∇ uh ]be ∥0,ωe  1/2   −1/2 1/2 ≤ Ch− ∥[k∇ uh ]∥0,e kK ∥k1/2 ∇ eh ∥0,K + hK kK ∥f ∥0,K , e K ∈ωe

which, together with (38) and the triangle inequality, implies



he We

∥[k∇ uh ]∥0,e ≤ C

  hK (∥fK ∥0,K + ∥f − fK ∥0,K ) ∥k1/2 ∇ eh ∥0,K + . 1/2 kK

K ∈ωe

Now (36) follows from the above inequality and (39).



Then the following efficiency bound is a straightforward consequence. Theorem 4.3. Let u ∈ H01 (Ω ) and uh ∈ Vh be the solutions of (1) and (10), respectively. There exists a positive constant C independent of the mesh size and the ratio kmax /kmin such that C ηR ≤ |||u − uh ||| +

 e∈Eh

We he

∥[uh ]∥20,e

1/2 +

 K ∈Th

h2K kK

∥f − fK ∥20,K

5. Flux recovery and the error estimator 5.1. Flux recovery Denote the local lowest order Raviart–Thomas ([34,35]) by RT0 (K ) = P0 (K )2 + xP0 (K )

1/2

.

(41)

538

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

where x = (x1 , x2 )t . Then the standard H (div; Ω ) conforming Raviart–Thomas space is defined by

Vh = {τ ∈ H (div; Ω ) : τ|K ∈ RT0 (K ), ∀K ∈ Th }. We recover the flux as follows: finding σ h ∈ Vh such that

(k−1 σ h , τ) = −(∇h uh , τ),

∀τ ∈ Vh .

(42)

Let ψ e be the nodal basis function on the edge e ∈ Eh . Define the explicit approximation σˆ h in Vh = span{ψ e : e ∈ Eh } by



σˆ h =

σe ψe ,

(43)

e∈Eh

where σe is the normal component of σˆh on the edge e ∈ Eh defined as a weighted average of the normal components of the approximated flux −k∇h uh , i.e.,

σe := {−k∇h uh · ne }ω

(44)

with weights ω1,i and ω2,i (i = 1, 2, 3) defined in (4). Hereafter we choose i = 2 or i = 3. Note that these weights satisfy the following inequality

(ω1,i )2 (ω2,i )2 max , 

k2



k1

1

≤

k1 + k2

=

2 We,1

,

i = 2, 3.

(45)

5.2. Recovery-based error estimators For any element K ∈ Th , define the following local a posteriori error estimator based on the solutions of (42) by

ησ2 ,K = ∥k1/2 σ h + k1/2 ∇ uh ∥0,K .

(46)

Let

ησ2 =



ησ2 ,K ,

ηJ2u =



ηJ2u ,K .

K ∈Th

K ∈Th

Now, we define the recovery-based a posteriori error estimators as follows:

 1/2 ξ = ησ2 + ηJ2u .

(47)

Theorem 5.1. Let u ∈ H01 (Ω ) and uh ∈ Vh be the solutions of (1) and (10), respectively. There exists a positive constant C independent of the mesh size and the ratio kmax /kmin such that

|||u − uh ||| ≤ C (ξ + Hf ),

(48)

with Hf =



1 2 2 k− K h K ∥ f ∥ 0 ,K

1/2

.

(49)

K ∈Th

Proof. To establish the validity of inequality (48), by (34), (25), P1 and the definition of the norm, it suffices to show that

∥k1/2 ∇α∥0,Ω ≤ C (ξ + Hf ).

(50)

Let eJ = α − Jh α . Using integrating by parts, Cauchy–Schwartz inequality and the fact ∇h · ∇h uh = 0, we get

(k∇h eh , ∇h eJ ) = (k∇ u + σ h , ∇ eJ ) − (σ h + k∇h uh , ∇h eJ ) ≤ (f − ∇ · σ h , eJ ) + C ησ ∥k1/2 ∇α∥0,Ω = (f , eJ ) − (∇h · (σ h + k∇h uh ), eJ ) + C ησ ∥k1/2 ∇α∥0,Ω ≤ C (ησ + Hf )∥k1/2 ∇α∥0,Ω . Combining with (29), (31) and (32), yields the inequality in (50). Thus completes the proof of the theorem.



L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

539

Lemma 5.2. For any element K ∈ Th , let ωK be the union of elements sharing a common edge with K . There exists a positive constant C independent of the mesh size and the ratio kmax /kmin such that

 1/2  h2 K ησ ≤ C ∥k1/2 ∇α∥20,ωK + ∥f − fK ∥20,K , K ∈Th

kK

∀K ∈ Th .

(51)

Proof. For any edge e ∈ ∂ K , without loss of generality, let ne be the outward unit vector normal to ∂ K , and denote by Ke the adjacent element with the common edge e. And set the weights by ω1,i and ω2,i (i = 2, 3) corresponding to element K and Ke respectively. Since τˆ h = −k∇h uh is piecewise constant, τˆ h |K may be represented in terms of the nodal basis function of Vh , {ψ e }e∈∂ K , as follows:

τˆ h |K =



τe,K ψe .

e∈∂ K

For any x ∈ K ,

σ h − τˆ h =



(σe − τe,K )ψe =

e∈∂ K

=−



(w1,i − 1)(τe,K − τe,Ke )ψe

e∈∂ K



w2,i [τˆ h ]ne · ψe .

e∈∂ K

Now, it follows from the triangle inequality, the fact that

 K

|ψe |2 dx ≤ C , and (45) that

1 ˆ h ∥20,K ησ2 ,K = ∥k−1/2 (σ h − τˆ h )∥20,K ≤ Ck− K ∥σ h − τ      he −1 2 2 2 ∥[k∇ uh ]∥0,e , ≤C kK (w2,i ) he [τˆ h ] ds ≤ C e

e∈∂ K

e∈∂ K

We,1

which together with (36) proves (51). This completes the proof of the Lemma.



The following Theorem is a straightforward consequence from Lemma 5.2. Theorem 5.3. Let u ∈ H01 (Ω ) and uh ∈ Vh be the solutions of (1) and (10), respectively. There exists a positive constant C independent of the mesh size and the ratio kmax /kmin such that

1/2   2 1/2    We hK ∥[uh ]∥20,e + ∥f − fK ∥20,K . ξ ≤ C |||u − uh ||| + e∈Eh

he

K ∈Th

kK

(52)

Our analysis can be easily extend to the interface problems with non-homogeneous Dirichlet boundary conditions. Indeed, this is the kind of boundary condition used in the following numerical experiments. 6. Numerical experiments In this section, we shall report several numerical results for the interface problems with continuous and discontinuous coefficients respectively. In all the examples, we choose δ = 1 for testing. Example 1 (L-shaped). Set k(x) ≡ 1 and Ω = [−1, 1]2 /[0, 1] × [−1, 0]. We choose an appropriate function uD such that u( r , θ ) = r

2/3

 sin

2 3

θ



in polar coordinates which satisfies:



−∇ · k∇ u = 0, u|∂ Ω = uD .

in Ω

(53)

This is a widely used problem with corner singularity. We want to show that our method is valid for this problem. We start with a coarse triangulation T0 , and then adopt the Dörfler’s bulk marking strategy [36] with θ = 0.2 to generate a sequence of refined mesh. The stopping criteria is chosen as

|||u − uh ||| ≤ tol. ∥k1/2 ∇ u∥0,Ω

540

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

a

b

Fig. 2. Error on adaptively generated meshes in Example 1 by ηR (a) and ξ (b).

a

b

Fig. 3. Final mesh for estimator ηR (a) and ξ (b).

Figs. 2 and 3 show that the norm of the error and the final refined mesh with tol = 0.01. Fig. 2 shows the results for residual-based estimator and recovery-based estimator, respectively. Since the exact solution u ∈ H 1+2/3−ϵ , ∀ϵ > 0, for the DG finite volume methods with linear polynomials, on uniformly bisecting mesh, we expect a convergence rate of order 23

for |||u − uh |||. Thus, theoretically, it is expected that the slope of the log(number of nodes)–log(relative error) is − 13 . It is observed, numerically the slope is very close to − 31 . Furthermore, for adaptive mesh, numerically we get a reduction with

slope − 12 of the log(Number of Nodes)–log(relative error) for all estimators, which indicates optimal rate. Both of the refined mesh shown in Fig. 3 is concentrated around the singularity located at the origin, which is the singularity. In addition, the value of ηR /|||u||| and ξ /|||u||| efficiently and accurately approximate the relative error of |||u − uh |||/|||u|||. Next we will consider another numerical example of the interface problem from [31], which is widely used by many authors. Example 2. For Ω = [−1, 1]2 , let u(r , θ ) = r γ µ(θ ) be the solution of (53), and set

 cos((π /2 − σ )γ ) · cos((θ − π /2 + ρ)γ ),   cos(ργ ) · cos((θ − π + σ )γ ), µ(θ ) = cos(σ γ ) · cos((θ − π − ρ)γ ),  cos((π /2 − ρ)γ ) · cos((θ − 3π /2 − σ )γ ),

if 0 ≤ 1 ≤ π /2 if π /2 ≤ 1 ≤ π if π ≤ 1 ≤ 3π /2 if 3π /2 ≤ 1 ≤ 2π

in the polar coordinates at the origin (r and θ denote the radial coordinate and the angular coordinate). Here the numbers γ , ρ, and σ satisfying the nonlinear relations

R := a /a = − tan((π /2 − σ )γ ) · cot(ργ ), 1 2    1/R = − tan(ργ ) · cot(σ ρ),   R = − tan(ργ ) · cot((π /2 − ρ)γ ), 0 < γ < 2    max{0, π γ − π } < 2γ ρ < min{π γ , π} max{0, π − π γ } < −2γ σ < min{π , 2π − π γ }.

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

541

Fig. 4. Refined mesh and convergence rate of ηR in Example 2 with tol = 0.1 for ω1,1 and ω2,1 .

Fig. 5. Refined mesh and convergence rate of ηR in Example 2 with tol = 0.1 for ω1,2 and ω2,2 .

Fig. 6. Refined mesh and convergence rate of ηR in Example 2 with tol = 0.1 for ω1,3 and ω2,3 .

Fig. 7. Refined mesh and convergence rate of ξ in Example 2 with tol = 0.1 for ω1,1 and ω2,1 .

We have tested the example in some worst-case scenarios, for example, γ = 0.1, which produces a very singular solution u that is barely in H 1.1 . Here R = a1 /a2 ≈ 161.4476, ρ = π /4, σ ≈ −14.9226, and k = a1 = R, in (0, 1)2 ∪ (−1, 0)2 and k = a2 = 1 in Ω \ ((0, 1)2 ∪ (−1, 0)2 ). Numerical results for tol = 0.1 are shown in Figs. 4–9. Figs. 4–6 show the final refinement and convergence results for residual-based a posteriori error estimator for three choices of weights ω1 and ω2 . Figs. 7–9 present results for recoverybased a posteriori error estimator for three choices of weights ω1 and ω2 . The results show that the mesh refinement based on all the estimators concentrates at the origin. Thus, these estimators could capture the interface correctly. Since the exact

542

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

Fig. 8. Refined mesh and convergence rate of ξ in Example 2 with tol = 0.1 for ω1,2 and ω2,2 .

Fig. 9. Refined mesh and convergence rate of ξ in Example 2 with tol = 0.1 for ω1,3 and ω2,3 .

Table 1 Error profiles for test problem 2 on adaptive meshes.

ηR i=1 ηR i=2 ηR i=3 ξi=1 ξi=2 ξi=3

l

N

||eh ||

||eh || = O(N k )

138 94 111 138 93 108

1699 4155 2364 1699 4155 2364

0.0565 0.0544 0.0558 0.0565 0.0544 0.0558

k k k k k k

= −0.5024 = −0.4290 = −0.4684 = −0.5024 = −0.4284 = −0.4689

η

η = O(N k )

0.1111 0.0434 0.0627 0.1111 0.0434 0.0627

k k k k k k

= −0.3142 = −0.1749 = −0.2200 = −0.3142 = −0.1750 = −0.2285

rel-err

eff-index

0.0999 0.0963 0.0988 0.0999 0.0963 0.0988

1.9683 0.7985 1.1229 1.9683 0.7985 1.1229

Table 2 Error profiles for test problem 3 on uniform meshes. h

||eh ||

ηR

ξ

1/4 1/8 1/16 1/32 1/64 1/128 1/256 Conv. Rate

6.7321e−01 3.3476e−01 1.6511e−01 8.1899e−02 4.0786e−02 2.0354e−02 1.0168e−02 1.0087

8.8667e−01 3.6739e−01 1.6580e−01 7.9376e−02 3.8985e−02 1.9343e−02 9.6372e−03 1.0769

4.4924e−01 2.6430e−01 1.4196e−01 7.3912e−02 3.8068e−02 1.9522e−02 9.9888e−03 9.2464e−01

solution u has regularity H 1.1 , theoretically, we expect a reduction rate of relative error in discrete H 1 -norm as 0.1. Thus, theoretically, it is expected that the slope of the log(number of nodes)–log(relative error) is −0.05. It is observed, numerically the slope is very close to −0.05. Furthermore, for adaptive mesh, numerically we get a reduction with slope − 12 of the log(number of nodes)–log(relative error) for all estimators, which indicates optimal rate. Additionally, the three choices of weights show different behaviors. We note that harmonic weight is the best option for testing, which shows consistent decay of error with respect to the number of nodes. Table 1 compares these estimators for tol = 0.1. Here l is the iteration steps; N is the number of nodes; ref-err= ||uh|| ; . It is observed that under the same stopping criteria: i = 2 need fewest iteration steps but most number err-index= Estimator ||eh || of nodes; i = 1 need maximum iteration steps but the minimum number of nodes; effective index of i = 3 gives the best approximation to 1. Indeed, all of these show that our method is effective and efficient. ||e ||

L. Mu, R. Jari / Journal of Computational and Applied Mathematics 255 (2014) 529–543

543

Fig. 10. The convergence test of Example 3.

Example 3. Let Ω = (0, 1)2 . The last example assumes that the problem has an exact solution u = sin(π x) sin(π y), and k = (x + 1)(y + 1). The homogeneous Dirichlet boundary condition is given for testing. In the experiment, we solve the test problem on uniform meshes. The sequences of meshes are generated by: (1) dividing Ω into n × n sub-rectangles; (2) dividing each sub-rectangle into two triangles by connecting the negative slope line. Denote h = 1/n. The error profile, the recovery-based estimator ξ , the residual-based estimator ηR , and rate of convergence are shown in Table 2 for i = 2 and plotted in Fig. 10. Our numerical results show that both estimators can give good estimate for the error in ||| · |||-norm. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36]

I. Babuška, W.C. Rheinboldt, Error estimates for adaptive finite element computations, SIAM J. Numer. Anal. 18 (1978) 736–754. I. Babuška, W.C. Rheinboldt, A posterior error estimates for the finite element methods, Int. J. Numer. Methods Eng. 12 (1978) 1597–1615. C. Bernardi, R. Verfürth, Adaptive finite element methods for elliptic equations with non-smooth coefficients, Numer. Math. 85 (4) (2000) 579–608. R. Verfürth, A Review of a Posteriori Error Estimation and Adaptive Mesh-refinement Techniques, Wiley-Teubner, Chichester, 1996. O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, 6th edition, Elsevier, 2005. M. Ainsworth, J.T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Wiley-Interscience, New York, 2000. I. Babuška, J. Whiteman, T. Strouboulis, Finite Elements: An Introduction to the Method and Error Estimation, Oxford University Press, 2011. P. Ladevèze, J.P. Pelle, Mastering Calculations in Linear and Nonlinear Mechanics, Springer, 2004. I. Babuška, T. Strouboulis, C.S. Upadhyay, S.K. Gangaraj, A model study of element residual estimators for linear elliptic problems: the quality of the estimators in the interior of meshes of triangles and quadrilaterals, Comput. Struct. 57 (6) (1995) 1009–1028. E. Dari, R.G. Durán, V. Vampa, A posteriori error estimators for nonconforming finite element methods, Math. Model. Numer. Anal. (RAIRO) 30 (4) (1996) 385–400. I. Babuška, T. Strouboulis, The Finite Element Method and its Reliability, Oxford University Press, 2001. R. Li, W.B. Liu, N.N. Yan, A posteriori error estimates of recovery type for distributed convex optimal control problems, J. Sci. Comput. 33 (2007) 155–182. N.N. Yan, A posteriori error estimators of gradient recovery type for elliptic obstacle problems, Adv. Comput. Math. 15 (1–4) (2001) 333–362. T. Tang, J.C. Xu, Adaptive Computations: Theory and Algorithms-Mathematics Monograph Series 6, Science publisher, 2007. M. Ainsworth, A. Craig, A posteriori error estimators in the finite element method, Numer. Math. 60 (1) (1991) 429–463. M. Ainsworth, J.T. Oden, A posteriori error estimation in finite element analysis, Comput. Methods Appl. Mech. Engrg. 142 (1997) 1–88. R.E. Bank, J.C. Xu, Asymptotically exact a posteriori error estimators, part I: grids with superconvergence, SIAM J. Numer. Anal. 41 (6) (2003) 2294–2312. N.N. Yan, A. Zhou, Gradient recovery type a posteriori error estimates for finite element approximations on irregular meshes, Comput. Methods Appl. Mech. Engrg. 190 (32–33) (2001) 4289–4299. Z. Zhang, J. Zhu, Analysis of the superconvergent patch recovery technique and a posteriori error estimator in the finite element method (I), Comput. Methods Appl. Mech. Engrg. 123 (1–4) (1995) 173–187. Z. Zhang, J. Zhu, Analysis of the superconvergent patch recovery technique and a posteriori error estimator in the finite element method (II), Comput. Methods Appl. Mech. Engrg. 163 (1–4) (1998) 159–170. O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery (SPR) and adaptive finite element refinement, Comput. Methods Appl. Mech. Engrg. 101 (1–3) (1992) 207–224. P. Katragadda, I.R. Grosset, A posteriori error estimation and adaptive mesh refinement for combined thermal-stress finite element analysis, Comput. Struct. 59 (6) (1996) 1149–1163. A. Bergam, Z. Mghazli, R. Verfurth, A posteriori estimators for the finite volume discretization of an elliptic problem, Numer. Math. 95 (2003) 599–624. J. Liu, L. Mu, X. Ye, An adaptive discontinuous finite volume method for elliptic problems, J. Comput. Appl. Math. 235 (18) (2011) 5422–5431. L. Mu, R. Jari, L. Conner, A recovery-based error estimator for finite volume metods of interface problems: conforming linear elements, Numer. Methods Partial Differential Equations 29 (2013) 131–143. D. Pietro, M. Vohralk, C. Widmer, An a posteriori error estimator for a finite volume discretization of the two-phase flow, in: Finite Volumes for Complex Applications VI Problems & Perspectives, vol. 4, 2011, pp. 341–349. J.C. Xu, Q.S. Zou, Analysis of linear and quadratic simplicial finite volume methods for elliptic equations, Numer. Math. 111 (3) (2009) 469–492. X. Ye, A posteriori error estimation for finite volume methods of the second order elliptic problem, Numer. Methods Partial Differential Equations 27 (2011) http://dx.doi.org/10.1002/num.20575. X. Ye, A new discontinuous finite volume method for elliptic problems, SIAM J. Numer. Anal. 42 (3) (2004) 1062–1072. Z.Q. Cai, X. Ye, S. Zhang, Discontinuous Galerkin finite element methods for interface problems: a priori and a posteriori error estimations, SIAM 49 (2011) 1761–1787. R.B. Kellogg, On the Poisson equation with intersecting interfaces, Appl. Anal. 4 (2) (1974) 101–129. D.N. Arnold, An interior penalty finite-element method with discontinuous elements, SIAM J. Numer. Anal. 19 (4) (1982) 742–760. M. Ainsworth, A posteriori error estimation for discontinuous Galerkin finite element approximation, SIAM J. Numer. Anal. 45 (2007) 1777–1798. P.A. Raviart, J.M. Thomas, A mixed finite element method for second order elliptic problems, math, aspects of the FEM, Lecture Notes in Math. 606 (1977) 292–315. F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer Verlag, New York, 1991. W. Döfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal. 33 (1996) 1106–1124.