Learning Adaptive Coarse Spaces of BDDC Algorithms for Stochastic Elliptic Problems with Oscillatory and High Contrast Coefficients

Abstract

In this paper, we consider the balancing domain decomposition by constraints (BDDC) algorithm with adaptive coarse spaces for a class of stochastic elliptic problems. The key ingredient in the construction of the coarse space is the solutions of local spectral problems, which depend on the coefficient of the PDE. This poses a significant challenge for stochastic coefficients as it is computationally expensive to solve the local spectral problems for every realization of the coefficient. To tackle this computational burden, we propose a machine learning approach. Our method is based on the use of a deep neural network (DNN) to approximate the relation between the stochastic coefficients and the coarse spaces. For the input of the DNN, we apply the Karhunen–Loève expansion and use the first few dominant terms in the expansion. The output of the DNN is the resulting coarse space, which is then applied with the standard adaptive BDDC algorithm. We will present some numerical results with oscillatory and high contrast coefficients to show the efficiency and robustness of the proposed scheme.

1. Introduction

Coefficient functions are often one of the main difficulties in modeling a real life problem. Due to the nature of the phenomena, these coefficients can be oscillatory and with high contrast features, all such conditions hinder the development of a fast and accurate numerical solution. Under some conditions, numerical methods like discontinuous Galerkin methods could give a robust scheme, as in [1,2]. However, more importantly, coefficients are not constant and they usually involve uncertainties and randomness. Therefore, a stochastic partial differential equation (SPDE) is considered for the modeling of these stochastic coefficient functions. There are several popular approaches for solving SPDE numerically, for example, the Monte Carlo simulation, stochastic Galerkin method, and the stochastic collocation method in [3,4,5]. All these methods require a high computational power in realistic simulations and cannot be generalized to similar situations. In the past few decades, machine learning techniques have been applied on many aspects, such as image classification, fluid dynamics, and solving differential equations in [6,7,8,9,10,11,12,13,14]. In this paper, instead of directly solving the SPDE by machine learning, we choose a moderate step that involves the advantage of machine learning and, also, the accuracy from the adaptive BDDC algorithm.

The combination of machine learning and domain decomposition methods has been discussed in [15,16,17] for different applications. A unified review on current approaches is presented in [10], which introduces three general approaches on this combination method: the first one uses machine learning techniques to improve the convergence properties and the computational efficiency of domain decomposition methods. For the second type of approach, DNN is used as discretization methods and solvers to solve the differential equations directly. The last one makes use of the properties and ideas of domain decomposition methods to improve the performance of neural network.

Although close to the first approach descried, our idea in this paper is not actually the same. We first apply the Karhunen–Loève expansion on the stochastic coefficient as in [18,19,20] to capture the characteristic with finite terms. After this expansion, the stochastic coefficient function of SPDE is reduced to a finite expansion with desired accuracy and can be solved with a balancing domain decomposition by constraints (BDDC) algorithm, which can prevent the computational problem caused by the refinement of spatial mesh resolution and give an accurate solution rather than directly using machine learning technique without suitably designed conditions. After obtaining sufficient number of random solutions from the BDDC method, we can derive the statistics of the original SPDE solution. However, the cost for forming and solving generalized eigenvalue problems in the BDDC algorithm with good accuracy is considerable especially for high dimensional problems. We note that such eigenvalue problems are considered to enrich the coarse space adaptivity in the BDDC algorithm. Therefore, we introduce neural network with a user-defined number of layers and neurons to compensate the increase in computational cost with unexpected times of BDDC solutions.

The adaptive BDDC algorithm with enriched primal unknowns is considered in this paper as its ability on oscillatory and high contrast coefficients has been shown in [21,22]. Among different domain decomposition methods, the considered adaptive BDDC does not require a strong assumption on coefficient functions and the subdomain partitions to achieve a good performance like the standard BDDC algorithm [23,24]. It is because the additional coarse space basis functions computed by the dominant eigenfunctions are related to the ill-conditioning in the standard BDDC algorithm with highly varying coefficients. The introduction of these dominant eigenfunctions, thus, could greatly improve the robustness of numerical scheme for problems with rough and high contrast coefficients. Moreover, an estimate of condition number of this adaptive BDDC algorithm can be shown to be only controlled by a given tolerance without any extra assumptions on the coefficients or subdomain partitions. For simplicity, in the remaining of the paper, we will call the resulting new numerical scheme combining the machine learning technique and adaptive BDDC algorithm as learning adaptive BDDC algorithm.

Finally, to end this section, we state the main idea of our proposed scheme and the model problem. In this paper, we integrate an artificial neural network with learning abilities into a BDDC algorithm with adaptively enriched coarse spaces for efficient solutions of the neural network approximation of stochastic elliptic problems. Let $\mathcal{D}$ be a bounded domain in ${\mathbb{R}}^d$, $d=2$, and $\mathsf{\Omega }$ be a stochastic probability space:

$$\begin{array}{cc}\hfill -\nabla ·\left(\rho \right(\mathit{x},\omega )\nabla u)& =f,\mathrm{in}\mathcal{D},\hfill \\ \hfill u& =0,\mathrm{on}\partial \mathcal{D},\hfill \end{array}$$

(1)

where $\mathit{x}\in \mathcal{D}$, $\omega \in \mathsf{\Omega }$, and $\rho (\mathit{x},\omega )$ is uniformly positive and is highly heterogeneous with very high contrast. Here, we only assume the two-dimensional spatial domain, however, the three-dimensional case is also applicable and is shown to be robust for deterministic elliptic problems in [21]. Moreover, as the model problem (1) can be extended to fluid flow problems, we seldom call the coefficient function $\rho (\mathit{x},\omega )$ as permeability function.

The rest of the paper is organized as follows. In Section 2.1, a brief formulation of the BDDC algorithm with adaptively enriched coarse problems is presented for two-dimensional elliptic problems. Then, the Karhunen–Loève expansion is introduced and the details of the artificial neural network used are clarified in Section 2.2. Two concrete examples with an explicit analytical expression of the Karhunen–Loève expansion, some network training, and testing parameters, are presented in the first half of Section 3 and the results of various numerical experiments are presented subsequently. Finally, a concluding remark is given in Section 4.

2. Materials and Methods

2.1. Adaptive BDDC Algorithm

The main feature of the adaptive BDDC is the local generalized eigenvalue problem defined on every subdomain interface, which introduces the adaptively enriched primal unknowns. In this subsection, we repeat the formulation presented by the first and second authors in [21] and give a brief overview of the adaptive BDDC algorithm. We refer to [21,22] for further details, where the analysis of the condition number estimation and the robustness of numerical scheme with oscillatory and high contrast coefficients can be found.

2.1.1. Local Linear System

We first introduce a discrete form of the model problem (1) in a deterministic fashion. Let $V_h$ be the space of conforming linear finite element functions with respect to a given mesh on $\mathcal{D}$ with the mesh size $h>0$ and with the zero value on $\partial \mathcal{D}$. We will then find the approximate solution $u\in V_h$ such that

$$a(u,v)=(f,v),\forall v\in V_h,$$

(2)

where

$$a(u,v)={\int }_{\mathcal{D}}\rho \left(\mathit{x}\right)\nabla u·\nabla vdx,(f,v)={\int }_{\mathcal{D}}fvdx.$$

(3)

We assume that the spatial domain $\mathcal{D}$ is partitioned into a set of N non-overlapping subdomains $\left\{{\mathcal{D}}_i\right\}$, $i=1,2,\cdots ,N$, so that $\mathcal{D}={\cup }_{i=1}^N{\mathcal{D}}_i$. We note that the subdomain boundaries do not cut triangles equipped for $V_h$. We allow the coefficient $\rho \left(\mathit{x}\right)$ to have high contrast jumps and oscillations across subdomains and on subdomain interfaces. Let $a_i(u,v)$ be the bilinear form of the model elliptic problem (2) restricted to each subdomain ${\mathcal{D}}_i$ defined as

$$a_i(u,v)={\int }_{{\mathcal{D}}_i}\rho \left(\mathit{x}\right)\nabla u·\nabla vdx,\forall u,v\in X_i,$$

where $X_i$ is the restriction of $V_h$ to ${\mathcal{D}}_i$.

In the BDDC algorithm, the original problem (2) is reduced to a subdomain interface problem and solved by an iterative method combined with a pre-conditioner. The interface problem can be obtained by solving a Dirichlet problem in each subdomain. After choosing dual and primal unknowns on the subdomain interface unknowns, the interface problem is then solved by utilizing local problems and one global coarse problem corresponding to the chosen sets of dual and primal unknowns, respectively. At each iteration, the residuals are multiplied by certain scaling factors to balance the errors across the subdomain interface regarding to the energy of each subdomain problem. The coarse problem aims to correct the global part of the error in each iteration and, thus, the choice of primal unknowns is important in obtaining a good performance as the number of subdomains increases. The basis for primal unknowns is obtained by the minimum energy extension for a given constraint at the location of primal unknowns and such a basis provides a robust coarse problem with a good energy estimate. We refer to [23,24,25,26] for general introductions to BDDC algorithm.

2.1.2. Notation and Preliminary Results

To facilitate our discussion, we first introduce some notation. Let $S^{\left(i\right)}$ be the Schur complement matrix obtained from the local stiffness matrix $A^{\left(i\right)}$ after eliminating unknowns interior to ${\mathcal{D}}_i$, where $A^{\left(i\right)}$ is defined by $a_i(u,v)=u^T{A}^{\left(i\right)}v$, for all $u,v\in X_i$. In the following, we will use the same symbol to represent a finite element function and its corresponding coefficient vector in order to simplify the notations.

Recall that $X_i$ is the restriction of the finite element space $V_h$ to each subdomain ${\mathcal{D}}_i$. Let $W_i$ be the restriction of $X_i$ to $\partial {\mathcal{D}}_i$. We then introduce the product spaces

$$X=\prod _{i=1}^N{X}_i,W=\prod _{i=1}^N{W}_i,$$

where we remark that the functions in X and W are totally decoupled across the subdomain interfaces. In addition, we introduce partially coupled subspaces $\tilde{X}$, $\tilde{W}$, and fully coupled subspaces $\widehat{X}$, $\widehat{W}$, where some primal unknowns are strongly coupled for functions in $\tilde{X}$ or $\tilde{W}$, while the functions in $\widehat{X}$, $\widehat{W}$ are fully coupled across the subdomain interfaces.

Next, we present the basic description of the BDDC algorithm; see [23,24,25,26]. For simplicity, the two-dimensional case will be considered. After eliminating unknowns interior to each subdomain, the Schur complement matrices $S^{\left(i\right)}$ are obtained from $A^{\left(i\right)}$, and $g_i$ is denoted as the corresponding right-hand side unknowns on the subdomain boundary. Summing up across the subdomains, the problem (2) is reduced in the BDDC algorithm, which is to find $\widehat{w}\in \widehat{W}$, such that

$$\sum _{i=1}^N{R}_i^T{S}^{\left(i\right)}{R}_i\widehat{w}=\sum _{i=1}^N{R}_i^T{g}_i,$$

(4)

where $R_i:\widehat{W}\to W_i$ is the restriction operator into $\partial {\mathcal{D}}_i$, and $g_i\in W_i$ depends on the source term f. Note, the solution $\widehat{w}$ to problem (4) is actually the restriction of u in (2) to the subdomain boundaries.

The BDDC pre-conditioner is built based on the partially coupled space $\tilde{W}$. Let ${\tilde{R}}_i:\tilde{W}\to W_i$ be the restriction into $\partial {\mathcal{D}}_i$, and let $\tilde{S}$ be the partially coupled matrix defined by

$$\tilde{S}=\sum _{i=1}^N{\tilde{R}}_i^T{S}^{\left(i\right)}{\tilde{R}}_i.$$

For the space $\tilde{W}$, we can express it as the product of the two spaces

$$\tilde{W}=W_{\Delta }\times {\widehat{W}}_{\Pi },$$

where ${\widehat{W}}_{\Pi }$ consists of vectors of the primal unknowns and $W_{\Delta }$ consists of vectors of dual unknowns, which are strongly coupled at the primal unknowns and decoupled at the remaining interface unknowns, respectively. We define $\widehat{R}:\widehat{W}\to \tilde{W}$ such that

$$\widehat{R}=\left(\begin{array}{c}{R}_{\Delta }\\ R_{\Pi }\end{array}\right),$$

where $R_{\Delta }$ is the mapping from $\widehat{W}$ to $W_{\Delta }$ and $R_{\Pi }$ is the restriction from $\widehat{W}$ to ${\widehat{W}}_{\Pi }$. We note that $R_{\Delta }$ is obtained as

$$R_{\Delta }=\left(\begin{array}{c}{R}_{\Delta }^{\left(1\right)}\\ R_{\Delta }^{\left(2\right)}\\ ⋮\\ R_{\Delta }^{\left(N\right)}\end{array}\right),$$

where $R_{\Delta }^{\left(i\right)}$ is the restriction from $\widehat{W}$ to $W_{\Delta }^{\left(i\right)}$ and $W_{\Delta }^{\left(i\right)}$ is the space of dual unknowns of ${\mathcal{D}}_i$.

The BDDC pre-conditioner is then given by

$$M_{BDDC}^{-1}={\widehat{R}}^T\tilde{D}{\tilde{S}}^{-1}{\tilde{D}}^T\widehat{R}$$

(5)

where $\tilde{D}$ is a scaling matrix of the form

$$\tilde{D}=\sum _{i=1}^N{\tilde{R}}_i^T{D}_i{\tilde{R}}_i.$$

Here, the matrices $D_i$ are defined for unknowns in $W_i$ and they are introduced to make the pre-conditioner robust to the heterogeneity in $\rho \left(\mathit{x}\right)$ across the subdomain interface. In more detail, $D_i$ consists of blocks $D_F^{\left(i\right)}$ and $D_V^{\left(i\right)}$, where F denotes an equivalence class shared by two subdomains, i.e., ${\mathcal{D}}_i$ and its neighboring subdomain ${\mathcal{D}}_j$, and V denotes the end points of F, respectively. We call such equivalence classes F and V as edge and vertex in two dimensions, respectively. In our BDDC algorithm, unknowns at subdomain vertices are included to the set of primal unknowns and adaptively selected primal constraints are later included to the set after a change of basis formulation. For a given edge F in two dimensions, the matrices $D_F^{\left(l\right)}$ and $D_V^{\left(l\right)}$ satisfy a partition of unity property, i.e., $D_F^{\left(i\right)}+D_F^{\left(j\right)}=I$ and ${\sum }_{l\in n\left(V\right)}{D}_V^{\left(l\right)}=1$, where $n\left(V\right)$ denotes the set of subdomain indices sharing the vertex V. The matrices $D_F^{\left(l\right)}$ and $D_V^{\left(l\right)}$ are called scaling matrices. As mentioned earlier, the scaling matrices help to balance the residual error at each iteration with respect to the energy of subdomain problems sharing the interface. For the case when $\rho \left(\mathit{x}\right)$ is identical across the interface F, $D_F^{\left(i\right)}$ and $D_F^{\left(j\right)}$ are chosen simply as multiplicity scalings, i.e., $1/2$, but for a general case when $\rho \left(\mathit{x}\right)$ has discontinuities, different choice of scalings, such as $\rho$-scalings or deluxe scalings, can be more effective. The scaling matrices $D_V^{\left(l\right)}$ can be chosen using similar ideas. We refer to [27] and references therein for scaling matrices.

We note that by using the definitions of $\widehat{R}$ and $\tilde{S}$, the matrix in the left hand side of (4) can be written as

$$\sum _{i=1}^N{R}_i^T{S}^{\left(i\right)}{R}_i={\widehat{R}}^T\tilde{S}\widehat{R}.$$

In the BDDC algorithm, the system in (4) is solved by an iterative method with the preconditioner (5). Thus, its performance is analyzed by estimating the condition number of

$$M_{BDDC}^{-1}{\widehat{R}}^T\tilde{S}\widehat{R}={\widehat{R}}^T\tilde{D}{\tilde{S}}^{-1}{\tilde{D}}^T\widehat{R}{\widehat{R}}^T\tilde{S}\widehat{R}.$$

(6)

2.1.3. Generalized Eigenvalue Problems

After the preliminaries, we are ready to state the generalized eigenvalue problems which introduce the adaptive enrichment of coarse components. For an equivalence class F shared by two subdomains ${\mathcal{D}}_i$ and ${\mathcal{D}}_j$, the following generalized eigenvalue problem is proposed in [28]:

$$\left({\left(D_F^{\left(j\right)}\right)}^T{S}_F^{\left(i\right)}{D}_F^{\left(j\right)}+{\left(D_F^{\left(i\right)}\right)}^T{S}_F^{\left(j\right)}{D}_F^{\left(i\right)}\right)v=\lambda \left({\tilde{S}}_F^{\left(i\right)}:{\tilde{S}}_F^{\left(j\right)}\right)v,$$

(7)

where $S_F^{\left(i\right)}$ and $D_F^{\left(i\right)}$ are the block matrices of $S^{\left(i\right)}$ and $D_i$ corresponding to unknowns interior to F, respectively. The matrix ${\tilde{S}}_F^{\left(i\right)}$ are the Schur complement of $S^{\left(i\right)}$ obtained after eliminating unknowns except those interior to F. In addition, for symmetric and semi-positive definite matrices A and B, their parallel sum $A:B$ is defined by, see [29],

$$A:B=B{(A+B)}^{+}A,$$

(8)

where ${(A+B)}^{+}$ is a pseudo inverse of $A+B$. We note that the problem in (7) is identical to that considered in [30] when $D_F^{\left(i\right)}$ are chosen as the deluxe scalings, i.e.,

$$D_F^{\left(i\right)}={(S_F^{\left(i\right)}+S_F^{\left(j\right)})}^{-1}{S}_F^{\left(i\right)}.$$

We note that a similar generalized eigenvalue problem is considered and extended to three-dimensional problems in the first and second authors’ work [21].

2.2. Learning Adaptive BDDC Algorithm

In addition to the BDDC algorithm, we perform the Karhunen–Loève (KL) expansion to decompose the stochastic permeability $\rho (\mathit{x},\omega )$ for the preparation of the learning adaptive BDDC algorithm. As shown by Theorem 2.7, Proposition 2.18 in [31] and the numerical examples in [20], a small number of terms in the KL expansion can approximate the stochastic permeability with a reasonable accuracy. Thus, corresponding computational cost can also be reduced efficiently. Once we obtain the KL series expansion, we can easily produce a certain amount of sample data as a training set for the artificial neural network. After the training process, the neural network learns the main characteristics of the dominant eigenvectors from the adaptive BDDC algorithm, and neural network approximation on coarse spaces can be obtained. In our proposed algorithm, we still need some resulting samples from the BDDC algorithm to train the neural network, which is a process of a relatively high computational cost. However, after this setup, we can apply this network to obtain an approximate solution directly without the use of other stochastic sampling methods, such as the Monte Carlo simulation, which involves even higher computational cost for the generation of large number of random samples or realizations.

In this subsection, we discuss the proposed algorithm in details. Moreover, our resulting neural network can be also applied on problems with similar stochastic properties. It thus saves the cost of training again and we can readily use the trained neural network for prediction. Some testing samples will be generated to test the network performance and to show its generalization capacity in the later section.

2.2.1. Karhunen–Loève Expansion

To start our learning adaptive algorithm, we first apply the KL expansion on the stochastic field $\rho (\mathit{x},\omega )$. To ensure the positivity of stochastic permeability field, a logarithmic transformation is considered as $K(\mathit{x},\omega )=log\left(\rho \right(\mathit{x},\omega \left)\right)$. Let $C_K(\mathit{x},\widehat{\mathit{x}})$ be the covariance function of $K(\mathit{x},\omega )$ at two locations $\mathit{x}$ and $\widehat{\mathit{x}}$. Since the covariance function $C_K(\mathit{x},\widehat{\mathit{x}})$ is symmetric and positive definite, it can be decomposed into:

$$\begin{array}{c}\hfill C_K(\mathit{x},\widehat{\mathit{x}})=\sum _{i=1}^{\infty }{\lambda }_i{f}_i\left(\mathit{x}\right)f_i\left(\widehat{\mathit{x}}\right),\end{array}$$

(9)

where ${\lambda }_i$ and $f_i$ are eigenvalues and eigenfunctions computed from the following Fredholm integral equation:

$$\begin{array}{c}\hfill {\int }_{\mathcal{D}}{C}_K(\mathit{x},\widehat{\mathit{x}})f\left(\mathit{x}\right)d\mathit{x}=\lambda f\left(\widehat{\mathit{x}}\right),\end{array}$$

(10)

where $\left\{f_i\left(\mathit{x}\right)\right\}$ are orthogonal and deterministic functions. It is noted that for some specially chosen covariance functions, we can find the $(\lambda ,f)$ eigenpair analytically. We will list some examples in the section of numerical results. After solving the Fredholm integral equation, we are ready to express the KL expansion of the log permeability:

$$\begin{array}{c}\hfill K(\mathit{x},\omega )=E\left[K\right]\left(\mathit{x}\right)+\sum _{i=1}^{\infty }\sqrt{{\lambda }_i}{\xi }_i{f}_i\left(\mathit{x}\right).\end{array}$$

(11)

Here $E\left[K\right]\left(\mathit{x}\right)$ is the expected value of $K(\mathit{x},\omega )$ and ${\xi }_i$ are independent and identically distributed Gaussian random variables with mean 0 and variance 1. As mentioned above, a few terms in KL expansion can give a reasonably accurate approximation, which can be explained now as we can just pick the dominant eigenfunctions in (11) and the other eigenvalues will decay rapidly.

Therefore, KL expansion is an efficient method to represent the stochastic coefficient and favors the subsequent process. Although analytical eigenfunctions and eigenvalues cannot always be found, there are various numerical methods to compute the KL coefficients in [31,32]. After the generation of some sets of $\left\{{\xi }_i\right\}$, which are the only stochastic variables in the permeability after KL expansion, as the input of neural network, we plug back KL expanded permeability coefficients into the adaptive BDDC algorithm, and the resulting dominant eigenfunctions in the coarse spaces will be the target output of the our designated neural network.

2.2.2. Neural Network

There are many types of neural networks for different usages. Since in our learning adaptive BDDC algorithm the desired neural network is to capture numeric feature from the data as a supervised learning, a fully connected feed-forward neural network is chosen. An illustration of the neural network structure is shown in Figure 1. In a general feed-forward neural network, there are L hidden layers between the input and output layer, and the arrows from layers to layers indicate that all the neurons in the previous layer are used to produce every neuron in the next layer. Computational cost increases with L, but significant improvements in predictions of our tests are not always resulted, $L=1$ is thus chosen for most of the numerical experiments.

We denote the number of neurons in the input layer to be R, which is the number of terms in the truncated KL expansion. The inputs are actually the Gaussian random variables $\left\{{\xi }_i\right\}$ in one realization. The output of the network is a column vector that consists of all dominant eigenvectors in the coarse space. The number of neurons in the output layer is denoted as O, which is the sum of length of all resulting eigenvectors obtained from the BDDC algorithm. This value O depends on many factors, such as the geometry of the spatial domain, the structure of the grid, and parameters of the BDDC algorithm.

In the training, we use the scaled conjugate gradient algorithm in [33] to do the task of minimization. Therefore, unlike the gradient descent method, we do not need to determine the value of learning rate at each step. However, the initial choice for the scaled conjugate gradient algorithm may still affect the performance of training. Before training, we have to generate training samples and decide some conditions on stopping criteria. When the minimization process is completed for all the training samples, we are said to finish one epoch. There are two stopping criteria that are usually considered. They are the tolerance of the cost function gradient and the maximum number of epochs trained. When either of these conditions are satisfied, that is, the cost function gradient is below the specified tolerance or the number of epochs reaches the maximum number, the training is stopped and the network will be tested for performance on a testing set.

In both training and testing results, in order to obtain an error measure less affected by scales of data set, a normalized root mean squared error (NRMSE) is used as an error of reference to see how well are the network estimation and prediction:

$$\begin{array}{cc}\hfill & RMSE=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^M\frac{1}{O}\left|\right|y\left({\mathit{\xi }}^{\left(i\right)}\right)-f_{NN}\left({\mathit{\xi }}^{\left(i\right)}\right){\left|\right|}_2^2}{M}}},\hfill \\ \hfill & NRMSE={\displaystyle \frac{RMSE}{{max}_{i=1,2,\cdots ,M}\left|\right|y\left({\mathit{\xi }}^{\left(i\right)}\right)-f_{NN}\left({\mathit{\xi }}^{\left(i\right)}\right){\left|\right|}_{\infty }}},\hfill \end{array}$$

where ${\mathit{\xi }}^{\left(i\right)}\sim N(0,1)\in {\mathbb{R}}^R$ is a column vector of random variables that follow the standard normal distribution in the i-th realization, $y(·)$ is the dominant eigenvectors obtained from the adaptive BDDC algorithm when the KL expanded stochastic permeability function is used, and $f_{NN}(·)$ is the regression function describing the neural network (NN). Here, O is the dimension of network output, and M is the number of samples. For simplicity, the considered NRMSE can be understood as the root mean squared error divided by the maximum element of the absolute difference of all the eigenvectors.

Overall, the proposed scheme can be concluded as follows:

• Step 1: Perform KL expansion on the logarithmic stochastic permeability function $K(\mathit{x},\omega )$;
• Step 2: Generate M realizations of $\left\{{\mathit{\xi }}^{\left(i\right)}\right\}$ and obtain the corresponding BDDC dominant eigenvectors $\left\{y\left({\mathit{\xi }}^{\left(i\right)}\right)\right\}$, which are the training data for the neural network;
• Step 3: Define training conditions and train the neural network;
• Step 4: Examine the network performance whether the NRMSE is satisfied, otherwise, go back to Step 3 and change training conditions.

3. Results and Discussion

In this section, supporting numerical results are presented to show the performance of our proposed learning adaptive BDDC algorithm. We will consider various choices of permeability coefficients $\rho (\mathit{x},\omega )$ with two major stochastic behaviors. As the concern of our numerical tests are on the learning algorithm, we fix all the parameters that are not related. In most of the experiments, $\mathcal{D}\in {\mathbb{R}}^2$ is chosen to be a unit square spatial domain, and it is partitioned into 16 uniform square subdomains with $H=4^{-1}$ as the coarse grid size. Each subdomain is then further divided into uniform grids with a fine grid size $h={32}^{-1}$. We will state clearly if other grid sizes are used. In the following, we will first describe the considered stochastic coefficients and training conditions.

3.1. Choices of Stochastic Coefficients

From the KL expansion expression (11), we know the eigenvalues and eigenfunctions computed from the Fredholm integral equation are related to the covariance function of $K(\mathit{x},\omega )$, which can be treated as the starting point of KL expansion. Throughout the experiments, two specially chosen covariance functions whose eigenvalues and eigenfunctions in the Fredholm integral equation can be computed analytically are considered with different expected values. For $\mathit{x}={(x_1,x_2)}^T$ and $\widehat{\mathit{x}}={({\widehat{x}}_1,{\widehat{x}}_2)}^T$, these two covariance functions are

• Brownian sheet covariance function:

  $$C_K(\mathit{x},\widehat{\mathit{x}})=min(x_1,{\widehat{x}}_1)min(x_2,{\widehat{x}}_2),$$

  where the corresponding eigenvalues and eigenfunctions are

  $$\begin{array}{cc}\hfill {\lambda }_k& ={\displaystyle \frac{16}{\left({(2i-1)}^2{\pi }^2\right)\left({(2j-1)}^2{\pi }^2\right)}},\hfill \\ \hfill f_k\left(\mathit{x}\right)& =2sin\left((i-{\displaystyle \frac{1}{2}})\pi x_1\right)sin\left((j-{\displaystyle \frac{1}{2}})\pi x_2\right).\hfill \end{array}$$
• Exponential covariance function:

  $$C_K(\mathit{x},\widehat{\mathit{x}})={\sigma }_K^{2}exp\left(-{\displaystyle \frac{|x_1-{\widehat{x}}_1|}{{\eta }_1}}-{\displaystyle \frac{|x_2-{\widehat{x}}_2|}{{\eta }_2}}\right),$$

  where ${\sigma }_K^2$ and ${\eta }_i$ are the variance and correlation length in the $x_i$ direction of the process, respectively. The corresponding eigenvalues and eigenfunctions are

  $$\begin{array}{cc}\hfill {\lambda }_k& ={\displaystyle \frac{4{\eta }_1{\eta }_2{\sigma }_K^2}{(r_{1,i}^2{\eta }_1^2+1)(r_{2,j}^2{\eta }_2^2+1)}},\hfill \\ \hfill f_k\left(\mathit{x}\right)& ={\displaystyle \frac{r_{1,i}{\eta }_{1}cos\left(r_{1,i}{x}_1\right)+sin\left(r_{1,i}{x}_1\right)}{\sqrt{(r_{1,i}^2{\eta }_1^2+1)/2+{\eta }_1}}}{\displaystyle \frac{r_{2,j}{\eta }_{2}cos\left(r_{2,j}{x}_2\right)+sin\left(r_{2,j}{x}_2\right)}{\sqrt{(r_{2,j}^2{\eta }_2^2+1)/2+{\eta }_2}}},\hfill \end{array}$$

  where $r_{1,i}$ is the i-th positive root of the characteristic equation

  $$(r_{1,i}^2{\eta }_1^2-1)sin\left(r_{1,i}\right)=2{\eta }_1{r}_{1,i}cos\left(r_{1,i}\right)$$

  in the $x_1$ direction. The same for $r_{2,j}$, which is in the j-th positive root of the characteristic equation in the $x_2$ direction. After arranging the roots $r_{1,i}$ and $r_{2,j}$ in ascending order, we can immediately obtain a monotonically decreasing series of ${\lambda }_k$, which favors our selection of dominant eigenfunctions in the truncated KL expansion with R terms.

For clarification, the positive indices i and j are mapped to index k so that ${\lambda }_k$ is monotonically decreasing. It is noted these eigenvalues and eigenfunctions are just the product of solutions of the one-dimensional process in the Fredholm integral equation, since the covariance functions are separable. First four terms of the eigenfunctions for the second covariance function with $({\sigma }_K^2,{\eta }_1,{\eta }_2)=(1,0.25,0.25)$ are illustrated in Figure 2.

3.2. Adaptive BDDC Parameters and Training Conditions

In the setting of BDDC parameters, we fix the number of dominant eigenvectors on every coarse edge to be 2 and deluxe scaling is always considered. We remark that we have fixed the number of coarse basis functions on each edge. In general, one should vary the number of coarse basis functions according to the local heterogeneities. If the number of coarse basis function is too large, then the coarse problem is more expensive to solve. On the other hand, if the number of coarse basis function is too small, then the number of iteration will increase. It is worth note that, the output layer of neural network is composed by the total dominant eigenvectors on all the coarse edges, so the output layer size O is determined by the coarse grid size H, fine grid size h, and also the number of dominant eigenvectors on every coarse edge. When solve the system by the conjugate gradient method, the iteration is stopped when the relative residual is below ${10}^{-10}$.

In the training, there are different sets of hyperparameter that could affect the training performance of the neural network. Unless specified, the following training conditions are considered:

• Number of truncated terms in KL expansion: $R=4$;
• Number of hidden layers: $L=1$;
• Number of neurons in the hidden layer: 10;
• Number of neurons in the output layer: $O=336$;
• Activation function in hidden layer: hyperbolic tangent function;
• Activation function in output layer: linear function;
• Stopping criteria:

  -

  Minimum value of the cost function gradient: ${10}^{-6}$; or

  -

  Maximum number of training epochs: 1,000,000;

• Sample size of training set: $M=$10,000;
• Sample size of testing set: $M=500$.

After we obtain an accurate neural network from the proposed algorithm, besides the NRMSE of testing set samples, we also present some characteristics of the pre-conditioner based on the approximate eigenvectors from the learning adaptive BDDC algorithm, which include the number of iterations required, minimum and maximum eigenvalues of the pre-conditioner. To show the performance of pre-conditioning on our predicted dominant eigenvectors, we will use the infinity norm to measure the largest difference and also the following symmetric mean absolute percentage error (sMAPE), proposed in [34], to show a relative error type measure with an intuitive range between 0% and 100%:

$$\begin{array}{c}\hfill sMAPE={\displaystyle \frac{1}{M}}\sum _{i=1}^M{\displaystyle \frac{|q_i-{\widehat{q}}_i|}{|q_i|+|{\widehat{q}}_i|}},\end{array}$$

where $q_i$ are some target quantities from pre-conditioning and ${\widehat{q}}_i$ is the same type of quantity as $q_i$ from pre-conditioning using our predicted eigenvectors. In the following subsections, all the computation and results were obtained using MATLAB R2019a with Intel Xeon Gold 6130 CPU and GeForce GTX 1080 Ti GPU in parallel.

3.3. Brownian Sheet Covariance Function

The first set of numerical tests bases on the KL expansion of Brownian sheet covariance function with the following expected functions $E\left[K\right]\left(\mathit{x}\right)$:

($\mathcal{A}$)

${10}^s$

($\mathcal{B}$)

$5sin\left(2\pi x_1\right)sin\left(2\pi x_2\right)+5$.

Expected function $\mathcal{A}$ is a random coefficient as $s\in [-1,1]$ is randomly chosen for each fine grid element, and expected function $\mathcal{B}$ is a smooth trigonometric function with values between $[0,10]$. The common property between these two choices is that the resulting permeability function $\rho (\mathit{x},\omega )$ is highly oscillatory with magnitude order about ${10}^{-2}$ to ${10}^5$, which are good candidates to test our method on oscillatory and high contrast coefficients. Here, we show the appearance of mean functions $\mathcal{A}$ and $\mathcal{B}$ and the corresponding permeability function $\mathcal{K}(\mathit{x},\omega )$ in Figure 3. We can observe that from the first row of Figure 3 to the second row, that is, from the appearances of expected function $E\left[K\right]\left(\mathit{x}\right)$ to the logarithmic permeability coefficient $K(\mathit{x},\omega )$, there are no significant differences except only a minor change on the values on each fine grid element. Nevertheless, these small stochastic changes cause a high contrast $\rho (\mathit{x},\omega )$ after taking exponential function as in the third row of Figure 3.

In

Table 1. Training records when expected functions $\mathcal{A}$ and $\mathcal{B}$ are considered.

Case	Epochs Trained	Training NRMSE	Training Time	Preparation Time for 1 Sample
$\mathcal{A}$	2.29 $\times {10}^5$	1.97$\times {10}^{-2}$	1.98 $\times {10}^3$ s	1.53 s
$\mathcal{B}$	5.39 $\times {10}^4$	3.77 $\times {10}^{-3}$	6.33 $\times {10}^2$ s	1.54 s

, the training results of these two expected functions $\mathcal{A}$ and $\mathcal{B}$ under the training conditions mentioned are presented. Although the random initial choice of conjugate gradient algorithm for minimization may have influences, in general, the neural network of expected function $\mathcal{A}$ usually needs more epochs until stopping criterion is reached and, thus, the training time is also longer. However, this increase in training epochs does not bring a better training NRMSE when compared with expected function $\mathcal{B}$. We can see the same phenomenon when the testing set is considered. The main reason is due to the smoothness of expected function $\mathcal{B}$, the neural network can capture its characteristics easier than a random coefficient.

In the testing results, we consider another set of data containing $M=500$ samples as a testing set. After obtaining the predicted eigenvectors from neural network, we plug it into the BDDC pre-conditioned solver and obtain an estimated pre-conditioner. For the estimated pre-conditioner, we report several properties such as the number of iteration needed for the iterative solver, the maximum and minimum eigenvalues of the pre-conditioned system. Therefore, to better show the performance of network, besides the testing error NRMSE, we also use the sMAPE and $l_{\infty }$ norm of difference of the iteration numbers, the maximum and minimum eigenvalues of the estimated pre-conditioner and the target pre-conditioner. The comparison results are listed in

Table 2. Testing records when expected functions $\mathcal{A}$ and $\mathcal{B}$ are considered.

Case	Testing NRMSE	sMAPE (${\mathit{l}}_{\mathit{\infty }}$ Error) in
		Iteration Number	${\mathit{\lambda }}_{\mathbf{min}}$	${\mathit{\lambda }}_{\mathbf{max}}$
$\mathcal{A}$	6.58$\times {10}^{-2}$	6.94 $\times {10}^{-2}$ (2)	2.58 $\times {10}^{-7}$ (4.02$\times {10}^{-6}$)	4.45 $\times {10}^{-3}$ (3.58 $\times {10}^{-2}$)
$\mathcal{B}$	6.73$\times {10}^{-3}$	6.93 $\times {10}^{-2}$ (1)	3.19 $\times {10}^{-5}$ (1.50$\times {10}^{-4}$)	2.01 $\times {10}^{-3}$ (4.94 $\times {10}^{-2}$)

below.

It is clear that the testing NRMSE of expected function $\mathcal{B}$ is also much smaller than expected function $\mathcal{A}$. One possible reason is because the function used is smooth, and the neural network can better capture the property of the resulting stochastic permeability function. We can see that values of the testing NRMSE of both expected function $\mathcal{A}$ and expected function $\mathcal{B}$ are just a bit larger than the training NRMSE. This suggests that the neural network considered is capable to give a good prediction on dominant eigenvectors in the coarse space even when the magnitude of coefficient function value changes dramatically across each fine grid element. Nevertheless, even though a more accurate set of eigenvectors can be obtained for the expected function $\mathcal{B}$, but its performance after the BDDC pre-conditioning is not always better than that of expected function $\mathcal{A}$. Obviously, $l_{\infty }$ errors in the aspects of ${\lambda }_{max}$ and ${\lambda }_{min}$ are larger for the expected function $\mathcal{B}$.

Detailed comparisons can be found in Figure 4, where the area of overlapping represents how close the specified quantity of the estimated pre-conditioner is compared to the target pre-conditioner. In the first row of figures, we can see the number of iteration required for the predicted eigenvectors is usually more than the target iteration number, however, the differences in the iteration number is not highly related to the difference in minimum and maximum eigenvalues. We can see more examples to verify it in this and later results. The bin width of the histogram of minimum eigenvalues for expected function $\mathcal{A}$ is about ${10}^{-7}$, which is in coherence with the $l_{\infty }$ error, but we eliminate this small magnitude effect and the sMAPE is considered. The high ratio of overlapping area for expected function $\mathcal{A}$ confirms the usage of sMAPE as a good measure, as the resulting sMAPE 2.58 $\times {10}^{-7}$ of expected function $\mathcal{A}$ is much smaller than 3.19$\times {10}^{-5}$ of expected function $\mathcal{B}$.

For the last column of errors in Table 2, the $l_{\infty }$ norm of maximum eigenvalues is larger for expected function $\mathcal{B}$. Nevertheless, it does not mean that the neural network of expected function $\mathcal{B}$ performs worse. From its histogram, we can see an extremely concentrated and overlapping area at the bin around 1.04, moreover, there are some outliers for the prediction results, which is one of the source for a larger $l_{\infty }$ norm. Therefore, the sMAPE can well represent the performance of preconditioning on the estimated results, and we can conclude that both neural networks show a good results on oscillatory and high contrast coefficients, which can be represented by NRMSE and sMAPE. In Section 3.4, besides an artificial coefficient, we will consider one that is closer to a real life case.

Although the $l_{\infty }$ norm seems not describe the characteristics of performance in a full picture, it is still a good and intuitive measure, for example, in the difference of iteration number, we can immediately realize how large will be the worst case.

Further Experiments on Different Training Conditions

To emphasize the changes of using different training conditions, we display the results of increasing the truncated terms R in KL expansion and using less training samples. The expected functions and stochastic property considered are the same as above to show the effect of changing conditions only. For clarification, there are two experiments each with only one condition is modified and other conditions are the same as stated in Section 3.2. We state the new conditions as the follows:

• Modified experiment 1: $R=9$;
• Modified experiment 2: Training sample is 5000.

When R is increased, we can see in

Table 3. Training records when modified experiment 1 is considered.

Case	Epochs Trained	Training NRMSE	Training Time	Preparation Time for 1 Sample
$\mathcal{A}$	2.17 $\times {10}^5$	2.74 $\times {10}^{-2}$	2.47 $\times {10}^3$ s	1.77 s
$\mathcal{B}$	1.15 $\times {10}^5$	4.02 $\times {10}^{-3}$	1.19 $\times {10}^3$ s	1.75 s

, the training time increases as the total number of parameters increases with the input layer. In

Table 4. Testing records when modified experiment 1 is considered.

Case	Testing NRMSE	sMAPE (${\mathit{l}}_{\mathit{\infty }}$ Error) in
		Iteration Number	${\mathit{\lambda }}_{\mathbf{min}}$	${\mathit{\lambda }}_{\mathbf{max}}$
$\mathcal{A}$	6.97 $\times {10}^{-2}$	5.85 $\times {10}^{-2}$ (1)	2.15 $\times {10}^{-6}$ (1.35 $\times {10}^{-5}$)	2.29 $\times {10}^{-2}$ (9.05 $\times {10}^{-2}$)
$\mathcal{B}$	5.40 $\times {10}^{-3}$	6.73 $\times {10}^{-2}$ (1)	2.99$\times {10}^{-5}$ (1.76 $\times {10}^{-4}$)	1.96 $\times {10}^{-3}$ (5.64 $\times {10}^{-2}$)

, the testing results of expected function $\mathcal{B}$ clearly have better performance than in Table 2, which suggests that a larger R could bring a smaller NRMSE.

The training and testing results of modified experiment 2 are shown in

Table 5. Training records when modified experiment 2 is considered.

Case	Epochs Trained	Training NRMSE	Training Time	Preparation Time for 1 Sample
$\mathcal{A}$	9.12 $\times {10}^4$	3.34 $\times {10}^{-2}$	8.94 $\times {10}^2$ s	1.53 s
$\mathcal{B}$	1.01 $\times {10}^5$	4.91 $\times {10}^{-3}$	1.12 $\times {10}^3$ s	1.54 s

and

Table 6. Testing records when modified experiment 2 is considered.

Case	Testing NRMSE	sMAPE (${\mathit{l}}_{\mathit{\infty }}$ Error) in
		Iteration Number	${\mathit{\lambda }}_{\mathbf{min}}$	${\mathit{\lambda }}_{\mathbf{max}}$
$\mathcal{A}$	7.82 $\times {10}^{-2}$	6.68 $\times {10}^{-2}$ (2)	7.15 $\times {10}^{-7}$ (6.59 $\times {10}^{-6}$)	3.82 $\times {10}^{-3}$ (1.57 $\times {10}^{-1}$)
$\mathcal{B}$	8.43 $\times {10}^{-3}$	6.91 $\times {10}^{-2}$ (1)	2.73 $\times {10}^{-5}$ (1.64 $\times {10}^{-4}$)	2.39 $\times {10}^{-3}$ (6.61 $\times {10}^{-2}$)

, respectively. Obviously, the testing results in Table 6 are worse than when 10,000 training samples are used. Moreover, as the training time does not have a distinct reduction, and with the aid of parallel computation, more training samples can be easily obtained. Therefore, using less samples is not always a good choice, we should use appropriate number of samples.

3.4. Exponential Covariance Function

To test our method on realistic highly varying coefficients, we consider expected functions that come from the second model of the 10th SPE Comparative Solution Project (SPE10), with the KL expanded exponential covariance function as the stochastic source. For clarification, in this set of experiment, we first use the modified permeability of Layer 35 in the x-direction of SPE10 data as the expected function $E\left[K\right]\left(\mathit{x}\right)$ with $({\sigma }_K^2,{\eta }_1,{\eta }_2)=(1,0.25,0.25)$ in the exponential covariance function $C_K(\mathit{x},\widehat{\mathit{x}})$. We then train the neural network, and the resulting neural network is used to test its generalization capacity on the following different situations:

• Different stochastic behavior:
  • All remain unchanged except $({\eta }_1,{\eta }_2)=(0.2,0.125)$ or $({\eta }_1,{\eta }_2)=(1,1)$;
• Different mean permeability function:
  • The expected function is changed to the Layers 1, 4,34 in the x-direction of SPE10 data, but the stochastic parameters are unchanged, i.e., $({\sigma }_K^2,{\eta }_1,{\eta }_2)=(1,0.25,0.25)$.

It is of note that, all the permeability fields used are modified into our computational domain and shown in Figure 5 for a clear comparison. In the Layer 35 permeability field, four sharp properties can be observed. There are one blue strip on the very left side, and a reversed but narrower red strip is next to it. In the top right corner, a low permeability area is shown and just below it, a small high permeability area is located in the bottom right corner. So a similar but with slightly different properties of Layer 34, and two Layers 1, 4 with more distinct differences are chosen to test the generalization capacity of our neural network.

Before showing the training and testing records of neural network, we further present the logarithmic permeability $K(\mathit{x},\omega )$, and permeability $\rho (\mathit{x},\omega )$ when $({\sigma }_K^2,{\eta }_1,{\eta }_2)=(1,0.25,0.25)$, $(1,0.2,0.125)$ and $(1,1,1)$ in Figure 6. The sets of $\left\{{\xi }_i\right\}$ used in these three rows of realizations are the same. We can thus see the change in parameters does cause the stochastic behavior to change, too. We present the training and testing records of neural network as below in

Table 7. Training record when Layer 35 permeability is considered.

Case	Epochs Trained	Training NRMSE	Training Time	Preparation Time for 1 Sample
Layer 35	2.62 $\times {10}^5$	1.52 $\times {10}^{-2}$	1.07 $\times {10}^3$ s	12.34 s

and

Table 8. Testing records when Layer 1, 4, 34, 35, 35 *, and 35 ** are considered.

Case	Testing NRMSE	sMAPE (${\mathit{l}}_{\mathit{\infty }}$ Error) in
		Iteration Number	${\mathit{\lambda }}_{\mathbf{min}}$	${\mathit{\lambda }}_{\mathbf{max}}$
Layer 1	8.02 $\times {10}^{-2}$	1.04 $\times {10}^{-1}$ (2)	3.74 $\times {10}^{-6}$ (3.71 $\times {10}^{-5}$)	4.00 $\times {10}^{-2}$ (1.36 $\times {10}^{-1}$)
Layer 4	6.95 $\times {10}^{-2}$	7.18 $\times {10}^{-2}$ (1)	1.11 $\times {10}^{-5}$ (1.36 $\times {10}^{-4}$)	1.61 $\times {10}^{-2}$ (6.85 $\times {10}^{-2}$)
Layer 34	5.08 $\times {10}^{-2}$	2.52 $\times {10}^{-3}$ (1)	3.25 $\times {10}^{-5}$ (1.74 $\times {10}^{-4}$)	3.06 $\times {10}^{-2}$ (8.47 $\times {10}^{-2}$)
Layer 35	2.48 $\times {10}^{-2}$	3.65 $\times {10}^{-2}$ (1)	7.52 $\times {10}^{-6}$ (7.50 $\times {10}^{-5}$)	1.36 $\times {10}^{-3}$ (6.66 $\times {10}^{-2}$)
Layer 35 *	2.27 $\times {10}^{-2}$	3.89 $\times {10}^{-2}$ (1)	7.57 $\times {10}^{-6}$ (7.17 $\times {10}^{-5}$)	1.27 $\times {10}^{-3}$ (6.59 $\times {10}^{-2}$)
Layer 35 **	2.75 $\times {10}^{-2}$	4.29 $\times {10}^{-2}$ (1)	8.52 $\times {10}^{-6}$ (7.13 $\times {10}^{-5}$)	1.81 $\times {10}^{-3}$ (8.10 $\times {10}^{-2}$)

, where we specify the expected functions considered in the first column.

In Table 8, we clarify that the results of Layer 35 * also use Layer 35 permeability as the expected function but with $({\eta }_1,{\eta }_2)=(0.2,0.125)$. Similarly, Layer 35 ** is with another set of parameters $({\eta }_1,{\eta }_2)=(1,1)$. Each row corresponds to the results of different testing sets, however, all the results are obtained using the same neural network of Layer 35 in Table 7. Therefore, the testing results of Layer 35 are used as a reference to decide whether the results of other testing sets are good or not.

We first focus on the testing results of Layer 35. From the corresponding histograms in the left column of Figure 7, besides the iteration number required for predicted eigenvectors are still usually larger by 1, we can observe a large area of overlapping, in particular, there is only 8.4% of data that are not overlapped in the histogram of maximum eigenvalues. Moreover, the small NRMSE and sMAPE again confirm with the graphical results. For generalization tests, we then consider other testing sets with different stochastic behaviors and expected functions.

When the stochastic behavior of Layer 35 is changed to $({\sigma }_K^2,{\eta }_1,{\eta }_2)=(1,0.2,0.125)$ and $({\sigma }_K^2,{\eta }_1,{\eta }_2)=(1,1,1)$, we list the corresponding sMAPE and $l_{\infty }$ norms of differences in iteration numbers, maximum eigenvalues and minimum eigenvalues of pre-conditioners in the last two rows of Table 8. We can see the values of errors are very similar to the row of Layer 35 and the NRMSE are very close. Therefore, we expect that the Layer 35 neural network can generalize the coefficient functions with similar stochastic properties, and give a good approximation on dominant eigenvectors for the BDDC pre-conditioner. It is then verified by the histograms in the middle column of Figure 7. Moreover, a large population of samples are concentrated around 0 in the histograms of difference between the estimated pre-conditioner and target pre-conditioner in Figure 8. All these show our method performs well with stochastic oscillatory and high contrast coefficients, and also coefficient functions with similar stochastic properties.

Finally, we compare the performances when the expected function is changed to Layer 1, 4, and 34. Note, the NRMSE, sMAPE and $l_{\infty }$ errors of the three layers are almost all much higher than those of Layer 35, which are crucial clues for a worse performance. It is clear that when the permeability fields do not have many similarities such as Layer 1 and 4, the testing NRMSE increases distinctly. Even though the permeability fields of Layer 34 and Layer 35 share some common properties and have a near magnitude, unlike previous results for Layer 35 and Layer 35 *, we can see in the right column of Figure 7, the prediction output and target output on maximum and minimum eigenvalues of pre-conditioners are even two populations with different centres and only a little overlapping area. The main reason is a minor difference in $K(\mathit{x},\omega )$ can already cause a huge change after taking exponential function. Therefore, when $E\left[K\right]\left(\mathit{x}\right)$ is changed to another entirely different mean permeability with just some common properties, the feed-forward neural network may not be able to capture these new characteristics, which will be part of our future research.

3.5. Exponential Covariance Function with Multiple Hidden Layers

It is noted that after taking exponential function, the permeability fields in previous experiments only show some global peaks, the single layer structure of neural network may then suffice to give an accurate prediction. Therefore, in this subsection, we illustrate an example with a high contrast coefficient cutting more interfaces and a realization of this new coefficient function $\widehat{\rho }(\mathit{x},\omega )$ is displayed in Figure 9. It is of note that, the stochastic properties considered are still the exponential covariance function with $({\sigma }_K^2,{\eta }_1,{\eta }_2)=(1,0.25,0.25)$.

To extract more features and details from the new coefficient, we consider a much finer grid sizes, which are $H=1/10$ and $h=1/600$, and the spatial domain is still the unit square ${[0,1]}^2$. With respect to this change, we change the following training conditions accordingly:

• Number of hidden layers: $L=1$ or $L=2$;
• Number of neurons in the hidden layers: $n^{\left(1\right)}=10$ and $n^{\left(2\right)}=0,5$ or 7;
• Number of neurons in the output layer: $O=$21,240;
• Sample size of training set: $M=2000$;
• Sample size of testing set: $M=100$.

Other conditions are kept to be unchanged as in Section 3.2. We show the new training and testing records in

Table 9. Training records when new coefficient function $\widehat{\rho }(\mathit{x},\omega )$ is considered.

$({\mathit{n}}^{\left(1\right)},{\mathit{n}}^{\left(2\right)})$	Epochs Trained	Training NRMSE	Training Time	Preparation Time for 1 Sample
(10,0)	1.00 $\times {10}^6$	3.81 $\times {10}^{-2}$	4.89 $\times {10}^4$ s	130.77 s
(10,5)	6.46 $\times {10}^5$	4.44 $\times {10}^{-2}$	2.65 $\times {10}^4$ s
(10,7)	1.00 $\times {10}^6$	4.11 $\times {10}^{-2}$	4.62 $\times {10}^4$ s

and

Table 10. Testing records when new coefficient function $\widehat{\rho }(\mathit{x},\omega )$ is considered.

$({\mathit{n}}^{\left(1\right)},{\mathit{n}}^{\left(2\right)})$	Testing NRMSE	sMAPE (${\mathit{l}}_{\mathit{\infty }}$ Error) in
		Iteration Number	${\mathit{\lambda }}_{\mathbf{min}}$	${\mathit{\lambda }}_{\mathbf{max}}$
(10,0)	5.47 $\times {10}^{-2}$	2.12 $\times {10}^{-2}$ (4)	1.86 $\times {10}^{-6}$ (2.98 $\times {10}^{-5}$)	1.42 $\times {10}^{-1}$ (4.34 $\times {10}^2$)
(10,5)	3.57 $\times {10}^{-2}$	7.81 $\times {10}^{-3}$ (5)	1.21 $\times {10}^{-6}$ (3.53 $\times {10}^{-5}$)	4.55 $\times {10}^{-2}$ (7.89 $\times {10}^1$)
(10,7)	2.30 $\times {10}^{-2}$	1.12 $\times {10}^{-2}$ (5)	1.13 $\times {10}^{-6}$ (2.50 $\times {10}^{-5}$)	9.18 $\times {10}^{-2}$ (7.01 $\times {10}^1$)

.

On one hand, when the number of neurons in the second hidden layer $n^{\left(2\right)}$ increases, the training results do not show a distinct difference, the three training NRMSE are all around $0.04$ with a much longer training time than previous experiments. However, on the other hand, the testing results do show some improvements when $n^{\left(2\right)}$ increases. Obviously, the testing NRMSE and pre-conditioner errors improve when a second hidden layer is introduced. Moreover, the testing NRMSE also decreases when $n^{\left(2\right)}$ increases. We can see the testing NRMSE of multiple hidden layers are slightly lower than the corresponding training NRMSE.

This indicates that for more complicated coefficient functions, introduction of multiple hidden layers could help to bring better performance. However, as it also increases the complexity of neural network and the training time, thus, a single layer of neural network should be considered first if the problem is not complicated. In Figure 10, comparisons of the iteration number, maximum and minimum eigenvalues are illustrated in histograms for easy observation. Unlike the histograms in Figure 7, we can see more iterations are needed for solving the problem and more outliers exist in the maximum eigenvalue comparison as the problem itself becomes more complicated than previous experiments.

4. Conclusions

A new learning adaptive BDDC algorithm is introduced and it consists of three main parts, which are the Gaussian random variables in Karhunen–Loève expansion, artificial neural network and the dominant eigenvectors obtained from the adaptive BDDC algorithm in the coarse spaces. The considered neural network acts as a fast computational tool, which turns the input Gaussian random variables into predicted dominant eigenvectors as output. In addition, the neural network in the proposed algorithm is suitable for other permeability coefficients with similar stochastic properties for generalization purpose. Numerical results confirm the efficiency of the proposed algorithm and the generalization abilities. On the other hand, we currently just use the simplest feed-forward neural network structure without retraining when new characteristics are entered. In our future plan, we will improve the neural network structure that matches with the computational domain and some adaptive BDDC features. Deep learning technique will also be included for higher accuracy in prediction, pre-conditioning, and also better generalization ability on wider aspects of applications.