Constraint Energy Minimizing Generalized Multiscale Discontinuous Galerkin Method
Siu Wun Cheung Eric T. Chung Wing Tat Leung




Constraint Energy Minimizing Generalized Multiscale Discontinuous Galerkin Method
    
Siu Wun CheungDepartment of Mathematics, Texas A&M University, College Station, TX 77843, USA ()
Eric T. ChungDepartment of Mathematics, The Chinese University of Hong Kong, Shatin, New Territories, Hong Kong SAR, China ()
Wing Tat LeungInstitute for Computational Engineering and Sciences, 
The University of Texas at Austin, Austin, Texas, USA ()

     
==================================================================================================================================================================================================================================================================================================================================================================




Numerical simulation of flow problems and wave propagation in heterogeneous media 
has important applications in many engineering areas. 
However, numerical solutions on the fine grid are often prohibitively expensive, 
and multiscale model reduction techniques are introduced to efficiently solve for 
an accurate approximation on the coarse grid.
In this paper, we propose an energy minimization based multiscale model reduction approach 
in the discontinuous Galerkin discretization setting.  
The main idea of the method is to extract the non-decaying component in the high conductivity regions by 
identifying dominant modes with small eigenvalues of local spectral problems, 
and define multiscale basis functions in coarse oversampled regions by constraint energy minimization problems. 
The multiscale basis functions are in general discontinuous on the coarse grid 
and coupled by interior penalty discontinuous Galerkin formulation. 
The minimal degree of freedom in representing high-contrast features 
is achieved through the design of local spectral problems, 
which provides the most compressed local multiscale space. 
We analyze the method for solving Darcy flow problem 
and show that the convergence is linear in coarse mesh size 
and independent of the contrast, provided that the oversampling size is appropriately chosen. 
Numerical results are presented to show the performance of the method for 
simulation on flow problem and wave propagation in high-contrast heterogeneous media.




§ INTRODUCTION




Many engineering applications require numerical simulation in heterogeneous media with multiple scales and high contrast.  
For example, Darcy flow equation in heterogeneous media is used to describe 
fluid flow in porous medium in reservoir simulation, 
and wave equation in heterogeneous media has been widely used for subsurface modeling. 
Numerical solutions on the fine grid are often prohibitively expensive in these complex multiscale problems. 


To this end, extensive research effort had been devoted to developing efficient methods 
for solving multiscale problems at reduced expense, 
for example, numerical homogenization approaches <cit.>
and multiscale methods, including
Multiscale Finite Element Methods (MsFEM) <cit.>, 
Variational Multiscale Methods (VMS) <cit.>, 
Heterogeneous Multiscale Methods (HMM) <cit.> and 
and Generalized Multsicale Finite Element Methods (GMsFEM) 
<cit.>. 
The common goal of these methods is to construct numerical solvers 
on the coarse grid, which is typically much coarser than the fine grid
which captures all the heterogeneities in the medium properties. 
In numerical homogenization approaches, effective properties are computed and 
the global problem is formulated and solved on the coarse grid. 
However, these approaches are limited to the cases when the medium properties possess scale separation.
On the other hand, multiscale methods construct of multiscale basis functions which are 
responsible for capturing the local oscillatory effects of the solution. 
Once the multiscale basis functions, coarse-scale equations are formulated.
Moreover, fine-scale information can be recovered by the coarse-scale coefficients and mutliscale basis functions. 


Many existing mutliscale methods, such as MsFEM, VMS and HMM, 
construct one basis function per local coarse region to handle the effects of local heterogeneities. 
However, for more complex multiscale problems, each local coarse region contains several high-conductivity regions and 
multiple multiscale basis functions are required to represent the local solution space. 
GMsFEM is developed to allow systematic enrichment of the coarse-scale space with fine-scale information 
and identify the underlying low-dimensional local structures for solution representation. 
The main idea of GMsFEM is to extract local dominant modes by carefully designed local spectral problems in coarse regions, 
and the convergence of the GMsFEM is related to eigenvalue decay of local spectral problems. 
For a more detailed discussion on GMsFEM, we refer the readers to 
<cit.> 
and the references therein. 
Our method developed in this work is motivated by GMsFEM and achieves spectral convergence. 
Through the design of local spectral problems, our method results in the minimal degree of freedom 
in representing high-contrast features. 


For typical mesh-based numerical discretization, such as the finite element method and finite difference method,
the most important issue of the solution accuracy is mesh convergence. 
However, for multiscale problems, it is difficult to adjust coarse-grid size based on scales and contrast, 
and it is desirable to have convergence independent of these physical parameters. 
it becomes non-trivial to derive multiscale methods with convergence on coarse mesh size independent of scales and contrast.  
Recently, several multiscale methods with mesh convergence are developed. 
<cit.>. 
This idea can be combined with the use of local spectral problems for 
achieving both spectral convergence and mesh convergence 
<cit.>. 

On the other hand, to overcome the difficulty of stability and conservation for 
convection-dominated problems and wave propagations in heterogeneous media, 
multiscale methods in the discontinuous Galerkin (DG) framework have been investigated 
<cit.>. 
In these approaches, unlike conforming finite element formulations, 
multiscale basis functions are in general discontinuous on the coarse grid, 
and stabilization or penalty terms are added to ensure well-posedness of the global problem. 
Our goal in this work is to develop a robust multiscale method in interior penalty discontinuous Galerkin formulation, 
which exhibits both spectral convergence and mesh convergence. 


In this paper, we present Constraint Energy Minimizing Generalized Multiscale
Discontinuous Galerkin Method (CEM-GMsDGM). 
There are two key ingredients of the presented approach.  
The first main ingredient is the local spectral problems in each coarse block for identification of auxiliary basis functions. 
The low-energy dominant modes, which are eigenvectors corresponding to small eigenvalues of local spectral problems, 
are used as auxiliary basis functions for further construction. 
The auxiliary basis functions possess the information related to high conductivity channels and 
it suffices to use the same number of auxiliary basis functions as the number of channels in a coarse block.
The second ingredient is the constraint energy minimization problems for definition of multiscale basis functions. 
Each of the auxiliary basis functions sets up an independent constraint and 
uniquely defines a corresponding multiscale basis function. 
The multiscale basis functions will then be used to span the multiscale space and 
used to solve the coarse-scale global problem in IPDG formulation.
We remark that the local spectral problems and the constraint energy minimization problems 
are carefully designed and supported by our analysis. 
Thanks to the design of local spectral problems, 
the auxiliary space is of minimal dimension for representing high-contrast features 
and obtaining a contrast-independent convergence. 
Due to the fact that the dimensions of the auxiliary space and the multiscale space are identical, 
the multiscale space is of minimum dimension as well. 
In the construction of multiscale basis functions, the constraints are responsible for 
handling non-decaying components represented by the auxiliary basis functions in the high conductivity regions 
and achieving linear convergence in coarse mesh size. 
On the other hand, the multiscale basis functions are supported in oversampled coarse regions 
and allowed to have discontinuity on the coarse grid. 
Therefore, the IPDG bilinear form is also used to define the energy term in the constraint energy minimization problems. 
The advantages of the method is verified both theoretically and numerically. 
We analyze the method for solving Darcy flow problem 
and establish a criterion for the oversampling size 
which is sufficient for linear coarse-mesh convergence independent of the contrast.
Numerical results are presented to show the performance of the method for 
simulation on flow problem and wave propagation in high-contrast heterogeneous media. 


The paper is organized as follows. In Section <ref>, we will introduce the notions of grids, 
and essential discretization details such as DG finite element spaces and IPDG formulation on the coarse grid. 
The details of the proposed method will be presented in Section <ref>. 
The method will be analyzed in Section <ref>.
Numerical results will be provided in Section <ref>.
Finally, a conclusion will be given in Section <ref>.



§ PRELIMINARIES


We consider the following high-contrast flow problem

    - div(κ∇ u) = f  in Ω,

subject to the homogeneous Dirichlet boundary condition u = 0 on ∂Ω, 
where Ω⊂ℝ^2 is the computational domain and f is a given source term. 
We assume that the permeability field κ is highly heterogeneous 
with very high contrast κ_0 ≤κ≤κ_1. 

Next, we introduce the notions of coarse and fine meshes. 
We start with a usual partition 𝒯^H of Ω into finite elements, 
which does not necessarily resolve any multiscale features. 
The partition 𝒯^H is called a coarse grid and 
a generic element K in the partition 𝒯^H is called a coarse element. 
Moreover, H > 0 is called the coarse mesh size.
We let N_c be the number of coarse grid nodes and 
N be the number of coarse elements. 
We also denote the collection of all coarse grid edges by ℰ^H.
We perform a refinement of 𝒯^H to obtain a fine grid 𝒯^h, 
where h > 0 is called the fine mesh size. 
It is assumed that the fine grid is sufficiently fine to resolve the solution. 
An illustration of the fine grid and the coarse grid and a coarse element are shown in Figure <ref>. 



We are now going to discuss the discontinuous Galerkin (DG) discretization 
and the interior penalty discontinuous Galerkin (IPDG) global formulation. 
For the i-th coarse block K_i, we denote the restriction of the Sobolev space H_0^1(Ω)
on K_i by V(K_i). We let V_h(K_i) be 
the conforming bilinear elements defined on the fine grid 𝒯^h in K_i, i.e.

    V_h(K_i) = { v ∈ V(K_i): v |_τ∈ℚ^1(τ)  for all τ∈𝒯^h  and τ⊂ K_i},

where ℚ^1(τ) stands for the bilinear element on the fine grid block τ. 
The DG approximation space is then given by the space of 
coarse-scale locally conforming piecewise bilinear 
fine-grid basis functions, namely 

    V_h = ⊕_i=1^N V_h(K_i).

We remark that functions in V_h are continuous within coarse blocks, 
but discontinuous across the coarse grid edges in general. 
The global formulation of IPDG method then reads:
find u_h ∈ V_h such that

    a_DG(u_h ,w) = ∫_Ω fw   dx  for all  w ∈ V_h,

where the bilinear form a_DG is defined by:

    a_DG(v,w) 
        = ∑_K ∈𝒯^H∫_K κ∇ v ·∇ w   dx 
    - ∑_E ∈ℰ^H∫_E {κ∇ v · n_E } w   dσ
          - ∑_E ∈ℰ^H∫_E {κ∇ w · n_E } v   dσ 
    + γh∑_E ∈ℰ^H∫_E κ v  w   dσ,

where γ > 0 is a penalty parameter and 
n_E is a fixed unit normal vector defined on the coarse edge E ∈ℰ^H. 
Note that, in (<ref>), the average and the jump operators 
are defined in the classical way. 
Specifically, consider an interior coarse edge E ∈ℰ^H and 
let K^+ and K^- be the two coarse grid blocks sharing the edge E, 
where the unit normal vector n_E is pointing from K^+ to K^-. 
For a piecewise smooth function G with respect to the coarse grid 𝒯^H, we define

    { G }    = 12(G^+ + G^-), 
     G     = G^+ - G^-,

where G^+ = G|_K^+ and G^- = G|_K^-. 
Moreover, on the edge E, we define κ = (κ_K^+ + κ_K^-)/2, 
where κ_K^± is the maximum value of κ over K^±. 
For a coarse edge E lying on the boundary ∂Ω, we define
{G}= G  =G, and κ=κ_K on E, 
where we always assume that n_E is pointing outside of Ω. 

First, we define the energy norm on the space V of coarse-grid piecewise smooth functions by

    w _a^2 = a_DG(w,w)  for all  w ∈ V.

We also define the DG-norm on V by 

    w _DG^2 = ∑_K ∈𝒯^H∫_K κ|∇ w |^2   dx 
    + γh∑_E ∈ℰ^H∫_E κ w ^2   dσ for all  w ∈ V.

The two norms are equivalent on the subspace of piecewise bi-cubic polynomials in V: 
there exists C_0 ≥ 1 such that 

    C_0^-1 w _a≤ w _DG≤ C_0  w _a.

The continuity and coercivity results of the bilinear form a_DG 
with respect to the DG-norm is ensured by a sufficiently large penalty parameter γ.
While the method works well for general highly heterogeneous field κ, 
we assume κ is piecewise constant on the fine grid 𝒯^h 
for the sake of simplicity in our analysis presented in Section <ref>. 



§ METHOD DESCRIPTION


In this section, we will present the construction of the multiscale basis functions. 
First, we will use the concept of GMsFEM to construct our auxiliary multiscale basis functions 
on a generic coarse block K in the coarse grid. 
We consider V_h(K_i) as the snapshot space in K_i and 
perform a dimension reduction through a spectral problem, 
which is to find a real number λ_j^(i) 
and a function ϕ_j^(i)∈ V_h(K_i) such that

    a_i(ϕ_j^(i), w) = λ_j^(i) s_i(ϕ_j^(i), w)  for all  w ∈ V_h(K_i),

where a_i is a symmetric non-negative definite bilinear operator 
and s_i is a symmetric positive definite bilinear operators defined on V_h(K_i) × V_h(K_i). 
We remark that the above problem is solved on the fine mesh in the actual computations. 
Based on our analysis, we can choose

    a_i(v,w)     = ∫_K_iκ∇ v ·∇ w   dx, 
    
    s_i(v,w)     = ∫_K_iκ v w   dx,

where κ̃ = ∑_j=1^N_cκ|∇χ_j^ms|^2 and 
{χ_j^ms}_j=1^N_c are the standard multiscale finite element (MsFEM) basis functions.
We let λ_j^(i) be the eigenvalues of (<ref>) 
arranged in ascending order in j, and use the first L_i eigenfunctions 
to construct our local auxiliary multiscale space

    V_aux^(i) = span{ϕ_j^(i): 1 ≤ j ≤ L_i}.

The global auxiliary multiscale space V^h_aux is then defined as 
the sum of these local auxiliary multiscale spaces

    V_aux = ⊕_i=1^N V_aux^(i).

For the local auxiliary multiscale space V^(i)_aux, the bilinear form s_i 
in (<ref>) defines an inner product 
with norm v_s(K_i) = s(v,v)^1/2. 
These local inner products and norms provide a natural definitions of inner
product and norm for the global auxiliary multiscale space V_aux, which are defined by

    s(v,w)     = ∑_i=1^N s_i(v,w)  for all  v,w ∈ V_aux, 
    v_s     = s(v,v)^1/2 for all  v ∈ V_aux.

We note that s(v, w) and v_s are also an inner product and norm for the space V_h.
Before we move on to discuss the construction of multiscale basis functions, 
we introduce some tools which will be used to describe our method and analyze the convergence.
We first introduce the concept of ϕ-orthogonality. 
For 1 ≤ i ≤ N and 1 ≤ j ≤ L_i, 
in coarse block K_i, given auxiliary basis function ϕ_j^(i)∈ V_aux,
we say that ψ∈ V_h is ϕ_j^(i)-orthogonal if

    s (ψ, ϕ_j'^(i')) = δ_i,i'δ_j,j' for all  1 ≤ j' ≤ L_i' and  1 ≤ i' ≤ N.

We also introduce a projection operator π: V_h → V_aux by 
π = ∑_i=1^N π_i, where

    π_i(v) = ∑_j=1^L_is_i(v,ϕ_j^(i))s_i(ϕ_j^(i),ϕ_j^(i))ϕ_j^(i) for all  v ∈ V_h,  for all  i = 1,2,…,N.

Next, we construct our global multiscale basis functions in V_h.
The global multiscale basis function ψ_j^(i)∈ V_h is defined as the solution of 
the following constrained energy minimization problem

    ψ_j^(i) = { a_DG(ψ, ψ) : 
    ψ∈ V_h  is ϕ_j^(i)-orthogonal}.

By introducing a Lagrange multiplier, 
the minimization problem (<ref>) is equivalent to the following variational problem: 
find ψ_j^(i)∈ V_h and μ_j^(i)∈ V_aux^(i) such that

    a_DG(ψ_j^(i), ψ) + s_i(ψ, μ_j^(i))     = 0  for all ψ∈ V_h, 
    
    s_i(ψ_j^(i) - ϕ_j^(i), μ)     = 0  for all μ∈ V_aux^(i).

Now we discuss the construction our localized multiscale basis functions. 
We first denote by K_i,m an oversampled domain formed by 
enlarging the coarse grid block K_i by m coarse grid layers. 
An illustration of an oversampled domain is shown in Figure <ref>.
We introduce the subspace V_h(K_i,m), which contains restriction of 
fine-scale basis functions in V_h on the oversampled domain K_i,m. 
We also define V_h,0(K_i,m) = V_h(K_i,m) ∩ H_0^1(K_i,m) by 
the subspace of functions in V_h(K_i,m) vanishing on the 
boundary of the oversampled domain K_i,m. 

Motivated by the construction of our global multiscale basis functions,
the method for construction of the localized multiscale basis functions are as follows:
The localized multiscale basis function ψ_j,ms^(i)∈ V_h(K_i,m) is defined as the solution of 
the following constrained energy minimization problem

    ψ_j,ms^(i) = { a_DG(ψ, ψ) : 
    ψ∈ V_h(K_i,m)  is ϕ_j^(i)-orthogonal}.

Using the method of Lagrange multiplier, 
the minimization problem (<ref>) is equivalent to the following variational problem: 
find ψ_j,ms^(i)∈ V_h(K_i,m) and μ_j^(i)∈ V_aux^(i) such that

    a_DG(ψ_j,ms^(i), ψ) + s_i(ψ, μ_j^(i))     = 0  for all ψ∈ V_h(K_i,m), 
    
    s_i(ψ_j,ms^(i) - ϕ_j^(i), μ)     = 0  for all μ∈ V_aux^(i).

We use the localized multiscale basis functions to construct
the multiscale DG finite element space, which is defined as

    V_ms = span{ψ_j,ms^(i) : 1 ≤ j ≤ L_i, 1 ≤ i ≤ N }.

We remark that the multiscale finite element space V_ms is a subspace of V_h. 
After the multiscale DG finite element space is constructed, 
the multiscale solution u_ms is given by: find u_ms∈ V_ms such that

    a_DG(u_ms,w) = ∫_Ω fw   dx  for all  w ∈ V_ms.




§ ANALYSIS


In this section, we will analyze the proposed method.
Besides the energy norm and the DG norm, we also define the s-norm on V by

    w _s^2 = ∑_K ∈𝒯^H∫_K κ| w |^2   dx.

Given a subdomain Ω' ⊆Ω formed by a union of coarse blocks K ∈𝒯^H, 
we also define the local s-norm by 

    w _s(Ω')^2 
    = ∑_K ⊆Ω'∫_K κ| w |^2   dx.


The flow of our analysis goes as follows. 
First, we prove the convergence using the global multiscale basis functions.
With the global multiscale basis functions constructed, 
the global multiscale finite element space is defined by

    V_glo = span{ψ_j^(i) : 1 ≤ j ≤ L_i, 1 ≤ i ≤ N },

and an approximated solution u_glo∈ V_glo is given by

    a_DG(u_glo,w) = ∫_Ω fw   dx  for all  w ∈ V_glo.

We remark that the construction of global multiscale basis functions 
motivates the construction of localized multiscale basis functions.
The approximated solution u_glo will also be used in our convergence analysis.
Next, we give an estimate of the difference between 
the global multiscale functions ψ_j^(i) and 
the localized multiscale basis functions ψ_j,ms^(i), 
in order to show that using the multiscale solution u_ms 
provide similar convergence results as the global solution u_glo.
For this purpose, we denote the kernel of 
the projection operator π by V_h. 
Then, for any ψ_j^(i)∈ V_glo, we have 

    a_DG(ψ_j^(i), w) = 0  for all  w ∈V_h,

which implies V_h ⊆ V_glo^⊥, 
where V_glo^⊥ is the orthogonal complement of 
V_glo with respect to the inner product a_DG(·, ·). 
Moreover, since dim(V_glo) = dim(V_aux), 
we have V_h = V_glo^⊥ and 
V_h = V_glo⊕V_h.



 §.§ Convergence result


The convergence analysis will start with the following lemma, 
which concerns about the convergence of the approximated solution 
by the global multiscale basis functions.


Let u_h ∈ V_h be the solution of (<ref>) 
and u_glo∈ V_glo be the solution of (<ref>)
with the global multiscale basis functions defined by 
the constrained energy minimization problem (<ref>).
Then we have u_h - u_glo∈V_h and

    u_h - u_glo_a ≤Λ^-1/2κ^-1/2 f _L^2(Ω),

where

    Λ = min_1 ≤ i ≤ Nλ_L_i+1^(i).

Moreover, if we replace the multiscale partition of unity {χ_j^ms} by 
the bilinear partition of unity, we have

    u_h - u_glo_a ≤ C H Λ^-1/2κ^-1/2 f _L^2(Ω).


By the definitions of u_h in (<ref>) and u_glo in (<ref>), we have

    a_DG(u_h,w)     = ∫_Ω fw   dx  for all  w ∈ V_h, 
    
    a_DG(u_glo,w)     = ∫_Ω fw   dx  for all  w ∈ V_glo.

Since V_glo⊆ V_h, this yields the Galerkin orthogonality property.

    a_DG(u_h - u_glo,w) = 0  for all  w ∈ V_glo,

which implies u_h - u_glo∈ V_glo^⊥ = V_h. 
In particular, if we take w = u_glo in (<ref>), together with (<ref>), we have

    u_h - u_glo_a^2
        = a_DG(u_h, u_h - u_glo) 
        = (f, u_h - u_glo)_0,Ω
       ≤κ^-1/2 f _L^2(Ω) u_h - u_glo_s.

Since u_h - u_glo∈V_h, we have π(u_h - u_glo) = 0. 
Furthermore, since K_i are disjoint, we have π_i(u_h - u_glo) = 0 for all i = 1,2,…,N.
This implies

    u_h - u_glo_s^2
        = ∑_i=1^N  u_h - u_glo_s(K_i)^2 
        = ∑_i=1^N (I - π_i) (u_h - u_glo) _s(K_i)^2

By the s_i-orthogonality of the eigenfunctions ϕ_j^(i), we have

    (I - π_i) (u_h - u_glo) _s(K_i)^2 
       ≤(λ^(i)_L_i+1)^-1 a_i(u_h - u_glo,u_h - u_glo) 
       ≤Λ^-1 a_i(u_h - u_glo,u_h - u_glo).

Therefore, we have

    u_h - u_glo_s^2
       ≤Λ^-1∑_i=1^N a_i(u_h - u_glo,u_h - u_glo) 
       ≤Λ^-1 u_h - u_glo_a^2.

Using (<ref>) and (<ref>), we obtain our desired result. 
The second part of the result follows from the property 
|∇χ_j | = O(H^-1) 
of the bilinear partition of unity.



The next step is to prove the global basis functions are indeed localizable. 
This makes use of the following lemma, which states some approximation properties 
of the projection operator π. 
In the analysis, we will make use of the Lagrange interpolation operator 
and a bubble function in the coarse grid. 
We define the Lagrange interpolation operator
I_h: C^0(Ω) ∩ H_0^1(Ω) → C^0(Ω) ∩ V_h by:
for all u ∈ C^0(Ω) ∩ H_0^1(Ω), 
the interpolant I_h u ∈ C^0(Ω) ∩ V_h is defined 
piecewise on each fine block τ∈𝒯^h by

    (I_h u)(x) = u(x)  for all vectices  x  of τ,

which satisfies the standard approximation properties: 
there exists C_I ≥ 1 such that for every u ∈ C^0(Ω) ∩ H_0^1(Ω), 

    κ^1/2 (u - I_h u) _L^2(τ) + h κ^1/2∇(u - I_h u) _L^2(τ)≤ C_I h κ^1/2∇ u _L^2(τ),

on each fine block τ∈𝒯^h. 
For any coarse grid block K,  we define a bubble function B on K, i.e.
B(x) = 0 for all x ∈∂ K and B(x) > 0 for all x ∈int(K). 
More precisely, we take B =  ∏_j χ_j^ms, where the product is taken 
over all the coarse grid nodes lying on the boundary ∂ K. 
We can then define the constant 

    C_π = sup_K ∈𝒯^H, μ∈ V_aux∫_K κμ^2∫_K κ B μ^2.

In our following analysis, we will assume 
the following smallness criterion on the fine mesh size h:

    C_π C_I (C_𝒯^2 + λ_max) Θ_L^∞(Ω)^1/2 h < 1,

where C_𝒯 is the maximum number of vertices over all coarse elements K ∈𝒯^H and

    λ_max    = max_1 ≤ i ≤ Nλ_L_i^(i), 
    Θ    = ∑_j |∇χ_j^ms|^2.




Assume the smallness criterion (<ref>) on the fine mesh size h. 
For any v_aux∈ V_aux, there exists a function v ∈ C^0(Ω) ∩ V_h such that 

    π(v) = v_aux,    v _a^2 ≤ D  v_aux_s^2,   supp(v) ⊆supp(v_aux),

where the constant D is defined by 

    D = (2C_π(1+C_I^2) ( C_𝒯^2 + λ_max)1 - C_π C_I ( C_𝒯^2+ λ_max) Θ_L^∞(K_i)^1/2 h)^2.


Let v_aux∈ V_aux^(i). 
We consider the following constraint minimization problem on the block K_i: 

    v = { a_DG(v,v):   v ∈ V_h,0(K_i), 
       s_i(v,ν) = s_i(v_aux, ν)  for all ν∈ V_aux^(i)}.

The minimization problem (<ref>) is equivalent to the following variational problem: 
find (v,μ) ∈ V_h,0(K_i) × V_aux^(i) such that 

    a_i(v, w) + s_i(w, μ)     = 0  for all  w ∈ V_h,0(K_i), 
    
    s_i(v - v_aux, ν)     = 0  for all ν∈ V_aux^(i).

The existence of solution of (<ref>) is based on an inf-sup condition: 

    inf_ν∈ V_aux^(i)sup_w ∈ V_h,0(K_i) ∖{0}s_i(w,ν)a_i(w,w)≥β,

where β > 0 is a constant independent to be determined. 
Pick any ν∈ V_aux^(i). 
We take w = I_h(B ν) ∈ C^0(Ω) ∩ V_h. Since ν∈ V_h(K_i) and 
B(x) = 0 for all vertices of K_i, we have w ∈ V_h,0(K_i). 
First, we see that 

    s_i(w,ν) 
        = ∫_K_iκ B ν^2 + ∫_K_iκ(I_h (B ν) - Bν)ν
       ≥ C_π^-1ν_s(K_i)^2 -  I_h (B ν) - Bν_s ν_s(K_i)
       ≥ C_π^-1ν_s(K_i)^2 - C_I h κ^1/2∇ (Bν) _L^2(K_i)ν_s(K_i)
       ≥ C_π^-1ν_s(K_i)^2 - C_I Θ_L^∞(K_i)^1/2 h
    κ^1/2∇ (Bν) _L^2(K_i)ν_s(K_i).

On the other hand, we observe that 

    a_i(w,w)    ≤ 2(κ^1/2∇ (Bν) _L^2(K_i)^2
    + κ^1/2∇ (Bν - I_h(Bν)) _L^2(K_i)^2)  
       ≤ 2(1+C_I^2) κ^1/2∇ (Bν) _L^2(K_i).

It remains to estimate the term κ^1/2∇ (Bν) _L^2(K_i). 
Since 0 ≤χ_j^ms≤ 1, we have 
0 ≤ B ≤ 1 and |∇ B |^2 ≤ C_𝒯^2 Θ. 
Using these facts together with ∇ (Bν)  = (∇ B) ν + B (∇ν), we imply 

    κ^1/2∇ (Bν) _L^2(K_i)^2
       ≤ C_𝒯^2 ν_s(K_i)^2 + a_i(ν,ν) 
       ≤( C_𝒯^2 + λ_max) ν_s(K_i)^2

By taking the inf-sup constant

    β = 1 - C_π C_I ( C_𝒯^2 + λ_max) Θ_L^∞(K_i)^1/2 h2C_π(1+C_I^2) ( C_𝒯^2 + λ_max),

we prove the inf-sup condition (<ref>) and therefore the existence 
of (v,μ) ∈ V_h,0(K_i) × V_aux^(i) in (<ref>). 
It is then direct to check that the solution v ∈ V_h,0(K_i) satisfies the desired properties. 



We remark that, without loss of generality, we can assume D ≥ C_0^2(1+C_I^2).
We are now going to establish an estimate of the difference between 
the global multiscale basis functions and localized multiscale basis functions. 
We will see that the global multiscale basis functions have a decay property, 
and their values are small outside a suitably large oversampled domain. 
We will make use of a cutoff function in our proof. 
For each coarse block K_i and M > m, 
the oversampling regions K_i,M and K_i,m 
define an outer neighborhood and an inner neighborhood respectively. 
We define χ_i^M,m∈span{χ_j^ms} 
such that 0 ≤χ_i^M,m≤ 1 and 

    χ_i^M,m = 1  in  K_i,m and χ_i^M,m = 0  in Ω∖ K_i,M.

Moreover, we define the following DG norm for w ∈ V 
on K_i,M∖ K_i,m. 

    w _DG(K_i,M∖ K_i,m)^2
    = ∑_K_k ⊂ K_i,M∖ K_i,m a_k(w,w) + 
    γh∑_E ∈ℰ^H(K_i,M∖ K_i,m)∫_E κ w ^2 dσ,

where ℰ^H(K_i,M∖ K_i,m) denotes the 
collection of all coarse grid edges in ℰ^H which lie within 
in the interior of K_i,M∖ K_i,m and 
the boundary of K_i,M. 
We remark that the definition also applies to a region Ω∖ K_i,m 
in the case when M is sufficiently large. 



Assume the smallness criterion (<ref>) on the fine mesh size h.
Suppose m > 2 is the number of coarse grid layers 
in the oversampled domain K_i,m extended from 
the coarse grid block K_i. 
Let ϕ_j^(i)∈ V_aux be a given auxiliary multiscale basis function. 
Let ψ_j^(i)∈ V_glo be the global multiscale basis function 
obtained from (<ref>), and 
ψ_j,ms^(i)∈ V_h(K_i,m) be the localized multiscale basis function 
obtained from (<ref>).
Then we have

    ψ_j^(i) - ψ_j,ms^(i)_a^2 ≤ E ϕ_j^(i)_s(K_i)^2,

where E = 40D^3 (1+Λ^-1) (1 + 6D^-2 (1 + Λ^-1/2)^-1)^1-m.

By the variational formulations (<ref>) and (<ref>), we have 

    a_DG(ψ_j^(i)-ψ_j,ms^(i),ψ) + s_i(ψ, μ_j^(i) - μ_j,ms^(i)) = 0  for all ψ∈ V_h(K_i,m).

By Lemma <ref>, there exists ϕ_j^(i)∈ V_h such that 

    π(ϕ_j^(i)) = ϕ_j^(i),   ϕ_j^(i)_a^2 ≤ D ϕ_j^(i)_s(K_i)^2,   supp(ϕ_j^(i)) ⊆ K_i.
 
We take η = ψ_j^(i) - ϕ_j^(i)∈ V_h and 
ζ =  ϕ_j^(i) - ψ_j,ms^(i)∈ V_h(K_i,m). 
By definition, we have π(η) = π(ζ) = 0 and therefore 
η, ζ∈V_h. 
Again, by Lemma <ref>, there exists ρ∈ V_h such that 

    π(ρ) = π(I_h(χ_i^m,m-1η)),   ρ_a^2 ≤ D π(I_h(χ_i^m,m-1η)) _s^2,   supp(ρ) ⊆ K_i,m∖ K_i,m-1.
 
Take τ = ρ - I_h(χ_i^m,m-1η) ∈ V_h. 
Again, π(τ) = 0 and hence τ∈V_h. 
Taking ψ = τ - ζ∈ V_h(K_i,m) in (<ref>) and 
making use of the fact τ - ζ∈V_h, we have

    a_DG(ψ_j^(i)-ψ_j,ms^(i),τ-ζ) = 0,
 
and therefore

    ψ_j^(i) - ψ_j,ms^(i)_a^2 
        = a_DG(ψ_j^(i)-ψ_j,ms^(i),η+ζ) 
        = a_DG(ψ_j^(i)-ψ_j,ms^(i),η+τ) 
       ≤ψ_j^(i) - ψ_j,ms^(i)_a 
    η+τ_a,

which in turn implies 

    ψ_j^(i) - ψ_j,ms^(i)_a^2
       ≤η+τ_a^2 
        =  I_h((1-χ_i^m,m-1)η)+ρ_a^2 
       ≤ 2(  I_h((1-χ_i^m,m-1)η) _a^2 + ρ_a^2) 
       ≤ 2( C_0^2  I_h((1-χ_i^m,m-1)η) _DG^2 + ρ_a^2) 
       ≤ 2( 2C_0^2 (1+C_I^2)  (1-χ_i^m,m-1)η_DG^2 + ρ_a^2).

For the first term on the right hand side of (<ref>), 
by using ∇((1-χ_i^m,m-1)η) = -∇χ_i^m,m-1η+(1-χ_i^m,m-1)∇η 
and 0 ≤ 1-χ_i^m,m-1≤ 1, we have 

    (1-χ_i^m,m-1)η_a^2 ≤ 2( 
    η_DG(Ω∖ K_i,m-1)^2 +
    η_s(Ω∖ K_i,m-1)^2).

For the second term on the right hand side of (<ref>), using the definition of ρ 
in (<ref>) and 0 ≤ 1-χ_i^m,m-1≤ 1, we obtain 

    ρ_a^2 ≤ D π (χ_i^m,m-1η) _s^2 ≤ D χ_i^m,m-1η_s^2 ≤ D η_s(Ω∖ K_i,m-1)^2.

Moreover, since η∈V_h, by the spectral problem (<ref>), we have 

    η_s(Ω∖ K_i,m-1)^2 ≤Λ^-1∑_K_k ⊂Ω∖ K_i,m-1 a_k(η,η).

Combining all these estimates, we obtain 

    ψ_j^(i) - ψ_j,ms^(i)_a^2
    ≤ 10D(1+Λ^-1) η_DG(Ω∖ K_i,m-1)^2.

Next, we will provide a recursive estimate for η in the number of oversampling layers m. 
We take ξ = 1- χ_i^m-1,m-2. Then 
0 ≤ξ≤ 1 and ξ = 1 in Ω∖ K_i,m-1. 
Using ∇(ξ^2 η) = ξ^2 ∇η + 2ξη∇ξ, 
for every K ∈𝒯^H, we have 

    ∫_K κ∇η·∇(ξ^2 η) 
    = ∫_K κ∇η·(ξ^2 ∇η + 2ξη∇ξ)
    = ∫_K κ|∇ (ξη) |^2 - ∫_K κ|∇ξ|^2 η^2.

In addition, using ∇(ξη) = ξ∇η + η∇ξ, 
for every E ∈ℰ^H, we have 

    -∫_E {κ∇η· n_E }ξ^2 η 
    -∫_E {κ∇ (ξ^2 η) · n_E }η + 
    γh∫_E κηξ^2 η
        = -∫_E {κ∇η· n_E }ξ^2 η 
    -∫_E {κ (ξ^2 ∇η + 2ξη∇ξ) · n_E }η + 
    γh∫_E κηξ^2 η
        = -2(∫_E {κξ∇η· n_E }ξη + 
    ∫_E {κη∇ξ· n_E }ξη) + 
    γh∫_E κξη^2  
        = -2∫_E {κ∇ (ξη) · n_E }ξη +
    γh∫_E κξη^2.

Summing over K ∈𝒯^H and E ∈ℰ^H, we obtain 

    ξη_a^2 
    
    ≤ a_DG(η,ξ^2η) + η_s(K_i,m-1∖ K_i,m-2)^2,

where we make use of the fact that ∇ξ = 0 outside K_i,m-1∖ K_i,m-2. 
We start with estimating the first term on the right hand side of (<ref>). 
For any coarse element K_k ∈Ω∖ K_i,m-1, 
since ξ = 1 in K_k and η∈V_h, we have 

    s(ξ^2 η, ϕ_j^(k)) = s(η, ϕ_j^(k)) = 0  for all  j = 1,2,…,L_k.

On the other hand, for any coarse element K_k ∈ K_i,m-2, 
since ξ = 0 in K_k, we have 

    s(ξ^2 η, ϕ_j^(k)) = 0  for all  j = 1,2,…,L_k.

Therefore, supp(π(I_h(ξ^2 η))) ⊂ K_i,m-1∖ K_i,m-2. 
By Lemma <ref>, there exists σ∈ V_h such that 

    π(σ) = π(I_h(ξ^2 η)),   γ_a^2 ≤ D π (I_h(ξ^2 η)) _s^2,   supp(σ) ⊂ K_i,m-1∖ K_i,m-2.

For any coarse element K_k ⊂ K_i,m-1∖ K_i,m-2, 
since 0 ≤ξ≤ 1 and π(η) = 0, we have 

    π(I_h(ξ^2 η)) _s(K_k)^2 
    ≤ I_h(ξ^2 η) _s(K_k)^2 ≤ I_h(η) _s(K_k)^2 
    = η_s(K_k)^2 ≤Λ^-1 a_k(η,η).

Summing over K_k ⊂ K_i,m-1∖ K_i,m-2, we obtain 

    π(I_h(ξ^2 η)) _s^2 ≤Λ^-1∑_K_k ⊂ K_i,m-1∖ K_i,m-2 a_k(η,η).

We take θ = I_h(ξ^2 η) - σ. 
Again, π(θ) = 0 and θ∈V_h, which yields

    a_DG(ψ_j^(i), θ) = 0.

On the other hand, supp(θ) ⊂Ω∖ K_i,m-2 and supp(ϕ_j^(i)) ⊂ K_i. 
Since θ and ϕ_j^(i) has disjoint supports, we have 

    a_DG(ϕ_j^(i), θ) = 0.

Therefore, we obtain 

    a_DG(η,θ) = a_DG(ψ_j^(i)-ϕ_j^(i), θ) = 0.

Recall from the definition that I_h(ξ^2 η) = θ + σ and 
supp(σ) ⊂ K_i,m-1∖ K_i,m-2.
Hence we have

    a_DG(η,I_h(ξ^2η)) 
        = a_DG(η,σ) 
       ≤ C_0η_DG(K_i,m-1∖ K_i,m-2)σ_a 
       ≤ C_0 D^1/2η_DG(K_i,m-1∖ K_i,m-2)π(I_h(ξ^2 η)) _s 
       ≤ D Λ^-1/2η_DG(K_i,m-1∖ K_i,m-2)^2.

On the other hand, making use of the fact that 
ξ^2 = 0 in K_i,m-2 and ξ^2 = 1 in Ω∖ K_i,m-1, 
we observe that ξ^2 η = I_h(ξ^2 η) outside K_i,m-1∖ K_i,m-2. 
Moreover, ξ^2 η - I_h(ξ^2 η) is globally continuous. 
Thus, we obtain 

    a_DG(η,ξ^2η - I_h(ξ^2η)) 
       ≤ C_0^2 η_DG(K_i,m-1∖ K_i,m-2)ξ^2η - I_h(ξ^2η) _DG(K_i,m-1∖ K_i,m-2)
       ≤ C_0^2 C_I η_DG(K_i,m-1∖ K_i,m-2)ξ^2η_DG(K_i,m-1∖ K_i,m-2)
       ≤D2(η_DG(K_i,m-1∖ K_i,m-2)^2 + 
    ξ^2η_DG(K_i,m-1∖ K_i,m-2)^2).

Again, using ∇(ξ^2 η) = ξ^2 ∇η + 2ξη∇ξ, we have 

    ξ^2η_DG(K_i,m-1∖ K_i,m-2)^2 ≤ 
    2 η_DG(K_i,m-1∖ K_i,m-2)^2 + 
    8 η_s(K_i,m-1∖ K_i,m-2)^2.

Combining (<ref>), (<ref>), (<ref>) and (<ref>), we arrive at 

    ξη_a^2  ≤ D ((32 + Λ^-1/2) 
    η_DG(K_i,m-1∖ K_i,m-2)^2 
    + 5η_s(K_i,m-1∖ K_i,m-2)^2).

Moreover, since 
π(η) = 0, we have 

    η_s(K_i,m-1∖ K_i,m-2)≤Λ^-1/2η_DG(K_i,m-1∖ K_i,m-2),

which implies 

    ξη_a^2  ≤ 6D (1 + Λ^-1/2) 
    η_DG(K_i,m-1∖ K_i,m-2)^2.

By the equivalence of norms, we have 

    η_DG(Ω∖ K_i,m-1)^2 ≤ C_0^2 ξη_a^2 
    ≤ 6D^2 (1 + Λ^-1/2) η_DG(K_i,m-1∖ K_i,m-2)^2.

We obtain the recurrence estimate 

    η_DG(Ω∖ K_i,m-2)^2 
        = η_DG(Ω∖ K_i,m-1)^2 + η_DG(K_i,m-1∖ K_i,m-2)^2 
       ≥(1 + 6D^-2 (1 + Λ^-1/2)^-1) η_DG(Ω∖ K_i,m-1)^2.

Inductively, we have 

    η_DG(Ω∖ K_i,m-1)^2 
       ≤(1 + 6D^-2 (1 + Λ^-1/2)^-1)^1-mη_DG(Ω∖ K_i,1)^2  
       ≤ D (1 + 6D^-2 (1 + Λ^-1/2)^-1)^1-mη_a^2.

Combining (<ref>) and (<ref>), we see that 

    ψ_j^(i) - ψ_j,ms^(i)_a^2 ≤ 10D^2 (1+Λ^-1)
    (1 + 6D^-2 (1 + Λ^-1/2)^-1)^1-mη_a^2

By the energy minimizing property of ψ_j^(i), we have 

    η_a ≤ψ_j^(i)_a + ϕ_j^(i)_a ≤ 2 ϕ_j^(i)_a ≤ 2D^1/2ϕ_j^(i)_s(K_i).

We obtain the desired result. 



Now, we are ready to establish our main theorem, 
which estimates the error between the solution u_h 
and the multiscale solution u_ms.



Let u_h ∈ V_h be the solution of (<ref>),
u_glo∈ V_glo be the solution of (<ref>)
with the global multiscale basis functions defined by (<ref>), 
and u_ms∈ V_ms be the multiscale solution of (<ref>)
with the localized multiscale basis functions defined on an oversampled domain 
with m > 2 coarse grid layers by (<ref>).
Then we have 

    u_h - u_ms_a ≤ C Λ^-1/2κ^-1/2 f _L^2(Ω) + 
    Cm^d E^1/2u_glo_s,

Moreover, if we let k = O(log(κ_1H)) and 
replace the multiscale partition of unity {χ_j^ms} by 
the bilinear partition of unity, we have

    u_h - u_ms_a ≤ C H Λ^-1/2κ^-1/2 f _L^2(Ω).


First, we write u_glo in the linear combination of the basis {ψ_k^(j)} 

    u_glo =  ∑_i=1^N ∑_j=1^L_iα_j^(i)ψ_j^(i).

and define u_ms∈ V_ms by 

    u_ms = ∑_i=1^N ∑_j=1^L_iα_j^(i)ψ_j,ms^(i).
 
From (<ref>) and (<ref>), we obtain the Galerkin orthogonality 

    a_DG(u_h - u_ms, w) = 0  for all  w ∈ V_ms,

which gives 

    u_h - u_ms_a≤ u_h - u_ms_a≤ u_h - u_glo_a +  u_glo - u_ms_a.

Using Lemma <ref>, we see that

    u_glo - u_ms_a^2
        = ∑_i=1^N ∑_j=1^L_iα_j^(i) (ψ_j^(i) - ψ_j,ms^(i)) _a^2 
       ≤ Cm^d ∑_i=1^N ∑_j=1^L_iα_j^(i) (ψ_j^(i) - ψ_j,ms^(i)) _a^2 
       ≤ C m^d E ∑_i=1^N ∑_j=1^L_iα_j^(i)ϕ_j^(i)_s^2 
        = C m^d E  u_glo_s^2,

where the last equality follows from the orthogonality of the eigenfunctions in (<ref>). 
Using the estimates (<ref>) and (<ref>) in (<ref>), we have

    u_h - u_ms_a≤Λ^-1/2κ^-1/2 f _L^2(Ω) + 
    C m^d/2 E^1/2 u_glo_s.

This completes the first part of the theorem. 
Next, we assume the partition of unity functions are bilinear, 
and we are going to estimate u_glo_s. 
Using the fact that |∇χ_k | = O(H^-1), we have

    u_glo_s^2 ≤ CH^-2κ_1  u_glo_L^2(Ω)^2.

Then, by Poincaré inequality, we have

    u_glo^2_L^2(Ω)≤ C κ_0^-1 u_glo_a^2.

By taking w = u_glo∈ V_glo in (<ref>), we obtain

    u_glo_a^2 = (f, u_glo)_0,Ω≤κ^-1/2 f _L^2(Ω) u_glo_s.

Combining these estimates, we have

    u_glo_s ≤ CH^-2κ_0^-1κ_1 κ^-1/2 f _L^2(Ω).

To obtain our desired result, we need

    H^-2κ_1 m^d/2E^1/2 = O(1).

Taking logarithm, we have

    log(H^-2) + log(κ_1) + d2log(m) + 1-m2log(1+Λ^-1/2) = O(1).

Thus, taking m = O(log(κ_1H)) 
completes the proof of the second result. 





§ NUMERICAL RESULTS



In this section, we will present numerical examples with high contrast media 
to demonstrate the convergence of our proposed method 
with respect to the coarse mesh size H and 
the number of oversampling layers m, 
and illustrate possible improvements in error robustness with respect to contrast by
employing the idea of constructing multiscale basis function by relaxation method introduced in <cit.>.
Lastly, we examine the performance of applying the method to the wave equation. 
In all the experiments, the IPDG penalty parameter in (<ref>) 
is set to be γ = 4, so as to ensure the coercivity of the bilinear form a_DG.



 §.§ Experiment 1: flow problem


In the first experiment, 
we consider a highly heterogeneous permeability field κ 
in Ω = [0.1]^2 as shown in Figure <ref>, 
with the background value is κ = 1 and 
the value in the channels and inclusions is 10^4. 
and the resolution is 400×400, i.e.
κ is piecewise constant on a fine grid with mesh size h = 1/400. 
The coarse mesh size varies from H = 1/80 to H = 1/10, 
and the number of oversampling layers varies from m = 3 to m = 6. 
In all these combinations, there are no more than 3 high conductivity channels in a 
coarse block K ∈𝒯^H. 
As a result, we have 3 small eigenvalues in a local spectral problem (<ref>), 
and it suffices to use 3 auxiliary basis functions per coarse block to construct the 
correspoding localized multiscale basis functions. 
The source function is taken as 

    f(x,y) = 2π^2 sin(π x) sin(π y)  for all  (x,y) ∈Ω.

Table <ref> records the error when 
we take the number of oversampling layer to be approximately
m ≈ 4 log(1/H)/ log(1/10). 
The results show that the method provides optimal convergence in energy norm, 
which agrees with our theoretical finding in Section <ref>,  
and the L^2 error converges with second order. 
Table <ref> records the error with various number of oversampling layers
and a fixed coarse mesh sizes H = 1/40.
It can be observed that increasing the number of oversampling layers 
improves the quality of approximations, but the decay in error is limited 
when the oversampling region is sufficiently large. 
This numerically verifies that the multiscale basis functions can indeed be localized. 
























































Next, we present the idea of the relaxed formulation of (<ref>).
Instead of using the method of Lagrange multiplier as in (<ref>), 
the ϕ-orthogonality is imposed weakly by a penalty formulation.













The localized multiscale basis function ψ_j,ms^(i)∈ V_h(K_i,m) is defined as the solution of 
the following relaxed constrained energy minimization problem

    ψ_j,ms^(i) = { a_DG(ψ, ψ) + s(π(ψ)-ϕ_j^(i),π(ψ)-ϕ_j^(i)) : 
    ψ∈ V_h(K_i,m)}.

The minimization problem (<ref>) is equivalent to the following variational problem: 
find ψ_j,ms^(i)∈ V_h(K_i,m) such that

    a_DG(ψ_j,ms^(i), ψ) + s(π(ψ_j,ms^(i)), π(ψ)) = s(ϕ_j^(i),π(ψ))  for all ψ∈ V_h(K_i,m).

The construction of multiscale finite element space and coarse-scale model 
then follow (<ref>) and (<ref>) respectively. 
We compare the performance of the multiscale method with 
multiscale basis functions constructed by method of Lagrange multiplier (<ref>) 
and the relaxation method (<ref>) at different contrast values,  
where the coarse mesh size is taken as H = 1/10 
and the number of oversampling layers as m = 4. 
In Table <ref>, we record the energy error and L^2 error 
with different contrast κ_1, 
where κ_1 ≫ 1 is the value of κ in the high conductivity channels.
It can be seen that the relaxation method is more robust with respect to contrast.



































 §.§ Experiment 2: wave propagation


In the second experiment, we consider a simple wave equation in the space-time domain Ω_T = [0,T] ×Ω, 
where [0,T] is the temporal domain and Ω is the spatial domain: 

    ∂^2 u∂ t^2  - div(κ∇ u) = f,

with homogeneous Dirichlet boundary condition u = 0 on [0,T] ×∂Ω. 
Here div and ∇ correspond to the divergence and gradient with respect to the spatial variable x, 
and κ is the bulk modulus which is assumed to be stationary and highly oscillatory. 
Using second-order central difference for temporal discretization 
on a uniform grid points t_k = kΔ t and 
IPDG formulation on the space V_h, 
the fully-discrete DG finite element scheme reads: find u^n+1_h ∈ V_h such that 

    u_h^n+1 - 2u_h^n + u_h^n-1Δ t^2 + a_DG(u_h^n, w) = ∫_Ω f^n w  for all  w ∈ V_h,

where the superscript k stands for the evaluation of a function at the time instant t_k, 
and the initial conditions u^0_h, u^1_h are given. 
Next, we use the CEM method in Section <ref> 
to construct a multiscale finite element space V_ms. 
The coarse-scale model then reads: find u^n+1_ms∈ V_ms such that 

    u_ms^n+1 - 2u_ms^n + u_ms^n-1Δ t^2 + a_DG(u_ms^n, w) = ∫_Ω f^n w  for all  w ∈ V_ms.

In this experiment, we take the bulk modulus on the spatial domain Ω = [0,1]^2 
as part of the Marmousi model as shown in Figure <ref>. 
The fine mesh size is taken as h = 1/256. 
The coarse mesh size varies from H = 1/32 to H = 1/8, 
and the number of oversampling layers varies from m = 4 to m = 6. 
In all these combinations, we use 4 auxiliary basis functions per coarse block to construct the 
correspoding localized multiscale basis functions. 
The source function f is taken as the Ricker wavelet

    f(t,x,y) = t-2/f_04h^2exp( -π^2 f_0^2 (t-2/f_0)^2 ) 
    exp((x-0.5)^2+(y-0.5)^24h^2)  for all  (t,x,y) ∈Ω_T.

where the central frequency is chosen as f_0 = 20.
We iteratively solve for the numerical solution at the final time T = 0.2 
with time step size Δ t = 10^-4. 
Table <ref> records the error of the final solution 
when we take the number of oversampling layer to be approximately
m ≈ 4 log(1/H)/ log(1/8). 
It can been observed that the method results in good accuracy and desired convergence in error. 
Figure <ref> depicts the numerical solutions 
by the fine-scale formulation and the coarse-scale formulation at the final time T = 2. 
The comparison suggests that the CEM method provides very good accuracy 
at a reduced computational expense. 











§ CONCLUSION


In this paper, we present CEM-GMsDGM, 
a local multiscale model reduction approach in the discontinuous Galerkin framework. 
The multiscale basis functions are defined in coarse oversampled regions by a constraint energy minimization problem, 
which are in general discontinuous on the coarse grid, and coupled by the IPDG formulation. 
Thanks to the definition of local spectral problems, 
the dimension of auxiliary space is minimal for sufficiently representing the high conductivity regions, 
and provides the most locally compressed multiscale space. 
In our analysis for the Darcy flow problem, we show that the method provides optimal convergence in the coarse mesh size, 
which is independent of the contrast, provided that the oversampling size is appropriately chosen. 
The convergence of the method for solving Darcy flow is theoretically analyzed and numerically verified. 
Numerical results for applying the method on a simple wave equation are also presented. 


plain




