A comparison of numerical methods for conservative and positive advection‐diffusion‐production‐destruction systems

In many simulations it is important to guarantee positivity of the numerical approximations. For standard ODE solvers positivity is usually achieved by restrictions of the time step size, in contrast, MPRK schemes are unconditionally positive, which allows for much faster computations. We compare MPRK schemes and standard ODE solvers applied to an advection‐diffusion‐produ system for which it is essential to avoid negative approximations .


Introduction
Modified Patankar-Runge-Kutta (MPRK) schemes are methods for the solution of positive and conservative productiondestruction systems (PDS), which ensure positivity and conservation of the numerical solution for any positive time step size. Based on an introduction of MPRK schemes in [1], a general definition and order analysis, together with families of MPRK schemes up to third order were recently presented in [2][3][4]. Furthermore, it was numerically shown that MPRK schemes are capable of solving stiff PDS like the Robertson problem in a highly efficient and accurate manner.
In this paper we demonstrate the benefits of the unconditionally positive MPRK schemes, when applied to an advectiondiffusion-production-destruction system (ADPDS), for which too large negative approximations lead to divergence of the numerical method. Thus, unconditional positivity preservation represents a crucial property to obtain reliable results. As a model problem for the mathematical description of viscous flow, we consider the ADPDS with advection coefficient a = 10 −2 and diffusion coefficient d = 10 −6 . The right hand sides of (1) form the positive and conservative PDS introduced in [5], which models the interaction between nutrients (u 1 ), phytoplankton (u 2 ), zooplankton (u 3 ) and detritus (u 4 ) based on Michaelis-Menten kinetics. In numerical simulations it is essential that approximations of u 1 are larger than the negative of the half-saturation constant, i. e. −0.01 < u 1 . If u 1 < −0.01 and u 2 > 0 then u ′ 1 may become negative, which may cause divergence of the method, as is the case for the classical schemes depicted in Figure 1.
To solve the ADPDS (1) we employ a Strang splitting and use an explicit finite-volume scheme for the advection-diffusion part. The focus is on the comparison of standard ODE solvers and MPRK schemes for the solution of the reaction step (2), for which it is indispensable to ensure positivity of the numerical approximations. In the case of standard schemes this requires significant time step restrictions compared to the maximum time step size allowed for the advection-diffusion step. In contrast, when applying MPRK schemes for the solution of (2), no additional time step restrictions are necessary, which considerably reduces computing time, while generating good approximations.
We also like to note, that MPRK schemes can be applied in a direct approach without splitting, as shown in [6,7], where MPRK schemes were derived from the SSP (strong stability preserving) formulation of Runge-Kutta schemes.
To solve the ADPDS (1) numerically, we use a finite volume discretization with an upwind discretization of the advective fluxes and a central discretization of the diffusive fluxes. We employ an equidistant grid with mesh width ∆x and assume periodic boundary conditions for simplicity.
The time integration is carried out by means of a Strang splitting, which separates advection-diffusion from the reactive PDS part (2). The advection-diffusion equation is integrated explicitly by Heun's method, which requires a time step size ∆t ad ≤ ∆x 2 2d + a∆x to ensure stability. For the time integration of the PDS (2) we compare the following second order accurate schemes.
• Heun's method: Explicit two-stage Runge-Kutta scheme , see [2], is defined as • ROS2: Linearly implicit two-stage Rosenbrock scheme, L-stable • SDIRK2: Implicit two-stage SDIRK scheme, L-stable The above methods will be used with fixed time steps. In addition, we also make comparisons with the adaptive solvers

Test configuration
Choosing a mesh width of ∆x = 1 100 for x ∈ [0, 1] the maximum time step size for the advection-diffusion step becomes where the CFL number CFL = 0.9 was chosen for safe-guarding. As mentioned above, negative approximations below the half-saturation rate in the reaction step will immediately result in divergence of the numerical scheme. Time step sizes ∆t r which guarantee positivity were determined numerically and are listed in Table 1. If t n is the current time level and ∆t r < ∆t ad several sub-steps of size ∆t r will be taken in the reaction part to traverse the time interval [t n , t n + ∆t ad ], where the last step size is potentially reduced in order to reach t n + ∆t ad exactly. The overall time span of interest is [0, 50].

Numerical results
First, we show numerical approximations of the PDS (2) on the time interval [0, 5] with initial values u = (8, 2, 1, 4) T and time step size ∆t = 0.02 in Figure 1. The MPRK22(1) scheme captures the solution well, the other three schemes compute defective approximations even for this small step size due to the occurrence of negative values of u 1 . We note that MPRK22(1) generates qualitatively correct approximations even for much larger time steps.  To avoid negative approximations the conditionally positive schemes must use much smaller time steps than the MPRK scheme.
The impact of the time step restrictions of the conditionally positive schemes with respect to computing time, becomes particularly apparent when solving the ADPDS (1) as described in sections 2 and 3 with initial conditions evident from Figure 2 a). Table 1 lists the time step sizes necessary to ensure positivity together with the elapsed wall clock time of the simulations. The MPRK scheme is at least 60 times faster than the other schemes, while generating meaningful approximations as shown in Figure 2.
Next, we compare the MPRK22(1) scheme to four adaptive time stepping schemes of the MATLAB ODE suite. Table 2 shows elapsed wall clock times for different absolute and relative tolerances. For coarse tolerances only ode15s with enabled non-negativity can compute valid approximations, but is still 10 times slower than MPRK22(1) with fixed time stepping. In addition, the positivity of the approximations is achieved by clipping which introduces a conservation error. The other adaptive solvers fail to compute sufficiently large nonnegative approximations for coarse tolerances and take at least 10 times the computing time of MPRK22(1) for fine tolerances.  11.29s 11.38s 11.6s 11.68s 11.6s 10 −1 11.46s 11.68s 11.83s 14.84s 10 −2 12.58s 12.50s 13.26s 10 −3 12.86s 13.59s 10 −4 14.44s www.gamm-proceedings.com