Accelerating Reservoir Simulation on Multi-core and Many-Core Architectures with Graph Coloring ILU(k)

Modern petroleum reservoir simulation serves as a primary tool for quantitatively managing reservoir production and planning new fields. It involves repeatedly solving the Jacobian of a set of strong nonlinear partial differential equations governing the mass and energy conduction and conservation. Most of the existing reservoir simulators adopt iterative solver with multiple stages of preconditioners, in which the incomplete LU (ILU) factorization is an outstanding universal smoother. However, it turns out that when the degree of freedom of each grid grows, ILU usually becomes the bottleneck of the solver. Moreover, ILU is difficult to parallelize due to its inherent data dependency. In this paper, we developed a sparse iterative solver with parallelized ILU and triangular solve using block-wise data structure. Compared with the state of art iterative solver on 14 industrial reservoir simulation matrices, the proposed ILU is 5.2x faster (on average) than the state of art iterative solver because of the block-wise data structure, which leads to 2.2x speedup on the total solver runtime. In addition, parallel ILU and triangular solve are developed to further accelerate the solver. To tackle the strong data dependency in ILU and triangular solve, we first partition the algorithm into separated tasks and construct a data flow graph to represent the data dependency. Then, tasks are scheduled in parallel according to the topological order of the data flow graph. On an 8-thread multicore architecture, we achieved another 3.6x speedup on ILU factorization, and 3.3x on triangular solve with good scalability.

A new approach to the application of direct method in petroleum reservoir simulation is introduced, This new approach is based on a modified form of the alternating diagonal ordering scheme that will rearrange each type of elements along different secondary diagonals in the coefficient matrix. In a typical 2-dimensional 3-phase petroleum reservoir simulation problem of 345 grid system, the new method is substantially faster than the band algorithm applied to the matrix resulting from the standard ordering. Introduction In petroleum reservoir simulation, a linear system of equations is obtained. The coefficient matrix of such a linear system takes different forms depending on the ordering of the grid system. The resulting system of equations can be solved either by direct methods or by iterative methods. Direct solutions are more reliable, but they normally require larger storage space and more computer time. Therefore, a development of a better direct method approach is needed. The ordering technique can be improved by the use of processes which are geared specifically to the direct solution. The method is the generate-and-solve algorithm used by Woo, et al. The drawback of this technique is the fact it requires large computer storage and overhead work to store and generate the nonlooping code. The yet known systematic optimal ordering scheme in petroleum reservoir simulation is the alternating diagonal scheme developed by Price and Coats. However, there is a major drawback to this ordering scheme. In the alternating diagonal scheme, the elements in the upper right and lower left corners of the coefficient matrix do not form perfect secondary diagonals. Thus, the solution of the resulting matrix requires band algorithm which treats the zeros within the band as nonzeros. The treatment of the zeros as nonzeros makes the band matrix technique require larger storage space and more computational labor. In this paper a modification to the alternating diagonal ordering is proposed. This modification will rearrange each type of elements along different secondary diagonals in the coefficient matrix. By this rearrangement no new element will be introduced in the process of factorization except in the lower right process of factorization except in the lower right corner of the composite matrix. Also the location of the original and new elements are known beforehand. This enables us to develop a general procedure to solve the system of equations in reservoir simulation. This general procedure eliminates the overhead calculations and extra storage needed for the generate-and-solve algorithm. Also this general procedure would eliminate the extra calculations required procedure would eliminate the extra calculations required on the zero elements in the band matrix algorithm. A NEW DIRECT SOLUTION APPROACH The new direct solution approach will be outlined to solve the following system of equations (1) The new approach begins with the definition of a new ordering scheme. Then the LU factorization will be applied to this new ordering. Finally the implementation of the new ordering scheme in typical reservoir simulators will be outlined. RESTRICTED ALTERNATING DIAGONAL ORDERING Here a modification to the alternating diagonal ordering that produces perfect secondary diagonals in the original A matrix is proposed. This ordering will be called restricted alternating diagonal (RAD). Figure 1 depicts the RAD ordering. Figure 2 shows the resulting coefficient matrix corresponding to RAD ordering scheme. The numbering of this scheme must go along the shortest diagonal to reduce the band width of the new nonzero elements that will be generated later in the process of factorization. process of factorization.

The parallelizable block ILU (incomplete LU) factorization preconditioners for a block-tridiagonal matrix have been recently proposed by the author. In this paper, we describe a parallelization of Krylov subspace methods with the block ILU factorization preconditioners on distributed-memory computers such as the Cray T3E, and then parallel performance results of a preconditioned Krylov subspace method are provided to evaluate the effectiveness and efficiency of the block ILU preconditioners on the Cray T3E.

This paper presents non-intrusive, efficient stochastic approaches for predicting uncertainties associated with petroleum reservoir simulations. The Monte Carlo simulation method, which is the most common and straightforward approach for uncertainty quantification in the industry, requires to perform a large number of reservoir simulations and is thus computationally expensive especially for large-scale problems. We propose an efficient and accurate alternative through the collocation-based stochastic approaches. The reservoirs are considered to exhibit randomly heterogeneous flow properties. The underlying random permeability field can be represented by the Karhunen-Loeve expansion (or principal component analysis), which reduces the dimensionality of random space. Two different collocation-based methods are introduced to propagate uncertainty of the reservoir response. The first one is the probabilistic collocation method that deals with the random reservoir responses by employing the orthogonal polynomial functions as the bases of the random space and utilizing the collocation technique in the random space. The second one is the sparse grid collocation method that is based on the multi-dimensional interpolation and high-dimensional quadrature techniques. They are non-intrusive in that the resulting equations have the exactly the same form as the original equations and can thus be solved with existing reservoir simulators. These methods are efficient since only a small number of simulations are required and the statistical moments and probability density functions of the quantities of interest in the oil reservoirs can be accurately estimated. The proposed approaches are demonstrated with a 3D reservoir model originating from the 9th SPE comparative project. The accuracy, efficiency, and compatibility are compared against Monte Carlo simulations. This study reveals that, compared to traditional Monte Carlo simulations, the collocation-based stochastic approaches can accurately quantify uncertainty in petroleum reservoirs and greatly reduce the computational cost.

Flux Corrected Transport (FCT) is a technique developed to reduce numerical diffusion around shock fronts which has not previously been used in reservoir engineering. The first section of this paper describes the application of FCT techniques to solving a one dimensional waterflooding problem, and shows how to include FCT in an IMPES simulator. A second feature of our program is the implicit timesplit treatment of capillary pressure effects. A comparison between laboratory data and simulator results is included, and good agreement is obtained.

Summary We describe a collection of block preconditioners for use in solving large, sparse, linear systems of equations by iterative methods, and we compare their performance with several point preconditioners in solving some systems arising in numerical reservoir simulation. We consider block preconditioners that handle either lines or planes in an implicit manner and pointwise incomplete LU (ILU) factorizations combined with partial elimination preprocessors. We conclude that the best of the pointwise methods are both more robust and faster, but that the best of the block methods are competitive for certain orderings of unknowns and require less storage.

In this paper three different schemes for the finite-difference approximation of ( ) terms with irregular spacings of the grid points are investigated. Several methods of selecting the grid spacing are considered, and numerical results are presented for some carefully selected examples. presented for some carefully selected examples. Results indicate that one of the most commonly used finite-difference schemes may unnecessarily result in large discretization errors; the other two schemes discussed result in greater accuracy and are no more difficult to use. Introduction Finite-difference methods are used widely for solving partial-differential equations encountered in petroleum reservoir simulation. In many practical problems it is necessary to refine the grid in certain problems it is necessary to refine the grid in certain parts of the reservoir, thus resulting in an irregular parts of the reservoir, thus resulting in an irregular grid spacing for the system. For example, the local refinement of grid is usually necessary in the study of the flow behavior near a well, i.e., the case of coning problems. On the other hand, it is usually advantageous to make the grid locally coarse over aquifers and large gas caps. In spite of the importance of irregular grids in reservoir simulation and other engineering problems, little is known about the effect of various schemes on the over-all discretization error. Most literature on the subject of irregular grid spacings is concerned with the treatment of irregular boundaries. Here, the properties of several difference approximations to properties of several difference approximations to the second-order operator Au = [ ] are investigated. Operator A is encountered in all diffusion-type equations, such as multiphase flow in porous media and heat conduction with temperature-dependent conductivity. For simplicity, we limit our study to one space variable in this study. All results reported here can be extended readily to any number of space variables.

Chapter 9 - Conjugate Gradient-Type Methods

Consider a preconditioner M based on an incomplete LU (or Cholesky) factorization of a matrix A. M−1, which represents an approximation of A−1, is applied by performing forward and back substitution steps; this can present a computational bottleneck. An alternative strategy is to directly approximate A−1 by explicitly computing M−1. Preconditioners of this kind are called sparse approximate inverse preconditioners. They constitute an important class of algebraic preconditioners that are complementary to the approaches discussed in the previous chapter. They can be attractive because when used with an iterative solver, they can require fewer iterations than standard incomplete factorization preconditioners that contain a similar number of entries while offering significantly greater potential for parallel computations.

We consider characteristics-based finite element methods for solving nonlinear, degenerate, advection-diffusion equations. These equations have applications in the simulation of petroleum reservoirs and groundwater aquifers and in the modeling of free boundary problems. Standard finite element Galerkin methods have been studied for these equations. In this paper, we analyze the characteristics-based finite element methods for them. The main difficulty in the analysis is that the equations are degenerate and the solution lacks regularity. Here we develop a technique that respects the degeneracy and the known minimal regularity. This technique is based on the Green operator for standard elliptic equations and is developed directly for the degenerate advection-diffusion equations. We concentrate our analysis on the modified method of characteristics (MMOC) and one of its variants, the modified method of characteristics with adjusted advection (MMOCAA), which conserves mass. We derive error estimates in various norms. The extension to other variants is discussed. The present technique is also applied to nondegenerate problems; error estimates previously obtained for the MMOC are derived under much weaker regularity assumptions on the solution, and the error estimates for the MMOCAA appear new even in the nondegenerate case. Finally, numerical results are presented to show the sharpness of the error estimates derived.

The use of graphics processing units (GPUs) in high-performance parallel computing continues to become more prevalent, often as part of a heterogeneous system. For years, CUDA has been the de facto programming environment for nearly all general-purpose GPU (GPGPU) applications. In spite of this, the framework is available only on NVIDIA GPUs, traditionally requiring reimplementation in other frameworks in order to utilize additional multi- or many-core devices. On the other hand, OpenCL provides an open and vendor-neutral programming environment and runtime system. With implementations available for CPUs, GPUs, and other types of accelerators, OpenCL therefore holds the promise of a write once, run anywhere ecosystem for heterogeneous computing. Given the many similarities between CUDA and OpenCL, manually porting a CUDA application to OpenCL is typically straightforward, albeit tedious and error-prone. In response to this issue, we created CU2CL, an automated CUDA-to-OpenCL source-to-source translator that possesses a novel design and clever reuse of the Clang compiler framework. Currently, the CU2CL translator covers the primary constructs found in CUDA runtime API, and we have successfully translated many applications from the CUDA SDK and Rodinia benchmark suite. The performance of the automatically translated applications via CU2CL is on par with their manually ported counterparts.

We present a study of the effectiveness of asynchronous incomplete LU factorization preconditioners for the time-implicit solution of compressible flow problems while exploiting thread-parallelism within a compute node. A block variant of the asynchronous fine-grained parallel preconditioner adapted to a finite volume discretization of the compressible Navier--Stokes equations on unstructured grids is presented, and convergence theory is extended to the new variant. Experimental (numerical) results on the performance of these preconditioners on inviscid and viscous laminar two-dimensional steady-state test cases are reported. It is found, for these compressible flow problems, that the block variant performs much better in terms of convergence, parallel scalability, and reliability than the original scalar asynchronous ILU preconditioner. For viscous flow, it is found that the ordering of unknowns may determine the success or failure of asynchronous block-ILU preconditioning, and an ordering of grid cells suitable for solving viscous problems is presented.

SummaryThis study explores the performance and scaling of a GMRES Krylov method employed as a smoother for an algebraic multigrid preconditioned Newton‐Krylov solution approach applied to a fully implicit variational multiscale finite element resistive magnetohydrodynamics formulation. In this context, a Newton iteration is used for the nonlinear system and a parallel MPI‐based Krylov method is employed for the linear subsystems. The efficiency of this approach is critically dependent on the scalability and performance of the parallel algebraic multigrid preconditioner for the linear solutions and the performance of the multigrid smoothers play a critical role. Krylov multigrid smoothers are considered in an attempt to reduce the time and memory requirements of existing robust smoothers based on additive Schwarz domain decomposition with incomplete LU factorization solves on each subdomain. Three time‐dependent resistive magnetohydrodynamics test cases are considered to evaluate the method. Compared with a domain decomposition incomplete LU smoother, the GMRES smoother can reduce the solve time due to a significant decrease in the preconditioner setup time and often a reduction in outer Krylov solver iterations, and requires less memory, typically 35% less memory.

Practical considerations for adaptive implicit methods in reservoir simulation

Nowadays, the efficient and probabilistic data-structure named Bloom filter is widely used in network processing applications. Moreover, the parallel nature of Bloom filter has made it suitable for multi/many-core architectures. In this paper, two schemes called multi-core architecture with shared queue and multi-core architecture with private queue (both employing Bloom filter cores) are proposed and implemented on FPGA. The inherent parallelism in querying of packets and different number of cores (such as 1, 2, 4, 8 and 16 cores) are considered. Experimental results show that the multi-core architecture with private queue achieves higher throughput than the latter one. Furthermore, Bloom filter is also implemented on GPU (as many-core architecture) and the results are compared to the CPU only version. When the number of packets in GPU memory is 16384, the speedup achieved by GPU implementations using CUDA is about 274 times compared to CPU implementation. However FPGA results outperform GPU, so that the throughput of first architecture (shared queue) and second architecture (private queue) with 16 cores are respectively almost 5.5 and 7.1 times higher than GPU throughput.