Prediction of Discretization of GMsFEM Using Deep Learning

Abstract

In this paper, we propose a deep-learning-based approach to a class of multiscale problems. The generalized multiscale finite element method (GMsFEM) has been proven successful as a model reduction technique of flow problems in heterogeneous and high-contrast porous media. The key ingredients of GMsFEM include mutlsicale basis functions and coarse-scale parameters, which are obtained from solving local problems in each coarse neighborhood. Given a fixed medium, these quantities are precomputed by solving local problems in an offline stage, and result in a reduced-order model. However, these quantities have to be re-computed in case of varying media (various permeability fields). The objective of our work is to use deep learning techniques to mimic the nonlinear relation between the permeability field and the GMsFEM discretizations, and use neural networks to perform fast computation of GMsFEM ingredients repeatedly for a class of media. We provide numerical experiments to investigate the predictive power of neural networks and the usefulness of the resultant multiscale model in solving channelized porous media flow problems.

1. Introduction

Multiscale features widely exist in many engineering problems. For instance, in porous media flow, the media properties typically vary over many scales and contain high contrast. Multiscale finite element methods (MsFEM) [1,2,3] and generalized multiscale finite element methods (GMsFEM) [4,5] are designed for solving multiscale problems using local model reduction techniques. In these methods, the computational domain is partitioned into a coarse grid ${\mathcal{T}}^H$, which does not necessarily resolve all multiscale features. We further perform a refinement of ${\mathcal{T}}^H$ to obtain a fine grid ${\mathcal{T}}^h$, which essentially resolves all multiscale features. The idea of local model reduction in these methods is based on idenfications of local multiscale basis functions supported in coarse regions on the fine grid, and replacement of the macroscopic equations by a coarse-scale system using a limited number of local multiscale basis functions. As in many model reduction techniques, the computations of multiscale basis functions, which constitute a small dimensional subspace, can be performed in an offline stage. For a fixed medium, these multiscale basis functions are reusable for any force terms and boundary conditions. Therefore, these methods provide substantial computational savings in the online stage, in which a coarse-scale system is constructed and solved on the reduced-order space.

However, difficulties arise in situations with uncertainties in the media properties in some local regions, which are common for oil reservoirs or aquifers. One straightforward approach for quantifying the uncertainties is to sample realizations of media properties. In such cases, it is challenging to find an offline principal component subspace which is able to universally solve the multiscale problems with different media properties. The computation of multiscale basis functions has to be performed in an online procedure for each medium. Even though the multiscale basis functions are reusable for different force terms and boundary conditions, the computational effort can grow very huge for a large number of realizations of media properties. To this end, building a functional relationship between the media properties and the multiscale model in an offline stage can avoid repeating expensive computations and thus vastly reduce the computational complexity. Due to the diversity of complexity of the media properties, the functional relationship is highly nonlinear. Modelling such a nonlinear functional relationship typically involves high-order approximations. Therefore, it is natural to use machine learning techniques to devise such complex models. In [6,7], the authors make use of a Bayesian approach for learning multiscale models and incorporating essential observation data in the presence of uncertainties.

Deep neural networks is one class of machine learning algorithm that is based on an artificial neural network, which is composed of a relatively large number of layers of nonlinear processing units, called neurons, for feature extraction. The neurons are connected to other neurons in the successive layers. The information propagates from the input, through the intermediate hidden layers, and to the output layer. In the propagation process, the output in each layer is used as input in the consecutive layer. Each layer transforms its input data into a little more abstract feature representation. In between layers, a nonlinear activation function is used as the nonlinear transformation on the input, which increases the expressive power of neural networks. Recently, deep neural network (DNN) has been successfully used to interpret complicated data sets and applied to tasks with pattern recognition, such as image recognition, speech recognition and natural language processing [8,9,10]. Extensive researches have also been conducted on investigating the expression power of deep neural networks [11,12,13,14,15].

Results show that neural networks can represent and approximate a large class of functions. Recently, deep learning has been applied to model reductions and partial differential equations. In [16], the authors studied deep convolution networks for surrogate model construction. on dynamic flow problems in heterogeneous media. In [17], the authors studied the relationship between residual networks (ResNet) and characteristic equations of linear transport, and proposed an interpretation of deep neural networks by continuous flow models. In [18], the authors combined the idea of the Ritz method and deep learning techniques to solve elliptic problems and eigenvalue problems. In [19], a neural network has been designed to learn the physical quantities of interest as a function of random input coefficients. The concept of using deep learning to generate a reduced-order model for a dynamic flow has been applied to proper orthogonal decomposition (POD) global model reduction [20] and nonlocal multi-continuum upscaling (NLMC) [21].

In this work, we propose a deep-learning-based method for fast computation of the GMsFEM discretization. Our approach makes use of deep neural networks as a fast proxy to compute GMsFEM discretizations for flow problems in channelized porous media with uncertainties. More specifically, neural networks are used to express the functional relationship between the media properties and the multiscale model. Such networks are built up in an offline stage. Sufficient sample pairs are required to ensure the expressive power of the networks. With different realizations of media properties, one can use the built network and avoid computations of local problems and spectral problems.

The paper is organized as follows. We start with the underlying partial differential equation that describes the flow within a heterogeneous media and the main ingredients of GMsFEM in Section 2. Next, in Section 3, we present the idea of using deep learning as a proxy for prediction of GMsFEM discretizations. The networks will be precisely defined and the sampling will be explained in detail. In Section 4, we present numerical experiments to show the effectiveness of our presented networks on several examples with different configurations. Finally, a conclusive discussion is provided in Section 5.

2. Preliminaries

In this paper, we are considering the flow problem in highly heterogeneous media

$$\begin{array}{cc}\hfill -\mathrm{div}\left(\kappa \nabla u\right)=f& \mathrm{in}\mathsf{\Omega },\hfill \\ \hfill u=0\mathrm{or}\frac{\partial u}{\partial n}=0& \mathrm{on}\partial \mathsf{\Omega },\hfill \end{array}$$

(1)

where $\mathsf{\Omega }$ is the computational domain, $\kappa$ is the permeability coefficient in $L^{\infty }(\mathsf{\Omega })$, and f is a source function in $L^2(\mathsf{\Omega })$. We assume the coefficient $\kappa$ is highly heterogeneous with high contrast. The classical finite element method for solving (1) numerically is given by: find $u_h\in V_h$ such that

$$\begin{array}{cc}\hfill a(u_h,v)={\int }_{\mathsf{\Omega }}\kappa \nabla u_h·\nabla vdx& ={\int }_{\mathsf{\Omega }}fvdx=(f,v)\mathrm{for}\mathrm{all}v\in V_h,\hfill \end{array}$$

(2)

where $V_h$ is a standard conforming finite element space over a partition ${\mathcal{T}}_h$ of $\mathsf{\Omega }$ with mesh size h.

However, with the highly heterogeneous property of coefficient $\kappa$, the mesh size h has to be taken extremely small to capture the underlying fine-scale features of $\kappa$. This ends up with a large computational cost. GMsFEM [4,5] serves as a model reduction technique to reduce the number of degrees of freedom and attain both efficiency and accuracy to a considerable extent. GMsFEM has been successfully extended to other formulations and applied to other problems. Here we provide a brief introduction of the main ingredients of GMsFEM. For a more detailed discussion of GMsFEM and related concepts, the reader is referred to [22,23,24,25,26].

In GMsFEM, we define a coarse mesh ${\mathcal{T}}^H$ over the domain $\mathsf{\Omega }$ and refine to obtain a fine mesh ${\mathcal{T}}^h$ with mesh size $h\ll H$, which is fine enough to restore the multiscale properties of the problem. Multiscale basis functions are defined on coarse grid blocks using linear combinations of finite element basis functions on ${\mathcal{T}}^h$, and designed to resolve the local multiscale behaviors of the exact solution. The multiscale finite element space $V_{\mathrm{ms}}$, which is a principal component subspace of the conforming finite space $V_h$ with $\dim (V_{\mathrm{ms}})\ll \dim (V_h)$, is constructed by the linear span of multiscale basis functions. The multiscale solution $u_{\mathrm{ms}}\in V_{\mathrm{ms}}$ is then defined by

$$a(u_{\mathrm{ms}},v)=(f,v)\mathrm{for}\mathrm{all}v\in V_{\mathrm{ms}}.$$

(3)

We consider the identification of dominant modes for solving (1) by multiscale basis functions, including spectral basis functions and simplified basis functions, in GMsFEM. Here, we present the details of the construction of multiscale basis functions in GMsFEM. Let $N_x=\{x_i|1\le i\le N_v\}$ be the set of nodes of the coarse mesh ${\mathcal{T}}^H$. For each coarse grid node $x_i\in N_x$, the coarse neighborhood ${\omega }_i$ is defined by

$${\omega }_i=\bigcup \{K_j\in {\mathcal{T}}^H;x_i\in {\overline{K}}_j\},$$

(4)

that is, the union of the coarse elements $K_j\in {\mathcal{T}}^H$ containing the coarse grid node $x_i$. An example of the coarse and fine mesh, coarse blocks and a coarse neighborhood is shown in Figure 1. For each coarse neighbourhood ${\omega }_i$, we construct multiscale basis functions ${\{{\varphi }_j^{{\omega }_i}\}}_{j=1}^{L_i}$ supported on ${\omega }_i$.

For the construction of spectral basis functions, we first construct a snapshot space $V_{\mathrm{snap}}^{(i)}$ spanned by local snapshot basis functions ${\psi }_{\mathrm{snap}}^{i,k}$ for each local coarse neighborhood ${\omega }_i$. The snapshot basis function ${\psi }_{\mathrm{snap}}^{i,k}$ is the solution of a local problem

$$\begin{array}{cc}\hfill -\mathrm{div}(\kappa \nabla {\psi }_{\mathrm{snap}}^{i,k})=0,& \mathrm{in}{\omega }_i\hfill \\ \hfill {\psi }_{\mathrm{snap}}^{i,k}={\delta }_{i,k},& \mathrm{on}\partial {\omega }_i.\hfill \end{array}$$

(5)

The fine grid function ${\delta }_{i,k}$ is a function defined for all $x_s\in \partial {\omega }_i$, where $\{x_s\}$ denote the fine degrees of freedom on the boundary of local coarse region ${\omega }_i$. Specifically,

$${\delta }_{i,k}(x_s)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}s=k\hfill \\ 0\hfill & \mathrm{if}s\ne k.\hfill \end{array}\right.$$

(6)

The linear span of these harmonic extensions forms the local snapshot space $V_{\mathrm{snap}}^{(i)}={\mathrm{span}}_k\{{\psi }_{\mathrm{snap}}^{i,k}\}$. One can also use randomized boundary conditions to reduce the computational cost associated with snapshot calculations [25]. Next, a spectral problem is designed based on our analysis and used to reduce the dimension of local multiscale space. More precisely, we seek for eigenvalues ${\lambda }_m^i$ and corresponding eigenfunctions ${\varphi }_m^{{\omega }_i}\in V_{\mathrm{snap}}^{(i)}$ satisfying

$$a_i({\varphi }_m^{{\omega }_i},v)={\lambda }_m^i{s}_i({\varphi }_m^{{\omega }_i},v),\forall v\in V_{\mathrm{snap}}^{(i)},$$

(7)

where the bilinear forms in the spectral problem are defined as

$$\begin{array}{cc}\hfill a_i(u,v)& ={\int }_{{\omega }_i}\kappa \nabla u·\nabla v,\hfill \\ \hfill s_i(u,v)& ={\int }_{{\omega }_i}\tilde{\kappa }uv,\hfill \end{array}$$

(8)

where $\tilde{\kappa }={\sum }_j\kappa {|\nabla {\chi }_j|}^2$, and ${\chi }_j$ denotes the multiscale partition of the unity function. We arrange the eigenvalues ${\lambda }_m^i$ of the spectral problem (7) in ascending order, and select the first $l_i$ eigenfunctions ${\{{\varphi }_m^{{\omega }_i}\}}_{m=1}^{l_i}$ corresponding to the small eigenvalues as the multiscale basis functions.

An alternative way to construct the multiscale basis function is by using the idea of simplified basis functions. This approach assumes the number of channels and position of the channelized permeability field are known. Therefore we can obtain multiscale basis functions ${\{{\varphi }_m^{{\omega }_i}\}}_{m=1}^{l_i}$ using these information and without solving the spectral problem [27].

Once the multiscale basis functions are constructed, the span of the multiscale basis functions will form the offline space

$$\begin{array}{cc}\hfill V_{ms}^{(i)}& =\mathrm{span}{\{{\varphi }_m^{{\omega }_i}\}}_{m=1}^{l_i},\hfill \\ \hfill V_{ms}& ={\oplus }_i{V}_{ms}^{(i)}.\hfill \end{array}$$

(9)

We will then seek a multlscale solution $u_{ms}\in V_{ms}$ satisfying

$$a(u_{ms},v)=(f,v)\mathrm{for}\mathrm{all}v\in V_{ms},$$

(10)

which is a Galerkin projection of the (1) onto $V_{ms}$, and can be written as a system of linear equations

$$A_c{u}_c=b_c,$$

(11)

where $A_c$ and $b_c$ are the coarse-scale stiffness matrix and load vector. If we collect all the multiscale basis functions and arrange the fine-scale coordinate representation in columns, we obtain the downscaling operator R. Then the fine-scale representation of the multiscale solution is given by

$$u_{ms}=R{u}_c.$$

(12)

3. Deep Learning for GMsFEM

In applications, there are uncertainties within some local regions of the permeability field $\kappa$ in the flow problem. Thousands of forward simulations are needed to quantify the uncertainties of the flow solution. GMsFEM provides us with a fast solver to compute the solutions accurately and efficiently. Considering that there is a large amount of simulation data, we are interested in developing a method utilizing the existing offline data and reducing direct computational effort later. In this work, we aim at using DNN to model the relationship between heterogeneous permeability coefficient $\kappa$ and the key ingredients of GMsFEM solver, i.e., coarse scale stiffness matrices and multiscale basis functions. When the relation is built up, we can feed the network any realization of the permeability field and obtain the corresponding GMsFEM ingredients, and further restore fine-grid GMsFEM solution of (1). The general idea of utilizing deep learning in the GMsFEM framework is illustrated in Figure 2.

Suppose that there are uncertainties for the heterogeneous coefficient in a local coarse block $K_0$, which we call the target block, and the permeability outside the target block remains the same. For example, for a channelized permeability field, the position, location and the permeability values of the channels in the target block can vary. The target block $K_0$ is inside three coarse neighborhoods, denoted by ${\omega }_1,{\omega }_2,{\omega }_3$. The union of the 3 neighborhoods, i.e.,

$${\omega }^{+}(K_0)={\omega }_1\cup {\omega }_2\cup {\omega }_3,$$

(13)

are constituted of by the target block $K_0$ and 12 other coarse blocks, denoted by ${\{K_l\}}_{l=1}^{12}$ A target block and its surrounding neighborhoods are depicted in Figure 3.

For a fixed permeability field $\kappa$, one can compute the multiscale basis functions ${\varphi }_m^{{\omega }_i}(\kappa )$ defined by (7), for $i=1,2,3$, and the local coarse-scale stiffness matrices $A_c^{K_l}(\kappa )$, defined by

$${\left[A_c^{K_l}(\kappa )\right]}_{m,n}^{i,j}={\int }_{K_l}\kappa \nabla {\varphi }_m^{{\omega }_i}(\kappa )·\nabla {\varphi }_n^{{\omega }_j}(\kappa ),$$

(14)

for $l=0,1,\dots ,12$. We are interested in constructing the maps $g_B^{m,i}$ and $g_M^l$, where

• $g_B^{m,i}$ maps the permeability coefficient $\kappa$ to a local multiscale basis function ${\varphi }_m^{{\omega }_i}$, where i denotes the index of the coarse block, and m denotes the index of the basis in coarse block ${\omega }_i$

  $$g_B^{m,i}:\kappa \mapsto {\varphi }_m^{{\omega }_i}(\kappa ),$$

  (15)
• $g_M^l$ maps the permeability coefficient $\kappa$ to the coarse grid parameters $A_c^{K_l}$ ($l=0,\cdots ,12$)

  $$g_M^l:\kappa \mapsto A_c^{K_l}(\kappa ).$$

  (16)

In this work, our goal is to make use of deep learning to build fast approximations of these quantities associated with the uncertainties in the permeability field $\kappa$, which can provide fast and accurate solutions to the heterogeneous flow problem (1).

For each realization $\kappa$, one can compute the images of $\kappa$ under the local multiscale basis maps $g_B^{m,i}$ and the local coarse-scale matrix maps $g_M^l$. These forward calculations serve as training samples for building a deep neural network for approximation of the corresponding maps, i.e.,

$$\begin{array}{cc}\hfill {\mathcal{N}}_B^{m,i}(\kappa )& \approx g_B^{m,i}(\kappa ),\hfill \\ \hfill {\mathcal{N}}_M^l(\kappa )& \approx g_M^l(\kappa ).\hfill \end{array}$$

(17)

In our networks, the permeability field $\kappa$ is the input, while the multiscale basis functions ${\varphi }_m^{{\omega }_i}$ and the coarse-scale matrix $A_c^{K_l}$ are the outputs. Once the neural networks are built, we can use the networks to compute the multiscale basis functions and coarse-scale parameters in the associated region for any permeability field $\kappa$. Using these local information from the neural networks together with the global information which can be pre-computed, we can form the downscale operator R with the multiscale basis functions, form and solve the linear system (11), and obtain the multiscale solution by (12).

3.1. Network Architecture

In general, an L-layer neural network $\mathcal{N}$ can be written in the form

$$\mathcal{N}(x;\theta )=\sigma (W_L\sigma (\cdots \sigma (W_2\sigma (W_{1}x+b_1)+b_2)\cdots )+b_L),$$

(18)

where $\theta :=(W_1,W_2,\cdots ,W_L,b_1,b_2,\cdots ,b_L)$, W’s are the weight matrices and b’s are the bias vectors, $\sigma$ is the activation function, x is the input. Such a network is used to approximate the output y. Our goal is then to find ${\theta }^{\ast }$ by solving an optimization problem

$${\theta }^{\ast }=\underset{\theta }{\mathrm{argmin}}\mathcal{L}(\theta ),$$

(19)

where $\mathcal{L}(\theta )$ is called loss function, which measures the mismatch between the image of the input x under the the neural network $\mathcal{N}(x,y;\theta )$ and the desired output y in a set of training samples $(x_j,y_j)$. In this paper, we use the mean-squared error metric to be our loss function

$$\mathcal{L}(\theta )=\frac{1}{N}\sum _{j=1}^N{\parallel y_j-\mathcal{N}(x_j;\theta )\parallel }_2^2,$$

(20)

where N is the number of the training samples. An illustration of a deep neural network is shown in Figure 4.

Suppose we have a set of different realizations of the permeability $\{{\kappa }_1,{\kappa }_2,\cdots ,{\kappa }_N\}$ in the target block. In our network, the input $x_j={\kappa }_j\in {\mathbb{R}}^d$ is a vector containing the permeability image pixels in the target block. The output $y_j$ is an entry of the local stiffness matrix $A_c^{K_l}$, or the coordinate representation of a multiscale basis function ${\varphi }_m^{{\omega }_i}$. We will make use of these sample pairs $(x_j,y_j)$ to train a deep neural network ${\mathcal{N}}_B^{m,i}(x;{\theta }_B^{\ast })$ and ${\mathcal{N}}_M^l(x;{\theta }_M^{\ast })$ by minimizing the loss function with respect to the network parameter $\theta$, such that the trained neural networks can approximate the functions $g_B^{m,i}$ and $g_M^l$, respectively. Once the neural is constructed, for some given new permeability coefficient ${\kappa }_{N+1}$, we use our trained networks to compute a fast prediction of the outputs, i.e., local multiscale basis functions ${\varphi }_m^{{\omega }_i,\mathrm{pred}}$ by

$${\varphi }_m^{{\omega }_i,\mathrm{pred}}({\kappa }_{N+1})={\mathcal{N}}_B^{m,i}({\kappa }_{N+1};{\theta }_B^{\ast })\approx g_B^{m,i}({\kappa }_{N+1})={\varphi }_m^{{\omega }_i}({\kappa }_{N+1}),$$

(21)

and local coarse-scale stiffness matrix $A_c^{K_l,\mathrm{pred}}$ by

$$A_c^{K_l,\mathrm{pred}}({\kappa }_{N+1})={\mathcal{N}}_M^l({\kappa }_{N+1};{\theta }_M^{\ast })\approx g_M^l({\kappa }_{N+1})=A_c^{K_l}({\kappa }_{N+1}).$$

(22)

3.2. Network-Based Multiscale Solver

Once the neural networks are built, we can assemble the predicted multiscale basis functions to obtain a prediction $R^{\mathrm{pred}}$ for the downscaling operator, and assemble the predicted local coarse-scale stiffness matrix $A_c^{K_l,\mathrm{pred}}$ in the global matrix $A_c^{\mathrm{pred}}$. Following (11) and (12), we solve the predicted coarse-scale coefficient vector $u_c^{\mathrm{pred}}$ from the following linear system

$$A_c^{\mathrm{pred}}{u}_c^{\mathrm{pred}}=b_c,$$

(23)

and obtain the predicted multiscale solution $u_{ms}^{\mathrm{pred}}$ by

$$u_{ms}^{\mathrm{pred}}=R^{\mathrm{pred}}{u}_c^{\mathrm{pred}}.$$

(24)

4. Numerical Results

In this section, we present some numerical results for predicting the GMsFEM ingredients and solutions using our proposed method. We considered permeability fields $\kappa$ with high-contrast channels inside the domain $\mathsf{\Omega }={(0,1)}^2$, which consist of uncertainties in a target cell $K_0$. More precisely, we considered a number of random realizations of permeability fields ${\kappa }_1,{\kappa }_2,{\kappa }_3,\cdots ,{\kappa }_{N+M}$. Each permeability field contained two high-conductivity channels, and the fields differ in the target cell $K_0$ by:

• in experiment 1, the channel configurations were all distinct, and the permeability coefficients inside the channels were fixed in each sample (see Figure 5 for illustrations), and
• in experiment 2, the channel configurations were randomly chosen among five configurations, and the permeability coefficients inside the channels followed a random distribution (see Figure 6 for illustrations).

In these numerical experiments, we assumed there were uncertainties in only the target block $K_0$. The permeability field in $\mathsf{\Omega }\setminus K_0$ was fixed across all the samples.

We followed the procedures in Section 3 and generated sample pairs using GMsFEM. Local spectral problems were solved to obtain the multiscale basis functions ${\varphi }_m^{{\omega }_i}$. In the neural network, the permeability field $x=\kappa$ was considered to be the input, while the local multiscale basis functions $y={\varphi }_m^{{\omega }_i}$ and local coarse-scale matrices $y=A_c^{K_l}$ were regarded as the output. These sample pairs were divided into the training set and the learning set in a random manner. A large number N of realizations, namely ${\kappa }_1,{\kappa }_2,\dots ,{\kappa }_N$, were used to generate sample pairs in the training set, while the remaining M realizations, namely, ${\kappa }_{N+1},{\kappa }_{N+2},\dots ,{\kappa }_{N+M}$ are used in testing the predictive power of the trained network. We remark that, for each basis function and each local matrix, we solved an optimization problem in minimizing the loss function defined by the sample pairs in the training set, and build a separate deep neural network. We summarize the network architectures for training local coarse scale stiffness matrix and multiscale basis functions as below:

• For the multiscale basis function ${\varphi }_m^{{\omega }_i}$, we built a network ${\mathcal{N}}_B^{m,i}$ using

  –

  Input: vectorized permeability pixels values $\kappa$,

  –

  Output: coefficient vector of multiscale basis ${\varphi }_m^{{\omega }_i}(\kappa )$ on coarse neighborhood ${\omega }_i$,

  –

  Loss function: mean squared error ${\displaystyle \frac{1}{N}\sum _{j=1}^N\left|\right|{\varphi }_m^{{\omega }_i}({\kappa }_j)-{\mathcal{N}}_B^{m,i}({\kappa }_j;{\theta }_B){\left|\right|}_2^2}$,

  –

  Activation function: leaky ReLu function,

  –

  DNN structure: 10–20 hidden layers, each layer have 250–350 neurons,

  –

  Training optimizer: Adamax.

• For the local coarse scale stiffness matrix $A_c^{K_l}$, we build a network ${\mathcal{N}}_M^l$ using

  –

  Input: vectorized permeability pixels values $\kappa$,

  –

  Output: vectorized coarse scale stiffness matrix $A_c^{K_l}(\kappa )$ on the coarse block $K_l$,

  –

  Loss function: mean squared error ${\displaystyle \frac{1}{N}\sum _{j=1}^N\left|\right|A_c^{K_l}({\kappa }_j)-{\mathcal{N}}_M^l({\kappa }_j;{\theta }_M){\left|\right|}_2^2}$,

  –

  Activation function: ReLu function (rectifier),

  –

  DNN structure: 10–16 hidden layers, each layer have 100–500 neurons,

  –

  Training optimizer: Proximal Adagrad.

For simplicity, the activation functions ReLU function [28] and Leaky ReLU function were used as they have the simplest derivatives among all nonlinear functions. The ReLU function proved to be useful in training deep neural network architectures. The Leaky ReLU function can resolve the vanishing gradient problem which can accelerate the training in some occasions. The optimizers Adamax and Proximal Adagrad are stochastic gradient descent (SGD)-based methods commonly used in neural network training [29]. In both experiments, we trained our network using Python API Tensorflow and Keras [30].

Once a neural network was built on training, it can be used to predict the output given a new input. The accuracy of the predictions is essential in making the network useful. In our experiments, we used M sample pairs, which were not used in training the network, to examine the predictive power of our network. On these sample pairs, referred to as the testing set, we compared the prediction and the exact output and computed the mismatch in some suitable metric. Here, we summarize the metric used in our numerical experiment. For the multiscale basis functions, we compute the relative error in $L^2$ and $H^1$ norm, i.e.,

$$\begin{array}{cc}\hfill e_{L^2}({\kappa }_{N+j})& ={\left({\displaystyle \frac{{\int }_{\mathsf{\Omega }}{\left|{\varphi }_m^{{\omega }_i}({\kappa }_{N+j})-{\varphi }_m^{{\omega }_i,\mathrm{pred}}({\kappa }_{N+j})\right|}^2}{{\int }_{\mathsf{\Omega }}{\left|{\varphi }_m^{{\omega }_i}({\kappa }_{N+j})\right|}^2}}\right)}^{\frac{1}{2}},\hfill \\ \hfill e_{H^1}({\kappa }_{N+j})& ={\left({\displaystyle \frac{{\int }_{\mathsf{\Omega }}{\left|\nabla {\varphi }_m^{{\omega }_i}({\kappa }_{N+j})-\nabla {\varphi }_m^{{\omega }_i,\mathrm{pred}}({\kappa }_{N+j})\right|}^2}{{\int }_{\mathsf{\Omega }}{\left|\nabla {\varphi }_m^{{\omega }_i}({\kappa }_{N+j})\right|}^2}}\right)}^{\frac{1}{2}}.\hfill \end{array}$$

(25)

For the local stiffness matrices, we computed the relative error in entrywise ${\ell }^2$, entrywise ${\ell }^{\infty }$ and Frobenius norm, i.e.,

$$\begin{array}{cc}\hfill e_{{\ell }^2}({\kappa }_{N+j})& ={\displaystyle \frac{\parallel A_c^{K_l}({\kappa }_{N+j})-A_c^{K_l,\mathrm{pred}}({\kappa }_{N+j}){\parallel }_2}{\parallel A_c^{K_l}({\kappa }_{N+j}){\parallel }_2}},\hfill \\ \hfill e_{{\ell }^{\infty }}({\kappa }_{N+j})& ={\displaystyle \frac{\parallel A_c^{K_l}({\kappa }_{N+j})-A_c^{K_l,\mathrm{pred}}({\kappa }_{N+j}){\parallel }_{\infty }}{\parallel A_c^{K_l}({\kappa }_{N+j}){\parallel }_{\infty }}},\hfill \\ \hfill e_F({\kappa }_{N+j})& ={\displaystyle \frac{\parallel A_c^{K_l}({\kappa }_{N+j})-A_c^{K_l,\mathrm{pred}}({\kappa }_{N+j}){\parallel }_F}{\parallel A_c^{K_l}({\kappa }_{N+j}){\parallel }_F}}.\hfill \end{array}$$

(26)

A more important measure of the usefulness of the trained neural network is the predicted multiscale solution $u_{ms}^{\mathrm{pred}}(\kappa )$ given by (23) and (24). We compared the predicted solution to $u_{ms}$ defined by (11) and (12), and computed the relative error in $L^2$ and energy norm, i.e.,

$$\begin{array}{cc}\hfill e_{L^2}({\kappa }_{N+j})& ={\left({\displaystyle \frac{{\int }_{\mathsf{\Omega }}{\left|u_{ms}({\kappa }_{N+j})-u_{ms}^{\mathrm{pred}}({\kappa }_{N+j})\right|}^2}{{\int }_{\mathsf{\Omega }}{\left|u_{ms}({\kappa }_{N+j})\right|}^2}}\right)}^{\frac{1}{2}},\hfill \\ \hfill e_a({\kappa }_{N+j})& ={\left({\displaystyle \frac{{\int }_{\mathsf{\Omega }}{\kappa }_j{\left|\nabla u_{ms}({\kappa }_{N+j})-\nabla u_{ms}^{\mathrm{pred}}({\kappa }_{N+j})\right|}^2}{{\int }_{\mathsf{\Omega }}{\kappa }_j{\left|\nabla u_{ms}({\kappa }_{N+j})\right|}^2}}\right)}^{\frac{1}{2}}.\hfill \end{array}$$

(27)

4.1. Experiment 1

In this experiment, we considered curved channelized permeability fields. Each permeability field contained a straight channel and a curved channel. The straight channel was fixed and the curved channel struck the boundary of the target cell $K_0$ at the same points. The curvature of the sine-shaped channel inside $K_0$ varied among all realizations. We generated 2000 realizations of permeability fields, where the permeability coefficients were fixed. Samples of permeability fields are depicted in Figure 5. Among the 2000 realizations, 1980 sample pairs were randomly chosen and used as training samples, and the remaining 20 sample pairs were used as testing samples.

For each realization, we computed the local multiscale basis functions and local coarse-scale stiffness matrix. In building the local snapshot space, we solved for harmonic extension of all the fine-grid boundary conditions. Local multiscale basis functions were then constructed by solving the spectral problem and multiplied the spectral basis functions with the multiscale partition of unity functions. With the offline space constructed, we computed the coarse-scale stiffness matrix. We used the training samples to build deep neural networks for approximating these GMsFEM quantities, and examined the performance of the approximations on the testing set.

Table 1. Percentage error of multiscale basis functions ${\varphi }_1^{{\omega }_i}$ in experiment 1.

Sample	${\mathit{\omega }}_1$		${\mathit{\omega }}_2$		${\mathit{\omega }}_3$
$\mathit{j}$	${\mathit{e}}_{{\mathit{L}}^{\mathbf{2}}}$	${\mathit{e}}_{{\mathit{H}}^{\mathbf{1}}}$	${\mathit{e}}_{{\mathit{L}}^{\mathbf{2}}}$	${\mathit{e}}_{{\mathit{H}}^{\mathbf{1}}}$	${\mathit{e}}_{{\mathit{L}}^{\mathbf{2}}}$	${\mathit{e}}_{{\mathit{H}}^{\mathbf{1}}}$
1	0.47%	3.2%	0.40%	3.6%	0.84%	5.1%
2	0.45%	4.4%	0.39%	3.3%	1.00%	6.3%
3	0.34%	2.3%	0.40%	3.1%	0.88%	4.3%
4	0.35%	4.2%	0.43%	5.4%	0.94%	6.6%
5	0.35%	3.3%	0.37%	3.9%	0.90%	6.1%
6	0.51%	4.7%	0.92%	12.0%	2.60%	19.0%
7	0.45%	4.1%	0.38%	3.2%	1.00%	6.4%
8	0.31%	3.4%	0.43%	5.5%	1.10%	7.7%
9	0.25%	2.2%	0.46%	5.6%	1.10%	6.2%
10	0.31%	3.5%	0.42%	4.5%	1.30%	7.6%
Mean	0.38%	3.5%	0.46%	5.0%	1.17%	7.5%

,

Table 2. Percentage error of multiscale basis functions ${\varphi }_2^{{\omega }_i}$ in experiment 1.

Sample	${\mathit{\omega }}_1$		${\mathit{\omega }}_2$		${\mathit{\omega }}_3$
$\mathit{j}$	${\mathit{e}}_{{\mathit{L}}^{\mathbf{2}}}$	${\mathit{e}}_{{\mathit{H}}^{\mathbf{1}}}$	${\mathit{e}}_{{\mathit{L}}^{\mathbf{2}}}$	${\mathit{e}}_{{\mathit{H}}^{\mathbf{1}}}$	${\mathit{e}}_{{\mathit{L}}^{\mathbf{2}}}$	${\mathit{e}}_{{\mathit{H}}^{\mathbf{1}}}$
1	0.47%	4.2%	0.40%	1.4%	0.32%	1.1%
2	0.57%	3.2%	0.31%	1.4%	0.30%	1.1%
3	0.58%	2.7%	0.31%	1.4%	0.33%	1.1%
4	0.59%	3.6%	0.13%	1.3%	0.32%	1.1%
5	0.53%	4.0%	0.51%	1.6%	0.27%	1.0%
6	0.85%	4.3%	0.51%	2.1%	0.29%	1.3%
7	0.50%	2.7%	0.22%	1.5%	0.29%	1.0%
8	0.43%	4.5%	0.61%	1.9%	0.35%	1.1%
9	0.71%	2.9%	0.14%	1.4%	0.27%	1.1%
10	0.66%	4.4%	0.53%	1.8%	0.26%	1.1%
Mean	0.59%	3.6%	0.37%	1.6%	0.30%	1.1%

and

Table 3. Percentage error of the local stiffness matrix $A_c^{K_0}$ in experiment 1.

Sample j	${\mathit{e}}_{{\mathit{\ell }}^2}$	${\mathit{e}}_{\mathit{F}}$
1	0.67%	0.84%
2	0.37%	0.37%
3	0.32%	0.38%
4	1.32%	1.29%
5	0.51%	0.59%
6	4.43%	4.28%
7	0.34%	0.38%
8	0.86%	1.04%
9	1.00%	0.97%
10	0.90%	1.08%
Mean	0.76%	0.81%

record the error of the prediction by the neural networks in each testing sample and the mean error measured in the defined metric. It can be seen that the predictions were of high accuracy. This is vital in ensuring the predicted GMsFEM solver is useful.

Table 4. Percentage error of multiscale solution $u_{ms}$ in experiment 1.

Sample j	${\mathit{e}}_{{\mathit{L}}^2}$	${\mathit{e}}_{\mathit{a}}$
1	0.31%	4.58%
2	0.30%	4.60%
3	0.30%	4.51%
4	0.27%	4.60%
5	0.29%	4.56%
6	0.47%	4.67%
7	0.39%	4.70%
8	0.30%	4.63%
9	0.35%	4.65%
10	0.31%	4.65%
Mean	0.33%	4.62%

records the error of the multiscale solution in each testing sample and the mean error using our proposed method. It can be observed that using the predicted GMsFEM solver, we obtained a good approximation of the multiscale solution compared with the exact GMsFEM solver.

4.2. Experiment 2

In this experiment, we considered sine-shaped channelized permeability fields. Each permeability field contained a straight channel and a sine-shaped channel. There were altogether five channel configurations, where the straight channel was fixed and the sine-shaped channel struck the boundary of the target cell $K_0$ at the same points. The curvature of the sine-shaped channel inside $K_0$ varied among these configurations. For each channel configuration, we generated 500 realizations of permeability fields, where the permeability coefficients followed random distributions. Samples of permeability fields are depicted in Figure 6. Among the 2500 realizations, 2475 sample pairs were randomly chosen and used as training samples, and the remaining 25 sample pairs were used as testing samples.

Next, for each realization, we computed the local multiscale basis functions and local coarse-scale stiffness matrix. In building the local snapshot space, we solved for harmonic extension of randomized fine-grid boundary conditions, so as to reduce the number of local problems to be solved. Local multiscale basis functions were then constructed by solving the spectral problem and multiplied the spectral basis functions with the multiscale partition of unity functions. With the offline space constructed, we computed the coarse-scale stiffness matrix. We used the training samples to build deep neural networks for approximating these GMsFEM quantities, and examined the performance of the approximations on the testing set.

Figure 7, Figure 8 and Figure 9 show the comparison of the multiscale basis functions in two respective coarse neighborhoods. It can be observed that the predicted multiscale basis functions were in good agreement with the exact ones. In particular, the neural network successfully interpreted the high conductivity regions as the support localization feature of the multiscale basis functions.

Table 5. Mean percentage error of multiscale basis functions ${\varphi }_m^{{\omega }_i}$ in experiment 2.

Basis	${\mathit{\omega }}_1$		${\mathit{\omega }}_2$		${\mathit{\omega }}_3$
$\mathit{m}$	${\mathit{e}}_{{\mathit{L}}^{\mathbf{2}}}$	${\mathit{e}}_{{\mathit{H}}^{\mathbf{1}}}$	${\mathit{e}}_{{\mathit{L}}^{\mathbf{2}}}$	${\mathit{e}}_{{\mathit{H}}^{\mathbf{1}}}$	${\mathit{e}}_{{\mathit{L}}^{\mathbf{2}}}$	${\mathit{e}}_{{\mathit{H}}^{\mathbf{1}}}$
1	0.55	0.91	0.37	3.02	0.20	0.63
2	0.80	1.48	2.17	3.55	0.27	1.51

and

Table 6. Percentage error of the local stiffness matrix $A_c^{K_0}$ in experiment 2.

	${\mathit{e}}_{{\mathit{\ell }}^2}$	${\mathit{e}}_{{\mathit{\ell }}^{\mathit{\infty }}}$	${\mathit{e}}_{\mathit{F}}$
Mean	0.75	0.72	0.80

record the mean error of the prediction by the neural networks, measured in the defined metric. Again, it can be seen that the prediction are of high accuracy.

Table 7. Percentage error of multiscale solution $u_{ms}$ in experiment 2.

	${\mathit{e}}_{{\mathit{L}}^2}$	${\mathit{e}}_{\mathit{a}}$
Mean	0.03	0.26

records the mean error between the multiscale solution using the neural-network-based multiscale solver and using exact GMsFEM. we obtain a good approximation of the multiscale solution compared with the exact GMsFEM solver.

5. Conclusions

In this paper, we develop a method using deep learning techniques for fast computation of GMsFEM discretizations. Given a particular permeability field, the main ingredients of GMsFEM, including the multiscale basis functions and coarse-scale matrices, are computed in an offline stage by solving local problems. However, when one is interested in calculating GMsFEM discretizations for multiple choices of permeability fields, repeatedly formulating and solving such local problems could become computationally expensive or even infeasible. Multi-layer networks are used to represent the nonlinear mapping from the fine-scale permeability field coefficients to the multiscale basis functions and the coarse-scale parameters. The networks provide a direct fast approximation of the GMsFEM ingredients in a local neighborhood for any online permeability fields, in contrast to repeatedly formulating and solving local problems. Numerical results are presented to show the performance of our proposed method. We see that, given sufficient samples of GMsFEM discretizations for supervised training, deep neural networks are capable of providing reasonably close approximations of the exact GMsFEM discretization. Moreover, the small consistency error provides good approximations of multiscale solutions.