**1. Introduction**

Geostatistical stochastic simulation is important for the research of geological phenomenon. In the past several decades, a large number of geostatistical methods have been developed based on the spatial covariance properties of the geological data. The traditional tool to quantify the spatial covariance is known as variogram, which measures the covariance among any two points separated by a certain distance [1]. Although variogram-based methods are successfully applied to multi-Gaussian system, they have limitations for the characterization of complex systems such as the curvilinear or long-range continuous facies [2–4]. An alternative known as multiple-

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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‐ tion using TIs is straightforward and smart [7].

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 largescale or pattern-rich simulation models.

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 models.

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 improvement of computational performance.
