Convolutional neural network based reduced order modeling for multiscale problems

The goal of RB methods is to find linear approximation spaces $V_{N}$ for which the worst best-approximation error for elements of $V_{h}$,

$$d_{V_{N}}(V_{h}):=\mathop{\sup}\limits_{u_{h}\in V_{h}}\mathop{\inf}\limits_{\nu\in V_{N}}\|u_{h}-\nu\|,$$

is near the Kolmogorov N-width of $V_{h}$

$$d_{N}(V_{h}):=\mathop{\inf}\limits_{\begin{subarray}{c}W\subseteq V_{h}\\ \dim W\leq N\end{subarray}}\mathop{\sup}\limits_{u_{h}\in V_{h}}\mathop{\inf}\limits_{\nu\in W}\|u_{h}-\nu\|.$$

The existence of the optimal subspace $V_{N}$ is guaranteed by the Theorem II.2.3 in [28].

In practical, the reduced basis functions $P_{N}$ of $V_{N}$ can be generated by either greedy algorithms or proper orthogonal decomposition (POD). We briefly recall the latter as we will use it in the numerical results for comparision. Given a snapshot matrix of data

$$\mathcal{U}=\{u_{h}(K_{1}),u_{h}(K_{2}),\ldots,u_{h}(K_{M})\}$$

the main idea of POD is to find a finite subspace $V_{N}\subseteq V_{h}$ with orthogonal basis $P_{N}:=[p_{1},p_{2},\ldots,p_{N}]$ to minimize the projection error,

$$P_{N}=\mathop{\arg\min}\limits_{\begin{subarray}{c}P^{T}P=I\\ p_{i}\in V_{N}\end{subarray}}\sum_{i=1}^{M}\|u_{h}(K_{i})-PP^{T}u_{h}(K_{i})\|_{2}^{2}.$$

By taking advantage of the singular value decomposition (SVD) of $\mathcal{U}$

$$U=P\Sigma Z^{T},$$

we can obtain the matrix $P_{N}$ whose columns are the $N$ orthogonal basis of $V_{N}$: that is the $N$ columns of $P$ corresponding to the largest diagonal elements in $\Sigma$.

We note that POD is a powerful tool to exploit the low-rank features of data. However, there are still two bottlenecks:

• •

  POD only provide the optimal linear representations of data in $l_{2}$ sense, and it is not guaranteed that they are the best basis functions in a broader sense.

• •

  RB methods, including POD-based RB methods, try to represent the solutions for different parameters in the same subspace of $V_{h}$, which means that a larger number of basis functions are needed to obtain accurate approximation of parametrized PDEs with high-dimensional parameter space. This will affect the efficiency and accuracy of RB methods.

Efficient implementation of RB-methods relies on the affine decompositions of stiffness matrix and load vectors, that is $A_{h}(u_{h};K)$ and $F_{h}(u_{h};K)$ can be written in parameter-seperable form, as follows

$$A_{h}(P_{N}^{T}u_{N};K)=\sum_{q=1}^{Q_{A}}c_{A}^{q}(u_{N},K)A_{h}^{q},\ F_{h}(P_{N}^{T}u_{N};K)=\sum_{q=1}^{Q_{F}}c_{F}^{q}(u_{N},K)F_{h}^{q}.$$

According to (3.2), we can obtain that

$$A_{N}(u_{N};K)=\sum_{q=1}^{Q_{A}}c_{A}^{q}(u_{N},K)A_{N}^{q},\ F_{N}(u_{N};K)=\sum_{q=1}^{Q_{F}}c_{F}^{q}(u_{N},K)F_{N}^{q},$$

where $A_{N}^{q}=P_{N}^{T}A_{h}^{q}P_{N},F_{N}^{q}=P_{N}^{T}F_{h}^{q}$. Precomputing $\{A_{N}^{q}\}_{q=1}^{Q_{A}}$ and $\{F_{N}^{q}\}_{q=1}^{Q_{F}}$ in the offline phase, we can avoid operations with a complexity dependent on $N_{h}$.

However, affine decompositions of bilinear forms (or linear forms) is a strong assumption, which is not satisfied in many cases. DEIM are used to provide an approximated parameter-seperable representations. During offline phase, POD is used to obtain parameter-independent functions $\{A_{N}^{q}\}_{q=1}^{Q_{A}}$ and $\{F_{N}^{q}\}_{q=1}^{Q_{F}}$. Then, during online phase, an interpolation problem is solved to obtain the coefficients $\{c_{N}^{q}(u_{N},K)\}_{q=1}^{Q_{A}}$ and $\{c_{N}^{q}(u_{N},K)\}_{q=1}^{Q_{F}}$ for each new parameters $K$.

DEIM techniques are powerful tool to deal with nonlinear and non-affine RB problems, however there are two main limitations:

• •

  For complex problems (i.e., high-dimensional and nonlinear problems), a large number of basis functions must be employed to obtain an accurate approximation, which prevents efficiency of online computation.

• •

  For nonlinear problems, DEIM can help to reduce the computation costs but can not avoid iteration (such as newton iteration methods) since it is an interpolation method.

In this subsection, we will introduce an activation function at the output layer to reconstruct the high-fidelity solution $u_{h}(K)$ from basis functions, that is, using $\hat{u}_{h}(K)=P_{N}(K)u_{N}(K)$ to approximate the high-fidelity solutions $u_{h}(K)$. Given $P_{N}(K)$ by Basis CNNs, we use Galerkin projection to obtain $u_{N}(K)$ by solving the algebra equation (3.1). Define the residual by

$$\mathcal{R}(u_{h},K)=A_{h}(K)u_{h}(K)-F_{h}(K).$$

Forcing the residual of $\hat{u}_{h}$ orthogonal to $V_{N}$, $u_{N}(K)$ satisfied

$$P_{N}(K)^{T}\mathcal{R}(\hat{u}_{h},K)=P_{N}(K)^{T}\Big{(}A_{h}(K)P_{N}(K)u_{N}(K)-F_{h}(K)\Big{)}=0.$$

We can then obtain the approximated solution $\hat{u}_{h}$ by solving above algebra equation, i.e.,

$$\hat{u}_{h}(K)=P_{N}(K)\Big{(}P_{N}(K)^{T}A_{h}(K)P_{N}(K)\Big{)}^{-1}P_{N}(K)^{T}F_{h}(K),$$

where $P_{N}(K)=\Phi_{\theta}^{Base}(K)$. We call it Galerkin-projection activation functions $\sigma_{G}$ and define it as a function of $P_{N}$, that is

$$\sigma_{G}(P_{N};A_{h},F_{h})=P_{N}(P_{N}^{T}A_{h}P_{N})^{-1}(P_{N}^{T}F_{h})=P_{N}A_{N}^{-1}F_{N},$$

where $A_{N}=P_{N}^{T}A_{h}P_{N},\ F_{N}=P_{N}^{T}F_{h}$.

Above all, We can write the forward propagation of Basis CNNs as

$$\hat{u}_{h}(K)=\sigma_{G}\circ\Phi_{\theta}^{Base}(K).$$

(4.2)

Given the dataset $\Big{\{}K_{i},u_{h}(K_{i})\Big{\}}_{i=1}^{M}$ and FEM bilinear forms (or linear forms) $\Big{\{}A_{h}(K_{i}),F_{h}(K_{i})\Big{\}}_{i=1}^{M}$, we define the loss function as follows

	$\displaystyle Loss_{RB}(\theta)$	$\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\|u_{h}(K_{i})-\hat{u}_{h}(K_{i})\|_{L_{2}}^{2}+\lambda_{G}\ cond_{F}\Big{(}A_{N}(K_{i})\Big{)}^{2},$		(4.3)
		$\displaystyle=\frac{1}{M}\sum_{i=1}^{M}\|u_{h}(K_{i})-\sigma_{G}\big{(}\Phi_{\theta}^{Base}(K_{i});A_{h}(K_{i}),F_{h}(K_{i})\big{)}\|_{L_{2}}^{2}$
		$\displaystyle+\lambda_{G}\ cond_{F}\Big{(}\big{(}\Phi_{\theta}^{Base}(K_{i})\big{)}^{T}A_{h}\Phi_{\theta}^{Base}(K_{i})\Big{)}^{2},$

where the $L_{2}$ is the $L_{2}$ norm. In practical, we approximate the $L_{2}$ norm in the loss function at the discrete points on the FEM mesh $\mathcal{M}$, and $cond_{F}$ is the condition number under Frobenious norm, i.e.,

$$cond_{F}(A)=\|A\|_{F}\|A\|_{F}^{-1}.$$

The first term of loss function (4.3) trains the neural networks to approximate the data. As shown in Remark 4.3, the training depends on the inverse of $A_{N}$. However, $A_{N}$ learned by neural networks may be ill-conditioned because it does not have good properties, such as orthogonality [33]. To guarantee the stability of training and the reduced-order model (3.1), we impose the second term in the loss function. In theoretical analysis, condition number under 2-norm are widely used, but it is intractable to compute during training. Therefore, we choose the condition number under Frobenious norm here. The stability parameter $\lambda_{G}$ is a hyperparameter.

In the fist test case, we consider the binomial point process $\kappa(x,\xi)$ defined on $\mathcal{P}=(\Omega,\mathcal{F},\mathbb{P})$ for any given $x\in\mathcal{S}$. Let any subset of $\mathcal{S}$ be $\mathcal{B}$ and $\{X_{i}\}_{i=1}^{n}$ be n $i.i.d$ points with the distribution $\mathbb{P}(\mathcal{S})$ and the density $p(x)$. Then the number of points in $\mathcal{\mathcal{B}}$ can be defined as the following random variable

$$\xi(\mathcal{B})=\sum_{i=1}^{n}\delta(X_{i}\in\mathcal{B}),$$

where $\delta(X_{i}\in\mathcal{B})$ is the indicator function and takes value as

$$\delta(X_{i}\in\mathcal{B})=\left\{\begin{aligned} &1\ \ \ \mathrm{if}\ X_{i}\in\mathcal{B},\\ &0\ \ \ \mathrm{otherwise}.\end{aligned}\right.$$

The probability of each random variables $X_{i}$ in $\mathcal{B}$ is the same and can be defined as follows,

$$p=\mathbb{P}(X_{i}\in\mathcal{B})=\frac{\int_{\mathcal{B}}p(x)dx}{\int_{\mathcal{S}}p(x)dx}.$$

Then $\xi(\mathcal{B})$ is distributed to a binomial distribution with $n$ and $p$. We introduce the set

$$F(\xi)=\cup_{i=1}^{n}\{x\in\mathbb{R}^{2}|\|x-x_{i}\|_{2}\leq r\},$$

where $x_{i}$ is the sample of $X_{i}\sim p(\mathcal{S})$ and $r$ is the radius of points. Then we can define the two-scale random permeability $\kappa(x,\xi)$, where points of permeability $\kappa^{(1)}$ are randomly distributed in the background media of conductivity $\kappa^{(0)}$. Specifically,

$$\kappa(x,\xi)=\left\{\begin{aligned} &\kappa^{(1)}\ \ \ \mathrm{if}\ x\in F(\xi),\\ &\kappa^{(0)}\ \ \ \mathrm{otherwise}.\end{aligned}\right.$$

Here, we consider the special case that $\mathcal{B}=\mathcal{S}$, $\mathbb{P}$ is the uniform distribution on $\mathcal{S}$ and takes $n=5,r=0.05,\kappa^{(1)}=1000,\kappa^{(0)}=1$. In the following experiments, we discretize $\kappa(x,\xi)$ over the uniform FEM mesh $\mathcal{M}$ with rectangular elements and $61\times 61$ nodes, denoted by $K$. Several instances of $K$ are given in Figure 5.2.

In this paper, we only approximate the solutions on the free nodes, i.e., the spatial resolution of inputs is $59\times 59$. The configuration of Basis CNNs are shown in Figure 4.2(b) with more details in Table 1. The 2nd-3th columns of Table 1 show the spatial resolution $H\times W$ of feature maps and the number of parameters of each layer. This configuration will be used in all experiments except that the spatial resolution of feature maps varies with different cases.

The Basis CNNs are trained with Adam optimizer with loss (4.3). The initial learning rate is 0.0001 and cosine decay, which is a built-in function of TensorFlow, is used. The batch size is set as 32 and the model is trained for 100 epochs. Each sample in the dataset is composed of four variables, they are permeability $K$, high-fidelity solution $u_{h}(K)$, stiffness matrix $A_{h}(K)$ and load vector $F_{h}(K)$, respectively (see equation (5.3)-(5.5)). We collect the data pairs for 12740 different samples of $K$, 10240 of which are selected as training data and the rest are used to test the generalization ability of CNN-based ROM.

The RB methods exploit low-rank features of the data, which is usually called modes, that characterize physical processes. Therefore, we show the RB basis functions learned by CNN-based ROM and POD-based RB methods in Figure 5.4. Note that for POD-based RB basis functions, the modes are distributed around the center and dominant modes alternate arrangement. In contrast, the shape of the CNN-based RB basis functions are very much like the data. Such input-different basis functions can help reduce the number of basis functions needed for approximation. However, they lack interpretability, which is the common problem of deep learning techniques.

Our training data are generated by FEM, which has approximation error. For instance, FEM solutions on a mesh with $31\times 31$ nodes are less accurate than those on a mesh with $61\times 61$ (see Figure 5.5). We use a Gaussian distribution with means of zero and standard deviations of 0.001 to simulate the errors caused by numerical methods, which is shown in the second column of Figure 5.6.

Using the noisy data, we train the Basis CNNs with $N=10$ and implement POD-based RB methods with $N=100,1000$ again. The numerical results are presented in Figure 5.6. CNN-based ROM shows a clear advantage and a de-noising effect.

In this subsection, we consider the random input in Figure 5.7, which simulates the hydraulic conductivity field of a channelized aquifer [36]. The total number of pixels is $2500\times 2500$. Channel and matrix materials are assigned hydraulic conductivity values of 1000 and 1, respectively. To introduce randomness, we extract 23104 small images of resolution $64\times 64$ from the large image by moving along both the horizontal and vertical directions with an interval of 16 pixels for each movement. To obtain enough training data, we flip the generated images horizontally (see Figure 5.7). We finally get 46208 images, among which 40960 images are collected as training dataset and the rest are used for test. Then Basis CNNs with architecture (Table 1) are trained using Adam optimizer for 50 epochs. The relative tests mean error is $2.16\%$ and the generalization ability of the proposed method is intuitively presented in Figure 5.8.

Next, we use the upper-left quarter of this background field to validate the ideas in Section 4.3, i.e., the total number of pixels is $1250\times 1250$. Here, we consider an extreme case that the oversampling domain covers the whole computational region. Linear boundary conditions defined in equation (4.5) are used for the oversampling domain. Four nodal basis functions are obtained with different boundary conditions. Rather than training four different CNNs, we fully utilize the rotational symmetry of the boundary conditions, i.e., each boundary condition can be rotated clockwise by 90 degrees to overlap with another boundary condition. This can be accomplished by rotating background field. More specifically, we use Basis CNNs to learn the following equation,

$$\left\{\begin{aligned} &-\nabla\cdot\Big{(}\kappa(x,\xi)\nabla\tilde{\Phi}_{ms}(x)\Big{)}=0,\ x\in\mathcal{S}:=[0,1]\times[0,1],\\ &\tilde{\Phi}_{ms}(0,1)=1,\tilde{\Phi}(0,0)=0,\tilde{\Phi}_{ms}(1,0)=0,\tilde{\Phi}_{ms}(1,1)=0.\end{aligned}\right.$$

To generate training and test dataset, we first obtain 5625 sub-image of resolution $65\times 65$ by moving along both the horizontal and vertical directions with an interval of 16 pixels for each movement. Then each sample are rotated clockwise by 90 degrees three times to acquire three additional samples (see Figure 5.9).

The relative test mean error is $1.38\%$ and several test instances are shown in Figure 5.10. Note that the basis functions learned by Basis CNNs are temporary, we construct the true multiscale finite element basis functions $\Phi_{ms}^{i},i=1,2,3,4$ in each coarse elements $\omega$ from the linear combination of $\tilde{\Phi}_{ms}^{i},i=1,2,3,4$, i.e.,

$$\Phi_{ms}^{i}=\sum_{k=1}^{4}c_{ik}\tilde{\Phi}_{ms}^{k},$$

where the coefficients $c_{ik}$ is determined by conditions $\Phi_{ms}^{i}(x_{j})=\delta_{ij}$. This leads to $N$ non-conforming basis functions which is dealt with by simply taking averages on the boundary because there only exists an $O(1)$ jump across $\partial{\omega}$ [5]. The numerical results depicted in Figure 5.10 demonstrate that the same basis functions can be used for different types of source terms and the numerical solution is convergent to the reference as the number of basis functions is large enough.

We first consider the case that the source term is nonlinear as follows,

$$\left\{\begin{aligned} &-\nabla\cdot\Big{(}\kappa(x,\xi)\nabla u(x)\Big{)}=\sin\Big{(}10\pi u(x)\Big{)}+\cos\Big{(}10\pi u(x)\Big{)},x\in\mathcal{S},\\ &u(x)=0,x\in\partial{\mathcal{S}}.\end{aligned}\right.$$

(5.9)

The weak form of above equation is that find $u\in V$, such that

$$\int_{\mathcal{S}}\kappa\nabla u\nabla v\ dx=\int_{\mathcal{S}}\Big{(}\sin(10\pi u)+\cos(10\pi u)\Big{)}\ v\ dx,\ \forall v\in V.$$

(5.10)

We solve above variational problems by Newton’s iteration methods on a fixed mesh with rectangular elements and $65\times 65$ nodes. Once we get the high-fidelity solution $u_{h}(K)$, we can obtain entries of stiffness matrix $A_{h}(K)\subset\mathbb{R}^{4225\times 4225}$ and the load vector $F_{h}(K)\subset\mathbb{R}^{4225\times 1}$ are

	$\displaystyle a_{h}(i,j)$	$\displaystyle=\int_{\mathcal{S}}K\nabla\varphi_{i}\nabla\varphi_{j}\ dx,$
	$\displaystyle f_{h}(i)$	$\displaystyle=\int_{\mathcal{S}}u_{h}(K)\nabla\varphi_{i}\ dx.$

where $\{\varphi_{i},\ldots,\varphi_{4225}\}$ is the bilinear basis functions of $V_{h}$. Here, we use element means to approximate the value of $u_{h}$ at the integral points (the same as Figure 5.1). The Basis CNNs with architecture (Table 1) are first trained using 5120 data pairs. Since $u_{h}(K)$ is unknown for test instances, Basis CNNs themselves can not realize online computation. We further use Coef CNNs to learn the mapping from $K$ to the coefficients $c(K)$. The specific architecture of Coef CNNs is presented in Table 3

The same dataset is used to train Coef RNN. The relative test mean error is $4.43\%$, which is estimated by 2500 test instances. Besides, several results are shown in Figure 5.12. Compared to traditional ROM, which need iterative methods on the online stage, we can directly obtain the prediction of solutions using CNN-based ROM.

In this section, we consider the p-Laplacian equations as follows,

$$\left\{\begin{aligned} &-\nabla\cdot\Big{(}\kappa(x,\xi)|u(x)|^{p-2}\nabla u(x)\Big{)}=1,x\in\mathcal{S},\\ &u(x)=0,x\in\partial{\mathcal{S}}.\end{aligned}\right.$$

(5.11)

The p-Laplacian operator in the left hand of equation (5.11) is derived from a nonlinear Darcy law and continuity equation. For p=2, equation (5.11) is equivalent to the equation in Section 5.1, which is linear. Here, we take $p=3$, and minimize the total potential energy of the flow, i.e.,

$$u=\mathop{\arg\min}\limits_{v\in V}\frac{1}{p}\int_{\mathcal{S}}\kappa\|\nabla v\|_{2}^{p}-\int_{\mathcal{S}}vdx.$$

We approximate $u$ on the mesh with triangular elements and $65\times 65$ nodes. The above optimization problems are solved for 12740 samples of the log-Gaussian random field with bandwidth 0.1 by an efficient algorithm [37], among which 10240 instances are split to training dataset and the rest are used for test. The stiffness matrix $A_{h}(K)\in\mathbb{R}^{4225\times 4225}$ is computed by assuming that $u(K)$ is accurately approximated by $u_{h}(K)$, which has entries

$$a_{h}(i,j)=\int_{\mathcal{S}}K|u_{h}(K)|^{p-2}\nabla\varphi_{i}\nabla\varphi_{j}\ dx,$$

where $\{\varphi_{i},\ldots,\varphi_{4225}\}$ is the piecewise linear basis functions of $V_{h}$. The definition of load vector $F_{h}(K)\in\mathbb{R}^{4225\times 1}$ is the same as equation (5.4). We train both Basis CNNs and Coef CNNs by Adam for 50 epochs. The relative test mean error is $1.49\%$ and several test instances are shown in Figure 5.13.