Introduction

Many mathematical models describing industrial problems are subject to uncertainty due to some limitation. Incorporating the uncertainty typically leads to more accurate representation of the world. However it is in the cost of solving a statistical problem, for example one can consider saturated flow in subsurface, or heat condition in Metal Matrix Composites. Such stochastic models require enormous computational effort, thereby requiring new fast algorithms to facilitate that need. In this paper we consider a scalar elliptic SPDE describing single phase flow throughout heterogeneous porous media. Although the approach here is not limited to this problem, we aim at computing the mean flux throughout saturated porous media with prescribed pressure drop and known distribution of the random coefficients.

One of the preferred and powerful methods for solving SPDE is Multilevel Monte Carlo algorithm (MLMC). The algorithm exploits a combination of fewer expensive computations with a plenty of cheap ones to compute expected values at significantly lower computational cost then standard Monte Carlo. A key component is the selection of the different levels of combination. Many different approaches exist. In [1] authors use the number of therms of Karhunen-Loewe expansion in order to define coarser levels. In [4], similarly to our approach, the coarser levels are defined on coarser grids via averaging of the coefficients in the PDE. Here we construct the levels by renormalization. For details we reefer to [5,3]. For generation of permeability (random field) we use circulant embedding algorithm [7].

Due to the rapidly advancing area of computer science, methods based on Monte Carlo sampling are of great interest. Such methods are suitable for parallelization and can compute the expected value in reasonable time. For realistic simulations usually a HPC implementation is need. Here we investigate the possible benefits of generating permeability on GPUs and study the reduction of time in overall computation of MLMC.

Model problem

In order to test the overall performance of MLMC in combination with CUDA generated permeability field, we consider simple model problem in a unit cube domain, steady state single phase flow in random porous media. This problem, illustrates well the challenges in solving stochastic PDEs.

−∇ • [k(x, ω)∇p(x, ω)] = 0 for x ∈ D = (0, 1) d , ω ∈ Ω ( 1 )

Subject to boundary conditions:

p x=0 = 1 p x=1 = 0 ∂ n p = 0 on other boundaries,(2)

with dimension d = {2, 3}, pressure p, scalar permeability k, and random vector ω. The quantify of interest is the mean (expected value) of the total flux through the unit cube:

Q(x, ω) := ˆx=0 k(x, ω)∂ n p(x, ω)dx (3)

Both the coefficient k(x, ω) and the solution p(x, ω) are subject to uncertainty, characterized by the random vector ω in a properly defined space.

Solving this equation can be broken into three sub-problems: generating permeability (random field), solving the deterministic problem and reducing the variance with MLMC method. We briefly discuss each of them.

Generating permeability field

Generating permeability fields is essential problem in solving the SPDE. Here we consider a practical covariance proposed by [8]:

C(x, y) = σ 2 exp(−||x − y|| p /λ), p = 2(4)

Where || • || p denotes the l p norm. in R d . which satisfies:

E[K(x, .)] = 0, E[K(x, .), K(y, .)] = C(x − y) = C(y − x) for x, y ∈ D and K(x, ω) = log(k(x, ω))

Several approaches has been developed of generating random permeability fields applicable to flow simulations. Here we use an algorithm based on forward and inverse Fourier transform over circulant covariance matrix to generate permeability. Realization of such field is governed by two parameters: standard deviation σ and correlation length λ. More details can be found in the papers [2] and [7].

Solving the deterministic problem

Literature provides different numerical schemes for solving PDEs. In [2], an Multi-scale Finite element methods is used. Here for solving the elliptic PDEs corresponding to each realization of permeability field, we use finite volume method on a cell centered grid. This method is mass conservative. More details can be found in [1,3].

Variance reduction

We shortly recall the idea proposed in [1]. Let {M l : l = 0 . . . L} ∈ N be increasing sequence of numbers called levels, with corresponding quantities {Q M l } L l=0 , and s ≥ 2 be coarsening factor, such that we have

M l = sM l−1 , for l = 0 . . . L. Defining Y l = Q M l − Q M l−1 and setting Y 0 = Q M 0 ,

we can use telescopic sum to write the following identitity for the expected value

E E[Q M ] = E[Q M 0 ] + L l=1 E[Q M l − Q M l−1 ] = L l=0 E[Y l ]( 5 )

The expectation on the finest level is equal to the expectation on the coarsest level plus sum of corrections, i.e. differences in the expectations on each pair of consecutive levels. The terms in (5) are approximated using standard MC independent estimators, with N l samples. For the mean square error we have:

e( Q ML M,N ) 2 = E[( Q ML M,N − E[Q]) 2 ] = L l=0 N −1 l V [Y l ] + (E[Q M − E[Q]) 2 (6)

Our goal is to have:

e( Q ML M,N ) 2 = L l=0 N −1 l V [Y l ] + (E[Q M − E[Q]) 2 ≤ 2ε 2(7)

Denote

v l = V [Y l ],

and let t l be the mean time computing difference Y l once, and T = L l=0 n l t l be the total time for the computation. Minimizing T under the above constraint, with Lagrangian multipliers and turning it to integer value gives us:

n l = α (v l /t l ) with Lagrangian multiplier α = 1 2 L l=0 (v l /t l )(8)

To define the levels, the resolution of the discretization is used, such that the number of square cells in D = 2 and cubic cells in D = 3 are exact power of 2. Thus on the finer level we have 4 times more cells than on the coarser level for D = 2, and 8 times more for D = 3. To approximate the random field, we employ heuristic technique, witch combines each 4 cells in 2D into one by combination of simple arithmetic, geometric and harmonic average. For details we refer to [3].

Simulation results

The computational time used by the setup of the generation of 100 permeability fields with σ = 2.0 and λ = 0.2 is shown on figure 2. Our underling grid is in 2D with size, approximately 1.7 * 10 7 . As it can be observed, the speed up that is achieved using only CPU cores for computation at all of the steps is significant up to around 8 cores. At this point it is approximately up to 6.755 times faster than calculation on a single core. Using more than 8 cores efficiency gradually decreases and with 24 cores the speed up achieved is approximately 15 times. The GPU generation is much faster: using single GPU, generation time is comparable to 11 CPU cores, and on 2 GPUs -approximately 22 CPU cores. This is expected since permeability generation is mainly a forward Fourier transform over a matrix, followed by multiplication with random numbers for each element and then inverse Fourier transform. Table 1 and 2 show the comparison of 3 level MLMC implementation using only CPU versus implementation, that uses GPU for permeability generation. The presented results are averaged over 10 runs of the algorithm. The fine grid is of size 2 10 * 2 10 and permeability generation parameters are σ = 1.5, λ = 0.1. Each problem is computed by own CPU core and all use shared GPU, thus in a given moment of the execution 12 or 24 concurrent calls to the GPUs may exist. This is more taxing situation for the GPU compared to a distribution where a group of processes solves a single problem. In table 1 one can observe that MLMC calculated on single GPU and 12 CPU cores is faster than same algorithm using only CPU cores and the execution of the generation step on the different levels is approximately 2 times faster while the execution of the same MLMC implementation on single GPU and 24 CPU cores shows that the overall time is slower than execution only on the same number of CPU cores. The performance is significantly improved with introduction of a second GPU. The generation times are notably lower than the case with a single GPU. The implementation of the algorithm, used for that test is written in C++. The method used for solving the PDE conjugate gradient preconditioned with AMG if provided by DUNE library [9]. For implementation of circulant embedding algorithm the fftw library is used for generation on CPU, and cufft library provided by NVIDIA is used for generation on GPU. The Multilevel Monte Carlo algorithm is implemented with pure MPI. All tests are performed on HybriLIT-education and testing cluster, Dubna Russia, on a GPU node with two NVIDIA TESLA K80 GPUs and Intel Xeon E5-2695 v2 processors.

#GPU E[Q] Time[s] RMS t gen l 0 [s] t

Conclusions

Generating random filed using modern GPU accelerated computing is promising possibility, witch can decrease the computational time of the generation step of the MLMC algorithm if the task is small enough to fit in the GPU memory. When the task is larger, however, the execution time using GPU can be larger than using only CPU cores. The idea of GPU acceleration is fairly new in the area of scientific computing and there are not many libraries implemented, using GPU accelerated approach. Further more, HPC systems with GPU enabled nodes are rather expensive and the number of large scientific computing clusters with GPU capabilities is relatively small. On the other hand, in the recent years, some large commercial clusters for GPU computing were built, inspired by the block-chain and neural networks trend. A further investigation on solving large scientific problems on such clusters is an interesting study.