A convergence analysis of Generalized Multiscale Finite Element Methods

TL;DR

This paper provides a comprehensive convergence analysis for the Generalized Multiscale Finite Element Method (GMsFEM), demonstrating its effectiveness for solving second-order elliptic equations with heterogeneous multiscale coefficients.

Contribution

The paper introduces a novel, general convergence analysis for GMsFEM that simplifies previous methods and extends to more sophisticated variants.

Findings

GMsFEM converges reliably for multiscale elliptic problems.

Error estimates depend on eigenvalues of unused eigenvectors.

Numerical experiments confirm theoretical convergence and stability.

Abstract

In this paper, we consider an approximation method, and a novel general analysis, for second-order elliptic differential equations with heterogeneous multiscale coefficients. We obtain convergence of the Generalized Multi-scale Finite Element Method (GMsFEM) method that uses local eigenvectors in its construction. The analysis presented here can be extended, without great difficulty, to more sophisticated GMsFEMs. For concreteness, the obtained error estimates generalize and simplify the convergence analysis of [J. Comput. Phys. 230 (2011), 937-955]. The GMsFEM method construct basis functions that are obtained by multiplication of (approximation of) local eigenvectors by partition of unity functions. Only important eigenvectors are used in the construction. The error estimates are general and are written in terms of the eigenvalues of the eigenvectors not used in the construction. The error analysis involve local and global norms that measure the decay of the expansion of the solution in terms of local eigenvectors. Numerical experiments are carried out to verify the feasibility of the approach with respect to the convergence and stability properties of the analysis in view of the good scientific computing practice.

Figures

Tables

Table 1. Table 1: Table of Notation

$\Omega$	pag.4	Definition domain of the elliptic problem
$H^1\left(\Omega \right)$	pag.9	Sobolev space of functions with continuous first derivatives on $\Omega$
$H_0^1\left(\Omega \right)$	pag.4	Sobolev space of functions with continuous first derivatives on $\Omega$ and vanishing on $\partial \Omega$
$w_i$	pag.5	Neighborhood of the node $y_i$ of a triangulation.
$w^K$	pag.5	Union of all neighborhood containing the element $K$
${\mu }_l^K$	pag.10	$l-th$ Dirichlet eigenvalue of the local eigenvalue problem on $K$
${\varphi }_l^K$	pag.10	$l-th$ Dirichlet eigenfunction of the local eigenvalue problem on $K$.
${\lambda }_l^{w_i}$	pag.12	$l-th$ Newmann eigenvalue of the local eigenvalue problem on $w_i$.
${\psi }_l^{w_i}$	pag.12	$l-th$ Newmann eigenfunction of the local eigenvalue problem on $w_i$.
${\Phi }_i$	pag.12	Generalased basis function associated to the neighborhood $w_i$
$𝒜$	pag.7	Elliptic operator on $\Omega$.
${𝒜}^K$	pag.11	Elliptic operator on the coarse block $K$.
${𝒜}^{w_i}$	pag.13	elliptic operator defined on the neighborhood $w_i$.
${𝒥}_L$	pag.9	Dirichlet projection operator truncated at $L$.
${𝒥}_L^K$	pag.13	Dirichlet projection operator on element $K$ truncated at $L$.
$I_L^{{\omega }_i}$	pag.18	Newmann projection operator in on neighborhood ${\omega }_i$ truncated at $L$.
$I_N$	pag.19	Generalized basis coarse interpolation operator
$J_D$	pag.19	Dirichlet basis coarse interpolation operator
$\tilde{V}\left({\omega }_i\right)$	pag.12	Space of functions in $H^1\left(w_i\right)$ which are zero on $\partial w_i\cap \partial \Omega$.
${\Vert \left|v\right|\Vert }_{s;D}$	pag.6	Norms based on the eigenvalue expansion decay on an set $D$.

Equations

$$-\text{mbox}div(\kappa (x)\nabla u)=f,$$

$$a(u,v)=f(v)\text{mbox}forallv\in H_0^1(\Omega ).$$

$$a(u,v)={\int }_{\Omega }\kappa (x)\nabla u(x)\nabla v(x)dx\text{mbox}forallu,v\in H_0^1(\Omega )$$

$$f(v)={\int }_{\Omega }f(x)v(x)dx\text{mbox}forallv\in H_0^1(\Omega ).$$

$${\omega }_i=\bigcup \{K\in {\mathcal{T}}^H;\text{leavevmode}\text{leavevmode}\text{leavevmode}{y}_i\in \overline{K}\}.$$

$${\omega }^K=\bigcup \{w_i\text{leavevmode}\text{leavevmode};\text{leavevmode}\text{leavevmode}\text{leavevmode}K\subset w_i\}.$$

$$a(u_0,v)=f(v),\text{for all}v\in \text{span}\{{\Phi }_i{\}}_{i=1}^{N_c}.$$

$$m(v,w)={\int }_{\Omega }\kappa vw\text{mbox}forallv,w\in H_0^1(\Omega ).$$

$$a(\varphi ,z)=\mu m(\varphi ,z)\text{mbox}forallz\in H_0^1(\Omega ).$$

$${\mu }_1\le {\mu }_2\le \cdots \le {\mu }_{\ell }\dots .$$

$$-\text{mbox}div(\kappa \nabla \varphi )=\mu \kappa \varphi$$

$$v=\sum _{\ell =1}^{\infty }m(v,{\varphi }_{\ell }){\varphi }_{\ell }$$

$$a(v,v)=\sum _{\ell =1}^{\infty }{\mu }_{\ell }m(v,{\varphi }_{\ell }){\varphi }_{\ell }$$

$$m(v,v)=\sum _{\ell =1}^{\infty }m(v,{\varphi }_{\ell }{)}^2.$$

$${\mu }_{\ell }m(u,\varphi )=a(u,\varphi )=f(\varphi ).$$

$$u=\sum _{\ell =1}^{\infty }\frac{1}{{\mu }_{\ell }}f({\varphi }_{\ell }){\varphi }_{\ell }.$$

$$|||v||{|}_{s;\Omega }^2=\sum _{\ell =1}^{\infty }{\mu }_{\ell }^{s}m(v,{\varphi }_{\ell }{)}^2.$$

$$|||v||{|}_{0;\Omega }^2=m(v,v)={\int }_{\Omega }\kappa v^2\text{mbox}and|||v||{|}_{1;\Omega }^2=a(v,v)={\int }_{\Omega }\kappa |\nabla v{|}^2.$$

$$|||u||{|}_{2;\Omega }^2=\sum _{\ell =1}^{\infty }{\mu }_{\ell }^{2}m(v,{\varphi }_{\ell }{)}^2.$$

$$\mathcal{A}u=\sum _{\ell =1}^{\infty }{\mu }_{\ell }m(u,{\varphi }_{\ell }){\varphi }_{\ell }.$$

$$|||\mathcal{A}u||{|}_{0;\Omega }^2=m(\mathcal{A}u,\mathcal{A}u)=\sum _{\ell =1}^{\infty }{\mu }_{\ell }^{2}m(u,{\varphi }_{\ell }{)}^2=|||u||{|}_{2;\Omega }^2<\infty .$$

$$a(u,v)=m(\mathcal{A}u,v)\text{mbox}forallv\in H_0^1(\Omega ).$$

$$\mathcal{A}u=-{\kappa }^{-1}\text{mbox}div(\kappa \nabla u).$$

$$|||u||{|}_{2;\Omega }^2=\Vert \mathcal{A}u{\Vert }_0^2={\int }_{\Omega }{\kappa }^{-1}|\text{mbox}div(\kappa \nabla u){|}^2.$$

$$u^N=\sum _{\ell =1}^{N}m(u,{\varphi }_{\ell }){\varphi }_{\ell }.$$

$${\kappa }^{-1}\text{mbox}div(\kappa \nabla u^N)=\sum _{\ell =1}^{N}m(u,{\varphi }_{\ell }){\mu }_{\ell }{\varphi }_{\ell }.$$

$$\Vert {\kappa }^{-1}\text{mbox}div(\kappa \nabla u^M)-{\kappa }^{-1}\text{mbox}div(\kappa \nabla u^N){\Vert }_0^2=\sum _{\ell =N+1}^{M}m(u,{\varphi }_{\ell }{)}^2{\mu }_{\ell }^2.$$

$${\int }_{\Omega }\kappa Uz=-\underset{N\to \infty }{\lim }{\int }_{\Omega }\kappa {\kappa }^{-1}\text{mbox}div(\kappa \nabla u^N)z={\int }_{\Omega }\kappa \nabla u\nabla z,$$

$${\int }_{\Omega }\kappa \mathcal{A}uz=m(\mathcal{A}u,v)=a(u,z)={\int }_{\Omega }\kappa \nabla u\nabla z.$$

Full text

††journal: Computer Physics Communications

1. 1 Introduction

2. 2 Preliminaries on multiscale finite element methods

3. 3 Global eigenvalue problem

4. 3.1 Global Norms based on eigenvalue expansion decay

5. 3.2 Approximation using global eigenvectors

6. 4 Dirichlet eigenvalue problem in coarse blocks

7. 5 Local Neumann eigenvalue problem in coarse neighborhoods

8. 6 GMsFEM space construction using local eigenvalue problems

9. 7 A technical assumption

10. 8 Approximation properties of the coarse space

11. 8.1 A coarse-scale interpolation operator

12. 8.2 Interpolation approximation

13. 9 Numerical experiments

14. 10 Discussions and final comments

1 Introduction

Approximation of partial differential equations posed on domains with multiscale and heterogeneous properties appear in variety of applications. For instance, when modeling subsurface flow scenarios, subsurface properties typically vary several orders of magnitude over multiple scales. In this case, the high-contrast in the properties such as permeability raises additional issues to be consider when constructing approximation of solutions. Several multiscale models to efficiently solve flow and transport processes have been considered. In despite of many contributions, the design and mathematical analysis of high-contrast multiscale problems continue being a challenging problem; See for instance [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. These approaches approximate the effects of the fine-scale features using a coarse mesh. They attempt to capture the fine scale effects on a coarse grid via localized basis functions. The main idea of Multiscale Finite Element Methods (MsFEMs) is to construct basis functions that are used to approximate the solution on a coarse grid. The accuracy of MsFEMs is found to be very sensitive to the particularities of the construction of the basis functions (e.g., boundary conditions of local problems). See for instance [13, 19, 20]).

It is known that the construction of the basis functions need to be carefully designed in order to obtain accurate coarse-scale approximations of the solution (e.g., [13]). In particular, the resulting basis functions need to have similar oscillatory behavior as the fine-scale solution. In classical multiscale methods, a number of approaches are proposed to construct basis functions, e.g., oversampling techniques or the use of limited global information (e.g., [1, 13]) that employs solutions in larger regions to reduce localization errors. Recently, a new and promising methodology was introduced for the construction of basis function. This methodology is referred as to Generalized Multiscale Finite Element Method (GMsFEM). The main goal of GMsFEMs is to construct coarse spaces for MsFEMs that result in accurate coarse-scale solutions. This methodology was first developed in [14, 15, 21] in connections with the robustness of domains decomposition iterative methods for solving the elliptic equation with heterogeneous coefficients subjected to appropriate boundary conditions

[TABLE]

where $\kappa(x)$ is a heterogeneous scalar field with high-contrast. In particular, it is assumed that $\kappa(x)\geq c_{0}>0$ (bounded below), while $\kappa(x)$ can have very large values.

A main ingredient in the construction was the use of local generalize eigenvalue problems and (possible multiscale) partition of unity functions to construct the coarse spaces. Besides using one coarse function per coarse node, in the GMsFEM it was proposed to use several multiscale basis functions per coarse node. These basis functions represent important features of the solution within a coarse-grid block and they are computed using eigenvectors of an eigenvalue problem. Then, in the works [22, 23, 24, 25], some studies of the coarse approximation properties of the GMsFEM were carried out. In these works and for applications to high-contrast problems, methodologies to keep small the dimension of the resulting coarse space were successfully proposed. The use of coarse spaces that somehow incorporates important modes of a (local) energy related to the problem motivated the general version of the GMsFEM. Thus,a more general and practical GMsFEM was then developed in [26] where several (more practical) options to compute important modes to be include in the coarse space was used. See also [27] for an earlier construction. It is important to mention that the methodology in [26] was designed for parametric and nonlinear problems and can be applied for variety of applications as it have been shown in recent developments not review here.

In this paper, we prove convergence of the GMsFEM method that uses local eigenvectors as developed in [23, 14, 15, 21]. The analysis presented here can be extended, without great difficulty, to more sophisticated GMsFEMs. Some convergence analysis of the GMsFEM, using local eigenvectors or approximation of them was obtained in [22]. The prove, as usual in finite element analysis, focuses on constructing interpolation operator to the coarse finite element space. The a priory error is obtained for square integrable right hand side $f$ in (2). Additionally, in [22] the authors make some assumptions concerning integrability of residuals and also concerning boundedness of the quotients of local energy norms with weight $\kappa$ and $\kappa\chi^{2}$ where $\chi^{2}$ is a especial partition of unity function. These assumptions are hard to verify in practice. Moreover, in the analysis they use a Caccioppoli inequality to write energy estimates from a region to a bigger region. Therefore, extensions of the analysis in [22] to other equations and/or different discretization is not straightforward.

In this paper, we substantially simplify the analysis of GMsFEM methods and remove the assumptions used in [22] to obtain convergence, yielding a general convergence proof and more suited for computational practice. We assume square integrability of the right hand side $f$. In order to obtain error bounds in terms of the decay of the eigenvalues used in the construction we assume that the problem is regular in the sense that the solution can be well approximated by local eigenvectors which in the case of smooth coefficients, square integrable right hand side and convex domains, is implied by the classical regularity of the problem.

It is worth to mention a main difference between the classical finite element analysis and the analysis of GMsFEM procedures for the case of heterogeneous multiscale coefficients. In the usual finite element analysis, to write the interpolation error estimates, it is assumed that the solution is smooth enough or regular enough in the classical (Sobolev) sense. This is done while using Hilbert norms (at least for elliptic problems). In the case of discontinuous multiscale coefficients, it is well know that solutions are not smooth in the classical sense. Then, the classical finite element analysis arguments do not work. In this paper, we are able to write interpolation error estimates using norms suitable for the problem at hand. In particular, to measure the “smoothness” of the solution we use the decay of the expansion of the solution in terms of global eigenvectors. This is motivated by the fact that, for a given elliptic operator, the eigenvectors are a good model for smooth functions in the scale of norms generated by powers of the operator. We define then global norms, using the decay of the expansion over global eigenvectors. We also define local norms using the decay of the expansion in terms of local eigenvectors (computed locally in a coarse node neighborhood). The main result of this paper is that we can compare the new local and global norms. With this new norms, we are able to write approximation results for the interpolation of functions that solve (2) with square integrable right hand side. We also prove error estimates in terms of the eigenvalues of the eigenvalue problem used in the construction.

The rest of the paper is organized as follow. In Section 2 we present some preliminaries on multiscale methods. In Section 3 we collect some facts on the global eigenvalue problem related to the problem. Here we introduce a scale of global norms used for the analysis. These norms measure the decay of the expansion in terms of global eigenvectors. In Sections 5 and 4 we study the local eigenvalue problems also using norms that measure the decay of the expansion in terms of the local eigenvectors. We also relate local norms to the boundary values of the eigenvalue problem. In Section 6 we review a very particular realization of the GMsFEM methodology that is the one analyzed in this paper. In Section 8 we obtain our interpolation error for the resulting method. We also write our convergence result. We present some numerical experiments in Section 9. Our numerical results verify our theoretical findings for smooth coefficients. We also consider a more practical case with heterogeneous multiscale coefficients. Finally, in Section 10 we draw some conclusions and make some final comments.

2 Preliminaries on multiscale finite element methods

In this section, we describe multiscale finite element method framework. In general terms, the MsFEMs compute the coarse-scale solution by using multiscale basis functions. It can be casted as a numerical upscaling procedure. Also as a numerical homogenization method where, instead of effective parameters representing small scale effects, basis functions are constructed that capture the small scale effects on solutions.

Multiscale techniques can be applied to variety of problems. In this paper, in order to fix ideas, we consider a second order elliptic problem with a possible multiscale high-contrast coefficient. More precisely, let $\Omega\subset\mathbb{R}^{2}$ (or $\mathbb{R}^{3}$) be a polygonal domain. We consider the elliptic equation with heterogeneous coefficients

[TABLE]

where $\kappa(x)$ is a heterogeneous scalar field with high-contrast. In particular, we assume that $\kappa(x)\geq c_{0}>0$ (bounded below), while $\kappa(x)$ can have very large values. We assume that $\kappa\ inL^{\infty}({\Omega})$ and therefore $\kappa$ might be discontinuous. The variational formulation of this problem is: Find $u\in H_{0}({\Omega})$ such that

[TABLE]

Here the bilinear form $a$ and the linear functional $f$ are defined by

[TABLE]

and

[TABLE]

Let $\mathcal{T}^{H}$ be a triangulation composed by elements $K$. We refer to the triangulation $\mathcal{T}^{H}$ as a coarse triangulation in the sense that does not necessarily resolve all the scales in the model (in our case that would be all variations and discontinuities of $\kappa$). We denote $\{y_{i}\}_{i=1}^{N_{v}}$ the vertices of the coarse mesh $\mathcal{T}^{H}$ and define the neighborhood of the node $y_{i}$ by

[TABLE]

and the neighborhood of an element $K$ by,

[TABLE]

Using the coarse mesh $\mathcal{T}^{H}$ we introduce coarse basis functions $\{\Phi_{i}\}_{i=1}^{N_{c}}$, where $N_{c}$ is the number of coarse basis functions. In our paper, the basis functions are supported in $\omega_{i}$; however, for $\omega_{i}$, there may be multiple basis functions. MsFEMs approximate the solution on a coarse grid as $u_{0}=\sum c_{i}\Phi_{i}$, where $c_{I}$ are determined from

[TABLE]

Once $c_{i}$’s are determined, one can define a fine-scale approximation of the solution by reconstructing via basis functions, $u_{0}=\sum_{i=1}^{N_{c}}c_{i}\Phi_{i}$.

3 Global eigenvalue problem

In this section, we recall some facts about the global eigenvalue problem associated to problem (1). We stress that the global eigenvalue problem is used in the analysis only and it is not use in the computations.

We start the presentation by introducing the global mass bilinear form. This is given by

[TABLE]

Note that we use the coefficient $\kappa$ in the mass matrix. The reason is that our main application in mind is on high-contrast problems and, as show in [25, 28, 14, 15], it is important to define the mass matrix with the coefficient $\kappa$. Moreover, more complicated bilinear forms can be also used as in recent developments in GMsFEM; see [26, 29].

We consider the eigenvalue problem (in weak form) that seeks to find eigenfunctions $\phi$ and scalars $\mu$ such that

[TABLE]

Denote it’s eigenvalues and eigenfunctions by $\{\mu_{\ell}\}$ and $\{\phi_{\ell}\}$, respectively. We order eigenvalues as

[TABLE]

We have $\mu_{1}>0$. The eigenvalue problem (4) is the weak form of the eigenvalue problem

[TABLE]

in ${\Omega}$ with homogeneous Dirichlet boundary condition on $\partial{\Omega}$.

We recall that the eigenvectors form a complete ($m$)-orthonormal system of $L^{2}({\Omega})$ that is also orthogonal with respect to the bilinear form $a$. Given any $v\in H^{1}_{0}(\Omega)$ we can write

[TABLE]

and compute the bilinear form $a$ as

[TABLE]

and the bilinear form $m$ as

[TABLE]

It is important to recall that the expansion of the solution can be explicitly given. In fact, from the weak form (2) we see that

[TABLE]

Then we have

[TABLE]

The eigenvector are the regular functions par excellence when working with the differential operator $-\mbox{div}(\kappa\nabla\cdot)$. In particular we stress the following fact.

Remark 1

We have that $\kappa\nabla\phi_{\ell}$ has a square integrable divergence. That is, $-\mbox{div}(\kappa\nabla\phi_{\ell})\in L^{2}({\Omega})$. This follows by observing that, in a generalize sense, $-\mbox{div}(\nabla\phi_{\ell})=\mu_{\ell}\kappa\phi_{\ell}$.

3.1 Global Norms based on eigenvalue expansion decay

In this section, we introduce a scale of norms that help measuring the decay of the expansions in terms of eigenvectors of the global eigenvalue problem. These norms are used in the a priori error estimates of our GMsFEM method. We note that, without assuming some sort of regularity of the solution of (2), it is difficult to measure the rate of the error in finite element approximations and give error estimates. For this paper, we only assume that the forcing term is square integrable in order to obtain approximation using global eigenvector. Later we consider the case of approximation using locally constructed basis functions with small support that employ local eigenvector in its design.

For any $v\in L^{2}({\Omega})$ written as $v=\sum_{\ell=1}^{\infty}m(v,\phi_{\ell})\phi_{\ell}$ and $s>0$, we introduce the norm $|||\cdot|||_{s;{\Omega}}$ defined by,

[TABLE]

We note that these norms depend on the bilinear forms $a$ and $m$ but, in order to make notation simpler, we do not stress this dependence in our notation. Note that

[TABLE]

In this paper, we mainly use the norm $|||\cdot|||_{s;{\Omega}}$ with $s=2$. We have that

[TABLE]

Then, if $|||u|||_{2;{\Omega}}<\infty$ we can define the operator $\mathcal{A}$ applied to $u$ by

[TABLE]

We readily have $\mathcal{A}u\in L^{2}({\Omega})$ since,

[TABLE]

Furthermore, if $|||u|||_{2;{\Omega}}<\infty$ we also have the following integration by parts relation, that can be verified by straightforward calculations,

[TABLE]

We now present a characterization of $\mathcal{A}$ using the divergence operator div and the coefficient $\kappa$. This implies that for even integer values of $s$ (in particular for $s=2$), the norm $|||\cdot|||_{s}$ is computed by subassembly of similar norms in subdomains.

Theorem 2

The operator $\mathcal{A}$ is a locally defined operator. More precisely, if $|||u|||_{2;{\Omega}}<\infty$ we have that $-\kappa^{-1}\mbox{div}(\kappa\nabla u)$ belongs to the space $L^{2}({\Omega})$ and we have

[TABLE]

Moreover, we have

[TABLE]

Proof. Recall that for $u\in L^{2}({\Omega})$, we have the expansion $u=\sum_{\ell=1}^{\infty}m(u,\phi_{\ell})\phi_{\ell}.$ For an integer $N$ define the truncated approximation of $u$ as,

[TABLE]

We construct the $-\kappa^{-1}\mbox{div}(\kappa\nabla u)$ as a limit in the $m-$norm of the sequence of rescaled divergences given by $\{\kappa^{-1}\mbox{div}(\kappa\nabla u^{N})\}_{N=1}^{\infty}$. To this end, we prove that the sequence $\{\kappa^{-1}\mbox{div}(\kappa\nabla u^{N})\}_{N=1}^{\infty}$ is a Cauchy sequence in the $m-$norm. Indeed, we have, by using the eigenvalue problem (6) and Remark 1, the following identity,

[TABLE]

So that, using the orthogonality of the eigenvectors, we conclude that for every $M>N$ we have,

[TABLE]

This implies the claim since the series $|||u|||^{2}_{2;{\Omega}}=\sum_{\ell=1}^{\infty}m(u,\phi_{\ell})^{2}\mu_{\ell}^{2}<\infty$. We conclude that there exist an $L^{2}({\Omega})$ function, denoted by $U$, such that we have $\|U+\kappa^{-1}\mbox{div}(\kappa\nabla u^{N})\|\to 0$ when $N\to\infty$. We also have that for any $z\in H^{1}_{0}({\Omega})$ it holds,

[TABLE]

which proves that $U=-\kappa^{-1}\mbox{div}(\kappa\nabla u^{N})z$.

Finally note that, by using (10), for every function $z\in H^{1}_{0}({\Omega})$ we have,

[TABLE]

Remark 3

Notice that if $u$ is the solution of (2) with $f\in L^{2}({\Omega})$, then, we have $\mathcal{A}u=\kappa^{-1}f$.

Lemma 4

Assume that $|||f|||_{s}^{2}<\infty$ and let $u$ be the solution of (2). Consider $t\geq 1$ such that $t-s-2\leq 0$ we have

[TABLE]

In particular, if $|||f|||_{0}<\infty$ we have that $|||u|||_{2}=|||\kappa^{-1}f|||_{0}=\int_{\Omega}\kappa^{-1}f^{2}$.

Proof. Using the explicit expansion in (9), the definition of the norm $\|\cdot\|_{t}$ and then increasing the order of eigenvalues we have

[TABLE]

This finishes the proof.

3.2 Approximation using global eigenvectors

In this section, we show how to obtain a priori error estimates if we use the space spanned by the first eigenvectors. Given an integer $L$ and $v\in H^{1}_{0}(\Omega)$, we define

[TABLE]

From (5), (7), and (8) it is easy to prove the following inequality

[TABLE]

When ${L}=1$ and $\kappa=1$ we obtain the usual Friedrichs’ inequality.

If $u$ is the solution of (2) and $|||f|||_{s}<\infty$ we have the following a priori estimate.

Lemma 5

Let $u$ be the solution of (2). If if $t-s-2<0$ we have

[TABLE]

Proof. Using the explicit expansion in (9), the definition of the norm $\|\cdot\|_{t}$ and then increasing oder of eigenvalues we have

[TABLE]

This finishes the proof.

We observe that, in particular, we have the following a priori error estimates.

[TABLE]

The space generated by the first eigenvalues gives good approximation spaces and the analysis becomes easy.

4 Dirichlet eigenvalue problem in coarse blocks

In this section, we study the local Dirichlet eigenvalue problem associated to problem (1). For any $K$, we define the following bilinear forms

[TABLE]

and

[TABLE]

We consider the eigenvalue problems that seek eigenfunctions $\phi\in H^{1}_{0}(K)$ and scalars $\mu$ such that

[TABLE]

and denote its eigenvalues and eigenvectors by $\{\mu_{\ell}^{K}\}$ and $\{\phi_{\ell}^{K}\}$, respectively. Note that the eigenvectors $\{\phi^{K}_{\ell}\}$ form an orthonormal basis of $L^{2}(K)$ with respect to the $m^{K}$ inner product. We order eigenvalues as

[TABLE]

The eigenvalue problem above corresponds to the approximation of the eigenvalue problem

[TABLE]

with homogeneous Dirichlet boundary condition on $\partial K$. These eigenvectors are the model of regular functions working with the differential operator $-\mbox{div}(\kappa\nabla\cdot)$. In particular, the operator is well defined and well behaved over these functions. We have that $\kappa\nabla\phi^{K}_{\ell}$ has a square integrable divergence. That is, $-\mbox{div}(\nabla\phi^{K}_{\ell})\in L^{2}(K)$. This follows by observing that, in a generalize sense, $-\mbox{div}(\nabla\phi_{\ell})=\mu_{\ell}\kappa\phi_{\ell}$.

Now we use the expansion in terms of local eigenvectors and define norms based on the decay of the eigenexpansion. These local norms are the ones that naturally appear in the local interpolation errors for our interpolation operator. A main issue is to compare this local norms with the global norms defined in Section 3.1 for $s>1$. For the case $s=2$, we prove in Theorem 6 that the local norms can be assembled to obtain a global norms equivalent to the norm $|||\cdot|||_{2,\Omega}$ only for functions that have zero value on block boundaries.

Given any $v\in L^{2}(K)$ we can write

[TABLE]

and compute the local energy bilinear form by

[TABLE]

We can also compute the local mass bilinear as,

[TABLE]

The local norm to measure the decay of the expansion is introduced as follows. We introduce the norm,

[TABLE]

Note that

[TABLE]

We consider the case $s=2$. If $|||u|||_{2,K}<\infty$ we can define the operator $\mathcal{A}^{K}$ by,

[TABLE]

which is square integrable since

[TABLE]

Additionally if $|||u|||_{2,K}<\infty$, we have the following local integration by parts relation that can be verified by direct calculations,

[TABLE]

We have the following result. This result reveals that the local norms $|||u|||_{2,K}$ is related to the $L^{2}$ integrability of $\mbox{div}(\kappa\nabla\cdot u)$. The proof of the following theorem follows the proof of Theorem 2 but we presented in the local setting in the interest of completeness.

Theorem 6

The operator $\mathcal{A}^{K}$ is a locally defined operator. More precisely, if $|||u|||^{2}_{2,K}<\infty$ we have that $-\kappa^{-1}\mbox{div}(\kappa\nabla u)$ belongs to the space $L^{2}(K)$ and we have

[TABLE]

Moreover, we have

[TABLE]

Proof. Since $u\in L^{2}(K)$ we have the expansion $u=\sum_{\ell=1}^{\infty}m^{K}(u,\phi_{\ell}^{K})\phi_{\ell}^{K}$. For any integer $N$, truncate this expansion to get,

[TABLE]

The sequence of rescaled divergences, $\{\kappa^{-1}\mbox{div}(\kappa\nabla u^{N})\}_{N=1}^{\infty}$, is a Cauchy sequence in the $m^{K}-$norm. Indeed, we have, by using the eigenvalue problem (12), the following identity,

[TABLE]

So that, using the orthogonality of the eigenvectors we conclude that for every $M>N$ we have,

[TABLE]

Then, there exists an $L^{2}(K)$ function, say $U$, such that $\|U+\kappa^{-1}\mbox{div}(\kappa\nabla\mu)\|_{0,{K}}\to 0$ when $N\to\infty$.

We also have that for any $z\in H^{1}_{0}(K)$ it holds,

[TABLE]

which proves that $U=-\kappa^{-1}\mbox{div}(\kappa\nabla\mu)$.

Finally note that for every function $z\in H^{1}_{0}(K)$ we have,

[TABLE]

Given an integer $L$ and $v\in H^{1}_{0}(K)$, we define

[TABLE]

From the analogous to (5), (7), and (8) it is easy to prove the following inequality

Lemma 7

Assume that $u\in H^{1}_{0}(K)$ and $|||u|||_{s,K}\leq\infty$ with $s>1$. We have for $1\leq s\leq t\leq 2$,

[TABLE]

In particular,

[TABLE]

5 Local Neumann eigenvalue problem in coarse neighborhoods

In this section, we study local eigenvalue problem associated to problem (1). For any $\omega_{i}$, we define the following bilinear forms

[TABLE]

and

[TABLE]

Define $\widetilde{V}(\omega_{i})=\{v\in H^{1}(\omega_{i})\leavevmode\nobreak\ :\leavevmode\nobreak\ v=0\mbox{ on }\partial\omega_{i}\cap\partial\Omega\}$ if $\partial\omega_{i}\cap\partial\Omega$ is non-empty and $\widetilde{V}(\omega_{i})=\{v\in H^{1}(\omega_{i})\leavevmode\nobreak\ :\leavevmode\nobreak\ \int_{\omega_{i}}v=0\}$ otherwise. We consider the eigenvalue problems that seek eigenfunctions $\psi\in\widetilde{V}(\omega_{i})$ and scalars $\lambda$ such that

[TABLE]

and denote its eigenvalues and eigenvectors by $\{\lambda_{\ell}^{\omega_{i}}\}$ and $\{\psi_{\ell}^{\omega_{i}}\}$, respectively. Note that the eigenvectors $\{\psi^{\omega_{i}}_{\ell}\}$ form an orthonormal basis of of $L^{2}(\omega_{i})$ with respect to the $m^{\omega_{i}}$ inner product. Note that $\lambda_{1}^{\omega_{i}}=0$ when $\partial\omega_{i}\cap\partial{\Omega}$ is empty, that is, when $\omega_{i}$ is a floating subdomain. We order eigenvalues as

[TABLE]

The eigenvalue problem above corresponds to the approximation of the eigenvalue problem

[TABLE]

with homogeneous Neumann boundary condition on $\partial\omega_{i}\cap\Omega$ and homogeneous Dirichlet boundary condition on $\partial\omega_{i}\cap\partial\Omega$ (when non-empty).

As mentioned before when studying the global eigenvalue problem, these eigenvectors are the model of regular functions working with the differential operator $-\mbox{div}(\kappa\nabla\cdot)$. In particular the operator is well defined and well behaved over these functions. We have that $\kappa\nabla\psi_{\ell}$ has a square integrable divergence. That is, $-\mbox{div}(\nabla\psi_{\ell})\in L^{2}(\omega_{i})$. This follows by observing that, in a generalize sense, $-\mbox{div}(\nabla\psi_{\ell})=\lambda_{\ell}\kappa\psi_{\ell}$.

Given any $v\in\widetilde{V}(\omega_{i})$ we can write

[TABLE]

and compute the local energy bilinear form by

[TABLE]

We can also compute the local mass bilinear as,

[TABLE]

The local norm to measure the decay of the expansion is introduced as follows. We introduce the semi-norm,

[TABLE]

Note that for $s=0$ we have a norm and for $s=1$ the semi-norms becomes a norm when restricted to non-constant functions on $\omega_{i}$, more precisely,

[TABLE]

We consider the case $s=2$. If $|||u|||_{2,\omega_{i}}<\infty$ we can define the operator $\mathcal{A}^{\omega_{i}}$ by,

[TABLE]

which is square integrable since

[TABLE]

Additionally, if $|||u|||_{2,\omega_{i}}<\infty$, we have the following local integration by parts relation (that can be verified directly by the series expansion of both sides),

[TABLE]

We have the following result.

Theorem 8

The operator $\mathcal{A}^{\omega_{i}}$ is a locally defined operator. More precisely, if $|||u|||^{2}_{2,\omega_{i}}<\infty$ we have that $-\kappa^{-1}\mbox{div}(\kappa\nabla u)$ belongs to the space $L^{2}(\omega_{i})$ and we have

[TABLE]

Moreover, we have

[TABLE]

Proof. Since $u\in L^{2}(\omega_{i})$ we have the expansion $u=\sum_{\ell=1}^{\infty}m^{\omega_{i}}(u,\psi_{\ell}^{\omega_{i}})\psi_{\ell}^{\omega_{i}}$. For any integer $N$, truncate this expansion to get,

[TABLE]

The sequence of rescaled divergences, $\{\kappa^{-1}\mbox{div}(\kappa\nabla u^{N})\}_{N=1}^{\infty}$, is a Cauchy sequence in the $m^{\omega_{i}}-$norm. Indeed, we have, by using the eigenvalue problem (13), the following identity,

[TABLE]

So that, using the orthogonality of the eigenvectors we conclude that for every $M>N$ we have,

[TABLE]

Then, there exists an $L^{2}(\omega_{i})$ function, say $U$, such that $\|U+\kappa^{-1}\mbox{div}(\kappa\nabla u_{N})\|_{0,{\omega_{i}}}\to 0$ when $N\to\infty$. We also have that for any $z\in H^{1}_{0}(\omega_{i})$ it holds,

[TABLE]

which proves that $U=-\kappa^{-1}\mbox{div}(\kappa\nabla u_{N})$. Finally, note that for every function $z\in H^{1}_{0}(\omega_{i})$ we have,

[TABLE]

Remark 9

In virtue of Theorem 8 and the equality (14) we see that a necessary condition for the $L^{2}$ integrability of $\mbox{div}(\kappa\nabla u)$ is that $\kappa\partial_{\eta}u=0$ on $\partial\omega_{i}$. Note that in general, if we are not sure $|||u|||^{2}_{2,\omega_{i}}<\infty$ and we do not assume $\partial_{\eta}u=0$ then, the integration by parts become

[TABLE]

Therefore, doing estimates about the eigenvalue decay is harder in this case. We also mention that in the analysis presented in [22], it is assume the square integrability of $\mbox{div}(\kappa\nabla u)$ that, as mentioned above, implies $\kappa\partial_{\eta}u=0$. Assuming that the solution has zero flux across boundaries neighborhoods is not a general assumption. A main contribution of this paper is to clarify this main assumption of [22] and to present an analysis valid for the general case where the solution $u$ does not have null fluxes across neighborhood boundaries. As we show in our numerical experiments, for the case of a solution that is not close to a function with null flux across neighborhood boundaries, the convergence rate of the GMsFEM as introduced in [22] is not optimal and additional basis functions constructed from local Dirichlet eigenvalues must be introduced to recover good convergence.

Lemma 10

Assume that $u\in H^{1}(\omega_{i})$, $\partial_{n}u=0$ and $|||u|||_{s,\omega_{i}}\leq\infty$ with $s>1$. We have for $1\leq s\leq t\leq 2$,

[TABLE]

In particular,

[TABLE]

6 GMsFEM space construction using local eigenvalue problems

In this section, we summarize the construction of coarse scale finite element spaces using a GMsFEM framework. In order to focus in the analysis of convergence we consider a particular case of the construction of spaces using the GMsFEM framework as introduced in [22, 28]. This construction evolved to the GMsFEM method as described in [26]. The method presented in this paper to obtain convergence can be also carried out for the constructions in [26] under appropriate assumptions of the local spectral problems used for the construction of coarse spaces.

We choose the basis functions that span the eigenfunctions corresponding to small eigenvalues. We note that $\{\omega_{i}\}_{y_{i}\in\mathcal{T}^{H}}$ is a covering of $\Omega$. Let $\{\chi_{i}\}_{i=1}^{N_{v}}$ be a partition of unity subordinated to the covering $\{\omega_{i}\}$ such that $\chi_{i}\in V^{h}(\Omega)$ and $|\nabla\chi_{i}|\leq\frac{1}{H}$, $i=1,\dots,N_{v}$. Define the set of coarse basis functions

[TABLE]

where $L_{i}$ is the number of eigenvalues that will be chosen for the node $i$; see [30, 31] for more details on the generalized finite element method using partitions of unity. Denote by $V_{0}$, as before, the local spectral multiscale space

[TABLE]

Define also

[TABLE]

Finally define,

[TABLE]

In practice, the computation of the multiscale basis functions have to be done. For instance, the computation of the multiscale basis functions can be performed in a fine grid (local to each region) that is sufficiently fine to resolve and represent the scales of the problem. In our case, this means that the fine-grid have to be sufficiently fine to represent the variations and discontinuities of the coefficient $\kappa$. The computations of the basis functions are local to each coarse region and can be done in a preprocessing step (taking advantage of parallel computations). For more details and related concepts, we refer to [26].

We define $u_{H}$ as the Galerkin approximation using the space $V_{0}$, that is,

[TABLE]

7 A technical assumption

In order to get convergence rates we assume that we can decompose the exact solution $u$ as $u\approx u_{D}+u_{N}$ where the $u_{D}$ can be approximated using the space $V_{D}$ and $u_{N}$ can be approximated using the space $V_{D}$. This requires that $u_{N}$ and $u_{D}$ have the right boundary condition on the coarse block edges. See Remark 9. More precisely we estate the following assumption.

Assumption 11

Let $u$ be the exact solution, that is $-\mbox{div}(\kappa\nabla u)=f$. We assume that there exists $u_{D}$, $u_{N}$ and $\epsilon\preceq H$ such that

We have

[TABLE] 2. 2.

We have the boundary data given by

[TABLE]

and

[TABLE] 3. 3.

We have the bounds,

[TABLE]

Note that 2. in the case $\epsilon=0$ implies that $u_{N}=u$ on $\partial K$ and that $\partial_{\eta}u_{D}=\partial_{\eta}u$ on $\partial K$ and therefore we should be able to split the Neumann and Dirichlet boundary data into different functions with regular divergece. This allows us to approximate each part by the rightly constructed subspace with eigenvalues with Nuemann and Dirichlet data. We can think of Assumption 11 as a natural extension (to the case of variable coefficient) of a regularity assumption. If fact we can give a following example for the case of regular coefficient.

Remark 12

In the case of regular coefficient $\kappa$ and regular right hand side $f$ it is known that, given any $\epsilon>0$ we can approximate the solution $u$ by $C^{1}$ finite elements defined on a sufficiently fine triangulation with square $H^{1}$-error smaller than $\epsilon^{2}\int_{\Omega}f^{2}$ so we can get the Assumption 11. In fact, for the case of $C^{1}$ (Hermite) finite element spaces defined on a triangulation we can relate $u_{D}$ to the derivative value degrees of freedom while $u_{N}$ will correspond to the nodal values of the function $u$.

Remark 13

In the case of regular coefficient, say $\kappa=1$, with regular right hand side it is known that the solution $u$ is regular $u\in H^{2}({\Omega})$. In this case, using standard regularity results we can show that Assumption 11 holds with $\epsilon=0$. We can construct $u_{D}$ and $u_{N}$ by solving a forth order problem. In fact, consider

[TABLE]

Define the global function $u_{D}$ by $u_{D}|_{K}=u_{K,D}$. Note that $u_{D}\in H^{1}_{0}(\Omega)$. Define also $u_{N}=u-u_{D}$. We have

[TABLE]

Note also that $u_{D}+u_{N}=u$.

8 Approximation properties of the coarse space

We mention that, in the presence of high-contrast multiscale coefficient $\kappa$, if $L$ is large enough then $\lambda_{L+1}^{\omega_{i}}$ is contrast independent and in this case we refer to (11) as a contrast independent weighted Poincaré inequality. The basis function encode information of the behavior of solutions due to the high-contrast in the multiscale coefficient and then allow us to compute using a coarse grid size $H$that does not need to resolve all discontinuities of the coefficient $\kappa$. This is a main motivation for the construction of the space presented above. For many recent developments using these ideas we refer the interested reader to [26] and references there in. In this paper, in order to focus in the convergence analysis, we work with the coarse space presented above. We also note that, more involved and sophisticated coarse space can be constructed as in [26]. The ideas developed in this paper also apply to variety of cases proposed in [26].

8.1 A coarse-scale interpolation

operator

Given an integer $L$, and $v\in V^{h}(\Omega)$, we define

[TABLE]

From Lemma 10 the following inequality holds for $1\leq s\leq 2$,

[TABLE]

and, if $\partial_{\eta}v=0$ on $\partial\omega_{i}$, we have

[TABLE]

Recall that we assume that $\partial_{\eta}v=0$ on $\partial\omega_{i}$, otherwise the term $\sum_{\ell=1}^{\infty}\int_{\partial\omega_{i}}\partial_{\eta}v\psi^{\omega_{i}}_{\ell}$ will come on the right hand side. These last term is harder to bound.

Define the coarse interpolation $I_{N}$ by

[TABLE]

where $I^{\omega_{i}}_{L_{i}}$ is defined in (19). Note that we have

[TABLE]

This interpolation was analyzed in [14, 15] for high-contrast problems, there, it was used to obtain a robust two level domain decomposition method and no approximation in energy norm was needed. Later, in [22] an analysis to obtain approximation results in $H^{1}$ norm was carried out. The analysis was rather complex and difficult to extend to other applications. The analysis we present in this paper simplifies that of [22]. Indeed and it is easier to extend and to combine with different techniques to analyze different realizations of the GMsFEM methodology.

In order to avoid the assumption of square integrable residuals in the approximation result as it is done in [22] (which implies assuming that the solution has zero flux across neighborhood boundaries - see Remark 9) we introduce an additional interpolation operator into $V_{D}$. Given $L_{K}$ and $v$ define

[TABLE]

and

[TABLE]

where we extend $J_{L_{K}}^{K}v$ by zero outside the block $K_{j}$.

8.2 Interpolation approximation

We note that the analysis presented in this section is closely related to the analysis in [14, 15, 22]. In particular we simplify and the analysis presented in [22].

Lemma 14

Consider $v\in H^{1}_{0}(\Omega)$. We have the following weighted $L^{2}$ approximation

[TABLE]

where $\lambda_{K,L+1}=\min_{y_{i}\in K}\lambda_{L_{i}+1}^{\omega_{i}}$ and therefore $\int_{\Omega}\kappa(v-I_{0}v)^{2}\preceq\frac{1}{\lambda_{L+1}}||v||_{H^{1}(\Omega)}^{2}$ where $\lambda_{L+1}=\min_{K}\lambda_{K,L+1}$.

Proof. First we prove (22). Using that $\chi_{i}\leq 1$ we have

[TABLE]

and using (20) to estimate the last term above, we obtain the result.

We now present the result in the $H^{1}$-norm. Here we are more explicit in the assumption than in the analogous result in [22] where they assume that in each neighborhood the residual is square integrable. See [22]. The proof is analogous to the one presented in [22] and we presented here for completeness.

Lemma 15

Assume that $f=\mbox{div}(\kappa\nabla u)\in L^{2}(w^{K})$ and also assume that for each $i$, $y_{i}\in K$ we have $\partial_{\eta}u=0$ on $\partial\omega_{i}$. Then, the following energy approximation holds,

[TABLE]

where $\lambda_{K,L+1}=\min_{y_{i}\in K}\lambda_{L_{i}+1}^{\omega_{i}}$ and $\omega^{K}$ is defined in (3).

Proof. We note that $\sum_{y_{i}\in K}\nabla\chi_{i}=0$ in $K$, and then we can fix $y_{j}\in K$ and write $\nabla\chi_{j}=-\sum_{y_{i}\in K\setminus\{y_{j}\}}\nabla\chi_{i}$. We obtain,

[TABLE]

which gives the following bound valid on $K$,

[TABLE]

From (24) we get

[TABLE]

To bound the first term above we use (20) as follows,

[TABLE]

The second term in (25) is estimated using (21)

[TABLE]

By combining (26), (27) and (25) we obtain (23).

A similar lemma for the case of Dirichlet boundary conditions on $K$ is presented next. This is direct consequence of Lemma 7.

Lemma 16

Assume that $f=\mbox{div}(\kappa\nabla u)\in L^{2}(\Omega)$ and also assume that for each $K\in\mathcal{T}^{H}$, we have $u=0$ on $\partial K$. Then, the following energy approximation holds,

[TABLE]

The idea is to use the Assumption 11 and then to apply $I_{N}$ to $u_{N}$ and the other part, that is $u_{D}$ will be approximated by a truncated expansion on $V_{D}$. Under Assumption 11 and for the solution $u$ define the coarse interpolation $I_{0}$ by

[TABLE]

Recall that $L_{i}$ is the number of Newmann eigenfunctions considered in the neighborhood $w_{i}$ and $L_{K}$ is the number of Dirichlet eigenfunctions considered on the element $K$. We finally present our main approximation result.

Theorem 17

Assume that $|||u|||_{2,\Omega}<\infty$ where $u$ is the solution of (1) and that Assumption 11 holds, the following approximation for the energy interpolation error holds,

[TABLE]

where $\lambda_{L+1}=\min_{K}\lambda_{K,L+1}^{\omega_{i}}$ and $\mu_{L+1}=\min_{K}\mu_{L+1}^{K}$.

Proof. By definition (28), our technical assumption and using the triangular inequality we have

[TABLE]

Using Lemma 15 and noting that $\bigcup K\subset\bigcup w^{K}$ we have for the firs term in (30)

[TABLE]

Now using Lemma 16 for the second term in (30) we get

[TABLE]

using (17) we obtain the bounds for the third term in (30)

[TABLE]

with (33) and using (18) in (31) and (32) we obtain (29) from (30).

With the tools we have at hand we can obtain the convergence of the GMsFEM method of Section 6. Under Assumption 11, by combining the Cea’s lemma with our interpolation approximation result (Lemma 17 ) and the estimates in Lemma 4 we obtain the following error estimates.

Theorem 18

Let $u$ be the solution of problem (2) with $f$ being square integrable and let $u_{H}$ the solution of (16) using the coarse basis functions constructed in (15). Suppose Assumption 11 holds. We have,

[TABLE]

where ${\lambda_{L+1}}$ is the minimum left-out eigenvalue and was introduced in Theorem 17

It is easy to see that if we map the local eigenvalue problem posed in $\omega_{i}$ (of diameter $H$) to a size one domain, then, the resulting eigenvalues scale with $H^{-2}$. If we re-scale all the eigenvalue problem to size one domains, we can then write the estimates in terms of eigenvalue problems posed in one size domains. In this way it is clear the $H$ dependence of our estimate. We have the following result.

Corollary 19

Under the assumptions of Theorem 18 we have,

[TABLE]

where the hidden constant involves eigenvalues of re-scaled eigenvalue problems posed on the unit square.

9 Numerical experiments

Our aim is to show that, when applying GMsFEM, including Dirichlet’s basis functions is necessary in some cases. So let’s consider problem (1) on a square domain $\Omega$ subject to homogeneous Dirichlet’s boundary conditions. Define a fine rectangular $512\times 512$ mesh $\mathcal{T}^{h}$ and a coarse rectangular $4\times 4$ mesh $\mathcal{T}^{H}$. Coarse basis functions are associated to each nodal point of $\mathcal{T}^{H}$ and supported on it’s $4$ adjacent rectangles. We apply the GMsFEM method to problem (1) with homogeneous medium ($\kappa\equiv 1$) considering $2$ different smooth sources.

Experiment $1$: Apply GMsFEM to problem (1) with source $f_{1}$ being a smooth function composed by product of polynomials and sines which correspond to the exact solution $u_{1}$

[TABLE]

Experiment $2$: Apply GMsFEM to problem (1) with source $f_{2}$ being a linear combination of sines that correspond to the exact solution $u_{2}$ also constructed by the linear combination of the $\mathbb{R}^{2}$ tonsorial product of functions $\{\mbox{sin}(4k\pi)\}_{k=1,2}$

[TABLE]

with $c_{ij}=1/(i+2(j-1))$ decreasing as $i$ and $j$ grow.

We vary the number of Newman function $N_{N}=1\cdots 40$ while fixing the number of Dirichlet function by element ($N_{D}=0\cdots 20$ in Experiment $1$ and $N_{D}=0\cdots 3$ in Experiment $2$), . The errors in the energy norm are presented in Figure 1 while in Figures 2 and 3 we can see the evolution of solution $u_{1}$ and $u_{2}$ respectively.

In the first example, including Dirichlet functions doesn’t have much impact in the approximation, the error in the energy norm decreases faster by including more Newmann functions. While in the second example by including Dirichlet functions the energy norm decreases faster than by including Newman functions as we can see in Figure 1

Now let us consider two different heterogeneous fields. The first is shown is composed by 3 chanels of high permeability in $64\times 64$ fictitious geological mesh shown in Figure 4. The second heterogeneous medium to be consider will be the last $64\times 64$ block of the geological permeability $SPE10$ porous medium taken from [32] (see Figure 5). This is a widely used heterogeneous porous medium for simulations (see for example [33])

Next we perform two numerical experiments to evaluate the performance of including Dirichlet basis when heterogeneity of the medium is present.

Experiment $3$: Let us consider a source taking alternating values $1$ and $-1$ on the $4\times 4$ coarse mesh $\mathcal{T}^{H}$ of $\Omega$

[TABLE]

with $K_{j}\in\mathcal{T}^{H}$ and $j\in\{1,..,16\}$. The coarse mesh basis are computed on the fine $512\times 512$ mesh $\mathcal{T}^{h}$ and approximated reference solution is computed in a $1024\times 1024$ finner mesh using classical finite element $\mathcal{Q}^{1}$. As we did in Examples $1$ and $2$ we fix the number of Dirichlet coarse basis in the expansion $N_{D}$ from $1$ to $20$ and vary the number of Newmann coarse basis $N_{N}$ from $1$ to $40$. We use the $3$-channels heterogeneous medium shown in Figure 4

Experiment $4$: Exactly the same parameters in Experiment $3$ except for the medium which will be taken from SPE10 as previously described and shown in Figure 5.

Both for Experiment $3$ and $4$ the approximated solution benefits from including Dirichlet basis in the expansion.

10 Discussions and final comments

In this paper, we obtained error estimates for GMsFEM approximation of high-contrast multiscale problems. This construction uses local Neumann eigenvectors on neighborhoods and Dirichlet eigenvectors on elements to construct finite element basis function. The analysis is based on eigenfunction expansions and the norms used for the error estimates measure the decay of the expansion of the solution in terms of local eigenfunctions. The norms in the interpolation error estimates can be bounded by the $L^{2}$ norm of a rescaling of the forcing term. For the analysis we assume that the solution can be approximated by a sum of two functions, one with zero flux across coarse blocks boundaries, and the other with zero value on coarse blocks boundaries. This assumption is easily verified for classical regular problems. The introduction of this assumption allowed us to extend and simplify the convergence analysis presented in [22]. The error estimates derived here can be applied to several situations that are under current research. For instance, it is possible to obtain error estimates for general bilinear forms using the analysis presented here combined with the construction of coarse spaces in [34].

Notation

In general we use the letters $\mu$ and $\phi$ to refer to Dirichlet eigenvalues and eigenvectors; and we use letters $\lambda$ and $\psi$ to refer to Newmann eigenvalues and eigenvectors. Letter $\Phi$ is reserved to generalized basis.

Acknowledgments

Eduardo Abreu thanks the FAPESP for support under grant 2016/23374-1. Juan Galvis wants to thank KAUST hospitality where part of this work was developed and also the discussion on coarse space approximations properties and related topics with several colleagues, among them, Joerg Willems, Marcus Sarkis, Raytcho Lazarov, Jhonny Guzmán, Chia-Chieh Chu, Florian Maris and Yalchin Efendiev.