An adaptive dynamically low-dimensional approximation method for multiscale stochastic diffusion equations

TL;DR

This paper introduces an adaptive low-dimensional approximation method combining multiscale basis functions and generalized polynomial chaos to efficiently solve complex multiscale stochastic diffusion equations.

Contribution

It develops a novel adaptive DyBO method with multiscale basis functions and online adaptation for improved efficiency and accuracy in multiscale stochastic PDEs.

Findings

Significant reduction in computational cost demonstrated.

High accuracy maintained with multiscale and online basis functions.

Method effectively handles multiscale and stochastic features.

Abstract

In this paper, we propose a dynamically low-dimensional approximation method to solve a class of time-dependent multiscale stochastic diffusion equations. A dynamically bi-orthogonal (DyBO) method was developed to explore low-dimensional structures of stochastic partial differential equations (SPDEs) and solve them efficiently. However, when the SPDEs have multiscale features in physical space, the original DyBO method becomes expensive. To address this issue, we construct multiscale basis functions within each coarse grid block for dimension reduction in the physical space. To further improve the accuracy, we also perform online procedure to construct online adaptive basis functions. In the stochastic space, we use the generalized polynomial chaos (gPC) basis functions to represent the stochastic part of the solutions. Numerical results are presented to demonstrate the efficiency of the proposed method in solving time-dependent PDEs with multiscale and random features.

Figures

Equations

$$\frac{\partial u^{\epsilon }}{\partial t}(x,t,\omega )={\mathcal{L}}^{\epsilon }{u}^{\epsilon }(x,t,\omega ),x\in \mathcal{D},t\in (0,T],\omega \in \Omega ,$$

$$u^{\epsilon }(x,t,\omega )\approx {\bar{u}}^{\epsilon }(x,t)+\sum _{i=1}^m{u}_i^{\epsilon }(x,t)Y_i(\omega ,t)={\bar{u}}^{\epsilon }(x,t)+\mathbf{U}(x,t){\mathbf{Y}}^T(\omega ,t),$$

$$\frac{\partial u^{\epsilon }}{\partial t}$$

$$u^{\epsilon }(x,0,\omega )$$

$$\displaystyle\mathcal{B}\big{(}u^{\varepsilon}(x,t,\omega)\big{)}$$

$${\tilde{u}}^{\epsilon }(x,t,\omega )={\bar{u}}^{\epsilon }(x,t)+\sum _{i=1}^m{u}_i^{\epsilon }(x,t)Y_i(\omega ,t)={\bar{u}}^{\epsilon }(x,t)+\mathbf{U}(x,t){\mathbf{Y}}^T(\omega ,t)\approx u^{\epsilon }(x,t,\omega ),$$

$$\mathbf{U}(x,t)=\left(u_1^{\epsilon }(x,t),\cdots ,u_m^{\epsilon }(x,t)\right)\text{and}\mathbf{Y}(\omega ,t)=\left(Y_1(\omega ,t),\cdots ,Y_m(\omega ,t)\right)$$

$$⟨{\mathbf{U}}^T,\mathbf{U}⟩=(⟨u_i,u_j⟩{\delta }_{ij})\text{and}\mathbb{E}\left[{\mathbf{Y}}^T\mathbf{Y}\right]={\mathbf{I}}_{m\times m},$$

$$\frac{\partial \bar{u}}{\partial t}$$

$$\frac{\partial \mathbf{U}}{\partial t}$$

$$\frac{\mathrm{d}\mathbf{Y}}{\mathrm{d}t}$$

$$\mathcal{Q}(\mathbf{M}):=\frac{1}{2}\left(\mathbf{M}-{\mathbf{M}}^T\right)\text{and}\tilde{\mathcal{Q}}(\mathbf{M}):=\mathcal{Q}(\mathbf{M})+\mathrm{diag}(\mathbf{M}),$$

$$\mathbf{C}-{\mathbf{\Lambda }}_{\mathbf{U}}^{-1}\tilde{\mathcal{Q}}\left({\mathbf{\Lambda }}_{\mathbf{U}}\mathbf{C}\right)$$

$$\mathbf{D}-\mathcal{Q}\left(\mathbf{D}\right)$$

$${\mathbf{D}}^T+\mathbf{C}$$

$${\int }_{-\infty }^{\infty }{H}_i(\xi )H_j(\xi )\rho (\xi )d\xi ={\delta }_{ij}.$$

$$Y_i(\omega ,t)\approx \sum _{\mathbf{\alpha }\in {\mathfrak{J}}_r^p}{\mathbf{H}}_{\mathbf{\alpha }}(\mathbf{\xi }(\omega ))A_{\mathbf{\alpha }i}(t)=\mathbf{H}\left(\mathbf{\xi }\right){\mathbf{A}}_i(t)i=1,2,\cdots ,m.$$

$$\mathbf{Y}(\omega ,t)=\mathbf{H}\left(\mathbf{\xi }(\omega )\right)\mathbf{A}(t)$$

$$\tilde{u}\approx \bar{u}+\mathbf{U}{\mathbf{A}}^T{\mathbf{H}}^T.$$

$$\frac{\partial \bar{u}}{\partial t}$$

$$\frac{\partial \mathbf{U}}{\partial t}$$

$$\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}t}$$

$$G_{\ast }(\bar{u},\mathbf{U},\mathbf{Y})={\mathbf{\Lambda }}_{\mathbf{U}}^{-1}⟨{\mathbf{U}}^T,\mathbb{E}\left[\tilde{\mathcal{L}}\tilde{u}\mathbf{Y}\right]⟩={\mathbf{\Lambda }}_{\mathbf{U}}^{-1}⟨{\mathbf{U}}^T,\mathbb{E}\left[\tilde{\mathcal{L}}\tilde{u}\mathbf{H}\right]⟩\mathbf{A}.$$

$$u^{\text{DyBO-gPC}}=\bar{u}+\mathbf{U}{\mathbf{A}}^T{\mathbf{H}}^T.$$

$$\tilde{a}(x,\omega )=\sum _{i=1}^r{a}_i(x){\xi }_i(\omega ),i=1,2\cdots ,r.$$

$$\frac{\partial \bar{u}}{\partial t}$$

$$\frac{\partial \mathbf{U}}{\partial t}$$

$$\frac{\mathrm{d}\mathbf{A}}{\mathrm{d}t}$$

$$\frac{\partial w}{\partial t}$$

$$w{|}_{t=0}$$

$$\frac{w^n-w^{n-1}}{\Delta t}=\nabla \cdot (\bar{a}\nabla w^n)+\mathcal{G}$$

$$-\nabla \cdot (\bar{a}\nabla w^n)+c{w}^n=\tilde{\mathcal{G}},$$

$$-\nabla \cdot (\bar{a}\nabla {\psi }_j^{(i)})$$

$${\psi }_j^{(i)}$$

$${\int }_{D_i}\bar{a}\nabla {\varphi }_j^{(i)}\cdot \nabla v={\lambda }_j^{(i)}{\int }_{D_i}\hat{a}{\varphi }_j^{(i)}v\forall v\in V_{snap}^{(i)},j=1,\cdots ,L_i,$$

$$-\nabla \cdot (\bar{a}\nabla {\chi }_i)$$

Full text

1 Introduction

Uncertainty arises in many real-world problems of scientific applications, such as heat propagation through random media or flow driven by stochastic forces. These kind of problems usually have multiple scale features involved in the spatial domain. For example, to simulate flows in heterogeneous porous media, the permeability field is often parameterized by random fields with multiple-scale structures.

Stochastic partial differential equations (SPDEs), which contain random variables or stochastic processes, play important roles in modeling complex problems and quantifying the corresponding uncertainties. Considerable amounts of efforts have been devoted to study SPDEs, see [babuska:04, Ghanem:91, WuanHou:06, Zabaras:09, Knio:01, matthies:05, Najm:09, Webster:08, Wan:06, Xiu:09] and references therein. These methods are effective when the dimension of solution space is not huge. However, when SPDEs have multiscale features, the SPDE problems become difficult since it requires tremendous computational resources to resolve the small scales of the SPDE solutions. This motivates us to develop efficient numerical schemes to solve these challenging problems.

In this paper, we shall consider the time-dependent SPDEs with multiscale coefficients as follows

[TABLE]

where suitable boundary and initial conditions are imposed, $\mathcal{D}\subset\mathbb{R}^{d}$ is a bounded spatial domain, $\Omega$ is a sample space, and $\mathcal{L}^{\varepsilon}$ is an elliptic operator that contains multiscale and random coefficient, where the smallest-scale is parameterized by $\varepsilon$.

The major difficulties in solving (1) come from two parts. In the physical space, we need a mesh fine enough to resolve the small-scale features. In the random space, we need extra degrees of freedom to represent the random features. Moreover, the problem (1) becomes more difficult if the dimension of the random input is high. To address these challenges, we shall explore low-dimensional structures hidden in the solution $u^{\varepsilon}(x,t,\omega)$. Specifically, if the solution $u^{\varepsilon}(x,t,\omega)$ is a second-order stochastic process at each time $t>0$, i.e., $u^{\varepsilon}(x,t,\omega)\in L^{2}(\mathcal{D}\times\Omega)$ for each $t$, one can approximate the solution $u^{\varepsilon}(x,t,\omega)$ by its $m$-term truncated Karhunen-Loève (KL) expansion [Karhunen1947, Loeve1978]

[TABLE]

where $\mathbf{U}(x,t)=(u_{1}^{\varepsilon}(x,t),\cdots,u_{m}^{\varepsilon}(x,t))$ and $\mathbf{Y}(\omega,t)=(Y_{1}(\omega,t),\cdots,Y_{m}(\omega,t))$. The KL expansion gives the compact representation of the solution. However, the direct computation of the KL expansion can be quite expensive since we need to form a covariance kernel and solve a large-scale eigenvalue problem.

In [ChengHouZhang1:13, ChengHouZhang2:13], a dynamically bi-orthogonal (DyBO) method was developed. This new method derives an equivalent system that can faithfully track the KL expansion of the SPDE solution. In other words, the DyBO method gives the evolution equations for $\bar{u}^{\varepsilon}$, $\mathbf{U}$, and $\mathbf{Y}$. The DyBO method can accurately and efficiently solve many time-dependent SPDEs, such as stochastic Burger’s equations and stochastic Navier-Stokes equations, with considerable savings over existing numerical methods. To explore the low-dimensional features of the solutions to the SPDEs, a dynamically orthogonal (DO) method was proposed [sapsis:09]. Later on, the equivalence of DO method and DyBO method has been shown in [choi2014equivalence] and the effectiveness of the DO and DyBO has also been discussed theoretically in [musharbash2015error]. This area is very active and highly demanded due to the latest advances in the UQ research.

If the SPDEs have multiscale features in the physical space (i.e., the smallest-scale parameter $\varepsilon$ is extremely small), however, the original DyBO method (as well as the DO method) becomes computationally expensive since one needs enormous degrees of freedom to represent the multiscale features in the physical space. To overcome this difficulty, we shall apply the generalized multiscale finite element methods (GMsFEM) [chung2016adaptive, efendiev2013generalized] to construct multiscale basis functions within each coarse grid block for model reduction in the physical space.

In the GMsFEM, we divide the computation into two stages: the offline stage and the online one. In the offline stage, we first compute global snapshot functions within each coarse neighborhood based on the given coarse and fine meshes and construct multiscale basis functions to represent the local heterogeneities. When the snapshot functions are computed, one can construct the multiscale basis functions in each coarse patch by solving some well-designed local spectral problems and identify the crucial multiscale basis functions to form the offline function space. In the online stage, we add more online multiscale basis functions that are constructed using the offline space. These online basis functions are computed adaptively in some selected spatial regions based on the current local residuals and their construction is motivated by the analysis in [chung2015online]. In general, the algorithm can guarantee that additional online multiscale basis functions will reduce the error rapidly if one chooses a sufficient number of offline basis functions. We should point out that there are many existing methods in the literature to solve multiscale problems though; see [Engquist2012, Bourgeat1984, Cruz1995, Dykaar1992, Engquist2003, HanZhangCMS:12, Hou1997, Juanes2005, Peterseim2014, Nolen2008] and references therein. Most of these methods are designed for multiscale problems with deterministic coefficients.

In our new method, we first derive the DyBO formulation for the multiscale SPDEs (1), which consists of deterministic PDEs for $\bar{u}^{\varepsilon}$ and $\mathbf{U}$ respectively and an ODE system for the stochastic basis $\mathbf{Y}$. For the deterministic PDEs (for $\bar{u}^{\varepsilon}$ and $\mathbf{U}$) in the formulation, we shall apply the GMsFEM to construct multiscale basis functions and use these multiscale basis functions to represent $\bar{u}^{\varepsilon}$ and $\mathbf{U}$. It leads to considerable savings over the original DyBO method. For the ODE system, the memory cost is relatively small and we shall apply a suitable ODE solver to compute the numerical solution. The GMsFEM enables us to significantly improve the efficiency of the DyBO method in solving time-dependent PDEs with multiscale and random coefficients.

The rest of the paper is organized as follows. In Section 2, we will introduce the framework of DyBO formulation. The GMsFEM and its online adaptive algorithm will be outlined in Section 3. The implementation issues of the algorithm and the numerical results will be given in Section 4. Finally, some concluding remarks will be drawn in Section LABEL:sec:conclusion.

2 The DyBO formulation for multiscale time-dependent SPDEs

In this paper, we consider a class of parabolic equations with multiscale and random coefficients

[TABLE]

where $\mathcal{D}\subset\mathbb{R}^{d}$ ($d=2,3$) is a bounded spatial domain, $(\Omega,\mathcal{F},\mathbb{P})$ is a probability space, and suitable boundary and initial conditions are imposed. The differential operator $\mathcal{L}^{\varepsilon}$ is defined as $\mathcal{L}^{\varepsilon}u^{\varepsilon}:=\nabla\cdot(a^{\varepsilon}(x,\omega))\nabla u^{\varepsilon}+f(x)$. The multiscale information is described by the parameter $\varepsilon$ and the force $f:\mathbb{R}^{d}\to\mathbb{R}$ is in $L^{2}(\mathcal{D})$. Assume that there exist two constants $a_{\max}\gg a_{min}>0$ such that $\mathbb{P}(\omega\in\Omega:a^{\varepsilon}(x,\omega)\in[a_{\min},a_{\max}],\text{a.e. }x\in\mathcal{D})=1$. Note that we are interested in the case that the coefficient $a^{\varepsilon}(x,\omega)$ has high contrast within the domain $\mathcal{D}$, where the mode reduction technique in physical space is necessary to reduce degrees of freedom in representing the solution.

2.1 An abstract framework for SPDEs

To make this paper self-contained, we briefly review the DyBO method [ChengHouZhang1:13, ChengHouZhang2:13]. We assume the solution $u^{\varepsilon}(x,t,\omega)$ to (3) satisfies $u^{\varepsilon}(\cdot,t,\cdot)\in L^{2}(\mathcal{D}\times\Omega)$ for each time $t\in(0,T]$. We consider the $m$-term truncated KL expansion with $m\in\mathbb{N}^{+}$

[TABLE]

as an approximation to the solution $u^{\varepsilon}(x,t,\omega)$. Here, $\bar{u}^{\varepsilon}(x,t)$ is the mean of the solution,

[TABLE]

are the spatial and stochastic modes (with zero mean), respectively. We omit the symbol $\varepsilon$ to simplify the notation. Next, by imposing the orthogonal conditions for $\mathbf{U}$ and $\mathbf{Y}$

[TABLE]

we obtain the evolution equations for $\bar{u}^{\varepsilon}$, $\mathbf{U}$ and $\mathbf{Y}$ as follows

[TABLE]

where $\left\langle\cdot,\,\cdot\right\rangle$ denotes the inner product in $L^{2}(\mathcal{D})$, $\mathbf{\Lambda}_{\mathbf{U}}=\mathrm{diag}(\left\langle\mathbf{U}^{T},\,\mathbf{U}\right\rangle)\in\mathbb{R}^{m\times m}$, and $\tilde{\mathcal{L}}\widetilde{u}=\mathcal{L}\widetilde{u}-\mathbb{E}\left[\mathcal{L}\widetilde{u}\right]$. Define two operators $\mathcal{Q}:\mathbb{R}^{k\times k}\rightarrow\mathbb{R}^{k\times k}$ and $\tilde{\mathcal{Q}}:\mathbb{R}^{k\times k}\rightarrow\mathbb{R}^{k\times k}$ as follows

[TABLE]

where $\mathbf{M}\in\mathbb{R}^{k\times k}$ is a square matrix and $\mathrm{diag}(\mathbf{M})$ is a diagonal matrix whose diagonal entries are equal to those of $\mathbf{M}$. Then, the matrices $\mathbf{C},\mathbf{D}\in\mathbb{R}^{m\times m}$ in (5) can be solved uniquely from the following linear system

[TABLE]

where the matrix is given by $G_{*}(\bar{u}^{\varepsilon},\mathbf{U},\mathbf{Y})=\mathbf{\Lambda}_{\mathbf{U}}^{-1}\left\langle\mathbf{U}^{T},\,\mathbb{E}\left[\tilde{\mathcal{L}}\widetilde{u}^{\varepsilon}\mathbf{Y}\right]\right\rangle\in\mathbb{R}^{m\times m}$.

In order to represent the stochastic modes $Y_{i}(\omega,t)$ in (5c), one can choose several different approaches including ensemble representations in sampling methods and spectral representations. In this paper, we use gPC basis functions to represent the stochastic modes $Y_{i}(\omega,t)$. Given two positive integers $r$ and $p$, we define $\mathfrak{J}_{r}^{p}:=\{\boldsymbol{\alpha}:\boldsymbol{\alpha}=(\alpha_{1},\cdots,\alpha_{r}),\alpha_{i}\in\mathbb{N}_{+},\lvert\alpha\rvert=\sum_{i=1}^{r}\alpha_{i}\leq p\}\backslash\{0\}$. Let $\{H_{i}(\xi)\}_{i=1}^{\infty}$ denote as an one-dimensional family of $\rho$-orthogonal polynomial, i.e.,

[TABLE]

If we write $\mathbf{H}_{\boldsymbol{\alpha}}(\boldsymbol{\xi})=\prod_{i=1}^{r}H_{\alpha_{i}}(\xi_{i})$ for $\boldsymbol{\alpha}\in\mathfrak{J}_{r}^{p}$ and $\boldsymbol{\xi}=(\xi_{i})_{i=1}^{r}$, then the Cameron-Martin theorem [Cameron1947] implies the stochastic modes $Y_{i}(\omega,t)$ in (4) can be approximated by

[TABLE]

Here, $\mathbf{H}\left(\boldsymbol{\xi}\right)=\left(\mathbf{H}_{\boldsymbol{\alpha}}\left(\boldsymbol{\xi}\right)\right)_{\boldsymbol{\alpha}\in\mathfrak{J}_{r}^{p}}\in\mathbb{R}^{1\times N_{p}}$, $\mathbf{A}_{i}(t)=(\mathbf{A}_{\boldsymbol{\alpha}i}(t))_{\boldsymbol{\alpha}\in\mathfrak{J}_{r}^{p}}\in\mathbb{R}^{N_{p}\times 1}$ and $N_{p}:=\lvert\mathfrak{J}_{r}^{p}\rvert$. We remark that for each $i=1,2,\cdots,m$, the coefficients $\{\mathbf{A}_{\boldsymbol{\alpha}i}(t)\}_{\boldsymbol{\alpha}\in\mathfrak{J}_{r}^{p}}$ represent the projection coefficients of the stochastic mode $Y_{i}(\omega,t)$ on the gPC basis functions $\mathbf{H}_{\boldsymbol{\alpha}}(\boldsymbol{\xi})$, $\boldsymbol{\alpha}\in\mathfrak{J}_{r}^{p}$. Moreover, $\{\mathbf{A}_{\boldsymbol{\alpha}i}(t)\}_{\boldsymbol{\alpha}\in\mathfrak{J}_{r}^{p}}$ change with respect to time. One may write

[TABLE]

where $\mathbf{A}(t)=(\mathbf{A}_{1}(t),\cdots,\mathbf{A}_{m}(t))\in\mathbb{R}^{N_{p}\times m}$. The KL expansion (4) now reads

[TABLE]

We can derive equations for $\bar{u}$, $\mathbf{U}$ and $\mathbf{A}$, instead of $\bar{u}$, $\mathbf{U}$ and $\mathbf{Y}$. Here and in the following, we have suppressed the variables $x$, $t$, and $\omega$ for notation simplicity. In other words, the stochastic modes $\mathbf{Y}$ are identified with a matrix $\mathbf{A}\in\mathbb{R}^{N_{p}\times m}$, which leads to the DyBO-gPC formulation of SPDE (3)

[TABLE]

where $\mathbf{C}(t)$ and $\mathbf{D}(t)$ can be solved from the linear system (6) with

[TABLE]

By solving the system (9), we have an approximate solution to (3)

[TABLE]

The condition $\mathbb{E}\left[\mathbf{Y}^{T}\mathbf{Y}\right]$ implies that the columns $(\mathbf{A}_{i})_{i=1}^{m}$ are orthonormal, i.e., $\mathbf{A}^{T}\mathbf{A}=\mathbf{I}_{m\times m}$. Note that $\mathbf{A}\mathbf{A}^{T}\in\mathbb{R}^{N_{p}\times N_{p}}$ in general is not an identity matrix as $m\ll N_{p}$ if the SPDE solution has a low-dimensional structure.

2.2 The DyBO formulation for the model problem

In this section, we shall derive the DyBO formulation for the model problem (3). Recall that the definition of the differential operator is $\mathcal{L}u=\nabla\cdot(a(x,\omega)\nabla u)+f(x)$ and we have omitted $\varepsilon$ for notation simplification. We assume that the coefficient $a(x,\omega)$ is of the form $a(x,\omega)=\bar{a}(x)+\tilde{a}(x,\omega)$, where $\bar{a}(x)=\mathbb{E}\left[a(x,\omega)\right]$ and $\tilde{a}(x,\omega)$ is the fluctuation, which can be parametrized as follows

[TABLE]

Here, $r\geq 1$ is a positive integer and $\{\xi_{i}(\omega)\}_{i=1}^{r}$ are independent identically distributed random variables assumed to be mean-zero. By substituting the expression of $\mathcal{L}u$ into (9), we obtain the DyBO-gPC formulation for the model problem (3) (see LABEL:cha:dybo-formulation-MsSPDE for the details of the derivation)

[TABLE]

where matrices $\mathbf{C}$ and $\mathbf{D}$ can be solved from (6) with $G_{*}=\mathbf{\Lambda}_{\mathbf{U}}^{-1}\left\langle\mathbf{U}^{T},\,\mathbb{E}\left[\tilde{\mathcal{L}}u\mathbf{H}\right]\right\rangle\mathbf{A}$.

Remark 2.1*.*

When the force function of the model problem (3) contains randomness, i.e., $f(x,\omega)$, one can derive the DyBO formulation accordingly without any difficulty.

Remark 2.2*.*

The boundary conditions and initial conditions for each physical component, and the initial condition for each stochastic component can be obtained by projection of the initial and boundary conditions of $u(x,t,\omega)$ on the corresponding components.

Remark 2.3*.*

As the system evolves, the norm of the mode $u_{i}$ (denoted as $\lambda_{i}$) in the KL expansion may change and some of them may get closer to each other. In this case, if the matrices $\mathbf{C}$ and $\mathbf{D}$ are still solved from (6), numerical errors will pollute the results. One may freeze $\mathbf{U}$ or $\mathbf{Y}$ temporarily for a short time and continue to evolve the system. At the end of this short period, the solution is recast into the bi-orthogonal form via the KL expansion. See [ChengHouZhang1:13, Section 4.2] for more details.

3 Multiscale model reduction using the GMsFEM

3.1 Motivations

Since the DyBO formulation (11) involves multiscale features in the physical space, one may consider an efficient solver to solve the problem in order to reduce the computational cost. As such, we shall apply the GMsFEM to discretize $\bar{u}$ and $\mathbf{U}$. Note that, Eq.(11a) and each component of Eq. (11b) have the following deterministic time-dependent PDE form

[TABLE]

In order to discretize the equation (12) in time, we apply the implicit Euler scheme with time step $\Delta t>0$ and obtain the discretization for each time $t_{n}=n\Delta t$, $n=1,2,\cdots,N$ ($T=N\Delta t$)

[TABLE]

where $w^{n}=w(t_{n})$ and the above equation is equivalent to the following

[TABLE]

where $c=1/\Delta t$ and $\tilde{\mathcal{G}}=cw^{n-1}+\mathcal{G}$. Hence, for each fixed $t_{n}>0$, we use the GMsFEM to solve the second order elliptic PDE (13) with multiscale coefficient $\bar{a}$.

3.2 The GMsFEM and the multiscale basis functions

Next, we present the framework of the GMsFEM for solving (13). We first introduce the notion of fine and coarse grids that we shall use in the method. Let $\mathcal{T}^{H}$ be a conforming partition of the spatial domain $\mathcal{D}$ with mesh size $H>0$. We refer to this partition as the coarse grid. Subordinate to $\mathcal{T}^{H}$, we define a fine grid partition denoted by $\mathcal{T}^{h}$, with mesh size $0<h\ll H$, by refining each coarse element in $\mathcal{T}^{H}$ into a connected union of fine elements. Assume the above refinement is performed such that $\mathcal{T}^{h}$ is a conforming partition of $\mathcal{D}$. Denote the interior nodes of $\mathcal{T}^{H}$ as $x_{i}$, $i=1,2,\cdots,N_{in}$, where $N_{in}$ is the number of interior nodes. The coarse elements of $\mathcal{T}^{H}$ are denoted as $K_{j}$, $j=1,2.,\cdots,N_{e}$, where $N_{e}$ is the number of the coarse elements. Define the coarse neighborhood of the node $x_{i}$ by $D_{i}:=\bigcup\{K_{j}\in\mathcal{T}^{H}:x_{i}\in\overline{K_{j}}\}$.

Once the coarse and fine grids are given, one may construct the multiscale basis functions for approximating the solution of (13). To obtain the multiscale basis functions, we first define the snapshot space. For each neighborhood $D_{i}$, define $J_{h}(D_{i})$ as the set of fine nodes of $\mathcal{T}^{h}$ lying on $\partial D_{i}$ and denote its cardinality as $L_{i}\in\mathbb{N}^{+}$. For each fine-grid node $x_{j}\in J_{h}(D_{i})$, define a fine-grid function $\delta_{j}^{h}$ on $J_{h}(D_{i})$ as $\delta_{j}^{h}(x_{k})=\delta_{jk}$. Next, for $j=1,\cdots,L_{i}$, define the snapshot function $\psi_{j}^{(i)}$ in coarse neighborhood $D_{i}$ as the solution to the following system

[TABLE]

The local snapshot space $V_{snap}^{(i)}$ corresponding to the coarse neighborhood $D_{i}$ is defined as follows $V_{snap}^{(i)}:=\text{snap}\{\psi_{j}^{(i)}:j=1,\cdots,L_{i}\}$ and the snapshot space reads $V_{snap}:=\bigoplus_{i=1}^{N_{in}}V_{snap}^{(i)}$.

The snapshot space defined above is usually of large dimension. Therefore, a dimension reduction is performed on $V_{snap}$ to archive a smaller space $V_{\text{off}}$, which contains the multiscale basis functions for simulation. This reduction is achieved by performing a spectral decomposition on each local snapshot space $V_{snap}^{(i)}$. The analysis in [Yalchin2011] motivates the following construction. For each $i=1,\cdots,N_{in}$, the spectral problem is to find $(\phi_{j}^{(i)},\lambda_{j}^{(i)})\in V_{snap}^{(i)}\times\mathbb{R}$ such that

[TABLE]

where $\hat{a}:=\bar{a}\sum_{i=1}^{N_{in}}H^{2}\lvert\nabla\chi_{i}\rvert^{2}$ and $\{\chi_{i}\}_{i=1}^{N_{in}}$ is a set of partition of unity satisfying the following system

[TABLE]

Assume that the eigenvalues obtained from (16) are arranged in ascending order and we may use the first $0<l_{i}\leq L_{i}$ (with $l_{i}\in\mathbb{N}^{+}$) eigenfunctions (related to the smallest $l_{i}$ eigenvalues) to form the local multiscale space $V_{\text{off}}^{(i)}:=\text{span}\{\chi_{i}\phi_{j}^{(i)}:j=1,\cdots,l_{i}\}$. The multiscale space $V_{\text{off}}$ is the direct sum of the local multiscale spaces, namely $V_{\text{off}}:=\bigoplus_{i=1}^{N_{in}}V_{\text{off}}^{(i)}$.

Once the multiscale space $V_{\text{off}}$ is constructed, we can find the GMsFEM solution $u_{\text{off}}^{n}\in V_{\text{off}}$ at time $t=t_{n}$, $n=1,\cdots,N$, by solving the following equation

[TABLE]

where $\mathcal{A}(u,v):=\int_{\mathcal{D}}\bar{a}\nabla u\cdot\nabla v$.

Remark 3.1*.*

The above derivation of $V_{\text{off}}$ is based on the (mean) coefficient $\bar{a}$ and the multiscale basis functions in $V_{\text{off}}$ are suitable for approximating $\bar{u}$. In this paper, we assume that the fluctuation of the coefficient is a small perturbation to the mean. Therefore, the multiscale space $V_{\text{off}}$ can also efficiently approximate $\mathbf{U}$.

3.3 Online adaptive algorithm

In order to achieve a rapid convergence in the GMsFEM, one may add some online basis functions to enrich the multiscale space $V_{\text{off}}$ based on local residuals. In this subsection, we briefly outline the online adaptive algorithm for GMsFEM.

Let $u_{\text{off}}^{n}\in V_{\text{off}}$ be the numerical solution obtained in (17) at time $t=t_{n}$. Given a coarse neighborhood $D_{i}$, we define $V_{i}:=H_{0}^{1}(D_{i})\cap V_{snap}$ equipped with the norm $\|v\|_{V_{i}}^{2}:=\int_{D_{i}}\bar{a}(x)\lvert\nabla v\rvert^{2}$. We also define the local residual operator $\mathcal{R}_{i}^{n}:V_{i}\to\mathbb{R}$ by

[TABLE]

where $u_{\text{f}}^{n-1}$ is the fine-scale solution at time $t=t_{n-1}$. The operator norm of $\mathcal{R}_{i}^{n}$, denoted by $\|\mathcal{R}_{i}^{n}\|_{V_{i}^{*}}$ gives a measure of the quantity of the residual. The online basis functions are computed during the time-marching process for a given fixed time $t=t_{n}$, contrary to offline basis functions that are pre-computed (defined in the Section 4.2).

Suppose that one needs to add an online basis $\phi$ into the space $V_{i}$. The analysis in [chung2015online] suggests that the required online basis $\phi\in V_{i}$ is the solution to the following equation

[TABLE]

We refer to $\tau\in\mathbb{N}$ as the level of enrichment and denote $u_{\text{off}}^{n,\tau}$ as the solution of (17) in $V_{\text{off}}^{n,\tau}$. Remark that $V_{\text{off}}^{n,0}:=V_{\text{off}}$ for time level $n\in\mathbb{N}$. Let $\mathcal{I}\subset\{1,2,\cdots,N_{in}\}$ be the index set over some non-overlapping coarse neighborhoods. For each $i\in\mathcal{I}$, we obtain an online basis $\phi_{i}\in V_{i}$ by solving (19) and define $V_{\text{off}}^{n,\tau+1}=V_{\text{off}}^{n,\tau}\oplus\text{span}\{\phi_{i}:i\in\mathcal{I}\}$. After that, solve (17) in $V_{\text{off}}^{n,\tau+1}$ to get $u_{\text{off}}^{n,\tau+1}$. Consequently, following the arguments in [chung2015online], we have at time $t=t_{n}$,

[TABLE]

where $C_{\text{err}}$ is a uniform constant and $\Lambda_{\text{min}}^{(\mathcal{I})}=\min_{i\in\mathcal{I}}\lambda_{l_{i}+1}^{(i)}$. Here, the norm is defined by $\|\cdot\|_{V}:=\sqrt{\mathcal{A}(\cdot,\cdot)}$. Inequality (20) shows that we can obtain a better accuracy by adding more online basis functions at each time $t=t_{n}$ and the rate of convergence depends on the constant $C_{\text{err}}$ and $\Lambda_{\text{min}}^{(\mathcal{I})}$.

3.4 The implementation of our new algorithm

We summarize the computational scheme for the problem in this section. Recall that the multiscale coefficient is $a(x,\omega)=\bar{a}(x)+\tilde{a}(x,\omega)$. The mean has high contrast in nature and the fluctuation part is $\tilde{a}(x,\omega)=\sum_{i=1}^{r}a_{i}\xi_{i}=a_{i}\xi_{i}$, where the Einstein notation is used. We rewrite the DyBO formulation for (9) as follows

[TABLE]

We assume that homogeneous boundary condition is imposed. Hence, the solutions $\bar{u}$ and $\mathbf{U}=(u_{1},\cdots,u_{m})$ will vanish on $\partial\mathcal{D}$. If the model problem (3) has inhomogeneous boundary condition, the boundary conditions for $\bar{u}$ and $\mathbf{U}$ can be obtained by taking expectation of $u$ on the corresponding components. The initial conditions for $\bar{u}$, $\mathbf{U}$ and $\mathbf{A}$ depend on the initial condition of $u$, which will be discussed in Section 4. For the details of implementation, see LABEL:app:implement.

4 Numerical experiments

In this section, we present some numerical examples to demonstrate the efficiency of our proposed method. The computational domain is $\mathcal{D}=(0,1)^{2}\subset\mathbb{R}^{2}$ and $T=1$. First, we divide the domain $\mathcal{D}$ into several equal square units with mesh size $H>0$ and refer to it as the coarse mesh $\mathcal{T}^{H}$. Next, we divide each coarse element into several equal square blocks with mesh size $h>0$ and refer to it as the fine mesh $\mathcal{T}^{h}$. Then, we discretize $\bar{u}$ and $\mathbf{U}$ in the DyBO formulation by using the GMsFEM to reduce degrees of freedom in representing the multiscale solutions. Thus, with the multiscale basis functions, we can represent $\bar{u}$ and $\mathbf{U}$ on the coarse mesh. In all examples, the number of initial local basis functions is $L_{i}=4$.

The multiscale coefficient is assumed to be $a(x,\omega)=\bar{a}(x)+\sum_{i=1}^{r}a_{i}(x)\xi_{i}(\omega)$, where $a_{i}(x)$ is a small (or multiscale) perturbation and $\{\xi_{i}(\omega)\}_{i=1}^{r}$ is a set of i.i.d. uniform-distributed random variables over $[-1,1]$. Moreover, we assume that there exist two constants $a_{\max}\gg a_{\min}>0$ such that $\mathbb{P}(\omega\in\Omega:a(x,\omega)\in[a_{\min},a_{\max}],~{}\text{a.e. }x\in\mathcal{D})=1$.

The initial condition of the solution is assumed to have the form of $m$-terms truncated KL expansion

[TABLE]

The stochastic basis $Y_{i}(\omega,t)$ can be expanded as $Y_{i}(\omega,t)=\sum_{j=1}^{N_{p}}H_{j}(\omega)A_{ji}(t)$ for each $i=1,\cdots,m$. Here, $\{H_{j}(\omega)\}_{j=1}^{N_{p}}$ is a set of tensor products of orthogonal polynomials in $\mathbb{R}$, $N_{p}=\frac{(p+r)!}{p!r!}-1$, and $p$ is the maximum degree of polynomial. Denote $\mathbf{A}(t)=(A_{ji}(t))_{N_{p}\times m}$. The initial condition of the matrix \mathbf{A}(t)|_{t=0}=\big{(}A_{ji}(0)\big{)}_{N_{p}\times m} should satisfy $\mathbb{E}\left[\mathbf{H}\mathbf{A}\right]=0$ and $\mathbf{A}^{T}(t)\mathbf{A}(t)|_{t=0}=\mathbf{I}_{m\times m}$.

For each function to be approximated (e.g. $\bar{u}$, $u_{i}$ or the variance function $\text{var}(u):=\sum_{i=1}^{m}u_{i}^{2}$), we define the following quantity at $t=t_{n}$ to measure the numerical error

[TABLE]

where $u_{\text{f}}^{n}$ is the reference solution and $u_{\text{approx}}^{n}$ is the approximation obtained by the proposed method. In the remaining part of this paper, we refer to this quantity as $L^{2}$-error.

Example 4.1**.**

We set the mesh size to be $H=\sqrt{2}/10$ and $h=\sqrt{2}/100$. The time step is $\Delta t=10^{-3}$. The multiscale fluctuation is parameterized by three independent random variables ($r=3$) and the number of terms in the KL expansion is $m=4$. Next, we set the coefficients $a_{i}$ ($i=1,2,3$) to be

[TABLE]

The mean $\bar{a}$ of the multiscale coefficient is of high-contrast (see Figure LABEL:fig:exp3_abar). The source function is chosen to be $f\equiv 1$. The initial conditions for the mean of the solution and the physical modes are given as follows

[TABLE]