Time adaptive conservative finite volume method

Journal Pre-proof Time adaptive conservative finite volume method

Patrick Jenny

PII:

S0021-9991(19)30772-7

DOI:

https://doi.org/10.1016/j.jcp.2019.109067

Reference:

YJCPH 109067

To appear in:

Journal of Computational Physics

Received date:

17 May 2019

Revised date:

23 October 2019

Accepted date:

24 October 2019

Please cite this article as: P. Jenny, Time adaptive conservative finite volume method, J. Comput. Phys. (2019), 109067, doi: https://doi.org/10.1016/j.jcp.2019.109067.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Highlights • The research highlight presented in this paper is an adaptive and still conservative time integration scheme for finite volume solvers. In combination with a new flux solver this sub-time stepping scheme is of almost second order (temporally and spatially), robust and much faster (in some of the cases more than 50 times) than conventional time stepping (at almost the same accuracy). This is demonstrated in the paper with numerical studies.

Journal of Computational Physics (2019)

Contents lists available at ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Time Adaptive Conservative Finite Volume Method Patrick Jenny∗ Swiss Federal Institute of Technology, Sonneggstrasse 3, 8092 Zurich, Switzerland

ARTICLE INFO

ABSTRACT

Article history: Received Received in final form Accepted Available online -

Finite volume methods for time dependent problems typically employ time integration schemes with constant time step sizes. The latter, in order to ensure stability and time accurate solutions, have to be chosen small enough, such that the Courant-Friedrichs-Lewy (CFL) number stays below a critical value everywhere. In many cases, however, this time step size limitation leads to huge numbers of very small time steps, which renders simulations very expensive. Motivated by this drawback of conventional time stepping schemes, various sub-time stepping algorithms have been proposed. However, most of them are inherently asynchronous, require small local CFL numbers or are not strictly conservative. In this paper a new adaptive time integration scheme for finite volume methods, which is conservative, of high spatial and temporal order, robust and easy to implement is presented. It relies on local subtime steps which are fractions of a global time step by powers of two, i.e., the grid cells proceed asynchronously in the order of their termination time, but since they all synchronize at the end of each global time step, it is possible to guarantee continuity of the mean fluxes and thus strict conservation at the global time step resolution. Numerical experiments with 1D and 2D test cases demonstrate that the adaptive conservative time integration (ACTI) scheme can achieve extreme speed-up factors over conventional time integration, while still maintaining high spatial and temporal accuracy. c 2019 Elsevier Inc. All rights reserved. 

Communicated by -

1. Introduction Almost all finite volume methods for time dependent problems rely on global time stepping schemes, were all cell values are updated synchronously. In order to guarantee stability with explicit time integration, but also to ensure an acceptable level of accuracy if an implicit scheme is employed, the global time step has to be small enough to honor the most restrictive Courant-Friedrichs-Lewy (CFL) criterion [1] among all grid cells. Consequently, a large fraction

∗ Corresponding

author: Tel.: +41-44-632-6987; e-mail: [email protected] (Patrick Jenny)

2

Patrick Jenny / Journal of Computational Physics (2019)

of grid cells may employ time steps much smaller than the one allowed based on the local CFL criterion. Motivated by the associated computational overhead, various local time stepping schemes have been developed. Such methods aim at exploiting the maximum possible local time step size and thus try to reduce the total number of steps and associated computational cost. However, such schemes necessarily are to some extent asynchronous, which introduces other challenges regarding spatial and temporal accuracy, conservation and robustness. Osher and Sanders [2] were among the first who proposed a locally varying time step size for 1D, non-linear scalar conservation laws. Their method, which is of first spatial order, is based local CFL number restrictions and relies on a conservative predictor-corrector scheme and on an arbitrary number of domains with fine and coarse temporal increments. To achieve conservation, intermediate values have to be stored at domain interfaces. The authors proved convergence for scalar monotone discretizations. Berger and Oliger [3] devised a sub-time stepping method for hyperbolic systems of conservation laws on discontinuous grids. Their approach can be seen as an extension of Osher and Sanders’ method [2]. They started their derivation from a discrete integral form of the conservation laws, which led to a conservative interface equation. As a result, they do not require a corrector step as Osher and Sanders’s method [2]. Later, Dawson and Kirby [4] introduced a similar method as the one by Osher and Sanders [2], which is of higher order, however. In addition, they rigorously proved a maximum principle. Kleb and Batina [5] introduced sub-time steps, which are fractions of the global time step by powers of two, such that the local CFL criterion is satisfied everywhere. As a result, the cells get synchronized at the end of each global time step and conservation as well as high temporal accuracy can be achieved. Cells with larger time steps get executed first, while those with smaller time steps have to catch up. This approach requires interpolation in time and storage of older values. Zhang et al. [6] proposed a sub-time stepping scheme with re-evaluation of the time step size after each global step in order to capture high wave speeds and Nutaro and Zeigler [7] devised a discrete event simulation algorithm for solving simple conservation laws. Constantinescu and Sandu [8] and Fox [9] devised multirate Runge-Kutta methods that preserve the stability properties of conventional time stepping. More recently, Seny et al. [10] presented a parallel implementation of the method by Constantinescu and Sandu. Shao et al.[11] devised an asynchronous parallel discrete event simulation algorithm based on local time stepping criteria, which is adapted for the hybrid finite-element nodecentred finite volume method, and they demonstrated impressive speedup for flow in fractured porous media. Another totally asynchronous local time stepping scheme was presented by Cavalcanti et al. [12]. Their method is based on a piecewise linear reconstruction in space and time in order to achieve high accuracy for variable time steps. They show results for scalar transport on unstructured grids and achieve speedup factors of about one order of magnitude. Carciopolo et al. [13] devised an implicit multirate method for hyperbolic problems, which is based on similar ideas as that by Kleb and Batina [5]. The local time steps are fractions of the global time step by powers of two and the cells with large time steps are processed first. To honor strict conservation, it is ensured that the fluxes used for the updates in the neighboring cells with smaller time steps are interpolated such that their sums are consistent. Since their scheme is implicit, it allows to employ very large CFL numbers locally, where temporal accuracy is not required. In [14], Carciopolo et al. employed their multirate method to multi-phase flow in porous media. In this paper, an adaptive conservative time integration (ACTI) scheme is devised, which is based on local time step sizes being smaller than the global time step by powers of two. Like in the scheme by Kleb and Batina [5] the cells get synchronized after each global time step, where conservation is exactly honored. But opposed to their scheme, here the cells with smaller time steps get processed first and their accumulated fluxes are then consistently used for the update of adjacent cells with larger time steps. This makes interpolations, interface equations and correction steps unnecessary and allows to employ established high order flux discretizations. The new ACTI scheme is strictly conservative and the local time steps can be dynamically adapted for each individual cell. It was combined with a flux solver, which is of second order in space and time. Numerical tests with linear and non-linear 1D and 2D test cases proved that the proposed ACTI scheme is very robust for large local CFL numbers, that it is strictly conservative and that it is of second order in time and space. In all cases the ACTI results are almost indistinguishable from those with conventional time stepping with the same maximum CFL number (maximum CFL numbers of 1.0 and 0.5 were employed in 1D and 2D, respectively), but with ACTI the required computational work could be reduced by factors beyond 50. The only additional overhead consists in sorting the cells with respect to their end times and time step sizes. It is very easy to implement the presented ACTI scheme and it should be possible to adopt it into existing finite volume codes without too much difficulty. The paper is structured as follows. In section 2 the ACTI scheme is introduced and discussed. Section 3 deals with spatially and temporally accurate numerical flux discretization. Although not the main focus of the paper, this is relevant, since for time dependent problems as considered here not only spatial, but also temporal accuracy is

Patrick Jenny / Journal of Computational Physics (2019)

3

important. An assessment of the new ACTI scheme with high order flux discretization is provided in section 4. Both spatial and temporal accuracy are investigated and compared to conventional time stepping. Most importantly, it is shown that compared to conventional time stepping dramatic speedup is possible with ACTI without compromizing neither accuracy nor robustness. Final discussions and conclusions are provided in section 5. 2. Adaptive Conservative Time Integration (ACTI) Here a new adaptive and conservative time integration scheme is presented, which is based on local time steps being fractions of a global time step by powers of two; similar as in the method by [5]. This choice of time step sizes allows to advance in such a way that all cell states get synchronized at the end of each global time step. Consequently, it is possible to use consistent numerical fluxes and thus obtain a conservative method. To explain this method, we consider scalar conservation laws of the form ∂u + ∇· f ∂t

= 0,

(1)

where u(x, t) is the spatially and temporally dependent scalar field and f (x, t) the flux. Note that bold symbols refer to vectors and subscript indexes to individual components. Further, throughout the paper all quantities are considered dimensionless. To compute numerical solutions, the finite volume method with explicit integration in time is considered, i.e., the computational domain Ω gets divided into N temporally fix disjunct cells ΩI (I ∈ {1, . . . , N}), which store the numerical cell values  1 U In ≈ u(x, tn ) dV(x) (2) |ΩI | ΩI at discrete times tn . To update these average numerical u-values, one can employ the scheme 1  n→n+1 F , U In+1 = U In + |ΩI | J∈A(I) J→I where A(I) is the index-set referring to all cells adjacent to ΩI and  tn+1  F n→n+1 ≈ f (x, t) · nJ→I dS (x) dt J→I tn

Γ JI

(3)

(4)

is the numerical flow across the interface Γ JI from Ω J to ΩI during the time t ∈ [tn , tn+1 ]. Further, nJ→I is the unit normal at Γ JI pointing from Ω J to ΩI . For simplicity, but without loss of generality, we now consider a hyperbolic problem with the approximation =

F n→n+1 J→I where

n f J→I

U In+1

n f J→I |Γ JI |Δt,

(5)

is the numerical flux and Δt the global time step size. Thus one obtains the forward Euler method =

U In+1 +

Δt  n f |Γ JI |, |ΩI | J∈A(I) J→I

(6)

to compute the numerical u-values at the new time tn+1 from those at the old time tn . In this case, in order to ensure stability, one has to choose Δt ≤

min (ΔtI )

I∈{0,...,N}

with the maximum allowed cell time step sizes   |ΩI | ΔtI = CFL min , J∈A(I) |Γ JI ||a J→I |

(7)

(8)

where a J→I is the wave speed at the interface ΓI J and CFL is case and scheme dependent. Note that while stability can be achieved with large CFL values, if an implicit time integration scheme is employed, they are still restricted in order to keep temporal truncation errors in reasonable bounds. In any case, the number of time steps required to complete

Patrick Jenny / Journal of Computational Physics (2019)

4

a simulation is determined by the cell with the smallest ΔtI , i.e., in many cells the maximum possible cell time step size cannot be exploited, which thus may be much larger than the global time step. Proceeding asynchronously, on the other hand, bears the potential to exploit the maximum allowed cell time steps everywhere thus resulting in an overall reduction of flow computations and variable updates. However, most such methods are not conservative and it is a challenge to maintain stability, robustness and high spatial and temporal accuracy. Here we present a sub-time stepping scheme, which is conservative, stable and spatially and temporally accurate. It relies on cell based sub-time steps, which are fractions of a global time step by powers of two. This choice allows to employ the accumulated flows from updates of adjacent cells with smaller local time steps, which guarantees exact balance and synchronization at the end of each coarse time step. A pseudo code of this adaptive conservative time integration (ACTI) method is presented in Algorithm. 1. First, the global time t is set to zero and a grid with N cells is Algorithm 1 : Pseudo code of the new sub-time stepping scheme · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · ·

define grid with N cells ΩI ; I ∈ {1, . . . , N} t←0 ∀I ∈ {1, . . . , N} { U I ← initialize tI ← 0 ∀J ∈ A(I) : F J→I ← 0 } while (t < tend ) { ∀I ∈ {1, . . . , N}: compute ΔtI  Δt ← min maxI∈{1,...,N} (ΔtI ), tend − t t ← t + Δt ∀I ∈ {1, . . . , N} { LI ← ceil(ln(ΔtI /Δt)/ ln(1/2)) ΔtI ← 2−LI Δt } while (minI (tI ) < t) { tnext,min = minI∈{1,...,N} (tI + ΔtI ) S ← (I ∈ {1, . . . , N}|tI + ΔtI = tnext,min ∧ ∀i < j ≤ size(S ) ⇒ LS i ≥ LS j ) for i from 1 to size(S ) { I ← Si tI ← tI + ΔtI ∀J ∈ A(I) { if (t J < tI ) { F J→I ← inflow from Ω J into ΩI during time t ∈ [tI − ΔtI , tI ] F I→J ← F I→J − F J→I } }  U I ← U I + |ΩI |−1 J∈A(I) F J→I ∀J ∈ A(I) : F J→I ← 0 } } }

defined. For each cell ΩI the dependent variable U I is initialized and the cell time tI and all flows across its interfaces are set to zero. Next, the global time step loop is entered, i.e., multiple global time steps are performed, until the global time t reaches the end time tend . At the beginning of each time step, the maximum allowed time step size ΔtI is determined for each cell ΩI and the global time step size Δt is chosen, whereas it has to be ensured that Δt does not exceed the remaining time till the end time. Then, after incrementing the global time with the global time step size, the levels LI are computed for each cell ΩI . The cell level is defined such that the maximum allowed cell time step is larger than or equal to 2−LI Δt and smaller than 2−(LI −1) Δt. The cell time step size ΔtI is then set to 2−LI Δt. Next, the sup-time step loop is entered and the following operations are performed until all cell times tI are equal to the global time t. At the beginning of each sub-time step iteration, a list S of cell indexes I, for which the end time tI + ΔtI is equal to the minimum end time of all cells, and which is sorted according to their levels LI (higher levels first), is created. The whole list is then traversed from the beginning till the end, while in each iteration first the cell time of the corresponding cell with index I is incremented by the cell time step size. Then, from each adjacent cell J the flow into ΩI is computed, if the cell time t J is smaller than tI . In that case the computed flow is subtracted from the

Patrick Jenny / Journal of Computational Physics (2019)

5

accumulation of flows from cell I into cell J. Otherwise, if t J is larger or equal to tI , there is no need to compute the flow, since this has already been done when cell J was updated earlier. In any case, once all flows across the faces of ΩI are available, U I gets updated and the flows in and out of ΩI are set to zero. For all studies presented in this paper the boundary conditions were treated the same way with ACTI as with conventional time stepping. An illustrtion of a concrete example is shown in Fig. 1. A domain with six cells is considered and the sub-figures show how the sub-steps in which the cells get updated until all reach the end of the global time step of size Δt. The sub-time step sizes for Ω1−6 are Δt, Δt, Δt/2, Δt/4, Δt/4 and Δt, respectively. During sub-step 1, the cells with earliest end times (end times are indicated by horizontal lines at the top of each cell bin) and the smallest sub-time step size get processed first. This means, the required flows across cell interfaces are computed (indicated by the vertical thick black lines along the cell bin boundaries) and subsequently the values in cells Ω4 and Ω5 are updated (indicated by the thick arrows pointing upwards). The vertical gray bars indicate up to which time the corresponding cells are already updated. The same two cells are updated again in sub-step 2, but in sub-step 3, Ω3 has the earliest end time. Note that the flow between Ω3 and Ω4 was already calculated, i.e., it is the accumulation of the flows at that interface from the first two sub-steps, and thus only the flow between Ω2 and Ω3 has to be computed in order to update cell Ω3 during sub-step 3. In sub-steps 4 and 5, once more Ω3 and Ω4 have the earliest end times and they are updated after computing the required flows across the cell interfaces. Now all cells have the same end times and the remaining cell processing order solely depends on the sub-time step size , i.e., cells with smaller sub-time step sizes are processed first. Concretely these are Ω4 and Ω5 in sub-step 5, Ω3 in sub-step 6 and Ω1 and Ω2 in sub-step 7. Note that flows only have to be computed at cell interfaces and times where indicated by black vertical lines, while gray vertical lines indicate where accumulated flows from previous sub-steps are used. Finally, all cells process times reach the end of the global time step. It is easy to figure out that in this concrete example 17 flow computations and 13 cell updates are necessary for ACTI to complete a global time step; versus 28 flow computations and 24 cell updates with conventional time stepping (using a time step size of Δt/4). This algorithm guarantees that all cell times are the same at the end of each global time step, i.e., they get synchronized, but in general they advance asynchronously in between. Moreover, since all computed flows are accounted for twice; once after being calculated to update the variable in the corresponding cell, and once as part of the accumulated flows used to update the neighboring cells. The computational gain compared to conventional time stepping is that most cells can proceed with time steps larger than the minimum cell time step size, which leads to fewer flow computations and fewer updates. Additional work is required to sort those cells with minimum end times according to their levels, but there exist very fast and established sorting algorithms of O(N  log(N  )) (note that N  refers to the number of cells with the same minimum end time and typically N   N). More on the cost for sorting is presented in section 4.3. 3. Numerical Flux Discretization Although not the main focus of this paper, four numerical flux schemes are briefly discussed. We consider hyperbolic scalar equations of the form ∂ f1 (u) ∂ f2 (u) ∂u + + ∂t ∂x1 ∂x2

=

0,

(9)

which are solved on Cartesian grids with cells ΩI1 ,I2 (I1 ∈ {1, . . . , N1 } and I2 ∈ {1, . . . , N2 }) and grid spacings h1I1 ,I2 and h2I1 ,I2 in x1 - and x2 -direction, respectively. Note that for convenience, in contrast to the previous section, two indexes are used here to identify grid cells. Further, |ΩI1 ,I2 | = h1I1 ,I2 h2I1 ,I2 . With conventional time stepping the numerical cell values are updated as ⎛ ⎞ ) − f1 (U In→n+1 ) ) − f2 (U In→n+1 ) ⎟⎟ f2 (U In→n+1 ⎜⎜⎜ f1 (U In→n+1 1 −1/2,I2 1 +1/2,I2 1 ,I2 −1/2 1 ,I2 +1/2 ⎟ n+1 n ⎜ ⎟⎟⎠ . U I1 ,I2 = U I1 ,I2 + Δt ⎜⎝ + (10) h1I1 ,I2 h2I1 ,I2 The numerical flux components are functions of u-estimates at the corresponding cell faces, e.g. the first flux component at the face ΓI1 +1/2,I2 is f1 (U In→n+1 ) (fluxes and u-values at the other cell faces are calculated analogously). The 1 +1/2,I2 scheme used here is based on a reconstruction of u at tn + Δ/2 around the interface. It employs numerical u-values at tn to obtain an initial interpolation, which is then advected with the estimated wave speed at the interface for half a time step. As a result one obtains an estimate of the local average u-distribution, but of interest for the numerical flux

Patrick Jenny / Journal of Computational Physics (2019)

6t

6t

6t

6t

6

t

t

sub-step 1

sub-step 2

11 12 13 14 15 16

11 12 13 14 15 16

t

t

sub-step 3

sub-step 4

11 12 13 14 15 16

11 12 13 14 15 16

t

t

sub-step 5

sub-step 6

11 12 13 14 15 16

11 12 13 14 15 16

t

t

sub-step 7

end of global time step

11 12 13 14 15 16

11 12 13 14 15 16

Fig. 1. Illustration of ACTI with a domain consisting of six cells, which get updated during seven sub-steps, until all cell times reach the end of the global time step of size Δt. The respective sub-time step sizes for Ω1−6 are Δt, Δt, Δt/2, Δt/4, Δt/4 and Δt. Vertical black and gray lines along the cell interfaces indicate where and when the flows have to be calculated and where accumulated flows from previous sub-steps are used, respectively. Gray bars indicate up to which time the individual cells are already updated and the thick upward pointing arrows indicate which cells are updated during the individual sub-steps.

Patrick Jenny / Journal of Computational Physics (2019)

7

calculation is only the value at the interface center at time tn+1/2 . For the time step from tn to tn+1 the computation at the vertical interface ΓI1 +1/2,I2 can be expressed as ⎧ σ1 σ1 σ2 ⎪ +β1 I21 ,I2 h1I1 ,I2 −Δt(β2 I21 ,I2 a1I1 +1/2,I2 + β3 I21 ,I2 a2I1 ,I2 ) if a1I1 +1/2,I2 > 0 U In1 ,I2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ σ1 σ1 σ2 ⎨ =⎪ U In→n+1 (11) U In1 +1,I2 −β1 I12+1,I2 h1I1 +1,I2 −Δt(β2 I12+1,I2 a1I1 +1/2,I2 + β3 I12+1,I2 a2I1 +1,I2 ) if a1I1 +1/2,I2 < 0 ⎪ 1 +1/2,I2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ n n ⎪ ⎪ ⎩ U I1 ,I2 +U I1 +1,I2 + β σ1I1 ,I2 h1I1 ,I2 −σ1I1 +1,I2 h1I1 +1,I2 − Δt β σ2I1 ,I2 a2I1 ,I2 +σ2I1 +1,I2 a2I1 +1,I2 else, 1 3 2 4 4 where a1 and a2 are the wave speed components in x1 - and x2 -directions, respectively (central discretization is used for the face normal component and simple upwind for the transverse component). Note that the upper and lower branches in Eq. (11) describe the scheme for positive and negative wave speeds in x1 direction, respectively. The computation of the numerical u-estimations at horizontal cell interfaces, e.g. between ΩI1 ,I2 and ΩI1 ,I2 +1 , is analogous to Eq. (11) (by swapping the two indexes) and is not explicitly shown. The parameters β1 , β2 and β3 were introduced to distinguish between different interpolation schemes. The one favored here, which is explained next for the face ΓI1 +1/2,I2 and a1I1 +1/2,I2 > 0 (upper branch in Eq. (11)), uses β1 = β2 = β3 = 1. In this case the reconstruction in ΩI1 ,I2 at time tn reads U˜ I1 ,I2 (x, tn )

= U In1 ,I2 + σ1I1 ,I2 (x1 − x1I1 ,I2 ) + σ2I1 ,I2 (x2 − x2I1 ,I2 )

(12)

and at t ∈ [tn , tn+1 ] one obtains U˜ I1 ,I2 (x, t)

= U In1 ,I2 + σ1I1 ,I2 (x1 − x1I1 ,I2 − (t − tn )a1I1 +1/2,I2 ) + σ2I1 ,I2 (x2 − x2I1 ,I2 − (t − tn )a2I1 ,I2 ).

(13)

One can easily see that for t = tn + Δt/2, x1 = x1I1 ,I2 + h1I1 ,I2 /2 and x2 = x2I1 ,I2 the last expression and the upper branch of Eq. (11) are identical. The slope component normal to the interface, in the upper branch of Eq. (11) this is σ1I1 ,I2 , is computed by the Koren limiter [15] as       s2 + 2 s1 2 s1 ,2 , ,0 , (14) σ1I1 ,I2 = σ (s1 , s2 ) = s2 max min min 3 s2 s2 where s1 = 2

U In1 +1,I2 − U In1 ,I2 h1I1 +1,I2 + h1I1 ,I2

and s2 = 2

U In1 ,I2 − U In1 −1,I2 h1I1 ,I2 + h1I1 −1,I2

.

(15)

On the other hand, for the slope in the cross direction (σ2I1 ,I2 in the upper branch of Eq. (11)) the very simple upwind based estimate ⎧ U In ,I −U In ,I −1 ⎪ ⎪ if a2I1 ,I2 > 0 2 h2 1 2 +h2 1 2 ⎪ ⎪ ⎪ I1 ,I2 I1 ,I2 −1 ⎪ ⎨ σ2I1 ,I2 = ⎪ (16) ⎪ ⎪ ⎪ U In ,I +1 −U In ,I ⎪ 1 2 ⎪ ⎩ 2 h 1 2 +h else 2 2 I1 ,I2 +1

I1 ,I2

is used. All other slopes are obtained analogously. As already mentioned, the parameters β1−3 determine the flux scheme, i.e., β1 = β2 = β3 = 0 corresponds to simple upwinding (flux scheme 1), β1 = 1 and β2 = β3 = 0 to a MUSCL upwind scheme (flux scheme 2) [16], β1 = β2 = 1 and β3 = 0 to MUSCL with slope advection (flux scheme 3) [17] and β1 = β2 = β3 = 1 to MUSCL with advection of inclined reconstruction (flux scheme 4). An illustration of scheme 4 is shown in Fig. 2. The blue surface represents the reconstruction of u in cell ΩI1 ,I2 (horizontal gray plane) at the beginning of the time step. This manifold is then translated by the wave speed aI1 +1/2,I2 obtained at the interface ΓI1 +1/2,I2 times the time step size Δt (green surface). Finally, the u value used for the flux calculation is the mean of U˜ I1 ,I2 (xI1 +1/2,I2 , tn ) and U˜ I1 ,I2 (xI1 +1/2,I2 , tn+1 ). To demonstrate the importance of not only spatial, but also temporal accuracy, advection with f1 = f2 = u and the initial condition  1 if x ∈ [0.2, 0.8] × [0.2, 0.8] u(x, t = 0) = (17) 0 else

Patrick Jenny / Journal of Computational Physics (2019)

8

u

~

n

U I 1, I 2 (x I 1+1/2, I 2,t )

aI

~

n+1

U I , I (x I +1/2, I ,t 1

+1/2, I2 6t 1

2

1

)

2

x2 x1 1I

1

1I

, I2

1

+1, I 2

Fig. 2. Illustration of flux scheme 4: Reconstruction of u in cell ΩI1 ,I2 (horizontal gray plane) at the beginning of the time step (blue surface) and reconstruction manifold translated by the wave speed aI1 +1/2,I2 times the time step size Δt (green surface). The u value used for the flux calculation is the mean of U˜ In +1/2,I and U˜ In+1 +1/2,I at the interface ΓI1 +1/2,I2 (vertical gray plane). 1

2

1

2

is simulated on a square domain Ω = [0, 1] × [0, 1] with periodic boundary conditions. Again, all dependent and independent variables are dimensionless. A Cartesian equidistant 100 × 100 grid is considered and results at t = 1 (when the exact solution matches the initial condition) obtained with the four flux schemes, all employing CFL=0.25, are compared. Figure 3 shows the grid, initial condition and numerical solutions with the four schemes. As expected, the simple upwind scheme (flux scheme 1) is very dissipative and inaccurate (Fig. 3c). Flux schemes 2-4 become identical for Δt → 0, but for larger CFL numbers (here CFL=0.25), which is of interest in order to save computational time, remarkable differences can be observed. While flux schemes 2 and 3 (Figs. 3d-e) show unacceptable distortions, flux scheme 4 (Fig. 3f) well preserves the square contours of the initial condition. Furthermore, flux scheme 4 proved to be much more stable than schemes 2 and 3, i.e., larger CFL numbers can be used. Therefore, all remaining numerical studies presented in this paper are based on flux scheme 4.

4. Assessment of the ACTI Scheme For demonstration purpose we first consider 1D advection and the 1D inviscid Burger’s equation with non-uniform grids. Then we study 2D advection on a non-uniform grid and finally tracer transport in a porous medium with a highly conductive fracture. In all these test cases either wave speed or cell size significantly vary throughout the computational domain, which makes them interesting targets for the new ACTI scheme. To assess the performance of ACTI compared to conventional time stepping, accuracy and computational efficiency are empirically investigated and compared to classical time stepping. Note that also here all dependent and independent variables are dimensionless. 4.1. 1D Advection For our first investigations of the presented ACTI scheme we consider the scalar advection equation ∂u ∂u + a ∂t ∂x

=

0

(18)

with a = 1 on the domain Ω = [0, 1] with periodic boundary conditions. The initial condition is u(x, t = 0) = sin(2πx) and the grid cells are defined as ΩI = [xI − hI /2, xI + hI /2] (I ∈ {1, . . . , N}), where hI = 1/N + cos(2πI/N)/(1.001N) (I ∈ {1, . . . , N}) is the cell size. Figure 4 depicts the grid spacing for N = 100; it can be observed that it is much smaller near the domain center than near the boundaries. Note that for CFL=1 the allowed time step size of a cell ΩI

Patrick Jenny / Journal of Computational Physics (2019)

9

1

0.8

0.8

0.8

0.6

0.6

0.6

x2

1

x2

1

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0.2

0.4

x1

0.6

0.8

0

1

0

0.2

(a)

0.4

x1

0.6

0.8

1

0

(b)

1

1

1

1

0.8

0.8

0.8

0.6

0.6

0.6

0.6

x2

x2

0.8

0.4

0.4

0.4

0.4

0.2

0.2

0.2

0.2

0

0

0.2

0.4

x1

0.6

0.8

1

0

0

0

0.2

(c)

0.4

x1

0.6

0.8

1

0

(d)

1

1

1

1

0.8

0.8

0.6

0.6

0.6

0.6

x2

0.8

x2

0.8

0.4

0.4

0.4

0.4

0.2

0.2

0.2

0.2

0

0

0.2

0.4

x1

(e)

0.6

0.8

1

0

0

0

0.2

0.4

x1

0.6

0.8

1

0

(f)

Fig. 3. Comparison of flux schemes 1-4 for 2D advection with constant wave speed a = (1, 1)T and periodic boundary conditions: (a) Cartesian equidistant 100 × 100 grid; (b) initial condition; (c) solution at t = 1 with flux scheme 1; (d) solution at t = 1 with flux scheme 2; (e) solution at t = 1 with flux scheme 3; (f) solution at t = 1 with flux scheme 4. For all flux schemes a CFL number of 0.25 was employed.

Patrick Jenny / Journal of Computational Physics (2019)

10

0.02 0.018 0.016 0.014

hI

0.012 0.01

0.008 0.006 0.004 0.002 0

0

10

20

30

40

50 I

60

70

80

90

100

Fig. 4. Grid spacing hI used for the 1D studies (linear advection and inviscid Burger’s equation).

is ΔtI = hI /a, i.e., if classical time stepping is applied, the time step size is dictated by the smallest cells. Figure 5 shows numerical solutions with ACTI (solid lines) and conventional time stepping (dashed lines) at t = 0.25 (left) and at t = 1 (right). The initial condition is represented by symbols. Both numerical solutions are almost identical

1

initial cond. ACTI conv. time stepping

0.5

0.5

0

0

u

u

1

-0.5

exact ACTI conv. time stepping

-0.5

-1

-1 0

0.2

0.4

x

t = 0.25

0.6

0.8

1

0

0.2

0.4

x

0.6

0.8

1

t=1

Fig. 5. 1D advection with a = 1 and periodic boundary conditions: Symbols represent the initial condition (for t = 1 the exact solution) and lines numerical solutions with ACTI (solid) and conventional time stepping (dashed); left at t = 0.25 and right at t = 1.

and both are in extremely good agreement with the analytical solution at t = 1. The interesting difference, however, is that conventional time stepping (with a maximum CFL number of one) requires 34 times more cell updates and 25 more flux calculations than ACTI. These ratios are also reported in table 2 together with those of the other test cases discussed below. The levels LI and CFL numbers ΔtI a/hI are shown in Figs. 6a and 6b, respectively. It can be seen that the highest level is 11, i.e., the smallest time step is 211 times smaller than the global one. Due to sub-time step adaptation, however, all cell based CFL numbers are in the range between 0.5 and 1.0, if ACTI is employed. Regarding  conservation it was verified that at the end of each global time step NI=1 (hI U I ) is the same as at the beginning of the simulation. A quantitative accuracy assessment of ACTI and conventional time stepping is presented in Fig. 7, which shows double-logarithmic convergence plots, i.e., the error 2-norm at t = 1 as a function of N, whereas the maximum CFL number was 1.0 in all cases. The slope for conventional time stepping is almost exactly −2 (compare with solid line), which indicates second order in space and time. For ACTI the slope is slightly worse, and at this point it can only be speculated about the drop in convergence order for larger N. It is encouraging, however, that ACTI is already very accurate for relatively small N and that its solutions are very close to those obtained with conventional time

Patrick Jenny / Journal of Computational Physics (2019)

11

12 1 10 0.8 0.6

CFLI

level LI

8 6

0.4

4

0.2

2 0

0

10

20

30

40

50 I

60

70

80

90

0

100

0

10

20

30

(a)

40

50 I

60

70

80

90

100

(b)

Fig. 6. 1D advection with a = 1 and periodic boundary conditions: (a) Levels LI and (b) CFL numbers ΔtI a/hI .

error 2-norm

0.1

-2

N ACTI conv. time stepping

0.01

0.001

25

N

50

100

Fig. 7. 1D advection with a = 1 and periodic boundary conditions: Error 2-norms at t = 1 as functions of N; the maximum CFL number was 1.0 for both ACTI and conventional time stepping.

stepping. 4.2. 1D Inviscid Burger’s Equation Similar conclusions as for the linear case can be drawn from numerical experiments with the inviscid Burger’s equation ∂u ∂u2 /2 + = 0. (19) ∂t ∂x The same domain, grid, global time step size, initial- and boundary conditions are considered. The wave speeds at the cell interfaces, e.g. between cells ΩI and ΩI+1 , is estimated as aI+1/2 = (U I + U I+1 )/2. Figure 8 shows the initial condition (cross symbols) together with numerical solutions for N = 100 at t ∈ {0.0625, 0.125, 0.1875, 0.25} computed with ACTI (solid lines) and conventional time stepping (dashed lines). The circle symbols depict the exact solutions u(xI , tn ) for tn ∈ {0.0625, 0.125}, which was obtained by finding the root of R(u, xI , tn )

= u − sin(2π(xI − u tn ))

for every cell mid-point, which can be achieved via the Newton scheme  −1 ∂R(u, xI , tn )  uν+1 = uν − R(uν , xI , tn )  uν ∂u

(20)

(21)

Patrick Jenny / Journal of Computational Physics (2019)

12

1

0.5

init. cond. ACTI conv. time stepping exact

0.5

0

u

u

1

init. cond. ACTI conv. time stepping exact

-0.5

0

-0.5

-1

-1 0

0.2

0.4

0.6

x

0.8

1

0

0.2

t = 0.0625

0.6

x

0.8

1

t = 0.125 1

init. cond. ACTI conv. time stepping

0.5

0.5

0

0

u

u

1

0.4

-0.5

init. cond. ACTI conv. time stepping

-0.5

-1

-1 0

0.2

0.4

x

t = 0.1875

0.6

0.8

1

0

0.2

0.4

x

0.6

0.8

1

t = 0.25

Fig. 8. 1D inviscid Burger’s equation with periodic boundary conditions: Cross symbols represent the initial condition and lines numerical solutions with ACTI (solid) and conventional time stepping (dashed) at four different times. The circle symbols depict the exact solutions for t ∈ {0.0625, 0.125}

Patrick Jenny / Journal of Computational Physics (2019)

error 2-norm

0.1

13

-2

N ACTI conv. time stepping

0.01

0.001

25

N

50

100

Fig. 9. 1D inviscid Burger’s equation with periodic boundary conditions: Error 2-norms at t = 0.125 as functions of N; the maximum CFL number was 1.0 for both ACTI and conventional time stepping.

10

10

8

8 level LI

12

level LI

12

6

6

4

4

2

2

0

0

10

20

30

40

50 I

60

t = 0.125

70

80

90

100

0

0

10

20

30

40

50 I

60

70

80

90

100

t = 0.25

Fig. 10. 1D inviscid Burger’s equation with periodic boundary conditions: Levels LI at two different times.

with a reasonable initial guess, e.g. u0 = sin(2πxI ). Once converged, uν+1 = uν is the analytical solution at xI and tn . Again, the numerical solutions obtained with the two methods are very accurate and almost indistinguishable from each other. They lead to a very sharp shock with a strength of two at t = 0.25. Grid convergence for t = 0.125 is shown in Fig. 9; again, almost second order is observed for both methods, while the accuracy of ACTI is slightly lower than with conventional time stepping. Regarding efficiency, it is impressive that for a simulation time of tend = 0.25 ACTI requires 68 times less cell updates and 50 times less flux calculations than conventional time stepping. Table 2 contains these ratios together with those for all other test cases presented in this paper. Figure 10 depicts the levels LI at t = 0.125 and t = 0.25 and the respective CFL numbers ΔtI U I /hI are shown in Fig. 11. Like in the previous 1D advection case, the smallest time step is 211 times smaller than the global one and most CFL numbers are in the range between 0.5 and 1.0. 4.3. 2D Advection To test the presented sub-time stepping scheme for multi-dimensional problems, first the 2D advection equation ∂u ∂u ∂u + a1 + a2 ∂t ∂x1 ∂x2

=

0

(22)

with a1 = a2 = 1, a domain Ω = [0, 1] × [0, 1] with periodic boundary conditions and a Cartesian grid with N1 × N2 cells are considered. The grid spacings in x1 - and x2 -directions are h1I1 ,I2 and h2I1 ,I2 , respectively, and the cells ΩI1 ,I2 =

Patrick Jenny / Journal of Computational Physics (2019)

14

0.8

0.8

0.6

0.6

CFLI

1

CFLI

1

0.4

0.4

0.2

0.2

0

0

10

20

30

40

50 I

60

70

80

90

100

0

0

10

20

30

t = 0.125

40

50 I

60

70

80

90

100

t = 0.25

Fig. 11. 1D inviscid Burger’s equation with periodic boundary conditions: CFL numbers ΔtI a/hI at two different times.

scheme ACTI conventional time stepping

CPU times [s]:

total 18.7 115.2

sorting 10.8

Table 1. Total CPU times for the 2D advection test case (N1 = N2 = 100) with and without ACTI. For ACTI also the CPU time for sorting is shown.

[x1I1 ,I2 − h1I1 ,I2 /2, x1I1 ,I2 + h1I1 ,I2 /2] × [x2I1 ,I2 − h2I1 ,I2 /2, x2I1 ,I2 + h2I1 ,I2 /2] (I1 ∈ {1, . . . , N1 } and I1 ∈ {1, . . . , N2 }) have a size of h1I1 ,I2 h2I1 ,I2 . The grid spacing in x1 -direction is the same as in the 1D test cases, i.e., h1I1 ,I2 = 1/N1 + cos(2πI1 /N1 )/(1.001N1 ) (N1 ∈ {1, . . . , N1 }); in x2 -direction, on the other hand, the spacing is uniform. Figure 12a showns the grid for N1 = N2 = 100. The initial condition u(x, t = 0) = sin(2πx1I1 ,I2 ) cos(2πx2I1 ,I2 ) shown in Fig. 12b is at the same time the exact reference for the numerical solutions at t = 1 shown in Figs. 12d and 12f obtained with ACTI and conventional time stepping, respectively. The corresponding errors are depicted in the maps of Figs. 12c and 12e. It can be observed that the errors produced by the two methods are comparable and in both cases the maximum error is less than 2.5% of the solution amplitude. Double logarithmic convergence plots showing the error 2-norms at t = 1 as functions of N1 = N2 are presented in Fig. 13; dashed line for ACTI and doted line for conventional time stepping. The maximum CFL number was 0.5 in both cases for all grids. Like for the 1D advection case the slope for conventional time stepping is almost exactly −2 (compare with solid line), which indicates second order in space and time. Again, the convergence rate is slightly worse for ACTI, which is expected, since in most cells much larger time steps are applied. Nevertheless, it is an encouraging finding that the corresponding solutions are in extremely close agreement; even for relatively coarse spatial resolutions and quite large CFL numbers. Regarding efficiency, it can be reported that conventional time stepping (with a maximum CFL number of 0.5) requires 35 times more variable updates and 30 times more flux calculations than ACTI. These ratios are also found in table 2 together with those for the other test cases. Although no effort has been made to optimize the implementation of ACTI, profiling of the code for this test case with N1 = N2 = 100 was performed and the total CPU times and that for sorting are shown in table 1; note, however, that more relevant at this stage is the potential efficiency gain, which is reflected by the reduction of flux computations and cell updates. It can be seen that with the current implementation of ACTI approximately 58% of the computational time is spent for sorting. Currently the sorting algorithm ”sort” from the C++ standard library is used, which is an implementation of QuickSort [18] and thus scales with O(N log(N)), but the way it is used in the code may be improved in the future. Nevertheless, for this case a speed up factor of almost six was achieved with ACTI in comparison with conventional time stepping.

Patrick Jenny / Journal of Computational Physics (2019)

1

0.8

0.8

0.6

0.6

1

0.5

0

x2

x2

1

15

0.4

0.4

0.2

0.2

0

0

0.2

0.4

x1

0.6

0.8

0

1

-0.5

0

0.2

(a)

0.4

x1

0.6

0.8

1

-1

(b)

1

0.025

1

1

0.02 0.8

0.015

0.8

0.5

0.01 0.005

x2

0

0.4

-0.005

0.6 0

x2

0.6

0.4

-0.01 0.2

-0.015

-0.5

0.2

-0.02 0

0

0.2

0.4

x1

0.6

0.8

1

-0.025

0

0

0.2

(c)

0.4

x1

0.6

0.8

1

-1

(d)

1

0.025

1

1

0.02 0.8

0.015

0.8

0.5

0.01 0.005

x2

0

0.4

-0.005

0.6 0

x2

0.6

0.4

-0.01 0.2

-0.015

-0.5

0.2

-0.02 0

0

0.2

0.4

x1

(e)

0.6

0.8

1

-0.025

0

0

0.2

0.4

x1

0.6

0.8

1

-1

(f)

Fig. 12. 2D advection with a = (1, 1)T and periodic boundary conditions: (a) Computational 100 × 100 grid, (b) initial condition (at the same time exact solution at t = 1), (d)&(f) numerical solutions at t = 1 with ACSI and conventional time stepping, respectively, (c)&(e) corresponding error maps.

Patrick Jenny / Journal of Computational Physics (2019)

16

error 2-norm

0.1

-2

N ACTI conv. time stepping

0.01

0.001

25

50 N1=N2

100

Fig. 13. 2D advection with a = (1, 1)T and periodic boundary conditions: Error 2-norms at t = 1 as functions of N1 = N2 ; the maximum CFL number was 0.5 for both ACTI and conventional time stepping.

4.4. Two Phase Flow in Porous Medium with Channels As a final test case, transport of a tracer u in a homogeneous porous medium with a single fracture is considered. Volumetric flux q is related to the pressure p via Darcy’s law as q

=

k − ∇p, μ

(23)

where k and μ are permeability of the medium and dynamic fluid viscosity; for simplicity we refer to the single parameter λ = k/μ. Incompressible fluid and incompressible porous material are assumed, which leads to the conservation equation ∇·q

=

s,

(24)

and together with Darcy’s law one obtains the elliptic pressure (or flow) equation ∇ · (λ∇p) =

−s.

(25)

Note that s is a volumetric source term, representing e.g. fluid injection (s > 0) and extraction (s < 0). Once the flow problem is solved, one can readily compute the tracer velocity a as a

= −

k ∇p, φμ

(26)

where φ is the porosity, which then can be employed to solve the tracer transport equation ∂u + ∇ · (a u) = ∂t

su . φ

(27)

Note that u is the tracer field and su a volumetric tracer source. These flow and transport equations are solved on a square domain Ω = [0, 1] × [0, 1] with no-flow Neumann boundary conditions everywhere (normal component of the pressure gradient is zero). A Cartesian equidistant 100 × 100 grid is considered with a maximum CFL number of 0.5. The source in the flow equation (25) is ⎧ ⎪ 1002 φ if x ∈ [0, 0.01] × [0, 0.01] (in cell Ω1,1 ) ⎪ ⎪ ⎨ 2 φ if x ∈ [0.99, 1] × [0.99, 1] (in cell Ω100,100 ) −100 s(x) = ⎪ (28) ⎪ ⎪ ⎩ 0 else and in the transport equation equation (27) it is ⎧ ⎪ 1002 φ U1,1 if x ∈ [0, 0.01] × [0, 0.01] ⎪ ⎪ ⎨ 2 φ U if x ∈ [0.99, 1] × [0.99, 1] −100 su (x) = ⎪ 100,100 ⎪ ⎪ ⎩ 0 else

(in cell Ω1,1 ) (in cell Ω100,100 )

(29)

Patrick Jenny / Journal of Computational Physics (2019)

17

production

fracture

homogeneous porous medium

tracer injection Fig. 14. Tracer transport in 2D porous medium with single fracture: Sketch of test case.

test case

end time tend

max. CFL number

1D advection 1D inviscid Burger 2D advection 2D fractured medium

1.0 0.25 1.0 0.5

1.0 1.0 0.5 0.5

cell updates without ACTI cell updates with ACTI

flux computations without ACTI flux computations with ACTI

34 68 35 20

25 50 30 20

Table 2. Ratios of required cell updates and flux computations without and with ACTI for the 1D advection-, 1D inviscid Burger-, 2D advection- and 2D fractured medium test cases.

with U1,1 = 1; U100,100 , on the other hand, results from the solution of the transport problem. The temporally constant parameter λ is defined as  100 if x ∈ [0.5, 0.51] × [0.25, 0.75] (representing a single fracture) λ(x) = (30) 1 else (representing a homogeneous porous medium) and the initial condition for the transport problem is u(x, t = 0) = 0. A sketch of the test case is shown in Fig. 14. Two-point central discretization is employed to estimate the pressure gradients at the cell interfaces and for λ at the cell interfaces harmonic averaging is employed. BiCGSTAB (Stabilized BiConjugate Gradient) from the Eigen library was used to solve the resulting linear system for the pressure, from which the wave speeds at the cell interfaces are calculated and then used for the transport problem. Figure 15 shows concentration maps (u-fields) for t ∈ {0.3, 0.4, 0.5}; on the left obtained with ACTI and on the right with conventional time stepping. At all times the differences are minute, which is also depicted in Fig. 16 for t ∈ {0.3, 0.5}. Also for this more relevant test case it was found that ACTI is much more efficient than conventional time stepping, i.e., 20 times fewer cell updates and flux calculations are required. The ratio of required computations without and with ACTI is also found in table 2 together with those for the two 1D cases and the 2D advection case. The reason for this speedup becomes evident from Fig. 17, which shows the cell based CFL numbers for the ACTI- and conventional time stepping simulations (left and right maps, respectively). While with ACTI almost all cell CFL numbers lie between 0.25 and 0.5, the majority of the cell CFL numbers are extremely tiny in the case of conventional time stepping. Figure 18 shows a surface plot of the levels used for ACTI. It can be observed that the ratio of global to smallest cell time step size is 29 . As expected, the levels are highest for those cells located in the fracture or which are near the injection or production sites. 5. Discussions and Conclusions A novel adaptive conservative time integration (ACTI) scheme for finite volume methods is presented, which proved to be of almost second temporal and spatial order. Results with ACTI of various linear and non-linear 1D

Patrick Jenny / Journal of Computational Physics (2019)

18

1

1

1 0.8

0.6

0.6

0.6 x2

0.8

x2

0.8

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0.2

0.4

x1

0.6

0.8

1

0

0

0

at t = 0.3 with ACTI

0.2

0.4

x1

0.6

0.8

1

at t = 0.3 without ACTI

1

1

1 0.8

0.6

0.6

0.6 x2

0.8

x2

0.8

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0.2

0.4

x1

0.6

0.8

1

0

0

0

at t = 0.4 with ACTI

0.2

0.4

x1

0.6

0.8

1

at t = 0.4 without ACTI

1

1

1 0.8

0.6

0.6

0.6 x2

0.8

x2

0.8

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0.2

0.4

x1

0.6

0.8

at t = 0.5 with ACTI

1

0

0

0

0.2

0.4

x1

0.6

0.8

1

at t = 0.5 without ACTI

Fig. 15. Tracer transport in 2D porous medium with single fracture: Shown are u-maps at t ∈ {0.3, 0.4, 0.5} obtained with ACTI (left) and without ACTI (right).

Patrick Jenny / Journal of Computational Physics (2019)

1

1

0.1 0.05

0.8

0.2

0.4

x1

0.6

0.8

1

-0.1

0.2

-0.15 0

-0.05

x2

x2 0

0

0.4

-0.1

0.2

0.05

0.6

-0.05

0.4

0.1

0.8

0

0.6

19

-0.2

0

-0.15 0

0.2

difference at t = 0.3

0.4

x1

0.6

0.8

1

-0.2

difference at t = 0.5

Fig. 16. Tracer transport in 2D porous medium with single fracture: Shown are differences in u between results obtained with ACTI and without ACTI at t ∈ {0.3, 0.5}.

1

0.5

1

0.45 0.8

0.4

0.8

0.35 0.6

0.3

0.6

0.4

x2

x2

0.25

0.4

0.2 0.15

0.2

0.2

0.1 0.05

0

0

0.2

0.4

x1

0.6

0.8

1

0

0

0

0.2

with ACTI

0.4

x1

0.6

0.8

1

without ACTI

Fig. 17. Tracer transport in 2D porous medium with single fracture: The maps show the cell CFL numbers for the simulations with ACTI (left) and with conventional time stepping (right).

9 8 7 6 5 4 3 2 1 1 0.8 0

0.2

0.6 0.4 x1

0.4 0.6

0.8

x2

0.2 1 0

Fig. 18. Tracer transport in 2D porous medium with single fracture: The surface plot depicts the levels for all cells as used for ACTI.

20

Patrick Jenny / Journal of Computational Physics (2019)

and 2D test cases are almost indistinguishable form the corresponding ones with conventional time stepping (with the same maximum CFL number), but much fewer cell updates and flux calculations were required. Since the extra overhead of ACTI only consists in sorting cell lists (these lists contain a fraction of all grid cells), this translates into dramatic computational gains. The new ACTI method is very simple to implement and proved to be very robust. Based on these encouraging findings, one may conclude that for many problems, e.g. for transport in sequentially coupled reservoir simulation [19, 20], ACTI renders implicit time stepping unnecessary. Moreover, sequentially coupled flow and transport problems are interesting candidates for multi-scale methods [21, 22, 23]; especially if heterogeneous coefficient distributions are involved. For example, for the flow problem the multi-scale finite volume (MSFV) method [24, 25, 26, 27, 28, 29, 30, 31] proved to be very efficient and transport can be tackled independently. To allow for arbitrarily large time steps to solve multi-phase transport problems, an unconditionally stable non-linear solver was devised [32]. While very effective to save computational work, it does not tackle the issue of large truncation errors. Instead, ACTI could be employed and especially in combination with the adaptive multi-scale method for transport by Lee et al. [33] further improvements of computational efficiency and accuracy would be possible. Acknowledgements The author thanks Piyush Panchal and Kristof Sarosi (ETH Zurich), who implemented an early version of the algorithm and performed preliminary studies. Further, the author thanks Daniel W. Meyer (ETH Zurich), Rasim Hasanzade and Hamdi Tchelepi (Stanford University) for their constructive feedback. References ¨ [1] R. Courant, K. Friedrichs, H. Lewy, Uber die partiellen Differenzengleichungen der mathematischen Physik, Mathematische Annalen 100 (1928) 32–74. [2] S. Osher, R. Sanders, Numerical approximations to nonlinear conservation laws with locally varying time and space grids, Math. Comp. 41 (1983) 321–321. [3] M. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, Journal of Computational Physics 53 (1984) 484–512. [4] C. Dawson, R. Kirby, High resolution schemes for conservation laws with locally varying time steps, SIAM J. Sci. Comput. 22 (2001) 2256—2281. [5] W. L. Kleb, J. T. Batina, M. H. Williams, Temporal adaptive Euler/Navier-Stokes algorithm involving unstructured dynamic meshes, AIAA Journal 30 (1992) 1980–1985. [6] X. Zhang, J. Trepanier, M. Reggio, R. Camarero, Time-accurate local time stepping method based on flux updating, AIAA Journal 32 (1994) 1926–1929. [7] J. Nutaro, B. P. Zeigler, R. Jammalamadaka, S. Akerkar, Discrete event solution of gas dynamics within the EVS framework, in: Computational Science ICCS 2003, 2003, pp. 319–328. [8] E. M. Constantinescu, A. Sandu, Multirate timestepping methods for hyperbolic conservation laws., J. Sci. Comput. 33 (2007) 239–278. [9] P. K. Fok, A linearly fourth order multirate Runge-Kutta method with error control, J. Sci. Comput. 66 (2015) 177–195. [10] B. Seny, J. Lambrechts, T. Toulorge, V. Legat, J. F. Remacle, An efficient parallel implementation of explicit multirate Runge-Kutta schemes for discontinuous Galerkin computations., Journal of Computational Physics 256 (2014) 135–160. [11] Q. Shao, S. Matthai, L. Gross, Efficient modelling of solute transport in heterogeneous media with discrete event simulation, Journal of Computational Physics 394 (2019) 134–155. [12] J. R. Cavalcanti, M. Dumbser, D. da Motta-Marques, C. R. Fragoso, A conservative finite volume scheme with time-accurate local time stepping for scalar transport on unstructured grids, Advances in Water Resources 86 (2015) 217–230. [13] L. D. Carciopolo, L. Bonaventura, A. Scotti, L. Formaggia, A conservative implicit multirate method for hyperbolic problems, Computational Geosciences 23 (2019) 647–664. [14] L. D. Carciopolo, A. Scotti, H. Hajibeygi, L. Formaggia, Conservative multirate multiscale simulation of multiphase flow in heterogeneous porous media, submitted to Journal of Computational Physics (2019). [15] B. Koren, A robust upwind discretisation method for advection, diffusion and source terms, in: C. B. Vreugdenhil, B. Koren (Eds.), Numerical Methods for Advection-Diffusion Problems, Braunschweig: Vieweg, 1993, p. 117. [16] B. van Leer, Towards the ultimate conservative difference scheme, a second order sequel to Godunov’s method, Journal of Computational Physics 32 (1979) 101–136. [17] R. J. Leveque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002. [18] C. A. R. Hoare, Algorithm 64: Quicksort, Communications of the ACM 4 (1961) 321. [19] A. Montcorge, H. A. Tchelepi, P. Jenny, Modified sequential fully implicit scheme for compositional flow simulation, Journal of Computational Physics 337 (2017) 98–115. [20] A. Montcorge, H. Tchelepic, P. Jenny, Sequential fully implicit formulation for compositional simulation using natural variables, Journal of Computational Physics 371 (2018) 690–711. [21] T. Hou, X. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, Journal of Computational Physics 134 (1997) 169–189.

Patrick Jenny / Journal of Computational Physics (2019)

21

[22] T. Arbogast, S. Bryant, Numerical subgrid upscaling for waterflood simulations, in: SPE 66375, presented at the SPE Symp. on Reservoir Simulation, 2001. [23] Y. Efendiev, T. Hou, X. Wu, Convergence of a nonconformal multiscale finite element method, SIAM J. Numer. Anal. 37 (2000). [24] P. Jenny, S. H. Lee, H. A. Tchelepi, Multi-scale finite-volume method for elliptic problems in subsurface flow simulation, Journal of Computational Physics 187 (2003) 47–67. [25] P. Jenny, S. H. Lee, H. Tchelepi, Adaptive multiscale finite-volume method for multi-phase flow and transport in porous media, SIAM Journal for Multiscale Modeling and Simulation 3 (2005) 50–64. [26] P. Jenny, I. Lunati, Multi-scale finite-volume method for elliptic problems with heterogeneous coefficients and source terms, PAMM - Proc. Appl. Math. Mech. 6 (2006) 485–486. [27] I. Lunati, P. Jenny, Multiscale finite-volume method for compressible multiphase flow in porous media, Journal of Computational Physics 216 (2006) 616–636. [28] P. Jenny, S. H. Lee, H. A. Tchelepi, Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media, Journal of Computational Physics 217 (2006) 627–641. [29] H. Hajibeygi, G. Bonfigli, M. A. Hesse, P. Jenny, Iterative multiscale finite-volume method, Journal of Computational Physics 227 (2008) 8604–8621. [30] H. Hajibeygi, P. Jenny, Multiscale finite-volume method for parabolic problems arising from compressible multiphase flow in porous media, Journal of Computational Physics 228 (2009) 5129–5147. [31] P. Jenny, I. Lunati, Modeling complex wells with the multi-scale finite-volume method, Journal of Computational Physics 228 (2009) 687–702. [32] P. Jenny, H. Tchelepi, S. Lee, Unconditionally convergent nonlinear solver for hyperbolic conservation laws with s-shaped flux functions, Journal of Computational Physics 228 (2009) 7497–7512. [33] S. Lee, H. Zhou, H. Tchelepi, Adaptive multiscale finite-volume method for nonlinear multiphase transport in heterogeneous formations, Journal of Computational Physics 228 (2009) 9036–9058.

 

 


 
 
 
 
 








  
 
 



 
 



 



 
 
 
 





  
 
















There is no conflict of interest.