Efficient simulation of unsaturated flow using exponential time integration
Introduction
We consider initial value problems resulting from semidiscrete formulations of the unsaturated Richards equation. These problems can be expressed in the following form (see Section 3):where and the number of nodes used in the discretisation determines the value of N.
The class of methods known as exponential integrators [17] solve (1) via the application of a function of the Jacobian matrix ∂g/∂u—either the matrix exponential or one of the closely related “φ functions” [1]. In recent times, these integrators have found application in the numerical integration of stiff problems; see for example the work by Hochbruck et al. [6], Minchev and Wright [17], Schulze et al. [22], Schmelzer and Trefethen [21] and the PhD thesis of Minchev [16]. A particular attraction of these methods is that they can perform very well without the need for preconditioning [17].
In this work we consider the simplest exponential integrator for (1), which we now describe. Employing a time-stepping strategy with un = u(tn) and setting τn = tn+1 − tn, the formulaprovides an approximate means for advancing the solution to (1) in time; where denotes the Jacobian matrix ∂g/∂u evaluated at u = un, gn = g(un) and . According to Minchev and Wright [17] the earliest reference to (2) seems to be the paper by Pope [18] published in 1960. Since then the scheme has been reinvented several times and published under various names, including the exponentially fitted Euler method [6] and the local linearisation method [7], [8], [13]. In our work we refer to the scheme as the exponential Euler method (EEM), in line with the review work of exponential integrators by Minchev and Wright [17].
Our aim is to compare the performance of EEM against the current state of the art for advancing the unsaturated Richards equation in time: variable-stepsize backward differentiation formulae (BDF) coupled with a Jacobian-free Newton–Krylov solver. These include both the backward Euler method [2], [3], [11], [14] and higher order BDF methods [9], [15], [24], [25]. Inherent in these numerical strategies is the solution of a nonlinear system at each time step. This nonlinear system is solved using Jacobian-free Newton–Krylov type methods [10] for the inner iteration, which rely heavily on the use of a good preconditioner to accelerate the convergence of the Newton (outer) iterations [12].
The remainder of this paper is organised as follows. In Section 2 we describe the Richards equation model for simulating unsaturated flow through heterogeneous porous media. In Section 3 we outline the use of the finite volume method to spatially discretise the model and thereby generate an initial value problem of the form (1). A brief overview of BDF methods is given in Section 4, by way of setting up the comparison with EEM. In Section 5 we describe our implementation of EEM in detail, and prove that despite the use of Krylov approximations to φ(τnJn)gn in (2), the scheme retains second order accuracy. A novel two-step scheme that allows efficient estimation of the local error is then proposed. The results of some numerical experiments comparing BDF and EEM are presented in Section 6, where it is shown that our implementation of EEM outperforms variable-order, variable-stepsize implementations of the backward differentiation formulae up to order 2 in terms of computational cost, and is competitive up to order 5 for large integrator tolerances. The main conclusions of the work are summarised in Section 7.
Section snippets
Richards equation for unsaturated flow
Pressure-driven flow in porous media can be modelled using Darcy’s law:where q = (qx, qz)T is the Darcy flux vector, h is pressure head, K is the hydraulic conductivity and ez is the unit vector in the vertical direction, oriented upwards.
Assuming incompressibility, conservation of mass requires thatwhere θ is the volumetric moisture content. Together, (3), (4) give the well-known Richards equation:To complete the model, a functional relationship
Spatial discretisation
We discretise in space using the finite volume method over a two-dimensional, rectangular mesh, such that an arbitrary control volume Vp with node situated at coordinate (xi, zj) has volume ΔVp = Δxi × Δzj. Integrating (4) over Vp, the following semidiscrete form is obtainedwhere indices on the components of the flux vector q denote the representative point of approximation on the control volume face. For example, recalling (3), we
Time integration with BDF
The backward differentiation formulae employed in this work are the variable-order, variable-stepsize BDFs, as implemented in the CVODE module of the Suite of Nonlinear and Differential/Algebraic Equation Solvers (SUNDIALS) [5]. A BDF of order q is an approximate relation between du/dt at t = tn+1, denoted , and the values ui for i = n − q + 1, … , n + 1, taking the form [5]:with stepsize τn = tn+1 − tn. The coefficients αn,i are uniquely determined by the order of the method and
Time integration with EEM
The exponential Euler method (2) is derived by linearising g(u) in (1) about t = tn. At each step this results in the linear initial value problem:to advance the solution from t = tn to t = tn+1. The exact solution of this problem at t = tn+1 determines the scheme (2), restated here for convenience:whereWe note that the method when applied to y′ = λy produces the exact solution . Hence EEM is comparable to
Test problems
We compare the proposed EEM integrator with a variable-order BDF integrator for two test problems. The first problem simulates the infiltration process into a region consisting of four different soils at dry initial conditions. This problem has featured in the work by Diersch and Perrochet [2] and initially Forsyth et al. [3]. Fig. 1 (left) presents a schematic view of the two-dimensional cross-section that forms the simulation domain. A no flux boundary condition is active on all boundaries,
Conclusion
In this paper we provided a practical variable-stepsize implementation of the exponential Euler method (EEM) for advancing semidiscrete solutions of the unsaturated Richards equation in time. In particular, we introduced a new second-order variant of the scheme that enables the local error to be estimated at the cost of a single additional function evaluation. Numerical experiments performed on two challenging test problems demonstrate that our implementation of EEM requires fewer time steps
References (26)
- et al.
On the primary variable switching technique for simulating unsaturated-saturated flows
Advances in Water Resources
(1999) - et al.
Robust numerical methods for saturated-unsaturated flow with dry initial conditions in heterogeneous media
Advances in Water Resources
(1995) - et al.
Dynamic properties of the local linearization method for initial-value problems
Applied Mathematics and Computation
(2002) - et al.
Rate of convergence of local linearization schemes for initial-value problems
Applied Mathematics and Computation
(2005) - et al.
Higher order time integration methods for two-phase flow
Advances in Water Resources
(2002) - et al.
Jacobian-free Newton-Krylov methods: A survey of approaches and applications
Journal of Computational Physics
(2004) - et al.
A higher order local linearization method for solving ordinary differential equations
Applied Mathematics and Computation
(2007) - et al.
A spatially and temporally adaptive solution of Richards’ equation
Advances in Water Resources
(2006) - et al.
Accurate and economical solution of the pressure-head form of Richards’ equation by the method of lines
Advances in Water Resources
(1997) - et al.
EXPINT–A MATLAB package for exponential integrators
ACM Transactions on Mathematical Software
(2007)
SUNDIALS: suite of nonlinear and differential/algebraic equation solvers
ACM Transactions on Mathematical Software
Exponential integrators for large systems of differential equations
SIAM Journal on Scientific Computing
Cited by (19)
A 3-variable PDE model for predicting fungal growth derived from microscopic mechanisms
2019, Journal of Theoretical BiologyCitation Excerpt :The exponential integrators method was first presented in 1960 by Certaine (Certaine, 1960) and has been well-known since late of 1990s by works of M. Hochbruck and others (Hochbruck et al., 1998; Hochbruck and Ostermann, 2010). Thanks to the approximation of the Krylov subspace to the matrix-function vector product φ(dtnJg(un))g(un), the exponential integrators method became very useful for large system of stiff equations (Carr et al., 2011; 2012; Michels et al., 2014). In this article, our work was inspired from the work of E. Carr, I.Turner and P. Perré (Carr et al., 2012) using the “variable-stepsize exponential Euler method” to solve the system of non-linear transport equations with great success.
New efficient substepping methods for exponential timestepping
2017, Applied Mathematics and ComputationCitation Excerpt :Working towards such an algorithm would be an interesting avenue of future research. We have extended the notion of recycling a Krylov subspace for increased accuracy in the sense of [26]. We have applied this new method to the first order ETD1 scheme and examined the effect of taking an arbitrary number of substeps (the parameter S).
The Meshfree Finite Volume Method with application to multi-phase porous media models
2017, Journal of Computational PhysicsThe extended distributed microstructure model for gradient-driven transport: A two-scale model for bypassing effective parameters
2016, Journal of Computational PhysicsEfficient simulation of geothermal processes in heterogeneous porous media based on the exponential Rosenbrock-Euler and Rosenbrock-type methods
2013, Advances in Water ResourcesCitation Excerpt :These approaches have been known since the 1960s, but have only in the last decade become a robust alternative in applications, due to recent developments resulting in improved efficiency in the computation of the involved matrix exponentials. Exponential integrators have been suggested as efficient and robust alternatives for the temporal discretization for several applications (see [8,17,41,45]). Rosenbrock-type methods have been intensively developed in the literature and used in a variety of settings (see [2,20,28,29] and reference therein).
Local Linearization-Runge-Kutta methods: A class of A-stable explicit integrators for dynamical systems
2013, Mathematical and Computer Modelling