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

TL;DR

This paper develops and compares a posteriori error estimators for the Enhanced Velocity Mixed Finite Element Method in incompressible Darcy flow, emphasizing pressure and velocity error detection to improve mesh adaptivity.

Contribution

It introduces new residual-based a posteriori error estimators for EVMFEM, including pressure postprocessing, and demonstrates their effectiveness through numerical tests.

Findings

Residual-based estimators outperform explicit residual estimators.

Pressure postprocessing enhances velocity error detection.

Numerical results confirm theoretical error bounds.

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. Reference. 1. Arbogast, T., & Chen, Z. (1995). On the implementation of mixed methods as nonconforming methods for second-order elliptic problems. Mathematics of Computation, 64(211), 943-972.

Figures

Equations

$$=-\mathbf{K}\nabla p\text{in}\Omega ,$$

$$\nabla \cdot \text{u}$$

$$p$$

$$k_{min}{\xi }^T\xi \le {\xi }^T\mathbf{K}(x)\xi \le k_{max}{\xi }^T\xi \forall x\in \Omega \forall \xi \in {\mathbb{R}}^d.$$

$$\left({\mathbf{K}}^{-1}\text{u},\mathbf{v}\right)-\left(p,\nabla \cdot \mathbf{v}\right)$$

$$\left(\nabla \cdot \text{u},w\right)$$

$$\Omega ={\left(\bigcup _{i=1}^{N_b}{\bar{\Omega }}_i\right)}^o,{\Gamma }_{i,j}=\partial {\Omega }_i\bigcap \partial {\Omega }_j,\Gamma ={\left(\bigcup _{i,j=1}^{N_b}{\bar{\Gamma }}_{i,j}\right)}^o,{\Gamma }_i={\Omega }_i\bigcap \Gamma =\partial {\Omega }_j\setminus \partial \Omega .$$

$${\mathbf{V}}_h(T)=\{\mathbf{v}=(v_1,v_2)\text{or}\mathbf{v}=(v_1,v_2,v_3):v_l={\alpha }_l+{\beta }_l{x}_l:{\alpha }_l,{\beta }_l\in \mathbb{R};l=1,..d\}\text{and}$$

$$W_h(T)=\{w=\text{constant}\}$$

$$W_{h}(\Omega)=\{w\in L^{2}(\Omega):w\biggr{\rvert}_{E}\in W_{h}(T),\forall T\in\mathcal{T}_{h}\}$$

$$\displaystyle\mathbf{V}_{h,i}=\{\mathbf{v}\in H({\rm div};\Omega_{i}):\mathbf{v}\biggr{\rvert}_{T}\in\mathbf{V}_{h}(T),\forall T\in\mathcal{T}_{h,i}\}\qquad i\in\{1,...n\}$$

$${\mathbf{V}}_h=\underset{i=1}{\overset{n}{⨁}}{\mathbf{V}}_{h,i}.$$

$${\mathbf{v}}_{T_k}\cdot \nu =\{\begin{array}{l}1,\mathrm{on}{e}_k\\ 0,\mathrm{other}\mathrm{edges}\end{array}$$

$${\mathbf{V}}_h^{\ast }=\underset{i=1}{\overset{n}{⨁}}{\mathbf{V}}_{h,i}^0⨁{\mathbf{V}}_h^{\Gamma }\cap H(\mathrm{div};\Omega ).$$

$$\left(K^{-1}{\text{u}}_h,\mathbf{v}\right)$$

$$\left(\nabla \cdot {\text{u}}_h,w\right)$$

$$e_{\mathbf{u}}=\mathbf{u}-{\mathbf{u}}_h\text{ and }{e}_p=p-p_h.$$

$$\mathcal{A}(\mathbf{u},p;\mathbf{v},w)=\left({\mathbf{K}}^{-1}\text{u},\mathbf{v}\right)-\left(p,\nabla \cdot \mathbf{v}\right)+\mu \left(\nabla \cdot \text{u},w\right)$$

$$L(\mathbf{v},w)=\mu (f,w)-⟨g,\mathbf{v}\cdot \nu ⟩$$

$$\mathcal{A}(\mathbf{u},p;\mathbf{v},w)=L(\mathbf{v},w)(\mathbf{v},w)\in \mathbf{V}\times W$$

$$\mathcal{A}({\mathbf{u}}_h,p_h;\mathbf{v},w)=L(\mathbf{v},w)(\mathbf{v},w)\in {\mathbf{V}}_h^{\ast }\times W_h$$

$$\mathcal{A}(e_{\mathbf{u}},e_p;\mathbf{v},w)=L(\mathbf{v},w)-\mathcal{A}({\mathbf{u}}_h,p_h;\mathbf{v},w)(\mathbf{v},w)\in \mathbf{V}\times W$$

$$\mathcal{A}(e_{\mathbf{u}},e_p;\mathbf{v},w)=0(\mathbf{v},w)\in {\mathbf{V}}_h^{\ast }\times W_h$$

$${\Vert p\Vert }_2\le C\left(\Vert f\Vert +{\Vert g\Vert }_{\frac{3}{2},{\Gamma }_D}\right)$$

$$(w-\hat{w},w_h)=0\forall w_h\in W_h$$

$${\Vert \mathbf{v}-\Pi \mathbf{v}\Vert }_T$$

$$\Vert w-\hat{w}\Vert$$

$$\forall T\in {\mathcal{T}}_h,e\in \partial T,{\Vert \varphi \Vert }_e$$

$$\forall T\in {\mathcal{T}}_h,e\in \partial T,{\Vert \varphi \Vert }_{1/2,e}$$

$$\forall T\in {\mathcal{T}}_h,e\in \partial T,{\Vert \mathbf{v}\cdot \nu \Vert }_e$$

$$2ab\le \epsilon a^2+\frac{1}{\epsilon }{b}^2$$

$$⟨{\Pi }^{\ast }\mathbf{q}\cdot \nu ,1{⟩}_e=⟨\mathbf{q}\cdot \nu ,1{⟩}_e$$

$$\left(\nabla \cdot ({\Pi }^{\ast }\mathbf{q}-\mathbf{q}),w\right)=0\forall w\in W_h$$

$${\Vert {\Pi }^{\ast }\mathbf{v}-\mathbf{v}\Vert }_T\le C{h}_T{\Vert \mathbf{v}\Vert }_{1,T}$$

Full text

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 $\mathbf{K}$, 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 $\mathbf{K}$ of $\Omega$ 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 $L^{2}$ 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 $p$ and the Darcy velocity u:

[TABLE]

where $\Omega\in\mathbb{R}^{d}(d=2$ or $3$) is multiblock domain, $f\in L^{2}(\Omega)$ and $\mathbf{K}$ is a symmetric, uniformly positive definite tensor representing the permeability divided by the viscosity with $L^{\infty}(\Omega)$ components, for some $0<k_{min}<k_{max}<\infty$

[TABLE]

Let $\Omega$ be divided into a series of small subdomains. For simplicity, the Dirichlet boundary condition is considered as zero, i.e. $g=0$. To formulate in mixed variational form, Sobolev spaces are exploited and the following space is defined for flux in $\mathbb{R}^{d}$ as usual to be $\textbf{V}=H({\rm div};\Omega)=\{\mathbf{v}\in\left(L^{2}(\Omega)\right)^{d}:\nabla\cdot\mathbf{v}\in L^{2}(\Omega)\}$ and equipped with the norm $\left\lVert\mathbf{v}\right\rVert_{V}=\left(\left\lVert\mathbf{v}\right\rVert^{2}+\left\lVert\nabla\cdot\mathbf{v}\right\rVert^{2}\right)^{\frac{1}{2}}$ and for the pressure the space is $W=L^{2}(\Omega)$ and the corresponding norm $\left\lVert w\right\rVert_{W}=\left\lVert w\right\rVert$.

We utilize standard notations. For subdomain $\zeta\subset\mathbb{R}^{d}$, the $L^{2}(\zeta)$ inner product (or duality pairing) and norm are denoted by $(\cdot,\cdot)_{S}$ and $\left\lVert\cdot\right\rVert_{\zeta}$, respectively, for scalar and vector valued functions. Let $W^{m,p}$ be the standard Sobolev space of $m$-differentiable functions in $L^{p}(\zeta)$. Let $\left\lVert\cdot\right\rVert_{m,\zeta}$ be norm of $H^{m}(\zeta)=W^{m,2}(\zeta)$ or $H^{m}(\zeta)$, where $\zeta$ and $m$ are omitted in case of $\zeta=\Omega$ and $m=0$ respectively, in other cases they are specified. We write $(\cdot,\cdot)$ for the $L^{2}(\zeta)$ or $\left(L^{2}(\zeta)\right)^{d}$ inner product, and $\langle\cdot,\cdot\rangle_{\partial\zeta}$ for duality pairing on boundaries and interfaces, where the pairing may be between two functions in $L^{2}$ or between elements of $H^{1/2}$ and $H^{-1/2}$, in either order.

Next, a weak variational form of the fluid flow problem $(\ref{eq:a})-(\ref{eq:c})$ is to find a pair $\textbf{u}\in\mathbf{V}$, $p\in W$

[TABLE]

where $\nu$ is the outward unit normal to $\partial\Omega$.

Discrete formulation

We consider

[TABLE]

This implies that the domain is divided into $N_{b}$ subdomains, the interface between $i^{th}$ and $j^{th}$ subdomains($i\neq j$), the interior subdomain interface for $i^{th}$ subdomain and union of all such interfaces, respectively.

Let $\mathcal{T}_{h,i}$ be a conforming, quasi-uniform and rectangular partition of $\Omega_{i}$, $1\leq i\leq N_{b}$, with maximal element diameter $h_{i}$. We then set $\mathcal{T}_{h}=\cup_{i=1}^{n}\mathcal{T}_{h,i}$ and denote $h$ the maximal element diameter in $\mathcal{T}_{h}$; note that $\mathcal{T}_{h}$ can be nonmatching as neighboring meshes $\mathcal{T}_{h,i}$ and $\mathcal{T}_{h,j}$ need not match on $\Gamma_{i,j}$. 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 $d=2$ and bricks for $d=3$. The $RT_{0}$ spaces are defined for any element $T\in\mathcal{T}_{h}$ by the following spaces:

[TABLE]

In fact, a vector function in $\mathbf{V}_{h}$ can be determined uniquely by its normal components $\mathbf{v}\cdot\nu$ at midpoints of edges (in 2D) or face (in 3D) of $T$. The degrees of freedom of $\mathbf{v}\in\mathbf{V}_{h}(T)$ were created by these normal components. The degree of freedom for a pressure function $p\in W_{h}(T)$ is at center of $T$ and piecewise constant inside of $T$. The pressure finite element approximation space on $\Omega$ is taken to be as

[TABLE]

We first construct a velocity finite element approximation space on $\Omega$, which is different from the velocity space of the Multiscale Mortar Mixed FEM. Let us formulate $RT_{0}$ space on each subdomain $\Omega_{i}$ for partition $\mathcal{T}_{h}$

[TABLE]

and then

[TABLE]

Although the normal components of vectors in $\mathbf{V}_{h}$ are continuous between elements within each subdomains, the reader may see $\mathbf{V}_{h}$ is not a subspace of $H({\rm div};\Omega)$, because the normal components of the velocity vector may not match on subdomain interface $\Gamma$. 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 $\Gamma$ to finer grids, which impose the strong flux continuity between subdomains. Let us define $\mathcal{T}_{h,i,j}$ as the intersection of the traces of $\mathcal{T}_{h,i}$ and $\mathcal{T}_{h,j}$, and let $\mathcal{T}^{\Gamma}_{h}=\bigcup_{1\leq i\leq j\leq N_{b}}\mathcal{T}_{h,i,j}$. We require that $\mathcal{T}_{h,i}$ and $\mathcal{T}_{h,j}$ need to align with the coordinate axes. Fluxes are constructed to match on each element $e\in\mathcal{T}^{\Gamma}_{h}$. We consider any element $T\in\mathcal{T}_{h,i}$ that shares at least one edge with the interface $\Gamma$, i.e., $T\cap\Gamma_{i,j}\neq\emptyset$, where $1\leq i,j\leq N_{b}$ and $i\neq j$. Then newly defined interface grid introduces a partition of the edge of $T$. This partition may be extended into the element $T$ as shown in Figure 2.

This new partitioning helps to construct fine-scale fluxes that is in $H(\textbf{ div},\Omega)$. So we represent a basis function $\mathbf{v}_{T_{k}}$ in the $\mathbf{V}_{h}(T_{k})$ space ($RT_{0}$) for given $T_{k}$ with the following way:

[TABLE]

i.e. a normal component $\mathbf{v}_{T_{k}}\cdot\nu$ equal to one on $e_{k}$ and zero on all other edges(faces) of $T_{k}$. Let $\mathbf{V}^{\Gamma}_{h}$ be span of all such basis functions defined on all sub-elements induced the interface discretization $\mathcal{T}_{h,i,j}$. Thus, the enhanced velocity space $\mathbf{V}^{*}_{h}$ is taken to be as

[TABLE]

where $\mathbf{V}^{0}_{h,i}=\{\mathbf{v}\in\mathbf{V}_{h,i}:\mathbf{v}\cdot\nu=0\text{ on }\Gamma_{i}\}$ is the subspace of $\mathbf{V}_{h,i}$. The finer grid flux allows to velocity approximation on the interface and then form the $H({\rm div},\Omega)$ conforming velocity space. Some difficulties arise, however, in analysis of method and implementation of robust linear solver for such modification of $RT_{0}$ velocity space at all elements, which are adjacent to the interface $\Gamma$. We now formulate the discrete variational form of equations $(\ref{eq:a})-(\ref{eq:c})$ as: Find $\mathbf{u}_{h}\in\mathbf{V}^{*}_{h}$ and $p_{h}\in W_{h}$ such that

[TABLE]

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

[TABLE]

We introduce the bilinear form $\mathcal{A}(\mathbf{u},p;\mathbf{v},w)$ defined as

[TABLE]

where $\mu=1$ or $\mu=-1$. We denote $\mathcal{A}(\cdot,\cdot)$ as $\mathcal{A}_{s}(\cdot,\cdot)$ when $\mu=-1$, which is a symmetric bilinear form and $\mathcal{A}(\cdot,\cdot)$ as $\mathcal{A}_{c}(\cdot,\cdot)$ when $\mu=1$, which is a nonsymmetric, but coercive, since $\mathcal{A}_{c}(\mathbf{v},w;\mathbf{v},w)=(\mathbf{K}^{-1}\mathbf{v},\mathbf{v})$. Linear functional $L(\mathbf{v},w)$ is defined as

[TABLE]

where $\mu=1$ or $\mu=-1$. Note that the solution does not depend on the choice of $\mu$. Therefore, weak variational form of (5)-(6) imply that $(\mathbf{u},p)\in\mathbf{V}\times W$ satisfy

[TABLE]

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

[TABLE]

Then using above notations we obtain the residual equation

[TABLE]

hold true for each pair $(e_{\mathbf{u}},e_{p})\in\mathbf{V}\times W$. Thus, (9) and (10) give the orthogonality condition

[TABLE]

Preliminaries and useful inequalities

We assume the model problem is $H^{2}$-regular, in other words, there exists a positive C that depends on $\mathbf{K}$ and $\Omega$ such that

[TABLE]

We refer to the reader to [30] for sufficient conditions for $H^{2}$-regularity.

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

[TABLE]

Furthermore, the following important approximation properties hold true. For all $T\in\mathcal{T}_{h}$, $e\in\mathcal{T}_{h,i}|_{\partial\Omega_{i}}$, and smooth enough $\mathbf{v}$ and $w$

[TABLE]

where $C_{P}=1/\pi$, if $T$ is convex. Above inequalities are standard $L^{2}$-projection approximation results and can be found in [31]. In the analysis below we will use of trace inequalities

[TABLE]

and Young’s inequality

[TABLE]

We know that from original work [25] the projection operator $\Pi^{*}$ 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 $\Pi^{*}$ the projection operator that maps $(H^{1}(\Omega))^{d}$ onto $\mathbf{V}^{*}_{h}$ that defined locally for any element $T\in\mathcal{T}_{h}$ and any $\mathbf{q}\in(H^{1}(T))^{d}$ such that for all $\mathbf{q}\in(H^{1}(T))^{d}$

[TABLE]

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

[TABLE]

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 $\mathbf{v}\in\left(H^{1}(T)\right)^{d}$ then $\exists$ C independent of $h$ such that

[TABLE]

Lemma 2**.**

Let $\mathbf{v}\in\left(H^{1}(T)\right)^{d}$ then $\exists$ C independent of $h$ such that

[TABLE]

Proof.

[TABLE]

∎

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 $g=0$.

3.1.1 Estimates for Pressure

Theorem 3**.**

There exists a constant $C$ independent of $h$ such that

[TABLE]

where, for all $T\in\mathcal{T}_{h}$

[TABLE]

Proof.

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

[TABLE]

By the elliptic regularity assumption 12 implies that

[TABLE]

Let $\mathbf{v}=-\mathbf{K}\nabla\varphi$ then

[TABLE]

Then

[TABLE]

By Green’s formula,

[TABLE]

we obtain

[TABLE]

We use the Cauchy-Schwarz inequality and approximation properties

[TABLE]

∎

3.1.2 Estimates for Velocity with Saturation Assumption

In order to get bounds on velocity error ($e_{\mathbf{u}}$), we employ a saturation assumption. Let be $\mathbf{V}^{*}_{h_{f}},W_{h_{f}}$ be the finite element approximation spaces which corresponds to refinement of $\mathcal{T}_{h}$, where $h_{f}=h/m$ for $m\geq 2$. This implies that $\mathbf{V}^{*}_{h}\subset\mathbf{V}^{*}_{h_{f}}$ and $W_{h}\subset W_{h_{f}}$. A priori error estimates from [25] allows us to employ the following saturation assumption.

Saturation assumption

There exist constant $\beta<1$,$\alpha<1$ and $h_{f}=h/m$ for $m\geq 2$, $m\in\mathbb{N}$ such that

[TABLE]

Lemma 4**.**

Let $\mathbf{u}\in\mathbf{H}\left({\rm div},\Omega\right)$ be the exact flux defined by Eqns $1-3$ with $g=0$. Let $\mathbf{u}_{h},\mathbf{u}_{h_{f}}\in\mathbf{L}^{2}\left(\Omega\right)$ be arbitrary and $h_{f}=h/m$ for $m\geq 2$, $m\in\mathbb{N}$. If there exist $\beta\in(0,1)$ such that

[TABLE]

then,

[TABLE]

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

We write $\mathbf{u}_{h_{f}}\in\mathbf{V}^{*}_{h_{f}}$ and $p_{h_{f}}\in W_{h_{f}}$ that are the enhanced velocity mixed finite element solution of equations (7) - (8) and let

[TABLE]

We have that $(\tilde{e}_{\mathbf{u}},\tilde{e}_{p})\in\mathbf{V}^{*}_{h,f}\times W_{h,f}$ which satisfy the residual equation

[TABLE]

and the orthogonality condition

[TABLE]

Theorem 5**.**

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

[TABLE]

where, for all $T\in\mathcal{T}_{h}$

[TABLE]

Proof.

It is enough to bound $\tilde{e}_{\mathbf{u_{h}}}$, since the saturation assumption gives the following bound

[TABLE]

[TABLE]

Using Green’s formula

[TABLE]

We treat three terms in the equation separately.

[TABLE]

Similarly for second term, we obtain

[TABLE]

[TABLE]

Combining all three terms for small enough $\varepsilon_{1}$ and $\varepsilon_{3}$ yields

[TABLE]

∎

3.2 Lower Bound

Theorem 6**.**

Assume $f$ is polynomial with degree $m$. There exists a constant $C$ independent of $h$ such that

[TABLE]

Proof.

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

[TABLE]

By applying the element bubble function technique we obtain

[TABLE]

From this we conclude Inequality (34). Note that $\left\lVert\llbracket p_{h}\rrbracket\right\rVert_{e}\leq\left\lVert p-p^{+}_{h}\right\rVert_{e}+\left\lVert p-p^{-}_{h}\right\rVert_{e}$. According to the trace inequality we get

[TABLE]

using above inequality for $\tilde{\zeta}_{P}$. ∎

3.3 Residual-based and Smoothing Estimators with Postprocessing

In this section, we show a general a posteriori error estimate in $L^{2}$ norm using suitable a polynomial functions $s$, which is obtained by the Arbogast and Chen postprocessing of $p_{h}$ [1]. This is discussed in Section 3.4. The advantage of using the postprocessed values ($s$) 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 $\tilde{p}_{h}\in H^{1}(\mathcal{T}_{h,i})$ and $s\in H^{1}(\Omega_{i})$. There exists a constant $C$ independent of $h$ such that

[TABLE]

where, for all $T\in\mathcal{T}_{h}$

[TABLE]

Proof.

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

[TABLE]

By the elliptic regularity assumption,

[TABLE]

Let $\mathbf{v}=-\mathbf{K}\nabla\varphi$ then

[TABLE]

Then

[TABLE]

By Green’s formula,

[TABLE]

we obtain after using $s\in H^{1}(\Omega_{i})$

[TABLE]

We use the Cauchy-Schwarz inequality and approximation properties

[TABLE]

We apply the Young inequality and (39) inequality complete the proof. ∎

Remark 1*.*

We note that in case of assumption on $\tilde{p}_{h}$ weakly continuous between elements and across interface, the jump terms would vanish.

Remark 2*.*

We note reconstruction of velocity can be useful in defining a posteriori error. There is recently proposed the improvement of velocity at interface in [33] which is a good candidate for evaluation of velocity error.

3.3.2 Estimates for Velocity without Saturation Assumption

Theorem 8**.**

Let $s\in H^{1}(\Omega_{i})$. There exists a constant $C$ independent of $h$ such that

[TABLE]

where, for all $T\in\mathcal{T}_{h}$

[TABLE]

Proof.

For velocity we can similarly consider the following set of equations

[TABLE]

[TABLE]

$k_{min}$ can be included as $\varepsilon$ term in the Young’s inequality which helps to cancel $\mathbf{K}$ terms values.

[TABLE]

We consider a last term above expression and use the trace inequality. We note $s$ is smooth in $\Omega_{i}$, jumps only at the interface $\Gamma$.

We recall that

[TABLE]

[TABLE]

by the discrete Poincare inequality or Cauchy-Schwarz inequality.

[TABLE]

[TABLE]

Therefore,

[TABLE]

∎

3.4 Practical construction of the post-processed pressure.

For given Enhanced Velocity finite element approximation, we briefly describe the construction procedure of $\tilde{p}_{h}$, $s_{h}$ and an illustration of implementation in the two dimensional case.

3.4.1 Construction of $\tilde{p}_{h}$.

Restricting Enhanced Velocity space, we can denote $\widehat{\mathbf{V}}_{h}$ be spaces omitting interface constraints $\mathbf{V}^{\Gamma}$, so $\widehat{\mathbf{V}}_{h,i}:=\bigoplus_{i=1}^{n}\mathbf{V}_{h,i}(T)$ and then $\widehat{\mathbf{V}}_{h}:=\bigoplus_{i=1}^{n}\widehat{\mathbf{V}}_{h,i}$. Let $\mathbf{u}_{h}$, $p_{h}$ be the solution of equations (7) - (8). We first compute Lagrange multipliers for each element. We define $\lambda_{h,T}\in\Lambda_{h}$, which is piecewise constant polynomials at edge or face,

[TABLE]

where the element $T\in\mathcal{T}_{h}$ and its side $e$. We note that basis functions are same for $\mathbf{V}^{\Gamma}$ and $\widehat{\mathbf{V}}_{h}$. We employ the $L^{2}$ projected velocity from the interface, which has a finer enhanced velocity approximation, to the edge or face of subdomain element and the formulation is provided in the next subsection. We denote polynomial space $\widetilde{W}_{h}$ in the following manner

[TABLE]

where $\mathbb{Q}_{m}$ is standard notation of space that is defined in [1, 18]. We next set the post-processed $\tilde{p}_{h}$ which is proposed in [1] and the construction is performed with the following properties, for each $T\in\mathcal{T}_{h}$

[TABLE]

3.4.2 Construction of $s_{h}$.

We propose to construct the $s_{h}$ in each subdomain $\Omega_{i}$ that has the conforming mesh in order to be an efficient in computation. Construction of $s_{h}$ involves the averaging operator $\mathcal{I}_{\rm av}:\mathbb{Q}_{k}(\mathcal{T}_{h})\rightarrow\mathbb{Q}_{k}(\mathcal{T}_{h})\cap H_{0}^{1}(\Omega_{i})$. For definition of $\mathbb{Q}_{m}$ we refer reader to [1, 18]. The operator is called Oswald operator and appeared in [18, 17, 34, 7] and the analysis can be found in [35, 34]. It is interesting to note that the mapping of the gradient of pressure through Oswald operator also considered in [10]. For given $\varphi_{h}\in\mathbb{Q}_{m}(\mathcal{T}_{h})$, we regard the values of $\mathcal{I}_{\rm av}(\varphi_{h})$ as being defined at a Lagrange node $V\in\Omega$ by averaging $\varphi_{h}$ values associated this node,

[TABLE]

where $|A|$ is cardinality of sets $A$ and $\mathcal{T}_{h}$ is all collection of $T\in\mathcal{T}_{h}$ for fixed $V$. One can see that $\mathcal{T}_{h}(V)=\varphi(V)$ at those nodes that are inside of given $T\in\mathcal{T}_{h}$. We set the value of $\mathcal{I}_{\rm av}(\varphi_{h})$ is zero at boundary nodes. For the convinience of the reader we restate the relevant lemma from [18, 35, 34] that is an important property of such construction.

Lemma 9**.**

Let $\mathcal{T}_{h}$ be shape regular, let $\varphi_{h}\in\mathbb{Q}_{m}(\mathcal{T}_{h})$ and let $\mathcal{I}_{\rm av}(\varphi_{h})$ be constructed as specified above. Then

[TABLE]

for all $T\in\mathcal{T}_{h}$ and $C$ depends only on the space dimension $d$, on the maximal polynomial degree $n$, and on the shape regularity parameter $\kappa_{\mathcal{T}}$. Now in our setting we define recovered pressure $s_{h}$ for the locally post-processed $\tilde{p}_{h}$ as follows.

[TABLE]

3.4.3 Implementation steps of construction

We provide a brief steps of numerical implementation of post-processed pressure in two dimensional case. Based on piecewise pressure and velocity from the lowest order Raviart-Thomas spaces over rectangles our aim to reconstruct smoother pressure $s_{h}$. For given element $T\in\mathcal{T}_{h}(\Omega_{i})$, the main steps are

Evaluate $\lambda_{h,T}$ at edge $e_{j}$, $j=1,..4$ based on $(\mathbf{u}_{h},p_{h})$, 2. 2.

Compute $\tilde{p}_{h}$ from known $\lambda_{h,T}$, and $p_{h}$ by using equation (40), 3. 3.

Based on $\tilde{p}_{h}$ compute $s_{h}$ equation (44) at Lagrange nodes in $\Omega_{i}$.

Step 1 is standard computation of Lagrange multiplier for each element. In step 2, we are relying on higher order polynomial, in our case, it is Span$\{1,x,y,x^{2},y^{2}\}$. It is sufficient to store coefficients of polynomials in the code. In step 3, we use Span$\{1,x,y,x^{2},y^{2},xy,x^{2}y,xy^{2},x^{2}y^{2}\}$ and 9 Lagrange nodes of rectangle elements that are four rectangle nodes, four midpoints at edge and center of rectangle. This case each node requires to find neighboring elements values to compute coefficients of $s_{h}$.

4 Numerical Examples

We conduct several numerical experiments to show a posteriori error bounds for two-dimensional flow problems. We set same domain $\Omega=\left(0,1\right)\times\left(0,1\right)$ and $H/h=2$, where $H$ is coarse subdomain discretization size and $h$ is fine subdomain discretization size, for all examples. Initial subdomains grids $\mathcal{T}_{h}$ are chosen in way that has a checkerboard pattern for subdomains. Examples of such discretization are shown in Figure 3. We focus on the flux error estimators which is key in flow and transport modeling and the comparision of actual error ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{u}-\mathbf{u}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$, which is defined as follows: ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{*}:=\left\lVert\mathbf{K}^{-\frac{1}{2}}\mathbf{v}\right\rVert^{2},\;\mathbf{v}\in L^{2}(\Omega).$ The localization of the element level computation brings efficiency in the many domains setting for EV scheme methodology. In the first numerical example, the residual-based error estimators and implicit error estimators are compared for pressure and velocity errors. Second example shows the post-processed error estimators (implicit) for heterogeneous porous media. In the third example, we demonstrate the advantages of estimator, $\eta_{P}$, as indicator for three by three subdomains.

Numerical Example 1

First example tests for uniform permeability, so $\mathbf{K}=\mathbf{I}$. We report the velocity error and corresponding a posteriori error estimators. We compute the source term and boundary conditions according to the analytical solution, which is taken as follows

[TABLE]

We want to illustrate in this example of the pressure estimator. However, the importance of flux is key in the mesh refinement, since the flow coupled with transport by the flux.

As can be seen from Figure 4, the pressure error was detected well by $\eta$ and $\zeta$, however, error estimator with postprocessing is closer to the actual error. We can see from Figure 5 the residual-based estimators of velocity error are not good in detection of velocity error, on the other hand, the estimator with postprocessing captures the actual error very well. To be specific, not only error within subdomains but also error at the interface are detected. The result indicates that the fine-grid subdomains error is also predicted for further refinement process, which might depends on chosen the marking strategy in adaptive setting. In this case, we note that the interface grids equals to the fine-grid subdomains mesh size and we set the coarser grid to show increased interface error. The residual-based estimators of velocity error are not good in detection of velocity error, see Figure 5.

From now we will refer as error estimator to the post-processed error estimators. As can be seen from Figure 6, we compare convergence rate of actual and estimated flux errors against to the coarse subdomain mesh size. The coarse subdomain discretization are $h=\{1/48,1/64,1/72,1/80,1/88\}$ and the fine subdomain discretization is two times smaller, i.e.

$h_{f}=\{1/96,1/128,1/144,1/160,1/172\}$. The convergence rate are almost the same and it implies that the effectivity indexes, resulting from the ratios of the estimate over the error, is decreasing slowly to constant number. By increasing number of degrees the estimator approaches the actual error at acceptable precision level.

We have also considered different components of flux estimator, i.e. $\eta_{P}$ and $\eta_{R}$. For visual representation of the convergence rates the reader refers to Figure 7. From this comparison we can see that the estimator $\eta_{P}$ and ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{u}-\mathbf{u}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$ converges similarly as $\mathcal{O}(h^{0.9})$ and $\mathcal{O}(h^{0.8})$, respectively. It is interesting to notice that $\eta_{P}$ is a good error indicator, since it is similar to the entire-domain error distribution. On the other hand, the estimator $\eta_{R_{h}}$ converges as $\mathcal{O}(h^{2})$, because $f$ is smooth. We expected that the estimator at interface $\eta_{EV}$ converges here $\mathcal{O}(h^{0.5})$ which is slower than $\mathcal{O}(h^{0.8})$.

Numerical Example 2

We consider the a diagonal heterogeneous permeability profiles of the medium

[TABLE]

We compute the source term and boundary conditions according to the analytical solution, which is taken as follows

[TABLE]

As shown in Figure 8, the actual spatial distribution of the flux error is predicted by the estimator well. The highly heterogeneous porous medium leads to increase of the estimate values. Some part is over predicted since the residual term has large values on high permeability regions compare the potential estimator term.

As can be seen from Figure 9, the comparison of actual and estimated flux error against to the coarse subdomain mesh size. The estimate is upper bound on the error as predicted by theory. The coarse subdomain discretization are $h=\{1/36,1/48,1/64,1/80,1/100\}$ and the fine subdomain discretization is two times smaller, i.e. $h=\{1/72,1/96,1/128,1/160,1/200\}$. The convergence rate are the same and it implies that the effectivity indices, resulting from the ratios of the estimate over the error, is decreasing quickly to constant number.

We have also considered different flux estimator components $\eta_{P}$, $\eta_{R_{h}}$ and $\eta_{EV}$. For visual representation of the convergence rate the reader refers to Figure 10. From this comparison we can see that the estimator $\eta_{P}$ and ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{u}-\mathbf{u}_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{*}$ converges as similar rate $\mathcal{O}(h^{0.9})$ and $\mathcal{O}(h^{0.8})$, respectively. It is interesting to notice that $\eta_{P}$ is a good error indicator if you look at the entire-domain error distribution. On the other hand, the estimator $\eta_{R_{h}}$ converges as $\mathcal{O}(h^{2.2})$ that is much better, because $f$ is smooth.

Numerical Example 3

In this numerical example, we want to show indicator advantages in the detecting flux error for three by three subdomains, see Figure 3. It is common to see in the subsurface flow problems $\eta_{R}$ term that might have impact near the well location due to source term. The rest of domain $\eta_{P}$ can be considered as indicator of the large error element. Therefore, we suggest to use $\eta_{P}$ as an estimator in the flow and transport problems to detect potential elements for refinement. Advantage of this estimator is a local evaluation of post-processed pressure approximation.

We consider the absolute permeability as $\mathbf{K}=\mathbf{I}$. The solution of the flow problem is

[TABLE]

and then corresponding boundary conditions and force term were computed. For visual comparison of error indicator and actual error and the $\eta_{P}$ the reader referred to Figure 11. As we discussed before, the numerical example confirms that the estimator $\eta_{P}$ is a good candidate as indicator in the adaptivity strategy.

5 Conclusions

The main goal of the present research was to determine a posteriori error estimators of Enhanced Velocity discretization for the incompressible flow problems. In this paper, we derived theoretically the explicit residual-based error estimators with the saturation assumption and the implicit error estimators, which is based on post-processed pressure, and confirmed them numerically. For pressure error, the residual-based error estimators indications are similar to the implicit error estimators. On the other hand, the flux error estimators were predicted well by the implicit error estimators. In addition, the numerical result indicates that the flux error can be detected by $\eta_{P}$. In the flow and transport setting, the flux is key component in coupling and therefore identifying regions of fine scale flux by using a posteriori error analysis is key in adaptivity strategy. The findings suggest that the provided error estimates could also be useful for reservoir simulator in industrial applications. In our future research we intend to concentrate on a posteriori error analysis of transient problems and the result of priori error study has been demonstrated in [27].

Acknowledgements

First author would like to acknowledge support of the Faculty Development Competitive Research Grant (Gran No. 110119FD4502), Nazarbayev University.