N many fields of science and engineering, such as fluid or structural mechanics, and nanotechnology, dynamical systems at different scale need to be simulated, optimized or controlled. They are often described by discretizations of systems of nonlinear partial differential equations yielding high-dimensional discrete phase spaces. For this reason, in recent decades, research has mainly focused on the development of sophisticated analytical and numerical (linear and nonlinear) tools to help understand the overall multiscale system behavior. Various models and numerical methods have been developed to simulate different physical processes at different scales. The choice of these methods will depend largely on the problem, the available computational resources and constitutive equations. Smoothed particle hydrodynamics (SPH) was developed a few decades ago to model inviscid fluid and gas flow dynamics in astrophysical problems. The SPH is an interpolation-based numerical technique that can be used to solve systems of partial differential equations (PDEs) using either Lagrangian or Eulerian descriptions. The nature of SPH method allows to incorporate different physical and chemical effects into the discretized governing equations with relatively small code-development effort. In addition, geometrically complex and/or dynamic boundaries, and interfaces can be handled without undue difficulty. The SPH numerical procedure of calculating state variables (i.e., density, velocity, and gradient of deformation) are computed as a weighted average of values in a local region. @en

Smoothed particle hydrodynamics (SPH) has been extensively used to model high and low Reynolds number flows, free surface flows and collapse of dams, study pore-scale flow and dispersion, elasticity, and thermal problems. In different applications, it is required to have a stable and accurate discretization of the elliptic operator with homogeneous and heterogeneous coefficients. In this paper, the stability and approximation analysis of different SPH discretization schemes (traditional and new) of the diagonal elliptic operator for homogeneous and heterogeneous media are presented. The optimum and new discretization scheme of specific shape satisfying thetwo-point flux approximation nature is also proposed. This scheme enhances the Laplace approximation (Brookshaw’s scheme (1985) and Schwaiger’s scheme (2008)) used in the SPH community for thermal, viscous, and pressure projection problems with an isotropic elliptic operator. The numerical results are illustrated by numerical examples, where the comparison between different versions of the meshless discretization methods are presented.@en

Commercial reservoir simulators must be very robust and fast. Moreover, current hardware requires the simulators to scale over multiple number of computing nodes and for a fixed (‘strong scalability’) as well as an increasing problem size per computing node (‘weak scalability’). In most current commercial reservoir simulators, due to the different geological structures and properties of hydrocarbon reservoirs and the use of enhanced oil recovery (EOR) techniques, the governing equations are strongly nonlinear and hard to solve. The Jacobian system is solved by FGMRES preconditioned by the two-level constrained pressure residual (CPR) preconditioner. The driving force of the CPR preconditioner is the solution of the pressure equation. The industry standard for solving the pressure equation is the algebraic multigrid (AMG) solver. AMG is well known for its ‘weak scalability’. However, in these applications, AMG has unfavorable ‘strong’ scalability properties. This degradation in scalability is due to the increased level of inter-processor communication in the algorithm. In this paper, a monotone non-Galerkin AMG (MNG-AMG) method is presented. The aim of the method is to reduce the overall communication in MNG-AMG by enforcing a predefined nonzero pattern and monotonicity property (i.e., M-matrices) on each multigrid level. This paper describes the application of the MNG-AMG method in the context of reservoir simulations. We will compare the parallel scalability of the default solver with the MNG-AMG solver and discuss the optimal values for the MNG-AMG solver for a variety of test cases based on full field reservoir simulations.@en

Commercial reservoir simulators must be very robust and fast. Moreover, current hardware requires the simulators to scale over multiple number of computing nodes and for a fixed (‘strong scalability’) as well as an increasing problem size per computing node (‘weak scalability’). In most current commercial reservoir simulators, due to the different geological structures and properties of hydrocarbon reservoirs and the use of enhanced oil recovery (EOR) techniques, the governing equations are strongly nonlinear and hard to solve. The Jacobian system is solved by FGMRES preconditioned by the two-level constrained pressure residual (CPR) preconditioner. The driving force of the CPR preconditioner is the solution of the pressure equation. The industry standard for solving the pressure equation is the algebraic multigrid (AMG) solver. AMG is well known for its ‘weak scalability’. However, in these applications, AMG has unfavorable ‘strong’ scalability properties. This degradation in scalability is due to the increased level of inter-processor communication in the algorithm. In this paper, a monotone non-Galerkin AMG (MNG-AMG) method is presented. The aim of the method is to reduce the overall communication in MNG-AMG by enforcing a predefined nonzero pattern and monotonicity property (i.e., M-matrices) on each multigrid level. This paper describes the application of the MNG-AMG method in the context of reservoir simulations. We will compare the parallel scalability of the default solver with the MNG-AMG solver and discuss the optimal values for the MNG-AMG solver for a variety of test cases based on full field reservoir simulations.@en

Commercial reservoir simulators must be very robust and fast. Moreover, current hardware requires the simulators to scale over multiple number of computing nodes and for a fixed (‘strong scalability’) as well as an increasing problem size per computing node (‘weak scalability’). In most current commercial reservoir simulators, due to the different geological structures and properties of hydrocarbon reservoirs and the use of enhanced oil recovery (EOR) techniques, the governing equations are strongly nonlinear and hard to solve. The Jacobian system is solved by FGMRES preconditioned by the two-level constrained pressure residual (CPR) preconditioner. The driving force of the CPR preconditioner is the solution of the pressure equation. The industry standard for solving the pressure equation is the algebraic multigrid (AMG) solver. AMG is well known for its ‘weak scalability’. However, in these applications, AMG has unfavorable ‘strong’ scalability properties. This degradation in scalability is due to the increased level of inter-processor communication in the algorithm. In this paper, a monotone non-Galerkin AMG (MNG-AMG) method is presented. The aim of the method is to reduce the overall communication in MNG-AMG by enforcing a predefined nonzero pattern and monotonicity property (i.e., M-matrices) on each multigrid level. This paper describes the application of the MNG-AMG method in the context of reservoir simulations. We will compare the parallel scalability of the default solver with the MNG-AMG solver and discuss the optimal values for the MNG-AMG solver for a variety of test cases based on full field reservoir simulations.@en

Preconditioning can be used to damp slowly varying error modes in the linear solver residuals, corresponding to extreme eigenvalues. Existing multiscale solvers use a sequence of aggressive restriction, coarse-grid correction and prolongation operators to handle low-frequency modes on the coarse grid. High-frequency errors are then resolved by employing a smoother on fine grid. In reservoir simulations, the Jacobian system is usually solved by FGMRES method with two-level Constrained Pressure Residual (CPR) preconditioner. In this paper, a parallel fully implicit smoothed particle hydrodynamics (SPH) based multiscale method for solving pressure system is presented. The prolongation and restriction operators in this method are based on a SPH gradient approximation (instead of solving localized flow problems) commonly used in the meshless community for thermal, viscous, and pressure projection problems. This method has been prototyped in a commercially available simulator. This method does not require a coarse partition and can be applied to general unstructured topology of the fine scale. The SPH based multiscale method provides a reasonably good approximation to the pressure system and speeds up the convergence when used as a preconditioner for an iterative fine-scale solver. In addition, it exhibits expected good scalability during parallel simulations. Numerical results are presented and discussed.@en

The reservoir simulation of the complex reservoirs with anisotropic permeability,which includes faults and non-orthogonal grids, with a fully discontinuous permeability tensor in the discretization is a major challenge. Several methods have already been developed and implemented within industry standard reservoir simulators for non-orthogonal grids (e.g., Multi-Point Flux Approximation (MPFA) “O” method). However, it has been noticed that some of the numerical methods for elliptic/parabolic equations may violate the maximum principle (i.e., lead to spurious oscillations), especially when the anisotropy is particularly strong. It has been found that the oscillations are closely related to the poor approximation of the pressure gradient in the flux computation. Therefore, proposed methods must correctly approximate underlying operators, satisfy a discrete maximum principle and have coercivity properties. Furthermore, the method must be robust and efficient. This paper presents the meshless multi-point flux approximation of second order elliptic operators containing a tensor coefficient. The method is based on a pressure gradient approximation commonly used in meshless methods (or Smoothed Particle Hydrodynamics method—SPH method). The proposed discretization schemes can be written as a sum of sparse positive semidefinite matrix and perturbation matrix. We show that convergence rates are retained as for finite difference methods O(hα), 1 ≤ α < 2, where h denotes the maximum particle spacing. The results are presented, discussed and future studies are outlined.@en

The reservoir simulation of the complex reservoirs with anisotropic permeability,which includes faults and non-orthogonal grids, with a fully discontinuous permeability tensor in the discretization is a major challenge. Several methods have already been developed and implemented within industry standard reservoir simulators for non-orthogonal grids (e.g., Multi-Point Flux Approximation (MPFA) “O” method). However, it has been noticed that some of the numerical methods for elliptic/parabolic equations may violate the maximum principle (i.e., lead to spurious oscillations), especially when the anisotropy is particularly strong. It has been found that the oscillations are closely related to the poor approximation of the pressure gradient in the flux computation. Therefore, proposed methods must correctly approximate underlying operators, satisfy a discrete maximum principle and have coercivity properties. Furthermore, the method must be robust and efficient. This paper presents the meshless multi-point flux approximation of second order elliptic operators containing a tensor coefficient. The method is based on a pressure gradient approximation commonly used in meshless methods (or Smoothed Particle Hydrodynamics method—SPH method). The proposed discretization schemes can be written as a sum of sparse positive semidefinite matrix and perturbation matrix. We show that convergence rates are retained as for finite difference methods O(hα), 1 ≤ α < 2, where h denotes the maximum particle spacing. The results are presented, discussed and future studies are outlined.@en