Computational multiscale methods for parabolic wave approximations in heterogeneous media

Eric Chung , Yalchin Efendiev , Sai-Mang Pun , and Zecheng Zhang Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong KongDepartment of Mathematics and Institute of Scientific Computing, Texas A&M University, College Station, TX 77843, USADepartment of Mathematics, Texas A&M University, College Station, TX 77843, USADepartment of Mathematics, Purdue University, West Lafayette, IN 47906, USA

Abstract

In this paper, we develop a computational multiscale to solve the parabolic wave approximation with heterogeneous and variable media. Parabolic wave approximation is a technique to approximate the full wave equation. One benefit of the method is that: one wave propagation direction can be taken as an evolution direction, and we then can discretize it using a classical scheme like Backward Euler. Consequently, we obtain a set of quasi-gas-dynamic (QGD) models with different heterogeneous permeability fields. Then, we employ constraint energy minimization generalized multiscale finite element method (CEM-GMsFEM) to perform spatial discretization for the problem. The resulting system can be solved by combining the central difference in time evolution. Due to the variable media, we apply the technique of proper orthogonal decomposition (POD) to further the dimension of the problem and solve the corresponding model problem in the POD space instead of in the whole multiscale space spanned by all possible multiscale basis functions. We prove the stability of the full discretization scheme and give the convergence analysis of the proposed approximation scheme. Numerical results verify the effectiveness of the proposed method.

1 Introduction

Parabolic wave approximation have been used to approximate wave equations with a preferred direction [2, 3]. To be more specific, we study the approximation of the full wave equation:

\displaystyle\rho\partial_{tt}u-\nabla\cdot(\mu\nabla u)=0.

(1)

Here, $\nabla$ denotes the gradient operator in $\mathbb{R}^{3}$ ; $\rho$ and $\mu$ are positive functions in $\mathbb{R}^{3}$ . We then consider the parabolic approximation of (1) in $D\times\Omega\subseteq\mathbb{R}^{3}$ : for $(z,x)\in D\times\Omega$ with the boundary $\Gamma:=\partial(D\times\Omega)$ , find $v=\sigma^{1/2}u$ such that

\displaystyle\begin{split}c^{-1}\partial_{tt}v+\partial_{t}(\partial_{z}v)-\frac{1}{2}\nabla_{x}\cdot\left(c\nabla_{x}v\right)&=0&\quad\text{in }D\times\Omega\times(0,T],\\ v(z,x,0)&=v_{0}(z,x)&\quad\text{in }D\times\Omega,\\ \partial_{t}v(z,x,0)&=v_{1}(z,x)&\quad\text{in }D\times\Omega,\\ v(z,x,t)&=g(z,x,t)&\quad\text{in }\Gamma\times(0,T].\end{split}

(2)

Here, $\nabla_{x}$ denotes the gradient operator defined in the bounded domain $\Omega\subset\mathbb{R}^{2}$ ; $D=[d_{1},d_{2}]\subset\mathbb{R}$ is a bounded domain; $v_{0}$ and $v_{1}$ are initial conditions; $g$ is boundary condition and $T>0$ is a given terminal time.

Wave equations of this type have a propagation direction $z$ which plays a role of time [4]. Many real-world problems can be solved by the parabolic approximation; in particular in the areas of geology [22], under water acoustics [5, 6, 43, 44] and optics [23, 32, 36]. There are two benefits for using the parabolic approximation instead of the full wave equation [2, 3]: (i) it is easier to be realized and computationally more efficient; and (ii) it makes it possible to use the approximation as an evolution equation in $z$ direction. The second property makes it possible to solve the problem as an equation evolving in $z$ direction; we hence can discretize ( $\partial_{z}v$ ) by applying classical difference scheme. Consequently, we obtain a quasi-gas-dynamic (QGD) equation which have been thoroughly studied in literature [11, 12, 13].

It is relatively common in the geological problems [22] that the medium under consideration is highly variable and heterogeneous; that is, the function field $c$ is non-homogeneous given a cross section in $z$ and is fast-changing in $z$ direction. This brings in two difficulties. The first to mention is the multiscale property brought by the heterogeneous field. The second difficulty is the intense work in solving a QGD model given a $z$ cross-section in the $z$ evolution. Throughout the work, we study (2) with variable and heterogeneous media. We will derive scheme in discretizing the time and $z$ evolution and also provide our approach to solve the problems mentioned before.

As we have discussed, the model has heterogeneous in $z$ direction. Directly solving this problem on fine mesh can capture the multiscale features; however, this is computationally intense and this issue becomes exacerbate when people are solving time-dependent problems. Therefore, many methods which solve the multiscale problems on coarser mesh have been proposed. These include homogenization-based approaches [8, 9, 10, 29, 40, 41, 47], multiscale finite element methods [31, 35, 37, 38], generalized multiscale finite element methods (GMsFEM) [14, 15, 18, 19, 21, 25, 30], constraint energy minimizing GMsFEM (CEM-GMsFEM) [16, 17], nonlocal multi-continua (NLMC) approaches [20], metric-based upscaling [45], heterogeneous multiscale method [1, 24], localized orthogonal decomposition (LOD) [33, 42], equation free approaches [46, 48, 49], computational continua [26, 27, 28], hierarchical multiscale method [7, 34, 50], and so on. Some of these approaches, such as homogenization-based approaches, are designed for problems with scale separation. In this work, we apply the CEM-GMsFEM [16, 17] and provide the convergence analysis of our proposed scheme based on the coarse mesh convergence results of the CEM-GMsFEM.

The second difficulty of the problem is the variable media. If we discretize $z$ evolution using some classical difference scheme, each $z$ level is a QGD model with heterogeneous. This model can be solved in the framework of CEM-GMsFEM [11]; however the coarse scale basis evaluation is time consuming; in particular, this process will be repeated for each level of $z$ . We hence proposed the proper orthogonal decomposition (POD) technique [39] which is target to find low dimensional subspace such that the error of the orthogonal projection is minimized in the sense of the norm induced by the inner product of the original space. To be more specific, the proposed method can be summarized as follows:

We proposed a method to solve a parabolic-wave model with variable and heterogeneous media. We first apply the backward Euler scheme to discretization the term $\partial_{z}v$ . This leads to a set of two dimensional (space) QGD models (3); we call this level of discretization the quasi-time scheme and an unconditional stability result is established. To further discretize the problem, we apply the central scheme to deal with the time derivative $v_{tt}$ , $v_{t}$ ; and then use CEM-GMsFEM in space on coarse scale to capture the heterogeneous brought by the media. We then prove that the full discretization scheme is stable in an energy norm under some CFL condition by using an inverse inequality in the multiscale space.

The key of the CEM-GMsFEM method is to construct the CEM basis. The standard procedure is first to build the auxiliary multiscale basis by solving local spectral problems in coarse mesh; we then can construct the CEM basis by evaluating a set of energy minimization problems. Due to the variable media, we need to construct a set of basis for each QGD model and the solution of the full discretized scheme is in the space of all multiscale basis. This is time consuming; and we hence can apply the POD technique to find the best set of orthogonal basis in the sense of $L_{2}$ minimization; that is, the projection error of the full discretized solution onto the POD basis is optimal in the norm induced by the original space. In practice, we collect CEM basis for some QGD models and then construct POD basis of the space spanned by all CEM basis; the POD models will finally be solved by using POD basis. A convergence analysis of the POD approximation is established and the numerical results prove the algorithm is successful.

The remainder of the paper is organized as follows. In Section 2, we present some preliminaries of the model problem and briefly overview the framework of proper orthogonal decomposition. Section 3 is devoted to the multiscale methods and we will briefly overview the construction of multiscale basis function within the framework of the CEM-GMsFEM. In Section 4, we present a complete analysis of the proposed computational multiscale method. We then present some numerical results to demonstrate the efficiency of the proposed method in Section 5. Concluding remarks are drawn in Section 6.

2 Preliminaries

In this section, we present some preliminaries of the model problem. For simplicity, we assume the homogeneous boundary condition $g=0$ is equipped in the model problem (2). The extension of inhomogeneous case is straightforward. The system (2) is the one we shall consider throughout the remainder of this paper and we shall develop computational multiscale method for efficiently simulating the problem. We remark that under appropriate regularity assumptions on initial and boundary conditions, the problem (2) has a unique solution such that

t\to v(z,x,t)\in W^{1,\infty}(0,T;L^{2}(D\times\Omega))\cap L^{\infty}(0,T;L^{2}(D\times\Omega)),

z\to v(z,x,t)\in L^{\infty}(D,H^{1}(0,T;L^{2}(\Omega))).

The last property enables us to consider $z$ as an evolution direction. The result is an application of the semigroup theory and the Hille-Yoshida theorem. See [3, Section 4] for more details.

Instead of the PDE formulation (2), we consider its corresponding variational formulation. In the following, we treat $z$ as an evolution direction and use backward Euler method to discretize the term $\partial_{z}v$ . We divide the domain $D=[d_{1},d_{2}]$ along $z$ -direction into $K$ pieces. We write, for $k=0,1,\cdots,K$ ,

v_{k}=v(z_{k})\quad\text{with}\quad z_{k}=d_{1}+k\Delta z\quad\text{and}\quad d_{2}=d_{1}+K\Delta z.

The quasi-time discretization of (2) reads: find $\theta_{k}\in V$ for $k=1,\cdots,K-1$ such that

\displaystyle\left(\ddot{\theta}_{k},w;z_{k}\right)_{c}+\frac{1}{\Delta z}\left(\dot{\theta}_{k}-\dot{\theta}_{k-1},w\right)+\frac{1}{2}a(\theta_{k},w;z_{k})=0\quad\text{ for all }w\in V:=H_{0}^{1}(\Omega)

(3)

or equivalently

\Delta z\left(\dot{\theta}_{k},w;z_{k}\right)_{c}+\left(\dot{\theta}_{k},w\right)+\frac{\Delta z}{2}a(\theta_{k},w;z_{k})=\left(\dot{\theta}_{k-1},w\right)\quad\text{ for all }w\in V,

where $(v,w;z_{k})_{c}:=\int_{\Omega}c^{-1}(z_{k})vw\leavevmode\nobreak\ dx$ , $(v,w):=\int_{\Omega}vw\leavevmode\nobreak\ dx$ , and $a(v,w;z_{k}):=\int_{\Omega}c(z_{k})\nabla v\cdot\nabla w\leavevmode\nobreak\ dx$ . Here, we denote $\dot{v}:=\partial_{t}v$ and $\ddot{v}:=\partial_{tt}v$ . We assume that $c\in L^{\infty}(D\times\Omega)$ and we denote $c_{\max}:=\left\|c\right\|_{L^{\infty}(D\times\Omega)}$ . Employing Galerkin’s method and the method of energy estimate, one can show the well-posedness of the variational formulation (3). Denote $\left\|v\right\|:=\sqrt{(v,v)}$ , $\left\|v\right\|_{c(z_{k})}:=\sqrt{(v,v;z_{k})_{c}}$ , and $\left\|v\right\|_{a(z_{k})}:=\sqrt{a(v,v;z_{k})}$ . We establish the following stability estimate for the quasi-time discretization (3).

Lemma 2.1.

Let $\{v_{k}\}_{k=0}^{K}\subseteq V$ solve the equation (3). Then, the following stability estimate holds:

\left\|\dot{v_{K}}\right\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\Delta z\sum_{k=1}^{K}\mathcal{E}_{k}(v_{k}(T))\lesssim\left\|\dot{v_{0}}\right\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\Delta z\sum_{k=1}^{K}\mathcal{E}_{k}(v_{k}(0)),

where we denote $\mathcal{E}_{k}(v):=\left\|\dot{v}\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v\right\|_{a(z_{k})}^{2}$ .

Remark.

The term $\mathcal{E}_{k}(v_{k}(0))=\left\|v_{1}(z_{k})\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{0}(z_{k})\right\|_{a(z_{k})}^{2}$ on the right-hand side of the inequality above is determined by the initial conditions.

Proof of Lemma 2.1.

Taking $w=\dot{v_{k}}$ , we have

\displaystyle\begin{split}&\quad\frac{\Delta z}{2}\frac{d}{dt}\left(\left\|\dot{v}_{k}\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}\right\|_{a(z_{k})}^{2}\right)+\left\|\dot{v}_{k}\right\|^{2}-(\dot{v}_{k-1},\dot{v}_{k})\\ &=\Delta z\left(\ddot{v}_{k},\dot{v}_{k};z_{k}\right)_{c}+(\dot{v}_{k}-\dot{v}_{k-1},\dot{v_{k}})+\frac{\Delta z}{2}a(v_{k},\dot{v_{k}};z_{k})=0.\end{split}

Thus, we have

\quad\frac{\Delta z}{2}\frac{d}{dt}\left(\left\|\dot{v}_{k}\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}\right\|_{a(z_{k})}^{2}\right)+\left\|\dot{v}_{k}\right\|^{2}=(\dot{v}_{k},\dot{v}_{k-1})\leq\left\|\dot{v}_{k}\right\|\left\|\dot{v}_{k-1}\right\|.

Multiplying by $2$ and integrating over $(0,T]$ , we have

\displaystyle\begin{split}&{\Delta z}\left(\left\|\dot{v}_{k}(T)\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}(T)\right\|_{a(z_{k})}^{2}\right)+2\int_{0}^{T}\left\|\dot{v}_{k}\right\|^{2}dt\\ &\leq 2\int_{0}^{T}\left\|\dot{v}_{k}\right\|\left\|\dot{v}_{k-1}\right\|dt+\Delta z\left(\left\|\dot{v}_{k}(0)\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}(0)\right\|_{a(z_{k})}^{2}\right)\\ &\leq\int_{0}^{T}\left\|\dot{v}_{k}\right\|^{2}dt+\int_{0}^{T}\left\|\dot{v}_{k-1}\right\|^{2}dt+\Delta z\left(\left\|\dot{v}_{k}(0)\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}(0)\right\|_{a(z_{k})}^{2}\right).\\ \end{split}

Therefore, we have

\displaystyle\begin{split}\left\|\dot{v}_{k}\right\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}&+{\Delta z}\left(\left\|\dot{v}_{k}(T)\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}(T)\right\|_{a(z_{k})}^{2}\right)\\ &\leq\left\|\dot{v}_{k-1}\right\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\Delta z\left(\left\|\dot{v}_{k}(0)\right\|_{c(z_{k})}^{2}+\frac{1}{2}\left\|v_{k}(0)\right\|_{a(z_{k})}^{2}\right).\end{split}

Summing over $k=1,\cdots,K$ , we obtain the desired result and this completes the proof. ∎

2.1 The Proper Orthogonal Decomposition

In this section, we briefly introduce the proper orthogonal decomposition (POD) method. This method aims to generate optimally ordered orthogonal basis functions in the least squares sense for a given set of theoretical, experimental, or computational data. Reduced-order models or surrogate models are then obtained by truncating this set of optimal basis functions, providing considerable computational savings over the original high-dimensional problems.

Let $X$ be a real Hilbert space endowed with inner product $(\cdot,\cdot)_{X}$ and norm $\left\|\cdot\right\|_{X}$ . We set $\mathcal{V}:=\operatorname{span}\{y_{1},y_{2},\cdots,y_{n}\}$ with each $y_{i}\in X$ for $i\in\{1,\cdots,n\}$ . We refer to $\mathcal{V}$ as ensemble consisting of the snapshots $\{y_{i}\}_{i=1}^{n}$ , at least one of which is assumed to be non-zero. Let $\{\psi_{k}\}_{k=1}^{\mathcal{N}}$ denote a set of orthonormal basis functions of $\mathcal{V}$ with $\mathcal{N}:=\text{dim}(\mathcal{V})\leq n$ . Then, each member of the ensemble can be expressed as

y_{j}=\sum_{k=1}^{\mathcal{N}}(y_{j},\psi_{k})_{X}\psi_{k}

for each $j\in\{1,\cdots,n\}$ . The POD method consists in choosing the orthonormal basis functions such that for every $\ell\in\{1,\cdots,\mathcal{N}\}$ the mean square error between the elements $y_{j}$ (for any $j\in\{1,\cdots,n\}$ ), and the corresponding $\ell$ -th partial sum is minimized on average:

\displaystyle\begin{split}&\min_{\{\psi_{k}\}_{k=1}^{\ell}}\frac{1}{n}\sum_{j=1}^{n}\left\|y_{j}-\sum_{k=1}^{\ell}(y_{j},\psi_{k})_{X}\psi_{k}\right\|_{X}^{2}\\ &\text{subject to}\quad(\psi_{k},\psi_{t})=\delta_{kt}\quad\text{for any }k,t\in\{1,2,\cdots,\ell\}.\end{split}

(4)

Here, $\delta_{kt}$ denotes the Kronecker-delta function. A solution $\{\psi_{k}\}_{k=1}^{\ell}$ to (4) is called a POD-basis of rank $\ell$ . We introduce the correlation matrix

K=\left(\frac{1}{n}(y_{j},y_{i})_{X}\right)\in\mathbb{R}^{n\times n}

corresponding to the snapshots $\{y_{j}\}_{j=1}^{n}$ . The matrix $K$ is positive semi-definite and has rank $\mathcal{N}$ . The minimization problem (4) can be reduced to an eigenvalue problem

\displaystyle Kv=\lambda v.

(5)

We sort all the positive eigenvalues in a decreasing order as $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{\mathcal{N}}>0$ and the associated eigenvectors are denoted by $v_{k}$ with $k=1,\cdots,\mathcal{N}$ . It can be shown that the POD-basis of rank $\ell\in\mathbb{N}^{+}$ with $\ell\leq\mathcal{N}$ is formed by

\displaystyle\varphi_{k}=\frac{1}{\sqrt{\lambda_{k}}}\sum_{j=1}^{n}(v_{k})_{j}y_{j}\quad\text{for }k=1,\cdots,\ell.

(6)

Here, $(v_{k})_{j}$ is the $j$ -th component of the eigenvector $v_{k}$ . The basis functions $\{\varphi_{k}\}_{k=1}^{\ell}$ form a POD-basis of rank $\ell$ and we have the following error formula.

Proposition 2.2.

Let $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{\mathcal{N}}>0$ be the positive eigenvalues of $K$ in (5) and $v_{1},\cdots,v_{\mathcal{N}}\in\mathbb{R}^{n}$ be the associated eigenvectors. Then, $\{\varphi_{k}\}_{k=1}^{\ell}$ given by (6) forms a set of POD-basis of rank $\ell$ with $\ell\leq\mathcal{N}$ . Moreover, we have the error formula

\frac{1}{n}\sum_{j=1}^{n}\left\|y_{j}-\sum_{k=1}^{\ell}(y_{j},\varphi_{k})_{X}\varphi_{k}\right\|_{X}^{2}=\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}.

In practice, we shall make use of the decay property of eigenvalues in $\lambda_{k}$ and choose the first $\ell$ dominant eigenvalues such that the ratio $\zeta:=\frac{\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}}{\sum_{k=1}^{\mathcal{N}}\lambda_{k}}$ is small enough to achieve an expected accuracy, for instance $\zeta=1\%$ . One would prefer the eigenvalues decays as fast as possible so that one can ensure high accuracy with few POD basis functions.

3 Multiscale Method

In this section, we develop the computational multiscale method in order to solve the parabolic wave approximation. For spatial discretization, we will apply the CEM-GMsFEM. In particular, for each node $z_{i}\in D$ along the $z$ -direction, we will construct a set of multiscale basis functions in the spirit of CEM-GMsFEM. To further reduce the dimension of the multiscale space, we will perform POD procedure related to these CEM basis functions. Once the multiscale space is constructed, one can use leapfrog scheme to discretize time derivatives and solve the resulting fully-discretized problem.

3.1 Spatial Discretization: CEM-GMsFEM

In this section, we briefly outline the framework of CEM-GMsFEM and present the construction of the multiscale space. First, we introduce fine and coarse grids for the computational domain. Let $\mathcal{T}^{H}=\{K_{i}\}_{i=1}^{N}$ be a conforming partition of the domain $\Omega$ with mesh size $H>0$ defined by

H:=\max_{K\in\mathcal{T}^{H}}\Big{(}\max_{x,y\in K}\lvert x-y\rvert\Big{)}.

We refer to this partition as coarse grid. We denote $N\in\mathbb{N}^{+}$ the total number of coarse elements. Subordinate to the coarse grid, we define the fine grid partition $\mathcal{T}^{h}$ (with mesh size $h\ll H$ ) by refining each coarse element $K\in\mathcal{T}^{H}$ into a connected union of finer elements. We assume the refinement above is performed such that $\mathcal{T}^{h}$ is also a conforming partition of the domain $\Omega$ . Denote $N_{c}$ the number of interior coarse grid nodes of $\mathcal{T}^{H}$ and $\{x_{i}\}_{i=1}^{N_{c}}$ the collection of interior coarse nodes in the coarse grid.

3.1.1 Spectral Decomposition

We present the construction of the auxiliary multiscale basis functions. Let $K_{i}\in\mathcal{T}^{H}$ be a coarse block. Define $V(K_{i})$ as the restriction of the abstract space $V$ on the coarse element $K_{i}$ . For each $z_{k}\in D$ , we consider a local spectral problem: Find $\sigma_{j}^{(i,k)}\in\mathbb{R}$ and $\phi_{j}^{(i,k)}\in V(K_{i})$ such that

\displaystyle a_{i}(\phi_{j}^{(i,k)},v;z_{k})=\sigma_{j}^{(i,k)}s_{i}(\phi_{j}^{(i,k)},v;z_{k})\quad\text{for all}\leavevmode\nobreak\ v\in V(K_{i}).

(7)

Here, $a_{i}:V(K_{i})\times V(K_{i})$ is a symmetric non-negative definite bilinear form and $s_{i}:V(K_{i})\times V(K_{i})$ is a symmetric positive definite bilinear form. We remark that the above problem is solved on the fine mesh in the actual computations. Based on the analysis, we choose

a_{i}(v,w;z_{k}):=\int_{K_{i}}c(z_{k})\nabla v\cdot\nabla w\leavevmode\nobreak\ dx,\quad s_{i}(v,w;z_{k}):=\int_{K_{i}}\tilde{\kappa}(z_{k})vw\leavevmode\nobreak\ dx,

where $\tilde{\kappa}:=\sum_{j=1}^{N_{c}}c(z_{k})\lvert\nabla\chi_{j,k}^{\text{ms}}\rvert^{2}$ . The functions $\{\chi_{j,k}^{\text{ms}}\}_{j=1}^{N_{c}}$ are the standard multiscale finite element basis functions which satisfy the partition of unity property. More precisely, $\chi_{j,k}^{\text{ms}}$ is the solution of the following system:

$\displaystyle\nabla\cdot(c(z_{k})\nabla\chi_{j,k}^{\text{ms}})$	$\displaystyle=0\quad$	$\displaystyle\text{in each }K\subset\omega_{j},$
$\displaystyle\chi_{j,k}^{\text{ms}}$	$\displaystyle=g_{j}\quad$	$\displaystyle\text{on }\partial K\setminus\partial\omega_{j},$
$\displaystyle\chi_{j,k}^{\text{ms}}$	$\displaystyle=0\quad$	$\displaystyle\text{on }\partial\omega_{j}.$

The function $g_{j}$ is continuous and linear along the boundary of the coarse element. We assume that the eigenvalues $\sigma_{j}^{(i,k)}$ are arranged in ascending order and we pick $\ell_{i,k}\in\mathbb{N}^{+}$ corresponding eigenfunctions to construct the local auxiliary space $V_{\text{aux}}^{(i,k)}:=\text{span}\{\phi_{j}^{(i,k)}:j=1,\cdots,\ell_{i}\}$ . We assume the normalization $s_{i}\left(\phi_{j}^{(i,k)},\phi_{j}^{(i,k)};z_{k}\right)=1$ . After that, we define the global auxiliary multiscale space $V_{\text{aux}}^{k}:=\bigoplus_{i=1}^{N}V_{\text{aux}}^{(i,k)}$ . We remark that the global auxiliary multiscale space is used to construct multiscale basis functions that are orthogonal to the auxiliary space with respect to the weighted $L^{2}$ inner product $s(\cdot,\cdot;z_{k})$ .

Note that the bilinear form $s_{i}(\cdot,\cdot;z_{k})$ defines an inner product with norm $\left\|\cdot\right\|_{s(K_{i};k)}:=\sqrt{s_{i}(\cdot,\cdot;z_{k})}$ in the local auxiliary space $V_{\text{aux}}^{(i,k)}$ . Based on these local inner products and norms, one can naturally define a new inner product and norm for the global auxiliary space $V_{\text{aux}}$ as follows: for all $v,w\in V_{\text{aux}}^{k}$ ,

\displaystyle s(v,w;z_{k}):=\sum_{i=1}^{N}s_{i}(v,w;z_{k})\quad\text{and}\quad\left\|v\right\|_{s(z_{k})}:=\sqrt{s(v,v;z_{k})}.

(8)

Note that if $\{\chi_{j,k}^{\text{ms}}\}_{j=1}^{N_{c}}$ is a set of bilinear partition of unity, then $\left\|v\right\|_{s(z_{k})}\leq H^{-1}(\max\{c\})^{1/2}\left\|v\right\|$ for any $v\in L^{2}(\Omega)$ . In addition, we define $\pi_{k}:L^{2}(\Omega)\to V_{\text{aux}}^{k}$ as the projection with respect to the inner product $s(\cdot,\cdot;z_{k})$ such that

\pi_{k}(u):=\sum_{i=1}^{N}\sum_{j=1}^{\ell_{i}}s_{i}(u,\phi_{j}^{(i,k)};z_{k})\phi_{j}^{(i,k)}\quad\text{for all }u\in L^{2}(\Omega).

3.1.2 The construction of multiscale basis functions

In this section, we present the construction of the multiscale basis functions. First, we define an oversampling region for each coarse element. Specifically, given a non-negative integer $m\in\mathbb{N}$ and a coarse element $K_{i}$ , we define the oversampling region $K_{i,m}\subset\Omega$ such that

K_{i,m}:=\left\{\begin{array}[]{lr}K_{i}&\text{if }m=0,\\ \displaystyle{\bigcup\{K:K_{i,m-1}\cap K\neq\emptyset\}}&\text{if }m\geq 1.\end{array}\right.

See Figure 1 for an illustration of oversampling region. For simplicity, we denote $K_{i}^{+}$ the oversampled region $K_{i,m}$ for some nonnegative integer $m$ .

Refer to caption — Figure 1: Oversampling region with $m=1$ .

Recall that $V(K_{i}^{+})$ is the restriction of $V$ on the coarse patch $K_{i}^{+}$ . Let $V_{0}(K_{i}^{+})$ be the subspace of $V(K_{i}^{+})$ with zero trace on the boundary $\partial K_{i}^{+}$ . For each eigenfunction $\phi_{j}^{(i,k)}\in V_{\text{aux}}^{k}$ , we define the multiscale basis $\psi_{j,\text{ms}}^{(i,k)}\in V_{0}(K_{i}^{+})$ to be the solution of the equation:

\displaystyle a(\psi_{j,\text{ms}}^{(i,k)},v)+s\left(\pi(\psi_{j,\text{ms}}^{(i,k)}),\pi_{k}(v)\right)=s(\phi_{j}^{(i,k)},v)\quad\text{for all }v\in V_{0}(K_{i}^{+}).

(9)

Then, the multiscale space is defined as $V_{\text{ms}}^{k}:=\text{span}\left\{\psi_{j,\text{ms}}^{(i,k)}:i=1,\cdots,N,\leavevmode\nobreak\ j=1,\cdots,\ell_{i,k}\right\}$ . By construction, we have $\text{dim}(V_{\text{ms}}^{k})=\text{dim}(V_{\text{aux}}^{k})$ . After that, we define

V_{\text{CEM}}:=\text{span}\left\{\psi_{j,\text{ms}}^{(i,k)}:i=1,\cdots,N,\leavevmode\nobreak\ j=1,\cdots,\ell_{i,k},\leavevmode\nobreak\ k=0,\cdots,K\right\}.

We will perform the procedure of POD for this space to construct the multiscale space that will be used in actual simulation.

3.2 Construction of multiscale reduced basis functions using POD

In this section, we present the construction of multiscale reduced basis functions using the POD technique. First, we define $\{v_{\text{CEM},k}\}_{k=0}^{K}\subset V_{\text{CEM}}$ to be the solution of the following equation:

\displaystyle\left(\ddot{v}_{\text{CEM},k},w;z_{k}\right)_{c}+\frac{1}{\Delta z}\left(\dot{v}_{\text{CEM},k}-\dot{v}_{\text{CEM},k-1},w\right)+\frac{1}{2}a(v_{\text{CEM},k},w;z_{k})=0\quad\text{ for all }w\in V_{\text{CEM}}.

(10)

In the abstract framework of POD, we set the space $X=L^{2}(0,T;H^{1}(\Omega))$ . We define the inner product associated with $X$ as $\left\langle u,v\right\rangle:=\int_{0}^{T}(\nabla u,\nabla v)+(u,v)\leavevmode\nobreak\ dt$ for $u,v\in X$ . The corresponding norm will then be defined as $\left\|\cdot\right\|_{L^{2}(0,T;H^{1}(\Omega))}:=\sqrt{\left\langle\cdot,\cdot\right\rangle}$ . Next, we define the snapshot space $\mathcal{V}$ as the collection of:

y_{j}:=v_{\text{CEM}}(t,z_{j-1}),\quad y_{j+K+1}:=\ddot{v}_{\text{CEM}}(t,z_{j-1}),\quad j=1,\cdots,K+1,

and

y_{j+2K+2}=\tilde{\partial}\dot{v}_{\text{CEM}}(t,z_{j}):=\frac{\dot{v}_{\text{CEM}}(t,z_{j})-\dot{v}_{\text{CEM}}(t,z_{j-1})}{\Delta z},\quad j=1,\cdots,K.

Then, we perform the POD procedure on the snapshot space $\mathcal{V}$ as described in Section 2.1. It should be noted that the correlation matrix is defined as $K_{ij}:=\frac{1}{3K+2}(y_{j},y_{i})_{H^{1}}$ . We denote $\{\psi_{k}\}_{k=1}^{\mathcal{N}}$ the corresponding POD basis functions with cardinality $\mathcal{N}\in\mathbb{N}^{+}$ . We then define $V_{\text{POD}}^{\ell}:=\operatorname{span}\{\psi_{k}\}_{k=1}^{\ell}$ for a given positive number $\ell\leq\mathcal{N}$ .

3.3 Fully Discretization

In this section, we present the fully discretization for the problem (3). We can further consider the discretization in time. We divide $(0,T]$ into $N$ pieces and write

v^{n}=v(t^{n})\quad\text{with}\quad t^{n}=n\Delta t\quad\text{and}\quad\Delta t=\frac{T}{N}.

Specifically, we use first- and second-order central difference schemes. We define $\tau:=\Delta t\Delta z$ . The fully discretization reads: for $k=1,\cdots,K-1$ and $n=1,\cdots,N-1$ , find $\{v_{k}^{n}\}\subset V_{\text{POD}}^{\ell}$ such that

\displaystyle\left(\frac{v_{k}^{n+1}-2v_{k}^{n}+v_{k}^{n-1}}{(\Delta t)^{2}},w;z_{k}\right)_{c}+\left(\frac{v_{k}^{n+1}-v_{k}^{n-1}}{2\tau},w\right)+\frac{1}{2}a(v_{k}^{n},w;z_{k})=\left(\frac{v_{k-1}^{n+1}-v_{k-1}^{n-1}}{2\tau},w\right)

(11)

for all $w\in V_{\text{POD}}^{\ell}$ . In short, we discretize implicitly in $z$ and explicitly in time. We show that the fully discretization is stable in time. To this aim, we first recall the inverse inequality for the multiscale space, which is proved in [11].

Proposition 3.1 (Lemma 4.5 in [11]).

Assume that $\{\chi_{j,k}^{\text{ms}}\}_{j=1}^{N_{c}}$ is a set of bilinear partition of unity. For any $v\in V_{\text{CEM}}$ , there is a constant $C_{\text{inv}}>0$ such that

\displaystyle\left\|\nabla v\right\|_{a}\leq C_{\text{inv}}H^{-1}\left\|v\right\|.

Then, we have the following stability result for (11).

Lemma 3.2.

Suppose that the following CFL condition

\frac{1}{c_{\max}}-\frac{1}{4}C_{\text{inv}}^{2}H^{-2}(\Delta t)^{2}\geq\delta,

holds for some $\delta>0$ . Then, we have the following stability estimate:

\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{K}^{n+1}-v_{K}^{n-1}\right\|^{2}+\sum_{k=1}^{K}\mathcal{E}_{N,k}\leq\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{0}^{n+1}-v_{0}^{n-1}\right\|^{2}+\sum_{k=1}^{K}\mathcal{E}_{0,k}.

Proof.

Define the energy as follows:

\displaystyle\mathcal{E}_{n,k}:=\left\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\|_{c}^{2}+\frac{1}{2}a(v_{k}^{n+1},v_{k}^{n}).

Using the inverse inequality for the multiscale functions, one can show that

	$\displaystyle\mathcal{E}_{n,k}$	$\displaystyle=\left\\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\\|_{c}^{2}+\frac{1}{2}a(v_{k}^{n+1},v_{k}^{n})$
		$\displaystyle=\left\\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\\|_{c}^{2}+\frac{1}{4}a(v_{k}^{n+1},v_{k}^{n+1})+\frac{1}{4}a(v_{k}^{n},v_{k}^{n})-\frac{1}{4}a(v_{k}^{n+1}-v_{k}^{n},v_{k}^{n+1}-v_{k}^{n})$
		$\displaystyle\geq\left\\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\\|_{c}^{2}+\frac{1}{4}\left(\left\\|v_{k}^{n+1}\right\\|_{a}^{2}+\left\\|v_{k}^{n}\right\\|_{a}^{2}\right)-\frac{1}{4}C_{\text{inv}}^{2}H^{-2}(\Delta t)^{2}\left\\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\\|^{2}$
		$\displaystyle\geq\big{(}\frac{1}{c_{\max}}-\frac{1}{4}C_{\text{inv}}^{2}H^{-2}(\Delta t)^{2}\big{)}\left\\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\\|^{2}+\frac{1}{4}\left(\left\\|v_{k}^{n+1}\right\\|_{a}^{2}+\left\\|v_{k}^{n}\right\\|_{a}^{2}\right)\geq 0.$

Next, taking $w=v_{k}^{n+1}-v_{k}^{n-1}$ in (11) and denoting $f^{n}_{k-1}=v_{k-1}^{n+1}-v_{k-1}^{n-1}$ , we have

\displaystyle\frac{1}{2\tau}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}+\left\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\|_{c}^{2}-\left\|\frac{v_{k}^{n}-v_{k}^{n-1}}{\Delta t}\right\|_{c}^{2}+\frac{1}{2}a(v_{k}^{n},v_{k}^{n+1}-v_{k}^{n-1})=\frac{1}{2\tau}\left(f^{n}_{k-1},v_{k}^{n+1}-v_{k}^{n-1}\right).

Then, it implies that

\displaystyle\begin{split}\frac{1}{2\tau}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}+\mathcal{E}_{n,k}-\mathcal{E}_{n-1,k}&=\frac{1}{2\tau}\left(f^{n}_{k},v_{k}^{n+1}-v_{k}^{n-1}\right)\\ &\leq\frac{1}{2\tau}\left\|v_{k-1}^{n+1}-v_{k-1}^{n-1}\right\|\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|\\ &\leq\frac{1}{4\tau}\left(\left\|v_{k-1}^{n+1}-v_{k-1}^{n-1}\right\|^{2}+\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}\right)\\ \implies\frac{1}{4\tau}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}+\mathcal{E}_{n,k}&\leq\frac{1}{4\tau}\left\|v_{k-1}^{n+1}-v_{k-1}^{n-1}\right\|^{2}+\mathcal{E}_{n-1,k}.\end{split}

Summing over $n=1,\cdots,N$ , we have

\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}+\mathcal{E}_{N,k}\leq\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{k-1}^{n+1}-v_{k-1}^{n-1}\right\|^{2}+\mathcal{E}_{0,k}.

Therefore, summing over $k=1,\cdots,K$ , we have

\frac{1}{4\tau}\sum_{n=1}^{N}\left(\sum_{k=1}^{K}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}-\sum_{k=0}^{K-1}\left\|v_{k}^{n+1}-v_{k}^{n-1}\right\|^{2}\right)+\sum_{k=1}^{K}\mathcal{E}_{N,k}\leq\sum_{k=1}^{K}\mathcal{E}_{0,k}

and

\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{K}^{n+1}-v_{K}^{n-1}\right\|^{2}+\sum_{k=1}^{K}\mathcal{E}_{N,k}\leq\frac{1}{4\tau}\sum_{n=1}^{N}\left\|v_{0}^{n+1}-v_{0}^{n-1}\right\|^{2}+\sum_{k=1}^{K}\mathcal{E}_{0,k}.

This completes the proof. ∎

Remark.

It should be noted that the CFL condition depends on the quantity $c_{\max}:=\left\|c\right\|_{L^{\infty}(D\times\Omega)}$ . It also relies on the constant $C_{\text{inv}}$ which comes from the inverse inequality (see Lemma 4.5 in [11]). From the proof of that lemma, the constant $C_{\text{inv}}$ is independent of the ratio of the contrast values in the permeability. We are using the CEM-GMsFEM, which is a coarse mesh method, hence the coarse mesh size $H$ will make the CFL much less restrictive. The idea to improve the CFL condition is to use the implicit scheme in time; more precisely, one can use $a(v_{k}^{n+1},w;z_{k})$ in (11). This will become one of our future works.

4 Convergence Analysis

In this section, we present the convergence analysis of the proposed POD-based multiscale method for the parabolic wave equation. Let $v_{\text{CEM}}\in V_{\text{CEM}}$ be the semi-discretized solution which solves the problem (10). Throughout this section, we assume that the coercivity and boundedness of this bilinear form hold. Specifically, there exist two constants $\gamma$ and $\beta$ , independent of any $z_{k}\in D$ , such that

\displaystyle\gamma\|u\|_{H_{1}(\Omega)}^{2}\leq a(u,u;z_{k})\quad\text{and}\quad a(u,w;z_{k})\leq\beta\|u\|_{H_{1}(\Omega)}\|w\|_{H_{1}(\Omega)}

(12)

for any $u,w\in V$ . We denote $(\cdot,\cdot)_{H^{1}}$ the $H^{1}$ -inner product such that $(u,v)_{H^{1}}:=(\nabla u,\nabla v)+(u,v)$ for any $u,v\in V$ with its associated norm being $\left\|\cdot\right\|_{H^{1}}:=\sqrt{(\cdot,\cdot)_{H^{1}}}$ .

We have the following lemma for the proper orthogonal decomposition.

Lemma 4.1 (cf. Proposition 1 in [39]).

Let $v_{\text{CEM}}\in V_{\text{CEM}}$ be the semi-discretized solution which solves the problem (10) and $K_{ij}:=\frac{1}{3K+2}(y_{j},y_{i})_{H^{1}}$ be the correlation matrix and $\lambda_{k}$ be the corresponding eigenvalues sorted ascending. For any integer $\ell$ with $0<\ell\leq\mathcal{N}$ , the following error formula holds:

\displaystyle\begin{split}\frac{1}{3K+2}&\left[\int_{0}^{T}\sum_{k=1}^{K+1}\left\|v_{\text{CEM}}(t,z_{k})-\sum_{j=1}^{\ell}\left(v_{\text{CEM}}(t,z_{k}),\psi_{j}\right)_{H^{1}}\psi_{j}\right\|_{H^{1}}^{2}\right.\\ &+\left.\sum_{k=1}^{K}\left\|\tilde{\partial}\dot{v}_{\text{CEM}}(t,z_{k})-\sum_{j=1}^{\ell}\left(\tilde{\partial}\dot{v}_{\text{CEM}}(t,z_{k}),\psi_{j}\right)_{H^{1}}\psi_{j}\right\|_{H^{1}}^{2}\right.\\ &\left.+\sum_{k=1}^{K+1}\left\|\ddot{v}_{\text{CEM}}(t,z_{k})-\sum_{j=1}^{\ell}\left(\ddot{v}_{\text{CEM}}(t,z_{k}),\psi_{j}\right)_{H^{1}}\psi_{j}\right\|_{H^{1}}^{2}dt\right]=\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}.\end{split}

(13)

Next, we define the Ritz-projection $P^{\ell}:V\rightarrow V_{\text{POD}}^{\ell}$ by:

\displaystyle(P^{\ell}v,\psi)_{H^{1}}=(v,\psi)_{H^{1}}\quad\text{ for all }\psi\in V_{\text{POD}}^{\ell},

(14)

for any $v\in V$ . We have the following estimate for the projection operator.

Lemma 4.2.

Let $\{z_{k}\}_{k=1}^{K+1}\subset D$ . For any $\ell\in\{1,2,...,\mathcal{N}\}$ and $v\in V$ , the projection operator $P^{\ell}$ satisfies the following estimate:

\displaystyle\frac{1}{K+1}\sum_{k=1}^{K+1}\int_{0}^{T}\left\|v_{\text{CEM}}(z_{k})-P^{\ell}v_{\text{CEM}}(z_{k})\right\|_{H^{1}}^{2}\leavevmode\nobreak\ dt\lesssim\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}.

(15)

Proof.

Let $v\in V$ be arbitrary. By the definition of the projection $P^{\ell}$ , we have

\displaystyle\begin{split}\left\|v-P^{\ell}v\right\|_{H^{1}}^{2}=(v-P^{\ell}v,v-P^{\ell}v)_{H^{1}}=(v-P^{l}v,v-\psi)_{H^{1}}\leq\left\|v-P^{\ell}v\right\|_{H^{1}}\left\|v-\psi\right\|_{H^{1}}\end{split}

for all $\psi\in V_{\text{POD}}^{\ell}$ . This implies that,

\displaystyle\left\|v-P^{\ell}v\right\|_{H^{1}}\lesssim\left\|v-\psi\right\|_{H^{1}}\quad\text{ for all }\psi\in V_{\text{POD}}^{\ell}.

(16)

The lemma follows immediately by Lemma 4.1 and (16). ∎

Clearly, from the proof of Lemma 4.2, the corollary below follows immediately.

Corollary 4.3.

Let $v_{\text{CEM}}\in V_{\text{CEM}}$ be the semi-discretized solution which solves (10).

Then, the following estimates hold:

\displaystyle\frac{1}{K}\sum_{k=1}^{K}\int_{0}^{T}\left\|\tilde{\partial}\dot{v}_{\text{CEM}}(z_{k})-P^{\ell}\tilde{\partial}\dot{v}_{\text{CEM}}(z_{k})\right\|_{H^{1}}^{2}\leavevmode\nobreak\ dt\lesssim\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}

(17)

and

\displaystyle\frac{1}{K+1}\sum_{k=0}^{K}\int_{0}^{T}\left\|\ddot{v}_{\text{CEM}}(z_{k})-P^{\ell}\ddot{v}_{\text{CEM}}(z_{k})\right\|_{H^{1}}^{2}dt\lesssim\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}.

(18)

Let $\{U_{k}\}_{k=0}^{K}\subset V_{\text{POD}}^{\ell}$ be the semi-discretized solution satisfying the following equation

\displaystyle(\ddot{U}_{k},\phi;z_{k})_{c}+(\tilde{\partial}\dot{U}_{k},\phi)+\frac{1}{2}a(U_{k},\phi;z_{k})=0

(19)

for all $\phi\in V_{\text{POD}}^{\ell}$ with appropriate initial condition $U_{0}$ . Here, we denote $\tilde{\partial}\dot{U}_{k}:=\frac{\dot{U}_{k}-\dot{U}_{k-1}}{\Delta z}$ . We are now able to analyze the error. Assume that $\{v_{k}^{n}\}\subset V_{\text{POD}}^{\ell}$ solves (11), $\{U_{k}\}_{k=0}^{K}\subset V_{\text{POD}}^{\ell}$ solves (19), and the function $v_{\text{CEM}}\in V_{\text{CEM}}$ solves (10). We decompose the error into three parts:

v_{k}^{n}-v_{\text{CEM}}(t_{n})=\underbrace{v_{k}^{n}-U_{k}(t_{n})}_{=:\mu_{k}^{n}}+\underbrace{U_{k}(t_{n})-P^{\ell}v_{\text{CEM}}(t_{n})}_{=:\nu_{k}^{n}}+\underbrace{P^{\ell}v_{\text{CEM}}(t_{n})-v_{\text{CEM}}(t_{n})}_{=:\rho_{k}^{n}}.

We denote $\mu_{k}(t)$ , $\nu_{k}(t)$ , and $\rho_{k}(t)$ the piecewise linear functions that interpolates $\{\mu_{k}^{n}\}$ , $\{\nu_{k}^{n}\}$ , and $\{\rho_{k}^{n}\}$ in time, respectively. Due to Lemma 4.2, we have

\displaystyle\frac{1}{K}\sum_{k=1}^{K}\left\|\rho_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}\lesssim\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}.

It should be noted that the term $\sum_{k=\ell+1}^{\mathcal{N}}\lambda_{k}$ comes from the method of POD and is a typical error term in the POD analysis. The error decay to $0$ very fast; and by the theory, this error is optimal [39] since the method of POD solves a minimization problem (4).

Next, using the notation $\tilde{\partial}\nu_{k}=\frac{\nu_{k}-\nu_{k-1}}{\Delta z}$ for all $k=1,...,K$ and the equation (10), we obtain

	$\displaystyle(\ddot{\nu}_{k},\psi;z_{k})_{c}+(\tilde{\partial}\dot{\nu}_{k},\psi)+\frac{1}{2}a(\nu_{k},\psi;z_{k})$	$\displaystyle=-(P^{\ell}\ddot{v}_{\text{CEM}},\psi;z_{k})_{c}-(\tilde{\partial}P^{\ell}\dot{v}_{\text{CEM}},\psi)-\frac{1}{2}a(v_{\text{CEM}},\psi;z_{k})$
		$\displaystyle=(\ddot{v}_{\text{CEM}}-P^{\ell}\ddot{v}_{\text{CEM}},\psi;z_{k})_{c}+(\tilde{\partial}\dot{v}_{\text{CEM}}-\tilde{\partial}P^{\ell}\dot{v}_{\text{CEM}},\psi)$

for any $\psi\in V_{\text{POD}}^{\ell}$ . Let us denote

h_{k}:=\tilde{\partial}\dot{v}_{\text{CEM}}-\tilde{\partial}P^{\ell}\dot{v}_{\text{CEM}}\quad\text{and}\quad w_{k}:=\ddot{v}_{\text{CEM}}-P^{\ell}\ddot{v}_{\text{CEM}}.

By Corollary 4.3, it follows that

\displaystyle\frac{1}{K}\sum_{k=1}^{K}\left\|h_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}+\frac{1}{K+1}\sum_{k=0}^{K}\left\|w_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}\lesssim\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}.

Using the same technique of showing Lemma 2.1, one can obtain an estimate for $\nu_{k}$ :

\displaystyle\|\dot{\nu}_{k}\|_{L^{2}(0,T;L^{2}(\Omega))}^{2}+\Delta z\sum_{k=1}^{K}\mathcal{E}_{k}(\nu_{k}(T))\lesssim\Delta z\sum_{k=1}^{K}\left\|h_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}+\left\|w_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}\lesssim\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}.

It remains to estimate the term $\mu_{k}$ . This error comes from the temporal discretization and it follows that

\left\|\mu_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}=\int_{0}^{T}\left\|\mu_{k}\right\|_{H^{1}}^{2}dt=\int_{0}^{T}O(\tau^{2})dt=O(\tau).

Denote $v_{k}(t)$ the piecewise linear function that interpolates $\{v_{k}^{n}\}$ in time. As a result, we have the following error estimate

\displaystyle\frac{1}{K}\sum_{k=1}^{K}\left\|v_{k}-v_{\text{CEM}}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}\lesssim O(\tau)+\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}.

(20)

We remark that the error between the solution $\theta$ (that solves (3)) and the solution $v_{\text{CEM}}$ has the relation: $\left\|\theta-v_{\text{CEM}}\right\|_{L^{2}(0,T;H^{1}(\Omega))}=O(H)$ . Therefore, the whole error estimate reads as follows:

\displaystyle\frac{1}{K}\sum_{k=1}^{K}\left\|\theta-v_{k}\right\|_{L^{2}(0,T;H^{1}(\Omega))}^{2}\lesssim O(H^{2})+O(\tau)+\sum_{i=\ell+1}^{\mathcal{N}}\lambda_{i}.

(21)

5 Numerical Experiments

In this section, we present some numerical results to demonstrate the efficiency of the proposed computational multiscale method. In the experiments below, we set the spatial domain to be $\Omega=(0,1)^{2}$ . We partition the spatial domain into $10\times 10$ uniform square elements and we refer this partition to the coarse grid with mesh size $H=\sqrt{2}/10$ . Further, we divide each coarse element into $10\times 10$ uniform square elements and the corresponding fine grid has resolution of the size $100\times 100$ . We refer this to the fine grid of the spatial domain with mesh size $h=\sqrt{2}/100$ . We compute the numerical solutions $v_{k}^{N}$ at the terminal time for all $k=1,\cdots,K$ and compare it with the reference solutions $v_{f}(z_{k},x,T)$ for all $k$ , which are computed by using the underlying fine grid. We measure the relative error in $L^{2}$ norm at each layer $z_{k}$ , which is defined as follows:

e_{2}^{k}:=\frac{\left\|v_{k}^{N}-v_{f}(z_{k},x,T)\right\|}{\left\|v_{f}(z_{k},x,T)\right\|}.

5.1 The First Experiment

In the first experiment, we set the initial condition to be $v(0,x,t)=\sin(\pi x_{1})\sin(\pi x_{2})\sin(t)$ where $x=(x_{1},x_{2})\in\Omega$ ; also we set homogeneous initial condition for all $z$ , i.e., $v(z,x,0)=0$ for all $(z,x)\in D\times\Omega$ . The temporal direction and $z$ -direction are discretized as in the scheme (11). We choose the time step to be $\Delta t=10^{-5}$ and the step size in $z$ -direction is $\Delta z=10^{-4}$ . We set $K=30$ , i.e., the $z$ direction is partitioned into $30$ levels and $D=[0,30\Delta z]$ .

In this experiment, the pattern of the permeability field $c(x)$ is given in Figure 2 (Left). In $x_{1}$ - and $x_{2}$ -directions, the permeability will follow the same pattern. However, along the $z$ -direction, the permeability field $c(x)$ has different values of contrast. See Figure 2 (Right) for an illustration for this layer-structured heterogeneous media. In particular, we partition the permeability field along the $z$ -direction into $3$ different layers; the value of contrast in the yellow region within each layer are equal to $10$ , $15$ , and $20$ , respectively. We remark that the patterns keep the same as demonstrated in Figure 2 (Left).

To obtain the set of POD basis functions, we first construct the CEM basis functions for each layer and then obtain the reduced basis functions by performing the POD procedure on the combined basis space. We choose 100 POD basis functions from the total $300$ CEM basis functions in this experiment. We remark that the number of the POD reduced basis functions depends on the decay of the eigenvalues $\{\lambda_{k}\}$ corresponding to the POD construction. In practice, one has to include sufficiently many reduced basis functions so that the tail of the eigenvalues is smaller than a given tolerance of accuracy.

The record of relative error $e_{2}^{k}$ for all $k\in\{1,\cdots,K\}$ at the terminal time is plotted in Figure 3 (Left). The relative error at each layer is around the level of magnitude of $10^{-3}$ . One can observe from the graph that the proposed numerical scheme is capable for accurately approximating the fine-scale solution.

5.2 Experiment 2

In the second experiment, we consider the so-called Marmousi heterogeneous field. More precisely, given $z_{k}\in D$ , $c(z_{k})$ is a Marmousi permeability field as shown in Figure 4.

The initial and boundary data are the same as in the first experiment. We set $T=5\times 10^{-3}$ and $\Delta t=5\times 10^{-6}$ . We set $\Delta z=10^{-4}$ and $D=[0,10\Delta z]$ with $K=10$ . For each $z_{k}\in D$ , the medium $c(z_{k})$ is obtained by taking the data from a very fine Marmousi field of size $600\times 600$ . To be more specific, we partition the fine field (which is of size $600\times 600$ ) into a coarse field of size $100\times 100$ whose coarse element size is $6\times 6$ . Then, in each coarse element, we randomly pick a value in one of the fine elements in this coarse element; and we formulate our computational fine grid ( $100\times 100$ ) as combining all selected elements accordingly.

To get the set of POD basis functions, we first calculate three sets of CEM basis functions with $k\in\{0,4,8\}$ and obtain the POD basis by performing the POD procedure on the combined set of CEM basis functions. We choose $50$ POD basis functions from the total $300$ CEM basis functions. The record of relative error for all $z_{k}$ at the terminal time is plotted in Figure 3 (Right). The relative error at each layer is around the level of magnitude of $10^{-4}$ in this case. In particular, the relative error $e_{2}^{k}$ is about $3.9\times 10^{-4}$ for $k=10$ . This demonstrates the efficiency of the proposed multiscale method for the simulation of parabolic wave formulation.

6 Conclusion

In this paper, we developed a computational multiscale method for simulating the parabolic wave formulation with highly heterogeneous media. For the spatial discretization, we employed the recently developed CEM-GMsFEM, which is proved to be efficient to reduce the dimension of the model in space. We then combined the technique of proper orthogonal decomposition to further reduce the dimension along the quasi time direction. A complete analysis of the proposed algorithm has been provided. Numerical results are provided to demonstrate the effectiveness and efficiency of the proposed multiscale method.

Acknowledgement

References

[1] A. Abdulle and B. Engquist. Finite element heterogeneous multiscale methods with near optimal computational complexity. Multiscale Modeling & Simulation, 6(4):1059–1084, 2007.
[2] A. Bamberger, B. Engquist, L. Halpern, and P. Joly. Higher order paraxial wave equation approximations in heterogeneous media. SIAM Journal on Applied Mathematics, 48(1):129–154, 1988.
[3] A. Bamberger, B. Engquist, L. Halpern, and P. Joly. Parabolic wave equation approximations in heterogenous media. SIAM Journal on Applied Mathematics, 48(1):99–128, 1988.
[4] A. Bamberger, B. Euggnist, L. Helpern, and P. Joly. The paraxial approximation for the wave equation: some new results. Adv. Comput. Methods Partial Differential Equations, pages 340–344, 1984.
[5] P. Blomgren, G. Papanicolaou, and H. Zhao. Super-resolution in time-reversal acoustics. The Journal of the Acoustical Society of America, 111(1):230–248, 2002.
[6] H. Brock, R. Buchal, and C. Spofford. Modifying the sound-speed profile to improve the accuracy of the parabolic-equation technique. The Journal of the Acoustical Society of America, 62(3):543–552, 1977.
[7] D. L. Brown, Y. Efendiev, and V. H. Hoang. An efficient hierarchical multiscale finite element method for stokes equations in slowly varying media. Multiscale Modeling & Simulation, 11(1):30–58, 2013.
[8] E. Cances, V. Ehrlacher, F. Legoll, and B. Stamm. An embedded corrector problem to approximate the homogenized coefficients of an elliptic equation. Comptes Rendus Mathematique, 353(9):801–806, 2015.
[9] J. Chen, S. Sun, and Z. He. Homogenize coupled stokes–cahn–hilliard system to darcy’s law for two-phase fluid flow in porous medium by volume averaging. Journal of Porous Media, 22(1), 2019.
[10] J. Chen, S. Sun, and X. Wang. Homogenization of two-phase fluid flow in porous media via volume averaging. Journal of Computational and Applied Mathematics, 353:265–282, 2019.
[11] B. Chetverushkin, E. Chung, Y. Efendiev, S.-M. Pun, and Z. Zhang. Computational multiscale methods for quasi-gas dynamic equations. arXiv preprint arXiv:2009.00068, 2020.
[12] B. Chetverushkin, A. Saveliev, and V. Saveliev. Compact quasi-gasdynamic system for high-performance computations. Computational Mathematics and Mathematical Physics, 59(3):493–500, 2019.
[13] B. N. Chetverushkin, N. D’Ascenzo, A. V. Saveliev, and V. Saveliev. Kinetic model and magnetogasdynamics equations. Computational Mathematics and Mathematical Physics, 58(5):691–699, 2018.
[14] E. Chung, Y. Efendiev, and T. Y. Hou. Adaptive multiscale model reduction with generalized multiscale finite element methods. Journal of Computational Physics, 320:69–95, 2016.
[15] E. Chung, Y. Efendiev, and C. S. Lee. Mixed generalized multiscale finite element methods and applications. Multiscale Modeling & Simulation, 13:338–366, 2014.
[16] E. Chung, Y. Efendiev, and W. T. Leung. Constraint energy minimizing generalized multiscale finite element method. Computer Methods in Applied Mechanics and Engineering, 339:298–319, 2018.
[17] E. Chung, Y. Efendiev, and W. T. Leung. Constraint energy minimizing generalized multiscale finite element method in the mixed formulation. Computational Geosciences, 22(3):677–693, 2018.
[18] E. T. Chung, Y. Efendiev, and W. T. Leung. Generalized multiscale finite element methods for wave propagation in heterogeneous media. Multiscale Modeling & Simulation, 12(4):1691–1721, 2014.
[19] E. T. Chung, Y. Efendiev, and W. T. Leung. Fast online generalized multiscale finite element method using constraint energy minimization. Journal of Computational Physics, 355:450–463, 2018.
[20] E. T. Chung, Y. Efendiev, W. T. Leung, M. Vasilyeva, and Y. Wang. Non-local multi-continua upscaling for flows in heterogeneous fractured media. Journal of Computational Physics, 372:22–34, 2018.
[21] E. T. Chung, W. T. Leung, and S. Pollock. Goal-oriented adaptivity for GMsFEM. Journal of Computational and Applied Mathematics, pages 625–637, 2015.
[22] J. F. Claerbout. Fundamentals of geophysical data processing, volume 274. Citeseer, 1976.
[23] J. D. Cole. Modern developments in transonic flow. SIAM Journal on Applied Mathematics, 29(4):763–787, 1975.
[24] W. E and B. Engquist. The Heterogeneous Multiscale Methods. Communications in Mathematical Sciences, 1(1):87–132, 2003.
[25] Y. Efendiev, J. Galvis, and T. Y. Hou. Generalized multiscale finite element methods (GMsFEM). Journal of Computational Physics, 251:116–135, 2013.
[26] D. Fafalis and J. Fish. Computational continua for linear elastic heterogeneous solids on unstructured finite element meshes. International Journal for Numerical Methods in Engineering, 115(4):501–530, 2018.
[27] J. Fish and S. Kuznetsov. Computational continua. International Journal for Numerical Methods in Engineering, 84(7):774–802, 2010.
[28] J. Fish and Z. Yuan. Multiscale enrichment based on partition of unity. International Journal for Numerical Methods in Engineering, 62(10):1341–1359, 2005.
[29] S. Fu, E. Chung, and G. Li. Edge multiscale methods for elliptic problems with heterogeneous coefficients. Journal of Computational Physics, 2019.
[30] K. Gao, S. Fu, R. L. Gibson Jr, E. T. Chung, and Y. Efendiev. Generalized multiscale finite-element method (GMsFEM) for elastic wave propagation in heterogeneous, anisotropic media. Journal of Computational Physics, 295:161–188, 2015.
[31] H. Hajibeygi, D. Kavounis, and P. Jenny. A hierarchical fracture model for the iterative multiscale finite volume method. Journal of Computational Physics, 230(4):8729–8743, 2011.
[32] A. Hasegawa and F. Tappert. Transmission of stationary nonlinear optical pulses in dispersive dielectric fibers. i. anomalous dispersion. Applied Physics Letters, 23(3):142–144, 1973.
[33] P. Henning, A. Målqvist, and D. Peterseim. A localized orthogonal decomposition method for semi-linear elliptic problems. ESAIM: Mathematical Modelling and Numerical Analysis, 48(5):1331–1349, 2014.
[34] V. Hoang and C. Schwab. High dimensional finite elements for elliptic problems with multiple scales. Multiscale Modeling & Simulation, 3:168–194, 2004.
[35] T. Y. Hou and X.-H. Wu. A multiscale finite element method for elliptic problems in composite materials and porous media. Journal of Computational Physics, 134:169–189, 1997.
[36] J. Hudson. A parabolic approximation for elastic waves. Wave Motion, 2(3):207–214, 1980.
[37] P. Jenny, S. Lee, and H. Tchelepi. Multi-scale finite volume method for elliptic problems in subsurface flow simulation. Journal of Computational Physics, 187:47–67, 2003.
[38] P. Jenny, S. Lee, and H. Tchelepi. Adaptive multi-scale finite volume method for multi-phase flow and transport in porous media. Multiscale Modeling & Simulation, 3:30–64, 2004.
[39] K. Kunisch and S. Volkwein. Galerkin proper orthogonal decomposition methods for parabolic problems. Numerische mathematik, 90(1):117–148, 2001.
[40] C. Le Bris, F. Legoll, and A. Lozinski. An MsFEM type approach for perforated domains. Multiscale Modeling & Simulation, 12(3):1046–1077, 2014.
[41] C. Le Bris, F. Legoll, and F. Thomines. Multiscale finite element approach for weakly random problems and related issues. ESAIM: Mathematical Modelling and Numerical Analysis, 48(3):815–858, 2014.
[42] A. Målqvist and D. Peterseim. Localization of elliptic multiscale problems. Mathematics of Computation, 83(290):2583–2603, 2014.
[43] S. T. McDaniel. Parabolic approximations for underwater sound propagation. The Journal of the Acoustical Society of America, 58(6):1178–1185, 1975.
[44] S. T. McDaniel. Propagation of normal mode in the parabolic approximation. The Journal of the Acoustical Society of America, 57(2):307–311, 1975.
[45] H. Owhadi and L. Zhang. Metric-based upscaling. Communications on Pure and Applied Mathematics, 60:675–723, 2007.
[46] A. Roberts and I. Kevrekidis. General tooth boundary conditions for equation free modeling. SIAM Journal on Scientific Computing, 29(4):1495–1510, 2007.
[47] A. Salama, S. Sun, M. F. El Amin, Y. Wang, and K. Kumar. Flow and Transport in Porous Media: A Multiscale Focus. Geofluids, 2017, 2017.
[48] G. Samaey, I. Kevrekidis, and D. Roose. Patch dynamics with buffers for homogenization problems. Journal of Computational Physics, 213(1):264–287, 2006.
[49] G. Samaey, D. Roose, and I. Kevrekidis. The gap-tooth scheme for homogenization problems. Multiscale Modeling & Simulation, 4(1):278–306, 2005.
[50] W. C. Tan and V. H. Hoang. High dimensional finite element method for multiscale nonlinear monotone parabolic equations. Journal of Computational and Applied Mathematics, 345:471–500, 2019.

	$\displaystyle\mathcal{E}_{n,k}$	$\displaystyle=\left\\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\\|_{c}^{2}+\frac{1}{2}a(v_{k}^{n+1},v_{k}^{n})$
		$\displaystyle=\left\\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\\|_{c}^{2}+\frac{1}{4}a(v_{k}^{n+1},v_{k}^{n+1})+\frac{1}{4}a(v_{k}^{n},v_{k}^{n})-\frac{1}{4}a(v_{k}^{n+1}-v_{k}^{n},v_{k}^{n+1}-v_{k}^{n})$
		$\displaystyle\geq\left\\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\\|_{c}^{2}+\frac{1}{4}\left(\left\\|v_{k}^{n+1}\right\\|_{a}^{2}+\left\\|v_{k}^{n}\right\\|_{a}^{2}\right)-\frac{1}{4}C_{\text{inv}}^{2}H^{-2}(\Delta t)^{2}\left\\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\\|^{2}$
		$\displaystyle\geq\big{(}\frac{1}{c_{\max}}-\frac{1}{4}C_{\text{inv}}^{2}H^{-2}(\Delta t)^{2}\big{)}\left\\|\frac{v_{k}^{n+1}-v_{k}^{n}}{\Delta t}\right\\|^{2}+\frac{1}{4}\left(\left\\|v_{k}^{n+1}\right\\|_{a}^{2}+\left\\|v_{k}^{n}\right\\|_{a}^{2}\right)\geq 0.$