**2. Methodology**

Currently, geostatistical simulations are processed with a large number of MPS algorithms using various techniques. Besides the difference in process definition and algorithmic imple‐ mentation, these algorithms share similarities of the fundamental elements [36].

#### **2.1. General framework**

point geostatistical (MPS) simulation was proposed to produce geologically realistic structure by using high-order spatial statistics based on conceptual training images [5, 6]. The concept of training image is introduced from explicit to represent geological structures with a numeral image generated from outcrops, expert sketching, or conditional simulations of variogrambased methods. Since the training images can incorporate additional information such as the expert guesses, prior database, and physical models, besides the spatial features, the simula‐

Due to the ability of reconstructing geological realistic, MPS methods are gaining popularity, and various algorithms have been proposed including pixel-based algorithms [5, 8, 9], patternbased algorithms [10–13], and optimal algorithms [14–16]. These algorithms have been applied to broad fields such as oil and gas industry [17–19], fluid prediction [20, 21], climate modeling [22]. However, some application suffers from the computational burden routinely. Since MPS methods need to scan the training image, abstract patterns, and reconstruct the patterns in the simulation grid, physical memory and running time are challenging or even unusable for large-

Many effects have been made to decrease the central processing units (CUP) and RAM expense. Approaches such as multiple grid technique [23], optimization of data storage with a list structure [24], hierarchical decomposing [25], Fourier space transform [26] have been intro‐ duced, while the computational burden remains heavy for very large grids and complex spatial

With the development of hardware, utilization of multiple-core central processing units (CPU), or graphic processing units (GPU), for parallel applications are increasing in popularity in various fields including geostatistical simulation. In 2010, Mariethoz investigated the possi‐ bility to parallelize the MPS simulation process on realization level, path-level, or node-level and proposed a general conflict management strategy [27]. This strategy has been implement‐ ed on a patch-based SIMPAT method [28]. Parallel implements for other geostatistical algorithms, such as the parallel two-point geostatistical simulation [29], parallel pixel-based

algorithms [30], and parallel optimal algorithms [31] have been proposed constantly.

In this article, we will present the parallel several schemes of MPS simulation on many-core and GPU architectures. The Compute Unified Device Architecture (CUDA) that provides access between CUPs and GPUs is used to illustrate the parallel strategies [32, 33]. Examples of the two general MPS algorithms known as SENSIM and DS are implemented and compared [34, 35] with the original algorithms to present the ability of pattern reproduction and the

Currently, geostatistical simulations are processed with a large number of MPS algorithms using various techniques. Besides the difference in process definition and algorithmic imple‐

mentation, these algorithms share similarities of the fundamental elements [36].

tion using TIs is straightforward and smart [7].

142 Modeling and Simulation in Engineering Sciences

scale or pattern-rich simulation models.

improvement of computational performance.

**2. Methodology**

models.

Generally, the overall framework of these algorithms is constructed as follows:


For each algorithm, one or some of the basic steps may have some identity. In this section, we will focus on the two algorithms following different theory, single normal equation simulation (SNESIM) and DS, to illustrate the application of parallel strategies on MPS simulations.

#### **2.2. SNESIM review**

The methodology of single normal equation simulation (SNESIM) is a sequential paradigm developed from the random concepts theory with the property of generating conditional distribution with local conditioning instead of global conditioning. The local conditional distribution is aggregated with the local data event. Considering a category property S with K possible states {s(k), k = 1, …, K}, a data template τ<sup>n</sup> is comprised of a geometry of *n* vectors {hα, α = 1, …, n}. The data event dn centered at u is defined with the template and the corre‐ sponding *n* values {s(u + hα), α = 1, …, n}. Consequently, the conditional probability of the occurrence of state sk denoted as Prob{S(u) = s*k*|S(u*α*), *α* = 1 … n} or Prob{S(u) = sk|dn} is defined following Bayes' theorem:

$$\text{Prob}\left\{\mathbf{S}(\mathbf{u}) \equiv \mathbf{s}\_{\mathbf{k}} \middle| \mathbf{d}\_{\mathbf{n}}\right\} = \frac{\text{Prob}\left\{\mathbf{S}(\mathbf{u}) \equiv \mathbf{s}\_{\mathbf{k}} \text{ and } \mathbf{d}\_{\mathbf{n}}\right\}}{\text{Prob}\left\{\mathbf{d}\_{\mathbf{n}}\right\}}\tag{1}$$

Where Prob{dn} and Prob{S(u) = sk and dn} are the probability of the occurrence of dn and the associated to a central value S(u) = s*<sup>k</sup>* respectively. In practice, this probability can be obtained by counting the number of replicates of the training image, which is calculated by:

$$\text{Prob}\left\{\mathbf{S}(\mathbf{u}) = \mathbf{s}\_k | \mathbf{d}\_n\right\} = \frac{\mathbf{c}\_k(\mathbf{d}\_n)}{\mathbf{c}(\mathbf{d}\_n)}\tag{2}$$

c{dn} is the number of replicates of the conditioning data dn, and c*k*{dn} is the number of replicates inferred from c*k*(dn) of which the central node S{u} has a value of s(k).

To reduce the computational burden of scanning the training image, SNESIM analyzes the training image and constructs a search tree to store all the patterns before sequential simula‐ tion. Then the same data template is used to aggregate the local distribution of the simulation grid for sampling values from the data base to the simulation grid.

The implementation of SNESIM algorithm is shown in **Figure 1**.

**Figure 1.** Flowchart for SNESIM process.

As described, the memory consumption is largely related to the size of data template and the pattern richness of the training image. At the same time, the pattern retrieving time also increases for the complex cases. To address this problem, a GPU-based parallel method is proposed in the following section.

#### **2.3. Direct Sampling review**

Direct Sampling is a novel MPS technique that borrows ideas from a sampling method introduced by Shannon totally abandoning the conditional probability approach. Instead of the conditional distribution, a measurement distance is used to calculate the similarity between the local conditioning and the data event got from the training image along a random path. The distance method enables the application to both categorical and continuous variables. During the sequential simulation, for each node to be simulated in the simulation grid, the training image is scanned along a circular unilateral path [37] and the distance is calculated. As long as the distance is smaller than the defined threshold, the current data event in the training image is sampled, and the central value is assigned to the node to be simulated in the simulation grid directly. The flowchart of Direct Sampling is shown in **Figure 2**.

**Figure 2.** Flowchart for Direct Sampling process.

To reduce the computational burden of scanning the training image, SNESIM analyzes the training image and constructs a search tree to store all the patterns before sequential simula‐ tion. Then the same data template is used to aggregate the local distribution of the simulation

As described, the memory consumption is largely related to the size of data template and the pattern richness of the training image. At the same time, the pattern retrieving time also increases for the complex cases. To address this problem, a GPU-based parallel method is

Direct Sampling is a novel MPS technique that borrows ideas from a sampling method introduced by Shannon totally abandoning the conditional probability approach. Instead of the conditional distribution, a measurement distance is used to calculate the similarity between

grid for sampling values from the data base to the simulation grid.

The implementation of SNESIM algorithm is shown in **Figure 1**.

144 Modeling and Simulation in Engineering Sciences

**Figure 1.** Flowchart for SNESIM process.

proposed in the following section.

**2.3. Direct Sampling review**

Since this method does not need to store patterns, it releases the memory intensity; on the other hand, parameter is the key factor for Direct Sampling.

There are three main input parameters:


Sensitivity analysis of parameters [38] shows trade-offs between the quality of realizations and the CPU times.
