Hybrid explicit-implicit learning for multiscale problems with time dependent source

The splitting method is a powerful method for solving partial differential equations. Various splitting methods have been designed to separate different physics, nonlinearities, and so on. Recently, a new splitting approach has been proposed where some degrees of freedom are handled implicitly while other degrees of freedom are handled explicitly. As a result, the scheme contains two equations, one implicit and the other explicit. The stability of this approach has been studied. It was shown that the time step scales as the coarse spatial mesh size, which can provide a significant computational advantage. However, the implicit solution part can be expensive, especially for nonlinear problems. In this paper, we introduce modified partial machine learning algorithms to replace the implicit solution part of the algorithm. These algorithms are first introduced in ‘HEI: Hybrid Explicit-Implicit Learning For Multiscale Problems’, where a homogeneous source term is considered along with the Transformer, which is a neural network that can predict future dynamics. In this paper, we consider time-dependent source terms which is a generalization of the previous work. Moreover, we use the whole history of the solution to train the network. As the implicit part of the equations is more complicated to solve, we design a neural network to predict it based on training. Furthermore, we compute the explicit part of the solution using our splitting strategy. In addition, we use Proper Orthogonal Decomposition based model reduction in machine learning. The machine learning algorithms provide computational saving without sacrificing accuracy. We present three numerical examples which show that our machine learning scheme is stable and accurate.

Many physical problems vary at multiple space and time scales. The solutions of these problems are difficult to compute as we need to resolve the time and space scales. For example, porous media flow and transport occur over many space and time scales. Many multiscale and splitting algorithms have been proposed to solve these problems [28, 18, 41]. In this paper, we would like to combine machine learning and splitting/multiscale approaches to predict the dynamics of complex forward problems. Our goal is to perform machine learning on a small part of the solution space, while computing “the rest” of the solution fast. Our proposed splitting algorithms allow doing it, which we will describe below.

Our approaches use multiscale methods in their construction to achieve coarser time stepping. Multiscale methods have extensively been studied in the literature. For linear problems, homogenization-based approaches [23, 39, 40], multiscale finite element methods [23, 31, 35], generalized multiscale finite element methods (GMsFEM) [12, 13, 14, 17, 21], Constraint Energy Minimizing GMsFEM (CEM-GMsFEM) [15, 16, 7], nonlocal multi-continua (NLMC) approaches [19], metric-based upscaling [44], heterogeneous multiscale method [20], localized orthogonal decomposition (LOD) [30], equation-free approaches [45, 46], multiscale stochastic approaches [33, 34, 32], and hierarchical multiscale methods [5] have been studied. Approaches such as GMsFEM and NLMC are designed to handle high-contrast problems. For nonlinear problems, we can replace linear multiscale basis functions with nonlinear maps [26, 27, 22].

Splitting methods are introduced when solving dynamic problems [3, 47, 29]. There is a wide range of applications for splitting methods, e.g. splitting physics, scales, and so on. Splitting methods allow simplifying the large-scale problem into smaller-scale problems. We couple them and obtain the final solution after we solve the smaller-scale problems. Recently, a partially explicit splitting method [28, 18, 24] is proposed for linear and nonlinear multiscale problems. In this approach, we treat some solution parts implicitly while treating the other parts explicitly. The stability analysis shows that the time step only depends on the explicit part of the solution. With the careful splitting of the space, we are able to make the time step scale as the coarse mesh size and independent of contrast.

To illustrate the main idea of our methods, we consider

In this paper, we combine machine learning and splitting approaches (such as the partially explicit method). In the splitting algorithm, we write the solution as $u=u_{1}+u_{2}$, where $u_{1}$ is computed implicitly and $u_{2}$ is computed explicitly. In our previous work, the Transformer has been used to solve this kind of problem in [25, 11]. In [25], we consider zero source term and use the Transformer to predict the future dynamics of the solution. In our current paper, we consider full dynamics and time-dependent source terms. The latter occurs in porous media applications and represents the rates of wells, which vary in time. The goal in these simulations is to compute the solutions at all time steps for a given time-dependent source term in an efficient way, which is not addressed before.

Many machine learning and multiscale methods are designed for solving multiscale parametric PDEs [11, 48, 49, 9, 43, 42]. The advantage of machine learning is that it treats PDEs as an input-to-output map. By training, and possibly, using the observed data, machine learning avoids computing solutions for each input (such as parameters in boundary conditions, permeability fields or sources). In this paper, we use machine learning to learn the difficult-to-compute part of the solution that is, in general, computed implicitly. This is a coarse-grid part of the solution, which can also be observed given some dynamic data. The splitting algorithm gives two coupled equations. The first equation for $u_{1}$ is implicit while the second equation for $u_{2}$ is explicit. We design a neural network to predict the first component of the solution. Then, we calculate the second component using the discretized equations. The inputs of our neural network are the parameters that characterize the source term of the equation. The outputs are the whole history of $u_{1}$, i.e. $u_{1}$ at all time steps. The output dimension of our neural network is the dimension of $u_{1}$ times the number of time steps, which is considerably large. Therefore, we use Proper Orthogonal Decomposition (POD) [48, 10, 38] to reduce the dimension of the solution space. We concatenate all solution vectors of every training input together to generate a large-scale matrix. After doing the singular value decomposition, we select several dominant eigenvectors to form the POD basis vectors. Thus, the output becomes the coordinate vectors based on the POD basis vectors and its dimension is greatly reduced, so is the number of parameters inside the neural networks.

In this paper, we show several numerical examples to illustrate our learning framework which is based on the partially explicit scheme. In our numerical examples, the diffusion term $f(u)$ is $-div(\kappa(x,u)\nabla u)$ and the reaction term is $g(x,t)$. $g(x,t)$ is a source term that is nonzero on several coarse blocks and is zero elsewhere. We choose such $g(x,t)$ as it simulates the wells in petroleum applications. By comparing with the solution which is obtained using fine grid basis functions, we find that the relative $L^{2}$ error is small and similar to the scheme where we calculate both components. Thus, the learning scheme can achieve similar accuracy while saving computation.

This paper is organized as follows. We present preliminaries and review the partially explicit splitting scheme in Section 2. In Section 3, we show details about our learning framework. We present numerical experiments in Section 4.

We consider the following parabolic partial differential equation

Let $V$ be a Hilbert space. In this equation, $f=\cfrac{\delta F}{\delta u}\in V^{\ast}$ with $F$ being the variational derivative of energy $F(u):=\int_{\Omega}E(u)$. In this paper, we use $(\cdot,\cdot)$ to denote $(\cdot,\cdot)_{V^{*},V}$ for simplicity. We assume that $f$ is contrast dependent (linear or nonlinear). Here, $g(x,t)$ is a time dependent source term with parameters. The weak formulation of Equation (2.1) is the following. We want to find $u\in V$ such that

The $(\cdot,\cdot)$ is the $L^{2}$ inner product. We assume that

The standard approach to solving this problem is the finite element method. Let $V_{H}$ be the finite element space and $u_{H}\in V_{H}$ is the numerical solution satisfying

The Backward Euler (implicit) time discretization for this Equation is

The Forward Euler (explicit) time discretization for Equation (2.4) is

One can find details about the stability analysis of the Backward and Forward Euler method in [18].

In Equation (2.1), we consider $g(x,t)$ as a source term with varying parameters representing temporal changes. For every parameter, we want to find the corresponding solution $u(x,t)$. For the traditional finite element methods, we need to construct meshes and basis functions, assemble the stiffness and mass matrices and solve the matrices equations, which makes the computation cost large, especially for the nonlinear PDEs as we have to assemble the stiffness matrix at every time step. Designing machine learning algorithms for the partially explicit scheme can alleviate this problem. Once the neural network is trained, the computation of solution becomes the multiplication of matrices and applying activation functions, which is more efficient in contrast to the finite element methods.

The Equation (2.7) is implicit for $u^{n+1}_{H,1}$, which makes the computations more complicated. We design a machine learning algorithm to predict $u^{n+1}_{H,1}$ and then we compute $u^{n+1}_{H,2}$ using Equation (2.8).

We present a flowchart in Figure 3.1 to show the workflow of our method. We first need to generate the dataset which will be used to train the neural network. Let $w_{m}=(w_{m1},w_{m2},\cdots,w_{mk})$ be the parameters in $g_{w}(x,t)$. We set $w_{mi}\in[L,U]$ for $i=1,2,\cdots,k$ with $L$ and $U$ being the lower bound and upper bound, respectively. We generate the training space $\{w_{m}\}_{m=1}^{M}$ by drawing samples from a uniform distribution over the interval $[L,U]$. For every $w_{m}$ in the training space, we calculate its corresponding partially explicit solution using Equations (2.7) and (2.8). Let $N$ be the number of time steps and $M$ be the size of our training sample space. We denote this solution space as $\{\{u_{H,m}^{n}\}_{n=2}^{N}\}_{m=1}^{M}$, where $u_{H,m}^{n}=u_{H,m,1}^{n}+u_{H,m,2}^{n}$. We remark that the reason why $n$ starts from $2$ is that we need $u^{0}$ and $u^{1}$ to compute $u^{2}$ in the partially explicit scheme. The $u^{1}$ is computed by the Backward Euler scheme using basis functions from both $V_{H,1}$ and $V_{H,2}$. We only need $\{\{u_{H,m,1}^{n}\}_{n=2}^{N}\}_{m=1}^{M}$ for our machine learning algorithm. Before we move on, we discuss a useful tool called Proper Orthogonal Decomposition [6].

For a given set of vectors $\{x_{i}\}_{i=1}^{c}$ with $x_{i}\in\mathbb{R}^{d}$, we want to construct a smaller dimensional vector space that can represent the space spanned by $\{x_{i}\}_{i=1}^{c}$ accurately, i.e., we want to find orthonormal vectors $\{y_{j}\}_{j=1}^{l}$ that solve

where $l\leq r$ and $r$ is the rank of $(x_{1},x_{2},\cdots,x_{c})$. We apply the singular value decomposition for this. Next, we briefly remind the POD procedure [6]. We first define a matrix $X=(x_{1},x_{2},\cdots,x_{c})$ whose dimension is $d\times c$. Then $X^{T}X$ and $XX^{T}$ are $c\times c$ and $d\times d$ dimensional matrices. Consider the eigenvalues and eigenvectors of $X^{T}X$ as well as $XX^{T}$,

Next, we define

We now have

where $\Sigma$ is the diagonal matrix whose main diagonal entries are $\sigma_{i}^{2}$. After doing rearrangement if necessary, we may assume

The solution of (3.1) is given by the $l$ eigenvectors of $XX^{T}$ corresponding to the largest $l$ eigenvalues. Thus, the desired $l$ dimensional subspace is the subspace spanned by eigenvectors corresponding to the largest $l$ eigenvalues.

There are $N\times M$ column vectors ($\in\mathbb{R}^{\dim(V_{H,1})}$) in the solution space $\{\{u_{H,m}^{n}\}_{n=2}^{N}\}_{m=1}^{M}$. We would like to output the solution $\{u_{H}^{n}\}_{n=2}^{N}$ at the same time so that we don’t have to do the iteration for every time step recursively. In this case, the output dimension is $\dim(V_{H,1})\times(N-1)$, which is too large ($\dim(V_{H,1})$ is $300$ in our examples). As a result, the number of parameters in the neural network will also be large, which increases the computation cost. We define

Hence, $\dim(S)=\dim(V_{H,1})\times((N-1)\cdot M)$. After the singular value decomposition to $S$, we define

where $p_{i}\in\mathbb{R}^{\dim(V_{H,1})}$ is the eigenvector corresponding to the $i$-th largest eigenvalue. We use $P$ to get a reduced dimensional coordinate vectors for every $u_{H,m,1}^{n}$,i.e.

Thus, the output dimension is reduced to $l\times(N-1)$. In the testing stage, after we get the output $\{u_{H,1}^{n,\prime}\}_{n=2}^{N}$, we use $P$ to get the coarse grid solution $\{u_{H,1}^{n}=Pu_{H,1}^{n,\prime}\}_{n=2}^{N}$.

As discussed above, in our neural networks, the inputs are the parameters $w$ in $g_{w}(x,t)$ and outputs are $\{u_{H,m,1}^{n,\prime}\}_{n}$. That is

where $\mathcal{N}$ is the neural network. To better illustrate,

where $Y_{1}$, $\phi_{1}$, $w$, $W_{1}$ and $b_{1}$ are the output of the first layer, the activation function, inputs, weight matrix and the biased term. Here, $Y_{1}$, $w$ and $b_{1}$ are row vectors. Similarly,

and

where $I$ is the number of layers and $i=1,2,\cdots,I-1$. Then, we can obtain $\{u_{H,m,1}^{n,\prime}\}_{n}$ by reshaping $Y_{I}$.

Next, we would like to prove two theoretical results. Under the setting in Section 3, we first define the map from the input to the output of our neural network as

where $W$ is a compact subset of $\mathbb{R}^{k}$ with norm $\|\cdot\|_{W}$. We will show that the neural network is capable of approximating $h$ precisely. We then show that when the loss function (mean squared error) goes to zero, the output of the neural network converges to the exact solution (which is computed via the partially explicit scheme). We state the following well-known theorem [36].

The universal approximation theorem relies on the continuity assumption of the mapping $\mathcal{H}$. More specifically, to apply the universal approximation theorem, we now need to show that $h$ is continuous.

The proof is given in Appendix A. With the continuity of $h$, we have the following theorem.

The theorem follows from the continuity of $h$ and the Universal Approximation Theorem. Finally, we present the estimation of the generalization error.

In this section, we present numerical examples. In all cases, the domain $\Omega$ is $[0,1]^{2}$. The permeability field $\kappa$ we use is shown in Figure 4.1. Note that there are high contrast channels in the permeability field. The fine mesh size is $\sqrt{2}/100$ and the coarse mesh size is $\sqrt{2}/10$.

We define the following notation.

• •

  $u_{H,l}$ is obtained by predicting $u_{H,1}$ and computing $u_{H,2}$ by Equation (2.8).

• •

  $u_{H,c}$ is obtained by computing $u_{H,1}$ and $u_{H,2}$ by Equation (2.7) and (2.8).

• •

  $u_{H,f}$ is obtained using the fine grid basis functions and Equation (2.5).

We also need to define the following relative error for better illustration. Let $\|\cdot\|$ be the $L^{2}$ norm. We know that $u_{H,l}=u_{H,l,1}+u_{H,l,2}$ and $u_{H,c}=u_{H,c,1}+u_{H,c,2}$.

• •

  $e_{1}=\cfrac{\|u_{H,l}-u_{H,f}\|}{\|u_{H,f}\|}\times 100\%$, $e_{2}=\cfrac{\|u_{H,c}-u_{H,f}\|}{\|u_{H,f}\|}\times 100\%$.

• •

  $e_{3}=\cfrac{\|u_{H,l,1}-u_{H,c,1}\|}{\|u_{H,c,1}\|}\times 100\%$, $e_{4}=\cfrac{\|u_{H,l}-u_{H,c}\|}{\|u_{H,c}\|}\times 100\%$.

$e_{1}$ and $e_{2}$ are used to measure the difference of accuracy between our learning scheme and the scheme in which we compute both $u_{H,1}$ and $u_{H,2}$. $e_{3}$ and $e_{4}$ are defined to measure how close the learning output is to their target. The neural network is fully connected and is activated by the scaled exponential linear unit (SELU). To obtain the solution of the minimization problem, we use the gradient-based optimizer Adam [37]. The loss function we use is the mean square error (MSE). Our machine learning algorithm is implemented using the open-source library Keras [8] with the Tensorflow [1] backend. The neural network we use is a simple neural network with $5$ fully connected layers. We remark that with this neural network, we can obtain high accuracy.

In the first two examples, we consider the following linear parabolic partial differential equations

where $\vec{n}$ is the outward normal vector and $u_{0}$ is the initial condition which is presented it in Figure 4.2. In our first example, the source term is defined as follows,

where $(w_{1},w_{2},w_{3},w_{4})$ are parameters. In this case, we let $1\leq w_{i}\leq 10$ for $i=1,2,3,4$. We show a simple sketch of $g_{w}(x,t)$ in Figure 4.2 to illustrate. This sketch only shows the position of the coarse blocks where the source term is nonzero. We choose this source, which has some similarities to the wells in petroleum engineering. The final time is $T=0.01$, there are $N=100$ time steps and the time step size is $\Delta t=T/N=10^{-4}$.

The training space contains $1000$ samples and the testing space contains $500$ samples. Every sample’s parameter $w$ is drawn from the uniform distribution over $[1,10]$. The number of POD basis vectors is $l=15$. We choose this $l$ based on the eigenvalue behavior in the process of singular value decomposition. We present $e_{1}$ and $e_{2}$ defined before in Figure 4.3. The error is calculated at each time step as the average over all samples in the testing space. From the plot in Figure 4.3, we find that the curves for $e_{1}$ and $e_{2}$ coincide, which means our learning algorithms can obtain similar accuracy as the scheme in which we compute $u_{H}$ by Equation (2.7) and (2.8). The average of $e_{3}$ and $e_{4}$ over time are $0.151\%$ and $0.155\%$, respectively, which means the output of our learning scheme nearly resembles the corresponding target (the solution we obtain by computing both equations).

In the next case, we test our algorithms using a more complicated source term with more parameters. The source term we use is,

where $(w_{1},w_{2},w_{3},\cdots,w_{10})$ are parameters. In this case, the parameters $w_{i}$ of every sample are drawn from the uniform distribution over $[1,20]$. We present an easy sketch of $g_{w}(x,t)$ in Figure 4.4 to illustrate the positions of the sources. The initial condition $u_{0}$ is also shown in Figure 4.4. The final time $T$, the number of time steps $N$ and the time step size are $0.01$, $100$ and $10^{-4}$, respectively.

The model is trained using 2000 samples and we test it with 1000 samples. The number of POD basis vectors is $l=25$. In Figure 4.5, the $e_{1}$ and $e_{2}$ are presented. By observing the error plot, we find that $e_{1}$ and $e_{2}$ are close to each other, which implies that our learning algorithms can obtain similar accuracy as the scheme where we calculate both $u_{H,1}$ and $u_{H,2}$ using Equation (2.7) and (2.8). The average of $e_{3}$ and $e_{4}$ over time are $0.673\%$ and $0.677\%$, respectively. Such small $e_{3}$ and $e_{4}$ indicate that the prediction is accurate.

In the next case, we consider the following equation

where $\vec{n}$ is the outward normal vector. This equation is nonlinear and its computation cost is much larger than the linear case because we have to assemble the stiffness matrix at every time step. We use the same initial condition $u_{0}$ and source term $g(x,t)$ as the second example. The final time is $T=0.001$, the number of time steps is $N=40$ and the time step size is $\Delta t=T/N=2.5\times 10^{-5}$. We train our neural network using $1000$ samples and the network is tested using $500$ samples. In Figure 4.6, we present $e_{1}$ and $e_{2}$ defined before. From the first plot, we find that the curves for $e_{1}$ and $e_{2}$ nearly coincide, which implies that our machine learning algorithms can have the same accuracy as we compute the solution by Equations (2.7) and (2.8). The average of $e_{3}$ and $e_{4}$ over time are $0.460\%$ and $0.450\%$, respectively, which means our neural network can approximate the target well and achieve a fantastic accuracy.

In this paper, we construct neural networks based on the partially explicit splitting scheme. We predict $u_{1}$, the implicit part of the solution, while solving for $u_{2}$, the explicit part of the solution. The output dimension resulting from $u_{1}$ is large and we apply the POD to reduce the dimension. The POD transformation greatly reduces the output dimension and at the same time the number of parameters in the neural networks. The POD error is small compared to the overall error. We discuss three numerical examples. We present a theoretical justification. We conclude that our machine learning algorithms can obtain similar accuracy as the partially explicit scheme.