Multiscale discontinuous Petrov--Galerkin method for the multiscale elliptic problems

TL;DR

This paper introduces a novel multiscale discontinuous Petrov--Galerkin method (MsDPGM) for elliptic problems with multiscale features, combining oversampling basis, Petrov--Galerkin, and discontinuous Galerkin techniques for improved accuracy and efficiency.

Contribution

It develops a new MsDPGM that eliminates resonance errors and reduces computational complexity, with detailed convergence analysis and numerical validation.

Findings

Eliminates resonance error in multiscale elliptic problems

Achieves accurate multiscale solution with reduced computational cost

Validated through numerical experiments with periodic and random coefficients

Abstract

In this paper we present a new multiscale discontinuous Petrov--Galerkin method (MsDPGM) for multiscale elliptic problems. This method utilizes the classical oversampling multiscale basis in the framework of Petrov--Galerkin version of discontinuous Galerkin finite element method, allowing us to better cope with multiscale features in the solution. The introduced MsDPGM takes advantages of the multiscale Petrov--Galerkin method (MsPGM) and discontinuous Galerkin method (DGM), which can eliminate the resonance error completely, and can decrease the computational complexity, allowing for more efficient solution algorithms. Upon the $H^2$ norm error estimate between the multiscale solution and the homogenized solution with the first order corrector, we give a detailed multiscale convergence analysis under the assumption that the oscillating coefficient is periodic. We also investigate the corresponding multiscale discontinuous finite element method (MsDFEM) which coupling the classical oversampling multiscale basis with DGM since it has not been studied detailedly in both aspects of error analysis and numerical tests in the literature. Numerical experiments are carried out for the multiscale elliptic problems with periodic and randomly generated log-normal coefficients to demonstrate the proposed method.

Tables

Table 1. Table 1: Compare different methods to show the accuracy of MsDPGM in periodic case given by ( 6.2 ). ρ = ε = 1 / 100 , d ~ = h = 1 / 32 , γ 0 = 20 formulae-sequence 𝜌 𝜀 1 100 ~ 𝑑 ℎ 1 32 subscript 𝛾 0 20 \rho=\varepsilon=1/100,\tilde{d}=h=1/32,\gamma_{0}=20 .

Relative error	$L^2$	$L^{\infty }$	Energy norm	CPU time(s)
				$T_1$	$T_2$
FEM	0.1150e-00	0.2311e-00	0.8790e-00	–	–
DFEM	0.2667e-00	0.2634e-00	0.5498e-00	–	–
MsPGM	0.7448e-01	0.7342e-01	0.2929e-00	–	–
OMsPGM	0.1430e-01	0.1521e-01	0.1641e-00	–	–
MsDFEM	0.1007e-01	0.1029e-01	0.1629e-00	1.300	0.028
MsDPGM	0.1266e-01	0.1395e-01	0.1631e-00	1.119	0.027

Table 2. Table 2: Convergence with respect to γ 0 subscript 𝛾 0 \gamma_{0} . ρ = ε = 1 / 100 , d ~ = h = 1 / 32 formulae-sequence 𝜌 𝜀 1 100 ~ 𝑑 ℎ 1 32 \rho=\varepsilon=1/100,\tilde{d}=h=1/32 .

Relative Error	$L^2$	$L^{\infty }$	Energy norm
${\gamma }_0=10$	0.1100e-01	0.1266e-01	0.1637e-00
${\gamma }_0=20$	0.1266e-01	0.1395e-01	0.1631e-00
${\gamma }_0=100$	0.1397e-01	0.1496e-01	0.1638e-00
${\gamma }_0=1000$	0.1426e-01	0.1519e-01	0.1641e-00
${\gamma }_0=10000$	0.1429e-01	0.1521e-01	0.1641e-00

Table 3. Table 3: Error with respect to h ℎ h . ρ = ε = 1 / 100 , d ~ = 1 / 32 , γ 0 = 20 formulae-sequence 𝜌 𝜀 1 100 formulae-sequence ~ 𝑑 1 32 subscript 𝛾 0 20 \rho=\varepsilon=1/100,\tilde{d}=1/32,\gamma_{0}=20 .

Relative error	$L^2$	$L^{\infty }$	Energy norm
$h=1/64$	0.1371e-01	0.1474e-01	0.1593e-00
$h=1/32$	0.1266e-01	0.1395e-01	0.1631e-00
$h=1/16$	0.1948e-01	0.2532e-01	0.1870e-00
$h=1/8$	0.5210e-01	0.8248e-01	0.2620e-00

Table 4. Table 4: Error with respect to δ 0 subscript 𝛿 0 \delta_{0} . ρ = ε = 1 / 100 , h = 1 / 32 , γ 0 = 20 formulae-sequence 𝜌 𝜀 1 100 formulae-sequence ℎ 1 32 subscript 𝛾 0 20 \rho=\varepsilon=1/100,h=1/32,\gamma_{0}=20 .

Relative error	$L^2$	$L^{\infty }$	Energy norm
${\delta }_0=1/32$	0.4304e-01	0.4423e-01	0.2184e-00
${\delta }_0=1/16$	0.2924e-01	0.3172e-01	0.1893e-00
${\delta }_0=1/8$	0.1969e-01	0.2156e-01	0.1728e-00
${\delta }_0=1/4$	0.1790e-01	0.2372e-01	0.1653e-00
${\delta }_0=1/2$	0.1531e-01	0.1766e-01	0.1642e-00
${\delta }_0=1$	0.1266e-01	0.1295e-01	0.1631e-00
${\delta }_0=2$	0.1197e-01	0.1398e-01	0.1631e-00

Table 5. Table 5: Relative errors in the L 2 superscript 𝐿 2 L^{2} , L ∞ superscript 𝐿 L^{\infty} and energy norm for the L–shaped problem with periodic coefficient ( 6.2 ). ρ = ε = 1 / 100 , d ~ = h = 1 / 16 , γ 0 = 20 formulae-sequence 𝜌 𝜀 1 100 ~ 𝑑 ℎ 1 16 subscript 𝛾 0 20 \rho=\varepsilon=1/100,\tilde{d}=h=1/16,\gamma_{0}=20 .

Relative error	$L^2$	$L^{\infty }$	Energy norm
MsPGM	0.7765e-02	0.3635e-01	0.2014e-00
OMsPGM	0.6285e-02	0.3277e-01	0.1035e-00
MsDFEM	0.3903e-02	0.2244e-01	0.9260e-01
MsDPGM	0.4654e-02	0.2299e-01	0.9275e-01

Table 6. Table 6: Relative errors in the L 2 superscript 𝐿 2 L^{2} , L ∞ superscript 𝐿 L^{\infty} and energy norm for the L–shaped problem with random coefficient σ 2 = 1.0 superscript 𝜎 2 1.0 \sigma^{2}=1.0 and l 1 = l 2 = 0.01 subscript 𝑙 1 subscript 𝑙 2 0.01 l_{1}=l_{2}=0.01 . d ~ = ρ = h = 1 / 16 , γ 0 = 20 formulae-sequence ~ 𝑑 𝜌 ℎ 1 16 subscript 𝛾 0 20 \tilde{d}=\rho=h=1/16,\gamma_{0}=20 .

Relative error	$L^2$	$L^{\infty }$	Energy norm
MsPGM	0.9074e-00	0.1290e+01	0.6601e+02
OMsPGM	0.9307e-02	0.3851e-01	0.1428e-00
MsDFEM	0.6504e-02	0.3718e-01	0.9931e-01
MsDPGM	0.8587e-02	0.3810e-01	0.1013e-00

Equations

$$\{\begin{array}{rlrl}-\nabla \cdot ({\mathbf{a}}^{\epsilon }(x)\nabla u_{\epsilon }(x))& =f(x)& & \text{in}\Omega ,\\ u_{\epsilon }(x)& =0& & \text{on}\partial \Omega ,\end{array}$$

$$\lambda |\xi {|}^2\le a_{ij}^{\epsilon }(x){\xi }_i{\xi }_j\le \Lambda |\xi {|}^2\forall \xi \in {\mathbf{R}}^n,x\in \bar{\Omega }$$

$$H^s({\mathcal{T}}_h)=\{v\in L^2(\Omega ):\forall K\in {\mathcal{T}}_h,v{|}_K\in H^s(K)\},$$

$$|||v|||_{H^{s}({\cal T}_{h})}=\Big{(}\sum_{K\in{\cal T}_{h}}||v||^{2}_{H^{s}(K)}\Big{)}^{1/2}.$$

$$\left\{v\right\}:=\frac{v{|}_{K_1^e}+v{|}_{K_2^e}}{2},\left[v\right]:=v{|}_{K_1^e}-v{|}_{K_2^e}\forall e=\partial K_1^e\cap \partial K_2^e.$$

$$\left\{v\right\}=\left[v\right]=v{|}_{K_1^e}\forall e=\partial K_1^e\cap \partial \Omega .$$

$${\int }_K{\mathbf{a}}^{\epsilon }\mathrm{NewA}{u}_{\epsilon }\cdot \mathrm{NewA}v\mathrm{d}x-{\int }_{\partial K}{\mathbf{a}}^{\epsilon }\mathrm{NewA}{u}_{\epsilon }\cdot {\mathbf{n}}_{K}v\mathrm{d}s={\int }_{K}fv.$$

$$\sum _{K\in {\mathcal{T}}_h}{\int }_{\partial K}{\mathbf{a}}^{\epsilon }\mathrm{NewA}{u}_{\epsilon }\cdot {\mathbf{n}}_{K}v\mathrm{d}s=\sum _{e\in {\Gamma }_h\cup \partial \Omega }{\int }_e\left[{\mathbf{a}}^{\epsilon }\mathrm{NewA}{u}_{\epsilon }\cdot \mathbf{n}v\right]\mathrm{d}s.$$

$$\sum _{K\in {\mathcal{T}}_h}{\int }_K{\mathbf{a}}^{\epsilon }\mathrm{NewA}{u}_{\epsilon }\cdot \mathrm{NewA}v\mathrm{d}x-\sum _{e\in {\Gamma }_h\cup \partial \Omega }{\int }_e\left\{{\mathbf{a}}^{\epsilon }\mathrm{NewA}{u}_{\epsilon }\cdot \mathbf{n}\right\}\left[v\right]\mathrm{d}s={\int }_{\Omega }fv,$$

$$\begin{array}{rl}a(u,v):& =\sum _{K\in {\mathcal{T}}_h}{\int }_K{\mathbf{a}}^{\epsilon }\mathrm{NewA}u\cdot \mathrm{NewA}v\mathrm{d}x-\sum _{e\in {\Gamma }_h\cup \partial \Omega }{\int }_e\left\{{\mathbf{a}}^{\epsilon }\mathrm{NewA}u\cdot \mathbf{n}\right\}\left[v\right]\mathrm{d}s\\ & +\beta \sum _{e\in {\Gamma }_h\cup \partial \Omega }{\int }_e\left[u\right]\left\{{\mathbf{a}}^{\epsilon }\mathrm{NewA}v\cdot \mathbf{n}\right\}\mathrm{d}s+\sum _{e\in {\Gamma }_h\cup \partial \Omega }\frac{{\gamma }_0}{\rho }{\int }_e\left[u\right]\left[v\right]\mathrm{d}s,\end{array}$$

$$a(u_{\epsilon },v)=(f,v)\forall v\in H^s({\mathcal{T}}_h).$$

$$-\nabla \cdot ({\mathbf{a}}^{\epsilon }\nabla {\psi }_i^S)=0\text{ in }S,{\psi }_i^S{|}_{\partial S}={\phi }_i^S.$$

$${\overline{{\psi }_i}}^K=\sum _{j=1}^3{c}_{ij}^K{\psi }_j^S{|}_K\text{ in }K,$$

$${\phi }_i^K=\sum _{j=1}^3{c}_{ij}^K{\phi }_j^S{|}_K\text{ in }K.$$

$${\Pi }_K\psi =c_i{\phi }_i^K\text{ if }\psi =c_i{\bar{\psi }}_i^K\in \text{OMS}(K).$$

$$\begin{array}{rl}{V}_{h,dc}^{ms}& =\{{\psi }_h\in L^2(\Omega ):{\psi }_h{|}_K\in \text{OMS}(K)\forall K\in {\mathcal{T}}_h\},\\ V_{h,dc}& =\{v_h\in L^2(\Omega ):v_h{|}_K\in P_1(K)\forall K\in {\mathcal{T}}_h\}.\end{array}$$

$$a({\tilde{u}}_h,v_h)=(f,v_h)\forall v_h\in V_{h,dc}^{ms}.$$

$$\Pi_{h}\psi_{h}|_{K}=\Pi_{K}\psi_{h}\ \ \mbox{for any $K\in{\cal T}_{h},\psi_{h}\in V_{h,dc}^{ms}$}.$$

$$a_h(u_h,v_h):=$$

$$-\sum _{e\in {\Gamma }_h\cup \partial \Omega }{\int }_e\left\{{\mathbf{a}}^{\epsilon }\mathrm{NewA}{u}_h\cdot \mathbf{n}\right\}\left[{\Pi }_h{v}_h\right]\mathrm{d}s$$

$$+\beta \sum _{e\in {\Gamma }_h\cup \partial \Omega }{\int }_e\left[{\Pi }_h{u}_h\right]\left\{{\mathbf{a}}^{\epsilon }\mathrm{NewA}{v}_h\cdot \mathbf{n}\right\}\mathrm{d}s$$

$$+J_0(u_h,v_h),$$

$$J_0(u_h,v_h):=$$

$$a_h(u_h,v_h)=(f,{\Pi }_h{v}_h)\forall v_h\in V_{h,dc}^{ms}.$$

$$\begin{split}\left\|v_{h}\right\|_{h,\Omega}:=&\Big{(}\sum_{K\in\mathcal{T}_{h}}\int_{K}\mathbf{a}^{\varepsilon}{\rm NewA}v_{h}\cdot{\rm NewA}v_{h}\,\mathrm{d}x+\sum_{e\in\Gamma_{h}\cup\partial\Omega}\frac{\rho}{\gamma_{0}}\int_{e}\left\{\mathbf{a}^{\varepsilon}{\rm NewA}v_{h}\cdot\mathbf{n}\right\}^{2}\,\mathrm{d}s\\ &\qquad\qquad+\sum_{e\in\Gamma_{h}\cup\partial\Omega}\frac{\gamma_{0}}{\rho}\int_{e}\left[\Pi_{h}v_{h}\right]^{2}\,\mathrm{d}s\Big{)}^{1/2}.\end{split}$$

$$\begin{split}E(v,v_{h}):=&\bigg{(}\sum_{K\in\mathcal{T}_{h}}\|(\mathbf{a}^{\varepsilon})^{1/2}{\rm NewA}(v-v_{h})\|_{L^{2}(K)}^{2}\\ &+\sum_{e\in\Gamma_{h}\cup\partial\Omega}\frac{\rho}{\gamma_{0}}\left\|\left\{\mathbf{a}^{\varepsilon}{\rm NewA}(v-v_{h})\cdot\mathbf{n}\right\}\right\|_{L^{2}(e)}^{2}\\ &+\sum_{e\in\Gamma_{h}\cup\partial\Omega}\frac{\gamma_{0}}{\rho}\left\|\left[v-\Pi_{h}v_{h}\right]\right\|_{L^{2}(e)}^{2}\bigg{)}^{1/2}\quad\forall v\in H^{2}(\Omega),v_{h}\in V_{h,dc}^{ms}.\end{split}$$

$$E(v,v_h)\lesssim E(v,w_h)+{\Vert w_h-v_h\Vert }_{h,\Omega }.$$

$$u_1=u_0+\epsilon {\chi }^j(x/\epsilon )\frac{\partial u_0}{\partial x_j},$$

$$-{\nabla }_y\cdot (\mathbf{a}(y){\nabla }_y{\chi }^j(y))={\nabla }_y\cdot (\mathbf{a}(y){\mathbf{e}}_j),j=1,\cdots ,n$$

$${\left|u_{\epsilon }-u_1\right|}_{H^2(\Omega )}\lesssim {\left|u_0\right|}_{H^2(\Omega )}+\frac{1}{\sqrt{\epsilon }}{\left|u_0\right|}_{W^{1,\infty }(\Omega )}+\epsilon {\left|u_0\right|}_{H^3(\Omega )}.$$

$${\psi }_i^S={\phi }_i^S+\epsilon {\chi }^j(x/\epsilon )\frac{\partial {\phi }_i^S}{\partial x_j}+\epsilon {\eta }^j(x)\frac{\partial {\phi }_i^S}{\partial x_j},$$

$$-\mathrm{NewA}\cdot ({\mathbf{a}}^{\epsilon }\mathrm{NewA}{\eta }^j)=0\text{ in }S,{\eta }^j{|}_{\partial S}=-{\chi }^j(x/\epsilon ).$$

Full text

1 Introduction

This paper considers the numerical approximation of second order elliptic problems with heterogeneous and highly oscillating coefficients. These problems arise in many applications such as flows in porous media or composite materials. The numerical simulation of such problems in heterogeneous media poses major mathematical and computational challenges. Standard numerical methods such as the finite element method (FEM) or the finite volume method (FVM) usually require the mesh size very fine. This leads to tremendous amount of computer memory and CPU time. In the past several decades, a number of multiscale numerical methods have been proposed to solve these problems, see e.g., multiscale finite element method (MsFEM) [41, 51, 52], heterogeneous multiscale method (HMM) [33, 34, 35], upscaling or numerical homogenization method [18, 32, 46, 47], variational multiscale method (or the residual-free bubble method) [14, 16, 55, 54, 63], wavelet homogenization techniques [30, 44], and multigrid numerical homogenization techniques [48, 61]. Most of them are presented on meshes that are coarser than the scale of oscillations. The small scale effect on the coarse scale is either captured by localized multiscale basis functions or modeled into the coarse scale equations with prescribed analytical forms.

In this paper, the framework of the MsFEM is used to propose a new method. Two main ingredients of the MsFEM are the global formulation of the method such as various finite element methods and the construction of basis functions. The key point of MsFEM is to construct the multiscale basis from the local solutions of the elliptic operator for finite element formulation. There have been many extensions and applications of the method in the past fifteen years (cf. [1, 2, 3, 19, 22, 24, 25, 29, 38, 39, 53, 56]). We refer the reader to the book [40] for more discussions on the theory and applications of MsFEMs.

It is shown that the oversampling MsFEM is a nonconforming FEM where the numerical solution has certain continuity across the inner-element boundaries, while its basis functions are discontinuous at the inner-element boundaries (see [51, 41]). Note that the discontinuous Garlerkin (DG) FEMs do not ask for any continuity, which comes the natural idea that using the DG FEM as the global formulation coupling with the oversampling multiscale bases (see [40]). DG FEMs for elliptic boundary value problems have been studied since late 1970s, and it is now an active research area (see [31, 7, 8]). Examples of the DG methods include the Local Discontinuous Galerkin (LDG) method [4, 17, 62], and the interior penalty discontinuous Galerkin (IPDG) methods [7, 8, 9, 10, 11, 15, 23]. In this paper we are concerned with the IPDG method, still named DGM. DG methods admit good local conservation properties of the state variable and also offer the use of very general meshes due to the lack of inter-element continuity requirement, e.g., meshes that contain several different types of elements or hanging nodes. These features are crucial in many multiscale applications (see [28, 65]).

There have been shown several multiscale methods related with DG methods in the past ten years. For instance, a multiscale model reduction technique in the framework of the DG FEM for the high–contrast problems, named Generalized Multiscale Finite Element method, was presented in [36]. The special multiscale basis functions of the DG approximation space to capture the singularity of the solutions were discussed in [66, 67, 68]. The variational multiscale method based on the DGM for the elliptic multiscale problems without any assumption on scale separation or periodicity were proposed in [42, 43]. Heterogeneous multiscale method based on DGM for homogenization or advection–diffusion problems were presented in [5, 6]. However, to our knowledge, the multiscale discontinuous finite element method (MsDFEM) which coupling the classical oversampling multiscale basis with the discontinuous Galerkin method has not been studied detailedly in both aspects of error analysis and numerical tests in the literature. To fill this gap, in this paper we provide the formulation and the corresponding error estimate of the MsDFEM. Our numerical experiments show that MsDFEM takes advantages of MsFEM and DGM, which can eliminate the resonance error and obtain more accurate results than the classical MsFEM.

Further, we notice that the Petrov–Galerkin (PG) method can decrease the computational complexity significantly, allowing for more efficient solution algorithms. Moreover, it has been found that the MsPGM can eliminate the resonance error by using the oversampling technique [51] and the conforming piecewise linear functions as test functions [53, 52]. Therefore, in this paper we try to use the discontinuous oversampling multiscale space in the framework of PG method, which couples both DGM and MsPGM, named multiscale discontinuous Petrov–Galerkin method (MsDPGM). There are two key issues of the MsDPGM to consider. The first one is how to define its bilinear form and prove the coercive condition of the bilinear form, which needs to use the transfer operator between the approximation space and the test function space. We emphasize that, compared to the MsDFEM, the bilinear form of MsDPGM is not just choosing the discontinuous piecewise linear function space as test function space. More delicate choice of the terms of bilinear form should be made. The second one is the error estimate of MsDPGM. We give the $H^{2}$ norm error estimate between the multiscale solution and the homogenized solution with the first order corrector, which plays an important role in the error estimate. The MsDPGM takes advantages of the MsPGM and DGM, which is expected to better approximate the multiscale solution than the standard MsPGM.

The proposed method is related with a combined finite element and oversampling multiscale Petrov–Galerkin method (FE-OMsPGM) [65]. The idea of FE-OMsPGM is to utilize the traditional FEM directly on a fine mesh of the problematic part of the domain and use the OMsPGM on a coarse mesh of the other part. The transmission condition across the FE-OMsPGM interface is treated by the penalty technique of DGM. In [65], they deal with the transmission condition by penalizing the jumps from linear function values as well as the fluxes of the finite element solution on the fine mesh to those of the oversampling multiscale solution on the coarse mesh. Compared to [65], in this paper, we develop and analyze MsDPGM for the multiscale elliptic problems. The basic idea is to use PG formulation based on the discontinuous multiscale approximation space. The jump terms across each inter-element are dealt with penalty technique. The penalty term of linear function values is taken as that of the FE-OMsPGM, while the penalty term of the fluxes is not needed here.

Although the error analysis is given under the assumption that the oscillating coefficient is periodic, our method is not restrict to the periodic case. The numerical results show that the introduced MsDPGM is very efficient for randomly generated coefficients. Recently, the multiscale methods on localization of the elliptic multiscale problems with highly varying (non-periodic) coefficients are studied in some papers. For instance, the new variational multiscale method is presented in [59]; a new oversampling strategy for the MsFEM is presented in [50]. In the future work, we will give more extensions and developments on our method with the new oversampling strategy.

The outline of this paper is as follows. In Section 2, we give the model problem and recall the DG variational formulation of the model problem in the broken Sobolev spaces. Section 3 is devoted to derive the MsDPGM. It includes the introduction of discontinuous oversampling multiscale approximation space and the derivation of the formulations of MsDFEM and MsDPGM. In Section 4, we review the homogenization results and give some preliminaries for the error analysis. In Section 5, we give the main results of our method. It includes the stability and a priori error estimate of the proposed method. In Section 6, we first give several numerical examples with periodic coefficients to demonstrate the accuracy of the method. Then we do the experiment to study how the size of oversampling elements affects the errors. Finally, we apply our method to multiscale problems on the L–shaped domain to demonstrate the efficiency of the method. Conclusions are given in the last section.

2 Model problem and DG variational formulation

In this section we introduce the multiscale model problem and give the DG variational formulation of the model problem. First we state some notations and conventions. Throughout this paper, the Einstein summation convention is used: summation is taken over repeated indices. Standard notation on Lebesgue and Sobolev spaces is employed. Subsequently $C,C_{0},C_{1},C_{2},\cdot\cdot\cdot$ denote generic constants, which are independent of $\varepsilon,h$, unless otherwise stated. We also use the shorthand notation $A\lesssim B$ and $B\lesssim A$ for the inequality $A\leq CB$ and $B\leq CA$. The notation $A\eqsim B$ is equivalent to the statement $A\lesssim B$ and $B\lesssim A$.

2.1 Model problem

Let $\Omega\subset\mathbf{R}^{n},n=2,3$ be a bounded polyhedral domain. Consider the following multiscale elliptic problem:

[TABLE]

where $\epsilon\ll 1$ is a parameter that represents the small scale in the physical problem, $f\in L^{2}(\Omega)$, and $\mathbf{a}^{\varepsilon}(x)=(a_{ij}^{\varepsilon}(x))$ is a symmetric, positive definite matrix:

[TABLE]

for some positive constants $\lambda$ and $\Lambda$.

2.2 DG variational Formulation

In this subsection, we derive the DG variational formulation of the model problem in the broken Sobolev spaces. Let ${\cal T}_{h}$ be a quasi-uniform triangulation of the domain $\Omega$. We define $h_{K}$ as diam $(K)$ and denote by $h=\text{max}_{K\in{\cal T}_{h}}h_{K}$.

We introduce the broken Sobolev spaces for any real number $s$,

[TABLE]

equipped with the broken Sobolev norm:

[TABLE]

Denote by $\Gamma_{h}$ the set of interior edges/faces of the ${\cal T}_{h}$. With each edge/face $e$, we associate a unit normal vector $\mathbf{n}$. If $e$ is on the boundary $\partial\Omega$, then $\mathbf{n}$ is taken to be the unit outward vector normal to $\partial\Omega$.

If $v\in H^{1}(\mathcal{T}_{h})$, the trace of $v$ along any side of one element $K$ is well defined. If two elements $K_{1}^{e}$ and $K_{2}^{e}$ are neighbors and share one common side $e$, there are two traces of $v$ along $e$. We define the average and jump for $v$. We assume that the normal vector $\mathbf{n}$ is oriented from $K_{1}^{e}$ to $K_{2}^{e}$:

[TABLE]

We extend the definition of jump and average to sides that belong to the boundary $\partial\Omega$:

[TABLE]

In the following, we assume that $s=2$. Multiplying (2.1) by any $v\in H^{s}({\cal T}_{h})$, integrating on each element $K$, and using integration by parts, we obtain

[TABLE]

We recall that $\mathbf{n}_{K}$ is the outward normal to $K$. Summing over all elements, and switching to the normal vectors $\mathbf{n}$, we have

[TABLE]

From the regularity of the solution $u_{\varepsilon}$, it follows that

[TABLE]

where we have used the formula $\left[vw\right]=\left\{v\right\}\left[w\right]+\left[v\right]\left\{w\right\}$ and the fact that

$\left[\mathbf{a}^{\varepsilon}{\rm NewA}u_{\varepsilon}\cdot\mathbf{n}\right]=0$.

We now define the DG bilinear form $a(\cdot,\cdot):H^{s}({\cal T}_{h})\times H^{s}({\cal T}_{h})\rightarrow\mathbf{R}:$

[TABLE]

where $\beta$ is a real number such as $-1,0,1$, $\gamma_{0}$ is called penalty parameter, and $\rho>0$ will be specified later.

The general DG variational formulation of the problem (2.1) is as follows: Find $u_{\varepsilon}\in H^{s}({\cal T}_{h})$, such that

[TABLE]

Remark 2.1**.**

It is easy to check that if the solution $u_{\varepsilon}$ of problem (2.1) belongs to $H^{2}(\Omega)$, then $u_{\varepsilon}$ satisfies the variational formulation (2.4). Conversely, if $u_{\varepsilon}\in H^{1}(\Omega)\cap H^{s}({\cal T}_{h})$ satisfies (2.4), then $u_{\varepsilon}$ is the solution of problem (2.1).

3 Multiscale discontinuous Petrov-Galerkin method

This section is devoted to the formulations of multiscale discontinuous methods for solving (2.1). In subsection 3.1, we introduce the oversampling multiscale approximation space defined on the triangulation ${\cal T}_{h}$. The formulations of the MsDFE and MsDPG methods are presented in subsection 3.2.

3.1 Oversampling multiscale approximation space

In this subsection, we introduce the oversampling multiscale approximation space defined on the triangulation ${\cal T}_{h}$ (cf. [20, 40, 51]). Here we only consider the case where $n=2$. For any $K\in{\cal T}_{h}$ with nodes $\{x_{i}^{K}\}_{i=1}^{3}$, let $\{\varphi_{i}^{K}\}^{3}_{i=1}$ be the basis of $P_{1}(K)$ satisfying $\varphi_{i}^{K}(x_{j}^{K})=\delta_{ij}$, where $\delta_{ij}$ stands for the Kroneckers symbol. For any $K\in{\cal T}_{h},$ we denote by $S=S(K)$ a macro-element which contains $K$ and $d_{K}=\text{dist}(\partial S,K)$. We assume that $d_{K}\geq\delta_{0}h_{K}$ for some positive constant $\delta_{0}$ independent of $h_{K}$. The minimum angle of $S(K)$ is bounded from below by some positive constant $\theta_{0}$ independent of $h_{K}$. In our later numerical experiments, for any coarse-grid element $K\in\mathcal{T}_{h}$ we put its macro-element $S(K)$ in such a way that their barycenters are coincide and their corresponding edges are parallel. See Figure 2 for an illustration.

Let $\psi_{i}^{S},i=1,2,3,$ with $\psi_{i}^{S}\in H^{1}(S)$, be the solution of the problem:

[TABLE]

Here $\{\varphi_{i}^{S}\}_{i=1}^{3}$ is the nodal basis of $P_{1}(S)$ such that $\varphi_{i}^{S}(x_{j}^{S})=\delta_{ij},i,j=1,2,3.$

The oversampling multiscale basis functions on $K$ are defined by

[TABLE]

with the constants so chosen that

[TABLE]

The existence of the constants $c_{ij}^{K}$ is guaranteed because $\{\varphi_{j}^{S}\}_{j=1}^{3}$ also forms the basis of $P_{1}(K)$. To illustrate the basis functions, we depict two examples of them in Figure 1 (cf. [65]) .

Let $\text{OMS}(K)=\mathrm{span\,}\{\bar{\psi_{i}}^{K}\}_{i=1}^{3}$ be the set of space functions on $K$. Define the projection $\Pi_{K}:\text{OMS}(K)\rightarrow P_{1}(K)$ as follows:

[TABLE]

Further, we introduce the discontinuous piecewise “OMS” approximation space and the discontinuous piecewise linear space:

[TABLE]

Here we use the abbreviated indexes ‘ms’, ‘dc’ for multiscale, discontinuous, respectively.

3.2 Formulation of the MsDFEM/MsDPGM

In this subsection we present the formulations of the MsDFEM and MsDPGM.

By use of the DG variational formulation (2.4) and the discontinuous piecewise “OMS” approximation space, we are now ready to define the MsDFE method: Find $\widetilde{u}_{h}\in V_{h,dc}^{ms}$ such that

[TABLE]

To define the discrete bilinear form for MsDPGM, we need the transfer operator $\Pi_{h}:V_{h,dc}^{ms}\rightarrow V_{h,dc}$ as following:

[TABLE]

Remark 3.1**.**

In general, the trial and test functions of PGM are not in the same space. For example, here we might use $V_{h,dc}^{ms}$ and $V_{h,dc}$ as the trial function and test function spaces respectively. However, it may result in a difficulty to prove the inf-sup condition of the corresponding bilinear form. Hence, in this paper we introduce the transfer operator $\Pi_{h}$ to connect the two spaces, and use it in the bilinear form which causes an easy way to establish the stability of the MsDPGM. The idea of connecting the trial function and test function spaces in the Petrov-Galerkin method through an operator was introduced in [26] (see also [60]).

The discrete bilinear form of MsDPGM on $V_{h,dc}^{ms}\times V_{h,dc}^{ms}$ is defined as:

[TABLE]

where $\beta$ is a real number such as $-1,0,1$, $\gamma_{0}$ is the penalty parameter, and $\rho$ will be specified later.

Remark 3.2**.**

It is well known that DG methods utilize discontinuous piecewise polynomial functions and numerical fluxes, which implies that the weak formulation subject to discretization must include jump terms across interfaces and that some penalty terms must be added to control the jump terms. Therefore, the methods need the restriction on the penalty parameter to ensure stability and convergence in some sense. In fact, the optional convergence is related with the penalty parameter (see [45]). In this paper, our MsDPGM takes advantage of the penalty technique, which is inevitable to consider the choice of the penalty parameter. In our theoretical analysis, the penalty parameter $\gamma_{0}$ is constrained by a large constant from below to ensure the coercivity of $a_{h}$. However, in practice, the penalty parameter is chosen through our experience. In later numerical tests, we try different choice of the penalty parameter to study its affection on the error.

Remark 3.3**.**

The parameter $\rho$ is chosen as $\epsilon$ in our later error analysis, while in practical computation, it may be chosen as the mesh size $h$.

Then, our MsDPG method is: Find $u_{h}\in V_{h,dc}^{ms}$ such that

[TABLE]

Remark 3.4**.**

The design of the last two terms in $a_{h}$ is tricky. As a matter of fact, we have tried numerically different possibilities of using $\Pi_{h}$ ( or not before each $u_{h}$ or $v_{h}$ ) before we found that the current form of $a_{h}$ is the best one and, most importantly, the corresponding MsDPGM can be analyzed theoretically. Indeed, our MsDPGM is some kind of pseudo Petrov-Galerkin formulation of the method that the test function space is formally the same as the solution space, however some terms involve a projection of the multiscale test function into a piecewise linear function space.

We denote the discrete norm for MsDPGM on $V_{h,dc}^{ms}$,

[TABLE]

Noting that the operator $\Pi_{h}$ is not defined for the exact solution $u_{\varepsilon}$, we introduce the following function to measure the error of the discrete solution:

[TABLE]

From the triangle inequality, it is clear that, for any $v\in H^{2}(\Omega),v_{h},w_{h}\in V_{h,dc}^{ms}$,

[TABLE]

4 Homogenization results and preliminaries

In this section we first review the results of classical homogenization theory, and give an important $H^{2}$ norm error estimate between the multiscale solution and the homogenized solution with the first order corrector. Then we recall some preliminaries for our following analysis.

4.1 The homogenization results

Hereafter, we assume that $\mathbf{a}^{\varepsilon}(x)$ has the form $\mathbf{a}(x/\varepsilon)$ and $a_{ij}(y)$ are sufficiently smooth periodic functions in $y$ with respect to a unit cube $Y$. For our analysis it is sufficient to assume that $a_{ij}(y)\in W^{2,p}(Y)$ with $p>n$.

For convenience sake, we take

[TABLE]

where $u_{0}$ is the homogenized solution, $\chi^{j}$ is the periodic solution of the following cell problem (cf. [12, 57]):

[TABLE]

with zero mean, i.e., $\int_{Y}\chi^{j}dy=0$, and $\mathbf{e}_{j}$ is the unit vector in the $j$th direction.

The following theorem gives the $H^{2}$ semi-norm estimate of the error $u_{\varepsilon}-u_{1}$, which plays a key role in the later error analysis.We arrange the proof in the Appendix A.

Theorem 4.1**.**

Assume that $u_{0}\in H^{3}(\Omega)$. Then the following estimate is valid:

[TABLE]

4.2 Preliminaries

In this subsection, we give some preliminaries for our later analysis. We first recall the definition of $\psi_{i}^{S},i=1,2,3$ (see (3.1)). By the asymptotic expansion (cf. [41, 53]), we know that

[TABLE]

with $\eta^{j}$ being the solution of

[TABLE]

Substituting (4.3) to (3.2), we see that $\bar{\psi_{i}}^{K}$ can be expanded as follows:

[TABLE]

Recall that $d_{K}=\text{dist}(\partial S,K)$, which satisfies: $d_{K}\geq\delta_{0}h_{K}$. Denote by $d=\text{min}_{K\in\mathcal{T}_{h}}d_{K}$.

Remark 4.1**.**

It has been shown in [51, 53] that the distance $d_{K}$ is determined by the thickness of the boundary layer of $\eta^{j}$. Numerically, it has been observed that the boundary layer of $\eta^{j}$ is about $O(\epsilon)$ thick (see [51]). It was also observed that $d_{K}=h_{K}(>\epsilon)$ is usually sufficient for eliminating the boundary layer effect. Therefore in our numerical tests we choose $h_{K}$ as the oversampling size in this paper. To study how the size of oversampling elements affects the errors, in Section 6 we include a numerical test which use a series of $d_{K}$ with different $\delta_{0}$ to compare the corresponding errors.

By the Maximum Principle we have

[TABLE]

which together with the interior gradient estimate (see [28, Lemma 3.6] or [41, Proposition C.1]) imply that

[TABLE]

Next, we give a trace inequality which will be used in this paper frequently (see [13, Theorem 1.6.6], [27]).

Lemma 1**.**

Let $K$ be an element of the triangulation ${\cal T}_{h}$. Then, for any $v\in H^{1}(K)$, we have

[TABLE]

The following lemma gives some approximation properties of the space OMS(K) (cf. [28, Lemma 4.1]).

Lemma 2**.**

Take $\phi_{h}^{K}=\sum\limits_{x_{i}^{K}\,\text{node of}\ K}u_{0}(x_{i}^{K})\bar{\psi}_{i}^{K}(x),\quad\forall K\in{\cal T}_{h}$. Then, the following estimates hold:

[TABLE]

Moreover, we recall the stability estimate for $\Pi_{K}$, which will be used in our later analysis (cf. Lemma 3.2 in [65]).

Lemma 3**.**

There exist positive constants $\gamma$, $\alpha_{1}$ and $\alpha_{2}$ which are independent of $h$ and $\varepsilon$ such that if $\varepsilon/h_{K}\leq\gamma$ for all $K\in\mathcal{T}_{h}$, then the following estimates are valid for any $v_{h}\in\text{OMS}(K)$,

[TABLE]

The following lemma gives an inverse estimate for the function in space $\text{OMS}(K)$ [28, Lemma 5.2].

Lemma 4**.**

Under the assumptions of Lemma 3, and assuming $\varepsilon\lesssim h_{K}\lesssim d_{K}$, we have

[TABLE]

5 Main results

In this section we only carry out the convergence analysis of MsDPGM. For the case of MsDFEM, similar results can be obtained by the same argument and are arranged in the Appendix B for convenience of the reader. For MsDPGM, we first show the stability of the bilinear form guaranteeing the existence and uniqueness of the solution, and then prove the error estimate with $\beta=-1,\rho=\varepsilon$. For other cases such as $\beta=0,1$, the analysis is similar and is omitted here.

5.1 Existence and uniqueness of the solution of the MsDPGM

We start by establishing the stability of the bilinear form of the MsDPGM.

Theorem 5.1**.**

We have

[TABLE]

Further, let the assumptions of Lemma 4 be fulfilled and $\gamma_{0}$ is large enough, then

[TABLE]

where $\kappa>0$ is a constant independent of $h,\varepsilon,\gamma_{0}$.

Proof.

From the definition of the norms, the Cauchy-Schwarz inequality and Lemma 3, it follows (5.1) immediately.

Next we prove (5.2). From (4.13), we get

[TABLE]

It is easy to see that,

[TABLE]

Hence, we obtain

[TABLE]

By use of Lemmas 1, 4 and $\varepsilon\lesssim h$, we have

[TABLE]

Therefore, from (5.3) and (5.4), we have

[TABLE]

where $\gamma_{0}$ is large enough such that $\frac{5C_{1}}{2\gamma_{0}}<\frac{C}{2}$. Then by choosing $\kappa=\text{min}(\frac{C}{2},\frac{1}{2})$, it follows (5.2). This completes the proof. ∎

Theorem 5.1 guarantees that there exists a unique solution to our MsDPGM. Now we establish an analogue of Céa lemma written in the following theorem:

Theorem 5.2**.**

For large enough $\gamma_{0}$, the following inequality holds:

[TABLE]

where the error function $E$ is defined in (3.6).

Proof.

It is clear that by Theorem 5.1 we have

[TABLE]

From (2.4), it follows that

[TABLE]

Then, using the facts that $\left[u_{\varepsilon}\right]=0$ and $\left[\mathbf{a}^{\varepsilon}{\rm NewA}u_{\varepsilon}\cdot\mathbf{n}\right]=0$, we have

[TABLE]

Therefore, we obtain

[TABLE]

which together with (3.7) yield

[TABLE]

The proof is completed. ∎

5.2 A priori error estimate of MsDPGM

We present the main result of the paper which gives the error estimate of the MsDPGM.

Theorem 5.3**.**

Let $u_{\varepsilon}$ be the solution of (2.1), and let $u_{h}$ be the numerical solution computed using MsDPGM defined by (3.5). Assume that $u_{0}\in H^{3}(\Omega),f\in L^{2}(\Omega)$, and that $\varepsilon\lesssim h\lesssim d$, and that the penalty parameter $\gamma_{0}$ is large enough. Then there exits a constant $\gamma$ independent of $h$ and $\varepsilon$ such that if $\varepsilon/h_{K}\leq\gamma$ for all $K\in\mathcal{T}_{h}$, the following error estimate holds:

[TABLE]

where $d=\text{min}_{K\in\mathcal{T}_{h}}d_{K}$.

Proof.

According to Theorem 5.2, the proof is devoted to estimating the interpolation error. To do this, we define $\psi_{h}$ by

[TABLE]

Clearly, $\psi_{h}\in V_{h,dc}^{ms}$. It is easy to see that

[TABLE]

where $I_{h}:H^{s}(\mathcal{T}_{h})\rightarrow V_{h,dc}$ is the Lagrange interpolation operator. Then we set $v_{h}$ as $\psi_{h}$. It is shown that in [28],

[TABLE]

Next, we estimate the term

[TABLE]

From (4.8), we have

[TABLE]

Therefore, it follows from Theorem 4.1, Lemma 2 and the assumption $\varepsilon\lesssim h\lesssim d$ that,

[TABLE]

where we have used $\frac{\varepsilon}{\sqrt{h}}<\sqrt{\varepsilon}$ and the Young’s inequality to derive the above inequality.

It remains to consider the term $\sum_{e\in\Gamma_{h}\cup\partial\Omega}\frac{\gamma_{0}}{\varepsilon}\left\|\left[u_{\varepsilon}-\Pi_{h}v_{h}\right]\right\|_{L^{2}(e)}^{2}.$ Noting that both $u_{\varepsilon}$ and $u_{0}$ are continuous functions, we have

[TABLE]

Then, by use of Lemma 1, we have

[TABLE]

which yields

[TABLE]

Hence, from (5.8)–(5.10), it follows (5.6) immediately. ∎

6 Numerical experiments

In this section, we present numerical experiments to confirm the theoretical results in Section 5. We show the numerical results of MsDPGM defined in (3.5), and also results of MsDFEM defined in (3.4) which show good performance as well as MsDPGM. In order to illustrate the accuracy of our methods, we also implement the standard MsFEM in Petrov–Galerkin formulation which is denoted as MsPGM, and the MsPGM which uses the classical oversampling multiscale basis (OMsPGM). We also show the results of the traditional linear finite element method (FEM) and discontinuous finite element method (DFEM) on the corresponding coarse grid to get a feeling for the accuracy of the multiscale methods. All numerical experiments are designed to show better performance of MsDPGM than the other MsPGMs.

In all tests, for simplicity, we use the standard triangulation which is constructed by first dividing the domain $\Omega$ into sub-squares of equal length $h$ and then connecting the lower-left and the upper-right vertices of each sub-square. For any coarse-grid element $K\in\mathcal{T}_{h}$ we put its macro-element $S(K)$ in such a way that their barycenters are coincide and their corresponding edges are parallel. The length of the horizontal and vertical edges of $S(K)$ is four times of the corresponding length of the edges of $K$. We assume that all right-angle sides of $S(K),K\in\mathcal{T}_{h}$ have the same length denoted by $h_{S}$. Recall the definition of the $d=\text{min}_{K\in\mathcal{T}_{h}}d_{K}$, define

[TABLE]

It is clear that $d\eqsim\tilde{d}$. See Figure 2 for an illustration.

In all of these computations, we have used finely resolved numerical solutions obtained using the traditional linear finite element method with mesh size $h_{f}=1/4096$ as the “exact” solutions which are denoted as $u_{e}$. Denoting $u_{h}$ as the numerical solutions computed by the methods considered in this section, we measure the relative error in $L^{2}$, $L^{\infty}$ and energy norms as following:

[TABLE]

where

[TABLE]

In all tests, the coefficient $\mathbf{a}^{\varepsilon}$ is chosen as the form $\mathbf{a}^{\varepsilon}=a^{\varepsilon}I$ where $a^{\varepsilon}$ is a scalar function and $I$ is the 2 by 2 identity matrix.

6.1 Application to elliptic problems with highly oscillating coefficients

We first consider the model problem (2.1) in the squared domain $\Omega=[0,1]\times[0,1]$. Assume that $f=1$ and the coefficient $\mathbf{a}^{\varepsilon}(x_{1},x_{2})$ has the following periodic form:

[TABLE]

where we fix $\varepsilon=1/100$.

In this test, we choose $h=1/32$ and report errors in the $L^{2},L^{\infty}$ and energy norms in Table 1. We can see that the MsDPGM and MsDFEM give more accurate results than the other multiscale methods considered here, while the FEM and DFEM give worse approximations to the gradient of solution. We also compare the CPU time $T_{1}$ and $T_{2}$ spent by the MsDFEM and MsDPGM to show the good performance of MsDPGM in computational complexity, where $T_{1}$ is the CPU time of assembling the stiffness matrix, and $T_{2}$ is the CPU time of solving the discrete system of algebraic equations. We can observe that the CPU time $T_{1}$ of our MsDPGM for assembling the stiffness matrix is shorter than that of MsDFEM since the Petrov–Galerkin method can decrease the computational complexity.

Secondly, we do an experiment to study how the penalty parameter $\gamma_{0}$ affects the errors. We fix $\rho=\varepsilon=1/100,\tilde{d}=h=1/32$ and choose a series of $\gamma_{0}$ in the test. The result is shown in Table 2. We observe that as $\gamma_{0}$ goes larger, the relative error is close to the error of the OMsPGM. It seems that MsDPGM converges to OMsPGM as the penalty parameter $\gamma_{0}$ goes to infinity (cf. [58]).

The third numerical experiment is to show the mesh size $h$ plays a role as that describing in Theorem 5.3. We fix $\tilde{d}=1/32$ and $\varepsilon=1/100$. Four kinds of mesh size are chosen: $h=1/64,1/32,1/16,1/8$. The results are shown in Table 3. Relative error in energy norm against the mesh size $h$ is clearly shown in Figure 3. It is easy to see that as $h$ goes larger, the relative error in energy norm goes larger, which is in agreement with the theoretical results in Theorem 5.3. We remark that the classical MsFEM suffers from the resonance error since the $H^{1}$–error estimate has the term $\epsilon/h$ due to the nonconforming error (see [52]). But for MsDPGM, the error estimate in Theorem 5.3, and the numerical results in Table 3 and Figure 3 show that the resonance error has been removed completely.

6.2 Affection of the size of the oversampling patches

In this subsection we study how the size of oversampling elements affects the error. The experiment to verify the inequality (4.7) for the model example with coefficient (6.2) has been done in [65]. The figures have been shown that $\left\|{\rm NewA}\eta^{j}\right\|_{L^{\infty}(K)}\cdot d_{K}$ are bounded by a constant which is consistent with (4.7) (see Figure 5 in [65]).

The following numerical experiment is to show how the oversampling size affects the error. Recalling the requirement of the oversampling size $d_{K}\geq\delta_{0}h_{K}$, we show the relative oversampling size $\delta_{0}$ against the error. Note the distance $d=\text{min}_{K\in\mathcal{T}_{h}}d_{K}$. Here it is equivalent to $d\geq\delta_{0}h$. We set $\rho=\varepsilon=1/100,h=1/32$. The result is shown in Table 4. We can see that as $\delta_{0}$ (equivalently $\tilde{d}$) goes larger, the relative error in energy norm goes smaller, which is coincided with the theoretical results in Theorem 5.3. We also notice that when $\tilde{d}$ is close to $\sqrt{\varepsilon}$, the errors begin to decrease very slowly. Recall that there is a homogenization error $\sqrt{\varepsilon}$ in the error estimate (5.6). We think that when $d$ is large enough, $\sqrt{\varepsilon}$ becomes the dominated error instead of $\varepsilon/d$.

6.3 Application to multiscale problems on L–shape domain

We consider the multiscale problem on the L–shaped domain of Figure 4 with Dirichlet boundary condition so chosen that the true solution is $u=r^{\frac{1}{3}}\sin(2\theta/3)$ in polar coordinates. It is known that the solution has the singular behavior around reentrant corners. So the classical finite element method fails to provide satisfactory result.

Firstly, we simulate the problem with coefficient given by (6.2). We fix $\varepsilon=1/100$ and choose $h=1/16$. The relative error is shown in Table 5. We observe that both MsDPGM and MsDFEM give better approximation than the other MsPG methods.

Secondly, we simulate the problem with the random log-normal permeability field $\mathbf{a}(x)$, which is generated by using the moving ellipse average [32] with the variance of the logarithm of the permeability $\sigma^{2}=1.0$, and the correlation lengths $l_{1}=l_{2}=0.01$ in $x_{1}$ and $x_{2}$ directions, respectively. One realization of the resulting permeability field is depicted in Figure 5, where $\frac{a_{\rm max}(x)}{a_{\rm min}(x)}=2.9642e+003$. In this test, we set $\rho=h=1/16$ since there is no explicit $\varepsilon$ in the example. The result is shown in Table 6. We can see that MsDPGM gives a better approximation than the other MsPG methods, while standard MsPGM gives the wrong approximation to the gradient of solution.

7 Conclusion

In this paper, we have proposed a new Petrov–Galerkin method based on the discontinuous multiscale approximation space for the multiscale elliptic problems. Under some assumptions on the coefficients, we give the error analysis of our method. The $H^{1}$–error is of the order

[TABLE]

which consists of the oversampling multiscale approximation error and the error contributed by the penalty. Note that the unpleasant resonance error does not appear. The reason is that our method uses discontinuous piecewise linear functions as test functions, which is only needed to estimate the interpolation error. Several numerical experiments have demonstrated the efficiency of MsDPGM. We also study the corresponding MsDFEM which coupling the classical oversampling multiscale basis with DGM. Our convergence analysis shows that MsDFEM can also eliminate the resonance error completely. That is the reason why MsDFEM is working as well as MsDPGM (even a little better). Furthermore, we can see that the CPU-time cost of MsDPGM for assembling the stiffness matrix is shorter than that of the MsDFEM due to its PG version. Therefore, we think that MsDPGM is a good choice when we need to take into consideration of the computational accuracy and the computer resource at the same time.

We emphasize that the proposed method is not restrict to the periodic case. The numerical experiments show that it is applicable to the random coefficient case very well. However, with the classical oversampling multiscale basis function space introduced in [51], the error estimate method is based on the classical homogenization theory, which needs the assumption that the oscillating coefficient is periodic. In the future work, we would like to combine the Petrov–Galerkin method with the new oversampling multiscale space [50] to consider the elliptic multiscale problems without any assumption on scale separation or periodicity. Besides, the introduced method may be inefficient for the multiscale problems which have some singularities, such as, the Dirac function singularities which stems from the simulation of steady flow transport through highly heterogeneous porous media driven by extraction wells [21], or high-conductivity channels that connect the boundaries of coarse-grid blocks [37]. To solve these problems, it needs some special definition of the multiscale basis functions around the channels such as the local spectral basis functions (see [37]), or local refinement of the elements near the channels (see [28]). We will couple these techniques with the introduced method in our future work. Finally, we remark that Generalized Multiscale Finite Element method coupling DGM was explored in [36]. The computation is divided into two stages: offline and online. In the offline stage, they construct a reduced dimensional multiscale space to be used for rapid computations in the online stage. In the online stage, they use the basis functions computed offline to solve the problem for current realization of the parameters. Similar to MsDPGM, in the online stage we can use the Petrov–Galerkin version of DGM to solve the problem with the basis functions computed offline, which leads to a kind of Generalized Multiscale Discontinuous Petrov-Galerkin method. The difficulty is the choice of the test function space and the proof of inf-sup condition, which is worth studying.

Acknowledgments

The authors would like to thank the referees for their carefully reading and constructive comments that improved the paper.

Appendix A Proof of Theorem 4.1

The following theorem plays an important role in our analysis (cf. [19, 20]).

Theorem A.1**.**

Assume that $u_{0}\in H^{2}(\Omega)\cap W^{1,\infty}(\Omega)$. There exists a constant $C$ independent of $u_{0},\varepsilon,\Omega$ such that

[TABLE]

where $\theta_{\varepsilon}$ denote the boundary corrector defined by

[TABLE]

We first estimate $\left|\varepsilon\theta_{\varepsilon}\right|_{H^{2}(\Omega)}$.

Lemma 5**.**

Assume that $u_{0}\in H^{2}(\Omega)\cap W^{1,\infty}(\Omega)$. Then the following estimate holds:

[TABLE]

Proof.

We only consider the case where $n=2$. For $n=3$, the proof is similar. Let $\xi\in C_{0}^{\infty}(\mathbf{R}^{2})$ be the cut-off function such that $0\leq\xi\leq 1,\xi=1$ in $\Omega\setminus\Omega_{\varepsilon/2},\xi=0$ in $\Omega_{\varepsilon}$, and $\left|{\rm NewA}\xi\right|\leq C/\varepsilon,\left|{\rm NewA}^{2}\xi\right|\leq C/\varepsilon^{2}$ in $\Omega$, where $\Omega_{\varepsilon}:=\{x:\text{dist}\{x,\partial\Omega\}\geq\varepsilon\}$. Then

[TABLE]

satisfies

[TABLE]

By use of Theorem 4.3.1.4 in [49] , and together with Theorem A.1, we have

[TABLE]

which implies

[TABLE]

This completes the proof. ∎

Proof of Theorem 4.1. It is shown that, for any $\varphi\in H_{0}^{1}(\Omega)$ (see [19, p.550] or [20, p.125]),

[TABLE]

where $\alpha^{k}(x/\varepsilon)=(\alpha_{ij}^{k}(x/\varepsilon))$ are skew-symmetric matrices which satisfy that (see [57, p.6])

[TABLE]

with

[TABLE]

From (A.5), it follows that,

[TABLE]

which combines the definition of $\theta_{\varepsilon}$ yield

[TABLE]

Thus, from Theorem 4.3.1.4 in [49], it follows that

[TABLE]

which combing (A.2) and Theorem A.1, yield (4.2) immediately.

Appendix B Theoretical Results of MsDFEM

We give some theoretical results of MsDFEM here for convenience of the reader. Detailed analysis can be found in the first author’s PHD thesis [64].

Lemma 6**.**

We have

[TABLE]

Further, let the assumptions of Lemma 4 be fulfilled and $\gamma_{0}$ is large enough, then

[TABLE]

Here

[TABLE]

Using the definition of the above norm, the Cauchy-Schwarz inequality and (4.8), Lemma 4, we can obtain (B.1) and (B.2) immediately. The proof is similar to Theorem 5.1 and is omitted here.

Theorem B.1**.**

Let $u_{\varepsilon}$ be the solution of (2.1), and let $\tilde{u}_{h}$ be the numerical solution computed by MsDFEM defined in (3.4). Assume that $u_{0}\in H^{3}(\Omega),f\in L^{2}(\Omega)$, $\varepsilon\lesssim h\lesssim d$, and that the penalty parameter $\gamma_{0}$ is large enough. Then there exits a constant $\gamma$ independent of $h$ and $\varepsilon$ such that if $\varepsilon/h_{K}\leq\gamma$ for all $K\in\mathcal{T}_{h}$, the following error estimate holds:

[TABLE]

where $d=\text{min}_{K\in\mathcal{T}_{h}}d_{K}$.

Proof.

By use of the Galerkin orthogonality of $a(\cdot,\cdot)$, we only need to estimate the interpolation error.

Take $v_{h}$ as $\psi_{h}$ (see (5.7)). The following two estimates of the error have been shown in the proof of Theorem 5.3:

[TABLE]

and

[TABLE]

It remains to consider the term $\sum_{e\in\Gamma_{h}\cup\partial\Omega}\frac{\gamma_{0}}{\varepsilon}\left\|\left[u_{\varepsilon}-v_{h}\right]\right\|_{L^{2}(e)}^{2}$. Noting that $\left[u_{\varepsilon}\right]=\left[u_{1}\right]=0$, then by use of the trace inequality (4.8) and Lemma 2, we have

[TABLE]

where we have used the assumption $\varepsilon\lesssim h\lesssim d$ and the Young’s inequality to derive the above inequality.

Hence, from (B.4), (B.5) and (B.6), it follows (B.3) immediately. This completes the proof. ∎