Adaptive Algebraic Multiscale Solver for Compressible Flow in Heterogeneous Porous Media

(1)

th

Mo A10

Adaptive Algebraic Multiscale Solver for

Compressible Flow in Heterogeneous Porous

Media

M. Ţene* (Delft University of Technology), H. Hajibeygi (Delft University of Technology), Y. Wang (Stanford University) & H.A. Tchelepi (Stanford University)

SUMMARY

An adaptive Algebraic Multiscale Solver for Compressible (C-AMS) flow in heterogeneous oil reservoirs is developed. Based on the recently developed AMS [Wang et al., 2014] for incompressible linear flows, the C-AMS extends the algebraic formulation of the multiscale methods for compressible (nonlinear) flows. Several types of basis functions (incompressible and compressible with and without accumulation) are considered to construct the prolongation operator. As for the restriction operator, C-AMS allows for both MSFV and MSFE methods. Furthermore, to resolve high-frequency errors, Correction Functions and ILU(0) are considered. The best C-AMS procedure is determined among these various strategies, on the basis of the CPU time for three-dimensional heterogeneous problems. The C-AMS is adaptive in all aspects of prolongation, restriction, and conservative reconstruction operators for time-dependent compressible flow problems. In addition, it is also adaptive in terms of linear-system update. Though the C-AMS is a conservative multiscale solver (i.e., only a few iterations are employed infrequently in order to maintain high-quality results), a benchmark study is performed to investigate its efficiency against an industrial-grade Algebraic Multigrid (AMG) solver, SAMG [Stuben, 2010]. This comparative study illustrates that the C-AMS is quite efficient for compressible flow simulations in large-scale heterogeneous 3D reservoirs.

(2)

Introduction

Accurate and efficient simulations of multiphase flow in large-scale highly heterogeneous formations has motivated the development of a multiscale (MS) family of numerical methods, e.g., Multiscale Finite Element (MSFE) (Hou and Wu (1997); Efendiev and Hou (2009)) and Finite Volume (MSFV) (Jenny et al. (2003, 2006)). The MS methods map a discrete fine-scale system to a much coarser space by using locally computed (and adaptively updated) basis functions. In the multigrid terminology (Trottenberg et al. (2001)), this map is considered as constructing a special prolongation operator, having locally-supported basis functions in its columns (Zhou and Tchelepi (2008)). The restriction operator is then defined based on a Finite Element (MSFE), Finite Volume (MSFV), or a combination of both (Cortinovis and Jenny (2014)).

The MSFV has been applied to a wide range of applications (See, e.g., Zhou and Tchelepi (2008, 2012); Lunati and Jenny (2008); Lee et al. (2009); Hajibeygi et al. (2008, 2011, 2012); Wolfsteiner et al. (2006); Lunati et al. (2011); Kunze et al. (2013); Tomin and Lunati (2013); Hajibeygi and Tchelepi (2014); Wang et al. (2014); Moyner and Lie (2014)). The combined developments have made it a promising framework for the next-generation reservoir simulators.

Most MS developments have focused on the incompressible problem. However, when compressibility effects are considered, the pressure equation becomes nonlinear, and its solution requires an iterative procedure involving a parabolic-type system of linear equations (Aziz and Settari (1979)).

The present study introduces an algebraic multiscale iterative solver for compressible flows in heteroge-nous porous media (C-AMS), along with a thorough study of its computational efficiency (based on CPU and number of iterations). In order to develop an efficient prolongation operator, several types of basis functions are considered. These basis functions differ from each other in the amount of compressibility involved in their formulation, ranging from incompressible elliptic (Lunati and Jenny (2006)), compress-ible elliptic (Zhou and Tchelepi (2008); Lee et al. (2008)), and a pressure-independent parabolic (incom-pressible advection with accumulation) method (Hajibeygi and Jenny (2009)). As for construction of an efficient restriction operator, both MSFE and MSFV operators are considered. For the second-stage solver, both the Correction Functions and ILU(0) methods are employed. The best C-AMS procedure is developed among these various strategies, on the basis of the CPU time for three-dimensional (3D) heterogeneous problems.

Though the C-AMS is a conservative multiscale solver (i.e., only a few iterations are performed in-frequently in order to maintain high-quality results), a benchmark study is performed to investigate its efficiency against an industrial-grade Algebraic Multigrid (AMG) solver, SAMG (Stuben (2010)). This comparative study for compressible problems is the first of its kind, and is made possible through our algebraic formulation. Numerical results illustrate that the C-AMS is quite efficient for the simulation of nonlinear compressible flow problems.

The paper is structured as follows. First, the Compressible Algebraic Multiscale Solver (C-AMS) is developed, where several options for prolongation and restriction operators as well as the second-stage solver are considered. Numerical results are then presented for a wide range of 3D heterogeneous test cases, aimed at determining the optimum strategy for the treatment of the heterogeneous transmissibility ﬁeld and nonlinear compressibility effects. Adaptive updating of C-AMS operators, and, importantly, the linear and nonlinear system of equations are also discussed. Finally, the C-AMS is compared with an industrial-grade Algebraic Multigrid Solver (SAMG) both in terms of the number of iterations and overall CPU time.

(3)

th

Compresible Flow in Heterogeneous Porous Media

Mass conservation for phaseα using Darcy’s law (without gravity and capillarity effects) can be written as:

∂

∂t(φραSα) − ∇ ·

ραλλλα∇p= ραqα ∀α ∈ {1,...,nph}, (1)

whereφ is the porosity, ρ density, S saturation, and q the source terms. Moreover, λλλ_α= KKKkrα/μα is the phase mobility with positive-deﬁnite permeability tensor, KKK, the relative permeability, kr_α, the phase viscosity, μ_α, and nph is the number of phases (and components). Note that mass exchange between phases (compositional effects) is neglected here.

After implicit (Euler-backward) time integration, the discrete equations are divided by the density at the new time step (n+ 1), then summed over all phases in order to obtain a pressure equation. This is achieved by using the constraint∑nph

α=1Snα+1= 1, leading to the nonlinear elliptic ﬂow problem φn+1 Δt − φn Δt nph

∑

α=1 ρn αSnα ρn+1 α − nph

∑

α=1 1 ρn+1 α ∇ · ρn+1 α λλλα∇pn+1= q, (2)

which is linearized as follows:

cν(pν+1− pν) − nph

∑

α=1 1 ρν α∇ · (ρ ν αλλλα· ∇pν+1) = bν, (3) where cν= 1 Δt ∂φ ∂ p ⏐ ⏐ ⏐ ⏐ ν − φn nph

∑

α=1 ∂ ∂ p 1 ρα ⏐_⏐ ⏐ ⏐ ν ρn αSnα (4) and bν= −φ ν Δt + φn Δt nph

∑

α=1 1 ρ_ανραnSnα+ q. (5)

The superscriptsν and ν + 1 denote the old and new iteration levels, respectively. Note that as ν → ∞, (3) converges to (2) and Cνpν+1− pν→ 0. The pressure-dependent coefﬁcient, c, is a by-product of the linearization lemma and only plays a role during iterations. This fact opens up the possibility to alter c, by computing it based on either pν (resulting in cν), or pn (corresponding to cn) - the pressure at the previous time-step. Each choice potentially can lead to a different convergence behavior.

Algebraically, (3) can be written as

Cν+ ˜Aνpν+1≡ Aνpν+1= rν≡ bν+Cνpν, (6)

where Cν is a diagonal matrix containing the values cν, and ˜A is the compressible ﬂow matrix. Compressible Algebraic Multiscale Solver (C-AMS)

Similar to AMS (Wang et al. (2014)), the C-AMS relies on the primal- and dual-coarse grids imposed on the ﬁne-scale grid (See Figure 1). There are Npand Nd coarse and dual-coarse grid cells in the domain of Nf ﬁne-grid cells.

The transfer operators between ﬁne-scale and coarse-scale are deﬁned using multiscale Restriction (R) and Prolongation (P). The former depends on the choice of discretization: either Finite Elements

(4)

Figure 1: Multiscale grids imposed on the given ﬁne-grid. Primal- and dual-coarse block are delimited using black lines and highlighted blue cells, respectively.The central primal-coarse block is drawn with red to show the overlap between the two coarse grids.

(MSFE), in which case R= PT, or Finite Volumes (MSFV), for which R corresponds to the volume integral over each primal-coarse block, i.e.,

R[i, j] =

dV , if ﬁne-cell j is contained in primal-coarse block i

0 , otherwise. (7)

The columns of P are populated using the values of the basis functions, which are computed on dual-coarse cells (see Figure 1), with simplified boundary conditions. In contrast to (incompressible) AMS, C-AMS can be formulated based on different choices of the basis functions, depending on the level of compressibility involved. The first two types, which are consistent with Eq. (3), are defined as

cnΦν+1_k_,h − nph

∑

α=1 1 ρν α∇(ρ ν αλλλα· ∇Φν+1k,h ) = 0 (8) − nph

∑

α=1 1 ρν α∇ · (ρ ν αλλλα· ∇Φν+1k_,h ) = 0, (9)

which change with pressure (through ρ_α), and their left-hand-side differs from that of the ﬁne-scale system only in the accumulation term, c. Alternatively, one can also change the advection term, leading to

cnΦk,h− ∇ · (λλλt· ∇Φk,h) = 0 (10)

−∇ · (λλλt· ∇Φk,h) = 0, (11)

which are both pressure independent, since cn is calculated based on pn. Here, λt = ∑nph

α=1λα. All of these equations are subject to reduced-problem boundary conditions along dual-coarse cell boundaries ∂ ˜Ωh_{(See, e.g., Jenny et al. (2006)). One can also obtain the equations for the corresponding four types of} correction functions,Ψh, by substituting the corresponding right-hand-side (RHS) in Eqs. (8)-(11). Note that Hajibeygi and Jenny (2009) used only one possible choice for their basis and correction functions. The basis functionsΦkare assembled over dual-coarse cells ˜Ωh ∀h ∈ {1,...,Nd}, i.e., Φk= ∑N_h₌₁d Φk_,h, and the correction functions can be assembled asΨ = ∑Nd

h₌₁Ψνh. Figure 2 illustrates two possible basis functions, from which it is clear that the presence, or absence, of the accumulation term plays a crucial role in determining the values of the (interpolation) basis functions. In particular, the basis functions do not add up to one if they account for compressibility effects.

(5)

th

The choices formulated above dictate how much computational effort is invested in constructing and updating the multiscale operators, i.e., while basis functions of Eqs. (8) and (9) depend on pressure (hence, updated adaptively when pressure changes), Eqs. (10) and (11) are pressure independent; thus, they do not get updated when the pressure changes. In choosing one over the other, one needs to consider the dynamics of the problem at hand. While the basis and correction functions from Eq. (10) were previously used by Hajibeygi and Jenny (2009), the other options are, as of yet, unexplored.

cnΦν+1− ∑Nph α=1ρ1_αν∇.(ρ_ανλα∇Φν+1) = 0 0 0.5 1 0 0.5 1 −∇.(λt∇Φ) = 0 0 0.5 1 0 0.5 1

Figure 2: Two choices of multiscale basis functions in a reference dual-coarse block (left), Summation of the basis functions over the dual-coarse block (right), i.e., partition of unity check.

The C-AMS approximates the ﬁne-scale solution pνby pν, i.e., Aνpν≈ Aνpν. Once the basis functions are computed and assembled, the C-AMS Prolongation operator PPP is then constructed, which is a matrix of size Nf × Nc, having basis functionΦk in its k− th column. Then, the map between the coarse ( ˆp) and ﬁne (p) solution reads

p= PPP ˆp. (12)

The coarse-scale system can be obtained by restricting the ﬁne-scale system (3) as

Ac,νpc,ν+1≡ (RRRAνPPP) ˆpν+1= RRR rν, (13) and its solution is prolonged to the ﬁne-scale using (12), i.e.,

pν+1= PPP(RRRAνPPP)−1RRR rν, (14) which in residual form reads as

(6)

Here pν+1= pν+ δ pν+1 andℜν = rν− Aνpν is the ﬁne-scale residual. Note that all different op-tions for basis funcop-tions are considered here to construct the prolongation operator, in order to ﬁnd the optimum procedure.

The C-AMS employs Eq. (14) as the global solver (for resolving low-frequency errors). In addition to the coarse-scale solver, an efficient convergent multiscale solver needs to include a second-stage solver at the fine-scale, in order to circumvent the simplified localization boundary conditions, the nonlinearity of the operator, and the complex RHS term. Among the choices for the second stage solver (block-, line-, or point-wise solvers)line-, the correction functions (CF) and ILU(0) are used in this work. The C-AMS procedure is summarized in Table 1.

Do until convergence (εν< e) achieved (See Eq. (17)) {

0. Update parameters, linear system matrix, and right-hand-side vector based on pν 1. Adaptively compute the Basis Functions using one of Eqs. (8)-(11)

2. Calculate residual:ℜν= rν− Aνpν

3. Pre-smoothing Stage (only when CF is used): Apply CF system onℜν and update Residual 4. Multiscale Stage: Solve (15) forδ pν+1/2

5. Post-smoothing Stage: iterate using the chosen smoother to improveδ pν+1/2, obtainingδ pν+1 6. pν+1= pν+1/2+ δ pν+1

}

Table 1: Outline of the C-AMS iteration procedure Numerical Results

The numerical experiments presented in this section are restricted to single-phase compressible ﬂuid ﬂow through the set of ‘patchy’ permeability realizations from Wang et al. (2014) (one of which is shown in Figure 3a), by imposing Dirichlet conditions pwest = 1 and peast = 0 on the west and east boundaries, respectively. For time-stepping, we use the non-dimensional time

t∗= μφL

2 ¯

K(pwest− peast),

(16) where ¯K is the average ﬁeld permeability, and L is the length of the domain (in our case, a cube with regular ﬁne-cells of size 1× 1 × 1). Figure 3b depicts the solution on a 64 × 64 × 64 grid-cell test case after 0.4t∗.

All implementations of this paper were designed to be single-threaded, and the CPU times were mea-sured on an Intel Xeon E5-1620 v2 quad-core system with 64GB RAM.

Convergence study

For our compressible multiscale solver (C-AMS), there are two ways to measure the error in our approx-imate solution. One approach is to use it in the nonlinear pressure Eq. (2), i.e.,

εν = q −φ_Δtn+1+φ_Δtn nph

∑

α=1 ρn αSnα ρn+1 α + nph

∑

α=1 1 ρn+1 α ∇ · ρn+1 α λλλα∇pn+1 (17)

and we will refer to εν 2 as the error norm, or in the linearized equation (3), which is equivalent to computing the residual norm, ℜν 2.

In order to determine a suitable stopping criterion, we considered a patchy domain with 64× 64 × 64 grid cells, for which the pressure equation is solved using the following solution strategy:

The error and residual norms were recorded after each iteration of the Richardson loop and are given in Figure 4a. Notice that the reduction of the residual norm beyond the ﬁrst few iterations does not

(7)

th

(a) log(permeability) (b) Pressure after 0.4t∗ Figure 3: One realization of a 64× 64 × 64 patchy permeability ﬁeld Do until ε 2< 10−6{

0. Update parameters, linear system matrix and right-hand-side vector based on pν

1. Solve linear system using a Richardson iterative scheme, preconditioned with one multigrid V-cycle until ℜ 2< 10−6

}

Table 2: Solution strategy used to determine a suitable stopping criterion

contribute to the reduction in the error norm. Therefore, one could speed up the solution scheme by monitoring the error norm and updating the linear system after its decrease gives stagnates. However, the computational cost of evaluating the nonlinear equation is roughly the same as that of a linear system update and, thus, orders of magnitude more than the evaluation of the residual norm.

0 20 40 60 80 100 120 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 101 102 Iteration Error norm Residual norm

(a) Full-solve of the linearized system

0 5 10 15 20 25 30 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 101 102 Iteration Error norm Residual norm

(b) Adaptive solve of the linearized system Figure 4: Error and residual norm histories for the 643grid-cell problem over a single time step of 0.4t∗. Figure 4a also reveals that the stagnation of the error norm seems to happen roughly after the residual norm has been reduced by 1/10 of its initial value (i.e. immediately after the linear system update). Figure 4b shows the convergence behaviour after implementing this heuristic. The two norms seem to be roughly in agreement; hence, in the following experiments we will evaluate ℜ

i 2 ℜ0

2 < 10 −1_after iteration i of the inner loop and εν 2< 10−6after iterationν of the outer loop. (see Table 2).

(8)

Adaptive updating of data structures

The previous study described the ﬁrst adaptivity method considered in this work, namely, updating the linear system only after the residual norm drops by an order of magnitude. This adaptive strategy leads to a 45% speedup, as shown in Fig. 5b. The basis functions need to be recomputed when the transmissibilities change signiﬁcantly. In addition to this, for the correction functions, one also needs to also monitor the corresponding right-hand-side terms. Figure 5c shows that adaptive update of the C-AMS basis functions lead to 85% speedup in terms of CPU time (Figure 5a).

Furthermore, the two adaptivity methods can be applied together, for maximum beneﬁt, as can be seen in Figure 5d. 0 5 10 15 20 25 0 5 10 15 20 25 Iteration

CPU time (sec)

515.8598 sec (a) No adaptivity 0 5 10 15 20 25 0 5 10 15 20 25 Iteration

CPU time (sec)

279.5869 sec

(b) Linear system update adaptivity

0 5 10 15 20 25 0 5 10 15 20 25 Iteration

CPU time (sec)

83.1234 sec

(c) Multiscale operator update adaptivity

0 5 10 15 20 25 0 5 10 15 20 25 Iteration

CPU time (sec)

74.2812 sec

(d) Both methods together

Multiscale solution Smoother solution Lin. sys. construction Basis functions Correction function

Figure 5: Effect of different types of adaptivity on the C-AMS performance for a 643grid-cell problem after a time step of 0.4t∗.

Local stage: choice of basis functions

The aim of this study is to determine the best choice of basis function to be used in the C-AMS algorithm. For this purpose we perform simulations on 20 realizations of a patchy reservoir with 64× 64 × 64 grid cells. The correction function is computed based on Eq. (8) in all cases (and, hence, updated adaptively with pressure), 20 iterations of ILU(0) are used for smoothing and all possibilities for the basis functions, i.e., Eqs. (8)-(11). Finally, there is a single time step in the simulation, which takes the initial solution at time 0 (p0= 0 everywhere) to time 0.4t∗.

(9)

th

top of each bar in Fig. 6), are measured. Results show that frequent updating of the basis functions does not impact the solution signiﬁcantly; hence, the best choice is to use the incompressible pressure independent basis functions, i.e., Eqs. (10) and (11). Also, the inclusion of the accumulation term or the type of Restriction (MSFE or MSFV) does not play an important role for this test case. Therefore, in the following experiments the “Incomp:FVM” is chosen as the C-AMS strategy for iterations, (i.e., FV-based restriction with incompressible basis functions). Note that for this patchy domain the performance of the “Incomp:FEM” (i.e., MSFE as the coarse-scale solver) and “Incomp:FVM” are comparable.

0 20 40 60 80

Comp+Accum:FVMComp+Accum:FEMComp:FVM Comp:FEM Incomp+Accum:FVMIncomp+Accum:FEMIncomp:FVM Incomp:FEM Average CPU time of CF + MS + 20 ILU with different types of Prolongation and Restriction

CPU time (sec)

29 30 34 20

30 39 34 21

Figure 6: Effect of the choice of basis function on the C-AMS performance for the 643grid-cell problem after a time step of 0.4t∗. Results are averaged over 20 statistically-the-same realizations. The number of iterations is shown on top of each bar.

Global stage: choice of correction function

After the specific choice for evaluating the basis functions, the same experiments are repeated for the correction functions. The CF is an independent stage; therefore, the possibility of eliminating it from the C-AMS algorithm is considered. The results in Fig. 7 confirm that eliminating it leads to overall speedup. Also Fig. 8 shows that the use of correction functions may lead to difficulty in convergence, a disadvantage that also supports eliminating it.

Global stage: number of smoother iterations

A ﬁnal variable in the C-AMS framework is the number of smoothing steps (ILU(0)) that should be applied in order to obtain the best tradeoff between convergence rate and CPU time. The results of our repeated experiments with the optimum choices found so far (MSFV with incompressible basis functions and without correction) and various numbers of ILU applications are illustrated in Fig. 9. It is clear that the optimum scenario lies between 5-10 ILU iterations per second-stage call. Note that all C-AMS runs (without correction functions) converged successfully.

CPU benchmark for different grid sizes

The results discussed so far indicate that for the patchy test case considered here C-AMS is optimum when MSFV (or MSFE) restriction operators with incompressible basis functions and 5 iterations of ILU are considered. In this subsection, the computational efﬁciency of C-AMS against the SAMG industrial-grad algebraic multigrid solver (Stuben (2010)), is presented for the patchy permeability ﬁelds

(10)

0 5 10 15 20 25 30 35

Comp+Accum:FVMComp+Accum:FEMComp:FVM Comp:FEM Incomp+Accum:FVMIncomp+Accum:FEMIncomp:FVMIncomp:FEMNo correction:FVMNo correction:FEM

Average CPU time of CF + MS + 20 ILU with different types of Correction and Restriction

CPU time (sec)

34 21 36 31 32 21 35 31 20 20 Multiscale solution Smoother solution Lin. sys. construction Basis functions Correction function

Figure 7: Effect of the choice of correction function on the CPU time of the multiscale solution on a set of 20 realizations of a 643grid-cell patchy reservoir after a time step of 0.4t∗. The number of iterations is shown on top of each bar. Only the last 2 bars on the right correspond to runs in which no correction function was used (i.e., MS + ILU).

0 20 40 60 80 100

Comp+Accum:FVMComp+Accum:FEMComp:FVMComp:FEMIncomp+Accum:FVMIncomp+Accum:FEMIncomp:FVMIncomp:FEMNo correction:FVMNo correction:FEM

Success rate of CF + MS + 20 ILU with different types of Correction and Restriction

% Runs

Converged Diverged

Figure 8: Effect of the choice of correction function on the success rate of the solver on a set of 20 realizations of the 643 grid-cell patchy reservoir after a time step of 0.4t∗. Only the last 2 bars on the right correspond to runs in which no correction function was used (i.e., MS + ILU).

(a realization is shown in Fig. 3). Figure 10 shows a typical solution after non-dimensioanl times of 0.4, 1.0, and 2.

Figures (11) and (12) show the number of iterations and CPU time at 0.4, 1.0, and 2 and different problem sizes of 64× 64 × 64, 128 × 128 × 128 and 256 × 256 × 256 ﬁne cells, respectively. Except for the ﬁrst time-step, when all the basis functions are fully computed, C-AMS has a slight edge over SAMG, mainly due to its adaptivity. The initialization cost of C-AMS is particularly high in the 2563 case, due to the large number of linear systems (solved with a direct solver) needed for the basis functions. By comparing Figs. 12(left) and 12(right) where the coarsening factor was changed from 8× 8 × 8 to 16×16×16, it is clear that, with larger primal-coarse blocks, C-AMS requires less setup time, but more

(11)

th 1 5 10 15 20 0 1 2 3 4 5 6 7 8 9 # Smoothing steps

CPU time (sec)

Average CPU time of MS + ILU with different number of Smoothing steps

217

51 30 23

20

Figure 9: Effect of the number of smoothing steps on the C-AMS performance (MSFV + ILU) for a set of 20 realizations of the 643grid-cell patchy reservoir after a time step of 0.4t∗. The number of iterations is shown on top of each bar.

iterations to converge.

Since reservoir simulators are typically run for many time-steps, the high initialization time of C-AMS is outweighed by the efﬁciency gained in subsequent steps.

(a) Pressure after 0.4t∗ (b) Pressure after 1.0t∗ (c) Pressure after 2.0t∗ Figure 10: Pressure solution on a realization of the 64× 64 × 64 patchy permeability ﬁeld Conclusions

In this paper we introduced the Algebraic Multiscale Solver for Compressible ﬂows (C-AMS) in het-erogeneous porous media. The new solver beneﬁts from adaptivity, both in terms of the infrequent updating of the linearized system and from the selective update of the basis functions used to construct the prolongation operator.

Our experiments on heterogeneous patchy reservoirs revealed that the most efﬁcient strategy is to use basis functions with incompressible advection terms, paired with 5 iterations of ILU(0).

Finally, in our benchmark study, we compared our algorithm to an industrial-grade multigrid solver (SAMG), and the results show that C-AMS is a competitive solver, especially in experiments that involve

(12)

AMS SAMG AMS SAMG AMS SAMG 0 1 2 3 4 5 6

CPU time (sec)

CPU bencmark on a 64 x 64 x 64 grid−cell reservoir

49 21 44 21 41 23 0.0 − 0.4t∗ _{0.4 − 1.0t}∗ _{1.0 − 2.0t}∗ Initialization Lin. sys. construction Solution

AMS SAMG AMS SAMG AMS SAMG

0 10 20 30 40 50 60 70

CPU time (sec)

Figure 11: CPU time comparison between the MSFV + 5 ILU and SAMG solvers on a 105 (left) and 106(right) grid-cell patchy reservoir over 3 succesive time-steps. The coarsening scheme of C-AMS is 8× 8 × 8 ﬁne-cells per coarse block. The number of iterations is given on top of each bar.

0 200 400 600 800 1000 1200 1400 1600 1800 2000

CPU time (sec)

33

32 28 32 26 32

0.0 − 0.4t∗ _{0.4 − 1.0t}∗ _{1.0 − 2.0t}∗

Initialization Lin. sys. construction Solution

0 100 200 300 400 500 600

CPU time (sec)

Figure 12: CPU time comparison between the MSFV + 5 ILU and SAMG solvers on a 107 grid-cell patchy reservoir over 3 succesive time-steps. The coarsening scheme of C-AMS is 8× 8 × 8 (left) and 16× 16 × 16 (right) ﬁne-cells per coarse block. The number of iterations is given on top of each bar. the simulation of a large number of time steps. The only drawback is the relatively high initialization time, which can be reduced by choosing an appropriate coarsening strategy, or, by running the basis function update on a parallel platform. With C-AMS only a few iterations per time step are required to obtain good quality results.

Acknowledgements

We appreciate ﬁnancial support from Chevron/Schlumberger Intersect Alliance Technology during Matei ¸

Tene’s scientiﬁc visit at TU Delft between November, 2013 - February, 2014. Since March 2014, Matei ¸

Tene is a PhD student at TU Delft sponsored by PI/ADNOC. References

Aziz, K. and Settari, A. [1979] Petroleum Reservoir Simulation. Blitzprint Ltd., Cagary, Alberta.

Cortinovis, D. and Jenny, P. [2014] Iterative galerkin-enriched multiscale ﬁnite-volume method. J. Comp. Phys., under review.

Efendiev, Y. and Hou, T.Y. [2009] Multiscale Finite Element Methods: Theory and Applications. Springer. Hajibeygi, H., Bonﬁgli, G., Hesse, M. and Jenny, P. [2008] Iterative multiscale ﬁnite-volume method. J. Comput.

Phys., 227, 8604–8621.

Hajibeygi, H. and Jenny, P. [2009] Multiscale ﬁnite-volume method for parabolic problems arising from com-pressible multiphase ﬂow in porous media. J. Comput. Phys., 228, 5129–5147.

(13)

th

Hajibeygi, H., Karvounis, D. and Jenny, P. [2011] A hierarchical fracture model for the iterative multiscale ﬁnite volume method. J. Comput. Phys., 230(24), 8729–8743.

Hajibeygi, H., Lee, S.H. and Lunati, I. [2012] Accurate and efficient simulation of multiphase flow in a hetero-geneous reservoir by using error estimate and control in the multiscale finite-volume framework. SPE Journal, 17(4), 1071–1083.

Hajibeygi, H. and Tchelepi, H.A. [2014] Compositional multiscale ﬁnite-volume formulation. SPE Journal, 19(2), 316–326.

Hou, T.Y. and Wu, X.H. [1997] A multiscale ﬁnite element method for elliptic problems in composite materials and porous media. J. Comput. Phys., 134, 169–189.

Jenny, P., Lee, S.H. and Tchelepi, H.A. [2003] Multi-scale ﬁnite-volume method for elliptic problems in subsur-face ﬂow simulation. J. Comput. Phys., 187, 47–67.

Jenny, P., Lee, S.H. and Tchelepi, H.A. [2006] Adaptive fully implicit scale ﬁnite-volume method for multi-phase ﬂow and transport in heterogeneous porous media. J. Comput. Phys., 217, 627–641.

Kunze, R., Lunati, I. and Lee, S.H. [2013] A multilevel multiscale ﬁnite-volume method. J. Comput. Phys., 225, 502–520.

Lee, S.H., Wolfsteiner, C. and Tchelepi, H.A. [2008] Multiscale finite-volume formulation for multiphase flow in porous media: black oil formulation of compressible, three-phase flow with gravity. Comput. Geosci., 12(3), 351–366.

Lee, S., Zhou, H. and Techelpi, H. [2009] Adaptive multiscale ﬁnite-volume method for nonlinear multiphase transport in heterogeneous formations. J. Comput. Phys., 228(24), 9036–9058.

Lunati, I. and Jenny, P. [2006] Multiscale ﬁnite-volume method for compressible multiphase ﬂow in porous media. J. Comput. Phys., 216(2), 616–636.

Lunati, I. and Jenny, P. [2008] Multiscale ﬁnite-volume method for density-driven ﬂow in porous media. Comput. Geosci., 12(3), 337–350.

Lunati, I., Lee, S. and Tyagi, M. [2011] An iterative multiscale ﬁnite volume algorithm converging to exact solution. J. of Comp. Phys., 230(5), 1849–1864.

Moyner, O. and Lie, K.A. [2014] The multiscale ﬁnite-volume method on stratigraphic grids. SPE Journal, in press, doi:http://dx.doi.org/10.2118/163649-PA.

Stuben, K. [2010] SAMG User’s Manual. Fraunhofer Institute SCAI.

Tomin, P. and Lunati, I. [2013] Hybrid multiscale ﬁnite volume method for two-phase ﬂow in porous media. J. Comput. Phys., 250(13), 293–307.

Trottenberg, U., Oosterlee, C. and Schueller, A. [2001] Multigrid. Elsevier Academic Press.

Wang, Y., Hajibeygi, H. and Tchelepi, H.A. [2014] Algebraic multiscale linear solver for heterogeneous elliptic problems. Journal of Computational Physics, 259, 284–303.

Wolfsteiner, C., Lee, S.H. and Tchelepi, H.A. [2006] Well modeling in the multiscale ﬁnite volume method for subsurface ﬂow simulation. SIAM Multiscale Model. Simul., 5(3), 900–917.

Zhou, H. and Tchelepi, H.A. [2008] Operator based multiscale method for compressible ﬂow. SPE Journal, 13(2), 267–273, doi:http://dx.doi.org/10.2118/106254-PA.

Zhou, H. and Tchelepi, H.A. [2012] Two-stage algebraic multiscale linear solver for highly heterogeneous reser-voir models. SPE J., SPE 141473-PA, 17(2), 523–539.