© J.-M. Gratien et al., published by IFP Energies nouvelles, 2016

Introduction

Reservoir simulation models are nowadays a key element for oil and gas companies in the development of fields. Reservoir simulation is used extensively to design advanced EOR processes like miscible displacement by CO₂ or nitrogen, chemical flooding (polymer, alkaline, surfactant, etc.) or thermal processes with in situ combustion or water steam injection, see Lake [1]. Such models with very complex physics are usually time-consuming and expensive as they lead to solve strongly non linear numerical models. To preserve the predictivity and the accuracy of the models with reasonable running simulation time, it is therefore crucial to improve the runtime performance of the simulators.

A posteriori estimate aims at controlling the error during the simulation and enables to develop adaptive algorithms leading to a significant computational saving. First rigorous approachs can be found in [2-4] suitable for unsteady nonlinear models, and [5-7] for degenerated problems. Recently, abstract frameworks for a posteriori estimates appeared, for two-phase flow [8, 9], multiphase compositional model [10], and the thermal multiphase compositional global model [11]. In these frameworks the dual norm of the residual augmented by a nonconformity evaluation term is controlled by fully computable estimators decomposed into space, time, linearization, and algebraic error components. This allows to formulate criteria for stopping the iterative algebraic solver and the iterative linearization solver when the corresponding error components do not affect significantly the overall error. Moreover, the spatial and temporal error components can be balanced by adaptations of both time step and space mesh. The aforementioned approach is promising. Significant gains in terms of the total number of algebraic solver iterations and important reductions in terms of the number of cells can be obtained in comparison with standard options of reservoir simulators. However, no straightforward comparison of the computational cost has been measured on reservoir engineering practices.

In this paper we propose an engineering approach of the a posteriori error estimates, provided by recent works [10, 11], in order to establish an efficient computation of the estimators in such a way that we can ensure important computational savings for realistic complex models. Our approach is based on providing a low-cost evaluation formula by avoiding flux reconstructions. As mentioned before we rely here on the general ideas of [10, 11], we apply a posteriori error estimators allowing to distinguish the different components of the global error, and we use them to formulate stopping criteria for the iterative algebraic and nonlinear solvers together with refinement/coarsening criteria for both the time step and the space mesh. The present work brings two major novelties: first, we propose a well-implemented evaluation of the estimators using a practical simplified formula making them usable in engineering practices; second, we apply the stopping criteria for the iterative algebraic solver in the context of a three-dimensional realistic reservoir simulation study.

The paper is organized as follows. In Section 1 we describe a model problem, the two-phase flow model, and we discuss its numerical resolution. In Section 2 we state the corresponding a posteriori error estimate. Then, we introduce in Section 3 a practical simplified formula for the evaluation of the estimators. In Section 4 we discuss the use of a posteriori error estimation tools to enhance the runtime reservoir simulator performance. Finally, in Section 5 we validate our approach on a real-life reservoir engineering example in three dimensions.

1 Reservoir model and numerical resolution

We introduce in this section a two-phase flow model, its finite volume discretization, the linearization of the resulting discrete system and the algebraic resolution. In [10, 11] a posteriori error estimates are driven for the thermal multiphase compositional general model. However, in this work we focus on the two-phase flow since it is the model of our chosen study case (SPE10 case, see Sect. 5), and for the clarity of the simplified evaluation process.

1.1 Model Problem

We consider an immiscible two-phase model where the phases are water and oil. Each of them is composed of only one component, water component and hydrocarbon component, respectively. As we have an immiscible model we will use the same index to represent the phase and its component. Consequently, we use the indices w,o to represent the water and the oil, respectively.

We have the following relation between the phase pressures(1)

where is the capillary pressure. Consequently, the unknowns of the model are the water pressure and the saturations , with . These unknowns are collected in the vector . We denote by the porosity of the medium and by the permeability tensor. For a phase is the relative permeability, is the dynamic viscosity and is the density. The system of governing equations is given by(2)

with no-flow boundary conditions. For all represents the average phase velocity given by Darcy’s law,(3)

where denotes the gravity vector acting in the negative z direction and g its Euclidian norm.

1.2 Discretization

To discretize our model we choose a popular scheme in the oil industry, the implicit in time cell centered finite volume scheme, that is fast and simple to implement and well-known by its stability.

In order to define this numerical scheme we start by giving some notations for the space-time mesh.

1.3 Discrete Setting

Let denote a sequence of positive real numbers corresponding to the discrete time steps such that . We consider the discrete times such that and, for ; then we define the time intervals . Let denote a family of meshes of the space domain superadmissible in the sense of Eymard et al. [12, Definition 3.1]. In what follows, for all , we denote by its -dimensional Lebesgue measure and by its diameter. For all ; the distance between the cell center and the face center , where denotes the faces of an element not lying on . Finally, we set the vertices of , and the elements from sharing the vertex .

1.4 An Implicit Finite Volume Scheme

The cell centered finite volume scheme is a two-point flux approximation scheme, where the unknowns of the model are discretized using one value per cell: For all we let , with . System (2) is then discretized as follows: For all , and , we require(4)with denoting the upstream mobility and the two-point finite volume approximation of the normal component of the average phase velocity on is given by(5)where and is an approximation of the mass density of the phase on the face . Boundary fluxes are set to zero to account for the homogeneous natural boundary condition.

1.5 Linearization and Algebraic Resolution

At this stage, we need to solve, at each time step, the system of nonlinear algebraic equations resulting from the discretization method of Section 1.4. First, for we apply the Newton linearization algorithm generating, for fixed, a sequence with solution to the following system of linear algebraic equations: For all and all ,(6)with, , we have

Second, for , and a given Newton iteration , to approximate the solution of System (6), we use an iterative algebraic solver generating, for fixed, a sequence solving Equation (6) up to the residuals, given for all and all by(7)with linearized phase fluxes(8)

2 A posteriori error estimation

The spirit of the a posteriori error estimates is to monitor the computational error. In [10, 11] we derive fully computable a posteriori error estimates of the dual norm of the residual augmented by a nonconformity evaluation term for the thermal multiphase compositional general model. Moreover, in [8] the authors bound an energy-type norm of immiscible incompressible two-phase flow.

In the previously mentioned articles, similar types of estimators (more simple for our model) are obtained, whether we choose to bound the energy-type norm or the dual norm of residual. Therefore, we will use the notation for the error norm between the known numerical approximation and the unknown exact solution, without giving more precision. The local-in-time error norm is denoted then by at time step , Newton iteration , and linear solver iteration .

As mention in the introduction, our main goal is to propose a simple evaluation of the estimators proposed in [10] in order to be usable for engineering practices. Accordingly, we will directly apply these a posteriori error estimates on our model in order to obtain an a posteriori error estimate distinguishing the space, time, linearization, and algebraic errors. To this end, first we have to recall the required reconstructions necessary to get the upper bound, see [10] for more details.

2.1 Flux and Pressure Reconstructions

We recall here the necessary reconstructions for the theory of the a posteriori error estimates [10]. For these reconstructions we need the Raviart–Thomas–Nédélec (RTN) finite-dimensional subspaces of . Since the meshes , used in the numerical experiments of Section 5 below, consist of rectangular parallelepipeds then we define the RTN subspaces bysee Brezzi and Fortin [13] for details. We recollect all the formulas since they are necessary for the evaluation process of Section 3.

2.1.1 Phase Pressure Postprocessings

For all we define such that, for all and all ,(9)with defined by (5) and on . Then, for each , we introduce the piecewise quadratic phase pressure such that, for all ,(10)and The space–time function is then continuous and piecewise affine in time, given by at the discrete times . We need also phase pressure reconstructions, , such that for all . These reconstructions are typically piecewise polynomial continuous in space and piecewise affine continuous in time.

2.1.2 Component Flux Reconstructions

The discretization flux reconstruction is such that, for all and all ,(11a)

with defined by (5), while on . We also define a linearization error flux reconstruction such that, for all and for all ,(11b)with defined by (8). Similarly, an algebraic error flux reconstruction is defined such that, for all and for all ,(11c)with defined by (7). To complete both (11b) and (11c), we set respectively and on . The total flux reconstruction is then given by(11d)

2.2 Local-in-Time a Posteriori Error Estimate

Following [10], under the assumption of the existence and the uniqueness of a weak solution of the considered problem and the reconstructions of phase pressures and fluxes evoked in Section 2.1.2, the local-in-time error measure can be bounded as follows(12)with

The local version of the estimators is defined as, for the spatial estimators (13)with (14) (15)the temporal estimators (16)the linearization and algebraic estimators (17)

In order to implement the different estimators in a reservoir simulation code, we propose here to simplify them. Our simplification is based on the following observations:

• our estimators are expressed in terms of norms of various discrete quantities. Therefore we only need to compute these norms and not necessarily the quantities themselves;

• all these norms are integrals. Therefore they can be computed/approximated on each adequate quadrature formulas by mesh element;

• we thus only need to know the values of the various quantities at the quadrature points;

• the implementation of RTN spaces and the physical flux reconstructions can be consequently avoided.

Let us now describe our simplification process in all details.

3 Evaluation of the estimators using a practical simplified formula

The a posteriori estimates of Section 2.2 use the reconstruction of the different fluxes and in the space (Sect. 2.1.2). Furthermore, the pressure reconstructions also involve spaces (Sect. 2.1.1) in a way that for all ,

Hence, a priori, implementing these estimators requires operations with RTN spaces as seen in the previous chapter. To the best of our knowledge, so far RTN spaces are not implemented in petroleum industrial codes. Therefore we search for a simplification of these estimators that avoids the use of RTN spaces.

3.1 A General Simplification Formula

Let a mesh cell be given and let . Then can be expressed by(18)where are the degrees of freedom on the face , i.e., the face fluxes , and are the lowest-order Raviart–Thomas–Nédélec basis functions, given by (for parallelepiped meshes such as the ones used in the test cases of Sect. 5 below)where is the barycenter of the face σ′ opposite to the face , and(19)with,

Consider now the following quadrature formula exact for polynomials of total degree three, on a cubic domain :

Here , and . We can also rewrite this formula with vector symbols by(20)with , and as given in Table 1.

To apply the formula (20) to the integration of a function over a cell of our parallelepiped mesh, we consider the transformation such that,(21a) (21b) (21c)

Using the quadrature formula (20) and the transformations in (21), we can evaluate the norm of the function given by (18) as follows(22)

As a conclusion, to compute an norm of reconstructed flux functions in the space , we just need to obtain the degrees of freedom , represented by the face normal fluxes. These can be obtained directly without any flux reconstruction. In what follows we use these simplifications to evaluate the estimators for the immiscible two-phase model.

3.2 Evaluation of the Estimators

We begin with the spatial estimators. Recall the reconstructions of the conservative fluxes and given by Equation (11d). For , one has

Owing to the reconstruction of in Equation (11a), in (11b), and in (11c), we satisfy the conservation propertyup to a neglected maladjustment from the practical construction of (Corollary 3.4). Then, the residual estimators and given by Equation (14) will be neglected. Note that, however, if we construct the algebraic flux following Equation (11c) then the residual estimators and are zero.

Consider now the second contribution of the phase spatial estimator, given by Equation (13). Following Equation (3), the velocity is given by,(23)

Then by the post-processing of the oil pressure given by Equation (10) and by the definition (11a) of the reconstructed flux we havethen

Consequently using formula (22) with (9) and formula (11a), one has(24)

The non conformity estimators are also evaluated in the same way. These estimators involve the difference of two terms: and In order to apply formula (22), these two terms must be expressed in the space. This is the case of the first term (by construction, see Eq. 10), but not the case of the second . Therefore, we lift this term into the space by preserving the corresponding normal fluxes over the faces: Let be the piecewise linear nodal basis functions of on the three-dimensional element . Then for all the continuous pressure can be expressed as

To alleviate the industrial implementation, we then approximate the value of the pressure in the nodes by(25)where stands for the cardinality of the set ; it is equal to 8 for interior vertices in the three-dimensional case. But the values of are preferable since they are more precise. The normal fluxes of over the faces of can then be expressed explicitly by: for all (26)

Then using formula (22) with relation (10) we evaluate the nonconformity estimators given by Equation (14), as(27)where

This leads finally to a simple evaluation of the local spatial estimators for all . Note that a more precise formula could be obtained by replacing the piecewise linear nodal basis functions of by the piecewise quadratic basis functions of .

Summarizing the above developments, we have:

Corollary 3.1 (simple formula to evaluate the local spatial estimators). We can approximate the spatial estimators (13) as (28) where and are given by (24), (27), respectively.

Finally to evaluate the global spatial estimator we approximate the time integral by a one-dimensional integration formula.

Now for the temporal estimators, recall that for all , the reconstructions of the phase pressures , are such that are in the RTN space. Thus, the evaluation of the local temporal estimators can be again done using formula (22).

Corollary. 3.2 (simple formula to evaluate the temporal estimators). By applying formula (22) to Equation (16) and using the reconstructions formulations (11a)–(11c), we can approximate the temporal estimators as (29) where and are given by (21) and (19),

For the linearization estimators, we remark that the fluxes which compose the estimators are the reconstructed fluxes in the RTN space: , given by Equation (11b). Thus evaluating of these local estimators is straightforward using formula (22):

Corollary 3.3 (simple formula to evaluate the linearization estimators). By applying formula Equation (22) in Equation (17) and using the reconstructions formulations (11b), we can approximate the local linearization estimators as (30) where and are given by (21) and (19), respectively.

Finally, for the algebraic estimators, the evaluation of the local norms is also straightforward using formula (22), as the contributions given by (11c), involved in these norms, are in the RTN space.

Following [10, Remark 3.3], we arrive at the following approximation formula:

Corollary 3.4 (computing practically the algebraic error). Following [10, Remark 3.3], to compute approximately the algebraic error, we perform j additional iterations of the algebraic solver from the stage (7) where j is a user-defined fixed number. Then, the local algebraic error estimators can be evaluated using formula (22) as follow:(31a) where and are given by (21) and (19), respectively.

As a conclusion, the evaluation of the different estimators can be carried out while avoiding the physical reconstructions of fluxes in the space. The key is the use of a quadrature formula for computing the norms, knowing the normal fluxes over all faces of any cell . Using this technique greatly simplifies the implementation of the estimators (in particular into industrial legacy codes) and yields an important computational saving compared to the previous technique where we need to build the RTN flux reconstructions.

4 A posteriori error estimate to enhance the simulator performance

The presented a posteriori error estimate framework has been implemented in a parallel reservoir prototype simulator [14], a thermal multi-purpose simulator with the capability to simulate various EOR processes like water injection, steam injection, etc. It is part of the next generation IFP Energies Nouvelle research simulators, based on Arcane framework [15]. To improve the runtime performance of our simulator, we elaborated two strategies: the first consists in improving our linear solver by reducing the number of iterations using adaptative linear stopping criteria based on the a posteriori error estimate framework; the second strategy consists in reducing the global number of degrees of freedom by reducing the number of grid blocks of the mesh combining the Arcane adaptive mesh refinement feature to an advanced local space a posteriori error estimator.

4.1 An Adaptative Linear Solver Stopping Criteria

We have implemented an adaptative linear solver stopping criteria which allows by a call back mechanism to evaluate regularly the a posteriori linear error estimate described in the previous section. This mechanism allows to stop iterative linear solvers once the relative linear stopping criteria has been reached, or when the linear error estimation is lower than the non linear error estimation by a given factor. We can, in that way, reduce the number of linear solver steps during the steps of the Newton iterative non linear solver.

4.2 An a Posteriori Space Error Estimate Criteria to Manage the Adaptative Mesh Refinement Tools

In our simulator, we developed an adaptative mesh refinement mechanism that allows to refine or coarsen cells based on a space error criteria. This mechanism is aimed at adapting the grid block size to optimize the total number of degrees of freedom for a given global discrete error criteria. The idea is to have coarse grid blocks where the local error criteria is low, and fine grid block where the local error criteria is high. As the performance of the simulator is directly linked to the total number of grid blocks, minimizing this number has a direct impact on the global simulator runtime performance. We implemented a space error criteria based on the a posteriori error estimate tools described in the previous section. The idea is to reduce the standard deviation of a computed error field, used to choose the cells to refine and the cells to coarsen. The adaptative mesh refinement mechanism is then applied to minimize the required number of degrees of freedom for a given error criteria. During the coarsening and refinement phases, we use some upscaling and interpolating mechanisms to ensure the mass balance. Geological data (rock model properties like porosity, permeability, etc.) are stored at the finest geological level as data on a cloud of space points. These data are used, during mesh adaptation by standard upscaling algorithms to compute the values of the rock model properties on new created cells [16].

5 Numerical experimentation

We validate our approach on the tenth SPE comparative solution project model [17]. It is an incompressible water-oil, black-oil, two-phase flow problem. It is built on a Cartesian regular geometry with no top structure or faults and simulates a 3D waterflood in a geostatistical model. In this five-spot production model one water injection well is located at the center of the reservoir and four production wells at the four corners of the reservoir grid. Producers are controlled by a bottom hole pressure target. We present some results obtained on the fine grids comparing the runtime performance between the standard options of the simulator and the adaptative linear solver stopping criteria presented in Section 4.1. Others runs experiment the adaptative mesh refinement feature with the space criteria presented in Section 4.2 and results are compared against fine and coarse grid simulations. The coarse grid is obtained from the fine reference grid after two levels of coarsening applied on all cells. To validate the simulation results, we compared production curves for different wells to ensure that the results are still valid from a reservoir engineering point of view.

5.1 Platform Description

We have run our test cases using 16 cores of a dual socket linux workstation with 2 Sandy Bridge 10-cores Intel^® Xeon^® E5-2680 CPU with a tact rate of 2.8 GHz, 25600 KB cache size.

5.2 Stopping Criteria Results

To have a better control of the linear solver, we applied here a stopping criteria similar to that of the adaptive algorithm proposed in [10] with the main difference being that we avoid the physical flux reconstructions and the evaluation of the estimators is based on the simplification proposed in Section 3. We consider here a simulation of 200 days. In Figure 1, we compare the oil production rate from both the resolution with the standard options and then with the stopping criteria. In the left part of Figure 1 we depict cumulative oil production for the first two production wells. Cumulative oil production for the third and fourth production wells are then presented in the right part of Figure 1. We observe that the stopping criteria does not have any effect on the accuracy of the well productions.

Regarding the numerical behavior of the previous two runs we have had the following results:

• with standard options the number of resolutions related to the time steps is 67, the total number of non linear steps is 272, the total number of linear solver steps is 82 053, and the global CPU time of this simulation is 4 653 seconds;

• with the stopping criteria the number of resolutions related to the time steps is 67, the total number of non linear steps is 272, the total number of linear solver steps is 56 732, and the global CPU time of this simulation is 2 754 seconds.

As it is obvious, significant saving is obtained in terms of the cumulated number of linear solver iterations. We further illustrate in Figure 2 a comparison in terms of the consumed CPU time as well as the cumulated CPU time at several time steps of the simulation. The results show an over-all gain in CPU time, with a speed-up factor reaching the value 2.

5.3 AMR Results

We evaluate the performance of our simulator using the AMR feature combined with the a posteriori space error criteria. We compare in Figure 3 the oil production rate of the producer wells resulting from the resolution on the fine mesh, the coarse mesh with and without the AMR option. The fine grid results are supposed to be the reference. We observe that for all wells except well1 the AMR simulation behaves identically to the fine grid and much better than the simulation on the coarse grid.

In Figure 4 we compare the CPU time and the number of active cells regarding the simulation time. We observe a gain of performance for each time step for the same accuracy than the resolution on the fine mesh, due to a reduce number of required degrees of freedoms intermediate between the number of cells of the fine and the coarse mesh. We can notice some peaks of CPU time. This is due to the fact that the time step management is not optimal with the AMR feature, there are more failed time steps than in the reference run. We need more investigation to better understand this difference of behaviour in our simulator time step management.

Conclusion

With our new a posteriori error estimation framework, we have improved the performance of our compositional reservoir simulator dedicated to EOR processes. We have first optimized our non linear solver by implementing an adaptative linear stopping criteria to avoid to over-solve the linear systems in the first non linear steps. We have implemented secondly a space criteria to better manage the use of the AMR feature by optimizing the number of grid blocks required for a given accuracy of the well production results. Our approach has been validated on the SPE10 benchmark, a 3D study case representative of real-live reservoir study cases. Indeed, we have obtained significant gain in CPU times without affecting the accuracy of the oil production forecast. However, we still need to study the effect of the AMR feature on our dynamic time step management that we may improve in future works, implementing new criteria based on time error estimation.

References

All Tables

All Figures