Monte Carlo sampling methods is a broad class of computational algorithms, that rely on repeated random sampling to obtain numerical results. The idea of such algorithms is to introduce randomness to solve problems, even in the deterministic case. Such algorithms are often used in physical and mathematical problems and are most useful when it is difficult or impossible to use other approaches due to limitations, such as cost of performing experiment or inability to take direct measures. The problems typically require solving a stochastic partial differential equations (SDPEs), where an uncertainty is incorporated in to the model. For example: as an input parameter, or initial boundary condition. Extensive efforts have been devoted to the development of accurate numerical algorithms, so that simulation predictions are reliable in since that the numerical errors are well understood and under control for practical problems. Multilevel Monte Carlo (MLMC) is a novel idea. Instead of sampling from the true solution, a sampling is done at different levels. Such approach is beneficial in terms of convergence rate. However for practical simulations a large number of problems has to be solved, with huge number of unknowns. Such computational restrictions naturally leads to challenging parallel algorithms. To overcome some of the limitations, here we consider a parallel implementation of MLMC algorithm for a model SPDE, that uses GPU acceleration for the permeability generation.
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 [
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. 2
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: 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: ˆ
x=0 Q(x, ω) :=
k(x, ω)∂np(x, ω)dx
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 fields is essential problem in solving the SPDE. Here
we consider a practical covariance proposed by [
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 [
Literature provides different numerical schemes for solving PDEs. In [
We shortly recall the idea proposed in [
Let {Ml : l = 0 . . . L} ∈ N be increasing sequence of numbers called levels, with corresponding quantities {QMl}lL=0, and s ≥ 2 be coarsening factor, such that we have Ml = sMl−1, for l = 0 . . . L. Defining Yl = QMl − QMl−1 and setting Y0 = QM0, we can use telescopic sum to write the following identitity for the expected value E
E[QM ] = E[QM0] +
E[QMl − QMl−1] =
E[Yl] (5) L l=1
L l=0 1 L 2 l=0 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 Nl samples. For the mean square error we have: nl = α (vl/tl) with Lagrangian multiplier α = (vl/tl) (8) L l=0 (6) (7) e(QMM,LN )2 = E[(QMM,LN − E[Q])2] =
Nl−1V [Yl] + (E[QM − E[Q])2 Our goal is to have: e(QMM,LN )2 =
Nl−1V [Yl] + (E[QM − E[Q])2 ≤ 2ε2
Denote vl = V [Yl], and let tl be the mean time computing difference Yl once, and T = lL=0 nltl be the total time for the computation. Minimizing T under the above constraint, with Lagrangian multipliers and turning it to integer value gives us:
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 [
(a) Permeability (b) Solution
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 ∗ 107 . 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. 0.1
Single Problem Performance
CPU Generation 1 GPU Generation 2 GPU Generation
The fine grid is of size 210 ∗ 210 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 [
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.