Error indicator for the incompressible Darcy flow problems using Enhanced Velocity Mixed Finite Element Method.

# Error indicator for the incompressible Darcy flow problems using Enhanced Velocity Mixed Finite Element Method.

Yerlan Amanbek Gurpreet Singh Oden Institute for Computational Engineering and Sciences, University of Texas at Austin Gergina Pencheva Oden Institute for Computational Engineering and Sciences, University of Texas at Austin Mary F. Wheeler Oden Institute for Computational Engineering and Sciences, University of Texas at Austin
July 2, 2019
###### Abstract

In the flow and transport numerical simulation, mesh adaptivity strategy is important in reducing the usage of CPU time and memory. The refinement based on the pressure error estimator is commonly-used approach without considering the flux error which plays important role in coupling flow and transport systems. We derive a posteriori error estimators for Enhanced Velocity Mixed Finite Element Method (EVMFEM) in the incompressible Darcy flow. We show numerically difference of the explicit residual based error estimator and implicit error estimators, where Arbogast and Chen post-processing procedure from [1] for pressure was used to improve estimators. A residual-based error estimator provides a better indicator for pressure error. Proposed estimators are good indicators in finding of the large error element. Numerical tests confirm theoretical results. We show the advantage of pressure postprocessing on the detecting of velocity error. To the authors’ best knowledge, a posteriori error analysis of EVMFEM has been scarcely investigated from the theoretical and numerical point of view.

Keywords.  a posteriori error analysis, enhanced velocity mixed finite element method, error estimates, adaptive mesh refinement.

## 1 Introduction

In subsurface problems, a computational saving of numerical simulation can be achieved by reducing degrees of freedom in the linear system. Appropriate reduction of degrees of freedom of the problem by controlling error of approximation is often handled via a posteriori error analysis. Such a special assessment of errors provides basis for mesh refinement or unrefinement strategy. A posteriori error estimator and Error Indicator of elements is essential to determine the large error element which can be considered for refinement to achieve the accurate and efficient subsurface simulations. By knowing the assessment of error provided by estimators, one can control of discretization error and achieve the anticipated quality of the numerical solution.

A posteriori error estimators for finite element method for elliptic boundary value problems started by the Babuska and Rheinbodt work in [2]. The main objective of a posteriori error estimation is to obtain the estimator that is close to the error in specific (e.g. energy) norm on each element [3]. The literature on a posteriori error estimates and adaptivity for mixed finite element approximation has highlighted several forms of estimates for mesh refinement strategy. To derive estimates for many finite element approximation, many researchers have proposed various methods of optimal a posteriori error estimate in [3, 4, 5, 6]. In particularly, conforming mixed finite element method were explored in [7, 8, 9, 10] as well as Discontinuous Galerkin (DG) method employing explicit error estimate in [11] and implicit error estimate in [12], goal-oriented Discontinuous Petrov-Galerkin(DPG) [13] with applications [14], for nonconforming FEM in [15, 16, 17] such as Multiscale Mortar Mixed FEM in [18, 19, 20, 21, 22]. However, to the best of our knowledge, a posteriori error analysis has not been conducted for Enhanced Velocity Mixed FEM, which is practical method in adaptive setting for many subsurface applications [23, 24, 25, 26, 27].

For large domain with heterogeneous permeability , which varies over scales from mm to hundreds kms, it is computationally intensive task. Direct discretization in fine scale of entire domain would involve high resolution permeability of that resulting a large and coupled system of equation. Such system solution would usually be CPU and Memory demanding.

Recent development for flow and transport in heterogeneous porous media, adaptive numerical homogenization using Enhanced Velocity Mixed Finite Element Method (EVMFEM) in [23, 28, 26], motivates us to study a posteriori error estimates in the adaptive mesh refinement strategy. The use of adaptive numerical homogenization captures fine scale features in efficient way for heterogeneous porous media. This method can be successfully used for a number of subsurface engineering applications to achieve efficient and accurate numerical solution.

The main purpose of this paper is to provide a posteriori error estimates for single phase flow using EVMFEM. We show theoretical derivations of estimators and then numerical results for selected indicators to confirm theoretical upper bounds. Velocity error is taken into account to increase accuracy of velocity that are important in transport equations. We consider nested version of mesh discretization at the interface for numerical experiments. First, we derive the explicit residual-based a posteriori error estimates for EVMFEM with saturation assumption. Second, we show the implicit error estimates without saturation assumption using a suitable post-processing of the finite element pressure approximation that is better indicator for mesh refinement. We show estimates in norm. For simplicity, the problem is considered with Dirichlet boundary conditions, but the result can be generalized.

The remainder of this paper is organized as follows. Section 2 outlines formulation of EVMFEM and preliminaries. Our a posteriori error estimates are presented in Section 3. In this section, we start by going through briefly preliminaries and useful inequalities, formulation and projections. Second, residual-based explicit error estimators were derived for the pressure and velocity errors. Then, residual-based estimators with smoothing were derived to improve the indicators of velocity error by using post-processed pressure such as the Arbogast and Chen postprocessing [1]. Section 4 shows computational results. The conclusion is reported in Section 5.

## 2 Model formulation

We start by giving the model formulation for the incompressible single-phase flow. For the convenience of reader we repeat the relevant material of domain decomposition method, discrete formulation with Enhanced Velocity from [25]. EVMFEM is a mass conservative and an efficient domain decomposition method which deals with non-matching grids or multiblock grids [25]. This method is extensively used in many complex multicomponenet, multicomponent, multiphase flow and transport processes in porous media [24]. EVMFEM approach is strongly mass conservative at the interfaces and impose strong continuity of fluxes between the subdomains. Earlier implementations [25, 24] employed a solution approach where only the coarse and fine domain contributions to the stiffness-matrix (or Jacobian matrix) were taken, neglecting interface contributions. The load vector (or residuals); however, contains contributions from both the coarse and fine subdomains as well as the interface. This resulted in an increase in the number of non-linear iterations to achieve convergence, for a given tolerance, even for a linear flow and transport problem. In this work, we use a fully coupled variant of the original EVMFEM approach wherein the interface terms are properly accounted for in the stiffness-matrix construction resulting in reduced non-linear iterations (one for a linear system).

### 2.1 Governing Equations of the Incompressible Flow

For convenience of analysis, we consider the incompressible single phase flow model for pressure and the Darcy velocity u:

 u =−K∇pinΩ, (1) ∇⋅u =finΩ, (2) p =gon∂Ω (3)

where or ) is multiblock domain, and is a symmetric, uniformly positive definite tensor representing the permeability divided by the viscosity with components, for some

 kminξTξ≤ξTK(x)ξ≤kmaxξTξ∀x∈Ω∀ξ∈Rd. (4)

Let be divided into a series of small subdomains. For simplicity, the Dirichlet boundary condition is considered as zero, i.e. . To formulate in mixed variational form, Sobolev spaces are exploited and the following space is defined for flux in as usual to be and equipped with the norm and for the pressure the space is and the corresponding norm .

We utilize standard notations. For subdomain , the inner product (or duality pairing) and norm are denoted by and , respectively, for scalar and vector valued functions. Let be the standard Sobolev space of -differentiable functions in . Let be norm of or , where and are omitted in case of and respectively, in other cases they are specified. We write for the or inner product, and for duality pairing on boundaries and interfaces, where the pairing may be between two functions in or between elements of and , in either order.

Next, a weak variational form of the fluid flow problem is to find a pair ,

 (K−1u,v)−(p,∇⋅v) =−⟨g,v⋅ν⟩∂Ω ∀v∈V (5) (∇⋅u,w) =(f,w) ∀w∈W (6)

where is the outward unit normal to .

#### Discrete formulation

We consider

 Ω=(Nb⋃i=1¯Ωi)o,Γi,j=∂Ωi⋂∂Ωj,Γ=(Nb⋃i,j=1¯Γi,j)o,Γi=Ωi⋂Γ=∂Ωj∖∂Ω.

This implies that the domain is divided into subdomains, the interface between and subdomains(), the interior subdomain interface for subdomain and union of all such interfaces, respectively.

Let be a conforming, quasi-uniform and rectangular partition of , , with maximal element diameter . We then set and denote the maximal element diameter in ; note that can be nonmatching as neighboring meshes and need not match on . We assume that all mesh families are shape-regular.

We narrow our work to the regularly used Raviart-Thomas spaces of lowest order on rectangles for and bricks for . The spaces are defined for any element by the following spaces:

 Vh(T)={v=(v1,v2)orv=(v1,v2,v3):vl=αl+βlxl:αl,βl∈R;l=1,..d}and Wh(T)={w=constant}

In fact, a vector function in can be determined uniquely by its normal components at midpoints of edges (in 2D) or face (in 3D) of . The degrees of freedom of were created by these normal components. The degree of freedom for a pressure function is at center of and piecewise constant inside of . The pressure finite element approximation space on is taken to be as

 Wh(Ω)={w∈L2(Ω):w∣∣∣E∈Wh(T),∀T∈Th}

We first construct a velocity finite element approximation space on , which is different from the velocity space of the Multiscale Mortar Mixed FEM. Let us formulate space on each subdomain for partition

 Vh,i={v∈H(div;Ωi):v∣∣∣T∈Vh(T),∀T∈Th,i}i∈{1,...n}

and then

 Vh=n⨁i=1Vh,i.

Although the normal components of vectors in are continuous between elements within each subdomains, the reader may see is not a subspace of , because the normal components of the velocity vector may not match on subdomain interface . To solve this issue, many researchers have proposed various methods such as Multiscale Mortar Mixed FEM [29], Enhanced Velocity Mixed FEM [25], etc. In Mortar Multiscale Mixed FEM , the mortar finite element space on coarse grid was introduced to connect subdomains together using Lagrange multipliers to enforce weak continuity for flux across subdomains. On the other hand, the Enhanced Velocity Mixed FEM modifies the degree of freedom on to finer grids, which impose the strong flux continuity between subdomains. Let us define as the intersection of the traces of and , and let . We require that and need to align with the coordinate axes. Fluxes are constructed to match on each element . We consider any element that shares at least one edge with the interface , i.e., , where and . Then newly defined interface grid introduces a partition of the edge of . This partition may be extended into the element as shown in Figure 2.

This new partitioning helps to construct fine-scale fluxes that is in . So we represent a basis function in the space () for given with the following way:

 vTk⋅ν={1,onek0,otheredges

i.e. a normal component equal to one on and zero on all other edges(faces) of . Let be span of all such basis functions defined on all sub-elements induced the interface discretization . Thus, the enhanced velocity space is taken to be as

 V∗h=n⨁i=1V0h,i⨁VΓh∩H(div;Ω).

where is the subspace of . The finer grid flux allows to velocity approximation on the interface and then form the conforming velocity space. Some difficulties arise, however, in analysis of method and implementation of robust linear solver for such modification of velocity space at all elements, which are adjacent to the interface . We now formulate the discrete variational form of equations as: Find and such that

 (K−1uh,v) =(ph,∇⋅v)−⟨g,v⋅ν⟩∂Ω ∀v∈V∗h (7) (∇⋅uh,w) =(f,w) ∀w∈Wh (8)

## 3 Methodology of Error Estimate

In this section, we derive error estimators for enhanced velocity mixed finite element discretization of elliptic problems. First, we briefly go through preliminaries and useful inequalities, formulation and projections. Second, residual based explicit error estimators were derived for pressure and velocity. Third, residual based estimators with smoothing were derived to improve the indicators of velocity error by using post-processed pressure such as the Arbogast and Chen postprocessing.

### Representation of error

We define

 eu=u−uh and ep=p−ph.

We introduce the bilinear form defined as

 A(u,p;v,w)=(K−1u,v)−(p,∇⋅v)+μ(∇⋅u,w)

where or . We denote as when , which is a symmetric bilinear form and as when , which is a nonsymmetric, but coercive, since . Linear functional is defined as

 L(v,w)=μ(f,w)−⟨g,v⋅ν⟩

where or . Note that the solution does not depend on the choice of . Therefore, weak variational form of (5)-(6) imply that satisfy

 A(u,p;v,w)=L(v,w)(v,w)∈V×W

and the discrete variational formulation in (7)- (8) implies that

 A(uh,ph;v,w)=L(v,w)(v,w)∈V∗h×Wh (9)

Then using above notations we obtain the residual equation

 A(eu,ep;v,w)=L(v,w)−A(uh,ph;v,w)(v,w)∈V×W (10)

hold true for each pair . Thus, (9) and (10) give the orthogonality condition

 A(eu,ep;v,w)=0(v,w)∈V∗h×Wh (11)

### Preliminaries and useful inequalities

We assume the model problem is -regular, in other words, there exists a positive C that depends on and such that

 (12)

We refer to the reader to [30] for sufficient conditions for -regularity.

We present below some of the approximation properties of the finite element spaces

 (w−^w,wh)=0∀wh∈Wh

Furthermore, the following important approximation properties hold true. For all , , and smooth enough and

 ∥v−Πv∥T ≤ChT∥v∥1,T (13) ∥w−^w∥ ≤CPhT∥w∥1,T (14)

where , if is convex. Above inequalities are standard -projection approximation results and can be found in [31]. In the analysis below we will use of trace inequalities

 ∀T∈Th,e∈∂T,∥ϕ∥e ≤C(h−1/2T∥ϕ∥+h1/2T∥∇ϕ∥)ϕ∈H1(T) (15) ∀T∈Th,e∈∂T,∥ϕ∥1/2,e ≤C∥ϕ∥1,Tϕ∈H1(T) (16) ≤Ch−1/2T∥v∥Tv∈V∗h (17)

and Young’s inequality

 2ab≤εa2+1εb2

We know that from original work [25] the projection operator was introduced and was utilized for a priori error analysis. For convenience of the reader, we repeat the relevant and brief definition. Thus, we denote by the projection operator that maps onto that defined locally for any element and any such that for all

 ⟨Π∗q⋅ν,1⟩e=⟨q⋅ν,1⟩e (18)

where is either any edge in 2D (or face in 3D) of not lying on or an edge in 2D (or face in 3D) of a sub-element, . Such projection is developed prior to conducting error analysis for a priori and a posteriori error estimates. As can be seen in Figure 2, has a common edge with the interface grid . According to divergence theorem, we have

 (∇⋅(Π∗q−q),w)=0∀w∈Wh (19)

We refer the reader to the original work [25] for more details. We want to scale it to local element. So by scaling argument we reach the following lemma

###### Lemma 1.

Let then C independent of such that

 ∥Π∗v−v∥T≤ChT∥v∥1,T (20)
###### Lemma 2.

Let then C independent of such that

 ∥(Π∗v−v)⋅ν∥∂T≤Ch1/2T∥v∥1,T (21)

### 3.1 Explicit Residual-based Error Estimators

In this section, we derive upper bounds on the local error. It is also called explicit estimators as they involve the input data and computed numerical solution without solving extra sub-problems. We are not interested in computing the constants in the error estimates and take the boundary condition as .

#### 3.1.1 Estimates for Pressure

###### Theorem 3.

There exists a constant independent of such that

 ∥ep∥2≤C{∑T∈Th(~ζP+~ζR)+~ζEV} (22)

where, for all

 ~ζP ~ζR,h =∥f−∇⋅uh∥2Th2T ~ζEV =∑e∈TΓh∥∥\llbracketph\rrbracket∥∥2ehT
###### Proof.

We consider a duality argument to derive bounds. Let be the solution of the auxiliary problem

 −∇⋅K∇φ =epinΩ, (23) φ =0on∂Ω. (24)

By the elliptic regularity assumption 12 implies that

 ∥φ∥2≤C∥ep∥

Let then

 As(v,φ;~v,~w)=(K−1v,~v)−(φ,∇⋅~v)−(∇⋅u,~w)=−(ep,~w)

Then

 ∥ep∥2 =−As(v,φ;eu,ep)=−As(eu,ep;v,φ)=−As(eu,ep;v−Π∗v,φ−^φ)= =−∑T∈Th{(K−1eu,v−Π∗v)T−(ep,∇⋅(v−Π∗v))T−(∇⋅eu,φ−^φ)T}= =−∑T∈Th{(K−1u,v−Π∗v)T−(K−1uh,v−Π∗v)T−(p,∇⋅(v−Π∗v))T +(ph,∇⋅v−Π∗v)T−(∇⋅u,φ−^φ)T+(∇⋅uh,φ−^φ)T}= =∑T∈Th{(K−1uh,v−Π∗v)T−(ph,∇⋅(v−Π∗v))T+(f,φ−^φ)T −(∇⋅uh,φ−^φ)T}

By Green’s formula,

 (ph,∇⋅(v−Π∗v))Ωi=−(∇ph,v−Π∗v)Ωi+⟨ph,(v−Π∗v)⋅ν⟩∂Ωi

we obtain

 ∥ep∥2 =∑T∈Th{(K−1uh+∇ph,v−Π∗v)T+(f−∇⋅uh,φ−^φ)T}−n∑i=1⟨ph,(v−Π∗v)⋅νi⟩∂Ωi

We use the Cauchy-Schwarz inequality and approximation properties

 ∥ep∥2 ≤∑T∈Th(∥∥K−1uh+∇ph∥∥T∥v−Π∗v∥T+∥f−∇⋅uh∥T∥φ−^φ∥T) −∑T∈Th⟨ph,(v−Π∗v)⋅ν⟩∂T ≤C∑T∈Th(∥∥K−1uh+∇ph∥∥ThT∥φ∥2,T+∥f−∇⋅uh∥ThT∥φ∥1,T) +C∑e∈TΓh,e∈∂Ti∪∂Tj∥∥\llbracketph\rrbracket∥∥eh12T∥v∥1,T ≤C∑T∈Th(∥∥K−1uh+∇ph∥∥ThT∥ep∥T+∥f−∇⋅uh∥ThT∥ep∥T) +C∑e∈TΓh,e∈∂Ti∪∂Tj∥∥\llbracketph\rrbracket∥∥eh12T∥v∥1,T ≤C∑T∈Th(∥∥K−1uh+∇ph∥∥ThT∥ep∥T+∥f−∇⋅uh∥ThT∥ep∥T) +C∑e∈TΓh,e∈∂Ti∪∂Tj∥∥\llbracketph\rrbracket∥∥eh12T∥ep∥T

#### 3.1.2 Estimates for Velocity with Saturation Assumption

In order to get bounds on velocity error (), we employ a saturation assumption. Let be be the finite element approximation spaces which corresponds to refinement of , where for . This implies that and . A priori error estimates from [25] allows us to employ the following saturation assumption.

### Saturation assumption

There exist constant , and for , such that

 ∥∥u−uhf∥∥ ≤β∥∥u−uh∥∥ (25) ∥∥p−phf∥∥ ≤α∥p−ph∥ (26)
###### Lemma 4.

Let be the exact flux defined by Eqns with . Let be arbitrary and for , . If there exist such that

 ∥∥u−uhf∥∥≤β∥u−uh∥, (27)

then,

 (28)

The proof of lemma 4 is straightforward by using the triangle inequality.

We write and that are the enhanced velocity mixed finite element solution of equations (7) - (8) and let

 ~eu=uhf−uh,~ep=phf−ph (29)

We have that which satisfy the residual equation

 A(~eu,~ep;~vh,~wh)=L(~vh,~wh)−A(uh,ph;~vh,~wh)∀(~vh,~wh)∈V∗hf×Whf (30)

and the orthogonality condition

 A(~eu,~ep;vh,wh)=0∀(vh,wh)∈V∗h×Wh (31)
###### Theorem 5.

Assume that the saturations assumptions (25) and (26) hold. Then there exists a constant independent of such that

 ∥eu∥2≤C⎛⎝∑T∈Th{ζP+ζR}+ζEV⎞⎠

where, for all

 ζP ζR,h =∥f−∇⋅uh∥2Th2T ζEV =∑e∈TΓh∥∥\llbracketph\rrbracket∥∥2eh−1e
###### Proof.

It is enough to bound , since the saturation assumption gives the following bound

 (32)
 ∥∥∥K−12~euh∥∥∥2 =A(~eu,~ep;~eu,~ep)=A(~eu,~ep;~eu−Π∗~eu,~ep) =L(~eu−Π∗~eu,~ep)−A(uh,ph;~eu−Π∗~eu,~ep) =−∑T∈Th{(K−1uh,~eu−Π∗~eu)T−(ph,∇⋅(~eu−Π∗~eu))T+(∇⋅uh,~ep)T} −(f,~ep)

Using Green’s formula

 ∥∥∥K−12~euh∥∥∥2= −n∑i=1⟨ph,(~eu−Π∗~eu)⋅νi⟩ΓiT3

We treat three terms in the equation separately.

 (33)

Similarly for second term, we obtain

 T2 ≤|(∇⋅uh−f,~ep−Ph~ep)T|≤∥∇⋅uh−f∥TChT∥∇~ep∥T ≤C∥∇⋅uh−f∥ThT(1+α)∥∇ep∥T ≤C∥∇⋅uh−f∥ThT∥∥−K−1u+K−1uh−K−1uh−∇ph∥∥T ≤C2∥∇⋅uh−f∥2Th2T+14∥∥K−1u−K−1uh∥∥2+14∥∥K−1uh+∇ph∥∥2T
 T3 ≤C∑e∈TΓh,T∩e≠∅∥∥\llbracketph\rrbracket∥∥eh−1/2∥~eu∥T ≤C24ε3∑e∈TΓh∥∥\llbracketph\rrbracket∥∥2eh−1+ε3∑T∈TΓh(Ω∗)∥~eu∥2T

Combining all three terms for small enough and yields

### 3.2 Lower Bound

###### Theorem 6.

Assume is polynomial with degree . There exists a constant independent of such that

 (34) ζEV≤C(∥p−ph∥+∥u−uh∥hT) (35)
###### Proof.

We have made use of a bubble function argument as it has been shown in [20, 32]

By applying the element bubble function technique we obtain

 ∥f−∇⋅uh∥ThT≤ChT∥∇⋅(u−uh)∥T≤C∥u−uh∥T=∥eu∥T

From this we conclude Inequality (34). Note that . According to the trace inequality we get

 ∥p−ph∥e ≤C(h−1/2T∥p−ph∥+h1/2T∥∇(p−ph)∥) ≤C(h−1/2T∥p−ph∥+h1/2T∥∥K−1uh+∇ph∥∥+h1/2T∥∥K−1(u−uh)∥∥) ≤C(h−1/2T∥p−ph∥+h1/2T∥∥K−1uh+∇ph∥∥+h1/2T∥∥K−1(u−uh)∥∥) ≤C(h−1/2T∥ep∥+h1/2T∥eu∥)

using above inequality for . ∎

### 3.3 Residual-based and Smoothing Estimators with Postprocessing

In this section, we show a general a posteriori error estimate in norm using suitable a polynomial functions , which is obtained by the Arbogast and Chen postprocessing of [1]. This is discussed in Section 3.4. The advantage of using the postprocessed values () lies in the fact that it leads to better indicator for pressure and flux error in the element.

#### 3.3.1 Estimates for Pressure

###### Theorem 7.

Let and . There exists a constant independent of such that

 ∥ep∥2≤C∑T∈Th{~ηP+~ηR+~ηNC} (36)

where, for all

 ~ηP =∥∥K−1uh+∇s∥∥2Th2T ~ηR,h =∥f−∇⋅uh∥2Th2T ~ηNC =∥∇(s−~ph)∥2Th2T
###### Proof.

We consider a duality argument to derive bounds. Let be the solution of the auxiliary problem

 −∇⋅K∇φ =epinΩ, (37) φ =0on∂Ω. (38)

By the elliptic regularity assumption,

 ∥φ∥2≤C∥ep∥ (39)

Let then

 A(v,φ;~v,~w)=−2∑i=1{(K−1v,~v)Ωi−(φ,∇⋅~v)Ωi+(∇⋅u,~w)Ωi}=−(ep,~w)

Then

 ∥ep