Fast Online Generalized Multiscale Finite Element Method using Constraint Energy Minimization

TL;DR

This paper introduces an improved online multiscale finite element method that adaptively constructs basis functions in oversampled regions, achieving significant error reduction with minimal iterations and providing an adaptive algorithm with controllable convergence.

Contribution

It presents a novel online basis construction approach using oversampling, enhancing error reduction and convergence in generalized multiscale finite element methods.

Findings

Achieves three orders of magnitude error reduction.

Error reduction is independent of physical parameters.

Adaptive algorithm with user-defined convergence rate.

Abstract

Local multiscale methods often construct multiscale basis functions in the offline stage without taking into account input parameters, such as source terms, boundary conditions, and so on. These basis functions are then used in the online stage with a specific input parameter to solve the global problem at a reduced computational cost. Recently, online approaches have been introduced, where multiscale basis functions are adaptively constructed in some regions to reduce the error significantly. In multiscale methods, it is desired to have only 1-2 iterations to reduce the error to a desired threshold. Using Generalized Multiscale Finite Element Framework, it was shown that by choosing sufficient number of offline basis functions, the error reduction can be made independent of physical parameters, such as scales and contrast. In this paper, our goal is to improve this. Using our recently proposed approach and special online basis construction in oversampled regions, we show that the error reduction can be made sufficiently large by appropriately selecting oversampling regions. Our numerical results show that one can achieve a three order of magnitude error reduction, which is better than our previous methods. We also develop an adaptive algorithm and enrich in selected regions with large residuals. In our adaptive method, we show that the convergence rate can be determined by a user-defined parameter and we confirm this by numerical simulations. The analysis of the method is presented.

Tables

Table 1. Table 1: Using source term f 1 subscript 𝑓 1 f_{1} and uniform enrichment.

Number of offline basis	online iteration	oversampling layers	$L_2$ error	energy error
3	0	2	0.37%	4.71%
3	1	2	6.75e-05%	1.28e-03%
3	2	2	1.57e-08%	2.64e-08%

Table 2. Table 2: Using source term f 1 subscript 𝑓 1 f_{1} and online adaptivity with θ = 0.95 𝜃 0.95 \theta=0.95 . Convergence rate is 0.9154 0.9154 0.9154 .

Number of offline basis $\omega$	DOF	oversampling layers	$L_2$ error	energy error
3	300	2	0.37%	4.71%
3	311	2	0.14%	2.21%
3	339	2	0.073%	1.12%
3	368	2	0.033%	0.57%

Table 3. Table 3: Using source term f 1 subscript 𝑓 1 f_{1} and online adaptivity with θ = 0.1 𝜃 0.1 \theta=0.1 . Convergence rate is 0.0589 0.0589 0.0589 .

Number of offline basis	DOF	oversampling layers	$L_2$ error	energy error
3	300	2	0.37%	4.71%
3	341	2	0.073%	1.09%
3	407	2	0.014%	0.21%
3	470	2	2.93e-03%	0.051%

Table 4. Table 4: Using source term f 2 subscript 𝑓 2 f_{2} and uniform enrichment.

Number of offline basis	online iteration	oversampling layers	$L_2$ error	energy error
3	0	2	1.06%	11.70%
3	1	2	6.43e-05%	1.51e-03%
3	2	2	1.57e-08%	4.25e-08%

Table 5. Table 5: Using source term f 2 subscript 𝑓 2 f_{2} and online adaptivity with θ = 0.95 𝜃 0.95 \theta=0.95 . Convergence rate is 0.9338 0.9338 0.9338 .

Number of offline basis	DOF	oversampling layers	$L_2$ error	energy error
3	300	2	1.06%	11.70%
3	309	2	0.13%	1.95%
3	324	2	0.062%	0.99%
3	347	2	0.031%	0.51%

Table 6. Table 6: Using source term f 2 subscript 𝑓 2 f_{2} and online adaptivity with θ = 0.1 𝜃 0.1 \theta=0.1 . Convergence rate is 0.09 0.09 0.09 .

Number of offline basis	DOF	oversampling layers	$L_2$ error	energy error
3	300	2	1.06%	11.70%
3	391	2	9.27e-03%	0.13%
3	513	2	3.65e-04%	6.44e-03%
3	578	2	5.31e-05%	9.89e-04%

Table 7. Table 7: Using source term f 3 subscript 𝑓 3 f_{3} and uniform enrichment.

Number of offline basis	online iteration	oversampling layers	$L_2$ error	energy error
3	0	2	30.01%	82.57%
3	1	2	0.0066%	0.0030%
3	2	2	4.45e-07%	1.22e-07%

Table 8. Table 8: Using source term f 3 subscript 𝑓 3 f_{3} and online adaptivity with θ = 0.1 𝜃 0.1 \theta=0.1 . Convergence rate is 0.0771 0.0771 0.0771 .

Number of offline basis	DOF	oversampling layers	$L_2$ error	energy error
3	300	2	30.01%	82.57%
3	356	2	8.68%	22.06%
3	378	2	4.87%	5.41%
3	392	2	4.46%	1.50%

Equations

$$-\mbox{div}\big{(}\kappa(x)\,\nabla u\big{)}=f,\quad\text{in}\;\Omega,$$

$$a(u,v)=(f,v),\forall v\in V,$$

$$a(u_{ms},v)=(f,v),\forall v\in V_{ms}.$$

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

$$V_{aux}(K_i)=\text{span}\{{\varphi }_j^{(i)}|1\le j\le J_i\}.$$

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

$$\pi (v)=\sum _{i=1}^N\sum _{j=1}^{J_i}{s}_i(v,{\varphi }_j^{(i)}){\varphi }_j^{(i)},\forall v\in V.$$

$$\tilde{V}=\{w\in V|\pi (w)=0\}.$$

$$a({\psi }_{j,ms}^{(i)},v)+s(\pi ({\psi }_{j,ms}^{(i)}),\pi (v))=s({\varphi }_j^{(i)},\pi (v)),\forall v\in V_0(K_i^{+}),$$

$$a({\psi }_j^{(i)},v)+s(\pi ({\psi }_j^{(i)}),\pi (v))=s({\varphi }_j^{(i)},\pi (v)),\forall v\in V.$$

$$r(v)=a(u_{ms},v)-{\int }_{\Omega }fv,\forall v\in V.$$

$$r_i(v)=r({\chi }_{i}v),\forall v\in V.$$

$$a({\beta }_{ms}^{(i)},v)+s(\pi ({\beta }_{ms}^{(i)}),\pi (v))=r_i(v),\forall v\in V_0({\omega }_i^{+}),$$

$$a({\beta }_{glo}^{(i)},v)+s(\pi ({\beta }_{glo}^{(i)}),\pi (v))=r_i(v),\forall v\in V.$$

$$a(u_{ms}^{(m)},v)=(f,v),\text{for all}v\in V_{ms}^{(m)}.$$

$$z_i(v)=a(u_{ms}^{(m)},v)-(f,v),\forall v\in V_0({\omega }_i).$$

$$\sum _{i=k+1}^N{\delta }_i^2<\theta \sum _{i=1}^N{\delta }_i^2.$$

$$a({\beta }_{ms}^{(i)},v)+s(\pi ({\beta }_{ms}^{(i)}),\pi (v))=r_i^{(m)}(v)\forall v\in V_0({\omega }_i^{+})$$

$$V_{ms}^{(m+1)}=V_{ms}^{(m)}+{\text{span}}_{1\le i\le k}\{{\beta }_{ms}^{(i)}\}.$$

$$\|u-u_{ms}^{(m+1)}\|^{2}_{a}\leq 3(1+\Lambda^{-1})\Big{(}C(\ell+1)^{d}E+2M^{2}\theta\Big{)}\|u-u^{(m)}_{ms}\|^{2}_{a}$$

$$\pi (v)=v_{aux},\Vert v{\Vert }_a^2\le D\Vert v_{aux}{\Vert }_s^2,\text{supp}(v)\subset \text{supp}(v_{aux}).$$

$$\|\psi_{j}^{(i)}-\psi_{j,ms}^{(i)}\|_{a}^{2}+\|\pi(\psi_{j}^{(i)}-\psi_{j,ms}^{(i)})\|_{s}^{2}\leq E\Big{(}\|\psi_{j}^{(i)}\|_{a}^{2}+\|\pi(\psi_{j}^{(i)})\|_{s}^{2}\Big{)}$$

$$\Vert \sum _{i=1}^N\sum _{j=1}^{J_i}{c}_j^{(i)}({\psi }_j^{(i)}-{\psi }_{j,ms}^{(i)}){\Vert }_a^2+\Vert \pi (\sum _{i=1}^N\sum _{j=1}^{J_i}{c}_j^{(i)}({\psi }_j^{(i)}-{\psi }_{j,ms}^{(i)})){\Vert }_s^2$$

$$\le$$

$$\|\beta_{glo}^{(i)}-\beta_{ms}^{(i)}\|_{a}^{2}+\|\pi(\beta_{glo}^{(i)}-\beta_{ms}^{(i)})\|_{s}^{2}\leq E\Big{(}\|\beta_{glo}^{(i)}\|_{a}^{2}+\|\pi(\beta_{glo}^{(i)})\|_{s}^{2}\Big{)}$$

$$\Vert \sum _{i=1}^{N_c}({\beta }_{glo}^{(i)}-{\beta }_{ms}^{(i)}){\Vert }_a^2+\Vert \pi (\sum _{i=1}^{N_c}({\beta }_{glo}^{(i)}-{\beta }_{ms}^{(i)})){\Vert }_s^2$$

$$\le$$

$$C_0=\underset{v\in V}{\sup }\frac{\Vert \pi (v){\Vert }_s^2}{\Vert v{\Vert }_a^2}.$$

$$\Vert u-u_{ms}^{(m+1)}{\Vert }_a\le \Vert u-w{\Vert }_a,\forall w\in V_{ms}^{(m+1)}.$$

$$a({\beta }_{glo}^{(i)},v)+s(\pi ({\beta }_{glo}^{(i)}),\pi (v))=r_i^{(m)}(v),\forall v\in V.$$

$$a(\sum _{i=1}^{N_c}{\beta }_{glo}^{(i)},v)+s(\pi (\sum _{i=1}^{N_c}{\beta }_{glo}^{(i)}),\pi (v))$$

Full text

1 Introduction

Many multiscale problems are prohibitively expensive to solve due to scale disparity and high contrast. These problems are often solved using some type of reduced-order models. These include numerical homogenization [7, 20], multiscale finite element methods [13, 11, 12], heterogeneous multiscale methods [8], variational multiscale methods [15], mortar multiscale methods [19], localized orthogonal decomposition [16], and so on. The main idea behind local reduced-order model reduction techniques is to compute multiscale basis functions in each coarse block. These basis functions are computed using solutions of local problems.

Various approaches have been developed for designing multiscale basis functions. One of the earlier works [11, 13, 12, 2] use harmonic extensions of standard finite element basis functions in computing multiscale basis functions. Because of “homogeneous” traces along coarse boundaries, these approaches can have large errors due to the mismatch between the fine-grid solution and multiscale solutions along the edges of coarse blocks. These approaches have been generalized by using oversampling ideas [13, 12], where one uses larger regions and solve local problems. The solutions of these local problems are then used in constructing boundary conditions for multiscale basis functions. These approaches reduce the errors due to boundary conditions.

In later works [9, 10, 2], the authors showed that in the presence of high-contrast, one needs multiple basis functions. In [10, 2], Generalized Multiscale Finite Element is introduced, where the authors propose a systematic way of computing multiscale basis functions. The multiscale basis functions are computed by solving spectral problems in each coarse patch and selecting the dominant eigenvectors. In particular, the eigenvalues are ordered in increasing order and the eigenvectors corresponding to small eigenvalues are selected. The spectral convergence rate $1/\Lambda$ has been derived for these approaches, where $\Lambda$ is the smallest eigenvalue (across all coarse blocks) whose corresponding eigenvector is not included in the coarse space. In [4], using oversampling ideas and localization ideas [18, 16, 17, 14], the authors propose an approach which provides both mesh-size dependent convergence and spectral convergence. The main idea of this approach, called CEM-GMsFEM, is to (1) compute some GMsFEM basis (2) use constrained energy minimization in oversampling domains to construct multiscale basis functions. As a result, we have a minimal number of basis functions and can show $H/\Lambda^{1/2}$ convergence rate.

The above approaches can be classified as offline methods because the construction of multiscale basis functions does not take into account the right hand side. The offline methods can be tuned in various ways to achieve smaller errors; however, the error decay slows down as we add basis functions after a certain number of basis functions are selected. This slow down is due to some slow decay after certain eigenvalue. To improve this, in [3, 6, 1, 5], the authors propose an online approach. The main idea of online approaches is to add multiscale basis functions using the residual information after computing the coarse-grid solution. These online multiscale basis functions are computed adaptively and are chosen to decrease the error the most. They are solutions of local problems. Our analysis in [3, 6] shows that the error decay is proportional to $1-C\Lambda$, where $C$ is the constant (independent of scales and contrast) that guarantees the positivity of this quantity. This indicates that the error is not reduced unless $\Lambda$ is sufficiently away from [math], i.e., we have suffcient number of offline basis functions. This was demonstrated analytically and numerically in our papers [3, 6]. Since the online procedure can be costly, our goal is to perform only 1-2 iterations.

In this paper, we would like to investigate online approaches for CEM-GMsFEM and show that one can significantly improve the existing online approaches for some cases. In the paper, first we present an online approach, which differs from our previous approach since CEM-GMsFEM uses oversampling. In particular, the online basis functions are formulated in the oversampled regions. Secondly, we present an analysis of the proposed method. Our analysis shows that the error decay by adding online basis functions can be significantly better compared to $1-\Lambda$ in online GMsFEM. The error decay can be made close to [math] (i.e., we obtain very accurate approximation in one iteration) by choosing larger oversampling regions provided we have sufficient number of offline basis functions. To our best knowledge, this is a first result of this kind. Moreover, the online approaches can be made adaptive and adaptive error indicators can be derived.

We present numerical result. In our numerical results, we consider high-contrast permeability fields and place the source term in different locations. All results show that the error drops 3 orders of magnitude, which is much better compared to previous online GMsFEM. We also present numerical results using adaptivity, which shows that by selecting only some (few) regions, one can achieve a significant error decay. Moreover, our adaptive algorithm allows one to input a parameter which specifies a desired convergence rate.

The paper is organized as follows. In Section 2, we present some preliminaries. In Section 3, we present the construction of offline basis functions, and in Section 4, we present our online adaptive enrichment algorithm. Section 5 is devoted to numerical results. The analysis of our method is presented in Section 6. Finally, we present some concluding remarks in Section 7.

2 Preliminaries

In this paper, we consider a class of multiscale problems of the form

[TABLE]

subject to the homogeneous Dirichlet boundary condition $u=0$ on $\partial\Omega$, where $\Omega\subset\mathbb{R}^{d}$ is the computational domain. We assume that $\kappa(x)$ is a heterogeneous coefficient with multiple scales and very high contrast. In solving (1), it is desirable to construct multiscale basis functions that can be computed locally and give coarse-mesh convergence rate independent of the heterogenities and contrast. In [4], such approach has been developed. When the source term $f$ does not belong to $L^{2}(\Omega)$ and when one needs to obtain solutions with more refined accuracy, it is desirable to construct online basis functions that capture properties unresolvable by offline basis functions. It is the purpose of this paper to do this.

Next, we introduce the notion of fine and coarse grids. We let $\mathcal{T}^{H}$ be a usual conforming partition of the computational domain $\Omega$ into $N$ finite elements (triangles, quadrilaterals, tetrahedra, etc.), and let $H$ be the mesh size of $\mathcal{T}^{H}$. 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. The fine grid partition will be denoted by $\mathcal{T}^{h}$, and is by definition a refinement of the coarse grid $\mathcal{T}^{H}$. Here, we use $h$ to denote the fine mesh size of $\mathcal{T}^{h}$. In Figure 1, we give an illustration of the fine grid, coarse grid, and oversampling domain. In the figure, the coarse grid is contained by a union of rectangular coarse elements, denoted generically by $K$. Each coarse element is a union of finer rectangular elements. Moreover, for each coarse grid node $x_{i}$, we define $\omega_{i}$ as the union of coarse elements having the vertex $x_{i}$. We also define $\omega_{i}^{+}$ as an oversampled region for $\omega_{i}$. Finally, we define $N_{c}$ as the number of coarse grid vertices.

We let $V=H_{0}^{1}(\Omega)$. The weak solution $u\in V$ of the problem (1) satisfies

[TABLE]

where $a(u,v)=\int_{\Omega}\kappa\nabla u\cdot\nabla v$ and $(f,v)=\int_{\Omega}fv$. We would like to find a multiscale solution $u_{ms}$ in a subspace of $V$, denoted as $V_{ms}$, satisfying

[TABLE]

We will introduce the construction of the multiscale finite element space $V_{ms}$ in the next section. We remark that the multiscale space $V_{ms}$ consists of two components, which are the offline part and the online part. For the offline part, we will construct multiscale basis functions in the offline stage, that is, before solving the problem (3). Note that these basis functions are independent of the source term $f$. For the online part, we will construct multiscale basis functions in the online stage using the residual of an approximate solution. Note that, these basis functions depend on the source term $f$. In Section 3, we present the construction of offline basis functions. In Section 4, we present the construction of online basis functions and an online adaptive enrichment algorithm.

3 Offline basis functions

In this section, we present the construction of the offline multiscale finite element space. The construction of our offline basis functions follows the framework in [4]. To construct the multiscale space, we will first construct the auxiliary space by solving local spectral problem for each coarse element $K$. Next, we will construct the multiscale basis functions by solving some local minimization problems using the auxiliary basis functions. Our multiscale finite element space will then be the span of these multiscale basis functions. Next, we will discuss the construction of both auxiliary space and multiscale space in detail.

3.1 Auxiliary basis functions

Now, we present the construction of the auxiliary basis functions. For each coarse element $K_{i}$, we define $V(K_{i})=H^{1}(K_{i})$, and solve the spectral problem: find $(\lambda^{(i)}_{j},\phi_{j}^{(i)})\in\mathbb{R}\times V(K_{i})$

[TABLE]

where $a_{i}(u,v)=\int_{K_{i}}\kappa\nabla u\cdot\nabla v$, $s_{i}(u,v)=\int_{K_{i}}\tilde{\kappa}uv$, $\tilde{\kappa}=\kappa\sum_{j=1}^{N_{c}}|\nabla\chi_{j}|^{2}$ and $\{\chi_{j}\}$ is a set of partition of unity functions with respect to the coarse grid. We remark that one can take $\{\chi_{j}\}$ as the standard multiscale basis functions or the standard piecewise linear functions. We assume that the eigenfunctions satisfy the normalized condition $s_{i}(\phi_{j}^{(i)},\phi_{j}^{(i)})=1$. We can assume that the eigenvalues are sorted in ascending order, that is, $\lambda_{1}\leq\lambda_{2}\leq\dots$. We then choose the first $J_{i}$ eigenfunctions with small eigenvalue, and define the local auxiliary space $V_{aux}(K_{i})$ as the span of these eigen-basis functions, which is

[TABLE]

Notice that, by construction, $\lambda_{J_{i}+1}^{(i)}=O(1)$. We define

[TABLE]

Finally, the auxiliary space $V_{aux}$ is defined as the sum of these local auxiliary spaces, that is, $V_{aux}=\oplus_{i}V_{aux}(K_{i})$.

Using the auxiliary space, we define a projection operator $\pi:V\rightarrow V_{aux}$ by

[TABLE]

We also denote the kernel of the operator $\pi$ as

[TABLE]

3.2 Multiscale basis functions

Now we present the construction of offline multiscale basis functions using the auxiliary space. For each coarse element $K_{i}$, we define an oversampled region $K_{i}^{+}$ by extending $K_{i}$ by $\ell$ coarse grid layers. See Figure 2 for an illustration of $K^{+}$ with $\ell=2$.

For each $\phi_{j}^{(i)}\in V_{aux}$, we will construct a basis function $\psi^{(i)}_{j,ms}$ whose support is $K_{i}^{+}$. Using the results in [4], the multiscale basis function $\psi^{(i)}_{j,ms}$ is constructed by solving the following local problem: find $\psi^{(i)}_{j,ms}\in V_{0}(K^{+}_{i})$ such that

[TABLE]

where $V_{0}(K_{i}^{+})=H^{1}_{0}(K_{i}^{+})$ and $s(u,v)=\sum_{i=1}^{N}s_{i}(u,v)$. We remark that the above definition is defined in the continuous space $V_{0}(K_{i}^{+})$. In our numerical simulations, we solve the above problem using the fine mesh defined in $K_{i}^{+}$ and an appropriate finite element method. Finally, the multiscale finite element space $V_{ms}$ is defined as the span of these multiscale basis functions, namely, $V_{ms}=\text{span}\{\phi^{(i)}_{j,ms}\}$. This method is called the CEM-GMsFEM.

We remark that our multiscale basis functions are used to approximate the related global basis functions. The global basis function $\psi^{(i)}_{j}$ is defined by solving the following problem: find $\psi^{(i)}_{j}\in V$ such that

[TABLE]

The global space is defined by $V_{glo}=\text{span}\{\psi_{j}^{(i)}\}$. We note that these global basis functions have an exponential decay property (see [4]), which motivates the definitions of the multiscale basis functions $\psi^{(i)}_{j,ms}$ having local supports. One important property of the global space is the orthogonal decomposition $V=V_{glo}\oplus\tilde{V}$ with respect to the inner product $a(u,v)$. This global space will be used in the analysis of the convergence result of our online adaptive enrichment method.

4 Online basis functions and adaptive enrichment

In this section, we will introduce an online enrichment method for this CEM. We will first present the construction of online basis functions. Then, we will give an error estimate for using the online enrichment method. Here, for online basis function, we mean the basis functions constructed in online stage by using the residual of the solution which contain the information of the source term. We will construct the online basis function in an iterative process. We remark that the error will decay rapidly such that the error will within a acceptable range in the first or two iterations.

To begin, we define a residual functional $r:V\rightarrow\mathbb{R}$. Let $u_{ms}\in V_{ms}$ be a numerical solution computed by solving (3). The residual functional $r$ is defined by

[TABLE]

For each coarse neighborhood $\omega_{i}$, we define the local residual functional $r_{i}:V\rightarrow\mathbb{R}$ by

[TABLE]

The local residual $r_{i}$ gives a measure of the error $u-u_{ms}$ in the coarse neighborhood $\omega_{i}$.

The construction of online basis function is related to the local residual $r_{i}$. Using the local residual $r_{i}$, we will construct an online basis function $\phi_{on}^{(i)}$ whose support is an oversampled region $\omega_{i}^{+}$. More precisely, the online basis function $\beta_{ms}^{(i)}\in V_{0}(\omega^{+}_{i})$ is obtained by solving

[TABLE]

where $V_{0}(\omega^{+}_{i})=H^{1}_{0}(\omega_{i}^{+})$. We can perform the above construction for each $r_{i}$, or for some selected $r_{i}$ (with $i\in I$ for an index set $I$) based on an adaptive criterion. We remark that the above online basis is obtained in the local region $\omega_{i}^{+}$. This is the result of a localization result for the corresponding global online basis function $\beta_{glo}^{(i)}\in V$ defined by

[TABLE]

After constructing the online basis functions, we can enrich our multiscale space by adding these online basis to the multiscale space, namely, $V_{ms}=V_{ms}+\text{span}_{i\in I}\{\beta_{ms}^{(i)}\}$. Using this multiscale finite element space, we can compute a new numerical solution by solving the equation (3). We can repeat the process to enrich our multiscale space until the residual norm is smaller than a given tolerance. Next, we present the precise online adaptive enrichment algorithm.

Online adaptive enrichment algorithm

We first choose an initial space $V_{ms}^{(1)}$. This is the space obtained by using the offline basis functions constructed in Section 3. We also choose a real number $\theta$ such that $0\leq\theta<1$. This number determines how many online basis functions are needed in each online iteration. Then, we will generate a sequence of spaces $V_{ms}^{(m)}$ and a sequence of multiscale solutions $u_{\text{ms}}^{(m)}$ obtained by solving (3).

For each $m=1,2,\dots$, we assume that $V_{ms}^{(m)}$ is given. We will preform the following procedures to obtain the new multiscale space $V_{ms}^{(m+1)}$.

1. Step 1:

Find the multiscale solution in the space $V^{(m)}_{ms}$. That is, find $u_{ms}^{(m)}\in V^{(m)}_{ms}$ such that

[TABLE] 2. Step 2:

Compute the local residuals $z_{i}(v)$ where

[TABLE]

Define $\delta_{i}=\|z_{i}\|_{a^{*}}$ where $\|z_{i}\|_{a^{*}}=\sup_{v\in V_{0}(\omega_{i})}\cfrac{r(v)}{\|v\|_{a}}$. We re-numerate the indices of $\omega_{i}$ such that $\delta_{1}\geq\delta_{2}\geq\cdots$. Choose the first $k$ regions so that

[TABLE] 3. Step 3:

Compute the local online basis functions. For each $1\leq i\leq k$ and coarse region $\omega_{i}$, find $\beta^{(i)}_{ms}\in V_{0}(\omega_{i}^{+})$ such that

[TABLE]

where $r^{(m)}_{i}(v)=a(u^{(m)}_{ms},\chi_{i}v)-\int_{\Omega}f\chi_{i}v$. 4. Step 4:

Enrich the multiscale space. Let

[TABLE]

In the next section, we will analyze the convergence rate for this online adaptive enrichment method. In particular, we will prove the following theorem.

Theorem 1**.**

Let $u$ be the solution of (1) and let $\{u_{ms}^{(m)}\}$ be the sequence of multiscale solutions obtained by our online adaptive enrichment algorithm. Then we have

[TABLE]

where $E=3(1+\Lambda^{-1})(1+2(1+\Lambda^{\frac{1}{2}})^{-1})^{1-\ell}$, $M$ is maximum number of overlapping subdomains and $C$ is a constant.

Remark 1**.**

We note that the convergence rate depends two terms $C(\ell+1)^{d}E$ and $2M^{2}\theta$. By using enough number of oversampling layers, the term $C(\ell+1)^{d}E$ tends to zero. Thus, the factor $2M^{2}\theta$ dominates the convergence rate. One can choose $\theta$ to obtain a desired convergence rate. We will also confirm this by some numerical examples. This is an improvement over the online method in [3], where the convergence rate is $(C_{1}+C_{2}\theta)$ with $0<C_{1}<1$.

5 Numerical Result

In this section, we present some numerical results to demonstrate the convergence of our proposed method. We take the computational domain $\Omega=(0,1)^{2}$. The medium parameter $\kappa$ in the equation (1) is chosen to be the function shown in Figure 3. We note that the medium $\kappa$ contains high contrast inclusions and channels. The fine mesh size $h$ is taken to be $1/200$, while the coarse mesh size $H$ in this example is $1/10$. In all our results, we take the number of oversampling layers $\ell=2$. We will illustrate the performance of our method by using two different source terms $f_{1}=((x-0.5)^{2}+(y-0.5)^{2})^{-\frac{1}{4}}$ and $f_{2}=((x-0.5)^{2}+(y-0.5)^{2})^{-\frac{3}{4}}$. We will test the performance by considering uniform enrichment and by using the online adaptive enrichment algorithm presented in Section 4.

In Table 1, we present the $L_{2}$ error and the energy error for the case $f_{1}$ with uniform enrichment, that is $\theta=0$. From the first two online iterations, we observe very fast convergence of the method. Next, we will consider some adaptive results for this case. In Table 2, we present the error decay by using our online adaptive enrichment algorithm with $\theta=0.95$. That is, we only add basis for regions which account for the largest $5\%$ of the residual. From the table, we observe that the convergence rate in the energy norm is $0.9154$, which is close to $0.95$. This results confirm our assertion that the convergence rate can be controlled by the user-defined parameter $\theta$. We remark that the convergence rate is computed by taking the maximum of all $\|u-u_{ms}^{(m+1)}\|_{a}^{2}/\|u-u_{ms}^{(m)}\|_{a}^{2}$. In Table 3, we present the adaptive result with $\theta=0.1$. That is, we add basis for regions which account for the largest $90\%$ of the residual. From the table, we observe that the convergence rate in the energy norm is $0.0589$, which is close to $0.1$. This result also confirms our prediction. Moreover, we note that the adaptive approach allows adding a very few online basis functions to reduce the error to $1\%$.

Now, we consider the second source term $f_{2}$. In Table 4, we present the error decay using uniform enrichement. We observe very fast decay in error from this table. Next, we test the performance using adaptivity. In Tables 5 and 6, we present the error decays with $\theta=0.95$ and $\theta=0.1$ respectively. We observe that the convergence rates in these two cases are $0.9338$ and $0.09$ respectively. This confirms that the user-defined parameter is useful in controlling the convergence rate of our adaptive method. Moreover, we note that the adaptive approach allows adding a very few online basis functions to reduce the error to $1\%$. Furthermore, in Figure 4, we show the number of online basis functions added in the computational domain. For $\theta=0.95$, we will add a small number of basis functions in each iteration. We observe that the basis functions are added near the singularity of the source $f_{2}$ and along the high contrast channel. For $\theta=0.1$, more basis functions are added throughout the domain with a faster convergence rate. We still observe that more basis are added near the singularity of $f_{2}$ and along the high contrast channel in $\kappa$.

Finally, we present a test case with a more singular source term $f_{3}=-\nabla\cdot(\kappa\nabla(xy))$, shown in Figure 5, where the reference solution is also presented. In Table 7, we present the error decay with uniform enrichment and observe the same type of exponential decay as the earlier examples. We also observe that the error is relatively large when no online basis function is used. In Table 8, we present the results with the online adaptive enrichment algorithm with $\theta=0.1$. We see that the numerically computed convergence rate is $0.0771$, which is close to the parameter $\theta$.

6 Convergence analysis

In this section, we analyze the convergence of the online adaptive enrichment algorithm presented in Section 4. First, we need some notations. We will define two different norms for the 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}$ respectively.

Next, we will recall a few theoretical results from [4] that are useful for our analysis. The first result is Lemma 1.

Lemma 1**.**

There is a constant $D$ such that for all $v_{aux}\in V_{aux}$ there exists a function $v\in V$ such that

[TABLE]

The second result is a localization result, saying that the global basis function defined in (5) has an exponential decay outside an oversampled region. This result motivates the local multiscale basis functions defined in (4).

Lemma 2**.**

We consider the oversampled domain $K_{i}^{+}$ obtained from $K_{i}$ by extending $\ell$ coarse grid layers with $\ell\geq 2$. 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 (4) and let $\psi_{j}^{(i)}$ be the global multiscale basis functions obtained in (5). Then we have

[TABLE]

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

Next, we will need the following lemma in our analysis. The proof is given in the Appendix.

Lemma 3**.**

Assume the same conditions in Lemma 2. We have

[TABLE]

In the next lemma, we see that the same localization result holds for the online basis function defined in (6). The proof of the following lemma is the same as that for Lemma 2 and Lemma 3, and is therefore omitted.

Lemma 4**.**

We consider the oversampled domain $\omega_{i}^{+}$ obtained from $\omega_{i}$ by extending $\ell$ coarse grid layers with $\ell\geq 2$. We let $\beta_{ms}^{(i)}$ be the online multiscale basis functions obtained in (6) and let $\beta_{glo}^{(i)}$ be the global online multiscale basis functions obtained in (7). Then we have

[TABLE]

where E=3(1+\Lambda^{-1})\Big{(}1+(2(1+\Lambda^{-\frac{1}{2}}))^{-1}\Big{)}^{1-\ell}. Furthermore, we have

[TABLE]

Finally, we define a constant $C_{0}\in\mathbb{R}$ as

[TABLE]

We remark that $C_{0}\leq\max\{\tilde{\kappa}\}C^{2}_{p}$, where $C_{p}$ is the Poincare constant defined by $\|w\|_{L^{2}(\Omega)}\leq C_{p}\|\nabla w\|_{L^{2}(\Omega)}$ for $w\in H^{1}_{0}(\Omega)$.

Now we are ready to prove Theorem 1.

6.1 Proof of Theorem 1

First, by using the Galerkin orthogonality, we have

[TABLE]

The proof is based on a suitable choice of $w\in V_{ms}^{(m+1)}$, and consists of $4$ steps.

Step 1:

In this step, we will give a representation of the error $u-u_{ms}^{(m)}$. For each $i=1,2,\cdots,N_{c}$, we construct a global online basis function $\beta_{glo}^{(i)}\in V$ such that

[TABLE]

Summing over all $i=1,2,\cdots,N_{c}$, we have

[TABLE]

By the definition of $r_{i}^{(m)}(v)$, we have

[TABLE]

Therefore, we have

[TABLE]

From the above relation, we see that

[TABLE]

Using the decomposition $V=V_{glo}\oplus\tilde{V}$, we have

[TABLE]

Hence, we obtain the representation

[TABLE]

where $c_{j}^{(i)}$ are some coefficients. We will use this representation in the next steps. We remark that, in Step 2 and Step 3, we will localize the terms $\psi_{j}^{(i)}$ and $\beta_{glo}^{(i)}$, and estimate the errors.

Step 2:

In this step, we will localize each $\psi_{j}^{(i)}$ in (12) and estimate the error. In particular, we will estimate $\|\sum_{i=1}^{N}\sum_{j=1}^{J_{i}}c_{j}^{(i)}\psi_{j}^{(i)}-\sum_{i=1}^{N}\sum_{j=1}^{J_{i}}c_{j}^{(i)}\psi_{j,ms}^{(i)}\|_{a}$. We define $\eta:=u-u_{ms}^{(m)}+\sum_{i=1}^{N_{c}}\beta_{glo}^{(i)}$. Using (12), we have

[TABLE]

where the last equality follows from (5). We let $v_{aux}^{(i)}=\sum_{j=1}^{J_{i}}c_{j}^{(i)}\phi_{j}^{(i)}\in V_{aux}(K_{i})$. By Lemma 1, there exists a function $q^{(i)}\in V_{0}(K_{i})$ such that $\pi(q^{(i)})=v_{aux}^{(i)}$ and

[TABLE]

Taking $v=q^{(i)}$ in (14), we have

[TABLE]

Thus, by the orthogonality of the eigenfunctions $\phi_{j}^{(i)}$ and the normalization condition $s_{i}(\phi_{j}^{(i)},\phi_{j}^{(i)})=1$, we have

[TABLE]

where the last inequality follows from the definition of the constant $C_{0}$. Recalling the definition of $\eta$, we have

[TABLE]

where the last inequality follows from (11). Finally, using Lemma 2 and Lemma 3, we have

[TABLE]

Step 3:

In this step, we will derive an estimate for $\beta_{glo}^{(i)}-\beta_{ms}^{(i)}$. By Lemma 4, we have

[TABLE]

Using the equation (7),

[TABLE]

Since $\pi$ is an orthogonal projection onto the space spanned by the eigenfunctions $\{\phi_{j}^{(i)}\}$ for $j=1,2,\cdots,J_{i}$ and $i=1,2,\cdots,N$, we have

[TABLE]

Therefore, we have

[TABLE]

Step 4:

In this final step, we will prove the required convergence. Let $I=\{1,\cdots,k\}$. From the adaptive enrichment algorithm, we add the online basis functions $\beta^{(i)}_{ms}$ for $i\in I$. We will take $w$ in (10) as

[TABLE]

Using (10) and (12), we have

[TABLE]

Using Step 2 and Step 3 as well as Lemma 4, we see that

[TABLE]

Next, we will estimate the remaining term $\|\sum_{i\notin I}\beta_{glo}^{(i)}\|^{2}_{a}$. We write $p:=\sum_{i\notin I}\beta_{glo}^{(i)}$. Then, by (7),

[TABLE]

Thus, we have

[TABLE]

where $\theta$ is defined in (9). Lastly, we will estimate $\|z_{i}\|_{a^{*}}$. By definition,

[TABLE]

Thus, we have

[TABLE]

Combining the above equations, we have

[TABLE]

This completes the proof.

7 Conclusions

In this paper, we develop an online adaptive enrichment algorithm for CEM-GMsFEM. The CEM-GMsFEM, developed in [4], provides a general methodology for constructing multiscale basis functions that give a mesh-dependent convergence rate, regardless the contrast and heterogeneties of the media. In some applications, there is a need to further reduce the error without refining the coarse mesh. In these cases, one needs to add basis functions. We propose a strategy to compute new basis functions in the online stage using local residuals. These online basis functions provide very fast (exponential) decay in error. Moreover, an adaptive strategy is proposed to enrich the basis in some selected regions with large residuals. This strategy is determined by a user-defined parameter. We show that this used-defined parameter relates directly to the convergence rate of the method. That is, one can determine the convergence rate by using this parameter correspondingly. Several numerical tests are shown to validate our estimates.

Appendix A Proof of Lemma 3

In this appendix, we give a proof for Lemma 3. Let $K_{i}$ be a coarse element and let $K_{i,n}$ be the oversampled region by enlarging $K_{i}$ by $n$ coarse grid layers. For integers $m>n$, we denote $\chi^{m,n}_{i}$ as the cutoff function used in [4]. In particular, this function satisfies $\chi^{m,n}_{i}=1$ in $K_{i,n}$ and $\chi^{m,n}_{i}=0$ in $\Omega\backslash K_{i,m}$. Next, we set $w^{(i)}$ as $\sum_{j=1}^{J_{i}}c^{(i)}_{j}(\psi_{j}^{(i)}-\psi_{j,ms}^{(i)})$ and $w$ as $\sum_{i=1}^{N}w^{(i)}$. Then, by the definitions (4) and (5), we have

[TABLE]

and

[TABLE]

since $\text{supp}\{\sum_{j=1}^{J_{i}}c^{(i)}_{j}\phi^{(i)}_{j,ms}\}\subset\text{supp}\{\sum_{j=1}^{J_{i}}c^{(i)}_{j}\psi^{(i)}_{j,ms}\}\subset K_{i,\ell}$. Therefore, by subtracting the above two equations, we have

[TABLE]

Now, we estimate the right hand side of the above. Notice that,

[TABLE]

Summing over $i=1,2,\cdots,N$, we have

[TABLE]

Hence, we conclude that

[TABLE]

This completes the proof of the lemma.