Mixed precision parallel preconditioner and linear solver for compositional reservoir simulation

Reservoir simulation role in value creation and strategic management decisions cannot be over emphasized. Simulation of complex challenging reservoirs with millions of grid blocks especially in compositional mode is very time-consuming even with fast modern computers. On the other hand, high price of cluster supercomputers prevents them for being commonly used for fast simulation of such reservoirs. In recent years, the development of many-core processors like cell processors, DSPs, and graphical processing units (GPUs) has provided a very cost-effective hardware platform for fast computational operations. However, programming for such processors is much more difficult than conventional CPUs, and new parallel algorithm design and special parallel implementation methods are needed. Using the computational power of CPUs, GPUs, and/or any other processing unit, Open Computing Language (OpenCL) provides a framework for programming for heterogeneous platforms. In this paper, OpenCL is used to employ the computational power of a GPU to build a preconditioner and solve the linear system arising from compositional formulation of multiphase flow in porous media. The proposed parallel preconditioner is proved to be quite effective, even in heterogeneous porous media. Using data-parallel modules on GPU, the preconditioner/solver runtime reduced at least 1 order of magnitude compared to their serial implementation on CPU.     

A GPU parallel computing strategy for the material point method

The material point method (MPM), which is a combination of the finite element and meshfree methods, suffers from significant computational workload due to the fine mesh that is required in spite of its advantages in simulating large deformations. This paper presents a parallel computing strategy for the MPM on the graphics processing unit (GPU) to boost the method’s computational efficiency. The interaction between a structural element and soil is investigated to validate the applicability of the parallelisation strategy. Two techniques are developed to parallelise the interpolation from soil particles to nodes to avoid a data race; the technique that is based on workload parallelisation across threads over the nodes has a higher computational efficiency. Benchmark problems of surface footing penetration and a submarine landslide are analysed to quantify the speedup of GPU parallel computing over sequential simulations on the central processing unit. The maximum speedup with the GPU used is ∼30 for single-precision calculations and decreases to ∼20 for double-precision calculations.     

A direct boundary element method for plane strain poroelasticity

This paper presents a direct boundary element method (BEM), formulated in the Laplace transform space, for plane strain poroelasticity. The paper expands on work by Cheng and Liggett by recasting the theoretical foundation of BEM within the framework of Rice and Cleary's formulation of the Biot theory of poroelasticity. Furthermore, the numerical algorithm is generalized to deal with both interior and exterior domain problems, and a method for indirectly calculating the Cauchy principal value of the singular integrals is presented. Formulae for the stress and flux inside the domain are also derived. Finally, the algorithm is validated by comparing the numerical results with the analytical solution of a borehole subject to a far-field deviatoric stress (exterior domain) and with the solution of Mandel's problem (interior domain). These two examples provide a critical test of the algorithm.     

A boundary element method for identifying orthotropic material parameters

A new approach is proposed for identifying elastic constants for orthotropic material bodies by using the boundary element method. The material parameters which must be determined are two Young's moduli, a shear modulus and two Poisson's ratios. The method is applied to geotechnical engineering problems.     

A linear constitutive model for unsaturated poroelasticity by micromechanical analysis

This paper extends the Biot theory of poroelasticity from the saturated to unsaturated case. The Biot phenomenological model uses parameters that are easily observable, such as the deformation of porous frame, total stress, pore pressure, and fluid specific discharge. Such model is preferred for engineering applications. At this macroscopic level, the extension of Biot theory from saturated to unsaturated is straightforward. The constitutive constants, however, are combined properties of solid, pore space, and fluids. In the unsaturated case, the constants are functions of the degree of saturation. Their measurements and tabulation over a range of saturation is generally not feasible for practical applications. In this work, a Biot-Willis type analysis is performed for the unsaturated case to provide a theory that the bulk material constants can be evaluated using laboratory measurable micromechanical constants under saturated condition, plus a capillary pressure curve (saturation versus suction pressure) typically available for unsaturated porous medium, without the need of measurement at each state of saturation. In particular, it is demonstrated that the surface energy contained in the meniscus interface manifests as a “capillary modulus,” given by the negative inverse slope of the capillary pressure curve. A rigorous analysis based on the thermodynamic variational energy approach is also conducted to lend theoretical support to the phenomenological approach. The presented model can bring a closure to the practical engineering modeling of the deformation of partially saturated porous medium that lacks the information of material constants over the range of saturation.     

一种高性能迭代预处理方法及其在岩土工程中的应用

岩土工程建设的发展极大地促进了三维数值模拟的应用。大规模三维有限元计算需要求解一系列大型线性方程组,这些线性方程组的求解直接影响着整个有限元计算的效率。复杂岩土工程问题通常涉及多相和多体耦合相互作用,各相之间或不同固体材料之间性质差别显著,可能导致Krylov子空间迭代法收敛缓慢,甚至求解失败。为了提高Krylov子空间迭代法的求解效率和可靠性,提出一种新的高效预处理技术,通过算例验证了所提出的分区块迭代预处理方法的有效性。     

A coupling of mixed and discontinuous Galerkin finite-element methods for poroelasticity

In this paper, we formulate a finite-element procedure for approximating the coupled fluid and mechanics in Biot’s consolidation model of poroelasticity. We approximate the flow variables by a mixed finite-element space and the displacement by a family of discontinuous Galerkin methods. Theoretical convergence error estimates are derived and, in particular, are shown to be independent of the constrained specific storage coefficient, c _o . This suggests that our proposed algorithm is a potentially effective way to combat locking, or the nonphysical pressure oscillations, which sometimes arise in numerical algorithms for poroelasticity.     

Online conservative generalized multiscale finite element method for highly heterogeneous flow models

Computational Geosciences - In this work, we consider an online enrichment procedure in the context of the Generalized Multiscale Finite Element Method (GMsFEM) for the two-phase flow model in...     

A supercoarsening multigrid method for poroelasticity in 3D coupled flow and geomechanics modeling

The Galerkin finite-element discretization of the force balance equation typically leads to large linear systems for geomechanical problems with realistic dimensions. In iteratively coupled flow and geomechanics modeling, a large linear system is solved at every timestep often multiple times during coupling iterations. The iterative solution of the linear system stemming from the poroelasticity equations constitutes the most time-consuming and memory-intensive component of coupled modeling. Block Jacobi, LSOR, and Incomplete LU factorization are popular preconditioning techniques used for accelerating the iterative solution of the poroelasticity linear systems. However, the need for more effective, efficient, and robust iterative solution techniques still remains especially for large coupled modeling problems requiring the solution of the poroelasticity system for a large number of timesteps. We developed a supercoarsening multigrid method (SCMG) which can be multiplicatively combined with commonly used preconditioning techniques. SCMG has been tested on a variety of coupled flow and geomechanics problems involving single-phase depletion and multiphase displacement of in-situ hydrocarbons, CO₂ injection, and extreme material property contrasts. Our analysis indicates that the SCMG consistently improves the convergence properties of the linear systems arising from the poroelasticity equations, and thus, accelerates the coupled simulations for all cases subject to investigation. The joint utilization of the two-level SCMG with the ILU1 preconditioner emerges as the most optimal preconditioning/iterative solution strategy in a great majority of the problems evaluated in this work. The BiCGSTAB iterative solver converges more rapidly compared to PCG in a number of test cases, in which various SCMG-accelerated preconditioning strategies are applied to both iterators.     

中国高放废物地质处置缓冲材料大型试验台架和热-水-力-化学耦合性能研究

缓冲材料作为高放废物深地质处置库中一道重要的人工屏障,与高放废物容器和处置库围岩直接接触,在高放废物衰变热、辐射作用和地下水等影响下产生复杂的热-水-力-化学耦合作用,为了验证缓冲材料是否能长期有效地发挥其屏障材料的作用,核工业北京地质研究院利用高庙子钠基膨润土组装并运行了模拟中国高放废物地质处置室 尺寸的大型缓冲材料膨润土试验台架（China-Mock-Up）。建立了缓冲材料试验台架的安装和试验方法,依据实测数据和理论分析,揭示了热-水-力-化学耦合作用条件下膨润土中的相对湿度是在加热器的热效应和外部供水的湿效应共同作用下发生变化的,压实膨润土中应力的变化主要是由于膨润土遇水膨胀和加热器的热效应引起的,试验验证了模拟高放废物地质处置室内加热器（废物罐）运行初期的位移过程,为缓冲材料和高放废物地质处置库的设计提供了重要的工程参数和理论依据。     

Stochastic computational modelling of highly heterogeneous poroelastic media with long‐range correlations

The compaction of highly heterogeneous poroelastic reservoirs with the geology characterized by long‐range correlations displaying fractal character is investigated within the framework of the stochastic computational modelling. The influence of reservoir heterogeneity upon the magnitude of the stresses induced in the porous matrix during fluid withdrawal and rock consolidation is analysed by performing ensemble averages over realizations of a log‐normally distributed stationary random hydraulic conductivity field. Considering the statistical distribution of this parameter characterized by a coefficient of variation governing the magnitude of heterogeneity and a correlation function which decays with a power‐law scaling behaviour we show that the combination of these two effects result in an increase in the magnitude of effective stresses of the rock during reservoir depletion. Further, within the framework of a perturbation analysis we show that the randomness in the hydraulic conductivity gives rise to non‐linear corrections in the upscaled poroelastic equations. These corrections are illustrated by a self‐consistent recursive hierarchy of solutions of the stochastic poroelastic equations parametrized by a scale parameter representing the fluctuating log‐conductivity standard deviation. A classical example of land subsidence caused by fluid extraction of a weak reservoir is numerically simulated by performing Monte Carlo simulations in conjunction with finite elements discretizations of the poroelastic equations associated with an ensemble of geologies. Numerical results illustrate the effects of the spatial variability and fractal character of the permeability distribution upon the evolution of the Mohr–Coulomb function of the rock. Copyright © 2004 John Wiley & Sons, Ltd.     

Computation of large-scale parameters for dispersion in fissured porous medium using finite-volume method

In this work, we are interested in solute transport in fissured porous media. The medium is considered as a special case of a two-region system, and a two-equation model previously obtained from a volume averaging technique is used to derive large-scale dispersion coefficients. These coefficients are obtained as solutions of a set of closure problems and the main objective of this work is to present an efficient method to solve these closure problems. The method makes use of an unstructured grid and special techniques to take into account the fissure network. Results are compared with other existing methods on simple fissured media. Finally, the technique is applied to a complex structure.     

Convergence analysis of fixed stress split iterative scheme for anisotropic poroelasticity with tensor Biot parameter

We perform a convergence analysis of the fixed stress split iterative scheme for the Biot system modeling coupled flow and deformation in anisotropic poroelastic media with tensor Biot parameter. The fixed stress split iterative scheme solves the flow subproblem with all components of the stress tensor frozen using a multipoint flux mixed finite element method, followed by the poromechanics subproblem using a conforming Galerkin method in every coupling iteration at each time step. The coupling iterations are repeated until convergence and Backward Euler is employed for time marching. The convergence analysis is based on studying the equations satisfied by the difference of iterates to show that the fixed stress split iterative scheme for anisotropic poroelasticity with Biot tensor is contractive. We also demonstrate that the scheme is numerically convergent using the classical Mandel’s problem solution for transverse isotropy.     

大型地下洞室考虑开挖卸荷效应的位移反分析

基于现场监测资料的位移反分析是地下工程动态监控、信息化施工的重要组成部分。以乌江彭水水电站大型地下厂房（开挖跨度为30 m,高度为78.5 m）为例,从围岩实测位移出发,建立了基于均匀设计－神经网络－遗传算法的围岩力学参数的系统反分析方法,反演考虑开挖卸荷效应的围岩力学参数。根据数值分析结果形成训练样本,利用BP人工神经网络映射围岩的变形与力学参数的关系,同时针对传统人工神经网络存在初始权值难以确定的问题,应用遗传算法优化神经网络的初始权值;利用现场监测的增量变形反演了围岩的力学参数;最后利用反演出的参数,进行地下厂房开挖预测分析。结果表明,预测位移与现场监测位移较为接近,进行统计检验结果为优,说明该参数反演方法是正确合理的。     

A coupling of mixed and continuous Galerkin finite element methods for poroelasticity II: the discrete-in-time case

In this paper, we formulate a finite element procedure for approximating the coupled fluid and mechanics in Biot’s consolidation model of poroelasticity. Here, we approximate the pressure by a mixed finite element method and the displacements by a Galerkin method. Theoretical convergence error estimates are derived in a discrete-in-time setting. Of particular interest is the case when the lowest-order Raviart–Thomas approximating space or cell-centered finite differences are used in the mixed formulation and continuous piecewise linear approximations are used for displacements. This approach appears to be the one most frequently applied to existing reservoir engineering simulators.     

A global method for coupling transport with chemistry in heterogeneous porous media

Modeling reactive transport in porous media, using a local chemical equilibrium assumption, leads to a system of advection–diffusion PDEs coupled with algebraic equations. When solving this coupled system, the algebraic equations have to be solved at each grid point for each chemical species and at each time step. This leads to a coupled non-linear system. In this paper, a global solution approach that enables to keep the software codes for transport and chemistry distinct is proposed. The method applies the Newton–Krylov framework to the formulation for reactive transport used in operator splitting. The method is formulated in terms of total mobile and total fixed concentrations and uses the chemical solver as a black box, as it only requires that one be able to solve chemical equilibrium problems (and compute derivatives) without having to know the solution method. An additional advantage of the Newton–Krylov method is that the Jacobian is only needed as an operator in a Jacobian matrix times vector product. The proposed method is tested on the MoMaS reactive transport benchmark.     

Micromechanical parameters in bonded particle method for modelling of brittle material failure

Bonded particle modelling (BPM) is nowadays being extensively used for simulating brittle material failure. In BPM, material is modelled as a dense assemblage of particles (grains) connected together by contacts (cement). This sort of modelling seriously depends on the mechanical properties of particle and contact, which are named here as micro‐parameters. However, a definite calibration methodology to obtain micro‐parameters has not been so far established; and many have reported some serious problems. In this research, a calibration procedure to find a unique set of micro‐parameters is established. To attain this purpose, discrete element code of UDEC is used to perform BPM. This code can be conveniently developed by the user. The proposed BPM is composed of rigid polygonal particles interacting at their contact points. These contacts can undergo a certain amount of tension, and their shear resistance is provided by cohesion and friction angle. The results demonstrate that each material macro‐property (i.e. Young's modulus, Poisson's ratio, internal friction angel, internal cohesion, and tensile strength) is directly originated from and distinctly related to the contact properties (i.e. normal and shear stiffness, friction angel, cohesion, and tensile strength). Copyright © 2010 John Wiley & Sons, Ltd.     

A parsimonious approach for large-scale tracer test interpretation

Hydrogeology Journal - Dye tracing is an efficient method for spring watershed delineation, but is also used in surface waters to assess pollution migration over several kilometers. The aim of this...