A full multigrid multilevel Monte Carlo method for the single phase subsurface flow with random coefficients

The subsurface flow is usually subject to uncertain porous media structures. In most cases, however, we only have partial knowledge about the porous media properties. A common approach is to model the uncertain parameters of porous media as random fields, then the statistical moments (e.g. expectation) of the Quantity of Interest(QoI) can be evaluated by the Monte Carlo method. In this study, we develop a full multigrid-multilevel Monte Carlo (FMG-MLMC) method to speed up the evaluation of random parameters effects on single-phase porous flows. In general, MLMC method applies a series of discretization with increasing resolution and computes the QoI on each of them. The effective variance reduction is the success of the method. We exploit the similar hierarchies of MLMC and multigrid methods and obtain the solution on coarse mesh Q^c_l as a byproduct of the full multigrid solution on fine mesh Q^f_l on each level l. In the cases considered in this work, the computational saving due to the coarse mesh samples saving is 20% asymptotically. Besides, a comparison of Monte Carlo and Quasi-Monte Carlo (QMC) methods reveals a smaller estimator variance and a faster convergence rate of the latter approach in this study.



There are no comments yet.


page 1

page 2

page 3

page 4


On the effective dimension and multilevel Monte Carlo

I consider the problem of integrating a function f over the d-dimensiona...

Enhanced Multi-Index Monte Carlo by means of Multiple Semi-Coarsened Multigrid for Anisotropic Diffusion Problems

In many models used in engineering and science, material properties are ...

Multilevel Monte Carlo Simulations of Composite Structures with Uncertain Manufacturing Defects

By adopting a Multilevel Monte Carlo (MLMC) framework, we show that only...

Solution of the 3D density-driven groundwater flow problem with uncertain porosity and permeability

As groundwater is an essential nutrition and irrigation resource, its po...

Probabilistic forecast of multiphase transport under viscous and buoyancy forces in heterogeneous porous media

In this study, we develop a probabilistic approach to map the parametric...

Constrained Ensemble Langevin Monte Carlo

The classical Langevin Monte Carlo method looks for i.i.d. samples from ...

Frozen Gaussian Sampling: A Mesh-free Monte Carlo Method For Approximating Semiclassical Schrödinger Equations

In this paper, we develop a Monte Carlo algorithm named the Frozen Gauss...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Flow and transport are the most fundamental phenomena in subsurface porous media associated with various physical processes, e.g., oil and gas flow in petroleum reservoir Terry et al. (2015), CO sequestration Zhang et al. (2015), water pollution dispersion Bear and Cheng (2010), etc. The numerical simulation and analyses of flow and transport in subsurface porous media are highly demanded in practical engineering and mechanism studies. However, the simulation results are always subject to the influence of uncertainties, mainly stemming from the inherent spatial heterogeneity of media properties caused by complex geological processes Boschan and Noetinger (2012). It has been widely recognized that in natural subsurface porous media, most properties, such as permeability, porosity, etc., exhibit an uneven spatial distribution. For example, the hydraulic conductivity can span several orders of magnitude in an aquifer or reservoir. How to quantificationally identify the influence of uncertainties of porous media properties on the flow and transport behaviors in subsurface physical processes has been a research hot spot in recent years.

Therefore, the uncertainty quantification is an essential task in the simulation of practical subsurface flows where porous media properties that unknown or partially known are taken as the input parameters. A possible way to deal with uncertainties of subsurface porous media is to treat porous media properties as random fields, then perform the stochastic simulation on the subsurface flow governing equations with random coefficients to evaluate the quality of interest (QoI). Among commonly-used stochastic simulation methods, e.g. Monte Carlo(MC) method, stochastic finite element method, stochastic collocation method, the MC method demonstrates apparent advantages such as it is a non-intrusive approach that only the realization of coefficients is needed while the original model code remains unchanged, and it is more easily to be implemented. In the standard MC method, the computer-generated (pseudo) random points are used, and in many cases, the computational efficiency is always unsatisfied for large-scale problems. The Quasi-Monte Carlo (QMC) method improve the demerit of MC by using deterministic quasi-random points. These points exhibit lower discrepancy and distribute more uniformly in the probability space. Moreover, to reduce the sample variance and further improve the computational efficiency, the multilevel Monte Carlo (MLMC) method was proposed and developed by Heinrich

Heinrich (2001) and Giles Giles (2008). It applies the control variates technique that a series of discretization is adopted with increasing resolution and computes the QoI on each of them, the success of which lies in the effective variance reduction sequentially.

It should be mentioned that in the particular case of subsurface flow with random coefficients, the problem is further aggravated where very detailed geological models are needed (a large number of cells) for an accurate description of the flow. To further alleviate the computational burden connected to the evaluation of random parameter effects on subsurface flow using the MLMC method, in this study, we exploit the similar hierarchies of MLMC and multigrid methods and proposed a full multigrid multilevel (quasi-) Monte Carlo (FMG-MLQMC) approach. In this proposed method, the solution on coarse mesh can be obtained as a byproduct of the full multigrid solution on fine mesh on each mesh level , instead of directly solving the equations on the coarse mesh as the standard MLMC does. The proposed FMG-MLQMC method saves the computation of the . There have been works coupling the multigrid solver with the multilevel framework, see Kumar et al. (2017); Robbe et al. (2018) for example. However, the FMG-MLQMC method we proposed saves the computational cost without modifying the MLMC framework. We exploit the implementation method for upscaling the random coefficient from fine mesh to neighboring coarse mesh. Although in this study we only focus on the simple single-phase subsurface flow with random coefficients, the proposed approach can be applied and extended naturally to multiphase flow and transfer in porous media and any other flow and transport problems associated with uncertainty effect.

The rest of the paper is organized as follows: In Section 2, we give a brief description of the single-phase subsurface flow and then introduce the proposed full multigrid multilevel (quasi-) Monte Carlo method in detail. The methodology on the upscaling method of random coefficients from the fine mesh to the coarse mesh, which preserves the random structure, is presented as well. In Section 3, we verify the effectiveness (a smaller estimator variance and faster convergence rate) of the presented method by comparing with standard MLMC method in two numerical experiments. Finally, in Section 4, we report the concluding remarks of this work along with a brief discussion of future directions.

2 Algorithms

2.1 Model problem and MLMC method

In this work, we consider the following elliptic problem,


where is the random, spatial-varying coefficient, is the computational domain, is a sample from the probability triple . and are Dirichlet and Neumann boundaries respectively. In single-phase flow context, Eq.1 corresponds the steady-state situation, when and prescribe the pressure and velocity of the fluid at the boundary, then the solution depicts the pressure in the domain .

In this work, we address the random elliptic problem using the multi-level algorithm. Basically, the MLMC method employs a series of control variates, which are often the discretized models with increasing resolution levels. Here, we associate each level with one mesh with given resolution. The approximations of quantity of interest(QoI) on these levels are denoted as , see Figure 1.

Figure 1: Multilevel Monte Carlo

We are interested in the approximation on the finest level . The MLMC method not only computes the solution on level itself, but also calculates the solutions on all the preceding meshes. The expectation of such quantity can be expressed by the following telescoping formula,


We can approximate each expectation using the Monte Carlo approach as follows,


Here, we associate level with the -th term in the telescoping formula (2). Notice that on each level , we use the same samples to calculate and . Then and are likely to correlate well, and the variance of will be small, see Eq.4.


where we denote as the variance of the on level .

Also notice that on different levels, independent samples are used so that the variance of the multilevel estimator is the summation of the variance on each level.


where we let . If we write and , then


be an unbiased estimator for


then the multilevel estimator becomes,


2.2 MLMC Complexity Theory

Let denote a quantity of interest, and denote the corresponding numerical approximation on -th mesh. If we assume that the weak error and the level variance decreases exponentially while the cost per sample on each level increases exponentially, there exist positive constants , and satisfying the following,


where is the cost per sample on level . With the mean square error less than a threshold,


the total computational cost satisfies,

as .

2.3 MLMC Algorithm

This subsection gives a MLMC algorithm initially proposed by M.GilesGiles (2008).

Start with , set the initial number of samples on level 0,1,2
while extra samples need to be evaluated do
       evaluate and , for
       update estimates for ,
       update optimal , compute the number of extra samples
       if  then
             and initialize ;
       end if
end while
Algorithm 1 MLMC Algorithm

In Algorithm 1, the variances are approximated by the sample variances on the run. The weak error is approximated by Richardson extrapolation .

And here we consider a equal split of the estimator variance and the approximation error, i.e.,


However, it is possible to determine the split factor in an optimal wayCollier et al. (2015).

The optimal number of samples can be obtained by solving a constrained optimization : minimizing the total computational cost subject to the constraint Eq.10.

MLMC algorithm will work under the following three conditions.


The sequence converges. Otherwise, the telescoping equation (2) does not yield a converging result.


and are estimated using the same underlying random sample in equation (3), and are thus well correlated. In this case, the estimator variance is significantly reduced.


The telescoping sum (2) introduces no bias error. Notice that in the telescoping equation, the term for appears twice. However, the two may be evaluated differently. If we denote the -th term in the telescoping equation and by and , respectively, which denote the fine mesh solution and coarse mesh solution on level , then the condition needs to be satisfied in order to introduce no bias error in equation (2). The expectation of fine mesh solution on level should be the same as that of the coarse mesh solution on level .

2.4 MLQMC Algorithm

The QMC method can solve integration problems as well. In contrast to the MC method, the QMC method replaces random points with deterministic points. Figure 2 gives an example of Monte Carlo points, lattice rule, and Sobol’ sequence.

(a) Monte Carlo
(b) Lattice rule
(c) Sobol’ sequence
Figure 2: An example of Monte Carlo points, Lattice rule and Sobol’ sequence in domain.

The QMC approximation is given by,

Notice that the points are deterministic. However, deterministic points yield biased estimates. In this work, we use the randomized QMC: random shift and digital scramble for lattice rule and Sobol’ sequence respectively. The interested readers may refer to Dick et al. (2013) for the above mentioned two randomization techniques.

There have been numerous studies combining the MLMC and QMC methodsGiles and Waterhouse (2009); Giles et al. (2016); Kuo et al. (2017). We follow the multilevel quasi-Monte Carlo (MLQMC) settings from these works and list the algorithm here.

Start with , set initial number of samples on level 0,1,2
while extra samples need to be evaluated do
       evaluate and , for
       update estimates for , and compute
       if  then
             select level such that and double
             if  then
                   and initialize
             end if
       end if
end while
Algorithm 2 MLQMC Algorithm

2.5 Multigrid

The multigrid method was originally introduced to solve elliptic boundary-value problems efficiently. It has since been developed to solve either linear or non-linear systems. Multigrid methods compute the solution on a sequence of grids. Figure. 3 gives an illustration of the full multigrid scheme.

We observe that, when the full multigrid solver is applied to the MLMC problem, based on the same level hierarchies, the solution on the coarse mesh can be obtained as a byproduct of the multigrid solution on fine mesh on each level . Thus, in our proposed FMG-MLMC method, we have saved the computation for . Also notice that at the red circles in Figure. 3 the solutions are exact, since they are the end point of each V-cycle.

Figure 3: An illustration of Full-Multigrid-Multilevel Monte Carlo method.

Recall that and are evaluated by the same underlying random sample. In level , we denote and as the coefficients of the fine and coarse models, respectively. Also recall the consistency condition, since we use the same numerical solver for and , we only require that and follow the same distribution.

can be generated by the matrix decomposition method, KL-expansion method or other random field generation methods, see Liu et al. (2019) for example. A way to generate is to coarsen . In order to prevent bias error, should satisfy the same distribution law as . Figure 4 shows a way of coarsening.

Figure 4: Coarsening

In this scheme, the value of the coefficient in the coarse grid is selected to be the corresponding block in the fine grid. We denote the blocks in fine grid and coarse grid by and respectively, then we have the following,

Here we give an algorithm to describe our proposed FMG-MLMC method.

Start with , set the initial number of samples on level 0,1,2
while extra samples need to be evaluated do
       coarsen to obtain
       use the multigrid solver to compute realizations , and obtain as a byproduct, for
       update estimates for ,
       update optimal , compute number of extra samples
       if  then
             and initialize ;
       end if
end while
Algorithm 3 FMG-MLMC Algorithm

Notice that in our scheme, no correlation is introduced into each levels. Thus the equation (5) still holds true.

3 Numerical Validation

3.1 Problem Statement

Recall the elliptic problem 1 in Section 2, now we focus the physical domain . In this work, we consider cases in two different boundary conditions and quantities of interest, as listed in Table 1. Case I is of the Dirichlet boundary type, with pointwise output quantity. Case II is of the mixed Dirichlet-Neumann boundary condition, whose output is the outflow at the east boundary.

Case Boundary Condition QoI
1 , , ,
2 , , ,
Table 1: Case Settings

(1) Discretization
The governing equation (1) is discretized by the finite-volume method on rectangular grids. On level

the degree of freedom is

(2) Random Fields
We choose the Matérn covariance function


where is the Euclidean distance of two points, is the correlation length, and controls the smoothness of the field.

The parameters of Matérn covariance for four different random fields are given in Table 2.

Random Field Parameters
Table 2: Random Field Parameters Settings

The random fields are generated using the KL-expansion method. The truncation term is determined when of the variability is captured, meaning that

where the summation of all eigenvalues satisfies the following

Readers of interest can see Ernst et al. (2009) for examples. is the random field region.
(3) QMC Method
In QMC method, the Lattice rule points are generated using the software from Kuo and Nuyens (2016), Sobol’ matrices from Joe and Kuo (2008)

are used to generate Sobol’ sequences. In both cases, 24 randomizations are applied. The confidence intervals are obtained by 10 sets of randomization.

3.2 Numerical results

In this subsection we present the numerical results of the two cases (Tab. 1) with four random field settings (Tab. 2). In each figure, the first and second row corresponds the cases and respectively, while the first column and second column corresponds the cases and . We will first give the results of the first case in the following.

The simulation starts by estimating the asymptotic rates in the assumptions (7). Figs. 5, 6 and 7 plot the variance of the QoI and against level . By comparison, the QMC method reduces in the variance not in the asymptotic variance convergence rate, but in the y-axis offsets.

Figure 5: Variance of and for 4 random fields

For the results of MLQMC-Lattice, the variance test (Fig. 6) is presented.

Figure 6: Variance of and for 4 random fields

For the results of MLQMC-Sobol’, the variance (Fig. 7) is presented.

Figure 7: Variance of and for 4 random fields

The QMC method does not affect the expectation value nor the computational cost. Here we skip the comparison of and , but present the final computational cost of the three methods, given . The results in 4 random fields are presented in Fig. 8.

Figure 8: Computational complexity of the three methods; the asymptotic rates are marked on the plot

For case II, we test the QMC convergence rate and then present the computational complexity.

Next the results of MLQMC-Lattice, the convergence test (Fig. 9) variance of level estimator against in each level. The offset between the lines reveals the variance decreases with the levels. The comparison shows Sobol’ sequence’s advantage in decreasing variance. The random parameter settings have small impact on the variances in this case.

Figure 9: Variance test of MLQMC-Lattice for all levels. The variance of the estimator is plotted as a function of the number of samples . The convergence rates are marked on the plot

Next, the results of MLQMC-Sobol’, the convergence test (Fig. 10)

Figure 10: Variance test of MLQMC-Sobol for all levels. The variance of the estimator is plotted as a function of the number of samples . The convergence rates are marked on the plot

Finally, we compare the computational cost of the three methods, given . The results in 4 random fields are presented in Fig. 11.

Figure 11: Computational complexity of the three methods

4 Conclusions

In this work, we combined the MLMC with a full multigrid method. We saved the computation for the coarse grid solution on each level without modifying MLMC hierarchy and introducing correlation between different levels. We applied the consistent coarsening approach such that no bias error was introduced in the telescoping formula (2).

We tested our FMG-MLQMC algorithm on 2-D elliptic PDE with random coefficients for two different types of boundary condition settings and QoIs. The random coefficients were modelled as lognormal random fields with the Matérn covariance function with various parameter settings. We observed that quasi-Monte Carlo approaches have better performance on smoother random fields and problems with more regularity. Also, the comparison of Monte Carlo and quasi-Monte Carlo methods (including Lattice rule and Sobol’ sequence) revealed that QMC outperforms MC due to a smaller estimator variance, and Sobol’ sequence performs slightly better than Lattice rule.

One future work could be the substitution of the geometric multigrid solver with the algebraic multigrid (AMG) solver. In the AMG scheme, the grids are not associated with physical meshes, rather, the grids are fully determined by the matrix entries algebraically.

Another work could be to extend the elliptic model to more sophisticated models, such as two-phase porous flow. The efficient sampling and fast simulation of multi-phase subsurface flow under heterogeneous media could produce practical values.

Further work could extend the multilevel model to a multiscale model, and multiscale meshes rather than geometric meshes would be used. In this case, one level could correspond to one scale, and sampling would be performed on each scale. The literature on multiscale modeling can be found in Jenny et al. (2003); Efendiev and Hou (2009).


The authors gratefully acknowledge the support from the National Natural Science Foundation of China (Nos. 51874262, 51904031) and the Research Funding from King Abdullah University of Science and Technology (KAUST) through the grants BAS/1/1351-01-01.


  • J. Bear and A. H. Cheng (2010) Modeling groundwater flow and contaminant transport. Vol. 23, Springer Science & Business Media. Cited by: §1.
  • A. Boschan and B. Noetinger (2012) Scale dependence of effective hydraulic conductivity distributions in 3d heterogeneous media: a numerical study. Transport in porous media 94 (1), pp. 101–121. Cited by: §1.
  • N. Collier, A. Haji-Ali, F. Nobile, E. von Schwerin, and R. Tempone (2015) A continuation multilevel monte carlo algorithm. BIT Numerical Mathematics 55 (2), pp. 399–432. External Links: ISSN 1572-9125 Cited by: §2.3.
  • J. Dick, F. Y. Kuo, and I. H. Sloan (2013) High-dimensional integration: the quasi-monte carlo way. Acta Numerica 22, pp. 133–288. Cited by: §2.4.
  • Y. Efendiev and T. Y. Hou (2009) Multiscale finite element methods: theory and applications. Vol. 4, Springer Science & Business Media. Cited by: §4.
  • O. G. Ernst, C. E. Powell, D. J. Silvester, and E. Ullmann (2009) Efficient solvers for a linear stochastic galerkin mixed formulation of diffusion problems with random data. SIAM Journal on Scientific Computing 31 (2), pp. 1424–1447. Cited by: §3.1.
  • M. B. Giles, F. Y. Kuo, and I. H. Sloan (2016) Combining sparse grids, multilevel mc and qmc for elliptic pdes with random coefficients. In International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing, pp. 265–281. Cited by: §2.4.
  • M. B. Giles and B. J. Waterhouse (2009) Multilevel quasi-monte carlo path simulation. Advanced Financial Modelling, Radon Series on Computational and Applied Mathematics 8, pp. 165–181. Cited by: §2.4.
  • M. B. Giles (2008) Multilevel monte carlo path simulation. Operations Research 56 (3), pp. 607–617. Cited by: §1, §2.3.
  • S. Heinrich (2001) Multilevel monte carlo methods. In LSSC, Cited by: §1.
  • P. Jenny, S. Lee, and H. A. Tchelepi (2003) Multi-scale finite-volume method for elliptic problems in subsurface flow simulation. Journal of Computational Physics 187 (1), pp. 47–67. Cited by: §4.
  • S. Joe and F. Y. Kuo (2008) Constructing sobol sequences with better two-dimensional projections. SIAM Journal on Scientific Computing 30 (5), pp. 2635–2654. Cited by: §3.1.
  • P. Kumar, C. W. Oosterlee, and R. P. Dwight (2017) A multigrid multilevel monte carlo method using high-order finite-volume scheme for lognormal diffusion problems. International Journal for Uncertainty Quantification 7 (1). Cited by: §1.
  • F. Kuo, R. Scheichl, C. Schwab, I. Sloan, and E. Ullmann (2017) Multilevel quasi-monte carlo methods for lognormal diffusion problems. Mathematics of Computation 86 (308), pp. 2827–2860. Cited by: §2.4.
  • F. Y. Kuo and D. Nuyens (2016) Application of quasi-monte carlo methods to elliptic pdes with random diffusion coefficients: a survey of analysis and implementation. Foundations of Computational Mathematics 16 (6), pp. 1631–1696. Cited by: §3.1.
  • Y. Liu, J. Li, S. Sun, and B. Yu (2019) Advances in gaussian random field generation: a review. Computational Geosciences. External Links: ISSN 1573-1499 Cited by: §2.5.
  • P. Robbe, D. Nuyens, and S. Vandewalle (2018) Recycling samples in the multigrid multilevel (quasi-) monte carlo method. arXiv preprint arXiv:1806.05619. Cited by: §1.
  • R. E. Terry, J. B. Rogers, and B. C. Craft (2015) Applied petroleum reservoir engineering. Pearson Education. Cited by: §1.
  • R. Zhang, P. H. Winterfeld, X. Yin, Y. Xiong, and Y. Wu (2015) Sequentially coupled thmc model for co2 geological sequestration into a 2d heterogeneous saline aquifer. Journal of Natural Gas Science and Engineering 27, pp. 579–615. Cited by: §1.