Second order fully semi-Lagrangian discretizations of advection-diffusion-reaction systems

We propose a second order, fully semi-Lagrangian method for the numerical solution of systems of advection-diffusion-reaction equations, which employs a semi-Lagrangian approach to approximate in time both the advective and the diffusive terms. Standard interpolation procedures are used for the space discretization on structured and unstructured meshes. The proposed method allows for large time steps, while avoiding the solution of large linear systems, which would be required by an implicit time discretization technique. Numerical experiments demonstrate the effectiveness of the proposed approach and its superior efficiency with respect to more conventional explicit and implicit time discretizations


Introduction
Large systems of advection-diffusion-reaction equations are responsible for most of the computational cost of typical environmental fluid dynamics models, such as are applied in climate modelling, water and air quality models and oceanic biogeochemistry [5], [12], [13]. Also in short and medium range weather forecasting, the number of interacting transported species is significant and drives the choice of optimal time discretization approaches towards methods that allow the use of large time steps [35]. Due to the potentially very large number of equations of this kind that have to be solved simultaneously in order to achieve a complete description of the relevant physical processes, even minor efficiency gains in the discretization of this very classical problem are of paramount practical importance. The standard ways to achieve such optimal efficiency are either the use of implicit schemes, or the application of semi-Lagrangian (SL) techniques, [14], [30] to the advection step, coupled to implicit methods for the diffusion and reaction step. As discussed in [12], [13], SL methods have the advantage that all the computational work that makes them computationally more expensive per time step than standard Eulerian techniques is indeed independent of the number of tracers, which allows to achieve easily a superior efficiency level in the limit of a large number of tracers.
In the recent papers [3], [4], a fully SL approach to both the advection and diffusion step was pursued, which combines the standard SL treatment of advection with SL-like schemes for diffusion, proposed, among others, in [7], [16], [22], [23], [24]. In particular, it was shown in [4] that, even for a single advection-diffusion equation, the fully SL approach can be more efficient than standard implicit techniques.
In the present work, we present a number of improvements to the fully SL approach of [3], [4]. In particular, we show how second order accuracy in time can be achieved. An improved treatment of Dirichlet boundary conditions is also discussed and analysed. The resulting approach yields a very efficient combination, which is validated on a number of classical benchmarks, both on structured and unstructured meshes. Numerical results show that the method yields good quantitative agreement with reference numerical solutions, while being superior in efficiency to standard implicit approaches and to approaches in which the SL method is only used for the advection term.
The outline of the paper is the following. In Section 2, the considered model problems are introduced. Section 3 describes the SL advection-diffusion solver. A stability and convergence analysis of the method is outlined in Section 4. The possible approaches to the treatment of boundary conditions are discussed in Section 5. A numerical validation of the proposed approach on both structured and unstructured meshes is presented in Section 6, while some conclusions and perspectives for future developments are outlined in Section 7.

The model problem
We consider as a model problem the advection-diffusion-reaction equation with Dirichlet boundary conditions c (x, 0) = c 0 (x) x ∈ Ω.
The unknown c : Ω × [0, T ] → R can be interpreted as the concentration of a chemical species that is transported through the domain Ω by the advection and diffusion processes, while undergoing locally a nonlinear evolution determined by the source term f (c). T is the final time, Ω ⊂ R 2 is an open bounded domain, u : Ω × [0, T ] → R 2 is a velocity field and b : ∂Ω × [0, T ] → R denotes the boundary value of the species c. While the proposed numerical method will be presented in this simpler case, the target for more realistic applications are systems of coupled advection-diffusion-reaction equations of the form x ∈ Ω, k = 1, . . . , S.
Here, D denotes a symmetric and positive semi-definite diffusivity tensor, possibly dependent on space and time. As remarked in Section 1, systems of this kind, with a possibly large number of species S, are responsible for the largest share of the computational cost of typical environmental fluid dynamics models, so that even minor increases in the efficiency of the discretization for this very classical problem are of great practical relevance.

Fully semi-Lagrangian methods
It can be observed that the evolution equation (1) is of the form where L denotes a linear differential operator. Consider then the homogeneous equationc associated to (1) with the same initial datumc(x, 0) = c 0 (x) and boundary conditions as in (1). If the evolution operator determining the solution of (3) is denoted by E t , so thatc(t) = E t [c 0 ], by formal application of the variation of constants formula, the solution of (1) can then be represented as If discrete time levels t n , n = 0, . . . , N are introduced, so that t n = n∆t and ∆t = T /N, a numerical method for the solution of equation (1) on the interval [t n , t n+1 ] can then be derived from the representation formula Discretizing the time integral by the trapezoidal rule, one obtains If the diffusion term is dropped in equation (1) and the evolution operator is approximated by a numerical approximation of the flow streamline and an interpolation at the departure point of the streamline, a numerical method based on formula (6) can be interpreted as a semi-Lagrangian extension of the trapezoidal rule with global truncation error of order 2. Semi-Lagrangian methods based on this formula have been used successfully in a large number of applications, see, among many others, [1], [10], [11], [31], [32], [34]. Often, a first order, off-centered version of the above formula is employed, defined for θ ∈ [1/2, 1] by The goal of this work is to extend this approach to the case of a semi-Lagrangian treatment of both advection and diffusion terms. In order to do this, we introduce a spatial discretization mesh G ∆x = {x i , x i ∈ Ω}, where ∆x denotes a measure of the mesh resolution. Notice that the only necessary restriction on the nature of the mesh is that it should be possible to define a polynomial interpolation operator I p [c](·) of degree p, constructed on the values of a grid function c defined on G ∆x . Assuming in the following that such operator is defined on structured and unstructured meshes (for the definition of such interpolation operator see for instance [27]), examples of applications using both meshes will be considered. We then denote the discrete approximations of the solution of (1) at the space-time mesh nodes by c n i and by some numerical approximations of the streamlines starting at x i and defined by the velocity field u over the interval [t n , t n+1 ]. In particular, in [3], [4] explicit Euler or Heun methods were employed to compute these approximations, coupled to a substepping approach along the lines of [9], [29]. More specifically, for the Euler method, given a positive integer m, a time substep was defined as ∆τ = ∆t/m and, for q = 0, . . . , m − 1, was computed, so that z n+1 i = y m . Following the approach outlined in [3] and applied in [4], a first order in time approximation of the solution of (1) can then be defined as where c n = (c n ) i . In (10), δ k is defined, for k = 1, . . . , 4, as δ k = ±δe j (e j (j = 1, 2) denoting the canonical basis in R 2 ), for all combination of both the sign and the index j = 1, 2; moreover, δ = √ 4ν∆t. Notice that, for simplicity, we neglect in (10) the treatment of boundary conditions. Possible approaches to handle nontrivial boundary conditions will be discussed in Section 5.
The method (10) will be denoted in what follows by SL1. This method inherits the same stability and convergence properties of the parent methods, as it will be discussed in Section 4. Notice that this approach can be extended to spatially varying diffusion coefficients and that, while only first order in time, its effective accuracy can be substantially superior to that of more standard techniques, if higher degree interpolation operators are used, as shown in [3].
In order to derive a method of second order in time, we follow the main steps of [16], [25]. The streamlines are interpreted as generalized characteristics, defined as the solutions of the stochastic differential equation: where W (s) denotes a standard 2-dimensional Wiener process. Applying the implicit weak method of order 2 defined in [19] for the approximation of the streamlines, we denote by the solutions of the nonlinear equations Here, e k denote the vectors: These vectors represent a realization of a vector random variable whose distribution is given by the probabilities It is to be remarked that method (12) can be rewritten in terms of the displacements δ n+1 thus yielding an implicit method that is a natural extension to stochastic differential equations of that introduced in [28] and commonly used in meteorological applications for the computation of streamlines in SL methods. A second order in time SL (SL2) scheme can then be defined by a Crank-Nicolson approach as Notice that, with respect to the simpler first order in time variant (10), nine interpolations at the foot of the streamlines must be computed, which clearly makes this approach substantially more expensive than the simpler, first order in time variant. In applications to systems of the form (2), the computational cost of scheme (14) can be marginally reduced by settingc and defining so as to reduce the number of the evaluations of a possibly costly nonlinear term. Furthermore, when the coupling of the diffusion and advection term is weak, it should be possible to decouple again the approximation of a single deterministic streamline from that of the diffusive displacements, which could be added at the end of each approximate streamline without increasing too much the error. The deterministic streamline can then be computed by a substepping procedure, applied either to explicit methods as in (9) or to implicit methods as in (12). For example, a decoupled substepping variant of (12) might be obtained by computing, for q = 0, . . . , m − 1, and setting z k = y m + √ 3∆tσe k . We will denote this decoupled variant with substepping by SL2s.
Notice that, in realistic problems, a bottleneck of the scheme (14) is the fact that the Crank-Nicolson method, while A-stable, is not L-stable, see e.g. [20]. Therefore, no damping is introduced by the method for very large values of the time step and spurious oscillations may arise, see also the discussion in [6]. In order to reduce the computational cost and to address the L-stability issue, different variants of the (14) scheme could also be introduced and compared, along the lines proposed in [33] for the pure advection case. However, this development goes beyond the scope of this paper and will not be pursued here.
Finally, even though achieving full second order consistency is quite complicated in the variable diffusion coefficient case, the previously introduced schemes can be nonetheless extended at least in the simpler configurations as suggested in [3] for the first order case, even though full second order accuracy is not guaranteed any more.

Convergence analysis
We present in this section a convergence analysis for the scheme (14). For simplicity, we assume a one-dimensional problem defined on R × [0, T ], with a time-independent drift term u: The multidimensional case, as well as the time dependence of u, require only small technical adaptations. On the other hand, the convergence analysis on bounded domains is still an open problem for high-order SL schemes, therefore we will not address this problem here. First, for i ∈ Z and n = 0, . . . , N − 1, we rewrite the scheme (14) with the shorthand notation where In one space dimension, the three discrete characteristics are defined by the equations with corresponding weights α + = α − = 1/6 and α 0 = 2/3. In what follows, we will use the symbol K to denote various positive constants, which do not depend on ∆t, x, t. We will also assume that: (H0) there exists a unique classical solution of (17);

Consistency
First, we derive a consistency result via a Taylor expansion. The same kind of result can be obtained by probabilistic arguments, see [23].
Proposition 1. Assume (H1)-(H3), and let c(x, t) be a smooth solution of (17). Then, for each (i, n) ∈ Z×{0, . . . , N −1} the consistency error of the scheme (14), defined as Proof. In what follows, we will omit the argument of functions computed at (x, t). Consider a smooth solution c of (17). Since assumption (H1) holds, by differentiating in time and space (17) we get that c is also solution of and hence, by differentiating again in space (20), of Using (20) and (21) in (19), we get : Define now U ± (x) = u (x) + u (z ± (x)). By a Taylor expansion of and, defining U 0 (x) = u (x) + u (z 0 (x)), Using (22), (23), (24) and the Taylor expansion Consider now the nonlinear reaction term. By assumption (H3), we have that Using (25) in (26), and taking into account that c(z By (27) and (25), we get the consistency error for the semi-discretization, Introducing the interpolation error and using assumptions (H2) and (H1), we finally prove the consistency error for the fully discrete scheme.

Stability
To prove stability, it is convenient to recast (18) in the matrix form where, for a vector c, f (c) denotes the vector obtained by applying f elementwise and the matrices B k (which represent the operation of interpolating c n at the points z k (x i )) have elements b k,ij defined by for a suitable basis of cardinal functions {ψ j }. The following proposition implies stability for the linear part of the scheme with respect to the 2-norm. For a formal definition of the basis functions ψ j we refer the reader to [17], [18].
Proposition 2. Assume (H2), and let the matrix B have elements defined by (29), with (ψ j ) basis functions for odd degree symmetric Lagrange or splines interpolation. Then, for each k, there exists a constant K B > 0 independent on ∆x, ∆t such that Proof. If we adopt a piecewise linear interpolation, the scheme is monotone and satisfies (30) both in the ∞-norm (with K B = 0) and in the 2-norm, see for example [17]. We will therefore focus on the case of high-order interpolations, for which the norm in (30) should be understood as 2-norm. Following [17], [18], we sketch the arguments to prove (30) for the cases of symmetric Lagrange and splines interpolation, which can be interpreted as Lagrange-Galerkin schemes with area-weighting. First, we make explicit the dependence of the points z k on x and ∆t. We recall that z k (x) = x + δ k (x), with δ k solving the equation: Expanding the term u(x + δ k (x)), we obtain therefore and hence, Dividing now by 1 + ∆t u /2, and using the fact that δ k = O( √ ∆t), we get, for ∆t → 0, Due to the term √ 3∆tσe k , the form of (31) does not coincide with that used in [17] for the points z k (x). However, when considering differences z k (x 1 ) − z k (x 2 ), this additional term is cancelled, so that As a consequence, the form (31) still satisfies the relevant properties used in the proof of (30). In particular, it satisfies the condition [17, for ∆t small enough, as well as the condition [17, Theorem 4] for a suitable positive constant K X . Then, a careful replica of the arguments used in [17] provides the estimate (30).

Convergence
We now present a convergence result in the discrete 2-norm.
Theorem 3. Assume (H0)-(H3). Let c(x, t) be the classical solution of (17), and c n be the solution of (18). Then, for any n such that t n ∈ [0, T ] and for (∆t, ∆x) → 0, where K T is positive constant depending on the final time T .
Proof. While a mere convergence proof could be carried out with weaker regularity assumptions, we will focus here on the error estimate above, which require the regularity assumed in (H0)-(H3).
Define the vectors γ n and ν n , so that γ n i = c(x i , t n ), and n = γ n − c n . Then, by Proposition 1, we get Subtracting (18) from (32), using the Lipschitz continuity of f and the triangle inequality in the form of a difference, we obtain from the left-hand side: Taking into account that k α k = 1, along with the bound (30), we also have from the right-hand side: Therefore, it turns out that Now, for ∆t small enough to have 1 − L f ∆t/2 > C > 0, we have that there exists a constant K T > 0 such that and hence, using this bound in (33), which, by standard arguments, implies that, for any n such that t n ∈ [0, T ], n 2 ≤ K T ∆t 2 + ∆x p ∆t .

Boundary conditions
The treatment of Dirichlet boundary conditions for this class of semi-Lagrangian methods has been considered in [24], where two methods are proposed. One approach has first order of consistency, but it does not seem possible to generalize it to multiple dimensions. The second approach has order of consistency 1/2. More recently, in [4], an easier treatment has been proposed for the scheme SL1 with timeindependent Dirichlet boundary condition, again with order of consistency 1/2. This approach has been extended in [2] to unstructured meshes. We propose here a new approach to obtain second order consistency for the scheme SL2 with Dirichlet boundary conditions. This technique is based on extrapolation, much in the spirit of the so-called ghost-point techniques, see e.g. [21].
In addition to the standard mesh G ∆x = {x i , x i ∈ Ω}, where the numerical solution is computed, we consider a second mesh G h = {v i , v i ∈ Ω} formed by a single layer of elements having their external side along the boundary of Ω. This second mesh is constructed with a size parameter h ∼ √ ∆t, and the degrees of freedom are chosen in order to allow a second-order interpolation. In Fig.1 we show, as an example, a square domain Ω = (−1, 1) × (−1, 1), for which the standard mesh G ∆x is formed by the blue triangular elements and the mesh G h is formed by the black rectangular elements. Note that the latter overlap at the corners. The asterisks in red are the nodes of G h , according to the standard Q 2 element. The values of the numerical solution on the nodes v i are obtained by interpolation at internal nodes, and by the Dirichlet boundary condition if the nodes lie on the boundary ∂Ω.
Let us call T ∆x a given triangulation, with G ∆x the set of the vertexes of the elements K ∈ T ∆x . Let us define the polygonal domain Ω ∆x := ∪ K∈τ ∆x K ⊂ Ω. If, for some i and k, z n+1 k,i / ∈ Ω ∆x , then its projection P (z n+1 k,i ) onto Ω ∆x is computed, defined as the point in Ω ∆x Figure 1: Unstructured computational mesh (blue triangular elements) together with boundary mesh G h (black rectangular elements and red asterisks nodes) at minimum distance from z n+1 k,i . The value of the numerical solution c n (z n+1 k,i ) is then approximated by a quadratic extrapolation operator Ψ 2 . This operator is constructed via the Q 2 associated to the element of G h to which the projection P (z n+1 k,i ) belongs: whereĉ n corresponds tô The method can be extended to more general domains, by considering triangular elements for G h . This technique will not be rigorously analysed here, but we will provide a numerical validation in Section 6.5.

Numerical results
A number of numerical experiments have been carried out, in order to assess the accuracy of the proposed methods on both structured and unstructured meshes. In the unstructured case, we have constructed a triangular mesh by the Matlab2019 function generateMesh, with a maximum mesh edge of ∆x, and used a P 2 space reconstruction. In the structured Cartesian case, the bicubic polynomial interpolation implemented in the Matlab2019 command interp2, has been used. Since the goal is to evaluate the accuracy of time discretization, both choices avoid to hide the time discretization error with the error introduced by a lower order space reconstruction.

Pure diffusion
In a first, basic test, we consider equation (1) in the pure diffusion case, i.e., with zero advection and reaction terms, on the square domain Ω = (−2, 2) × (−2, 2), with T = 1 and ν = 0.05. Based on the test case proposed in [26], we assume a Gaussian initial datum centered in (0, 0), with σ = 0.1, so that the exact solution in an infinite plane would be .
For this test case, we only consider structured meshes with constant steps ∆x = ∆y = 4/N in both directions. Following [3], we consider different time step values ∆t, which correspond to different values of the parabolic stability parameter µ = ∆tν/∆x 2 . We compare method SL1 (10) and method SL2 (14), and collect the results in Table 1. It can be observed that the expected convergence rates are recovered. Furthermore, it is apparent that scheme SL2 yields a substantial accuracy improvement, without an excessive increase in computational cost. Indeed, the SL2 runs require between 30% and 60% more CPU time, depending on the resolution, while leading to corresponding error reductions between 140% and 730%. As a comparison, a standard second order discretization in space coupled to an explicit second order method in time yields at the finest resolution an error 5 times larger than that of method SL2 at approximately the same computational cost.
It can be observed that the expected convergence rates with respect to the time discretization error are recovered, in the constant ∆x, constant C or constant µ convergence studies. It can also be observed that the decoupled variant SL2s, in spite of the loss of second order convergence, does indeed improve the results with respect to the SL1 method and is competitive with the full second order method SL2. As a comparison, a standard centered finite difference, second order discretization in space coupled to an explicit second order method in time yields at the finest resolution an error analogous to that of method SL2 but requires approximately three times its CPU time.
In the unstructured case, the quadratic polynomial interpolation naturally associated to P 2 finite elements was employed and only the SL2s and SL2 methods were considered. The triangular mesh used was chosen with maximum triangle size ∆x approximately equal to the corresponding structured meshes. The results are reported in Table 3.
While the behaviour of the SL2 scheme is entirely analogous to that of the structured mesh case, the SL2s method shows in this case little error reduction when the spatial resolution is kept fixed.

Reaction-diffusion equations
Following [15], we consider the Allen-Cahn equation on the domain Ω = (0, 1) × (0, 1), with periodic boundary conditions and for t ∈ [0, 2]. As in [15], we take the initial datum c 0 (x, y) = sin (2πx) sin (2πy) and a reference solution is computed by a pseudo-spectral Fourier discretization in space, see e.g. [8], and a fourth order Runge-Kutta scheme in time with a very large number of time steps. The results

Resolution
Relative error Convergence rates ∆x λ µ l 2 l ∞ p 2 p ∞ 0.04 16   are reported in Table 4, for the values ν = 0.01 and ν = 0.05 of the diffusion parameter, respectively. In this case, only unstructured meshes were considered and the reference solution was interpolated onto the unstructured mesh nodes using a higher order interpolation procedure. Both tests show a quadratic order of convergence.
which represent two coupled Lotka-Volterra prey-predator systems. As initial datum for c 1 , c 3 , the function was considered, while the initial datum for c 2 , c 4 , was taken to be equal to 3c 0 . In this test, only a structured mesh was considered with constant steps ∆x = ∆y = 1/20. A reference solution is computed by a pseudo-spectral Fourier discretization in space and a fourth order Runge-Kutta scheme in time, using a very large number of time steps. The reference solution is reported for two sample components in Figure 2, while the absolute error distributions obtained for the same components with the second order method SL2 (14) using cubic interpolation, using a timestep corresponding to λ ≈ 7 and µ ≈ 1/2, are shown in Figure 3. As a reference, the errors for a second order finite difference approximation of (36) using a second order Runge-Kutta scheme in time with a time step 20 times smaller are shown in Figure 4, while the errors obtained using a fourth order finite difference approximation for the advection term in (36) with a third order Runge-Kutta scheme in time are displayed in Figure 5, again computed with a time step 20 times smaller than that used for the SL2 method. It can be seen that the SL2 method allows to achieve errors of the same order of magnitude as those of the third order Runge-Kutta in time, while allowing for a much larger time step without solving large algebraic systems.

Advection-diffusion equation, nonhomogeneous boundary conditions
In this last set of numerical experiments, we consider nonhomogeneous, possibly time-dependent Dirichlet boundary conditions in four cases: pure diffusion, constant advection-diffusion, solid body rotation with diffusion and advection-diffusion on a nonconvex domain. In all these tests, we have used an unstructured mesh. In the first three cases, we consider Ω = (−1, 1) × (−1, 1), final time T = 1 and an initial condition in the form of a Gaussian centered at (x 0 , y 0 ) = (0.5, 0), with σ = 0.1. In Fig.1, we show the space meshes G ∆x and G h corresponding to the steps ∆x = 0.04, h = 0.5, which were used to compute the results in the first two rows of Tables 5-7. In order to have a reference solution to compare with, we compute the exact solution on the whole of R 2 and enforce its values at the boundary as boundary conditions, so that b(x, y, t) = c(x, y, t) for (x, y) ∈ ∂Ω, t ∈ [0, T ]. For all the three cases, we have set ν = 0.05 and T = 1. In the second and third test, the advection field has been chosen as u = (1, 0), and u = (−2πy, 2πx), respectively. Tables 5-7
The velocity field u (x, y) is given by where we set u 0 = 0.2 and r 2 = (x − x 0 ) 2 + (y − y 0 ) 2 . In Fig. 6, we show the domain Ω, discretized using a Delaunay mesh G ∆x with ∆x = 0.1, refined around the circular hole. In Fig. 7, we show the numerical solution computed with SL2 with time step ∆t = 0.005 for time t = 0.5, 1, 2, 3. The nonhomogeneus boundary condition are computed by extrapolation with an extra grid G h with h = 1.5 √ ∆t. In this case, the additional mesh G h has been built around the circular Figure 6: Unstructured mesh for the non-convex problem hole, as well as along the external rectangular boundary. Note that, in spite of the discontinuity of the initial configuration and the sharp boundary layer around the hole, the boundary condition is smoothly propagated in the interior of the domain.

Conclusions
A family of fully semi-Lagrangian approaches for the discretization of advection-diffusion-reaction systems has been proposed, which extend the methods outlined in [3], [4] to full second order accuracy. The stability and convergence of the basic second order method has been analyzed. The proposed methods have been validated on a number of classical benchmarks, on both structured and unstructured meshes. Numerical results show that these methods yield good quantitative agreement with reference numerical solutions, while being superior in efficiency to standard implicit approaches and to approaches in which the SL method is only used for the advection term. In future developments, the proposed method will be extended to higher order discontinuous finite element discretizations along the lines of [33] and will be applied to the development of second order fully semi-Lagrangian methods for the Navier-Stokes equations along the lines of [2], [4].