Constraint Energy Minimizing Generalized Multiscale Finite Element Method

TL;DR

This paper develops a multiscale finite element method with basis functions designed to achieve contrast-independent convergence, using local spectral problems and constraint energy minimization within oversampling domains.

Contribution

It introduces a new GMsFEM approach that ensures contrast-independent convergence with minimal basis functions through a novel constraint energy minimization technique.

Findings

Convergence rate is independent of contrast when auxiliary space is correctly chosen.

Oversampling size weakly depends on contrast, improving efficiency.

Numerical results confirm theoretical predictions and effectiveness of the method.

Abstract

The main goal of this paper is to design multiscale basis functions within GMsFEM framework such that the convergence of method is independent of the contrast and linearly decreases with respect to mesh size if oversampling size is appropriately chosen. We would like to show a mesh-dependent convergence with a minimal number of basis functions. Our construction starts with an auxiliary multiscale space by solving local spectral problems. In auxiliary multiscale space, we select the basis functions that correspond to small (contrast-dependent) eigenvalues. These basis functions represent the channels (high-contrast features that connect the boundaries of the coarse block). Using the auxiliary space, we propose a constraint energy minimization to construct multiscale spaces. The minimization is performed in the oversampling domain, which is larger than the target coarse block. The constraints allow handling non-decaying components of the local minimizers. If the auxiliary space is correctly chosen, we show that the convergence rate is independent of the contrast (because the basis representing the channels are included in the auxiliary space) and is proportional to the coarse-mesh size (because the constrains handle non-decaying components of the local minimizers). The oversampling size weakly depends on the contrast as our analysis shows. The convergence theorem requires that channels are not aligned with the coarse edges, which hold in many applications, where the channels are oblique with respect to the coarse-mesh geometry. The numerical results confirm our theoretical results. In particular, we show that if the oversampling domain size is not sufficiently large, the errors are large. To remove the contrast-dependence of the oversampling size, we propose a modified construction for basis functions and present numerical results and the analysis.

Figures

Tables

Table 1. Table 1: Numerical results with varying coarse grid size H 𝐻 H for the test case 1 1 1 .

Number basis per element	H	# oversampling coarse layers	$L_2$error	energy error
3	1/10	4	2.62%	15.99%
3	1/20	6 (log(1/20)/log(1/10)*4=5.20)	0.51%	7.04%
3	1/40	7 (log(1/40)/log(1/10)*4=6.41)	0.11%	3.31%
3	1/80	8 (log(1/80)/log(1/10)*4=7.61)	0.0015%	0.17%

Table 2. Table 2: Comparison of various number of oversampling layers and different contrast values for test case 1 1 1 .

Layer $\backslash$ Contrast	1e+3	1e+4	1e+5	1e+6
3	48.74%	73.22%	87.83%	91.14%
4	19.12%	15.99%	26.68%	58.67%
5	3.70%	4.19%	7.17%	19.24%

Table 3. Table 3: Numerical results with varying coarse grid size H 𝐻 H for the test case 2 2 2 .

Number basis per element	H	# oversampling coarse layers	$L_2$error	energy error
4	1/10	6	1.55%	11.10%
4	1/20	8 (log(1/20)/log(1/10)*6=7.8062)	0.02%	0.59%
4	1/40	10 (log(1/40)/log(1/10)*6=9.61)	0.0042%	0.23%

Table 4. Table 4: Comparison of various number of oversampling layers and different contrast values for test case 2 2 2 .

Layer $\backslash$ Contrast	1e+3	1e+4	1e+5	1e+6
4	21.88%	47.58%	65.94%	81.47%
5	6.96%	22.49%	33.09%	47.19%
6	1.82%	4.33%	6.49%	11.09%

Table 5. Table 5: Numerical result for the test case 1 1 1 with the relaxed method.

Number basis per $K$	H	# oversampling coarse layers	$L_2$error	energy error
3	1/10	3	0.33%	3.73%
3	1/20	4 (log(1/20)/log(1/10)*3=3.9031)	0.047%	1.17%
3	1/40	5 (log(1/40)/log(1/10)*3=4.8062)	0.010%	0.47%
3	1/80	6 (log(1/40)/log(1/10)*3=5.7093)	0.0015%	0.15%

Table 6. Table 6: Comparison for the test case 1 1 1 with the relaxed method.

Layer $\backslash$ Contrast	1e+4	1e+6	1e+8	1e+10
3	3.73%	3.89%	11.99%	65.19%
4	3.72%	3.72%	3.73%	5.14%
5	3.72%	3.72%	3.73%	3.72%

Table 7. Table 7: Numerical result for the test case 2 2 2 with the relaxed method.

Number basis per element	H	# oversampling coarse layers	$L_2$error	energy error
4	1/10	4	0.11%	1.50%
4	1/20	6 (log(1/20)/log(1/10)*4=5.2041)	0.021%	0.57%
4	1/40	7 (log(1/40)/log(1/10)*4=6.4082)	0.0042%	0.23%

Table 8. Table 8: Comparison for the test case 2 2 2 with the relaxed method.

Layer $\backslash$ Contrast	1e+6	1e+8	1e+10
3	5.00%	41.11%	83.26%
4	1.50%	1.50%	11.88%
5	1.50%	1.50%	1.50%

Table 9. Table 9: Using various numbers of basis functions for the test case 2 2 2 .

Number basis per element	H	# oversampling coarse layers	$L_2$error	energy error
1	1/10	4	77.30%	87.07%
2	1/10	4	30.21%	49.66%
3	1/10	4	24.27%	44.46%
4	1/10	4	0.11%	1.50%
5	1/10	4	0.08%	1.26%

Equations

$$-\mbox{div}\big{(}\kappa(x)\,\nabla u\big{)}=f\quad\text{in}\quad\Omega\subset\mathbb{R}^{d},$$

$$a(u,v)={\int }_{\Omega }fv\text{ for all }v\in V$$

$$a(u_{ms},v)={\int }_{\Omega }fv\text{ for all }v\in V_{ms}.$$

$$\mathcal{L}({\eta }_l^{(i)})=0\text{in}{\omega }_i$$

$$\mathcal{L}({\eta }_l^{(i),+})=0\text{in}{\omega }_i^{+}$$

$$\begin{array}{r}{a}_{\Omega }(u-u_H,u-u_H)⪯\sum _{\omega }{a}_{\omega }(u^{\omega }-u_H^{\omega },u^{\omega }-u_H^{\omega }),\end{array}$$

$$a_{\omega }(\eta -{\eta }_0,\eta -{\eta }_0)⪯\delta s_{\omega }(\eta -{\eta }_0,\eta -{\eta }_0),$$

$$a_i({\varphi }_j^{(i)},w)={\lambda }_j^{(i)}{s}_i({\varphi }_j^{(i)},w),\forall w\in V(K_i),$$

$$a_i(v,w)={\int }_{K_i}\kappa \nabla v\cdot \nabla w,s_i(v,w)={\int }_{K_i}\tilde{\kappa }vw,$$

$$s(v,w)=\sum _{i=1}^N{s}_i(v,w),\Vert v{\Vert }_s=s(v,v{)}^{\frac{1}{2}},\forall v\in V_{aux}.$$

$$s(\psi ,{\varphi }_j^{(i)})=1,s(\psi ,{\varphi }_{j^{\prime }}^{(i^{\prime })})=0,\text{if }{j}^{\prime }\ne j\text{ or }{i}^{\prime }\ne i.$$

$${\pi }_i(u)=\sum _{j=1}^{l_i}\frac{s_i(u,{\varphi }_j^{(i)})}{s_i({\varphi }_j^{(i)},{\varphi }_j^{(i)})}{\varphi }_j^{(i)},\forall u\in V.$$

$$\pi (u)=\sum _{i=1}^N\sum _{j=1}^{l_i}\frac{s_i(u,{\varphi }_j^{(i)})}{s_i({\varphi }_j^{(i)},{\varphi }_j^{(i)})}{\varphi }_j^{(i)},\forall u\in V.$$

$$\psi_{j,ms}^{(i)}=\text{argmin}\Big{\{}a(\psi,\psi)\;|\;\psi\in V_{0}(K_{i,m}),\quad\psi\text{ is }\phi_{j}^{(i)}\text{-orthogonal}\Big{\}}$$

$$V_{ms}=\text{span}\Big{\{}\psi_{j,ms}^{(i)}\;|\;1\leq j\leq l_{i},\;1\leq i\leq N\Big{\}}.$$

$$\psi_{j}^{(i)}=\text{argmin}\Big{\{}a(\psi,\psi)\;|\;\psi\in V,\quad\psi\text{ is }\phi_{j}^{(i)}\text{-orthogonal}\Big{\}}.$$

$$V_{glo}=\text{span}\Big{\{}\psi_{j}^{(i)}\;|\;1\leq j\leq l_{i},\;1\leq i\leq N\Big{\}}.$$

$$a({\psi }_j^{(i)},v)=0,\forall v\in \tilde{V}.$$

$$a(u_{glo},v)={\int }_{\Omega }fv\text{ for all }v\in V_{glo}.$$

$$\Vert u-u_{glo}{\Vert }_a\le {\Lambda }^{-\frac{1}{2}}\Vert {\tilde{\kappa }}^{-\frac{1}{2}}f{\Vert }_{L^2(\Omega )}$$

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

$$\Vert u-u_{glo}{\Vert }_a\le CH{\Lambda }^{-\frac{1}{2}}\Vert {\kappa }^{-\frac{1}{2}}f{\Vert }_{L^2(\Omega )}.$$

$$a(u,v)$$

$$a(u_{glo},v)$$

$$a(u-u_{glo},v)=0,\forall v\in V_{glo}.$$

$$a(u-u_{glo},u-u_{glo})$$

$$\le \Vert {\tilde{\kappa }}^{-\frac{1}{2}}f{\Vert }_{L^2(\Omega )}\Vert u-u_{glo}{\Vert }_s.$$

$$\Vert u-u_{glo}{\Vert }_s^2=\sum _{i=1}^N\Vert (u-u_{glo}){\Vert }_{s(K_i)}^2$$

$$\sum _{i=1}^N\Vert (I-{\pi }_i)(u-u_{glo}){\Vert }_{s(K_i)}^2$$

$$a(u-u_{glo},u-u_{glo})=(f,u-u_{glo})=(f-\pi f,u-u_{glo}),$$

$$\Vert u-u_{glo}{\Vert }_a\le CH{\Lambda }^{-\frac{1}{2}}\Vert {\kappa }^{-\frac{1}{2}}(f-\pi f){\Vert }_{L^2(\Omega )}.$$

Full text

1 Introduction

Many practical applications contain multiple scales and high contrast. These include flows in fractured media, processes in channelized porous media and so on. Due to scale disparity and the contrast, some type of coarse-grid models are used to solve these problems. The coarse grid is typically much larger than the fine-grid size and it (the coarse grid) contains many heterogeneities and high contrast. In modeling and simulations of multiscale problems, it is difficult to adjust coarse-grid sizes based on scales and contrast. Thus, it is important that the numerical performance is independent of these physical parameters.

There have been many existing approaches in the literature to handle multiscale problems. In this paper, we focus on Darcy flow equation in heterogeneous media. These multiscale approaches include homogenization approaches [4, 14, 11], numerical upscaling methods [12, 7, 35, 13], multiscale finite element methods [20], variational multiscale methods [22, 5, 23, 29, 24], heterogeneous multiscale methods [15, 1, 28, 2, 18], mortar multiscale methods [33, 3, 32], localized orthogonal decomposition methods [27], equation-free approaches [25, 34, 31, 26], generalized multiscale finite element methods [16, 8, 10] and so on. Some of these approaches are based on homogenization methods and compute effective properties. Once the effective properties are computed, the global problem is solved on the coarse grid. Our methods are in the class of multiscale finite element methods, where we seek multiscale basis functions to represent the local heterogeneities. In multiscale methods, one constructs multiscale basis functions that can capture the local oscillatory behavior of the solution.

Our approaches are based on Generalized Multiscale Finite Element Method (GMsFEM), [16, 8, 10]. This approach, as MsFEM, constructs multiscale basis functions in each coarse element via local spectral problems. Once local snapshot space is constructed, the main idea of the GMsFEM is to solve local spectral problems and identify multiscale basis functions. These approaches share some common elements with multi-continuum approaches and try to identify high-contrast features that need to be represented individually. These non-local features are typically channels (high-contrast regions that connect the boundaries of the coarse grid) and need separate (individual) basis functions. These observations about representing channels separately are consistent with multi-continuum methods; however, GMsFEM provides a general framework for deriving coarse-grid equations. We note that the localizations of channels are not possible, in general, and this is the reason for constructing basis functions for channels separately as discussed in [17, 19]. These ideas are first used in designing optimal preconditioners [19]. In GMsFEM, the local spectral problems and snapshots, if identified appropriately, correctly identify the necessary channels without any geometric interpretation.

It was shown that the GMsFEM’s convergence depends on the eigenvalue decay [9]. However, it is difficult to show a coarse mesh dependent convergence without using oversampling and many basis functions. In this paper, we would like to show a mesh-dependent convergence with a minimal number of basis functions. The convergence analysis of the GMsFEM suggests that one needs to include eigenvectors corresponding to small eigenvalues in the local spectral decomposition. We note that these small eigenvalues represent the channelized features, as we discussed above. To obtain a mesh-dependent convergence, we use the ideas from [30, 27, 21]111We learned about [21] in IPAM workshop (April 2017), which is similar to Section 3 (and Section 4) and done independently and earlier by Tom Hou and Pengchuan Zhang., which consists of using oversampling domains and obtaining decaying local solutions. For high-contrast problems, the local solutions do not decay in channels and thus, we need approaches that can take into account the information in the channels when constructing the decaying local solutions.

The proposed approach starts with auxiliary multiscale basis functions constructed using the GMsFEM in each coarse block. This auxiliary space contains the information related to channels and the number of these basis functions is the same as the number of the channels, which is a minimal number of basis functions required representing high-contrast features. This auxiliary space is used to take care of the non-decaying component of the oversampled local solutions, which occurs in the channels. The construction of multiscale basis functions is done by seeking a minimization of a functional subject to a constraint such that the minimizer is orthogonal (in a certain sense) to the auxiliary space. This allows handling non-decaying component of the oversampled local solutions. The resulting approach contains several basis functions per element and one can use an adaptivity ([9]) to define the basis functions. This construction allows obtaining the convergence rate $H/\Lambda$, where $\Lambda$ is the minimal eigenvalue that the corresponding eigenvector is not included in the space. Our analysis also shows that the size of the oversampling domain depends on the contrast weakly (logarithmically). To remove the contrast dependence of the oversampling domain size, we propose a modified algorithm. In this algorithm, we use the same auxiliary space; however, the minimization is done by relaxing the constraint.

In the paper, we present numerical results for two heterogeneous permeability fields. In both cases, the permeability fields contain channels and inclusions with high conductivity values. We select auxiliary basis functions such that to include all channelized features (i.e., the eigenvectors corresponding to very small (contrast-dependent) eigenvalues). Our numerical results show that the error decays as we decrease the coarse-mesh size; however, this is sensitive to the oversampling domain size. We present numerical results that show that if there is not sufficient oversampling, the errors are large and contrast dependent. Furthermore, we also present the numerical results for our modified algorithm and show that this contrast dependence is removed and the oversampling domain sizes are less sensitive to the contrast.

The paper is organized as follows. In Section 2, we present some preliminaries. In Section 3, we present the construction of multiscale basis functions. We present the analysis of the approach in Section 4. In Section 5, we present numerical results. In Section 6, we discuss an extension and show a modified basis construction. The conclusions are presented in Section 7.

2 Preliminaries

We consider

[TABLE]

where $\kappa$ is a high-contrast with $\kappa_{0}\leq\kappa(x)\leq\kappa_{1}$ and is a multiscale field. The above equation is subjected to the boundary conditions $u=0$ on $\partial\Omega$. Next, the notions of fine and coarse grids are introduced. Let $\mathcal{T}^{H}$ be a conforming partition of $\Omega$ into finite elements. Here, $H$ is the coarse-mesh size and this partition is called coarse grid. We let $N_{c}$ be the number of vertices and $N$ be the number of elements in the coarse mesh. We assume that each coarse element is partitioned into a connected union of fine-grid blocks and this partition is called $\mathcal{T}^{h}$. Note that $\mathcal{T}^{h}$ is a refinement of the coarse grid $\mathcal{T}^{H}$ with the mesh size $h$. It is assumed that the fine grid is sufficiently fine to resolve the solution. An illustration of the fine grid, coarse grid, and oversampling domain are shown in Figure 1.

We let $V=H^{1}_{0}(\Omega)$. Then the solution $u$ of (1) satisfies

[TABLE]

where $a(u,v)=\int_{\Omega}\kappa\nabla u\cdot\nabla v$. We will discuss the construction of multiscale basis functions in the next section. We consider $V_{ms}$ to be the space spanned by all multiscale basis functions. Then the multiscale solution $u_{ms}$ is defined as the solution of the following problem, find $u_{ms}\in V_{ms}$ such that

[TABLE]

Notice that, in order to show the performance of our method, we will compute the solution of (2) on a fine mesh, which is fine enough to resolve the heterogeneities of the true solution $u$. Moreover, the construction of the multiscale basis functions is also performed on the fine mesh, even though the definition of the multiscale basis functions is constructed in the space $V$. We will give the details in the next section.

The computation of the multiscale basis functions is divided into two stages. The first stage consists of constructing the auxiliary multiscale basis functions by using the concept of generalized multiscale finite element method (GMsFEM). The next step is the construction of the multiscale basis functions. In this step, a constrained energy minimizing is performed in the oversampling domain. The construction of the multiscale basis function will be discussed in Section 3 In the next subsection, we will first introduce the basic concepts of GMsFEM.

2.1 The basic concepts of Generalized Multiscale Finite Element Method

Generalized Multiscale Finite Element Method (GMsFEM) uses two stages: offline and online. In the offline stage, a small dimensional finite element space is constructed to solve the global problem for any input parameter, such as a right-hand side or boundary condition, on a coarse grid.

The snapshot space, $V_{H,\text{snap}}^{(i)}$ is constructed for a generic domain $\omega_{i}$ or $K_{i}$. For simplicity, we use the notation $\omega_{i}$, though multiscale basis functions in Section 3 will be constructed in $K_{i}$. The snapshot solutions are used to compute multiscale basis functions. The appropriate snapshot space (1) can provide a faster convergence, (2) can provide problem relevant restrictions on the coarse spaces (e.g., divergence free solutions) and (3) can reduce the cost associated with constructing the offline spaces.

One can use various snapshot spaces (see [8]), which are (1) all fine-grid functions; (2) harmonic snapshots; (3) oversampling harmonic snapshots; and (4) foce-based snapshots. Here, we briefly discuss harmonic snapshots in oversampling domain.

We briefly discuss the snapshot space that consists of harmonic extensions of fine-grid functions that are defined on the boundary of $\omega_{i}$. For each fine-grid function, $\delta_{l}^{h}(x)$, which is defined by $\delta_{l}^{h}(x_{k})=\delta_{l,k},\,\forall x_{k}\in\textsl{J}_{h}(\omega_{i})$, where $\textsl{J}_{h}(\omega_{i})$ denotes the set of fine-grid boundary nodes on $\partial\omega_{i}$, we obtain a snapshot function $\eta_{l}^{(i)}$ by

[TABLE]

with the boundary condition, $\eta_{l}^{(i)}=\delta_{l}^{h}(x)$ on $\partial\omega_{i}$, and $\delta_{l,k}=1$ if $l=k$ and $\delta_{l,k}=0$ if $l\neq k$. We remark that the snapshot functions can be computed in the oversampling region $\omega_{i}^{+}$. In this case, for each fine-grid function, $\delta_{l}^{h}(x)$, $\delta_{l}^{h}(x_{k})=\delta_{l,k},\,\forall x_{k}\in\textsl{J}_{h}(\omega_{i}^{+})$, where $\textsl{J}_{h}(\omega_{i}^{+})$ denotes the set of fine-grid boundary nodes on $\partial\omega_{i}^{+}$, we obtain a snapshot function $\eta_{l}^{(i),+}$ by

[TABLE]

with $\eta_{l}^{(i),+}=\delta_{l}^{h}(x)$ on $\partial\omega_{i}^{+}$. Finally, we remark that one can use randomized boundary conditions to reduce the computational cost associated with the snapshot construction [8, 6].

The offline space, $V_{ms}^{(i)}$ is computed for each $\omega_{i}$ (with elements of the space denoted $\psi_{l}^{(i)}$). We perform a spectral decomposition in the snapshot space and select the dominant (corresponding to the smallest eigenvalues) to construct the offline (multiscale) space. The convergence rate of the resulting method is proportional to $1/\Lambda_{*}$, where $\Lambda_{*}$ is the smallest eigenvalue that the corresponding eigenvector is not included in the multiscale space. We would like to select local spectral problem such that we can remove many small eigenvalues with fewer multiscale basis functions.

The spectral problem depends on the analysis. In the analysis, the error is decomposed into coarse subdomains. The energy functional corresponding to the domain $\Omega$ is denoted by $a_{\Omega}(u,u)$, e.g., $a_{\Omega}(u,u)=\int_{\Omega}\kappa\nabla u\cdot\nabla u$. Then,

[TABLE]

where $\omega$ are coarse regions ($\omega_{i}$), $u^{\omega}$ is the localization of the solution. The local spectral problem is chosen to bound $a_{\omega}(u^{\omega}-u_{H}^{\omega},u^{\omega}-u_{H}^{\omega})$. We seek the subspace $V_{ms}^{\omega}$ such that for any $\eta\in V_{H,\text{snap}}^{\omega}$, there exists $\eta_{0}\in V_{ms}^{\omega}$ with,

[TABLE]

where $s_{\omega}(\cdot,\cdot)$ is an auxiliary bilinear form. The auxiliary bilinear form needs to be chosen such that the solution is bounded in the corresponding norm. Below, we will use a bilinear form defined using the mass matrix.

3 The construction of the multiscale basis functions

In this section, we will present the construction of the auxiliary multiscale basis functions. As we mentioned before, we will use the concept of GMsFEM to construct our auxiliary multiscale basis functions, which will be constructed for each coarse block $K$ in the coarse grid. Let $K_{i}$ be the $i$-th coarse block and let $V(K_{i})$ be the restriction of $V$ on $K_{i}$. Recall that for (5), we need a local spectral problem, which is to find a real number $\lambda^{(i)}_{j}$ and a function $\phi_{j}^{(i)}\in V(K_{i})$ such that

[TABLE]

where $a_{i}$ is a symmetric non-negative definite bilinear operator and $s_{i}$ is a symmetric positive definite bilinear operators defined on $V(K_{i})\times V(K_{i})$. We remark that the above problem is solved on the fine mesh in the actual computations. Based on our analysis, we can choose

[TABLE]

where $\widetilde{\kappa}=\sum_{j=1}^{N_{c}}\kappa|\nabla\chi_{j}^{ms}|^{2}$ and $\{\chi_{j}^{ms}\}_{j=1}^{Nc}$ are the standard multiscale finite element (MsFEM) basis functions (see [20]), which satisfy the partition of unity property. We let $\lambda_{j}^{(i)}$ be the eigenvalues of (6) arranged in ascending order. We will use the first $l_{i}$ eigenfunctions to construct our local auxiliary multiscale space $V_{aux}^{(i)}$, where $V_{aux}^{(i)}=\text{span}\{\phi_{j}^{(i)}|\;j\leq l_{i}\}$. The global auxiliary multiscale space $V_{aux}$ is the sum of these local auxiliary multiscale space, namely $V_{aux}=\oplus_{i=1}^{N}V^{(i)}_{aux}$. This space is used to construct multiscale basis functions that are $\phi$-orthogonal to the auxiliary space as defined above. The notion of $\phi$-orthogonality will be defined next.

For the local auxiliary multiscale space $V_{aux}^{(i)}$, the bilinear form $s_{i}$ in (3) defines an inner product with norm $\|v\|_{s(K_{i})}=s_{i}(v,v)^{\frac{1}{2}}$. These local inner products and norms provide a natural definitions of inner product and norm for the global auxiliary multiscale space $V_{aux}$, which are defined by

[TABLE]

We note that $s(v,w)$ and $\|v\|_{s}$ are also an inner product and norm for the space $V$. Using the above inner product, we can define the notion of $\phi$-orthogonality in the space $V$. Given a function $\phi_{j}^{(i)}\in V_{aux}$, we say that a function $\psi\in V$ is $\phi_{j}^{(i)}$-orthogonal if

[TABLE]

Now, we let $\pi_{i}:V\rightarrow V_{aux}^{(i)}$ be the projection with respect to the inner product $s_{i}(v,w)$. So, the operator $\pi_{i}$ is given by

[TABLE]

In addition, we let $\pi:V\rightarrow V_{aux}$ be the projection with respect to the inner product $s(v,w)$. So, the operator $\pi$ is given by

[TABLE]

Note that $\pi=\sum_{i=1}^{N}\pi_{i}$.

We next present the construction of our multiscale basis functions. For each coarse element $K_{i}$, we define an oversampled domain $K_{i,m}\subset\Omega$ by enlarging $K_{i}$ by $m$ coarse grid layers, where $m\geq 1$ is an integer. We next define the multiscale basis function $\psi_{j,ms}^{(i)}\in V_{0}(K_{i,m})$ by

[TABLE]

where $V(K_{i,m})$ is the restriction of $V$ in $K_{i,m}$, and $V_{0}(K_{i,m})$ is the subspace of $V(K_{i,m})$ with zero trace on $\partial K_{i,m}$. Our multiscale finite element space $V_{ms}$ is defined by

[TABLE]

The existence of the solution of the above minimization problem will be proved in Lemma 2, where a more general version is considered. We remark that the linear independence of the above basis functions $\psi_{j,ms}^{(i)}$ is obvious.

The following are the main ideas behind this multiscale basis function construction.

• •

The $\phi$-orthogonality of the multiscale basis functions allows a spatial decay, which is one of the contributing factors of a mesh-size dependent convergence.

• •

The multiscale basis functions minimize the energy, which is important and, in particularly, for the decay.

• •

We note if we do not choose an appropriate auxiliary space, this will affect the convergence rate and the decay rate in terms of the contrast.

In Figure 2, we illustrate the importance of the auxiliary space on the decay of multiscale basis function. We consider a high-contrast channelized medium as shown in the left plot in Figure 2. In the middle plot of Figure 2, we show a multiscale basis with the use of only one eigenfunction in the auxiliary space. We see that the basis function has almost no decay. On the other hand, in the right plot of Figure 2, we show a multiscale basis with the use of $4$ eigenfunctions, and we see clearly that the basis function has very fast decay outside the coarse block.

Remark 1**.**

The local multiscale basis construction is motivated by the global basis construction as defined below. The global basis functions are used in the convergence analysis. We will present the construction of the global basis functions. The global multiscale basis function $\psi_{j}^{(i)}\in V$ is defined by

[TABLE]

Our multiscale finite element space $V_{glo}$ is defined by

[TABLE]

This global multiscale finite element space $V_{glo}$ satisfies a very important orthogonality property, which will be used in our convergence analysis. In particular, we define $\tilde{V}_{h}$ as the null space of the projection $\pi$, namely, $\tilde{V}=\{v\in V\;|\;\pi(v)=0\}$. Then for any $\psi_{j}^{(i)}\in V_{glo}$, we have

[TABLE]

Thus, $\tilde{V}\subset V_{glo}^{\bot}$, where $V_{glo}^{\bot}$ is the orthogonal complement of $V_{glo}$ with respect to the inner product defined using the bilinear form $a$. Since $\text{dim}(V_{glo})=\text{dim}(V_{aux})$, we have $\tilde{V}=V_{glo}^{\bot}$. Thus, we have $V=V_{glo}\oplus\tilde{V}$. In Figure 3, we illustrate the decay of the basis function.

4 Analysis

In this section, we will prove the convergence of our proposed method. Before proving the convergence of the method, we need to define some notations. We will define two different norms for the finite element space $V$. One is the $a$-norm $\|\cdot\|_{a}$ where $\|u\|_{a}^{2}=\int_{\Omega}\kappa|\nabla u|^{2}$. The other is $s$-norm $\|\cdot\|_{s}$ where $\|u\|_{s}^{2}=\int_{\Omega}\tilde{\kappa}u^{2}$. For a given subdomain $\Omega_{i}\subset\Omega$, we will define the local $a$-norm and $s$-norm by $\|u\|_{a(\Omega_{i})}^{2}=\int_{\Omega_{i}}\kappa|\nabla u|^{2}$ and $\|u\|_{s(\Omega_{i})}^{2}=\int_{\Omega_{i}}\tilde{\kappa}u^{2}$.

To prove the convergence result of the proposed method, we will first show the convergence result of using the global multiscale basis functions. Next, we will give an estimate of the difference between the global basis functions and the multiscale basis functions to show the convergent rate of the proposed method is similar to using global basis functions. The approximate solution $u_{glo}\in V_{glo}$ obtained in the global multiscale space $V_{glo}$ is defined by

[TABLE]

The convergence analysis will start with the following lemma.

Lemma 1**.**

Let $u$ be the solution in (2) and $u_{glo}$ be the solution of (9). We have $u-u_{glo}\in\tilde{V}$ and

[TABLE]

where

[TABLE]

Moreover, if we replace the multiscale partition of unity $\{\chi_{j}^{ms}\}$ by the bilinear partition of unity, we have

[TABLE]

Proof.

By the definitions of $u$ and $u_{glo}$, we have

[TABLE]

Combining these two equations, we obtain

[TABLE]

So, we have $u-u_{glo}\in V_{glo}^{\bot}=\tilde{V}$. Using this orthogonality property and (2), we have

[TABLE]

Since $u-u_{glo}\in\tilde{V}_{h}$, we have $\pi(u-u_{glo})=0$. By the fact that the coarse blocks $K_{i}$ are disjoint, we also have $\pi_{i}(u-u_{glo})=0$ for all $i=1,2,\cdots,N$. Therefore, we have

[TABLE]

By using the orthogonality of the eigenfunctions $\phi_{j}^{(i)}$ of (6), we have

[TABLE]

The proof for the second part follows from the fact that $|\nabla\chi_{j}|=O(H^{-1})$ when $\{\chi_{j}\}$ is the set of bilinear partition of unity functions. ∎

We remark that, by using the fact that

[TABLE]

we can, under sufficient regularity assumption on $f$, improve the above result to

[TABLE]

After proving the above lemma, we have the convergence of the method for using global basis functions. Next, we are going to prove these global basis functions are localizable. The estimate of the difference between the global basis functions and the multiscale basis functions will require the following lemma. For each coarse block $K$, we define $B$ to be a bubble function with $B(x)>0$ for all $x\in\text{int}(K)$ and $B(x)=0$ for all $x\in\partial K$. We will take $B=\Pi_{j}\chi_{j}^{ms}$ where the product is taken over all vertices $j$ on the boundary of $K$. Using the bubble function, we define the constant

[TABLE]

We also define

[TABLE]

The following lemma considers the following minimization problem defined on a coarse block $K_{i}$:

[TABLE]

for a given $v_{aux}\in V_{aux}^{(i)}$ with $\|v_{aux}\|_{s(K_{i})}=1$, where $v_{aux}^{\perp}\subset V_{aux}^{(i)}$ is the orthogonal complement of $\text{span}\{v_{aux}\}$ with respect to the inner product $s_{i}$. We note that the minimization problem (10) is a more general version of (7).

Lemma 2**.**

For all $v_{aux}\in V_{aux}$ there exists a function $v\in V$ such that

[TABLE]

We write $D=C_{\mathcal{T}}C_{\pi}(1+\lambda_{max})$, where $C_{\mathcal{T}}$ is the square of the maximum number of vertices over all coarse elements.

Proof.

The proof consists of two steps. In the first step, we will show that the problem (10) has a unique solution. In the second step, we will prove the desired result of the lemma.

Step 1:

Let $v_{aux}\in V_{aux}^{(i)}$. The minimization problem (10) is equivalent to the following variational problem: find $v\in V_{0}(K_{i})$ and $\mu\in V_{aux}^{(i)}$ such that

[TABLE]

Note that, the well-posedness of the minimization problem (10) is equivalent to the existence of a function $v\in V_{0}(K_{i})$ such that

[TABLE]

where $C$ is a constant to be determined.

Note that $v_{aux}$ is supported in $K_{i}$. We let $v=Bv_{aux}$. By the definition of $s_{i}$, we have

[TABLE]

Since $\nabla(Bv_{aux})=v_{aux}\nabla B+B\nabla v_{aux}$, $|B|\leq 1$ and $|\nabla B|^{2}\leq C_{\mathcal{T}}\sum_{j}|\nabla\chi_{j}^{ms}|^{2}$, we have

[TABLE]

Finally, using the spectral problem (6), we have

[TABLE]

This completes the first step.

Step 2:

From the above proof, we see that the minimization problem (10) has a unique solution $v\in V_{0}(K_{i})$. So, we see that $v$ and $v_{aux}$ satisfy (11)-(12). From (12), we see that $\pi_{i}(v)=v_{aux}$. We note that the other two conditions in the lemma follow from the above proof.

∎

Notice that, we can assume $D\geq 1$ in Lemma 2.

Before estimate the difference between the global and multiscale basis function, we need some notations for the oversampling domain and the cutoff function with respect to these oversampling domains. For each $K_{i}$, we recall that $K_{i,m}\subset\Omega$ is the oversampling coarse region by enlarging $K_{i}$ by $m$ coarse grid layers. For $M>m$, we define $\chi_{i}^{M,m}\in\text{span}\{\chi^{ms}_{j}\}$ such that $0\leq\chi_{i}^{M,m}\leq 1$ and

[TABLE]

Note that, we have $K_{i,m}\subset K_{i,M}$. Moreover, $\chi_{i}^{M,m}=1$ on the inner region $K_{i,m}$ and $\chi_{i}^{M,m}=0$ outside the outer region $K_{i,M}$.

The following lemma shows that our multiscale basis functions have a decay property. In particular, the multiscale basis functions are small outside an oversampled region specified in the lemma.

Lemma 3**.**

We consider the oversampled domain $K_{i,k}$ with $k\geq 2$. That is, $K_{i,k}$ is an oversampled region by enlarging $K_{i}$ by $k$ coarse grid layers. Let $\phi_{j}^{(i)}\in V_{aux}$ be a given auxiliary multiscale basis function. We let $\psi_{j,ms}^{(i)}$ be the multiscale basis functions obtained in (7) and let $\psi_{j}^{(i)}$ be the global multiscale basis functions obtained in (8). Then we have

[TABLE]

where E=8D^{2}(1+\Lambda^{-1})\Big{(}1+\cfrac{\Lambda^{\frac{1}{2}}}{2D^{\frac{1}{2}}}\Big{)}^{1-k}.

Proof.

For the given $\phi_{j}^{(i)}\in V_{aux}$, using Lemma 2, there exist a $\tilde{\phi}_{j}^{(i)}\in V_{h}$ such that

[TABLE]

We let $\eta=\psi_{j}^{(i)}-\tilde{\phi}_{j}^{(i)}$. Note that $\eta\in\tilde{V}_{h}$ since $\pi(\eta)=0$. By using the resulting variational forms of the minimization problems (7) and (8), we see that $\psi_{j}^{(i)}$ and $\psi_{j,ms}^{(i)}$ satisfy

[TABLE]

and

[TABLE]

for some $\mu_{j}^{(i)},\mu_{j,ms}^{(i)}\in V_{aux}$. Subtracting the above two equations and restricting $v\in\tilde{V}_{0}(K_{i,k})$, we have

[TABLE]

Therefore, for $v\in\tilde{V}_{0}(K_{i,k})$, we have

[TABLE]

since $-\psi_{j,ms}^{(i)}+\tilde{\phi}_{j}^{(i)}\in\tilde{V}_{0}(K_{i,k})$. So, we obtain

[TABLE]

for all $v\in\tilde{V}_{0}(K_{i,k}).$

Next, we will estimate $\|\psi_{j}^{(i)}-\psi_{j,ms}^{(i)}\|_{a}$. We consider the $i$-th coarse block $K_{i}$. For this block, we consider two oversampled regions $K_{i,k-1}$ and $K_{i,k}$. Using these two oversampled regions, we define the cutoff function $\chi_{i}^{k,k-1}$ with the properties in (13)-(14), where we take $m=k-1$ and $M=k$. For any coarse block $K_{j}\subset K_{i,k-1}$, by (13), we have $\chi_{i}^{k,k-1}\equiv 1$ on $K_{j}$. Since $\eta\in\tilde{V}$, we have

[TABLE]

From the above result and the fact that $\chi_{i}^{k,k-1}\equiv 0$ in $\Omega\backslash K_{i,k}$, we have

[TABLE]

Using Lemma 2, for the function $\pi(\chi_{i}^{k,k-1}\eta)$, there is $\mu\in V$ such that $\text{supp}(\mu)\subset K_{i,k}\backslash K_{i,k-1}$ and $\pi(\mu-\chi_{i}^{k,k-1}\eta)=0$. Moreover, also from Lemma 2,

[TABLE]

where the last inequality follows from the fact that $\pi$ is a projection. Hence, taking $v=\mu-\chi_{i}^{k,k-1}\eta$ in (16), we have

[TABLE]

Next, we will estimate the two terms on the right hand side. We will divide the proof in four steps.

Step 1: We will estimate the term $\|(1-\chi_{i}^{k,k-1})\eta\|_{a}$ in (18). By a direct computation, we have

[TABLE]

Note that, we have $1-\chi_{i}^{k,k-1}\leq 1$. For the second term on the right hand side of the above inequality, we will use the fact that $\eta\in\tilde{V}_{h}$ and the spectral problem (6). Thus, we conclude that

[TABLE]

We will estimate the right hand side in Step 3.

Step 2: We will estimate the term $\|\mu\|_{a(K_{i,k}\backslash K_{i,k-1})}$ in (18). By (17) and $|\chi_{i}^{k,k-1}|\leq 1$, we have

[TABLE]

Combining Step 1 and Step 2, we obtain

[TABLE]

Step 3: Finally, we will estimate the term $\|\eta\|_{a(\Omega\backslash K_{i,k-1})}$. We will first show that the following recursive inequality holds

[TABLE]

where $k-2\geq 0$. Using (20) in (19), we have

[TABLE]

By using (20) again in (21), we conclude that

[TABLE]

By the definition of $\eta$ and the energy minimizing property of $\psi_{j}^{(i)}$, we have

[TABLE]

where the last inequality follows from (15).

Step 4. We will prove the estimate (20). Let $\xi=1-\chi_{i}^{k-1,k-2}$. Then we see that $\xi\equiv 1$ in $\Omega\backslash K_{i,k-1}$ and $0\leq\xi\leq 1$ otherwise. Then we have

[TABLE]

We estimate the first term in (23). For the function $\pi(\xi^{2}\eta)$, using Lemma 2, there exist $\gamma\in V$ such that $\pi(\gamma)=\pi(\xi^{2}\eta)$ and $\text{supp}(\gamma)\subset\text{supp}(\pi(\xi^{2}\eta))$. For any coarse element $K_{m}\subset\Omega\backslash K_{i,k-1}$, since $\xi\equiv 1$ on $K_{m}$, we have

[TABLE]

On the other hand, since $\xi\equiv 0$ in $K_{i,k-2}$, we have

[TABLE]

From the above two conditions, we see that $\text{supp}(\pi(\xi^{2}\eta))\subset K_{i,k-1}\backslash K_{i,k-2}$, and consequently $\text{supp}(\gamma)\subset K_{i,k-1}\backslash K_{i,k-2}$. Note that, since $\pi(\gamma)=\pi(\xi^{2}\eta)$, we have $\xi^{2}\eta-\gamma\in\tilde{V}$. We note also that $\text{supp}(\xi^{2}\eta-\gamma)\subset\Omega\backslash K_{i,k-2}$. By (15), the functions $\tilde{\phi}_{j}^{(i)}$ and $\xi^{2}\eta-\gamma$ have disjoint supports, so $a(\tilde{\phi}_{j}^{(i)},\xi^{2}\eta-\gamma)=0$. Then, by the definition of $\eta$, we have

[TABLE]

By the construction of $\psi_{j}^{(i)}$, we have $a(\psi_{j}^{(i)},\xi^{2}\eta-\gamma)=0$. Then we can estimate the first term in (23) as follows

[TABLE]

For all coarse element $K\subset K_{i,k-1}\backslash K_{i,k-2}$, since $\pi(\eta)=0$, we have

[TABLE]

Summing the above over all coarse elements $K\subset K_{i,k-1}\backslash K_{i,k-2}$, we have

[TABLE]

To estimate the second term in (23), using the spectral problem (6),

[TABLE]

Hence, by using the above results, (23) can be estimated as

[TABLE]

By using the above inequality, we have

[TABLE]

This completes the proof.

∎

The above lemma shows the global basis is localizable. Next, we use the above result to obtain an estimate of the error between the solution $u$ and the multiscale solution $u_{ms}$.

Theorem 1**.**

Let $u$ be the solution of (2) and $u_{ms}$ be the multiscale solution of (3). Then we have

[TABLE]

where $u_{glo}\in V_{glo}$ is the multiscale solution using global basis. Moreover, if $k=O(log(\cfrac{\max\{\kappa\}}{H}))$ and $\chi_{i}$ are bilinear partition of unity, we have

[TABLE]

Proof.

We write $u_{glo}=\sum_{i=1}^{N}\sum_{j=1}^{l_{i}}c_{j}^{(i)}\psi_{j}^{(i)}$. Then we define $v=\sum_{i=1}^{N}\sum_{j=1}^{l_{i}}c_{j}^{(i)}\psi_{j,ms}^{(i)}\in V_{ms}$. So, by the Galerkin orthogonality, we have

[TABLE]

Recall that the basis functions $\psi_{j,ms}^{(i)}$ have supports in $K_{i,k}$. So, by Lemma 3,

[TABLE]

where the equality follows from the orthogonality of the eigenfunctions in (6). We also note that, in the above estimate, we apply Lemma 3 to the function $\sum_{j=1}^{l_{i}}c_{j}^{(i)}(\psi_{j}^{(i)}-\psi_{j,ms}^{(i)})$. By using Lemma 1, we obtain

[TABLE]

This completes the proof for the first part of the theorem.

To proof the second inequality, we need to estimate the $s$-norm of the global solution $u_{glo}$. In particular,

[TABLE]

Since $u_{glo}$ satisfies (9), we have

[TABLE]

Therefore, we have

[TABLE]

Thus, to obtain the second inequality, we need to show that $\max\{\tilde{\kappa}\}Ck^{d}E^{\frac{1}{2}}$ is bounded. Assuming the partition of unity functions are bilinear, then we have

[TABLE]

Taking logarithm,

[TABLE]

Thus, we need $k=O(\log(\cfrac{\max\{\kappa\}}{H}))$. This completes the proof. ∎

We remark that the decay rate of the basis functions depends on the factor $1+\frac{\Lambda^{\frac{1}{2}}}{2D^{\frac{1}{2}}}$. Since $\Lambda$ is not small as eigenfunctions with small eigenvalues are used in the construction of auxiliary space, we see that, when $D$ is not large, the decay is exponential.

5 Numerical Result

In this section, we will present two numerical examples with two different high contrast media to demonstrate the convergence of our proposed method. We take the domain $\Omega=(0,1)^{2}$. For the first numerical example, we consider the medium parameter $\kappa$ as shown in Figure 4 and assume that the fine mesh size $h$ to be $1/400$. That is, the medium $\kappa$ has a $400\times 400$ resolution. In this case, we consider the contrast of the medium is $10^{4}$ where the value of $\kappa$ is large in the red region. The convergence history with various coarse mesh sizes $H$ are shown in Table 1. In all these simulations, we take the number of oversampling layer to be approximately $4\log(1/H)/\log(1/10)$. Form Table 1, we can see the energy norm error converges in first order with respect to $H$ and the $L^{2}$ norm error converges in second order with respect to $H$. The first order convergence in the energy norm matches our theoretical bound.

We emphasize that, in this example, we use $3$ basis functions per coarse region. The reason of this is that the eigenvalue problem on each coarse region has $3$ small eigenvalues, and according to our theory, we need to include the first three eigenfunctions in the auxiliary space. As our theory shows, for a certain contrast value, one needs to use a large enough oversampling size in order to obtain the desired convergence order. Moreover, for a fixed contrast value, the results can be improved as the oversampling size increases. On the other hand, for a fixed oversampling size, the performance of the scheme will deteriorate as the medium contrast increases. This is confirmed by our numerical evidence shown in Table 2.

In our second test case, we consider the medium parameter $\kappa$ defined in Figure 5. In this case, the medium has a contrast value of $10^{6}$, and the fine grid size $h$ is $1/200$. We will show the convergence history in Table 3 using different choices of coarse mesh sizes. For all simulations, the number of oversampling layer is approximately $6\log(1/H)/\log(10)$, and we use $4$ multiscale basis functions per coarse block since there are $4$ small eigenvalues for some coarse blocks. Form Table 3, we can see that the method achieves the theoretical convergence rate. In Table 4, we compare the performance of the method with different choices of oversampling layers and contrast values. We see that, for a fixed choice of oversampling layer, the error increases moderately with respect to the contrast value. On the other hand, for a fixed contrast value, the error will improve as the number of oversampling layers increases. We also see that, the error will be small once an enough number of oversampling layer is used.

6 Relaxed constraint energy minimizing generalized multiscale finite element method

In this section, we consider a relaxed version of our method. In particular, we relax the constraint in the minimization problem (7). Instead of (7), we solve the following un-constrainted minimization problem: find $\psi_{j,ms}^{(i)}\in V_{0}(K_{i,m})$ such that

[TABLE]

This minimization problem is equivalent to the following variational formulation

[TABLE]

The global multiscale basis function $\psi_{j}^{(i)}\in V$ is defined in a similar way, namely,

[TABLE]

which is equivalent to the following variational form

[TABLE]

We remark that the spaces $V_{glo}$, $V_{ms}$ and $\tilde{V}$ are defined as before.

Given a function $\tilde{v}\in\tilde{V}_{h}$, we have $\pi(\tilde{v})=0$ by the definition. Therefore $a(\psi_{j}^{(i)},\tilde{v})=0,$ that is $\tilde{V}_{h}\subset V_{glo}^{\bot}$ and since $\text{dim}(V_{glo})=\text{dim}(V_{aux})$, we have $\tilde{V}_{h}=V_{glo}^{\bot}$. Thus, we have $V_{h}=V_{glo}\oplus\tilde{V}_{h}$. Using this property, the solution $u_{glo}$ of (9) satisfies Lemma 1.

6.1 Analysis

In this section, we analyze the convergence of this method. In the following, we prove a result similar to the one in Lemma 3, with the aim of estimating the difference between $\psi_{j,ms}^{(i)}$ and $\psi_{j}^{(i)}$.

Lemma 4**.**

We consider the oversampled domain $K_{i,k}$ with $k\geq 2$. That is, $K_{i,k}$ is an oversampled region by enlarging $K_{i}$ by $k$ coarse grid layers. Let $\phi_{j}^{(i)}\in V_{aux}$ be a given auxiliary multiscale basis function. We let $\psi_{j,ms}^{(i)}$ be the multiscale basis functions obtained in (24) and let $\psi_{j}^{(i)}$ be the global multiscale basis functions obtained in (26). Then we have

[TABLE]

where E=3(1+\Lambda^{-1})\Big{(}1+(2(1+\Lambda^{-\frac{1}{2}}))^{-1}\Big{)}^{1-k}.

Proof.

By the definitions $\psi_{j,ms}^{(i)}$ and $\psi_{j}^{(i)}$ in (25) and (27), we have

[TABLE]

Subtracting the above two equations, we have

[TABLE]

for all $v\in V_{0}(K_{i,k})$. Taking $v=w-\psi_{j,ms}^{(i)}$ with $w\in V_{0}(K_{i,k})$ in the above relation, we have

[TABLE]

Let $w=\chi^{k,k-1}_{i}\psi_{j}^{(i)}$ in the above relation, we have

[TABLE]

Next, we will estimate these two terms on the right hand side of (28). We divide the proof into four steps.

Step 1: We will estimate the term $\|(1-\chi_{i}^{k,k-1})\psi_{j}^{(i)}\|^{2}_{a}$ in (28). By the definition of the norm $\|\cdot\|_{a}$ and the fact that $\text{supp}(1-\chi_{i}^{k,k-1})\subset\Omega\backslash K_{i,k-1}$, we have

[TABLE]

We note that for each $K\in\mathcal{T}^{H}$, we have

[TABLE]

Therefore, we have

[TABLE]

Step 2: We will estimate the term \|\pi\Big{(}(1-\chi_{i}^{k,k-1})\psi_{j}^{(i)}\Big{)}\|^{2}_{s} in (28). Notice that

[TABLE]

By using (29), we have

[TABLE]

From the above two steps, we see that (28) can be estimated as

[TABLE]

Next we will estimate the right hand side of (30).

Step 3: We will estimate $\|\psi_{j}^{(i)}\|_{a(\Omega\backslash K_{i,k-1})}^{2}+\|\pi(\psi_{j}^{(i)})\|_{s(\Omega\backslash K_{i,k-1})}^{2}$. We will show that this term can be estimated by $\|\psi_{j}^{(i)}\|^{2}_{(a({K_{i,k-1}\backslash K_{i,k-2}})}+\|\pi(\psi_{j}^{(i)})\|^{2}_{s({K_{i,k-1}\backslash K_{i,k-2}})}$. By using the variational form (27) and using the test function $(1-\chi_{i}^{k-1,k-2})\psi_{j}^{(i)}$, we have

[TABLE]

where the last equality follows from the facts that $\text{supp}(1-\chi_{i}^{k-1,k-2})\subset\Omega\backslash K_{i,k-2}$ and $\text{supp}(\phi_{j}^{(i)})\subset K_{i}$. Note that

[TABLE]

so we have

[TABLE]

Consequently, we have

[TABLE]

Next, we note that, since $\chi_{j}^{k-1,k-2}\equiv 0$ in $\Omega\backslash K_{i,k-1}$, we have

[TABLE]

Thus, we have

[TABLE]

Finally, summing (32) and (33) and using (31), we have

[TABLE]

where the last inequality follows from (29).

Step 4: We will show that $\|\psi_{j}^{(i)}\|_{a(\Omega\backslash K_{i,k-1})}^{2}+\|\pi(\psi_{j}^{(i)})\|_{s(\Omega\backslash K_{i,k-1})}^{2}$ can be estimated by $\|\psi_{j}^{(i)}\|_{a(\Omega\backslash K_{i,k-2})}^{2}+\|\pi(\psi_{j}^{(i)})\|_{s(\Omega\backslash K_{i,k-2})}^{2}$. This recursive property is crucial in our convergence estimate. To do so, we note that By using this, we have

[TABLE]

where we used (34) in the last inequality. Using the above inequality recursively, we have

[TABLE]

∎

Finally, we state and prove the convergence.

Theorem 2**.**

Let $u$ be the solution of (2) and $u_{ms}$ be the multiscale solution of (3). Then we have

[TABLE]

where $u_{glo}\in V_{glo}$ is the multiscale solution using global basis. Moreover, if $k=O(log(\cfrac{\max\{\kappa\}}{H}))$ and $\chi_{i}$ are bilinear partition of unity, we have

[TABLE]

Proof.

The proof follows the same procedure as the proof of Theorem 1. We write $u_{glo}=\sum_{i=1}^{N}\sum_{j=1}^{l_{i}}c_{ij}\psi_{j}^{(i)}$ and define $v=\sum_{i=1}^{N}\sum_{j=1}^{l_{i}}c_{ij}\psi_{j,ms}^{(i)}$. It suffices to estimate $\|u_{glo}-v\|_{a}$. By Lemma 4,

[TABLE]

since $\|\phi_{j}^{(i)}\|_{s}=1$. Notice that, in the above, we use Lemma 4 to the function $\sum_{j=1}^{l_{i}}c_{ij}(\psi_{j}^{(i)}-\psi_{j,ms}^{(i)})$.

Note that, we have

[TABLE]

So, we obtain

[TABLE]

Using the variational problem (27), we have

[TABLE]

Let $b_{lk}=s(\pi u_{glo},\phi_{k}^{(l)})$ and $\vec{b}=(b_{lk})$. We have

[TABLE]

where $A\in\mathbb{R}^{p\times p}$ is the matrix representation of the bilinear form \Big{(}s(\pi(\psi_{j}^{(i)}),\pi(\psi_{k}^{(l)}))+a(\psi_{j}^{(i)},\psi_{k}^{(l)})\Big{)} where $p=\sum_{i=1}^{N}l_{i}$, and $\vec{c}=(c_{ij})$. We will derive a bound for the largest eigenvalue of $A^{-1}$. For a given $\vec{c}\in\mathbb{R}^{p}$, we define an auxiliary function $\phi=\sum_{i=1}^{N}\sum_{j=1}^{l_{i}}c_{ij}\phi_{j}^{(i)}\in V_{aux}$. Using the variational problem (27), there is $\psi\in V$ such that

[TABLE]

and $\psi=\sum_{i=1}^{N}\sum_{j=1}^{l_{i}}c_{ij}\psi_{j}^{(i)}$. Using the given $\phi\in V_{aux}$, by Lemma 2, there is $z\in V$ such that

[TABLE]

Taking $w=z$ in (36),

[TABLE]

Notice that $s(\phi,\pi z)=s(\phi,\phi)=\|\vec{c}\|_{2}^{2}$. Thus,

[TABLE]

From the above, we see that the largest eigenvalue of $A^{-1}$ is $(1+D)^{\frac{1}{2}}$. So, we can estimate (35) as

[TABLE]

Using (37), we conclude that

[TABLE]

The rest of the proof follows from the proof of Theorem 1.

∎

Finally, we remark that the above theorem provides an improved bound compared with Theorem 1, since the convergence rate in Theorem 2 is independent of the constant $D$. This is confirmed by our numerical results presented next.

6.2 Numerical Result

In this section, we present numerical results to show the performance of the relaxed version of the method. As predicted by the theory, the relaxed version of the method is more robust with respect to the contrast. We will consider two test cases, which are the same as those considered in Section 5. First, in Table 5, we show the errors for the first test case using different choices of coarse mesh sizes. We clearly see that the method gives the predicted convergence rate since we have included enough eigenfunctions in the auxiliary space. More importantly, by comparing to the similar test case in Table 1, we see that the relaxed version is able to improve the number of oversampling layers. In particular, we see that one needs fewer oversampling layers and obtains much better results.

In Table 6, we show the performance of the relaxed version with respect to the relation between contrast values and number of oversampling layers. From the results, we see that the relaxed version needs a much smaller number of oversampling layers in order to achieve a good result. In particular, for a given oversampling layer, it can handle a much larger contrast value. This confirms the theoretical estimates. We performed a similar computation for the test case $2$ and obtain the same conclusion. For the numerical results, see Table 7 and Table 8

Finally, we test the performance with different choices of eigenfunctions in the auxiliary space, and the results are shown in Table 9. We consider the second medium parameter $\kappa$. For this medium, there are $4$ high-contrast channels in some coarse blocks, and therefore we need to use $4$ eigenfunctions in the auxiliary space. As predicted by our theory, using less eigenfunctions will result in a poor decay in the multiscale basis functions, the hence poor performance of the scheme. This fact is confirmed by using one, two or three basis functions per coarse blocks. We see that using $4$ basis functions will significantly improve the performance. We also note that, using more than $4$ basis functions will not necessarily improve the result further.

7 Conclusions

In this paper, we propose Constraint Energy Minimizing GMsFEM for solving flow equations in high-contrast media. The proposed method first constructs an auxiliary space, which uses eigenvectors corresponding to small eigenvalues in the local spectral problem. This space contains the subgrid information, which can not be localized and it is a minimal dimensional space that one needs for preconditioning and obtaining the errors that do not depend on the contrast. Next, using local constraint energy minimizing construction in the oversampled domain, we construct multiscale basis functions. The constraint consists of some type of orthogonality with respect to the auxiliary space. The choice of the auxiliary space is important to guarantee that the method converges as we decrease the mesh size and the convergence is independent of the contrast if the oversampling domain size is appropriately selected. Our main theorem shows that the convergence depends on the oversampling domain size that depends on the contrast. We note that this is achieved with a minimal dimensional coarse space. To remove this effect, we propose a relaxation in imposing the constraint. In our numerical results, we vary the number of oversampling layers, the contrast, the coarse-mesh size, and the number of auxiliary multiscale basis functions. Our numerical results show that one needs a minimum number of auxiliary basis functions to provide a good accuracy, which does not depend on the contrast. The numerical results are presented, which confirm our theoretical findings.

8 Acknowledgements

EC and YE would like to acknowledge the support of Hausdorff Institute of Mathematics and Institute for Pure and Applied Mathematics for hosting their long-term visits. WTL would like to acknowledge the support of Institute for Pure and Applied Mathematics for his long-term visits.