A Multilevel, Hierarchical Sampling Technique for Spatially Correlated Random Fields

TL;DR

This paper introduces a scalable, hierarchical sampling method for spatially correlated random fields using a multilevel decomposition of stochastic reaction-diffusion equations, improving efficiency for large-scale uncertainty propagation.

Contribution

It presents a novel multilevel, hierarchical sampling technique that overcomes computational limitations of traditional methods like KL expansion for large-scale spatial fields.

Findings

Method is scalable and efficient for large problems.

Numerical experiments demonstrate effective multilevel Monte Carlo simulations.

Applicable to subsurface porous media flow problems.

Abstract

We propose an alternative method to generate samples of a spatially correlated random field with applications to large-scale problems for forward propagation of uncertainty. A classical approach for generating these samples is the Karhunen-Lo\`{e}ve (KL) decomposition. However, the KL expansion requires solving a dense eigenvalue problem and is therefore computationally infeasible for large-scale problems. Sampling methods based on stochastic partial differential equations provide a highly scalable way to sample Gaussian fields, but the resulting parametrization is mesh dependent. We propose a multilevel decomposition of the stochastic field to allow for scalable, hierarchical sampling based on solving a mixed finite element formulation of a stochastic reaction-diffusion equation with a random, white noise source function. Numerical experiments are presented to demonstrate the scalability of the sampling method as well as numerical results of multilevel Monte Carlo simulations for a subsurface porous media flow application using the proposed sampling method.

Figures

Equations

$$\mathbf{q}(\mathbf{x})+k(\mathbf{x})\nabla p(\mathbf{x})$$

$$\mathrm{div}\mathbf{q}(\mathbf{x})$$

$$-\Delta \theta (\mathbf{x},\omega )+{\kappa }^2\theta (\mathbf{x},\omega )=\mathcal{W}(\mathbf{x},\omega )$$

$$\mathrm{cov}(\mathbf{x},\mathbf{y})=\frac{{\sigma }^2}{2^{\nu -1}\Gamma (\nu )}(\kappa \Vert \mathbf{y}-\mathbf{x}\Vert {)}^{\nu }{K}_{\nu }(\kappa \Vert \mathbf{y}-\mathbf{x}\Vert ),$$

$$\mathrm{cov}(\mathbf{x},\mathbf{y})={\sigma }^2{e}^{-\kappa \Vert \mathbf{y}-\mathbf{x}\Vert }.$$

$$\mathcal{C}v(\mathbf{x})={\int }_D\mathrm{cov}(\mathbf{x},{\mathbf{x}}^{\prime })v({\mathbf{x}}^{\prime })d{\mathbf{x}}^{\prime }.$$

$$(q,\mathcal{C}v{)}_D=\lambda (q,v{)}_D\forall q\in L^2(D),$$

$$\theta (\mathbf{x},\omega )=\sum _{i=0}^{\infty }{\xi }_i(\omega )\sqrt{{\lambda }_i}{v}_i(\mathbf{x}),{\xi }_i(\omega )\sim N(0,{\sigma }^2)\text{i.i.d.}$$

$${\theta }_h(\omega )=\sum _{i=0}^{M_h}{\xi }_i(\omega )\sqrt{{\lambda }_{h,i}}{v}_{h,i}(\mathbf{x}),{\xi }_i(\omega )\sim N(0,{\sigma }^2)\text{i.i.d},$$

$$(q_h,{\mathcal{C}}_h{v}_{h,i}{)}_h={\lambda }_{h,i}(q_h,v_{h,i}{)}_h$$

$$({\kappa }^2-\Delta {)}^{\alpha /2}\theta (\mathbf{x},\omega )=g\mathcal{W}(\mathbf{x},\omega )\mathbf{x}\in {\mathbb{R}}^d,\alpha =\nu +\frac{d}{2},\kappa >0,\nu >0,$$

$$g=(4\pi {)}^{d/4}{\kappa }^{\nu }\sqrt{\frac{\Gamma \left(\nu +d/2\right)}{\Gamma (\nu )}},$$

$$({\kappa }^2-\Delta )\theta (\mathbf{x},\omega )=g\mathcal{W}(\mathbf{x},\omega ).$$

$$(\theta ,q)={\int }_D\theta qd\mathbf{x}\text{mbox}and(\mathbf{u},\mathbf{v})={\int }_D\mathbf{u}\cdot \mathbf{v}d\mathbf{x}.$$

$$\begin{array}[]{lcr}\mathbf{R}=H(\operatorname{div};D):=\left\{\mathbf{u}\in\mathbf{L}^{2}(D)\mid\operatorname{div}\,\mathbf{u}\in L^{2}(D)\text{ and }\mathbf{u}\cdot\mathbf{n}=0\text{ on }\partial D\right\}\text{ and }\\ \Theta=L^{2}(D).\end{array}$$

$$\begin{array}[]{lll}m(\mathbf{u},\mathbf{v})&:=(\mathbf{u},\mathbf{v})&\forall\,\mathbf{u},\mathbf{v}\in\mathbf{R},\\ w(\theta,q)&:=(\theta,q)&\forall\,\theta,q\in\Theta,\\ b(\mathbf{u},q)&:=(\operatorname{div}\,\mathbf{u},q)&\forall\,\mathbf{u}\in\mathbf{R},q\in\Theta,\\ \end{array}$$

$$F^{\mathcal{W}}(q):=(\mathcal{W},q)\forall q\in \Theta .$$

$$\begin{array}[]{ll}m(\mathbf{u}_{h},\mathbf{v}_{h})+b(\mathbf{v}_{h},\theta_{h})=0&\forall\mathbf{v}_{h}\in\mathbf{R}_{h}\\ b(\mathbf{u}_{h},q_{h})-\kappa^{2}\,w(\theta_{h},q_{h})=-g\,F^{\mathcal{W}}(q_{h})&\forall q_{h}\in\Theta_{h}.\end{array}$$

$$\mathbb{E}[(q_i,\mathcal{W})]$$

$$\mathrm{cov}\left((q_i,\mathcal{W}),(q_j,\mathcal{W})\right)$$

$$f_h=W_h^{\frac{1}{2}}{\xi }_h(\omega ),$$

$$\mathcal{A}_{h}U_{h}=\begin{array}[]{lcr}\begin{bmatrix}M_{h}&B_{h}^{T}\\ B_{h}&-\kappa^{2}W_{h}\end{bmatrix}\begin{bmatrix}\mathbf{u}_{h}\\ \theta_{h}\end{bmatrix}=\begin{bmatrix}0\\ -g\,f_{h}(\omega)\end{bmatrix}=F_{h},\par\end{array}$$

$$m({\nabla }_h{\theta }_h,{\nabla }_h{q}_h)+{\kappa }^{2}w({\theta }_h,q_h)=g{F}^W(q_h),\forall q_h\in {\Theta }_h.$$

$${\mathcal{A}}_h=\left[\begin{array}{cc}{A}_h& B_h^T\\ 0& -{\kappa }^2{W}_h\end{array}\right]\left[\begin{array}{cc}I& 0\\ -{\kappa }^{-2}{W}_h^{-1}{B}_h& I\end{array}\right],$$

$$A_h{\mathbf{u}}_h=-g{\kappa }^{-2}{B}_h^T{f}_h(\omega ),$$

$${\theta }_h={\kappa }^{-2}{W}_h^{-1}\left(B_h{\mathbf{u}}_h+g{f}_h(\omega )\right).$$

$$a({\mathbf{u}}_h,{\mathbf{v}}_h):=({\mathbf{u}}_h,{\mathbf{v}}_h)+{\kappa }^{-2}(\mathrm{div}{\mathbf{u}}_h,\mathrm{div}{\mathbf{v}}_h),$$

$$P_{\theta }:{\Theta }_{\ell +1}\to {\Theta }_{\ell }\text{for}\ell =0,\dots ,L-1.$$

$$P_{\mathbf{u}}:{\mathbf{R}}_{\ell +1}\to {\mathbf{R}}_{\ell }\text{for}\ell =0,\dots ,L-1.$$

$$\mathcal{P}=\left[\begin{array}{cc}{P}_{\mathbf{u}}& 0\\ 0& P_{\theta }\end{array}\right],$$

$$\left[\begin{array}{c}{\mathbf{u}}_{\ell }\\ {\theta }_{\ell }(\omega )\end{array}\right]={\mathcal{A}}_{\ell }^{-1}\left[\begin{array}{c}0\\ -g{W}_{\ell }^{1/2}{\xi }_{\ell }(\omega )\end{array}\right]$$

Full text

{AMS}

65C05, 60H15, 35R60, 65N30, 65M75, 65C30

1 Introduction

Generating spatially correlated Gaussian random fields with specific statistical properties is an important component and active topic of research in a diverse array of application areas such as ecology, meteorology and geology [10, 15, 22]. Gaussian random fields are collections of random variables indexed by elements from a multidimensional space with the property that their joint distribution is Gaussian. These types of fields are specified by expectations and positive semi-definite covariance functions and are generally good models for many phenomena; see, e.g., [5].

In this paper we focus on geophysics applications where the spatially correlated random field represents the permeability of a porous medium. Specifically, we consider subsurface flow in a specified domain, $D\subset\mathbb{R}^{d}$, governed by Darcy’s law,

[TABLE]

where ${\bf q}({\bf x})$ is the fluid velocity, $p({\bf x})$ is the fluid pressure and $k({\bf x})$ the hydraulic conductivity, measuring the transmissibility of the fluid through the porous medium.

In practice, only small portions of the hydraulic conductivity are known via noisy data measurements and this contributes a large source of uncertainty into the Darcy model. Understanding the effects that such uncertainties introduce into the model, or forward propagation uncertainty quantification (UQ) [21, 37], is quite important. To quantify this uncertainty we consider Monte Carlo methods, in particular multilevel Monte Carlo methods (MLMC), where equations (1)-(2) are solved repeatedly using different realizations of the hydraulic conductivity field. In multilevel Monte Carlo methods, more accurate (and expensive) simulations are run with fewer samples, while less accurate (and inexpensive) simulations are run with a larger number of samples. At the end of these simulations, quantities of interest depending on the velocity and/or pressure, such as the effective permeability, can be computed. Since a large number of Monte Carlo simulations are needed to produce accurate quantities of interest, it is essential to have efficient algorithms to rapidly generate these different realizations.

A standard way of modeling the hydraulic conductivity is as a log-normal field, $k({\bf x},\omega)=e^{\theta({\bf x},\omega)}$, where $\theta({\bf x},\omega)$ is a random field with a specified covariance structure [13, 18]. A classical, and well-studied approach for generating these samples with the desired statistical properties is the Karhunen-Loève decomposition. This sampling technique amounts to solving an eigenvalue problem with the dense covariance matrix. The samples are then different (random) linear combinations of the eigenvectors. However, sampling in this manner suffers from a significant computational cost and high memory requirement. In fact, the computational complexity for solving a dense eigenvalue problem grows cubically with the size of the covariance matrix and the memory necessary to store all the eigenvectors in the Karhunen-Loève decomposition scales quadratically, making such methods prohibitively expensive for large-scale problems of interest. New approaches based on randomized methods and hierarchical semi-separable matrices can drastically reduce the cost of solving the eigenvalue problem (see e.g. [36]), however they allow to compute only the dominant eigenmodes of the Karhunen-Loève expansion and, therefore, introduce bias in the sampling.

A different approach for generating random field samples is by solving a stochastic partial differential equation (SPDE) with a white noise source function [39, 40, 29]. To use the inverse of an elliptic differential operator as covariance function is a common approach for the solution of large-scale Bayesian inverse problem governed by PDE forward models, see e.g. [38, 11], as it allows for efficient evaluation of the covariance operator using a fast and scalable multigrid solver. The authors in [29] provide a link between Karhunen-Loève sampling from a Matérn distribution and SPDE sampling, further motivating the SPDE approach. For example, one can solve the following reaction-diffusion equation

[TABLE]

where $\theta({\bf x},\omega)$ is the logarithm of the hydraulic conductivity, $\kappa^{2}$ is a constant depending on the correlation length, and $\mathcal{W}({\bf x},\omega)$ is a white noise function. This sampling approach has two significant benefits: the first is avoiding the computational cost incurred by solving a dense eigenvalue problem and the second is that optimal solution methods for solving sparse linear systems arising from the finite element discretization of (3) can be applied. Despite the benefits of this approach, it is not without imperfections as generated samples contain artificial boundary effects.

Compared to the work in [29], three new ideas are introduced in this paper: i) a mixed discretization of the SPDE, ii) a hierarchical version of the sampler, and iii) mitigation of boundary artifacts using embedded domains. Specifically, given a general unstructured fine grid, we construct a hierarchy of algebraically coarsened grids and finite element spaces using the element-based algebraic multigrid techniques (AMGe) presented in [27, 26], which provide coarse spaces with improved approximation properties than the ones from the original work [34]. Then we generate random field samples at each level of the hierarchy by solving linear systems arising from the mixed discretization of (3). This allows us to compute our (piecewise constant) samples in a hierarchical fashion (as needed in MLMC simulations), while leveraging existing scalable methods and software for solving deterministic PDEs. We remark that while a hierarchical sampler can be constructed in a similar way for the primal formulation of the SPDE using geometric multigrid hierarchies, our framework offers more flexibility with respect to the geometry of the physical domain (since it does not require a sequence of nested grids). In addition, the computational advantages of resorting to the mixed formulation are twofold. First, the finite element discretization of the SPDE (3) requires the assembly of the square root of a mass matrix. In [29], the use of a continuous Galerkin (CG) finite element space in the discretization of the primal formulation leads to a non-diagonal mass matrix and mass lumping is used to make computation of the square-root tractable at the cost of accuracy. On the contrary, the finite element pair of lowest Raviart-Thomas and piecewise constant functions in the mixed formulation leads to a diagonal mass matrix for the variable $\theta$. Second, in the AMGe framework, the construction of $L^{2}$ and $H(\operatorname{div})$-conforming spaces (required in the mixed formulation) on agglomerated meshes is simpler than the one of $H^{1}$-conforming spaces (required by the primal formulation), and also leads to sparser, i.e. faster to apply, grid transfer operators.

Finally, we can not stress enough that performing realistic simulations with MLMC requires a considerable computational cost. To make the method more feasible in practice, parallelism must be fully exploited. Our work focuses solely on parallelism across the spatial domain in computing realizations of the input random field, then performing the subsequent solve of the model of interest. Further parallelism could be added using the scheduling approaches suggested in [20], where the authors investigate the complex task of scheduling parallel tasks within and across levels of MLMC.

The paper is structured as follows. In Section 2 we give an overview of both the classical KL expansion and the SPDE sampling approach. The hierarchical SPDE sampling procedure is introduced and examined in Section 3. Section 4 contains a brief review of the multilevel Monte Carlo method (MLMC), a scalable alternative to standard Monte Carlo methods that uses the solution of the PDE on a hierarchy of grid to effectively reduce the variance of the estimator. Here, we also present our numerical results when our proposed sampling technique is used inside the MLMC method. Lastly, Section 5 contains our concluding remarks.

2 Sampling from Gaussian Random Fields

In this section we review two ways of generating samples of a log-normal random field. These are collections of random variables, $\{\theta({\mathbf{x}},\omega):~{}{\mathbf{x}}\in D,~{}\omega\in\Omega\}$, where $\Omega$ is the sample space for the probability space $(\Omega,\mathcal{F},\mathbb{P})$. For a fixed point, ${\bf x}_{0}\in D$, $\theta({\mathbf{x}}_{0},\omega)$ is a random variable. For a fixed $\omega_{0}\in\Omega$, $\theta({\mathbf{x}},\omega_{0})$ is a deterministic function called a realization or sample of the random field.

Random fields with a particular covariance functions are used to model the correlation between points in spatial data. In geostatistical applications, the logarithm of the hydraulic conductivity, ${\theta({\mathbf{x}},\omega)=\log{(k({\bf x},\omega))}}$ is modeled as a random field with a specific class of covariance function.

In particular, we consider a stationary isotropic Gaussian field with the widely used class of Matérn covariance functions [32], given by

[TABLE]

where $\sigma^{2}$ is the marginal variance, $\nu>0$ determines the mean-square differentiability of the underlying process, $\kappa>0$ is a scaling factor inversely proportional to the correlation length, $\Gamma(\nu)$ is the gamma function, and $K_{\nu}$ is the modified Bessel function of the second kind. When $\nu=\frac{1}{2}$ (4) reduces to the common exponential covariance function, given by

[TABLE]

In this section we consider two ways for generating such realizations or samples of random fields. The first method we consider is the classical Karhunen-Loève expansion. The second is based on solving a particular stochastic reaction-diffusion equation with a white-noise source function.

2.1 The Karhunen-Loève Expansion

The Karhunen-Loève (KL) expansion of a second-order random field 222 A random field $\theta(\mathbf{x},\omega)$ is second-order if for each $\mathbf{x}\in D$ the random variable $\theta(\mathbf{x})$ has finite variance. provides a series representation using the orthonormal basis provided by the eigenfunctions of the underlying covariance operator [30]. For a bounded regular domain $D$ we define the convolution operator

[TABLE]

Since $\mathcal{C}$ is a compact operator, the eigenvalue problem

[TABLE]

admits a countable sequence of eigenpairs $(\lambda_{i},v_{i})$ where $\lim_{i\rightarrow\infty}\lambda_{i}=0.$ A Gaussian random field can be expanded as

[TABLE]

In practice, a discrete version of the eigenvalue problem (4) is computed. A triangulation of the domain $\mathcal{T}_{h}\subset D$ is generated with the discrete function space $\Theta_{h}\subset L^{2}(D)$ of piecewise constant functions. Then, we can compute realizations of the Gaussian field with a Karhunen-Loève expansion truncated after $M_{h}$ terms as

[TABLE]

where the pairs $(\lambda_{i},v_{i})$ solve the generalized eigenvalue problem:

[TABLE]

with inner product defined as $(r,s)_{h}=\int_{\mathcal{T}_{h}}r\,s\,d\mathbf{x}$ for $r,s\in\Theta_{h}.$

One of the main reasons for sampling using this approach is that any truncation of the KL expansion gives a sample with the minimal mean square error, see e.g. [24]. However, this method is not computationally feasible due to the high computational cost for solving a dense eigenvalue problem. On a mesh with $N$ degrees of freedom, the cost to factor a dense covariance matrix is $\mathcal{O}(N^{3})$. This makes generating samples in this manner impractical for finely resolved meshes where the number of degrees of freedom could be of the order of millions. This motivates other sampling methods that do not suffer from such a complexity cost.

2.2 Stochastic PDE Sampling

An important link between Gaussian fields and Gaussian Markov random fields is established in [29], where a random process on $\mathbb{R}^{d}$ with a Matérn covariance function can be obtained as the solution of a particular stochastic partial differential equation (SPDE). This provides an alternative method for computing the desired samples of a Gaussian random field via a SPDE, as opposed to the computationally intensive KL expansion.

Realizations of a Gaussian random field with an underlying Matérn covariance, $\theta({\bf x},\omega)$, solve the following linear stochastic PDE:

[TABLE]

for some other realization of the standard Gaussian white noise function $\mathcal{W}$ with scaling factor $g$, [39, 40]. The scaling factor is chosen to be

[TABLE]

to impose unit marginal variance of $\theta(\mathbf{x},\omega)$ [10].

It is worth noting that if $\nu=1$ in two dimensions, and $\nu=\frac{1}{2}$ in three dimensions, then (8) reduces to the following standard reaction-diffusion equation,

[TABLE]

Additionally in three dimensions, realizations of a Gaussian random field with exponential covariance are solutions of (9). Thus a scalable sampling alternative is equivalent to efficiently solving the stochastic reaction-diffusion equation given by (9) where we are able to leverage existing scalable solution strategies. We remark that, as explained in [29], samples from Matérn distributions with $\nu=2k+1$ (in two spatial dimensions) and from Matérn distributions with $\nu=2k+\frac{1}{2}$ (in three spatial dimensions), can be obtained by recursively solving $(\kappa^{2}-\Delta)\theta_{i+1}(\mathbf{x},\omega)=\theta_{i}(\mathbf{x},\omega)$ for $i=1,\ldots,k$, where $\theta_{0}(\mathbf{x},\omega)$ is the solution of (9).

Equation (9) is discretized using a mixed finite element method, [9, 17].

Following standard notation, for scalar functions $\theta,q\in L^{2}(D)$ and vector functions $\mathbf{u},\mathbf{v}\in\mathbf{L}^{2}(D)=[L^{2}(D)]^{d}$, we define the inner products:

[TABLE]

We also define the functional spaces $\mathbf{R}$ and $\Theta$ as

[TABLE]

Finally, we introduce the bilinear forms

[TABLE]

and the linear form

[TABLE]

Following standard finite element techniques, let $\mathbf{R}_{h}\subset\mathbf{R}$ denote the lowest order Raviart-Thomas finite element space and $\Theta_{h}\subset\Theta$ denote the finite element space of piecewise constant functions. Then, a mixed finite element discretization of (9) reads

Problem 2.1**.**

Find $(\mathbf{u}_{h},\theta_{h})\in\mathbf{R}_{h}\times\Theta_{h}$ such that

[TABLE]

*with essential boundary conditions $\mathbf{u}_{h}\cdot\mathbf{n}=0$.

Remark**.**

*The choice of a low order finite element discretization is optimal with respect to the regularity of the solution. For example, in 3D space and for $\nu=\frac{1}{2}$ the realizations of a Gaussian random field with Matérn covariance are only almost surely Hölder continuous with any exponent $\beta<\frac{1}{2}$, see e.g. [12]. *

In the following, we denote the discrete linear algebra representations of the bilinear forms $m$, $b$, and $w$ with the matrices $M_{h}$, $W_{h}$, and $B_{h}$ where $M_{h}$ is the mass matrix for the space $\mathbf{R}_{h}$, $W_{h}$ is the mass matrix for the space $\Theta_{h}$, and $B_{h}$ stems from the divergence operator. Particular care must be taken for the linear algebra representation of the stochastic right hand side $F^{\mathcal{W}}(q_{h})$. This requires recalling the following two properties of Gaussian white noise defined on a domain $D$. For any set of test functions $\left\{q_{i}\in L^{2}(D),i=1,\dots,n\right\},$ the expectation and covariance measures are given by

[TABLE]

By taking $q_{i},q_{j}$ as piecewise constants so that $q_{i},q_{j}\in\Theta_{h}$, the second equation implies that the covariance measure over a region of the domain is equal to the area of that region [29].

As a result of these properties, the computation of the discrete stochastic linear functional amounts to computing

[TABLE]

where the coefficients of $\xi_{h}(\omega)$ in the finite element expansion form a random vector drawn from $\mathcal{N}(0,I)$. It should be noted that the mass matrix $W_{h}$ for the space $\Theta_{h}$ is diagonal, hence its square root can be computed cheaply. This is not the case in the original primal formulation presented in [29] since they employ piecewise linear continuous elements. The authors suggest then using mass lumping to make computation feasible at the cost of accuracy.

Then the discrete mixed finite element problem can be written as the linear system,

[TABLE]

where $f_{h}(\omega)\sim\mathcal{N}(0,W_{h}).$ An efficient iterative solution method of (13), required to scalably generate the desired Gaussian field realizations, is described in the following section.

It is worth noticing that the lowest order Raviart-Thomas spaces allow us to define a discrete gradient $\nabla_{h}:\Theta_{h}\mapsto\mathcal{R}_{h}$ such that the identity $w(\nabla_{h}\theta_{h},\mathbf{v}_{h})=-b(\mathbf{v}_{h},\theta_{h})$ holds for all $\theta_{h}\in\Theta_{h}$ and all $\mathbf{v}_{h}\in\mathcal{R}_{h}$ with zero normal trace on $\partial D$. Specifically, we define $\nabla_{h}=-M_{h}^{-1}B^{T}_{h}$. Then, the Schur complement of (13) with respect to $\theta_{h}$ can be viewed as a non-local discontinuous Galerkin (interior penalty) discretization of the original PDE (9), cf., [35]. Namely, we have

[TABLE]

The above formulation allows us to apply the theory in [29] to show convergence of $\theta_{h}$ to $\theta$, however a detailed analysis of the convergence of the mixed system (10) is outside the scope of this work.

2.3 SPDE Sampler Numerical Solution

The sparse large matrix $\mathcal{A}_{h}$ admits the following block UL decomposition

[TABLE]

where the blocks $A_{h}=M_{h}+\kappa^{-2}B_{h}^{T}W_{h}^{-1}B_{h}$ and $-\kappa^{-2}W_{h}^{-1}B_{h}$ are both sparse since $W_{h}$ is diagonal. Then solving the large sparse linear system  (13) amounts to first finding $\mathbf{u}_{h}$ such that

[TABLE]

and then to set

[TABLE]

It is worth noting that $A_{h}$ in (15) is symmetric positive definite and stems from the matrix representation of the weighted $H(\operatorname{div})$ inner product,

[TABLE]

thanks to the particular choice of the spaces $\mathbf{R}_{h}$ and $\Theta_{h}$, see e.g. [6, 7]. Then to solve the linear system (15), we use the conjugate gradient (CG) method preconditioned by the Auxiliary Space AMG preconditioner for $H(\operatorname{div})$ problems in [25], which ensures a mesh independent convergence and robustness with respect to the choice of the correlation length $\kappa$. In particular, for our numerical results, we use hypre’s HypreADS preconditioner [1].

This approach allows for scalable sampling, but the parametrization of $\theta_{h}(\omega)$ is mesh-dependent. This is not suitable for MLMC where the same realization of the Gaussian field needs to be computed at a fine and a coarse spatial resolution. In Section 3, we detail our proposed method that allows for hierarchical sampling based on a multilevel decomposition of the stochastic field.

2.4 Boundary artifacts and embedded domains

As noted in [29], the SPDE sampling method introduces errors on the boundary resulting in a larger marginal variance of the random field close to the boundary of the domain. These errors are due to the introduction of artificial boundary conditions where the SPDE is defined globally, yet must be discretized on a finite domain. To mitigate this issue we embed the original mesh into a larger mesh. Equation (8) is discretized on the larger domain $\overline{D}$ where a corresponding linear system of the form (13) is solved and the corresponding random field realization is restricted to the original domain $D$. A rule of thumb proposed in [28] suggests the boundary effect is negligible at a distance equal to the correlation length from the boundary. Figure 1 shows the sample marginal variance of 2000 Gaussian field samples with correlation length $0.01$, marginal variance $\sigma^{2}=1$ generated by the SPDE sampler. In Figure 1a, the samples were computed using the original circular mesh with diameter equal $0.2$, whereas in Figure 1c the SPDE given by (8) is solved on the original mesh embedded in a larger square domain with sides of length $0.3$ shown in Figure 1b. Thus, mesh embedding alleviates the issue of the error on the boundary.

3 Multilevel, Hierarchical Sampling

In this section we describe our proposed hierarchical sampling technique based on a multilevel decomposition of the stochastic field for sampling Gaussian fields based on solving the SPDE given by (9). We first briefly describe the agglomeration process of a fine grid into a sequence of coarser levels, introducing the necessary finite element spaces and inter-level operators that will be used. After introducing the multilevel structure of the stochastic field, the implementation details of the proposed method are described, followed by numerical results demonstrated the scalability of the method.

3.1 Multilevel Structure

Using methodology from element-based algebraic multigrid (AMGe), we are able to construct operator-dependent coarse spaces with guaranteed approximation properties on general, unstructured grids, see  [26, 27, 34] for further details.

Let $\mathcal{T}_{0}$ denote a fine grid discretization of the domain $D$. This fine grid is agglomerated into a hierarchy of $L$ coarser algebraic levels, $\{\mathcal{T}_{\ell}\}_{\ell=1}^{L}$ where $L$ denotes the coarsest grid. Agglomerates are formed by grouping together fine-grid elements. Based upon this hierarchy, we have the sequence of spaces $\mathbf{R}_{\ell},\Theta_{\ell}$ for $\ell=0,\dots,L$, that are the discrete analogues for $H(\operatorname{div},D)$ and $L^{2}(D)$ respectively on the discretizations $\{\mathcal{T}_{\ell}\}_{\ell=0}^{L}$. The space $\mathbf{R}_{0}$ is discretized by lowest order Raviart-Thomas finite elements and $\Theta_{0}$ is discretized by piecewise constant finite elements.

In addition, we will define the following interpolation operators from the coarse space $\Theta_{\ell+1}$ to fine space $\Theta_{\ell}$ where,

[TABLE]

We also define the operators from the coarse space $\mathbf{R}_{\ell+1}$ to the fine space $\mathbf{R}_{\ell}$ as

[TABLE]

Details on the construction and properties of the interpolation operators $P_{\theta}$ and $P_{\mathbf{u}}$ using AMGe techniques are given in [34, 27, 26]. Here we limit our discussion to observe that such operators reduce to the canonical interpolation operators of geometric multigrid when a nested hierarchy of uniformly refined meshes is given (in our case of constant PDE coefficients). In addition, we define

[TABLE]

as the block interpolation operator for the matrix $\mathcal{A}_{\ell}.$ Using this hierarchical notation, a discrete realization of a Gaussian random field on level $\ell$ corresponding to the $\mathcal{T}_{\ell}$ discretization will be denoted as $\theta_{\ell}$.

3.2 Hierarchical SPDE Sampler

We now propose a multilevel, decomposition of the random field and show that the coarse representation is a realization of the Gaussian random field on the coarse space $\Theta_{\ell+1}.$

Proposition 3.1**.**

The Gaussian random field $\theta_{\ell}(\omega)$ given by

[TABLE]

admits the following two-level decomposition:

[TABLE]

where $\theta_{\ell+1}(\omega)$ is a coarse representation of a Gaussian random field from the same distribution, and

[TABLE]

with the block expressions given by

[TABLE]

Proof 3.2**.**

We start with a block LU factorization of the operator in (20)

[TABLE]

where $I_{\ell}$ and $I_{\ell+1}$ are the identity matrices at level $\ell$ and $\ell+1$. Here the Schur complement operator $\mathcal{S}$ is defined as

[TABLE]

where, by definition of the prolongation operator, $\mathcal{A}_{\ell+1}$ is the coarse operator stemming from the discretization of the SPDE operator at level $\ell+1$. Recalling the definition of $U_{\ell}$, the above LU factorization implies

[TABLE]

To show that $\theta_{\ell+1}(\omega)$ is a realization of a Gaussian random field on the coarse space $\Theta_{\ell+1},$ we observe that $\theta_{\ell+1}(\omega)$ satisfies

[TABLE]

where $W_{\ell+1}=P_{\theta}^{T}W_{\ell}P_{\theta}$ is the mass matrix on the coarser level $\ell+1$. The random forcing term is defined as $\xi_{\ell+1}(\omega):=W_{\ell+1}^{-1/2}P_{\theta}^{T}W_{\ell}^{1/2}\xi_{\ell}(\omega)$. We note that $\xi_{\ell+1}(\omega)$ is a multivariate Gaussian vector with zero mean and covariance matrix

[TABLE]

where we have exploited the fact that $\xi_{\ell}(\omega)\sim\mathcal{N}(0,I)$.

*The two-level decomposition (19) then follows from (21)3 by noticing that $\theta_{\ell+1}(\omega)$ is the solution of the SPDE (9) discretized at level $\ell+1$ with random forcing term $\xi_{\ell+1}(\omega)\sim\mathcal{N}(0,I)$. *

Thus given $\xi_{\ell}(\omega)\sim\mathcal{N}(0,I)$ we are able to efficiently sample the Gaussian random field on both the fine and upscaled discretization using the SPDE sampler.

3.3 Hierarchical SPDE Sampler Numerical Solution

We now describe the solution procedure that we employ for the hierarchical SPDE sampler. Starting with an unstructured mesh, a coarse problem is constructed by grouping together fine-grid elements using the graph partitioner METIS [23] and the coarse finite element spaces are computed as described in detail in [27, 26]. The same procedure is applied recursively so that a nested hierarchy of agglomerated meshes and coarse spaces is constructed. We will present results on general unstructured meshes, however the technique immediately translates to the case of nested refined meshes generated with uniform refinement by choosing the canonical interpolators.

Given $\xi_{\ell}(\omega)$, we compute the realizations of the Gaussian field at levels $\ell$ and $\ell+1$ as follows. First we compute the sample for level $\ell+1$ by solving the saddle point system

[TABLE]

using the methodology described in Section 2.3. Then we compute the sample for level $\ell$ by iteratively solving (13) with $\mathcal{P}U_{\ell+1}$ as the initial guess.

We conclude this section with a remark on the independence of different realizations of the random field for a particular spatial resolution generated by the hierarchical SPDE sampler. Specifically, two realizations $\theta_{\ell}(\omega_{i})$ and $\theta_{\ell}(\omega_{j})$ are independent if and only if $\xi_{\ell}(\omega_{i})$ and $\xi_{\ell}(\omega_{j})$ are independent. This is of extreme importance for MLMC and requires the use of a quality random number generator. In our numerical experiments, we use Tina’s Random Number Generator Library [4] which is a pseudo-random number generator with dedicated support for parallel, distributed environments [8].

3.4 Numerical Results

The numerical results in this section are used to demonstrate the effectiveness of our SPDE sampling approach. For the simulations, an absolute stopping criteria of $10^{-12}$ and a relative stopping criteria of $10^{-6}$ is used for the linear solver. For the results presented in this and the following sections, we use the C++ finite element library MFEM [2] to assemble the discretized problems.

First we consider the performance of the sampler in two space dimensions. The physical domain $D=(0,1200)\times(0,2200)$ is embedded in a larger domain $\overline{D}=(-100,1300)\times(-100,2300)$ to decrease the effect of the artificial Neumann boundary condition, see Section 2.4. The computational domain is then discretized using a structured quadrilateral mesh with 294,400 elements in the physical domain. Then a hierarchy of unstructured agglomerated meshes is constructed using the graph partitioner METIS [23] with a coarsening ratio of 8 elements per agglomerate, see Figure 2. Note that, on coarse levels, agglomerated elements have irregular shapes and an arbitrary number of neighboring elements. Coarse faces are also not flat. Figure 3 shows a sequence of computed Gaussian random field samples $\theta(\mathbf{x},\omega)$ on different levels with correlation length $b=100$. Observe the similarity in the fields generated on the algebraically coarsened levels. Locations of the essential features of the fine grid sample are preserved, however contours are blurred and boundaries are spread out due to the irregular shape of the agglomerated elements.

We now investigate the performance of the hierarchical sampler under weak scaling, i.e. when the number of mesh elements is proportional to the number of processes. Structured coarsening is used to create the hierarchy of agglomerated meshes with a coarsening ratio of 4 elements per agglomerate, that is, the original mesh is uniformly refined to build the hierarchy of levels. The code was executed on Sierra at Lawrence Livermore National Laboratory consisting of a total of 1,944 nodes where each node has two 6-core Xeon EP X5660 Intel CPUs (2.8 Ghz), and 24GB of memory. We use 8 MPI processes per node. Figure 4 shows the number of MPI processes versus the average time to generate a realization of the Gaussian field for 1000 samples. The number of MPI processes ranges from 8 to 512 and size of the stochastic dimension on the fine grid ranges from 3.7 to 235 million. The proposed sampling method exhibits the optimal near linear scaling as the number of processes increases for all levels in the hierarchy.

As a second example, we consider the generation of large-scale 3D spacially correlated random fields. The computational domain is inspired by one from the SAIGUP benchmark, [16], and represents a realistic geometry of a shallow marine hydrocarbon reservoirs. The computational domain roughly covers an area of $3000\times 9000$ square meters, and the difference (in depth) between the highest and the lowest point in the domain is $930$ meters. The original mesh consists of $78720$ hexahedral elements and it is uniformly refined several times to build the hierarchy of levels. We note that in this example, we have not used embedding of the computational domain into a larger domain, since at this point we are only interested in the scalability of computing a single sample which is done by a scalable multilevel method.

Figure 5 shows a sequence of computed Gaussian random field samples $\theta(\mathbf{x},\omega)$ on different levels (the finer level is obtained by uniformly refining the original mesh twice). To model $\theta(\mathbf{x},\omega)$ we prescribe the exponential covariance function (5) with the correlation length $b=100$ meters and unitary marginal variance. Observe the similarity in the fields generated at the different levels. The weak scalability of the proposed sampling method is demonstrated in Figure 6. The number of MPI processes ranges from 48 to 3072 and the total number of degree of freedom (i.e. the number of unknowns in the mixed system (13)) on the fine grid ranges from 2.1 to 134 million. The proposed sampling method exhibits fairly scalable behavior as the number of processors increases.

4 Multilevel Monte Carlo Methods

In this section we briefly describe the Monte Carlo method, and then its extension to the multilevel Monte Carlo method. Finally, we present results with MLMC simulations for subsurface porous media flow using our hierarchical SPDE sampling method.

In Monte Carlo methods, one is interested in approximating the expected value of some quantity of interest, $Q(\omega)=\mathcal{B}\mathcal{U}(\mathbf{x},\omega)$ where $\mathcal{U}({\bf x},\omega)$ is the solution of a PDE with random input coefficients. In our case, $Q$ represents some functional of the solution $({\mathbf{q}},p)$ to (1)-(2). In general, the quantity of interest, $Q$, is inaccessible so an approximation, $Q_{h}$, is computed. The standard Monte Carlo estimator for the quantity of interest is then

[TABLE]

where $Q_{h}^{i}$ is the $i^{\text{th}}$ sample of $Q_{h}$ and $N$ is the number of (independent) samples.

The mean square error of the method is given by

[TABLE]

For the root mean square error (RMSE) to be below a prescribed tolerance, $\varepsilon$, both terms should be smaller than $\varepsilon^{2}/2$.

The first term is the estimator variance and the second term is the estimator bias. The estimator bias measures the discretization error and is controlled by the spatial resolution of the approximate solution. For sufficiently fine spatial discretizations, the estimator bias is small and reducing the mean square error amounts to reducing the estimator variance. The estimator variance is then reduced by increasing the number of samples, $N$. Thus, the standard MC method RMSE converges in $\mathcal{O}(1/\sqrt{N})$. This is favorable as the convergence rate in independent of the stochastic dimension of the problem, yet when high accuracy is necessary the number of samples required can be a prohibitive expense as the samples must be computed with a fine mesh. For this reason, standard MC methods are not a scalable approach for the forward propagation of uncertainties when the forward problem is a PDE, due to the high cost of computing samples with a fine spatial discretization.

This motivates the multilevel Monte Carlo (MLMC) method, see e.g. [14, 19], which is an effective variance reduction technique based on hierarchical sampling which aims to alleviate the burden of standard MC by computing samples on a hierarchy of spatial discretizations.

Consider the sequence of spatial discretizations $\ell=\{0,1,\ldots L\}$, with level $\ell=0$ denoting the finest spatial discretization, and level $\ell=L$ the coarsest used to approximate $Q_{\ell}$ where $Q_{0}=Q_{h}$.

The idea behind MLMC is to estimate the correction with respect to the next coarser level, $\mathbb{E}[Q_{\ell}-Q_{\ell+1}]$ rather than estimating $\mathbb{E}[Q_{\ell}]$ directly. Linearity of the expectation gives the following expression for $\mathbb{E}[Q_{h}]$ on the finest level,

[TABLE]

where we have defined $Y_{\ell}=Q_{\ell}-Q_{\ell+1}$ for $i=0,\ldots L-1$ and $Y_{L}=Q_{L}$.

Similarly, one computes the estimator for $Y_{\ell}$,

[TABLE]

and then the multilevel Monte Carlo estimator is defined as

[TABLE]

Consequently, the mean square error for the MLMC method becomes

[TABLE]

The three terms in the right hand side of (27) represent, respectively, the variance on the coarsest level, the variance of the correction with respect to the next coarser level, and lastly the discretization error. For a prescribed level of accuracy, the number of realizations at the coarsest level, $N_{L}$, still needs to be large, but samples are much cheaper to obtain on the coarser level, and the number of realizations required for levels $(\ell<L)$ given by $N_{\ell}$ is much smaller, since $\mathbb{V}[Q_{\ell}-Q_{\ell+1}]\rightarrow 0$ as $h_{\ell}\rightarrow 0$. Thus, fewer samples are needed for the finest, most computationally expensive level. To minimize the overall cost of the MLMC algorithm (e.g. the computational time to reach a desired MSE), the optimal number of samples of each level $N_{\ell}$ is given by

[TABLE]

where $C_{\ell}$ is the cost of computing one sample at level $\ell$. We refer to [14] for additional details.

4.1 Numerical Experiments

In this section we include standard results from MLMC computations using our proposed SPDE sampler. The forward model is the mixed Darcy equations given by

[TABLE]

with homogeneous Neumann boundary conditions $\mathbf{q}\cdot\mathbf{n}=0$ on $\Gamma_{N}$ and Dirichlet boundary conditions $p=p_{D}$ on $\Gamma_{D}$. Here, $\Gamma_{N}\in\partial D$, $\Gamma_{D}\in\partial D$ are a non overlapping partition of $\partial D$, and $\mathbf{n}$ denotes the unit normal vector to $\partial D$.

We use the mixed finite element method to discretize the model problem (28), specifically we choose the lowest order Raviart-Thomas element for the flux $\mathbf{q}$ and piecewise constant functions for the pressure $p$. Then, for each input realization $k(\mathbf{x},\omega)$, we write the resulting discretized saddle point problem as

[TABLE]

Here $\mathbf{f}_{h}$ stems from the discretization of the Dirichlet boundary condition $p=p_{D}$ on $\Gamma_{D}$. The large sparse indefinite linear system (29) is solved using MINRES preconditioned with a block diagonal preconditioner based on the pressure Schur Complement, namely the $L^{2}-H^{1}$ preconditioner described in [31]. Specifically, we consider the symmetric positive definite preconditioner

[TABLE]

where $H_{k,h}=\text{diag}(M_{k,h})$ and $\widetilde{\Sigma}_{k,h}=B_{h}H^{-1}_{k,h}B_{h}^{T}$.

It is well-known that $\mathcal{N}_{k,h}$ is a robust preconditioner for (29) as long as $M_{k,h}$ is not too anisotropic. In fact, as shown in [33], $\mathcal{N}_{k,h}$ is an optimal preconditioner for $\widetilde{\mathcal{A}}_{k,h}$,

[TABLE]

and $\widetilde{\mathcal{A}}_{k,h}$ is spectrally equivalent to $\mathcal{A}_{k,h}$, since the Raviart-Thomas finite element matrix $M_{k,h}$ is spectrally equivalent to its diagonal $H_{k,h}$.

In the computations, we use BoomerAMG from hypre [1] to precondition the Schur complement $\widetilde{\Sigma}_{k,h}$ which is explicitly available and sparse. It should be noted that the AMG preconditioner of $\widetilde{\Sigma}_{k,h}$ is recomputed for each input realization. For the simulations, an absolute stopping criteria of $10^{-12}$ and a relative stopping criteria of $10^{-6}$ is used for the linear solver.

4.1.1 Top Layer of SPE10 Dataset

First we show experiments incorporating data from the Tenth SPE Benchmark (SPE10) [3]. We consider a 2D slice of the dataset of dimension $1200\times 2200\,{\rm ft}^{2}$ divided into cells of size $20\times 10\,{\rm ft}^{2}$ resulting in a mesh with $60\times 220$ quadrilateral elements. The PDE coefficient on each slice is a scalar function. The original $60\times 220$ quadrilateral mesh corresponds to the coarsest one in our MLMC experiments. To produce the other (finer) levels, we uniformly refine the initial 2D mesh several times.

We have $D=(0,1200)\times(0,2200)$ and assume the random conductivity coefficient $k(\mathbf{x},\omega)$ is modeled as a log-normal random field. A realization of $k(\mathbf{x},\omega)$ is generated by computing the exponential of a realization of a Gaussian random field. In particular, we assume that the mean of the Gaussian random field is the logarithm of the top horizontal slice from the SPE10 dataset.

Figure 7c shows a particular realization of the random conductivity coefficient $k(\mathbf{x},\omega)$ modeled as a log-normal field where $k(\mathbf{x},\omega)=\operatorname{exp}\left[{\operatorname{log}[k_{SPE10slice}(\mathbf{x})]+\theta(\omega)}\right]$ where $\theta(\omega)$ is a realization of the Gaussian random field generated with our sampler shown in Figure 7b and $k_{SPE10slice}(\mathbf{x})$ is shown in Figure 7a.

We solve (28) with the following boundary conditions:

[TABLE]

The quantity of interest we use is the expected value of the effective permeability, that is the (horizontal) flux through the “outflow” part of the boundary, defined as

[TABLE]

Figure 8 contains four subplots relating to the multilevel estimator and performance of the multilevel Monte Carlo method with hierarchical, SPDE sampling. The target mean square error is $\epsilon^{2}=8.73e$-5 for Figures 8a-8c. The first figure, Figure 8a, displays the multilevel estimator where the blue line with circles represents the expectation at each level $\mathbb{E}[Q_{\ell}]$ and the green dashed line represents the expectation of the difference in levels, $\mathbb{E}[Q_{\ell}-Q_{\ell+1}]$. Figure 8b illustrates the multilevel variance reduction. This plot contains two lines, the blue line represents the variance of the particular level, whereas the green dashed line represents the variance of the difference in levels. The plot shows the effectiveness of the MLMC method at reducing the variance as the number of unknowns increases. The average sampling time to generate the required Gaussian field realizations and solve the forward model for each level is shown in Figure 8c. This plot indicates near optimal scaling of the MLMC method with the proposed hierarchical sampler. Figure 8d shows the number of samples required at each level of the MLMC method for different prescribed mean square error tolerances. The plot clearly shows that more samples are generated on the coarse levels (fewer degrees of freedom) than on the finest levels (many degrees of freedom). This merely confirms the MLMC theory with our proposed hierarchical sampling technique.

4.1.2 Unit Cube

In this section we present similar experiments for estimating the expectation of the effective permeability, but for the unit cube domain $D=(0,1)\times(0,1)\times(0,1)$.

We solve (28) with the following boundary conditions:

[TABLE]

The quantity of interest is the expected value of the effective permeability defined in (30). The original mesh consists of 64 hexahedral elements and is uniformly refined several times to build the hierarchy of levels. We examine the performance of the multilevel estimator for the unit cube in Figure 9 which contains four subplots. The first figure, Figure 9a, displays the multilevel estimator where the blue line with circles represents the expectation at each level $\mathbb{E}[Q_{\ell}]$ and the green dashed line represents the expectation of the difference in levels, $\mathbb{E}[Q_{\ell}-Q_{\ell+1}]$.

Figure 9b illustrates the multilevel variance reduction of the method, while Figure 9c shows the average sampling time required to generate the Gaussian random field realizations and solve the forward model problem for each level. The method with hierarchical, SPDE sampling exhibits near optimal scaling for the 3D problem formulation. The number of samples required at each level of the MLMC method for different prescribed mean square error tolerances is shown in 9d.

For both of the examined computational domains, the hierarchical, SPDE sampler yields the expected results for MLMC variance reduction and displays the desired scaling properties for the possibility of large-scale MLMC simulations.

5 Conclusions

Multilevel Monte Carlos simulations for PDEs with uncertain input coefficients employ a hierarchy of spatial resolutions as a variance reduction technique for the approximation of expected quantities of interest. A key component in the multilevel Monte Carlo method is the ability to generate samples of a random field at different spatial resolutions. The Karhunen-Loève expansion provides a parametrization independent of the spatial discretization. However, both the computation and the memory requirements become infeasible at large-scale as the expansion requires the ability to compute and store eigenpairs of a large, dense covariance matrix. We suggest a sampling method based on the solution of a particular stochastic PDE. This method is highly scalable, but the parametrization is mesh dependent. We have proposed a multilevel decomposition of the stochastic field to allow for scalable, hierarchical stochastic PDEs samplers. Numerical results are provided that suggest the method possesses the desired scalability as the method leverages existing scalable solvers. We also have applied the new sampling technique to MLMC simulations of subsurface flow problems with over 10 million parameters in the stochastic dimension.