An analysis of the NLMC upscaling method for high contrast problems

# An analysis of the NLMC upscaling method for high contrast problems

Lina Zhao111Department of Mathematics,The Chinese University of Hong Kong, Hong Kong Special Administrative Region. (lzhao@math.cuhk.edu.hk)   Eric T. Chung222Department of Mathematics,The Chinese University of Hong Kong, Hong Kong Special Administrative Region. (tschung@math.cuhk.edu.hk)

Abstract: In this paper we propose simple multiscale basis functions with constraint energy minimization to solve elliptic problems with high contrast medium. Our methodology is based on the recently developed non-local multicontinuum method (NLMC). The main ingredient of the method is the construction of suitable local basis functions with the capability of capturing multiscale features and non-local effects. In our method, each coarse block is decomposed into various regions according to the contrast ratio, and we require that the contrast ratio should be relatively small within each region. The basis functions are constructed by solving a local problem defined on the oversampling domains and they have mean value one on the chosen region and zero mean otherwise. Numerical analysis shows that the resulting basis functions can be localizable and have a decay property. The convergence of the multiscale solution is also proved. Finally, some numerical experiments are carried out to illustrate the performances of the proposed method. They show that the proposed method can solve problem with high contrast medium efficiently. In particular, if the oversampling size is large enough, then we can achieve the desired error.

Keywords: Contraint energy minimization, Upscaling, Non-local multicontinuum method, High contrast

## 1 Introduction

In this paper we consider

 −∇⋅(κ∇u)=finΩ,u=0on∂Ω, (1.1)

where is the computational domain and is a high contrast with and is a multiscale field. The proposed method can be extended to 3D easily.

If the coefficient is rough, then the solution to (1.1) will also be rough; to be specific, will not in general be in and may not be in for any . For this kind of low regularity, standard analysis usually fails. Moreover, the classical polynomial based finite element methods could perform arbitrary badly for such problems, see, e.g., [4]. To resolve this issue, various numerical methods have been proposed and analyzed, and among all the methods we mention in particular the special finite element methods [2, 3], the upscaled models [12, 30] and the multiscale methods [20, 18, 21, 19, 1, 10, 7, 6, 23, 24, 16, 17].

The concept of non-local upscaling has been successfully applied to problems in porous media, see, e.g., [14, 11, 13]. Motivated by the work given in [15], the nonlocal multicontinua (NLMC) upscaling technique was initially introduced for flows in heterogeneous fractured media in [9], and have been successfully applied to different problems under application [25, 26, 27, 28]. The main idea of NLMC upscaling technique is to construct the multiscale basis functions over the oversampling domain via an energy minimization principle. Note that the constraint should be chosen properly in order to make the localization possible. One distinctive feature of the method is that it allows a systematic upscaling for processes in the fractured porous media, and provides an effective coarse scale model whose degrees of freedom have physical meaning.

Inspired by the work given in [9, 25, 26], the goal of this paper is to extend the idea of nonlocal multicontinua to problem (1.1). For our approach, we start with decomposing the coarse block into different regions and the criterion used for the decomposition is to have relatively small contrast ratio within each region. Then, we define the constraint energy minimzation problem in the oversampling domain, where the restriction for the basis functions is defined such that they have mean value one in the chosen region and zero mean otherwise, in addition the basis functions vanish on the boundary of the oversampling domain. We remark that the vanishing property is important for the localization of the multiscale basis functions and the localization idea has also been exploited in [22] to solve problems with heterogeneous and highly varying coefficients. Next, we can solve the local minimization problem by using the equivalent saddle point formulation to achieve the multiscale basis functions. The resulting multiscale basis functions have decay property, in addition, it can capture the fine-grid information well provided proper number of overampling layers are chosen. With the multiscale basis functions, we can solve the upscaled equation to obtain the upscaled coarse grid solution. It is worth mentioning that in our method the number of basis function is relatively small and it is equal to the number of scales over the domain. We also analyze the convergence of the proposed method. For this, we first compare the difference between the multiscale basis functions and the global basis functions, combining this with the convergence of the global solution, then we can prove the convergence of the multiscale solution in norm and weighted energy norm. The analysis indicates that the convergence rate only depends on the local contrast ratio, namely, the contrast ratio within each region. With proper number of oversampling layers, the first order convergence measured in energy norm can be obtained. Some numerical experiments are also carried out. The numerical experiments show that with the fixed coarse mesh size, the oversampling layers should be selected properly to achieve the desired error, in addition, for a fixed oversampling size, the performance of the scheme will deteriorate as the medium contrast increases.

The rest of the paper is organized as follows. In the next section, we present the construction of the proposed method for (1.1). The convergence analysis for the multiscale solution is proposed in Section 3. Then, some numerical experiments are investigated in Section 4 to confirm the theoretical results. Finally, the conclusions are given in Section 5.

## 2 Preliminaries

### 2.1 Description of NLMC method

The solution of (1.1) satisfies

 a(u,v)=(f,v)∀v∈H10(Ω), (2.1)

where .

Next, the notations of the fine grids and coarse grids are introduced. Let be a coarse-grid of the domain and be a conforming fine triangulation of . We assume that is a refinement of , where and represent the fine and coarse mesh sizes, respectively. Let be the -th coarse block and let be the corresponding oversampled region obtained by enlarging the coarse block by coarse grid layers (See Figure 1 for an illustration). We let be the number of elements in . Furthermore, each coarse block is decomposed into different regions and is the number of regions within coarse block . In addition, we require that within each region , should satisfy and the contrast ratio should be relatively small. In addition, we define for any . We remark that each region is a continuum.

Consider an oversampling region of the coarse block , then the multiscale basis function is constructed by minimizing subject to the following conditions

 1|Knl|∫Knlψ(j)i=δliδnj∀Knl⊂Ki,m,

where is the Dirac delta function and denotes the area of . We can see that has mean value on the -th region within the coarse block and mean in other regions inside the oversampling domain.

We remark that the above minimization problem is implicit, to solve it explicitly, we can write down the following equivalent variational formulation over each :

 a(ψji,v)+∑Knl⊂Ki,mλnl∫Knlvdx =0∀v∈H10(Ki,m), (2.2) ∫Knlψjidx =∫Knlδliδnjdx∀Knl⊂Ki,m, (2.3)

where and is a piecewise constant function with respect to each region of , and denotes restricted to . An illustration of the multiscale basis functions can be found in Figure 2.

Then we obtain our multiscale space

 Vms=span{ψ(j)i,ms}.

The resulting coarse grid equation can be written as

 a(¯u,v)=(f,v)∀v∈Vms.

The construction of the local multiscale basis function is motivated by the global basis construction as defined below, and in the subsequent analysis we will exploit the global basis functions to show the convergence analysis. The global basis function is defined by

 ψ(j)i=argmin{a(q(j)i,q(j)i)|q(j)i∈H10(Ω),1|Knl|∫Knlq(j)idx=δliδnj,∀Knl⊂Ω}. (2.4)

Out multiscale finite element space is defined by

 Vglo=span{ψ(j)i|1≤i≤N,1≤j≤li}.

For later analysis, we define to be the projection which is defined for each region as

 πij(v)=1|Kji|∫Kjivdx∀v∈L2(Ω)

and

 π(v)=N∑i=1li∑j=1πij(v).

In addition, we define as the null space of the projection , namely, . Then for any , we have

 a(ψ(j)i,v)=0∀v∈~V.

We remark that and interested readers can refer to [8] for the explanations.

The approximate solution obtained in the global multiscale space is defined by

 a(uglo,v)=(f,v)∀v∈Vglo. (2.5)

For later analysis, we define . In addition, for a given subdomain , we define the local -norm by .

### 2.2 Computational issue

For the convenience of the readers, we write down the implementation of the proposed method as follows.

1. Calculate the multiscale basis functions by solving (2.2)-(2.3) for each region .

2. Generate the projection matrix

 RT=[ψ(1)1⋯,ψ(l1)1,ms,⋯,ψ(1)N,ms,⋯,ψ(lN)N,ms],

where is a column vector using its representation in the fine grid.

3. Construct the coarse grid system

 RART¯u=Rb

and solve the above equation to get .

Note that the downscale solution can be defined by . Our coarse grid solutions have physical meaning, which is the average value of the solution on each region .

## 3 Error analysis

In this section, we will carry out the error analysis for the proposed method. We first show the convergence of the global basis function defined in (2.4), then we show the decay property of the local multiscale basis function, using which the convergence of the multiscale solution can be obtained.

### 3.1 Convergence

This subsection presents the convergence of the approximate solution obtained in (2.1) as stated in the next lemma.

###### Lemma 3.1.

Let be the solution in (2.1) and be the solution in (2.5), then we have

 ∥u−uglo∥a ≤CHC1/2ratio∥κ−1/2f∥0.
###### Proof.

By the definitions of and , we have

 a(u,v) =(f,v)∀v∈H10(Ω), a(uglo,v) =(f,v)∀v∈Vglo.

Combining these two equations, we can get

 a(u−uglo,v)=0∀v∈Vglo.

So, we have . It then follows that

 a(u−uglo,u−uglo)=a(u,u−uglo)=(f,u−uglo)≤∥κ−1/2f∥0∥κ1/2(u−uglo)∥0,

Since , the Poincaré inequality yields

 ∫Kji(u−uglo)2≤CH2∫Kji|∇(u−uglo)|2.

Therefore, the preceding arguments reveal that

 ∥u−uglo∥2a≤CHC1/2ratio∥κ−1/2f∥0∥u−uglo∥a,

which gives the desired estimate.

### 3.2 Decay property of the multiscale basis functions

This section aims to proving the global basis functions are localizable. To this end, for each coarse block , we define to be a bubble function and , where is barycentric coordinates and denotes the fine grids restricted to , and more information regarding the bubble function can be found in [29].

The next lemma considers the following minimization problem defined on a coarse block :

 v(j)i=argmin{a(q(j)i,q(j)i)|q(j)i∈H10(Ki),πil(q(j)i)=vaux∀l,=1⋯,li} (3.1)

for a given .

###### Lemma 3.2.

For all , there exists a function such that

 π(v)=vaux,∥v∥2a≤D∥κ1/2vaux∥20,supp(v)⊂supp(vaux).
###### Proof.

Let . The minimization problem is equivalent to the following variational problem: find and such that

 ai(v(j)i,w)+∑Kli⊂Kiμl∫Kliwdx =0∀w∈H10(Ki), (3.2) ∫Kliv(j)idx =∫Klivauxdx∀l=1,⋯,li. (3.3)

Let . Note that, by the mixed finite element theory (cf. [5]), the well-posedness of the minimization problem is equivalent to the existence of a function such that

 si(v,vaux)≥C∥vaux∥20,Ki,∥v∥a(Ki)≤C∥vaux∥0,Ki.

Note that is supported in . We let . By the definition of , we have

 si(v,vaux)=∑Kli⊂Ki∫KliBv2aux≥C∥vaux∥20,Ki.

 ∥v∥2a(Ki)=∥Bvaux∥2a(Ki)≤C∥v∥a(Ki)∥κ1/2vaux∥0,Ki,

Thus

 ∥v∥a(Ki)≤C∥κ1/2vaux∥0,Ki

and the minimization problem (3.1) has a unique solution . Therefore, and satisfy (3.2)-(3.3). From (3.3), we can obtain . The assertion follows.

The rest of this section attempts to estimating the difference between the global and multiscale basis functions. For this purpose, we first introduce some notations used for the subsequent analysis. We define the cutoff function with respect to these oversampling domains. For each , we recall that is the oversampling coarse region by enlarging by coarse grid layers. For , we define such that and

 χM,mi =1inKi,m, (3.4) χM,mi =0inΩ∖Ki,M. (3.5)

Note that we have and are the standard multiscale finite element (MsFEM) basis functions (cf. [18]).

The next lemma shows the difference between the global and multiscale basis functions, which will play an important role in the proof of the convergence of the multiscale solution.

###### Lemma 3.3.

We consider the oversampled domain with . That is, is an oversampled region by enlarging by grid layers. Let be the Dirac delta function. We let be the multiscale basis functions obtained in (2.2)-(2.3) and let be the global multiscale basis functions obtained in (2.4). Then we have

 ∥ψ(j)i−ψ(j)i,ms∥2a≤CE∥κ1/2δlj∥0,Ki∀l,j=1⋯,li

and

 E=D2(1+CratioH2)(1+12D1/2HC1/2ratio)1−k. (3.6)
###### Proof.

For the given , by Lemma 3.2, there exists a such that

 πil(~ϕ(j)i)=δlj,∥~ϕ(j)i∥2a≤D∥κ1/2δlj∥20andsupp(~ϕ(j)i)⊂Ki. (3.7)

We let , then we have . Therefore, . We see that and satisfy

 a(ψ(j)i,v)+∑Kli⊂Ωμ(l)i∫Klivdx=0∀v∈H10(Ω) (3.8)

and

 a(ψ(j)i,ms,v)+∑Kli⊂Ki,kμ(l)i,ms∫Klivdx=0∀v∈H10(Ki,k) (3.9)

for some , . Subtracting the above two equations and restricting , we have

 a(ψ(j)i−ψ(j)i,ms,v)=0∀v∈~V0(Ki,k).

Here, we have . Therefore, for , we can get

 ∥ψ(j)i−ψ(j)i,ms∥2a =a(ψ(j)i−ψ(j)i,ms,ψ(j)i−ψ(j)i,ms) =a(ψ(j)i−ψ(j)i,ms,ψ(j)i−~ϕ(j)i−ψ(j)i,ms+~ϕ(j)i)=a(ψ(j)i−ψ(j)i,ms,η−v),

where . Thus, we obtain

 ∥ψ(j)i−ψ(j)i,ms∥a≤∥η−v∥a. (3.10)

Now, we will estimate . We consider the th coarse block . For this block, we consider two oversampled regions and . Using these two overampled regions, we define the cutoff function with the properties in (3.4)-(3.5), where we take and . For any coarse block by (3.4), we have on . Since , we have

 ∑Knj⊂Kj∫Knjχk,k−1iη=∑Knj⊂Kj∫Knjη=0.

From the above result and the fact that in , we have

 supp(π(χk,k−1iη))⊂Ki,k∖Ki,k−1.

By Lemma 3.2, for the function , there is such that and . Moreover, it also follows from Lemma 3.2, the definition of and the Cauchy-Schwarz inequality that

 ∥μ∥a(Ki,k∖Ki,k−1)≤D1/2∥κ1/2π(χk,k−1iη)∥0,Ki,k∖Ki,k−1≤D1/2∥κ1/2χk,k−1iη∥0,Ki,k∖Ki,k−1, (3.11)

Hence, taking in (3.10), we can obtain

 ∥ψ(j)i−ψ(j)i,ms∥a≤∥η−v∥a≤∥(1−χk,k−1i)η∥a+∥μ∥a(Ki,k∖Ki,k−1). (3.12)

Next, we will estimate the two terms on the right hand side of (3.12).

Step 1: We first estimate the first term in (3.12). By a direct computation, we have

 ∥(1−χk,k−1i)η∥2a≤2(∫Ω∖Ki,k−1κ(1−χk,k−1i)2|∇η|2+∫Ω∖Ki,k−1κ|∇χk,k−1i|2η2).

Note that, we have . For the second term on the righ hand side of the above inequality, we will use the fact that and the Poincaré inequality

 ∥(1−χk,k−1i)η∥2a≤2(1+H2Cratio)∫Ω∖Ki,k−1κ|∇η|2.

We will estimate the right hand side in Step 3.

Step 2: We will estimate the second term on the right hand side of (3.12). By (3.11), the fact that and the Poincaré inequality, we have

 ∥μ∥2a(Ki,k∖Ki,k−1)≤D∥κ1/2χk,k−1iη∥20,Ki,k∖Ki,k−1≤DH2Cratio∫Ki,k∖Ki,k−1κ|∇η|2.

Combining Steps 1 and 2, we obtain

 ∥ψ(j)i−ψ(j)i,ms∥2a≤2D(1+CratioH2)∥η∥2a(Ω∖Ki,k−1). (3.13)

Step 3: Finally, we will estimate the term . We will first show that the following recursive inequality holds

 ∥η∥a(Ω∖Ki,k−1)≤(1+12HD1/2C1/2ratio)−1∥η∥2a(Ω∖Ki,k−2), (3.14)

where . Using (3.14) in (3.13), we can get

 ∥ψ(j)i−ψ(j)i,ms∥2a≤2D(1+CratioH2)(1+12HD1/2C1/2ratio)−1∥η∥2a(Ω∖Ki,k−2). (3.15)

By using (3.14) again in (3.15), we can obtain

 ∥ψ(j)i−ψ(j)i,ms∥2a ≤2D(1+CratioH2)(1+12HD1/2C1/2ratio)1−k∥η∥2a(Ω∖Ki) ≤2D(1+CratioH2)(1+12HD1/2C1/2ratio)1−k∥η∥2a.

By employing the definition of , the energy minimizing property of and Lemma 3.2, we have

 ∥η∥2a=∥ψ(j)i−~ϕ(j)i∥a≤2∥~ϕ(j)i∥a≤2D1/2∥κ1/2δlj∥0,Ki∀l,j=1,⋯,li.

Step 4: We will prove the estimate (3.14). Let . Then we see that in and otherwise. Then we have

 ∥η∥2a(Ω∖Ki,k−1)≤∫Ωκξ2|∇η|2=∫Ωκ∇η⋅∇(ξ2η)−2∫Ωκξη∇ξ∇η. (3.16)

We estimate the first term in (3.16). For the function , using Lemma 3.2, there exists such that and . For any coarse elements , since on , we have for any

 sm(ξ2η,ϕ(n)m)=0∀n=1,…,lm.

On the other hand, since in , we have

 sm(ξ2η,ϕ(n)m)=0∀n=1,…,lm,∀Km⊂Ki,k−2.

From the above two conditions, we see that and consequently . Note that, since , we have . We also note that . By (3.7), the functions and have disjoint supports, so . Then, by the definition of , we have

 a(η,ξ2η−γ)=a(ψ(i)j,ξ2η−γ).

By the construction of , we have . Then we can estimate the first term in (3.16) by the Cauchy-Schwarz inequality and Lemma 3.2

 ∫Ωκ∇η⋅∇(ξ2η) =∫Ωκ∇η⋅∇γ ≤D1/2∥η∥a(Ki,k−1∖Ki,k−2)∥κ1/2π(ξ2η)∥0,Ki,k−1∖Ki,k−2.

For all coarse elements and assume that within , since , we have from the Poincaré inequality that

 ∥κ1/2π(ξ2η)∥20,K≤κ1∥ξ2η∥20,K≤CratioH2∫Kκ|∇η|2.

Summing the above over all coarse elements , we have

 ∥κ1/2π(ξ2η)∥0,Ki,k−1∖Ki,k−2≤C1/2ratioH∥η∥a(Ki,k−1∖Ki,k−2).

To estimate the second term in (3.16), we have from the Poincaré inequality

 2∫Ωκξη∇ξ⋅∇η≤2∥κ1/2η∥0,Ki,k−1∖Ki,k−2∥η∥a(Ki,k−1∖Ki,k−2)≤2HC1/2ratio∥η∥2a(Ki,k−1∖Ki,k−2).

Hence, the preceding arguments yield the upper bound for (3.16)

 ∥η∥2a(Ω∖Ki,k−1)≤2C1/2ratioD1/2H∥η∥2a(Ki,k−1∖Ki,k−2).

Thus

 ∥η∥2a(Ω∖Ki,k−2)=∥η∥2a(Ω∖Ki,k−1)+∥η∥2a(Ki,k−1∖Ki,k−2)≥(1+12D1/2HC1/2ratio)∥η∥2a(Ω∖Ki,k−1).

###### Lemma 3.4.

With the same assumptions as in Lemma 3.3, we can obtain

 ∥N∑i=1(ψ(j)i−ψ(j)i,ms)∥2a ≤C(k+1)2N∑i=1∥ψ(j)i−ψ(j)i,ms∥2a.
###### Proof.

Let . By the constructions in (2.2)-(2.3) and (2.4) and Lemma 3.2, there is such that

 π(zi)=π((1−χk+1,ki)w),supp(zi)⊂Ki,k+1∖Ki,k,∥zi∥a≤D∥κ1/2π((1−χk+1,kiw))∥0.

It then follows from (3.8) and (3.9) that

 a(ψ(j)i−ψ(j)i,ms,v)+∑Kli⊂Ki,k(μ(l)i−μ(l)i,ms)∫Klivdx=0∀v∈H10(Ki,k). (3.17)

Putting in (3.17), we can obtain

 a(ψ(j)i−ψ(j)i,ms,((1−χk+1,k