Python framework for hp-adaptive discontinuous Galerkin methods for two-phase flow in porous media

Applied Mathematical Modelling 67 (2019) 179–200

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Python framework for hp-adaptive discontinuous Galerkin methods for two-phase flow in porous media Andreas Dedner a,∗, Birane Kane b, Robert Klöfkorn c,1, Martin Nolte d a

Mathematics Institute, Zeeman Building, University of Warwick, Coventry CV4 7AL, UK Institute of Applied Analysis and Numerical Simulation, University of Stuttgart, Pfaffenwaldring 57 70569 Stuttgart, Germany c International Research Institute of Stavanger, P.O. Box 8046 N-4068 Stavanger, Norway d Department for Applied Mathematics, Albert-Ludwigs-University Freiburg, Hermann-Herder-Str. 10, 79104 Freiburg, Germany b

a r t i c l e

i n f o

Article history: Received 18 April 2018 Revised 24 September 2018 Accepted 10 October 2018 Available online 23 October 2018 Keywords: Discontinuous Galerkin hp-adaptivity Porous media two-phase flow IMPES Dune Python

a b s t r a c t In this paper we present a framework for solving two-phase flow problems in porous media. The discretization is based on a Discontinuous Galerkin method and includes local grid adaptivity and local choice of polynomial degree. The method is implemented using the new Python frontend Dune-FemPy to the open source framework Dune. The code used for the simulations is made available as Jupyter notebook and can be used through a Docker container. We present a number of time stepping approaches ranging from a classical IMPES method to a fully coupled implicit scheme. The implementation of the discretization is very flexible allowing to test different formulations of the two-phase flow model and adaptation strategies. © 2018 Elsevier Inc. All rights reserved.

1. Introduction Simulation of multi-phase flows and transport processes in porous media requires careful numerical treatment due to the strong heterogeneity of the underlying porous medium. The spatial discretization requires locally conservative methods in order to be able to follow small concentrations [1]. Discontinuous Galerkin (DG) methods, Finite Volume methods and Mixed Finite Element methods are examples of discretization techniques achieving local conservation at the element level [2]. Application of DG methods to incompressible two-phase flow started within the framework provided by a decoupled approach called Implicit Pressure Explicit Saturation (IMPES) where first a pressure equation is solved implicitly and then the saturation is advanced by an explicit time stepping scheme. Upwinding, slope limiting techniques, and sometimes H(div)projection were required in order to remove unphysical oscillations and to ensure convergence to a solution. In the fully implicit and fully coupled approach, the mass balances are usually discretized in time by the implicit Euler method, resulting in a fully coupled system of nonlinear equations that has to be solved at each time step. The main ∗

Corresponding author. E-mail addresses: [email protected] (A. Dedner), [email protected] (B. Kane), [email protected] (R. Klöfkorn), [email protected] (M. Nolte). URL: https://warwick.ac.uk/fac/sci/maths/people/staff/andreas_dedner/ (A. Dedner), http://www.ians.uni-stuttgart.de/nmh/kane (B. Kane), https://aam.uni-freiburg.de/mitarb/ehemalige/nolte/index.html (M. Nolte) 1 www.iris.no https://doi.org/10.1016/j.apm.2018.10.013 0307-904X/© 2018 Elsevier Inc. All rights reserved.

180

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

advantage of a fully implicit scheme is the possibility of using significantly larger time step sizes, which can be crucial in view of long-term scenarios like atomic waste disposal. Commonly, rather simple, yet robust space-discretization schemes, like cell-centered or vertex-centered finite volume schemes, are used [1,3]. Fully implicit DG schemes have been proposed in [4] and [5], where the schemes are formulated in two space dimensions for incompressible fluid phases and numerical tests are performed without any kind of adaptivity. Bastian introduced in [6] a fully coupled symmetric interior penalty DG scheme for incompressible two-phase flow based on a wetting-phase potential and capillary potential formulation. Discontinuity in capillary-pressure functions is taken into account by incorporating the interface conditions into the penalty terms for the capillary potential. Heterogeneity in absolute or intrinsic permeability is treated by weighted averages. A higher-order diagonally implicit Runge–Kutta method in time is used and there is neither post processing of the velocity nor slope limiting. Only piecewise linear and piecewise quadratic functions are employed and no adaptive method is considered. A general abstract framework allowing for an a-posteriori estimator for porous-media two-phase flow problem was introduced by Vohralik and Wheeler [7]. This paved the way for an h-adaptive strategy for homogeneous two-phase flow problems [8]. However, it has not been applied to DG methods so far. Finally, Darmofal et al. introduced recently a space-time discontinuous Galerkin h-adaptive framework for 2d reservoir flows. Implicit estimators are derived through the use of dual problems [9,10] and a higher-order discretization is performed on anisotropic, unstructured meshes. Application to 3d problems and hp-adaptive strategies haven’t been considered yet. In this paper, we implement and evaluate numerically interior penalty DG methods for incompressible, immiscible, twophase flow. We consider strongly heterogeneous porous media, anisotropic permeability tensors and discontinuous capillarypressure functions. We write the system in terms of a phase-pressure/phase-saturation formulation. Adams–Moulton schemes of first or second order in time are combined with an Interior Penalty DG discretization in space. This implicit space time discretization leads to a fully coupled nonlinear system requiring to build a Jacobian matrix at each time step for the Newton-Raphson method. A thorough analysis of the numerical schemes providing stability estimates and a-priori error estimates with respect to the L2 (H1 ) norm for the pressure and the L∞ (L2 ) ∩ L2 (H1 ) norm for the saturation is carried out in [11,12]. This paper extends our previous work in [13,14] and [15]. We consider here new hp-adaptive strategies and compare the fully implicit scheme with the iterative IMPES scheme and the implicit iterative scheme. The iterative IMPES belongs to a class of schemes based on the L-scheme (or stabilised Picard), see [16]. It has been already applied to two-phase flow in porous media in [17,18]. Also, especially in combination with higher order DG in [19]. We refer to [11,20,21] for similar iterative IMPES schemes. The implicit iterative scheme is based on the iterative IMPES approach presented in [22] and treats the capillary pressure term implicitly to ensure stability. We also provide a more comprehensive model framework allowing to conveniently implement and compare various two-phase flow formulations. The implementation is based on the open-source PDE software framework Dune-FemPy, which is a Python frontend for Dune-Fem [23] based on the new Dune-Python module [24] and which adds support for the Unified Form Language [25]. It allows for a compact, legible presentation of the different discretizations under consideration. We combine Dune-FemPy with Jupyter [26] and Docker [27] to ensure reproducibility of our numerical experiments. The adaptive grid implementation is based on Dune-Alugrid [28] and parts of the stabilization mechanisms used are provided by Dune-Fem-DG [29]. The rest of this document is organised as follows. In Section 2, we provide some more details on the software package used to perform the experiments in this paper. In Section 3, we describe the two-phase flow model. The DG discretization is introduced in Section 4. Numerical examples are provided in Section 5. Conclusions are drawn in the last section.

2. Python frontend to

Dune

The experiments presented in this paper are based on a newly developed Python frontend to the well established software framework Dune. Dune is a open source software package for solving partial differential equations using grid based numerical schemes. The hallmarks of Dune are its modularity and flexibility. It is based around a set of core modules which form the building blocks for user modules. In most cases these user modules will be based on a discretization module which contains interfaces and realizations of common features of numerical schemes, such as grid functions, discrete functions and operators. The results shown here are based on the discretization module Dune-Fem [23] and an extension module Dune-Fem-DG [29] which provides additional functionality specifically for using discontinuous Galerkin methods. Dune is written in C++ using static polymorphism to implement the interfaces and realizations. It uses many advanced features and has proven to have a steep learning curve especially for users not used to high level C++ programming. To mitigate this, a Python frontend is being implemented to allows high level access to the all Dune features [24]. This frontend is designed to maintain the full flexibility and modularity of Dune while at the same time making it straightforward to control pre and post processing and general program flow using the scripting language Python. Flexibility is maintained by just in time compilation of C++ code and performance is achieved by only adding a thin layer to the C++ objects when passing them to Python which is removed when these objects are passed back into a C++ algorithm. Experiments have shown that there is no performance impact when using the high level Dune-Fem objects through the Python module Dune-FemPy [30].

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

181

In addition to implementing bindings to all the high level discretization interfaces provided by Dune-Fem, Dune-Fem-DG, Dune-FemPy also allows to use the UFL to describe the mathematical model, i.e., the bilinear forms. The code presented in the appendices of this paper make heavy use of UFL for describing not only the PDE model but also other parts of the algorithm, e.g., the adaption indicators. 3. Mathematical model This section introduces the mathematical formulation of a two-phase porous-media flow. In all that follows, we assume that the flow is immiscible and incompressible with no mass transfer between phases. 3.1. Two-phase flow formulation Let  be a polygonal bounded domain in Rd , d ∈ {2, 3}, with Lipschitz boundary ∂  and let T ∈ R+ . The flow of the wetting phase (e.g. water) and the nonwetting phase (e.g. oil, gas) is described by Darcy’s law and the continuity equation (i.e. balance of mass) for each phase α ∈ {w, n}[31]. In all that follows, we denote with subscript w the wetting phase and with subscript n the nonwetting phase. The unknown variables are the phase pressures pw , pn :  × (0, T ) → R and the phase saturations sw , sn :  × (0, T ) → R. For each phase α ∈ {w, n}, the Darcy velocity vα :  × (0, T ) → Rd is given by

vα = −λα K(∇ pα − ρα g )

in  × (0, T )

(1)

where λα :  × (0, T ) → R is the phase mobility, K :  → Rd×d is the absolute or intrinsic permeability tensor of the porous medium, ρα :  × (0, T ) → R is the phase density, and g ∈ Rd is the constant gravitational vector. Phase mobilities λα :  × (0, T ) → R are defined by

λα =

kr α

μα

,

α ∈ {w, n},

(2)

where μα is the constant phase viscosity and krα :  × (0, T ) → R is the relative permeability of phase α . The relative permeabilities are functions that depend nonlinearly on the phase saturation (i.e. krα = krα (sα )). Models for the relative permeability are the Van-Genuchten model [32] and the Brooks–Corey model [33]. For example, in the Brooks–Corey model,

krw (sn,e ) = (1 − sn,e )

2+3θ

θ

,

krn (sn,e ) = (sn,e )2 (1 − (1 − sn,e )

2+θ

θ

),

(3)

where the effective saturation sα , e is

sα ,e =

sα − sα ,r , 1 − sw,r − sn,r

∀α ∈ {w, n}.

(4)

Here, sα , r , α ∈ {w, n}, are the phase residual saturations. The parameter θ ∈ [0.2, 3.0] is a result of the inhomogeneity of the medium. For each phase α ∈ {w, n}, the balance of mass yields the saturation equation

φ

∂ ( ρα s α ) + ∇ · ( ρα v α ) = ρα q α , ∂t

(5)

where φ :  → R is the porosity, qα :  × (0, T ) → R is a source or sink term (e.g. wells located inside the domain in the case of a reservoir problem). In addition to (1) and (5), the following closure relations must also be satisfied:

sw + sn = 1,

(6)

pn − pw = pc ( sn ),

(7)

where pc (sn ) :  × (0, T ) → R is the capillary pressure, a Lipschitz continuous function of the phase saturation. For the Brooks–Corey formulation,

pc (sn,e ) = pd (1 − sn,e )−1/θ .

(8)

Here, pd ≥ 0 is the constant entry pressure needed to displace the fluid from the largest pore. 3.2. Model A: Wetting-phase-pressure/nonwetting-phase-saturation formulation Considering the phases are incompressible (i.e. the densities ρ α are constant), we get a total fluid conservation equation by summing the two mass balance equations from (5),

φ

∂ ( sn + sw ) + ∇ · ( vn + vw ) = qn + qw . ∂t

182

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

Thanks to relation (6),

∇ · ( vn + vw ) = qn + qw . From relation (1) we have

−∇ · (λn K(∇ pn − ρn g ) + λw K(∇ pw − ρw g )) = qn + qw . The last closure relation (7) allows to write

−∇ · (λn K(∇ pc + ∇ pw − ρn g ) + λw K(∇ pw − ρw g )) = qn + qw . Finally,





−∇ · (λw + λn )K∇ pw + λn K∇ pc − (ρw λw + ρn λn )Kg = qw + qn . To complete our system, we consider as second equation the nonwetting phase conservation relation

φ

∂ sn + ∇ · vn = qn . ∂t

Using relation (1) and (7) yields

φ

    ∂ sn − ∇ · λn K(∇ pw − ρn g ) − ∇ · λn K∇ pc = qn . ∂t

We get therefore a system of two equations with two unknowns pw and sn ,





−∇ · (λw + λn )K∇ pw + λn K∇ pc − (ρw λw + ρn λn )Kg = qw + qn ,

φ

    ∂ sn − ∇ · λn K(∇ pw − ρn g ) − ∇ · λn K∇ pc = qn . ∂t

(9)

(10)

Substituting ∇ pc = p c (sn )∇ sn for ∇ pc as in [31,34], the system (9)-(10) becomes





−∇ · (λw + λn )K∇ pw + λn p c K∇ sn − (ρw λw + ρn λn )Kg = qw + qn ,

φ

    ∂ sn − ∇ · λn K(∇ pw − ρn g ) − ∇ · λn p c K∇ sn = qn . ∂t

(11)

(12)

In order to have a complete system, we add appropriate boundary and initial conditions. Thus, we assume that the boundary of the system is divided into disjoint sets such that ∂  = D ∪ N . We denote by ν the outward normal to ∂  and set

pw (·, 0 ) = p0w (· ),

sn (·, 0 ) = s0n (· ),

pw = pw,D , vα · ν = Jα ,

sn = sD ,  Jt = Jα ,

in , on D × (0, T ), on N × (0, T ).

α ∈{w,n}

Here, Jα ∈ R, α ∈ {w, n}, is the inflow, s0n , p0w , sD , and pw, D are real numbers. In order to make pw uniquely determined, we require D = ∅. 3.3. General model framework We provide here a unified model framework allowing for the representation of the models introduced in the previous sections,





−∇ · A pp (s )∇ p + A ps (s )∇ s + G p (s ) = q p ,

(13)

  ∂t s − ∇ · Asp (s )(∇ p − Pg ) + Ass (s )∇ s + Gs (s ) = qs .

(14)

The model is described once the physical parameter functions A, G, and Pg are known. For Model A (i.e. (11) and (12)), we have for example p = pw , s = sn and

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

183

A ps (s ) = λn (s ) p c (s )K,

A pp (s ) = (λn (s ) + λw (s ))K,

Ass (s ) = λn (s ) p c (s )K,

Asp (s ) = λn (s )K,

G p (s ) = −(ρw λw (s ) + ρn λn (s ))Kg,

Gs ( s ) = 0, Pg = ρn g, q p = qw + qn ,

qs = qn .

4. Discretization In this section, we provide a discretization framework for a two-phase flow in a strongly heterogeneous and anisotropic porous medium. 4.1. Space discretization Let Th = {E } be a family of non-degenerate, quasi-uniform, possibly non-conforming partitions of  consisting of Nh elements (quadrilaterals or triangles in 2d, tetrahedrons or hexahedrons in 3d) of maximum diameter h. Let h be the union of the open sets that coincide with internal interfaces of elements of Th . Dirichlet and Neumann boundary interfaces are collected in the set Dh and Nh , respectively. Let e denote an interface in h shared by two elements E− and E+ of Th ; we associate with e a unit normal vector ν e directed from E− to E+ . We also denote by |e| the measure of e. The discontinuous finite element space is Dr (Th ) = {v ∈ L2 () : v|E ∈ PrE (E ) ∀E ∈ Th }, with r = (rE )E∈Th , PrE (E ) denotes QrE (resp. PrE ) the space of polynomial functions of degree at most rE ≥ 1 on E (resp. the space of polynomial functions of total degree rE ≥ 1 on E). We approximate the pressure and the saturation by discontinuous polynomials of total degrees r p = (r p,E )E∈Th and rs = (rs,E )E∈Th , respectively. For any function q ∈ Dr (Th ), we define the jump operator · and the average operator { · } over the interface e: ∀e ∈ h , q := qE− νe − qE+ νe , {q} := 12 qE− + 12 qE+ , ∀e ∈ ∂ , q := qE− ν , {q} := qE− . In order to treat the strong heterogeneity of the permeability tensor, we follow [35] and introduce a weighted average operator { · }ω : ∀e ∈ h , {q}ω = ωE− qE− + ωE+ qE+ , ∀e ∈ ∂ , {q}ω = qE− . +

−

The weights are ωE− = k+k+k− , ωE+ = k+k+k− with k− = νeT KE− νe and k+ = νeT KE+ νe . Here, KE− and KE+ are the permeability tensors for the elements E− and E+ , respectively. The derivation of the semi-discrete DG formulation is standard (see [6,35,36]). First, we multiply each equation of (13) and (14) by a test function and integrate over each element, then we apply the Green formula to obtain the semidiscrete weak DG formulation. The bulk integrals are thus given by:

B p (( p, s ), ϕ ; s¯) =

 E∈Th

Bs (( p, s ), ϕ ; s¯) =



 E

 E

E∈Th



A pp (s¯)∇ p + A ps (s¯)∇ s · ∇ϕ +



 E∈Th

E

Asp (s¯)(∇ p − Pg ) + Ass (s¯)∇ s · ∇ϕ +

G p (s¯) · ∇ϕ −

E∈Th

 E∈Th



E

E

Gs (s¯) · ∇ϕ −

q pϕ , 

E∈Th

E

qs ϕ .

The consistency terms on the skeleton are



C p ( ( p, s ), ϕ ; s ) =

e∈ h ∪ Dh ∪ Nh



Cs ( ( p, s ), ϕ ; s ) =

e∈ h ∪ Dh ∪ Nh

 

e

e

{App (s )∇ p + Aps (s )∇ s + G p (s )}ω · ϕ , {Asp (s )(∇ p − Pg ) + Ass (s )∇ s + Gs (s )}ω · ϕ 

with a given s¯ coming from the linearization. To stabilize the scheme, we define interior penalty terms on the skeleton with σ > 0 a given constant:

S p ( p, ϕ ) = σ Ss (s, ϕ ) = σ

 

e∈ h ∪ Dh

e

e∈ h ∪

e

  h D

γep  p · ϕ , γes s · ϕ .

We follow the suggestions from [37] and choose σ = r (r + 1 ) where r is the highest polynomial degree of the discrete p spaces. For Model A, the penalty terms γe and γes are given by

184

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

2k+ k− |e| × , k+ + k− min(|E+ |, |E− | ) 2k+ k− |e| = max(δs+ , δs− ) + × k + k− min(|E+ |, |E− | )

γep = max(δ +p , δ −p ) γes

where the terms δ p and δ s are the largest eigenvalues of App (0.5) and Ass (0.5), respectively. The two bilinear forms thus are

Fp (( p, s ), ϕ ; s¯) = B p (( p, s ), ϕ ; s¯) − C p (( p, s ), ϕ ; s¯) + S p ( p, ϕ ), Fs (( p, s ), ϕ ; s¯) = Bs (( p, s ), ϕ ; s¯) − Cs (( p, s ), ϕ ; s¯) + Ss ( p, ϕ ). 4.2. Time stepping Denoting with (pi , si ) the approximation to the solution in the discrete function space at some point in time ti , we use a simple one step scheme to advance the solution (pi , si ) at time ti to ( pi+1 , si+1 ) at the next point in time t i+1 = t i + τ based on

Fp (( pi+1 , si+1 ), ϕ ; s¯) = 0, 

(si+1 − si )ϕ + τ Fsα (( pi+1 , si+1 ), ϕ ; s¯) = 0,

(15)

(16)

defining for a given constant α ∈ [0, 1] the bilinear form

Fsα (( p, s ), ϕ ; s¯) = (1 − α )Fs (( pi , si ), ϕ ; si ) + α Fs (( pi+1 , si+1 ), ϕ ; s¯).

(17)

The starting points of the iteration (p0 , s0 ) are taken as L2 projections of the functions given by the initial conditions into the discrete space. In our tests we have always used α = 1 since we were more interested in investigating the influence of s¯. We also used a fixed time step τ throughout the whole course of the simulation although varying time steps can be easily used as well. Different choices for s¯ lead to different approaches for handling the nonlinearities in the pressure. We tested five different approaches described in the following: Linear For this approach we simply take s¯ = si leading to a forward Euler time stepping scheme. Implicit Taking s¯ = si+1 leads to a backward Euler scheme. The resulting fully coupled system is solved iteratively using a Newton method. Iterative This is similar to the previous approach, replacing the Newton method by an outer fixed point iteration to solve the system. We define s¯k = si+1,k with s¯0 = si and in each step of the iteration we therefore solve for k ≥ 0:

Fp (( pi+1,k+1 , si+1,k+1 ), ϕ ; si+1,k ) = 0, 

(si+1,k+1 − si )ϕ + τ Fsα (( pi+1,k+1 , si+1,k+1 ), ϕ ; si+1,k ) = 0.

(18)

(19)

IMPES-iterative This results in an iterative scheme using an IMPES approach. This is similar to the previous approach except that in each step of the iteration we solve

Fp (( pi+1,k+1 , si+1,k ), ϕ ; si+1,k ) = 0, 

(si+1,k+1 − si )ϕ + τ Fsα (( pi+1,k+1 , si+1,k+1 ), ϕ ; si+1,k ) = 0.

(20)

(21)

We note that the difference between the IMPES-iterative and the Iterative is the evaluation of the gradient of the saturation in Fp . For (18) we evaluate it at the step k + 1 and for (20) at the step k, respectively. IMPES Finally, we use a classical IMPES approach, which is similar to the previous without carrying out the iteration. The saturation in the pressure equation is taken explicitly and the new pressure is used in the saturation equation (in contrast to the first approach where the old pressure is used):

Fp (( pi+1 , si ), ϕ ; si ) = 0, 

(si+1 − si )ϕ + τ Fsα (( pi+1 , si+1 ), ϕ ; si ) = 0.

(22)

(23)

With the exception of the first and the last approach, all methods use an iteration to obtain a fixed point to the fully implicit equation

Fp (( pi+1 , si+1 ), ϕ ; si+1 ) = 0,

(24)

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200



(si+1 − si )ϕ + τ Fsα (( pi+1 , si+1 ), ϕ ; si+1 ) = 0.

185

(25)

In the third and the fourth method, this is achieved using an outer iteration (based on the first or the last method, respectively) while the second method uses a Newton method. To make the approaches easier to compare, we use the same stopping criteria for the iteration in all three cases. We take ( pi+1 , si+1 ) = ( pi+1,l , si+1,l ) with l such that

si+1,l − si+1,l−1 L2 () < toliter sl−1 L2 () .

(26)

We use a value of toliter = 3 · 10−2 to stop the iteration when the relative change between two steps is less then three percent. 4.3. Adaptivity Different adaptive strategies are possible depending on how elements are refined/coarsened; whether the elements should be p-refined or h-refined; when should the refinement process be stopped (e.g. maximum level of refinement, stopping criterion). Keeping this in focus, we provide in this section a brief introduction to different adaptive strategies implemented and tested in this work. In all that follows, the parameters maxpoldeg and maxlevel refer respectively to the maximum polynomial degree and the maximum level of refinement allowed. 4.3.1. Error indicator In the sequel, we implement an explicit estimator originally designed for non-steady convection-diffusion problems. A thorough analysis is available in [38]. Applying the estimator to the phase conservation Eq. (14) yields:

ηE2 = h2E Rvol 2L2 (E ) + +

 e∈∂ E∩∂ 



 1  1 2 2 he Re2 L2 (e ) + Re1 L2 (e ) 2 h e h e∈

he Re2 L2 (e ) + 2

1 Re1 2L2 (e) he





.

(27)

Here Rvol is the interior residual indicating how accurate the discretized solution satisfies the original PDE at every interior point of the domain,

Rvol = qs − φ

  ∂s + ∇ · Asp (s )(∇ p − Pg ) + Ass (s )∇ s + Gs (s ) . ∂t

The term Re1 is the numerical zero order inter-element (resp. Dirichlet boundary condition) residual depending on the jump of the discrete solution at the elements interfaces (resp. at the Dirichlet boundary), hence reflecting the regularity of the DG approximation (resp. the accuracy of the approximation on the Dirichlet boundary),

 σ γes s Re1 = σ γes (sD − s )

if e ∈ h , if e ∈ D .

The term Re2 is the first order numerical inter-element residual (resp. Neumann boundary condition residual) depending on the jump of numerical approximation of the normal flux at the elements interfaces (resp. at the Neumann boundary). It also allows to assess the regularity of the DG approximation (resp. the accuracy of the approximation on the Neumann boundary),



Re2 =

Asp (s )(∇ p − Pg ) + Ass (s )∇ s + Gs (s ) · νe

Jn + (Asp (s )(∇ p − Pg ) + Ass (s )∇ s + Gs (s ) ) · νe

if e ∈ h , if e ∈ N .

4.3.2. Adaptive strategies The indicator presented above will be used to drive adaptive algorithms. The h-adaptive algorithm is depicted in Algorithm 1. Given the error indicator ηEr,n defined in Eq. (27) for a polynomial degree r in time step n for each element E, we refine each element whose error indicator is greater than a refinement threshold value hT olEn and we coarsen elements where the indicator is smaller than the coarsening threshold 0.01 × hT olEn . In order to automatically compute the tolerance for refinement hT olEn used in each timestep during the simulation we choose the following approach. We pre-describe a tolerance for the initial adaptation such that the resulting refined grid looks satisfactory. Then we apply an equi-distribution strategy which aims to equally distribute the error contribution over all time steps and grid elements. As a result we compute hT olEn based on the initially computed error indicator,

186

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

Algorithm 1 h-adapt. Let ηEr,n be given ∀E ∈ Th for all E ∈ Th do 3: hE = diam(E ) if ηEr,n > hT olEn AND maxl evel > l evel (E ) then 4: 1:

2:

5: 6: 7: 8: 9: 10: 11:

h

hnew := 2E E else if ηEr,n < 0.01 × hT olEn AND l evel (E ) > 0 then hnew := 2hE E else hnew := hE E end if end for

hT olEn := tT ol

|

τn

Thn

|

with

tT ol :=

1  r,0 ηE . T

(28)

E∈Th

In the following we use ηEr = ηEr,n as abbreviation for ease of reading. The choice between increasing or decreasing the local polynomial order depends heavily on the value of an indicator ςE (ηEr , ηEr−1 ) where ηEr , E ∈ Th is a given error indicator and ηEr−1 is the same indicator evaluated for the L2 projection of the solution into a lower order polynomial space. The derivation of this L2 projection is quite straightforward due to the hierarchical aspect of the modal DG bases implemented. We considered a marking strategy based on these indicators

ςE := |ηEr − ηEr−1 |.

(29)

When this difference on a given element is non zero we expect the higher order to contribute to the accuracy of the scheme and keep or increase the given polynomial on that element, otherwise the polynomial order is decreased. Details are provided in Algorithm 2. Algorithm 2 p-adapt: markpDiff. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19:

Let ςE = |ηEr − ηEr−1 | be given ∀E ∈ Th for all E ∈ Th do rE := poldeg(E ) if ςE < ptol then if rE > 1 then rEnew := rE − 1 else rEnew := rE end if else if ςE > 100 × ptol then if rE < maxpoldeg then rEnew := rE + 1 else rEnew := rE end if else rEnew := rE end if end for

4.4. Stabilization Although due to the presence of the capillary pressure terms strong shocks do not occur in the numerical experiments carried out in this paper, the DG schemes needs stabilization to avoid unphysical values, such as negative saturation which would lead to an undefined state in Eq. (8). We follow the approach from [39] which has been initially proposed by Zhang and Shu in [40]. The general idea is to scale each polynomial on each element such that a constraint on minimum and maximum values of the saturation is respected. We define the following projection operator s : Dr (Th ) −→ Dr (Th ) with

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

187

Table 1 2d problem parameters.

 [-] Swr [-] Snr [-] θ [-] pd [Pa]

lens

\lens

0.39 0.1 0.00 2.0 50 0 0

0.40 0.12 0.00 2.70 755

Fig. 1. Geometry and boundary conditions for the DNAPL infiltration problem. The purple line in the right picture is described by x(σ ) from Eq. (33). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

 

s [s] · ϕ =

 

s˜ · ϕ

∀ϕ ∈ Dr

(30)

where on each element E of the grid we define a scaled saturation





s˜(x ) := χE s(x ) − s¯ + s¯

(31)

with s¯ being the mean value of s on element E. Note that by construction this leads to a mass conservative scheme since the scaled part has mean value zero and therefore s˜ has the same mean value as s. The scaling factor is

χE := min{1, |(s¯ − smin )/(s¯ − s(x ))|, |(smax − s¯)/(s¯ − s(x ))|} x∈E

(32)

for the combined set of all quadrature points E used for evaluation of the bilinear forms defined earlier, i.e. interior and surface integrals. The scaling limiter is applied after each Newton iteration for the implicit scheme and after each iteration of the iterative schemes. 5. Numerical experiments This section provides different numerical experiments aiming to demonstrate the effectiveness and robustness of the DG discretization of porous media flow models. All test are implemented with the hp-adaptive DG method described in the previous section using the different approaches for the time step computation. The maximal grid level was fixed to three and the maximal polynomial degree was also three. The main components of the code are described in some detail in Appendix B and provided as part of a Docker container as explained in Appendix A. We also show results based on some alternative approaches for example for the underlying model or for the adaptive strategy. These modifications to the python code are also provided in Appendix B. They demonstrate the flexibility of the Python code. 5.1. Problem setting A container is filled with two kinds of sand and saturated by water with density ρw = 10 0 0 kg m−3 and viscosity μw = 1 × 10−3 kg/ms. The dense non-aqueous phase liquid (DNAPL) considered in the experiment is Tetrachloroethylene with density ρn = 1460 kg/ms − 3 and viscosity μn = 9 × 10−4 kg/ms. We consider a two-dimensional DNAPL infiltration problem with different sand types and anisotropic permeability tensors. The material properties are detailed in Table 1. The bottom of the reservoir is impermeable for both phases. Hydrostatic conditions for the pressure pw and homogeneous Dirichlet conditions for the saturation sn are prescribed at the left and right boundaries. A flux of Jn = −5.137 × 10−5 ms−1 of the DNAPL is infiltrated into the domain from the top. Detailed boundary conditions are specified in Fig. 1 and Table 2. Initial conditions where the domain is fully saturated with water and a hydrostatic pressure distribution is considered (i.e. p0w = (0.65 − y ) · 9810, s0n = 0). The permeability tensor K\lens of the domain \lens is



K\lens =

10−10 −5 × 10−11



−5 × 10−11 m2 10−10

188

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200 Table 2 2d problem boundary conditions.

IN

Jn = −5.137 × 10−5 , Jw = 0

N S E ∪ W

Jn = 0.00, Jw = 0.00 Jw = 0, Jn = 0.00 pw = (0.65 − y ) · 9810, sn = 0

Fig. 2. Evolution of the non wetting saturation sn at times t = 20 0, 40 0, 60 0, and t = 800 (top) and the corresponding adaptive grid structure (bottom).

and the permeability tensor Klens of the lens lens is



Klens =

6 × 10−14 0



0 m2 . 6 × 10−14

The coarsest (macro) mesh consists of 60 quadrilateral elements globally refined everywhere to the finest level would result in 3840 elements. The final time is T = 800 s. For visualization, we later plot the solution of s over the line

x(σ ) = (1 − σ )(0.25, 0.65 )T + σ (0.775, 0.39 )T

(33)

with σ ∈ [0, 1]. Snapshots of the evolution of the resulting flow and the grid structure are shown in Fig. 2. 5.2. Time step stability In this section, we compare the various splitting and solution strategy described in Section 4.2. We compare three implicit and iterative coupling schemes and two loosely coupled schemes, one of them the classical IMPES scheme. Stability analyses of the IMPES scheme are carried in [41–43]. There are usually strong restrictions with respect to the choice of the time step size as shown by the cross-comparison of Todd’s [44], Coats’ [43] and classical Courant–Friedrichs–Lewy’s [45] criteria provided in [46]. Thus, it is quite likely that using one of those conditions would lead to a smaller time step for the IMPES method reducing its competitiveness. Regarding the fully implicit scheme, it is worth mentioning the work of Bastian [6] where the time step size was chosen adaptively by the algorithm depending on the convergence of the Newton method. In order to provide a fair minded comparison of the different methods, i.e., the fully implicit and the partially explicit ones, we use here a straightforward ”trial and error” approach by setting the same time step size for all the methods. For a time step τ > 3, only the fully coupled schemes are able to produces reasonable solutions. This is illustrated in Fig. 3. In Fig. 4 we compare the solution of the various fully coupled schemes for different time step sizes. If the solution converges then a correct solution profile is produced. The stability of the implicit scheme is influenced by the fact that the stabilization operator is only applied before and after the Newton solver. In principle, the explicit coupling schemes work fine for small time steps and fail to produce a valid solution for larger time steps. Here, the implicit schemes show their strength allowing for faster computation once the time step is chosen sufficiently large. 5.3. Different model: Model B The aim of this section is to demonstrate the flexibility of our code, by comparing our original model formulation with a description where the two-phase flow problem is modeled as a system of equations with two unknowns p¯ and sn : p¯ = pw + 12 pc .



−∇ · (λw + λn )K∇ p¯ +

(λn − λw ) 2



p c K∇ sn − (ρw λw + ρn λn )Kg = qw + qn on  × (0, T ),

(34)

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

0.6 0.5

implicit iterative impesIterative impes linear

0.5 0.4 s (x(σ),t=800s)

0.4 s (x(σ),t=800s)

0.6

implicit iterative impesIterative impes linear

0.3 0.2

0.3 0.2

0.1

0.1

0

0

-0.1

189

-0.1 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Fig. 3. All schemes for τ = 3 (left) and τ = 5 (right). All schemes are able to capture the solution characteristics for times steps smaller and up to τ = 3. For a time step τ > 3 only the fully coupled schemes are able to produce reasonable solutions, while the loosely coupled schemes do no longer capture the front position correctly.

0.6

0.6

τ=3 τ=5 τ=7

0.5

0.5 0.4 s (x(σ),t=800s)

0.4 s (x(σ),t=800s)

τ=3 τ=5 τ=7 τ=7 τ=9 τ=13

0.3 0.2

0.3 0.2

0.1

0.1

0

0

-0.1

-0.1 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Fig. 4. Left, the solution for fully coupled implicit scheme for τ = 3, 5, 7. Right the solution for the IMPES-iterative scheme for τ = 3, 5, 7, 9, 11, 13.

φ

  1   ∂ sn − ∇ · λn K(∇ p¯ − ρn g ) − ∇ · λn p c K∇ sn = qn ∂t 2

on  × (0, T ).

The required changes to the original code are shown in Appendix C.1. To complete the system, we add appropriate boundary and initial conditions.

p¯ (·, 0 ) = p0w (· ) + p¯ = pw,D +

1 pc (s0n (· )), 2

1 pc ( sD ), 2

vα · ν = Jα ,

sn (·, 0 ) = s0n (· ),

in ,

sn = sD ,  Jt = Jα ,

on D × (0, T ), on N × (0, T ).

α ∈{w,n}

Here, Jα ∈ R, α ∈ {w, n} is the inflow, s0n , p0w , sD , and pw, D are real numbers. Following the general description of the problem given in Section 3 we have p = p¯ , s = sn and

λn (s ) − λw (s )

p c (s )K,

A pp (s ) = (λn (s ) + λw (s ))K,

A ps (s ) =

Asp (s ) = λn (s )K,

Ass (s ) =

Gs ( s ) = 0,

G p (s ) = −(ρw λw (s ) + ρn λn (s ))Kg,

λn (s ) 2

2 p c (s )K,

Pg = ρn g, q p = qw + qn ,

qs = qn .

The required changes to the Python code are again minimal and described in the Appendix C.1.

(35)

190

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

0.6

0.6

implicit iterative

0.5

0.5 0.4 s (x(σ),t=800s)

s (x(σ),t=800s)

0.4 0.3 0.2

0.3 0.2

0.1

0.1

0

0

-0.1

-0.1 0.2

0.3

0.4

0.5

0.6

0.6

0.7

0.8

0.2

0.3

0.4

0.5

0.6

implicit iterative

0.5

0.6

0.7

0.8

implicit iterative

0.5 0.4 s (x(σ),t=800s)

0.4 s (x(σ),t=800s)

implicit iterative

0.3 0.2

0.3 0.2

0.1

0.1

0

0

-0.1

-0.1 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Fig. 5. Results using modelB. Top row: τ = 1, 3, bottom row: τ = 5, 7.

In the following, we investigate the stability of the different methods with respect to the time step size when applied to modelB. We perform the same investigation described in the previous section where we used modelA. The results are summarized in Fig. 5. We only investigated the stability of the three methods implicit, iterative, and impes-iterative. For modelB the splitting introduced in the impes type approach failed even for τ = 1 while the other two methods produce results in line with the results produced with modelA, although for higher values of τ , the iterative methods produce a discontinuity at the right most front as can be seen in the plots on the bottom row of Fig. 5. For τ > 5, the implicit method fails, making modelB a less stable choice for this scheme. On the other hand, the iterative approach produced results also for larger time steps τ = 9, 11, 13, 15 (not shown here) but in each case the solution showed the same type of discontinuity. Taking all the approaches into account, it is clear that modelA is the more stable representation of the problem. But our results also indicate that the stability of the iterative scheme does not seem to depend so much on the choice of the model (at the least for the two versions tested) and produces very similar results in both cases.

5.4. P-adaptivity Here, we compare different approaches for the indicator used to set the local polynomial degree. In the following, we always use the implicit method with τ = 5, h-adaptivity with a maximum level of three and also a maximum level of three for the polynomial order. In addition to the approach used previously, we test a version without p-adaptivity and an indicator based on determining the smoothness of the solution. In regions where the indicator detects a reduction in smoothness the polynomial order is reduced but only if the grid has been refined to the maximum allowed level. The smoothness indicator is based on ςE = Algorithm 3. The changes to the code are described in Appendix C.2.

ηEr and we set ptol = 1. Details are provided in ηEr−1

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

191

Algorithm 3 p-adapt: markpFrac. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19:

Let ςE =

ηEr be given ηEr−1

∀E ∈ Th

for all E ∈ Th do rE := poldeg(E ) if ςE < 0.01 × ptol then if rE < maxpoldeg then rEnew := rE + 1 else rEnew := rE end if else if ςE > ptol then if rE > 1 then rEnew := rE − 1 else rEnew := rE end if else rEnew := rE end if end for

Fig. 6. Comparison of different approaches for choosing the local polynomial degree. The top row shows the distribution of the polynomial order (left: original approach, right: modified indicator). Red color refers to r = 3 and blue refers to r = 1. The bottom row shows the grid level used for the three simulation. Red refers to l = 3 and blue refers to l = 0. Left to right: original approach, modified indicator, with uniform polynomial degree of three. White lines show 20 contour levels between sn = 0 and sn = 0.55. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 shows the distribution of the polynomial order for the two p-adaptive approaches. The local grid adaptivity is of course also influenced by the choice of indicator for the polynomial degree because the values of the residuals change. This can also be seen in Fig. 6. As expected for the approach discussed in Section 4.3 the polynomial order is reduced to the smallest admissible value (r = 1) in the regions where sn is constant. The approach based on the smoothness indicator described in this section clearly leads to a reduction in the polynomial order at the interface of the plume. This is expected since here the solution can be said to have lower regularity. This demonstrates that the indicator works as expected. On the left in Fig. 7 the solution for the different approaches along the same line given in Eq. (33) is shown. While the solution for the original adaptive method and the solution without adaptivity are indistinguishable, the solution with the smoothness indicator given above shows clearly an increase in numerical diffusion at the interfaces - especially at the entry point to the lens which also leads to an increase within the lens. The number of degrees of freedom during the course of the simulation depends on the number of elements and the local distribution of the polynomial degree used. The number of elements increases in time and so does the number of degrees of freedom. The right plot in Fig. 7 shows the number of degrees of freedom as a function of time. Clearly, the method with maximal polynomial degree on all elements requires the most degrees of freedom. Using the original indicator of p-adaptivity reduces the number of degrees of freedom to about 66% at the beginning of the simulation and still to 75% at the final time. The smoothness indicator given in this section only leads to a reduction of 20% at the beginning of the simulation and by only 6% at the final time.

192

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200 18000

reference smoothness indicator original no p adaptivity

0.6

number of degrees of freedom

0.5

original smoothness indicator no p adaptivity

16000

0.4

0.3

0.2

14000 12000 10000 8000 6000 4000

0.1 2000 0

0 0 0.2

0.3

0.4

0.5

0.6

0.7

100

200

0.8

300

400 time

500

600

700

800

Fig. 7. Comparison of different approaches for choosing the local polynomial degree. Shown sn over the line in Eq. (33).

Fig. 8. Evolution of the non wetting saturation sn for the isotropic weak lens setting at times t = 80 0, 160 0, 240 0, and t = 320 0 (top) and the corresponding adaptive grid structure (bottom).

Overall the indicator described in Section 4.3 does seem to lead to a better distribution of the polynomial degree with negligible influence of the actual solution. As can be seen from Fig. 6, there is only little reduction of the order in the actual plume and the intermediate order p = 2 is hardly used anywhere in the domain. Both these observations indicate that further research into p-adaptivity for this type of problem is required. 5.5. Isotropic flow over a weak lens In the final section we study a second test case. The setup is the same as in the previous example but the permeability tensors is isotropic and the lens is weaker:



K\lens =

10−10 0



0 m2 , 10−10



Klens =

10−12 0



0 m2 . 10−12

On the Python side the problem class the definition of K has to be modified accordingly as shown in Appendix C.3. Snapshots of the evolution of the resulting flow are shown in Fig. 8. The symmetry is clearly visible both in the solution and in the grid refinement. The grid is locally refined at the interface and around the lens once the flow reaches that point. Due to the weaker heterogeneity, flow also passes through the lens. Overall our tests indicate that the conclusions obtained from our previous tests also apply to this setting. 6. Conclusion We have presented a framework that allows to study hp-adaptive schemes for two-phase flow in porous media. The presented approach allows to easily study different adaptive strategies and time stepping algorithms. The change from one algorithm to another is easily implemented. All python based implementation is in the end forwarded to C++ based implementations to ensure performance of the applications.

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

193

Furthermore, the prototypes build on very mature implementations from the Dune community to allow for a short transition from prototypes to production codes. Parallelization and extension to 3d are straightforward with the prototype presented. We focused in this paper mostly on the stability of different approaches for evolving the solution from one time step to the next. Our results indicate that IMPES type splitting schemes can be used when smaller time steps are acceptable while for larger time steps schemes solving the fully coupled equations should be used. The simple fixed point iteration seems to work quite well and seems to be quite stable with respect to the algorithmic approaches and models tested. It turns out that the method does show convergence issues when too small values for the stopping tolerance are used. The Newton method has more difficulty with convergence for large time steps. This suggests a combined approach where first the iterative scheme is used to reach a reasonable starting point for the Newton solver. We tested this approach and could obtain reasonable results up to time steps of size τ = 25 effectively tripling the maximum time steps achievable when using only the Newton scheme. In future work, we will focus more on this approach. In addition, we will investigate 3d examples and the possible extension to polyhedral cells which are widely used in industrial applications. Preliminary work has been carried out, for example, in [47]. The deployment of higher order adaptive schemes is an essential tool for capturing reactive flows for applications such as polymer injections for improved oil recovery or CO2 sequestration. Here, improved numerical algorithms help to reduce uncertainty for predictions and thus ultimately improve decision making capabilities for involved stakeholders.

Acknowledgments Birane Kane acknowledges the German Research Foundation through the Cluster of Excellence in Simulation Technology (EXC310) at the University of Stuttgart for financial support. Robert Klöfkorn acknowledges the Research Council of Norway and the industry partners, ConocoPhillips Skandinavia AS, Aker BP ASA, Eni Norge AS, Maersk Oil; a company by Total, Statoil Petroleum AS, Neptune Energy Norge AS, Lundin Norway AS, Halliburton AS, Schlumberger Norge AS, Wintershall Norge AS, and DEA Norge AS, of The National IOR Centre of Norway for support.

Appendix A. Reproducing the results in a Docker container To easily reproduce the results of this paper we provide a Docker [27] image containing the presented code in a Jupyter notebook [26] and all necessary software to run it. Once Docker is installed, the following shell command will start the Jupyter server within a Docker container:

Notice that all user data will be put into and kept in the Docker volume named dune for later use. This volume should not exist prior to the first run of the above command. Open your favorite web browser and connect to 127.0.0.1:8888 and log in; the password is dune. The notebook twophaseflow contains the code used to obtain the results in this paper. Note that the cells in the notebook depend on each other and thus can not be executed individually. To run the notebook use the Run all option from the Cell menu item.

Appendix B. Main code structure In this section we show parts of the python script used in the simulation. The snippets are not self contained but should provide enough information to understand the overall structure. The full code which can be used to produce the simulations in Section 5.2 is available as a jupyter notebook (see Appendix A) for details. In the following, we describe parts of the code following the overall structure of Sections 3 and 4.

B1. Model A: Wetting-phase-pressure/nonwetting-phase-saturation formulation The model description is decomposed into two parts: this first part consists of a problem class containing a static function for the pressure law pc , the permeability tensor K, boundary, and initial data, and the further constants needed to fully describe the problem:

194

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

The Brooks–Corey pressure law is given by a function taking a problem class as first argument and the value non wetting phase sn :

The actual PDE description requires three vector valued coefficient functions one for the solution on the new time level (u), one for the solution on the previous time level (solution_old), and one for the intermediate state s¯ used in the iterative approaches (intermediate). The vector valued test function is v. Furthermore, τ , β are constants used used for the time step size, the penalty factor, respectively. These can be set dynamically during the simulation:

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

195

B2. Space discretization Given the expressions defined previously the bulk integrals for the bilinear forms can now be easily defined (compare Section 4):

Next, we describe the skeleton terms required for the DG formulation. We use some geometric terms defined for the skeleton of the grid and also the weighted average:

As shown in Section 4 it is straightforward to construct the required penalty and consistency terms

The factors for the penalty terms for the DG discretization depend on the model and are given by

196

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

B3. Time stepping The final bilinear forms used to carry out the time stepping depend on the actual schemes used. We first need to distinguish between the three schemes linear, implicit, iterative that are based on the full coupled system and the two schemes impes, iterative-impes which are based on a decoupling of the pressure and saturation equation. In the first case the final bilinear form is simply

while in the second case we define a pair of scalar forms:

Finally, we need to fix s¯, i.e., intermediate according to the scheme used. In the case of the implicit scheme we have intermediate=u, for linear and impesintermediate=solution_old, while for the other two schemes intermediate is an independent function used during the iteration. The following code demonstrates how the evolution of the solution from ti to t i+1 is carried out:

where the stopping criteria is given by

The implementation of the iterative-impes method looks almost the same

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

197

B4. Stabilization Note how we apply the limiting operator directly after the next iterate has been computed. The stabilization projection operator is available as limit(solution).

B5. Adaptivity The estimator is given as a form taking vector valued solution u with a scalar test function v0. This will later be used to generate an operator taking the solution and mapping into a piece wise constant scalar space with the value ηE on each element:

and since we want to use the estimator for the fully coupled implicit problem independent of the actual time stepping approach used, we add

The actual grid adaptivity is then carried out by calling:

where the marking function is

Finally the p-adaptivity requires calling

where the marking function

markp is given by

198

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

Appendix C. Code modifications C1. Different model: Model B Changing the formulation of the two-phase flow model requires redefining the terms for the bulk integrals and the penalty factor for the DG stabilization. The adaptation indicators and other DG terms do not need to be modified:

C2. P-adaptivity To change the marking strategy for the p-adaptivity the function

markp needs to be redefined:

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

199

C3. Isotropic flow over weak lens Changing the set up of the problem requires modifying the static components of the problem class, i.e., for the isotropic setting with the weaker lens we need to change permeability tensors:

References [1] P. Bastian, Numerical computation of multiphase flow in porous media, habilitationsschrift Univeristät Kiel, 1999. Ph.D. thesis. [2] D.A. Di Pietro, A. Ern, Mathematical Aspects of Discontinuous Galerkin Methods, 69, Springer Science & Business Media, 2011. [3] P. Bastian, R. Helmig, Efficient fully-coupled solution techniques for two-phase flow in porous media: Parallel multigrid solution and large scale computations, Adv. Water Resour. 23 (3) (1999) 199–216, doi:10.1016/S0309-1708(99)0 0 014-7. [4] Y. Epshteyn, B. Rivière, Fully implicit discontinuous finite element methods for two-phase flow., Appl. Numer. Math. 57 (4) (2007) 383–401, doi:10. 1016/j.apnum.20 06.04.0 04. [5] Y. Epshteyn, B. Rivière, Analysis of hp discontinuous Galerkin methods for incompressible two-phase flow., J. Comput. Appl. Math. 225 (2) (2009) 487–509, doi:10.1016/j.cam.2008.08.026. [6] P. Bastian, A fully-coupled discontinuous Galerkin method for two-phase flow in porous media with discontinuous capillary pressure, Comput. Geosci. 18 (5) (2014) 779–796, doi:10.1007/s10596- 014- 9426- y. [7] M. Vohralík, M.F. Wheeler, A posteriori error estimates, stopping criteria, and adaptivity for two-phase flows, Comput. Geosci. 17 (5) (2013) 789–812, doi:10.1007/s10596- 013- 9356- 0. [8] C. Cancès, I.S. Pop, M. Vohralík, An a posteriori error estimate for vertex-centered finite volume discretizations of immiscible incompressible two-phase flow, Math. Comput. 83 (285) (2014) 153–188, doi:10.1090/S0025- 5718- 2013- 02723- 8. [9] R. Becker, R. Rannacher, A feed-back approach to error control in finite element methods: Basic analysis and examples, East-West J. Numer. Math. 4 (1996) 237–264. [10] R. Becker, R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numerica 10 (2001) 1–102. [11] J. Kou, S. Sun, Convergence of discontinuous Galerkin methods for incompressible two-phase flow in heterogeneous media, SIAM J. Numer. Anal. 51 (6) (2013) 3280–3306. [12] B. Kane, Adaptive higher order discontinuous Galerkin methods for porous-media multi-phase flow with strong heterogeneities, 2018. Ph.D. thesis, University of Stuttgart http://dx.doi.org/10.18419/opus-9863. [13] B. Kane, Using dune-fem for adaptive higher order discontinuous Galerkin methods for strongly heterogenous two-phase flow in porous media, Arch. Numer. Softw. 5 (1) (2017), doi:10.11588/ans.2017.1.28068. [14] B. Kane, R. Klöfkorn, A. Dedner, Adaptive discontinuous Galerkin methods for flow in porous media, in: Proceedings of the 12th European conference on numerical mathematics and advanced applications, ENUMATH, Voss, Norway, 2018. In press. [15] B. Kane, R. Klöfkorn, C. Gersbacher, hp–adaptive discontinuous Galerkin methods for porous media flow, in: C. Cancès, P. Omnes (Eds.), Finite Volumes for Complex Applications VIII - Hyperbolic, Elliptic and Parabolic Problems, Springer International Publishing, Cham, 2017, pp. 447–456, doi:10.1007/ 978- 3- 319- 57394- 6_47. [16] F. List, F.A. Radu, A study on iterative methods for solving Richards? equation, Comput. Geosci. 20 (2) (2016) 341–353. [17] F.A. Radu, J.M. Nordbotten, I.S. Pop, K. Kumar, A robust linearization scheme for finite volume based discretizations for simulation of two-phase flow in porous media, J. Comput. Appl. Math. 289 (2015) 134–141. [18] F.A. Radu, K. Kumar, J.M. Nordbotten, I.S. Pop, A robust, mass conservative scheme for two-phase flow in porous media including hölder continuous nonlinearities, IMA J. Numer. Anal. 38 (2) (2017) 884–920. [19] S. Karpinski, I.S. Pop, F.A. Radu, Analysis of a linearization scheme for an interior penalty discontinuous Galerkin method for two-phase flow in porous media with dynamic capillarity effects, Int. J. Numer. Methods Eng. 112 (6) (2017) 553–577. [20] J. Kou, S. Sun, On iterative impes formulation for two phase flow with capillarity in heterogeneous porous media, Int. J. Numer. Anal. Model. Ser. B 1 (1) (2010) 20–40. [21] S. Lacroix, Y. Vassilevski, J. Wheeler, M. Wheeler, Iterative solution methods for modeling multiphase flow in porous media fully implicitly, SIAM J. Sci. Comput. 25 (3) (2003) 905–926. [22] A. Kvashchuk, A robust implicit scheme for two-phase flow in porous media, University of Bergen, 2015 Master thesis. [23] A. Dedner, R. Klöfkorn, M. Nolte, M. Ohlberger, A generic interface for parallel and adaptive scientific computing: abstraction principles and the DUNEFEM module, Computing 90 (3–4) (2010) 165–196, doi:10.10 07/s0 0607- 010- 0110- 3. [24] A. Dedner, M. Nolte, The dune-python module, Arch. Numer. Softw. (2018). In press. [25] M.S. Alnæs, A. Logg, K.B. Ølgaard, M.E. Rognes, G.N. Wells, Unified form language: a domain-specific language for weak formulations of partial differential equations, ACM Trans. Math. Softw. 40 (2) (2014) 9:1–9:37, doi:10.1145/2566630. [26] T. Kluyver, B. Ragan-Kelley, F. Pérez, B. Granger, M. Bussonier, J. Frederic, K. Kelley, J. Hamrick, J. Grout, S. Corlay, P. Ivanov, D. Avila, S. Abdalla, C. Willing, Jupyter notebooks – a publishing format for reproducible computational workflows, in: F. Loizides, B. Schmidt (Eds.), Positioning and Power in Academic Publishing: Players, Agents and Agendas, IOS Press, 2016, pp. 87–90. [27] C. Boettinger, An introduction to docker for reproducible research, with examples from the r environment, SIGOPS Oper. Syst. Rev. 49 (1) (2015) 71–79. [28] M. Alkämper, A. Dedner, R. Klöfkorn, M. Nolte, The DUNE-ALUGrid module., Arch. Numer. Softw. 4 (1) (2016) 1–28, doi:10.11588/ans.2016.1.23252. [29] A. Dedner, S. Girke, R. Klöfkorn, T. Malkmus, The DUNE-FEM-DG module, Arch. Numer. Softw. 5 (1) (2017), doi:10.11588/ans.2017.1.28602. [30] A. Dedner, L. Connelian, M. Nolte, The dune-fempy module, Arch. Numer. Softw. (2018). In press. [31] R. Helmig, et al., Multiphase Flow and Transport Processes in the Subsurface: A Contribution to the Modeling of Hydrosystems., Springer-Verlag, 1997. [32] M.T. Van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J. 44 (5) (1980) 892–898. [33] R.H. Brooks, A.T. Corey, Hydraulic properties of porous media and their relation to drainage design, Trans. ASAE 7 (1) (1964) 26–0028. [34] R. Helmig, R. Huber, Comparison of Galerkin-type discretization techniques for two-phase flow in heterogeneous porous media, Adv. Water Resour. 21 (8) (1998) 697–711. [35] A. Ern, I. Mozolevski, L. Schuh, Discontinuous Galerkin approximation of two-phase flows in heterogeneous porous media with discontinuous capillary pressures, Comput. Methods Appl. Mech. Eng. 199 (23) (2010) 1491–1501. [36] W. Klieber, B. Rivière, Adaptive simulations of two-phase flow by discontinuous Galerkin methods, Comput. Methods Appl. Mech. Eng. 196 (1) (2006) 404–419. [37] M. Ainsworth, R. Rankin, Constant free error bounds for nonuniform order discontinuous Galerkin finite-element approximation on locally refined meshes with hanging nodes, IMA J. Numer. Anal. (2009), doi:10.1093/imanum/drp025.

200

A. Dedner, B. Kane and R. Klöfkorn et al. / Applied Mathematical Modelling 67 (2019) 179–200

[38] S. Sun, M.F. Wheeler, L2 (H1 ) norm a posteriori error estimation for discontinuous Galerkin approximations of reactive transport problems, J. Sci. Comput. 22 (1) (2005) 501–530, doi:10.1007/s10915- 004- 4148- 2. [39] Y. Cheng, F. Li, J. Qiu, L. Xu, Positivity-preserving DG and central DG methods for ideal MHD equations, J. Comput. Phys. 238 (2013) 255–280, doi:10. 1016/j.jcp.2012.12.019. [40] X. Zhang, C.-W. Shu, On positivity-preserving high order discontinuous Galerkin schemes for compressible Euler equations on rectangular meshes, J. Comput. Phys. 229 (23) (2010) 8918–8934, doi:10.1016/j.jcp.2010.08.016. [41] T. Russell, et al., Stability analysis and switching criteria for adaptive implicit methods based on the CFL condition, in: SPE Symposium on Reservoir Simulation, Society of Petroleum Engineers, 1989. [42] K. Coats, et al., Impes stability: the CFL limit, in: SPE Reservoir Simulation Symposium, Society of Petroleum Engineers, 2001. [43] K. Coats, et al., Impes stability: selection of stable timesteps, SPE J. 8 (02) (2003) 181–187. [44] M. Todd, P. O’dell, G. Hirasaki, et al., Methods for increased accuracy in numerical reservoir simulators, Soc. Petroleum Eng. J. 12 (06) (1972) 515–530. [45] R. Courant, K. Friedrichs, H. Lewy, Über die partiellen differenzengleichungen der mathematischen physik, Mathematische annalen 100 (1) (1928) 32–74. [46] J. Franc, P. Horgue, R. Guibert, G. Debenest, Benchmark of different CFL conditions for impes, Comptes Rendus Mécanique 344 (10) (2016) 715–724. [47] R. Klöfkorn, A. Kvashchuk, M. Nolte, Comparison of linear reconstructions for second-order finite volume schemes on polyhedral grids, Comput. Geosci. (2017) 1–11, doi:10.1007/s10596- 017- 9658- 8.