Jostein R. Natvig

Open-source MATLAB implementation of consistent discretisations on complex grids

Accurate geological modelling of features such as faults, fractures or erosion requires grids that are flexible with respect to geometry. Such grids generally contain polyhedral cells and complex grid-cell connectivities. The grid representation for polyhedral grids in turn affects the efficient implementation of numerical methods for subsurface flow simulations. It is well known that conventional two-point flux-approximation methods are only consistent for K-orthogonal grids and will, therefore, not converge in the general case. In recent years, there has been significant research into consistent and convergent methods, including mixed, multipoint and mimetic discretisation methods. Likewise, the so-called multiscale methods based upon hierarchically coarsened grids have received a lot of attention. The paper does not propose novel mathematical methods but instead presents an open-source Matlab® toolkit that can be used as an efficient test platform for (new) discretisation and solution methods in reservoir simulation. The aim of the toolkit is to support reproducible research and simplify the development, verification and validation and testing and comparison of new discretisation and solution methods on general unstructured grids, including in particular corner point and 2.5D PEBI grids. The toolkit consists of a set of data structures and routines for creating, manipulating and visualising petrophysical data, fluid models and (unstructured) grids, including support for industry standard input formats, as well as routines for computing single and multiphase (incompressible) flow. We review key features of the toolkit and discuss a generic mimetic formulation that includes many known discretisation methods, including both the standard two-point method as well as consistent and convergent multipoint and mimetic methods. Apart from the core routines and data structures, the toolkit contains add-on modules that implement more advanced solvers and functionality. Herein, we show examples of multiscale methods and adjoint methods for use in optimisation of rates and placement of wells.

Well-balanced finite volume schemes of arbitrary order of accuracy for shallow water flows

Many geophysical flows are merely perturbations of some fundamental equilibrium state. If a numerical scheme shall capture such flows efficiently, it should be able to preserve the unperturbed equilibrium state at the discrete level. Here, we present a class of schemes of any desired order of accuracy which preserve the lake at rest perfectly. These schemes should have an impact for studying important classes of lake and ocean flows.

Solving the euler equations on graphics processing units

The paper describes how one can use commodity graphics cards (GPUs) as a high-performance parallel computer to simulate the dynamics of ideal gases in two and three spatial dimensions. The dynamics is described by the Euler equations, and numerical approximations are computed using state-of-the-art high-resolution finite-volume schemes. These schemes are based upon an explicit time discretisation and are therefore ideal candidates for parallel implementation.

Visual simulation of shallow-water waves

A commodity-type graphics card (GPU) is used to simulate nonlinear water waves described by a system of balance laws called the shallow-water system. To solve this hyperbolic system we use explicit high-resolution central-upwind schemes, which are particularly well suited for exploiting the parallel processing power of the GPU. In fact, simulations on the GPU are found to run 15‐30 times faster than on a CPU. The simulated cases involve dry-bed zones and non-trivial bottom topographies, which are real challenges to the robustness and accuracy of the discretization.

Fast computation of multiphase flow in porous media by implicit discontinuous Galerkin schemes with optimal ordering of elements

We present a family of implicit discontinuous Galerkin schemes for purely advective multiphase flow in porous media in the absence of gravity and capillary forces. To advance the solution one time step, one must solve a discrete system of nonlinear equations. By reordering the grid cells, the nonlinear system can be shown to have a lower triangular block structure, where each block corresponds to the degrees-of-freedom in a single or a small number of cells. To reorder the system, we view the grid cells and the fluxes over cell interfaces as vertices and edges in a directed graph and use a standard topological sorting algorithm. Then the global system can be computed by processing the blocks sequentially using a standard Newton-Raphson algorithm for the degrees-of-freedom in each block. Decoupling the system offers greater control over the nonlinear solution procedure and reduces the computational costs, memory requirements, and complexity of the scheme significantly. In particular, the first-order version of the method may be at least as efficient as modern streamline methods when accuracy requirements or the dynamics of the flow allow for large implicit time steps.

Simulation and visualization of the Saint-Venant system using GPUs

We consider three high-resolution schemes for computing shallow-water waves as described by the Saint-Venant system and discuss how to develop highly efficient implementations using graphical processing units (GPUs). The schemes are well-balanced for lake-at-rest problems, handle dry states, and support linear friction models. The first two schemes handle dry states by switching variables in the reconstruction step, so that bilinear reconstructions are computed using physical variables for small water depths and conserved variables elsewhere. In the third scheme, reconstructed slopes are modified in cells containing dry zones to ensure non-negative values at integration points. We discuss how single and double-precision arithmetics affect accuracy and efficiency, scalability and resource utilization for our implementations, and demonstrate that all three schemes map very well to current GPU hardware. We have also implemented direct and close-to-photo-realistic visualization of simulation results on the GPU, giving visual simulations with interactive speeds for reasonably-sized grids.

Flow-based coarsening for multiscale simulation of transport in porous media

Geological models are becoming increasingly large and detailed to account for heterogeneous structures on different spatial scales. To obtain simulation models that are computationally tractable, it is common to remove spatial detail from the geological description by upscaling. Pressure and transport equations are different in nature and generally require different strategies for optimal upgridding. To optimize the accuracy of a transport calculation, the coarsened grid should generally be constructed based on a posteriori error estimates and adapt to the flow patterns predicted by the pressure equation. However, sharp and rigorous estimates are generally hard to obtain, and herein we therefore consider various ad hoc methods for generating flow-adapted grids. Common for all is that they start by solving a single-phase flow problem once and then continue to form a coarsened grid by amalgamating cells from an underlying fine-scale grid. We present several variations of the original method. First, we discuss how to include a priori information in the coarsening process, e.g. to adapt to special geological features or to obtain less irregular grids in regions where flow-adaption is not crucial. Second, we discuss the use of bi-directional versus net fluxes over the coarse blocks and show how the latter gives systems that better represent the causality in the flow equations, which can be exploited to develop very efficient nonlinear solvers. Finally, we demonstrate how to improve simulation accuracy by dynamically adding local resolution near strong saturation fronts.

Numerical Solution of the Polymer System by Front Tracking

The paper describes the application of front tracking to the polymer system, an example of a nonstrictly hyperbolic system. Front tracking computes piecewise constant approximations based on approximate Riemann solutions and exact tracking of waves. It is well known that the front tracking method may introduce a blowup of the initial total variation for initial data along the curve where the two eigenvalues of the hyperbolic system are identical. It is demonstrated by numerical examples that the method converges to the correct solution after a finite time, and that this time decreases with the discretization parameter.For multidimensional problems, front tracking is combined with dimensional splitting, and numerical experiments indicate that large splitting steps can be used without loss of accuracy. Typical CFL numbers are in the range 10–20, and comparisons with Riemann free, high-resolution methods confirm the high efficiency of front tracking.The polymer system, coupled with an elliptic pressure equation, models two-phase, three-component polymer flooding in an oil reservoir. Two examples are presented, where this model is solved by a sequential time stepping procedure. Because of the approximate Riemann solver, the method is non-conservative and CFL numbers must be chosen only moderately larger than unity to avoid substantial material balance errors generated in near-well regions after water breakthrough. Moreover, it is demonstrated that dimensional splitting may introduce severe grid orientation effects for unstable displacements that are accentuated for decreasing discretization parameters.

A Multiscale Mixed Finite Element Solver for Three Phase Black Oil Flow

Previous research has shown that multiscale methods are robust and capable of providing more accurate solutions than traditional upscaling methods. Multiscale methods solve the pressure equation on a coarse grid, but capture the effects from fine-scale heterogeneities through basis functions computed numerically from local single-phase problems on the underlying geocellular grid. Published results have so far been limited to simple Cartesian grids and/or incompressible flow. Here, we present a multiscale mixed finite-element method for three-phase black-oil flow on geomodels with industry-standard complexity. In particular, we discuss which effects can be incorporated in the multiscale basis functions and which effects should be modeled only on the coarsened simulation grid. Moreover, we describe how to handle degenerate hexahedral cells and non-matching interfaces that occur across faults. Finally, we present results of flow simulations on models of industry-standard complexity and demonstrate how multiscale methods can be used to simulate three-phase black-oil flow directly on high-resolution geomodels. The multiscale methods presented herein enable varying resolution and provide a systematic procedure for coarsening or refining the simulation model. Introduction For the oil industry to succeed in increasing oil recovery there is a growing trend for model-based decisions. New and exciting developments are seen in a variety of areas such as real-time reservoir management, uncertainty quantification, integrated operations, closed-loop management, and production optimization. Common to all these fields of endeavor is the requirement for fast flow simulation in which the simulation model is tightly coupled to the geology and dynamic data sources. However, there is a significant, and increasing, gap between the level of detail seen in geological models and the capabilities of contemporary reservoir simulators. Mature fields have a large amount of geological and geophysical data that can be used to create static models, and sizes of high-resolution geological models range from a few million and up to a billion cells. Contemporary reservoir simulators typically operate on model sizes from tens of thousands to a few million cells. Similarly, mature fields usually have a lot of dynamic data (pressure tests, production data, 4-D seismics, etc) that could be used to calibrate and history match the high-resolution geological models. Unfortunately, instead of focusing on understanding the physical characteristics of reservoirs and the economic consequences of different developments, a lot of valuable human resources is diverted to upscaling (and downscaling) and its negative consequences for the representation of heterogeneities and fluid flow. Upscaling is a costly process which additionally wastes much of the information inherent in high-resolution geological models since local flow structures are only preserved in an average sense on the upscaled grid. Enabling the oil industry to make a step-change in its work processes therefore calls for a radical speedup of flow simulation and for simulators that are equipped to utilize both static data and the vast amount of dynamic data that becomes available. As an example, it would be highly attractive if reservoir simulation could be performed at seismic resolution in order to use 4-D seismics to history-match simulation models. There are several technological developments that can contribute to a radical speedup of flow simulation: advances in hardware, parallel algorithms, improved (non)linear solvers, and alternative formulations (streamlines, operator splitting), to name a few. Another important contribution may come from multiscale methods, as will be discussed herein. Generally speaking, multiscale methods are numerical methods and strategies that aim to describe physical phenomena on coarse grids while accounting for the influence of fine-scale structures in the porous media. However, unlike traditional upscaling techniques, multiscale methods often provide a mechanism to recover an approximate fine-scale solutions. Multiscale modeling of flow and transport in porous media has become a hot research topic in recent years. A quite comprehensive overview of current developments is found in a recent issue of the Computational Geosciences journal (Juanes and Tchelepi 2008). Common for all these methods is that they seek efficient solutions of elliptic (or parabolic) equations with rough coefficients in the absence of scale separation, which is often assumed in many other multiscale methods. In the race for making a

Constrained pressure residual multiscale (CPR-MS) method for fully implicit simulation of multiphase flow in porous media

We develop the first multiscale method for fully implicit (FIM) simulations of multiphase flow in porous media, namely CPR-MS method. Built on the FIM Jacobian matrix, the pressure system is obtained by employing a Constrained Pressure Residual (CPR) operator. Multiscale Finite Element (MSFE) and Finite Volume (MSFV) methods are then formulated algebraically to obtain efficient and accurate solutions of this pressure equation. The multiscale prediction stage (first-stage) is coupled with a corrector stage (second-stage) employed on the full system residual. The converged solution is enhanced through outer GMRES iterations preconditioned by these first and second stage operators. While the second-stage FIM stage is solved using a classical iterative solver, the multiscale stage is investigated in full detail. Several choices for fine-scale pre- and post-smoothing along with different choices of coarse-scale solvers are considered for a range of heterogeneous three-dimensional cases with capillarity and three-phase systems. The CPR-MS method is the first of its kind, and extends the applicability of the so-far developed multiscale methods (both MSFE and MSFV) to displacements with strong coupling terms.