Reservoir Simulator Runtime Enhancement Based on a Posteriori Error Estimation Techniques

— In this work, we show how the a posteriori error estimation techniques proposed in [Di Pietro et al. (2014) Computers & Mathematics with Applications 68 , 2331-2347] can be ef ﬁ ciently employed to improve the performance of a compositional reservoir simulator dedicated to Enhanced Oil Recovery (EOR) processes. This a posteriori error estimate allows to propose an adaptive mesh re ﬁ nement algorithm leading to signi ﬁ cant gain in terms of the number of cells in mesh compared to a ﬁ ne mesh resolution, and to formulate criteria for stopping the iterative algebraic solver and the iterative linearization solver without any loss of precision. The emphasis of this paper is on the computational cost of the error estimators. We introduce an ef ﬁ cient computation using a practical simpli ﬁ ed formula that can be easily implemented in a reservoir simulation code. Numerical results for a real-life reservoir engineering example in three dimensions show that we obtain a signi ﬁ cant gain in CPU times without affecting the accuracy of the oil production forecast. Résumé


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 2 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][3][4] suitable for unsteady nonlinear models, and [5][6][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.

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.

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 where P c o;w is the capillary pressure.Consequently, the unknowns of the model are the water pressure P w and the saturations S p , with p 2 fo; wg.These unknowns are collected in the vector X ¼ ðP w ; S w ; S o Þ.We denote by / the porosity of the medium and by K the permeability tensor.For a phase p 2 o; w f g; k r;p is the relative permeability, l p is the dynamic viscosity and q p is the density.The system of governing equations is given by with no-flow boundary conditions.For all p 2 o; w f g; m p represents the average phase velocity given by Darcy's law, where g denotes the gravity vector acting in the negative z direction and g its Euclidian norm.

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.

Discrete Setting
Let ðs n Þ 1 n N denote a sequence of positive real numbers corresponding to the discrete time steps such that t F ¼ P N n¼1 s n .We consider the discrete times ðt n Þ 0 n N such that t 0 :¼ 0 and, for 1 n N ; t n :¼ P n i¼1 s i ; then we define the time intervals I n :¼ ðt nÀ1 ; t n Þ.Let ðM n Þ 0 n N denote a family of meshes of the space domain X superadmissible in the sense of Eymard et al. [12,Definition 3.1].In what follows, for all M 2 M n , we denote by jM j its ddimensional Lebesgue measure and by h M its diameter.For all r 2 E i;n M ;d M ;r the distance between the cell center x M and the face center x r , where E i;n M denotes the faces of an element M 2 M n not lying on oX.Finally, we set V n M the vertices of M , and M n V the elements from M n sharing the vertex V .

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 1 n N we let 2) is then discretized as follows: For all 1 n N ; M 2 M n , and p 2 fw; og, we require with m " p ðX n M Þ denoting the upstream mobility and F p;M ;r ðX n M Þ the two-point finite volume approximation of the normal component of the average phase velocity on r is given by where a K :¼ K K d Kr 8K 2 fM; Lg; and q n p;r is an approximation of the mass density of the phase p on the face r.Boundary fluxes are set to zero to account for the homogeneous natural boundary condition.

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 1 n N we apply the Newton linearization algorithm generating, for X n;0 M fixed, a sequence ðX n;k M Þ k!1 with X n;k M solution to the following system of linear algebraic equations: For all p 2 fw; og and all M 2 M n , with, 8p 2 fw; og; 8M 2 M n , we have Second, for 1 n N, and a given Newton iteration k ! 1, to approximate the solution of System (6), we use an iterative algebraic solver generating, for X n;k;0 M fixed, a sequence ðX n;k;i M Þ i!1 solving Equation ( 6) up to the residuals, given for all p 2 fw; og and all M 2 M n by with linearized phase fluxes 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 N 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 N n;k;i at time step n, Newton iteration k ! 1, and linear solver iteration i ! 1.
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.

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 Hðdiv; XÞ.Since the meshes M n , used in the numerical experiments of Section 5 below, consist of rectangular parallelepipeds then we define the RTN subspaces by see Brezzi and Fortin [13] for details.We recollect all the formulas since they are necessary for the evaluation process of Section 3.

Phase Pressure Postprocessings
For all p 2 fw; og we define with F p;M ;r defined by ( 5) and C n;k;i p;h Á n X ¼ 0 on oX.Then, for each p 2 fw; og, we introduce the piecewise quadratic phase pressure P n;k;i p;h such that, for all M 2 M n , and ðP n;k;i p;h ;1Þ M jM j ¼ P n;k;i p;M : The space-time function P n;k;i p;hs is then continuous and piecewise affine in time, given by P n;k;i p;h at the discrete times t n .We need also phase pressure reconstructions, P p;hs; p 2 fw; og, such that P p;hs; 2 X for all p 2 fw; og.These reconstructions are typically piecewise polynomial continuous in space and piecewise affine continuous in time.

Component Flux Reconstructions
The discretization flux reconstruction H n;k;i dis;p;h 2 RTNðM n Þ is such that, for all M 2 M n and all r 2 E i;n M , with F p;M ;r defined by (5), while H n;k;i dis;p;h Á n X; ¼ 0 on oX.We also define a linearization error flux reconstruction H n;k;i lin;p;h 2 RTNðM n Þ such that, for all M 2 M n and for all r 2 E i;n M , with F n;k;i p;M ;r defined by (8).Similarly, an algebraic error flux reconstruction H n;k;i alg;p;h 2 RTNðM n Þ is defined such that, for all M 2 M n and for all r 2 E i;n M , with R n;k;i p;M defined by (7).To complete both (11b) and (11c), we set respectively H n;k;i lin;p;h Á n X ¼ 0 and H n;k;i alg;p;h Á n X ¼ 0 on oX.The total flux reconstruction H n p;h is then given by 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 The local version of the estimators is defined as, for p 2 fw; og the spatial estimators Àr ÁH 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.

EVALUATION OF THE ESTIMATORS USING A PRACTICAL SIMPLIFIED FORMULA
The a posteriori estimates of Section 2. 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.

A General Simplification Formula
Let a mesh cell M be given and let D h 2 RTNðM Þ.Then D h can be expressed by where c r are the degrees of freedom on the face r, i.e., the face fluxes ðD h Án M ; 1Þ r , and K r 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 x r 0 ¼ ðx r 0 ; y r 0 ; z r 0 Þ is the barycenter of the face r 0 opposite to the face r, and with, e w r;r 0 :¼ " 1 w r ;w r 0 ; w 2 x; y; z f g: " otherwise Consider now the following quadrature formula exact for polynomials of total degree three, on a cubic domain K ¼ ðÀ1; 1Þ 3 : We can also rewrite this formula with vector symbols by with W k ¼ 1 8 ; k 2 f1; 2; ::; 8g, and P k ¼ ðP k x ; P k y ; P k z Þ as given in Table 1.
To apply the formula (20) to the integration of a function over a cell Using the quadrature formula (20) and the transformations in (21), we can evaluate the ½L 2 ðM Þ 3 norm of the function Á h 2 RTNðM Þ given by (18) as follows As a conclusion, to compute an ½L 2 3 norm of reconstructed flux functions in the space RTNðM n Þ, we just need to obtain the degrees of freedom c r , 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.

Evaluation of the Estimators
We begin with the spatial estimators.Recall the reconstructions of the conservative fluxes H n;k;i w;h and H n;k;i o;h given by Equation (11d).For H n;k;i p;h ; p 2 fw; og, one has Owing to the reconstruction of H n;k;i dis; p;h in Equation (11a), H n;k;i lin; p;h in (11b), and H n;k;i alg; p;h in (11c), we satisfy the conservation property up to a neglected maladjustment from the practical construction of H n;k;i alg; p;h (Corollary 3.4).Then, the residual estimators g n;k;i R;M ;p 0 and g n;k;i R;M ;o given by Equation ( 14) will be neglected.Note that, however, if we construct the algebraic flux H n;k;i alg; p;h following Equation (11c) then the residual estimators g n;k;i R;M ;w ; and g n;k;i R;M ;o are zero.
Consider now the second contribution of the phase spatial estimator, given by Equation (13).Following Equation (3), the velocity m p ðP n;k;i p;hs Þ is given by, m p ðP n;k;i p;hs Þ ¼ ÀK rP n;k;i p;h À q p ðP n;k;i p;h Þg ð23Þ Then by the post-processing of the oil pressure P n;k;i o;h given by Equation (10)  Consequently using formula ( 22) with ( 9) and formula (11a), one has The non conformity estimators are also evaluated in the same way.These estimators involve the difference of two terms: ðÀKrP n;k;i p;hs Þj M and ÀKrP n;k;i p;hs j M : In order to apply formula (22), these two terms must be expressed in the RTNðM n Þ space.This is the case of the first term ðÀKrP n;k;i p;hs Þj M (by construction, see Eq. 10), but not the case of the second ðÀKrP n;k;i p;hs Þj M .Therefore, we lift this term ðÀKrP p;hs Þj M into the RTNðM n Þ space by preserving the corresponding normal fluxes over the faces: Let fu M i g i2f1;2;::;8g be the piecewise linear nodal basis functions of Q 1 on the three-dimensional element M .Then for all M 2 M n the continuous pressure P n;k;i p;hs can be expressed as To alleviate the industrial implementation, we then approximate the value of the pressure in the nodes by where N M n V stands for the cardinality of the set M n V ; it is equal to 8 for interior vertices in the three-dimensional case.
But the values of P n;k;i p;hs are preferable since they are more precise.The normal fluxes of ðÀKrP n;k;i p;hs Þj M over the faces of M can then be expressed explicitly by: for all r 2 E i;n M ðÀKrP n;k;i p;hs Á n X ; 1Þ r ¼ Then using formula (22) with relation (10) we evaluate the nonconformity estimators g n;k;i NC;M ;p ; given by Equation ( 14), as where This leads finally to a simple evaluation of the local spatial estimators g n;k;i sp;M ;w ; g n;k;i sp;M ;o for all M 2 M n .Note that a more precise formula could be obtained by replacing the piecewise linear nodal basis functions of Q 1 by the piecewise quadratic basis functions of Q 2 .
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 where g n;k;i F;M ;p and g n;k;i NC;M ;p 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 M 2 M n , the reconstructions of the phase pressures P n;k;i p;hs j M ; p 2 fw; og, are such that KrP n;k;i p;hs j M are in the RTNðM n Þ 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 where f k and E r;r 0 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ðM n Þ space: H n;k;i lin; p;h ; p 2 fw; og, 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 where f k and E r;r 0 are given by (21) and (19), respectively.Finally, for the algebraic estimators, the evaluation of the local ½L 2 3 norms is also straightforward using formula (22), as the contributions H n;k;i alg; p;h ; p 2 fw; og given by (11c), involved in these norms, are in the RTNðM n Þ 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: where f k and E r;r 0 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 RTNðM n Þ space.The key is the use of a quadrature formula for computing the L 2 norms, knowing the normal fluxes over all faces of any cell M 2 M n .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.

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.

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.

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].

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.

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.

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. Cumulative oil production of wells during the simulation.

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.

Figure 3
Figure 3Cumulative oil production of wells during the simulation.
Figure 4 CPU time (left) and active cell number (right) of the simulation.

Figure 2 CPU
Figure 2 CPU time (left) and total CPU time (right) of the simulation.