Constraint Energy Minimizing Generalized Multiscale Finite Element Method in the Mixed Formulation

TL;DR

This paper introduces a mass-conservative mixed multiscale finite element method for flow in heterogeneous porous media, achieving contrast-independent first-order convergence by constructing specialized basis functions.

Contribution

It develops a novel multiscale method with auxiliary spaces and local basis functions that ensure contrast-independent accuracy in heterogeneous media.

Findings

First-order convergence rate achieved

Convergence independent of media contrast

Method effectively captures high-conductivity channels

Abstract

This paper presents a novel mass-conservative mixed multiscale method for solving flow equations in heterogeneous porous media. The media properties (the permeability) contain multiple scales and high contrast. The proposed method solves the flow equation in a mixed formulation on a coarse grid by constructing multiscale basis functions. The resulting velocity field is mass conservative on the fine grid. Our main goal is to obtain first-order convergence in terms of the mesh size which is independent of local contrast. This is achieved, first, by constructing some auxiliary spaces, which contain global information that can not be localized, in general. This is built on our previous work on the Generalized Multiscale Finite Element Method (GMsFEM). In the auxiliary space, multiscale basis functions corresponding to small (contrast-dependent) eigenvalues are selected. These basis functions represent the high-conductivity channels (which connect the boundaries of a coarse block). Next, we solve local problems to construct multiscale basis functions for the velocity field. These local problems are formulated in the oversampled domain taking into account some constraints with respect to auxiliary spaces. The latter allows fast spatial decay of local solutions and, thus, allows taking smaller oversampled regions. The number of basis functions depends on small eigenvalues of the local spectral problems. Moreover, multiscale pressure basis functions are needed in constructing the velocity space. Our multiscale spaces have a minimal dimension, which is needed to avoid contrast-dependence in the convergence. The method's convergence requires an oversampling of several layers. We present an analysis of our approach. Our numerical results confirm that the convergence rate is first order with respect to the mesh size and independent of the contrast.

Figures

Tables

Table 1. Table 1: Numerical results for the first example with 3 3 3 basis functions per element.

Number basis per element	H	# oversampling coarse layers	$e_p$	$e_v$
3	1/8	3	12.2392%	3.4897%
3	1/16	4	2.7549%	0.9931%
3	1/32	5	0.8292%	0.3227%
3	1/64	6	0.2811%	0.1098%

Table 2. Table 2: Numerical results for the first example with 1 1 1 basis function per element.

Number basis per element	H	# oversampling coarse layers	$e_p$	$e_v$
1	1/8	3	28.5354%	75.6655%
1	1/16	4	12.5646%	18.2546%
1	1/32	5	5.6206%	3.0618%
1	1/64	6	2.3959%	0.8156%

Table 3. Table 3: Numerical results for the second example with 4 4 4 basis functions per element.

Number basis per element	H	# oversampling coarse layers	$e_p$	$e_v$
4	1/8	3	58.8673%	67.3385%
4	1/16	4	7.5104%	17.9436%
4	1/32	5	2.5633%	6.5846%
4	1/64	6	0.7907%	1.9459%

Table 4. Table 4: Numerical results for the second example with 1 1 1 basis functions per element.

Number basis per element	H	# oversampling coarse layers	$e_p$	$e_v$
1	1/8	3	86.5779%	105.0276%
1	1/16	4	211.2798%	89.6828%
1	1/32	5	44.3964%	50.7899%
1	1/64	6	6.4084%	8.3319%

Table 5. Table 5: Numerical results for the third example.

Number basis per element	H	# oversampling coarse layers	$e_p$	$e_v$
4	1/10	3	8.5204%	12.7005%
4	1/20	4 (log(1/20)/log(1/10)*3=3.90)	3.7952%	7.2393%
4	1/40	5 (log(1/40)/log(1/10)*3=4.8062)	1.4931%	4.2122%

Equations

$$\text{mbox}divv=f,{\kappa }^{-1}v=-\nabla p,\text{in}\Omega ,$$

$$\begin{array}{rl}a(v,w)-b(w,p)& =0,\forall w\in V_0,\\ b(v,q)& =(f,q),\forall q\in Q,\end{array}$$

$$a(v,w)={\int }_{\Omega }{\kappa }^{-1}v\cdot w,b(w,q)={\int }_{\Omega }q\text{div}w.$$

$$\Vert q{\Vert }_{L^2(\Omega )}\le C_0\underset{v\in V_0}{\sup }\frac{b(v,q)}{\Vert v{\Vert }_{H(\text{div};\Omega )}}$$

$$\begin{array}{rl}a({\varphi }_j^{(i)},v)-b(v,p_j^{(i)})& =0,\forall v\in V_0(K_i),\\ b({\varphi }_j^{(i)},q)& ={\lambda }_j^{(i)}{s}_i(p_j^{(i)},q),\forall q\in Q(K_i),\end{array}$$

$$s_i(p,q)={\int }_{K_i}\tilde{\kappa }pq,\tilde{\kappa }=\kappa \sum _{j=1}^{N_c}|\nabla {\chi }_j{|}^2$$

$$Q_{aux}(K_i)$$

$$\pi (q)=\sum _{i=1}^N\sum _{j=1}^{J_i}\frac{s_i(q,p_j^{(i)})}{s_i(p_j^{(i)},p_j^{(i)})}{p}_j^{(i)},\forall q\in Q.$$

$$a({\psi }_{j,ms}^{(i)},v)-b(v,q_{j,ms}^{(i)})$$

$$b({\psi }_{j,ms}^{(i)},q)-s({\mu }_{j,ms}^{(i)},\pi q)$$

$$s(\pi q_{j,ms}^{(i)},\gamma )$$

$$a({\psi }_{j,ms}^{(i)},v)-b(v,q_{j,ms}^{(i)})$$

$$s(\pi q_{j,ms}^{(i)},\pi q)+b({\psi }_{j,ms}^{(i)},q)$$

$$a({\psi }_j^{(i)},v)-b(v,q_j^{(i)})$$

$$s(\pi q_j^{(i)},\pi q)+b({\psi }_j^{(i)},q)$$

$$a(\psi ,v)-b(v,r)$$

$$s(\pi r,\pi q)+b(\psi ,q)$$

$$\begin{array}{rl}a(u,v)-b(v,p)& =0,\forall v\in V_0,\\ b(u,q)& =s(p_{aux},q),\forall q\in Q,\end{array}$$

$$a({\psi }_{null},v)-b(v,q_{null})$$

$$s(\pi q_{null},\pi q)+b({\psi }_{null},q)$$

$$b(v,q_{null})=0,\forall v\in V_0.$$

$$b(v,q_{null}-{\overline{q}}_{null})=0,\forall v\in V_0$$

$$s(q_{null},q)=s(p_{null},q).$$

$$a(u_{ms},v)-b(v,p_{ms})$$

$$b(u_{ms},q)$$

$$\Vert v{\Vert }_V^2$$

$$\Vert v{\Vert }_a^2$$

$$\Vert q{\Vert }_s^2$$

$$\Vert v{\Vert }_{V(D)}^2$$

$$\Vert v{\Vert }_{a(D)}^2$$

$$\Vert q{\Vert }_{s(D)}^2$$

$$\Lambda =\underset{1\le i\le N}{\min }{\lambda }_{J_i+1}^{(i)}.$$

$$a(u_{snap},v)+b(v,p_{snap})$$

$$b(u_{snap},q)$$

$$\Vert u-u_{snap}{\Vert }_a^2+\Vert p-p_{snap}{\Vert }_{L^2(\Omega )}^2\le \frac{1}{\Lambda }\Vert (I-\pi )({\tilde{\kappa }}^{-1}f){\Vert }_s^2.$$

$$\Vert u-u_{snap}{\Vert }_V^2\le (1+\frac{1}{\Lambda })\Vert (I-\pi )({\tilde{\kappa }}^{-1}f){\Vert }_s^2.$$

Full text

1 Introduction

Multiscale problems occur in many subsurface applications. Multiple scales in permeabilities (subsurface properties) can span a large range and the variations in the permeability can be of several orders of magnitude. Moreover, the subsurface problems are often solved multiple times with different source terms and boundary conditions, e.g., in optimization of well locations. These models often admit a reduced representation and both local and global model reduction techniques have been developed in many papers. In local model reduction techniques, a reduced-order model on a coarse grid (with the grid size much larger than the heterogeneities) is sought. One of the commonly used methods includes upscaling of the permeability field [18, 9, 39]. This is done by solving local problems in each coarse block and computing the effective properties. Similar ideas have been employed for upscaling relative permeabilities in multi-phase flow simulations [5, 19]. The resulting effective medium is, then, used in a mass-conservative simulation.

An alternative approach to upscaling methods for simulations of single and multi-phase flows is the multiscale method [27, 24, 29, 23, 12, 30]. In multiscale methods, the local solutions on a coarse grid are used as basis functions in simulating flow equations. Because subsurface applications require mass conservative velocity fields, multiscale finite volume [30, 26, 32, 17, 31], mixed multiscale finite element methods [10, 1, 2, 3, 7, 8, 15], mortar multiscale methods [38, 4, 37, 14], and various post-processing methods are proposed to guarantee mass conservation [6, 34]. Obtaining a mass-conservative velocity in multiscale simulations requires special formulations.

One of the mass-conservative multiscale approaches that has been developed and widely used is based on a mixed finite element formulation [10, 1, 2, 3]. In a mixed formulation, the flow equation is formulated as two first-order equations for the pressure and velocity field. The multiscale basis functions are mainly sought for the velocity field, which is used to approximate the transport equations and requires mass conservation. Previous approaches developed multiscale basis functions for the velocity field using local solutions with Neumann boundary conditions and used piecewise constant basis functions for the pressure equations. In this case, support of the pressure basis function consists of a single coarse block, while support of the velocity basis functions consists of two coarse blocks that share a common face. These multiscale approaches have been widely developed in [10, 1, 2, 3] and applied to various multiphase applications. In these approaches, one multiscale basis function is computed for each edge (that two coarse blocks share) and spatially distributed fluxes on an internal face is imposed to reduce the subgrid effects.

In some recent works [12, 21], a Generalized Multiscale Finite Element Method has been proposed in a mixed formulation [16, 25]. In these papers, the multiscale basis functions for the velocity are computed, as before, in two coarse blocks that share a common face and the pressure basis functions are computed in each coarse block. The main difference compared to previous approaches is that a systematic approach is designed to compute multiple velocity basis functions. These velocity basis functions are computed by constructing the snapshot spaces in two coarse-block regions and then solving local spectral problems in the snapshot spaces. The typical snapshot functions consist of local solutions with prescribed Neumann boundary conditions on the common interface between two coarse blocks. The local spectral decomposition (based on the analysis) identifies appropriate local spaces. In particular, by ordering the eigenvalues in increasing order, we select the eigenvectors corresponding to small eigenvalues. The convergence analysis of these approaches is performed in [16], which shows a spectral convergence $1/\Lambda$, where $\Lambda$ is the smallest eigenvalue corresponding to the eigenvector that is not included in the coarse space. However, this convergence estimate does not include a coarse-mesh dependent convergence. The latter is our goal in this paper.

This paper follows some of the main concepts developed in [13] to design a mixed multiscale method, which has a convergence rate proportional to $H/\Lambda$, where $H$ is the coarse-mesh size and $\Lambda$ is the smallest eigenvalue that the corresponding eigenvector is not included in the coarse spaces. To obtain first order convergence in terms of the coarse-mesh size, we use some ideas from [36, 33, 35, 28]. The proposed approach requires some oversampling and a special basis construction. Our approach and the analysis significantly differ from our earlier work [13], where we designed a continuous Galerkin technique. Moreover, the proposed method provides a mass conservative velocity field. Because of the high-contrast in the permeability field, we first need to identify non-local information that depends on the contrast. This is done by solving local spectral problems for the pressure field and Constructing an auxiliary space. The auxiliary space contains non-local information represented by channels (high-conductivity regions that connect the boundaries of the coarse block). It is important to separately include this information in the coarse space (e.g., [20, 22]). To construct multiscale basis functions, we solve local problems in a large oversampled region with an appropriate orthogonality condition. The latter allows fast decay of the “remaining” information and takes into account the non-decaying part of the local solution in the oversampled domain. The proposed method provides a convergence rate $H/\Lambda$ under certain conditions on the oversampled region. Here $\Lambda$ is the minimal eigenvalue as discussed above.

We present numerical results and consider three permeability fields. In the first two cases, the permeability fields have distinct channels and the number of channels are selected to have a certain number of minimum basis functions in each coarse block. Our numerical results show that with an appropriate number of basis functions, we can achieve first-order convergence with respect to the coarse-mesh size. In our third numerical example, we consider SPE 10 Benchmark tests [11]. In this case, the number of channels is not explicitly identified and we use an auxiliary space with several basis functions. Our numerical results show that one can achieve less than $5$% error in the velocity field using four dimensional auxiliary space and oversampling.

The paper is organized as follows. In Section 2, we present some preliminaries and discuss coarse and fine grids. Section 3 is devoted to describing the multiscale method. The analysis is presented in Section 4. We present the numerical results in Section 5.

2 Preliminaries

Motivated by the importance of mass conservation for flow problems in heterogeneous media, we consider a class of high-contrast flow problems in the following mixed form

[TABLE]

subject to the homogeneous boundary condition $v\cdot\vec{n}=0$ on $\partial\Omega$, where $\Omega\subset\mathbb{R}^{d}$ is the computational domain, $d$ is the dimension and $\vec{n}$ is the outward unit normal vector on $\partial\Omega$. Note that the source function $f$ satisfies $\int_{\Omega}f=0$. We assume that $\kappa(x)$ is a heterogeneous coefficient with multiple scales and very high contrast. We assume that $1\leq\kappa(x)\leq B$, where $B$ is large. Let $V=H(\text{div};\Omega)$ and $Q=L^{2}(\Omega)$. We define $V_{0}=\{v\in V\;|\;v\cdot\vec{n}=0,\text{ on }\partial\Omega\}$. Then we can formulate the problem (1) in the following weak formulation: find $v\in V_{0}$ and $p\in Q$ such that

[TABLE]

where

[TABLE]

We remark that the above problem (2) is solved together with the condition $\int_{\Omega}p=0$. We also remark that the following inf-sup condition holds: for all $q\in L^{2}(\Omega)$ with $\int_{\Omega}q=0$, we have

[TABLE]

where $C_{0}$ is independent of $B$.

We next introduce the notion of fine and coarse grids. We let $\mathcal{T}^{H}$ be a usual conforming partition of the computational domain $\Omega$ into finite elements (triangles, quadrilaterals, tetrahedra, etc.). We refer to this partition as the coarse grid and assume that each coarse element is partitioned into a connected union of fine grid blocks, see Figure 1, where a coarse mesh is shown. The fine grid partition will be denoted by $\mathcal{T}^{h}$, and is by definition a refinement of the coarse grid $\mathcal{T}^{H}$. We let $N$ be the number of coarse elements and let $N_{c}$ be the number of coarse grid nodes of $\mathcal{T}^{H}$. We remark that our basis functions are constructed based on the coarse mesh. In order to compute numerically the basis functions on coarse elements, we need the underlying fine mesh. In the next section, we will discuss in detail the construction of basis functions.

3 The Multiscale Method

In this section, we will present the construction of our method. The construction consists of two general steps. In the first step, we will construct a multiscale space for the pressure $p$. In the second step, we will use the multiscale space for pressure to construct a multiscale space for the velocity $v$. We remark that the basis functions for pressure are local, that is, the supports of pressure basis are the coarse elements. On the other hand, for each pressure basis function, we will construct a related velocity basis function, whose support is an oversampled region containing the support of the pressure basis functions. This oversampled region is usually obtained by enlarging a coarse element by several coarse grid layers, see Figure 1. This localized feature of the velocity basis function is the key to our method.

3.1 Pressure basis functions

In this section, we will present the construction of the pressure multiscale basis functions. We will construct a set of auxiliary multiscale basis functions for each coarse element $K_{i}$ using a spectral problem. First we define some notations. For a general set $S$, we define $Q(S)=L^{2}(S)$, and $V_{0}(S)=\{v\in H(\text{div};S)\;|\;v\cdot n=0,\text{ on }\partial S\}$. Next, we define the required spectral problem. For each coarse element $K_{i}$, we solve the eigenvalue problem: find $(\phi_{j}^{(i)},p_{j}^{(i)})\in V_{0}(K{}_{i})\times Q(K_{i})$ and $\lambda^{(i)}_{j}\in\mathbb{R}$ such that

[TABLE]

where

[TABLE]

and $\{\chi_{j}\}$ is the set of standard multiscale basis functions, which satisfy the partition of unity property. One can also use other types of partition of unity functions, for example, the set of standard bilinear basis functions. Assume that $s_{i}(p_{j}^{(i)},p_{j}^{(i)})=1$ and that we arrange the eigenvalues of (5) in non-decreasing order $0=\lambda_{1}\leq\lambda_{2}\leq\dots$. To define the required multiscale basis functions, we pick the first $J_{i}$ eigenfunctions $\{p_{j}^{(i)}\}$ and define the local auxiliary multiscale finite element space $Q_{aux}(K_{i})$ by

[TABLE]

The (global) auxiliary multiscale finite element space $Q_{aux}$ is defined by $Q_{aux}=\oplus_{i}Q_{aux}(K_{i})$. Notice that, the space $Q_{aux}$ will be used as the approximation space for the pressure.

3.2 Velocity basis functions

In this section, we will present the construction of the velocity basis functions. For each of the pressure basis function in $Q_{aux}$, we will construct a corresponding velocity basis function, whose support is an oversampled region containing the support of the pressure basis function (see Figure 1).

First of all, we define the interpolation operator $\pi:Q\rightarrow Q_{aux}$ by

[TABLE]

Note that $\pi$ is the $L^{2}$ projection of $Q$ onto $Q_{aux}$ with respect to the inner product $s(p,q)=\sum_{i=1}^{N}s_{i}(p,q)$, where $s_{i}(p,q)$ is defined in (6).

Next, we give two types of velocity multiscale basis functions. These two types are motivated by the two types of constructions as in [13]. Let $p^{(i)}_{j}\in Q_{aux}$ be a given pressure basis function supported in $K_{i}$. We let $K_{i}^{+}$ be the corresponding oversampled region, see Figure 1. We will define a velocity basis function $\psi_{j,ms}^{(i)}\in V_{0}(K_{i}^{+})$ by solving the following problem. The multiscale space is defined as $V_{ms}=\text{span}\{\psi_{j,ms}^{(i)}\}$. Note that the basis function is supported in $K_{i}^{+}$, which is a union of connected coarse elements and contains $K_{i}$. We define $L_{i}$ as the set of indices such that if $k\in L_{i}$, then $K_{k}\subset K_{i}^{+}$. We also define $Q_{aux}(K_{i}^{+})=\text{span}\{p_{j}^{(k)}\;|\;1\leq j\leq J_{k},k\in L_{i}\}$.

• •

Type 1 basis functions

We find $\psi_{j,ms}^{(i)}\in V_{0}(K_{i}^{+})$, $q_{j,ms}^{(i)}\in Q(K_{i}^{+})$ and $\mu_{j,ms}^{(i)}\in Q_{aux}(K_{i}^{+})$ such that

[TABLE]

• •

Type 2 basis functions

We find $\psi_{j,ms}^{(i)}\in V_{0}(K_{i}^{+})$ and $q_{j,ms}^{(i)}\in Q(K_{i}^{+})$ such that

[TABLE]

We remark that the Type 1 basis functions are motivated by the constraint energy minimization problem defined in [13] (more precisely, equation (7) in [13]). One can show that the variational formulation of the corresponding mixed formulation of the minimization problem is exactly given by (7)-(9). The Type 2 basis functions are motivated by the unconstraint energy minimization problem in [13] (more precisely, equation (24) in [13]). One can show that the variational formulation of the corresponding mixed formulation of the minimization problem is exactly given by (10)-(11).

In the following, we will consider only the Type 2 basis functions, due to its better performance as shown in [13]. Notice that, we can define the basis function $\psi_{j,ms}^{(i)}$ in the local region $K_{i}^{+}$ is because of a localization property of the related global basis function $\psi_{j}^{(i)}$. The global basis function $\psi_{j}^{(i)}\in V_{0}$ is constructed by solving the following problem. The global multiscale space is defined as $V_{glo}=\text{span}\{\psi_{j}^{(i)}\}$.

• •

Type 2 basis functions (global)

We find $\psi_{j}^{(i)}\in V_{0}$ and $q_{j}^{(i)}\in Q$ such that

[TABLE]

Notice that the system (12)-(13) defines a mapping $G$ from $Q_{aux}$ to $V_{glo}\times Q$. In particular, given $p_{aux}\in Q_{aux}$, the image $G(p_{aux})=(\psi,r)\in V_{glo}\times Q$ is defined by

[TABLE]

We remark that the mapping $\psi=G_{1}(p_{aux})$ is surjective.

Next, we give a characterization of the space $V_{glo}$ in the following lemma.

Lemma 1**.**

Let $V_{glo}$ be the global multiscale space for the velocity. For each $p_{aux}\in Q_{aux}$ with $s(p_{aux},1)=0$, there is a unique $u\in V_{glo}$ such that $(u,p)\in V_{0}\times Q$ is the solution of

[TABLE]

and $\int_{\Omega}p=0$.

Proof.

We define the space $\widehat{V}\subset V_{0}$ using the following construction. We can define a mapping $p_{aux}\mapsto u$ using the system (16), and then define $\widehat{V}$ as the image of this map. It suffices to show that $V_{glo}=\widehat{V}$.

It is clear that $V_{glo}\subset\widehat{V}$. By the construction of the above map, we have $\text{dim}(\widehat{V})=\text{dim}(Q_{aux})-1$. Thus, we will prove next that $\text{dim}(V_{glo})=\text{dim}(Q_{aux})-1$. To do so, we will show that the null space of the mapping $G_{1}$ has dimension one. Assume that $p_{null}\in Q_{aux}$ is in the null space of $G_{1}$. Let $(\psi_{null},q_{null})=G(p_{null})$ and $\psi_{null}=G_{1}(p_{null})$. By definition of $G$,

[TABLE]

By assumption, $\psi_{null}=0$. By (17), we have

[TABLE]

Since $v\in V_{0}$, we have

[TABLE]

where $\overline{q}_{null}$ is the mean value of $q_{null}$. By the inf-sup condition 4, we conclude that $q_{null}=\overline{q}_{null}$, which implies that $q_{null}$ is a constant function. Using (18) and $Q_{aux}$ contains constant functions, we have

[TABLE]

Thus, $p_{null}$ is a constant function. Therefore, we conclude that the dimension of the null space of $G_{1}$ is one.

∎

3.3 The method

From the above, we have the multiscale spaces $Q_{aux}$ and $V_{ms}$ for the approximation of pressure and velocity. The multiscale solution $(v_{ms},p_{ms})\in V_{ms}\times Q_{aux}$ is obtained by solving

[TABLE]

To analyze the method, we will first define some norms for our spaces. We define the norms $\|\cdot\|_{V}$ and $\|\cdot\|_{a}$ for $V_{0}$ by

[TABLE]

Next, we define a norm $\|\cdot\|_{s}$ for the space $Q$ by

[TABLE]

For a given subregion $D\subset\Omega$, we will define a local version of the norms $\|\cdot\|_{a(D)}$, $\|\cdot\|_{V(D)}$ and $\|\cdot\|_{s(D)}$ by

[TABLE]

Finally, $\|\cdot\|_{L^{2}(D)}$ denotes the standard $L^{2}$ norm on $D$.

4 Analysis

The analysis consists of several steps. In Section 4.1, we will analyze the stability and the convergence of using the global basis functions. Then in Section 4.2, we will analyze the well-posedness of the construction of the multiscale basis functions. We will in Section 4.3 prove a decay property of the global multiscale basis functions, and justify the localization procedure. Finally, in Section 4.4, we will analyze the stability and the convergence of using the multiscale basis functions.

For our analysis, we define

[TABLE]

Moreover, we define the snapshot solution $(u_{snap},p_{snap})\in V_{0}\times Q$ by

[TABLE]

The well-posedness of (21)-(22) is standard. Note that $s(\pi(\tilde{\kappa}^{-1}f),q)=(f,\pi q)$ for all $q\in Q$. We remark that $(u_{snap},p_{snap})$ is considered as the reference solution is our analysis, and we will estimate the differences $u_{snap}-u_{ms}$ and $p_{snap}-p_{ms}$. In the following lemma, we will show that the differences $u-u_{snap}$ and $p-p_{snap}$ are small. That is, the snapshot error is small.

Lemma 2**.**

Let $(u_{snap},p_{snap})$ be the snapshot solution defined in (21)-(22) and let $(u,p)$ be the exact solution defined in (2). Then we have

[TABLE]

In addition, we have

[TABLE]

Proof.

First, we have

[TABLE]

Taking $v=u-u_{snap}$ and $q=p-p_{snap}$ in the above system, and summing the two equations, we have

[TABLE]

For each coarse element $K_{i}$, we define $w_{i}$ be the restriction of $(I-\pi)(p-p_{snap})$ on $K_{i}$. Then, by the definition of $\pi$, we can write

[TABLE]

We define $\mu_{i}\in V_{0}(K_{i})$ as

[TABLE]

Letting $v=\mu_{i}$ in the first equation of (25), we have

[TABLE]

Using the spectral problem (5), we have

[TABLE]

In addition, we have $a(u-u_{snap},\mu_{i})\leq\|u-u_{snap}\|_{a}\,\|\mu_{i}\|_{a}$. Then using (5), we have

[TABLE]

Combining the above results, we obtain

[TABLE]

Hence we proved (23).

Using the inf-sup condition (4) and (25), we have

[TABLE]

This proves (23).

Finally, using (26), we have

[TABLE]

This proves (24).

∎

The above lemma shows that, when the partition of unity functions $\{\chi_{i}\}$ is chosen as the standard piecewise bilinear functions, we obtain the convergence rate $\|u-u_{snap}\|_{V}+\|p-p_{snap}\|_{L^{2}(\Omega)}=O(H)$, which is independent of the contrast.

4.1 Stability and convergence of using global basis functions

In this section, we consider the use of the global basis functions to solve the problem. In particular, we approximate the problem (1) using the space $Q_{aux}$ for pressure and $V_{glo}$ for velocity. So, we define the solution $(u_{glo},p_{glo})\in V_{glo}\times Q_{aux}$ using global basis function by the following

[TABLE]

Note that $s(\pi(\tilde{\kappa}^{-1}f),q)=(f,q)$ for all $q\in Q_{aux}$. We will prove the stability and the convergence of (28)-(29).

First, we prove the following inf-sup condition.

Lemma 3**.**

For every $p_{aux}\in Q_{aux}$ with $s(p_{aux},1)=0$, there is $w\in V_{glo}$ such that

[TABLE]

Proof.

Let $p_{aux}\in Q_{aux}$ with $s(p_{aux},1)=0$. Similar to (16), we consider the following problem

[TABLE]

where the solution $w\in V_{glo}$. Taking $q=p_{aux}$ in the second equation of (31), we have $b(w,p_{aux})=\|p_{aux}\|_{s}^{2}$. Next, taking $v=w$ in the first equation of (31) and $q=p$ in the second equation of (31), we have

[TABLE]

Using the inf-sup condition (4), we have $\|p\|_{L^{2}(\Omega)}\leq C_{0}\|w\|_{a}$. This completes the proof. ∎

From the above inf-sup condition, we obtain the existence of solution of the system (28)-(29). We remark that, the constant $C_{glo}=O(\max_{x\in\Omega}\tilde{\kappa}(x))$, which depends on the contrast of the coefficient $\kappa$. We will see later that this fact leads to the conclusion that the oversampling width is the logarithm of the constrast.

Next, we will prove the following property.

Lemma 4**.**

Let $(u_{glo},p_{glo})$ be the solution of (28)-(29) and let $(u_{snap},p_{snap})$ be the snapshot solution of (21)-(22). Then we have $u_{glo}=u_{snap}$, $p_{glo}=\pi p_{snap}$ and

[TABLE]

Proof.

Note that, by Lemma 1, we see that $u_{snap}\in V_{glo}\subset V_{0}$. We note also that (28)-(29) is a conforming approximtion to the system (21)-(22), and that $\text{div}\,V_{glo}\subset\tilde{\kappa}Q_{aux}$. Thus, we conclude that $u_{snap}=u_{glo}$. Next, using (28) and (21), we have

[TABLE]

Since $\text{div}\,V_{glo}\subset\tilde{\kappa}Q_{aux}$, we conclude that $p_{glo}=\pi p_{snap}$. Finally,

[TABLE]

where the last inequality follows from a proof similar to that of (27). Finally, we note that

[TABLE]

where the first term on the right can be estiamted using (23) and the second term on the right can be estimated as $\|u\|_{a}\leq C_{0}\|f\|_{L^{2}(\Omega)}$.

∎

In next section, we will prove the existence of the global basis functions.

4.2 Well-posedness of global and multiscale basis functions

In this section, we consider the well-posedness of finding the global basis functions (defined in (12)-(13)) and the multiscale basis functions (defined in (10)-(11)). Let $S$ be a region which is a union of connected coarse elements. We will take $S=K_{i}^{+}$ for multiscale basis and take $S=\Omega$ for global basis functions. Consider the problem of finding $\psi_{j}^{(i)}\in V_{0}(S)$ and $q_{j}^{(i)}\in Q(S)$ such that

[TABLE]

We will show that the above problem has a solution. To do so, we prove the following inf-sup condition.

Lemma 5**.**

For every $q\in Q(S)$, there is $v\in V_{0}(S)$ such that

[TABLE]

Proof.

Let $q\in Q(S)$. Recall that $S$ is a union of coarse elements. Consider a coarse element $K_{i}\subset S$ and $q^{(i)}=q|_{K_{i}}$. Note that, we can write $q^{(i)}=\sum_{j}c^{i}_{j}p_{j}^{(i)}$, since the set of eigenfunctions $\{p^{(i)}_{j}\}$ forms a basis for the space $Q(K_{i})$. Using $\{\phi_{j}^{(i)}\}$ from (5), we define

[TABLE]

and define $v=\sum v^{(i)}$, where the sum is taken over all $K_{i}\subset S$. Next, we will show the required property (35). Note that, by the orthogonality of eigenfunctions,

[TABLE]

On the other hand, by (5), we have

[TABLE]

So, we have

[TABLE]

Thus, by the above and the definition of the norm $\|\cdot\|_{V}$, we have

[TABLE]

Next, by (36), we have

[TABLE]

which implies

[TABLE]

Finally, we have

[TABLE]

This completes the proof.

∎

We remark that the above lemma implies the existence of solution of (33)-(34), see Appendix A.

Finally, we give an estimate of $\|\psi_{j}^{(i)}-\psi_{j,ms}^{(i)}\|_{V}$ and $\|q_{j}^{(i)}-q_{j,ms}^{(i)}\|_{s}$.

Lemma 6**.**

Let $(\psi_{j}^{(i)},q_{j}^{(i)})$ be the solution of (12)-(13) and let $(\psi_{j,ms}^{(i)},q_{j,ms}^{(i)})$ be the solution of (10)-(11). Then we have the following approximation property

[TABLE]

for all $(v,q)\in V_{0}(K_{i}^{+})\times Q(K_{i}^{+})$.

Proof.

Using (12)-(13) and (10)-(11), we have

[TABLE]

Then, for all $(v,q)\in V_{0}(K_{i}^{+})\times Q(K_{i}^{+})$, using (38)-(39), we have

[TABLE]

Note that the first two terms on the right hand side of (40) can be estimated easily. For the other two terms on the right hand side of (40), we observe that

[TABLE]

We will next estimate the two terms on the right hand side of (41).

• •

For the first term on the right hand side of (41), we apply the Cauchy-Schwarz inequality and the triangle inequality,

[TABLE]

Note that $q-q_{j,ms}^{(i)}\in Q(K_{i}^{+})$. By the inf-sup condition (35), there is $w\in V_{0}(K_{i}^{+})$ such that

[TABLE]

Note that $\|\pi(q-q_{j,ms}^{(i)})\|_{s}\leq\|q-q_{j}^{(i)}\|_{s}+\|\pi(q_{j}^{(i)}-q_{j,ms}^{(i)})\|_{s}$. In addition, by (38), we have

[TABLE]

Using the above, we see that (42) becomes

[TABLE]

This gives the required estimate for the first term on the right hand side of (41).

• •

For the second term on the right hand side of (41), using the Cauchy-Schwarz inequality, we have

[TABLE]

We note that (39) also holds for any test function $r\in Q$, that is,

[TABLE]

since both $q_{j,ms}^{(i)}$ and $\psi_{j,ms}^{(i)}$ are zero outside $K_{i}^{+}$. Taking $z=\tilde{\kappa}^{-1}\nabla\cdot(\psi_{j}^{(i)}-\psi_{j,ms}^{(i)})\in Q$,

[TABLE]

So, we have

[TABLE]

This gives the required estimate for the second term on the right hand side of (41).

The rest of the proof follows from the Young’s inequality.

∎

4.3 Decay property of global basis functions

In this section, we will prove a decay property of the global basis functions $\psi_{j}^{(i)}$ and $q_{j}^{(i)}$. In particular, we will show that the differences $\|\psi_{j}^{(i)}-\psi_{j,ms}^{(i)}\|_{V}$ and $\|q_{j}^{(i)}-q_{j,ms}^{(i)}\|_{s}$ are small when the oversampled region $K_{i}^{+}$ is large enough.

Let $K_{i}$ be a given coarse element. We define $K_{i,m}$ as the oversampling coarse neighborhood of enlarging $K_{i}$ by $m$ coarse grid layer. See Figure 1 for an illustration of $m=2$. For $M>m$, we define $\chi_{i}^{m,M}\in\text{span}\{\chi_{i}\}$ such that

[TABLE]

We also assume $\|\nabla\chi_{i}^{M,m}\|_{L^{\infty}(\Omega)}=O(H^{-1})$. Note that, we have $|\nabla\chi_{i}^{M,m}|^{2}\leq C\sum_{j=1}^{N_{c}}|\nabla\chi_{j}|^{2}$. In the following lemma, we will prove the decay property using $K_{i}^{+}=K_{i,l}$, that is, the basis functions are constructed in a region which is a $l$ coarse grid layer extension of $K_{i}$, with $l\geq 2$. (See Figure 1).

Lemma 7**.**

Let $(\psi_{j}^{(i)},q_{j}^{(i)})$ be the solution of (12)-(13) and let $(\psi_{j,ms}^{(i)},q_{j,ms}^{(i)})$ be the solution of (10)-(11). For $K_{i}^{+}=K_{i,l}$ with $l\geq 2$, we have

[TABLE]

where $E=C(1+\cfrac{1}{\Lambda})(1+C^{-1}(1+\frac{1}{\Lambda})^{-\frac{1}{2}})^{1-l}$.

Proof.

By Lemma 6,

[TABLE]

for all $(v,q)\in V_{0}(K_{i,l})\times Q(K_{i,l})$. Next we will choose $v$ and $q$ in order to obtain the required bound. We define

[TABLE]

Then we see that $v\in V_{0}(K_{i,l})$ and $q\in Q(K_{i,l})$.

Step 1: In the first step, we will show that

[TABLE]

Using (43), we need to estimate $\|\psi_{j}^{(i)}-v\|_{V}^{2}$ and $\|q_{j}^{(i)}-q\|_{s}^{2}$. By the property of $\chi_{i}^{l,l-1}$, we see that

[TABLE]

Similarly, we have

[TABLE]

Notice that

[TABLE]

By (13), $\tilde{\kappa}^{-\frac{1}{2}}\nabla\cdot\psi_{j}^{(i)}=\tilde{\kappa}^{\frac{1}{2}}\pi(q_{j}^{(i)})$ outside $K_{i}$. So,

[TABLE]

By assumption on $\chi_{i}^{l,l-1}$, we have $|\nabla\chi_{i}^{l,l-1}|^{2}\leq C\sum_{i=1}|\nabla\chi_{i}|^{2}$. So,

[TABLE]

This completes the proof.

Step 2: In the second step, we will show that

[TABLE]

To do so, we note that $\|q_{j}^{(i)}\|_{s(\Omega\backslash K_{i,l-1})}^{2}=\|\pi q_{j}^{(i)}\|_{s(\Omega\backslash K_{i,l-1})}^{2}+\|(I-\pi)q_{j}^{(i)}\|_{s(\Omega\backslash K_{i,l-1})}^{2}$. For each coarse element $K_{m}\subset\Omega\backslash K_{i,l-1}$, we let $r_{m}$ be the restriction of $(I-\pi)q_{j}^{(i)}$ on $K_{m}$. Then we can write $r_{m}=\sum_{n>J_{m}}d_{n}^{(m)}p_{n}^{(m)}$. We define $u_{m}=\sum_{n>J_{m}}(\lambda_{n}^{(m)})^{-1}d_{n}^{(m)}\phi_{n}^{(m)}$. Notice that, by the spectral problem (5), we have $\tilde{\kappa}^{-\frac{1}{2}}\nabla\cdot\phi_{n}^{(m)}=\tilde{\kappa}^{\frac{1}{2}}\lambda_{n}^{(m)}p_{n}^{(m)}$. Thus, by the orthogonality of eigenfunctions,

[TABLE]

where $z_{m}$ is the restriction of $q_{j}^{(i)}$ on $K_{m}$. Summing the above over all $K_{m}\subset\Omega\backslash K_{i,l-1}$, we have

[TABLE]

where $u=\sum_{K_{m}\subset\Omega\backslash K_{i,l-1}}u_{m}\in V_{0}(\Omega\backslash K_{i,l-1})\subset V_{0}$. Using (12), we have

[TABLE]

To estimate $\|u\|_{a(\Omega\backslash K_{i,l-1})}$, for any $K_{m}\subset\Omega\backslash K_{i,l-1}$, using the spectral problem (5), we obtain

[TABLE]

Since $\|p_{n}^{(m)}\|_{s(K_{m})}=1$, we have

[TABLE]

Summing the above over all $K_{m}\subset\Omega\backslash K_{i,l-1}$, we have

[TABLE]

This completes the proof.

Step 3: We remark that, by using a proof similar to that of Step 2, we obtain the following result.

[TABLE]

for any $m\geq 1$.

Step 4: In the fourth step, we will prove the following inequality

[TABLE]

To do so, we let $\xi=1-\chi_{i}^{l-1,l-2}$. Then, using (12),

[TABLE]

In addition, using (13),

[TABLE]

Adding the above two equations, we have

[TABLE]

where the last step follows from Step 3.

Step 5: In the final step, we will prove the desired estimate. We notice that

[TABLE]

Using the above inequality, we obtain

[TABLE]

The proof of the lemma is completed by using (12)-(13). ∎

4.4 Stability and convergence of using multiscale basis functions

In this section, we prove the stability and the convergence of the multiscale method (19)-(20). For the results in this section, we assume the size of oversample region $K_{i,l}$ is large enough so that $E$ is small. First we prove the following inf-sup condition for the scheme (19)-(20).

Lemma 8**.**

Assume that $Cl^{2d}E(1+C_{glo}^{2})\leq\alpha<1$. For any $q_{0}\in Q_{aux}$ with $s(q_{0},1)=0$, there is $u\in V_{ms}$ such that

[TABLE]

where $C_{ms}$ is a constant.

Proof.

Let $q_{0}\in Q_{aux}$ such that $s(q_{0},1)=0$. By using the proof of (30), there is $w\in V_{glo}$ and $p\in\text{span}\{q_{j}^{(i)}\}$ such that $\|w\|_{V}\leq C_{glo}\|q_{0}\|_{s}$,

[TABLE]

and $s(p,1)=0$. Note that, since $p\in\text{span}\{q_{j}^{(i)}\}$, we can write

[TABLE]

Using (12)-(13), we can see that

[TABLE]

The above motivates the definition of $u\in V_{ms}$, which is given by

[TABLE]

We call $u$ the projection of $w\in V_{glo}$ into the space $V_{ms}$.

Next, we note that

[TABLE]

By the Cauchy-Schwarz inequality, we have

[TABLE]

By the definition of $u$ and $w$, we have

[TABLE]

By Lemma 7,

[TABLE]

Using the system (12)-(13) and the definitons of $p$ and $w$, we see that $p$ and $w$ satisfy the following

[TABLE]

where $\tilde{p}=\sum_{i=1}^{N}\sum_{j=1}^{J_{i}}d_{j}^{(i)}p_{j}^{(i)}$. Since $w\in V_{glo}$ and $V_{glo}\subset V_{0}$, the above can be written as

[TABLE]

Let $q=\tilde{p}$ in the second equation of system (48), we have

[TABLE]

By (48), $\pi p\in Q_{aux}$ and the inf-sup condition (30), $\|\pi p\|_{s}\leq C_{glo}\|w\|_{V}$. Combining the above, we have

[TABLE]

We note that

[TABLE]

Finally, using (46), we have

[TABLE]

Next, using (45), we have

[TABLE]

On the other hand, by (46),

[TABLE]

Hence, we take

[TABLE]

This completes the proof.

∎

Finally, we state and prove the main convergence theorem.

Theorem 1**.**

Assume that $(1+C_{ms})^{2}\alpha=O(1)$. Let $u_{ms}$ be the multiscale solution of (19)-(20) and let $u_{snap}$ be the snapshot solution of (21)-(22). Then we have

[TABLE]

When $\{\chi_{i}\}$ is the standard bilinear functions, we have

[TABLE]

Proof.

Using (21)-(22) and (19)-(20), we have

[TABLE]

By the construction of the multiscale basis (11), we have $\nabla\cdot V_{ms}\subset\tilde{\kappa}Q_{aux}$. Therefore, we have

[TABLE]

Thus, for any $v\in V_{ms}$, we have

[TABLE]

where we used (49) and (51). Next, taking $q=\pi(p_{snap})-p_{ms}$ in (50), we have

[TABLE]

Hence, (52) becomes

[TABLE]

Consequently, we have

[TABLE]

Using the inf-sup condition (44), there is $w\in V_{ms}$ such that

[TABLE]

Therefore,

[TABLE]

Since $u_{snap}=u_{glo}$, we have

[TABLE]

We take $v$ to be the projection of $u_{glo}$ into the space $V_{ms}$ (c.f. Lemma 8), we have

[TABLE]

For the estimate for pressure, we use (53), we have

[TABLE]

∎

We remark that, to obtain $O(H)$ convergence, we need to choose the size of the oversampling domain such that

[TABLE]

So, we see that the size of the oversampling domain $l=O(\log(\kappa/H^{2}))$.

5 Numerical results

In this section, we will present some numerical results of our proposed method. In our simulations, we take the computational domain $\Omega$ as $\Omega=[0,1]^{2}$. In the first example, we consider the medium parameter $\kappa$ shown in Figure 2 and the fine grid size $h=1/256$. There are some high-contrast inclusions and high-contrast channels with the contrast ratio of $10^{4}$. Also, the spectral problem (5) has at most $3$ small eigenvalues over all coarse elements. In the second example, we consider the medium parameter $\kappa$ shown in Figure 3 and the fine grid size $h=1/256$. In this case, there are also some high-conductivity channels with the contrast ratio of $10^{6}$, and also at most $4$ small eigenvalues for the spectral problem (5). In the third example, we consider the medium parameter $\kappa$ shown in Figure 5 and the fine grid size $h=1/200$. This is a standard SPE 10 test case [11]. To illustrate the performance of our method, we use the following error quantities

[TABLE]

where $(v_{h},p_{h})$ is the approximate solution of (2) on the fine grid $\mathcal{T}^{h}$ using the standard mixed finite element method.

In Table 1, we present the results for the first example. In this case, the source function $f$ is defined with respect to the grid with $H=1/8$. We will take $f=1$ on the coarse element $[0,1/8]\times[7/8,1]$ and $f=-1$ on the coarse element $[7/8,1]\times[0,1/8]$, and take $f=0$ otherwise. Since there are at most $3$ small eigenvalues for the spectral problem (5) over all coarse elements, our theory suggests that $3$ basis functions are needed for each coarse element. From this table, we see clearly that our method gives a convergence rate proportional to the coarse mesh size $H$. We remark that the number of coarse grid layers is predicted by our theoretical results. We have observed that the convergence rate is independent of the contrast. Moreover, one can observe that the velocity error is small even for relatively large coarse-mesh sizes ($H=1/8$).

In Table 2, we present the results for the first example with the use of only one basis function per element. In this case, our theory suggests that the decay of the global basis functions is relatively slow. Thus, with the use of only a few layers in the oversampling region, the performance of the method is worse compared with the use of $3$ basis functions per element.

In Table 3, we present the numerical results for the second example, where the source function $f$ is the same as the source used in the first example. Since there are at most $4$ small eigenvalues for the spectral problem (5) over all coarse elements, our theory suggests that $4$ basis functions are needed for each coarse element. From this table, we see that our proposed method gives a convergence rate proportional to the coarse mesh size $H$. We remark that the number of coarse grid layers is predicted by our theoretical results. In addition, we note that, the constrast ratio for this case is $100$ times larger than that of the first example. Nevertheless, we observe that the performance of our method is robust with respect to the contrast ratio of the medium parameter. Moreover, one can obtain a very good velocity approximation with $H=1/32$. In Table 4, we present the numerical results with a fewer basis functions and one can observe large errors. We also observe that there is no convergence with respect to $H$ when the coarse mesh size is relatively large. In Figure 4, we present the effect of localization with respect to the number of basis functions. We see that, when only one basis function is used, our theory predicts that the decay rate of basis function is slow. On the other hand, when all eigenfunctions with small eigenvalues ($4$ in this example) are included in the construction of basis functions, our theory predicts that the decay of basis function is fast and this is confirmed by this figure.

In our third example, we consider a SPE 10 test case shown in Figure 5, where the constrast ratio is approximately $10^{7}$. In this case, there are no clear separate channels. The corresponding results are shown in Table 5. We select the source function $f$ with respect to the grid with $H=1/10$. In particular, $f=0$ except that $f=1$ on $[1,1/10]\times[9/10,1]$ and $f=-1$ on $[9/10,1]\times[1,1/10]$. We determine the oversampling layers by using our theoretical results. We note that in this case, there are not many small eigenvalues and the number of basis functions per coarse element in the auxiliary space is chosen to be 4. In this case, we observe a good agreement in the velocity field for $H=1/20$. Moreover, from this table, we see that we obtain the convergence rate predicated by our theoretical estimates.

6 Conclusion

In this paper, we propose a new multiscale method for a class of high-contrast flow problems in the mixed formulation. We construct multiscale basis functions for both the pressure and the velocity. For the basis functions for pressure, we use some dominant eigenfunctions on each coarse elements. For the velocity basis functions, we solve some cell problems on oversampled domains using the pressure basis functions. In particular, for each pressure basis function supported on a coarse element, we will construct the corresponding velocity basis function on an oversampled region containing the coarse element. We show that this localized construction is justifed by showing an exponential decay of a corresponding global construction. Furthermore, we show that, when eigenfunctions with small eigenvalues are included in the basis, the convergence rate depends on the mesh size and is independent of the contrast of the media. Some numerical results are shown to validate our theoretical statements.

Appendix A Existence of global and multiscale basis functions

In this appendix, we prove the existence of the problem (33)-(34). It suffices to prove the following. There exists a constant $\widetilde{C}$ such that for all $(u,p)\in V_{0}(K^{+}_{i})\times Q(K^{+}_{i})$, there exist a pair $(v,q)\in V_{0}(K^{+}_{i})\times Q(K^{+}_{i})$ such that

[TABLE]

where

[TABLE]

and

[TABLE]

First, it is clear that

[TABLE]

By the infsup condition (35), there exist a $w$ such that $\|w\|_{V}=\|(I-\pi)p\|_{s}$ and

[TABLE]

Therefore,

[TABLE]

and

[TABLE]

Next, we take $(v,q)=(u+\beta w,p+\tilde{\kappa}^{-1}\nabla\cdot u)$, so

[TABLE]

Finally, we choose $\beta=C$, we have

[TABLE]

and

[TABLE]

This completes the proof.