In this work we develop numerical methods for the scalar nonlinear convection-diffusion equation
with . We focus on methods that work well for regimes ranging from the hyperbolic setting () to the diffusion-dominated setting.
A great deal of research has been devoted to the development of numerical methods for hyperbolic conservation laws that are accurate and preserve qualitative solution properties, such as guaranteeing a priori bounds on the solution or avoiding spurious oscillations. Among the most useful schemes available are high-order finite volume methods, which tend to exhibit much less oscillation than generic schemes, but can still produce small over- or undershoots near discontinuities. For some applications, these violations of the bounds are unacceptable and more robust schemes are thus required.
Certain broad theoretical results establish the difficulty of guaranteeing properties like positivity or a maximum principle. These often involve a tradeoff between accuracy and robustness. Godunov’s theorem dictates that any linear PDE discretization that does not generate new local extrema must either be nonlinear or only first-order accurate. Similarly, landmark results by Bolley & Crouzeix bolley1978 and by Spijker spijker1983 establish that any general linear method for initial value ODEs that is guaranteed to satisfy a positivity or monotonicity property can be only first-order accurate.
High-order bound-preserving discretizations
High-order discretizations that satisfy positivity, monotonicity, or total variation diminishing (TVD) properties must therefore fall outside the usual straightforward discretizations. In the context of hyperbolic PDEs, such methods have constituted a major area of research for several decades. Well established techniques include second-order methods based on TVD limiters, as well as flux corrected transport (FCT) methods. Each of these imposes a local bound on solution updates, but such methods are at most second order accurate osher1984high ; zhang2010maximum , due to being too restrictive around smooth extrema.
Some approaches, like weighted essentially non-oscillatory (WENO) methods, give up on providing strict preservation of the maximum principle in favor of achieving better than second-order accuracy. Within the context of finite elements, high-order Bernstein polynomials can be employed with FCT-like methods to impose local bounds; see for example anderson2017high ; hajduk2021monolithic ; lohmann2017flux . To recover high-order accuracy in lohmann2017flux , the authors use smoothness indicators to deactivate the limiters around smooth extrema. Therefore, small violations of the maximum principle could occur.
High order accuracy can be achieved by instead imposing the global bounds
Note that (2) is essentially a one-step discrete form of the maximum principle, which states that a solution must remain bounded by its maximum (and/or minimum) value at the beginning of the simulation. We will refer to schemes that satisfy (2) as maximum principle preserving (MPP). Higher than second-order full discretizations that strictly satisfy (2) are a more recent development and include zhang2011maximum ; xiong2015high ; chen2016third .
In the present work, our starting point is a spatial discretization, like WENO or the method presented in lohmann2017flux for finite elements, that uses limiters to achieve a considerable reduction of oscillations without degrading the formal order of accuracy of the solution (but may violate (2)). We modify this discretization by performing an extra flux limiting step that enforces (2) strictly. Thus the resulting method employs a two-tier limiting strategy.
Bound-preserving time discretization
Most of the works cited above are based on method-of-lines finite volume or (dis)continuous Galerkin finite element discretizations. A key difficulty in this area is to find a time integration scheme that preserves the boundedness properties of the semi-discrete scheme. This difficulty is commonly solved by applying a strong stability preserving (SSP) Runge-Kutta time discretization. This means that the schemes are limited to 4th- or 6th-order accuracy, depending on whether an explicit or implicit method is used SSPbook . Furthermore, existing high-order implicit SSP methods are not A-stable, so that such schemes will (when applied to (1)) be subject to severe time step restrictions even if an implicit time discretization is used.
In arbogast2020third , the authors take a different approach by combining the backward Euler method with a third-order fully-implicit Runge-Kutta method in order to have L-stability. However, the spatial discretization is based on WENO reconstruction, which leads to a scheme that does not strictly satisfy the maximum principle. Herein, we provide a general technique that allows the use of any high-order Runge-Kutta method with a spatial discretization based on WENO reconstruction. The time discretization need not be SSP, can be of arbitrarily high-order, and can be chosen to be diagonally implicit. To obtain a full discretization that satisfies the maximum principle, we combine the high-order method with a low-order MPP scheme based on backward Euler with local Lax-Friedrichs numerical fluxes. Therefore, we obtain anti-diffusive fluxes (or flux corrections) that contain corrections to the spatial and temporal components of the low-order scheme. The FCT method has been used before with fluxes that combine corrections in space and time; see for instance lee2010multistep ; lohmann2017flux ; xiong2015high ; yang2016high . However, those references are based on combining explicit schemes. As a result, their flux limited update is MPP only under a restricted time step. It is also worth mentioning feng2019time , wherein the authors use continuous Galerkin finite elements in space and discontinuous Galerkin finite elements in time to obtain a full discretization that is then modified via FCT to obtain an MPP method. The baseline discretization in feng2019time resembles an implicit scheme.
Schemes for problems with parabolic terms
Much work in this area has focused on the purely hyperbolic setting, since parabolic terms tend to have a smoothing effect and may make it less challenging to achieve discrete boundedness properties. Some recent works have focused on convection-diffusion applications where convection is dominant and it is still difficult to avoid oscillations or strictly satisfy the maximum principle zhang2012maximum ; chen2016third ; xiong2015high ; yang2016high . These methods are explicit, and will be subject to tight restrictions on the time step when the diffusion coefficient is not small. In addition, strict preservation of the maximum principle is fulfilled only if a time step restriction is satisfied. In this work, we consider implicit schemes. As a result, we obtain a high-order full discretization that is MPP with no time step restriction.
1.1 Our contribution
The techniques proposed here are closely related to those of exGMCL , which provided fully discrete explicit schemes for hyperbolic problems. Here we extend the approach to problems that include diffusion and to implicit time integration. We make use of two ideas based on general techniques that have long been used in this area. The first is that of combining a low-order method that satisfies the desired property with a high-order method that may not, in such a way that the high-order "correction" is guaranteed not to break the property. The second is reminiscent of Harten’s Theorem harten1974method ; harten1997high , and consists of writing a scheme as a sum of updates, each of which is proportional to the difference between the current state and a neighboring state. By bounding the neighboring states and the proportions, one obtains bounds on the updated solution.
The main contributions of our new schemes are: i) strict enforcement of the maximum principle (2); ii) arbitrarily high-order accuracy in time and space; iii) application not only to hyperbolic problems but also to convection-diffusion equations; and iv) linear stability and MPP under arbitrarily large time step sizes. The result is a framework to obtain an arbitrarily high-order full discretization for the convection-diffusion equation that is strictly MPP for time steps of any size.
The rest of this manuscript is organized as follows. In Section 2, we review a general and well-known framework for discretizations of convection-diffusion problems based on limiting flux corrections. This methodology relies on low- and high-order schemes, which we present in Sections 3 and 4, respectively. Afterwards, in Sections 5 and 6, we use two flux limiting techniques that guarantee the scheme is MPP. In Section 7, we discuss how to impose the MPP on the intermediate solutions within the stages of the Runge-Kutta method. In Section 8, we present four one-dimensional tests and four two-dimensional tests. In our numerical examples, we use a 5th-order WENO discretization and a 5th-order singly diagonally implicit Runge-Kutta (SDIRK) method. Concluding remarks are given in Section 9. For completeness, in A, we discuss how we solve the algebraic equations associated with the implicit discretizations.
2 Flux correction
We are interested in a finite volume spatial discretization of the convection-diffusion equation (1). For simplicity, we assume the domain is a hyperrectangle and prescribe periodic boundary conditions on . The initial condition is given by
We partition into cells , where . Let denote the boundary of , denote the face shared between cells and , and denote the volume and the area of and , respectively, and denote the set of indices of neighbors of that share a face with it. To obtain a finite volume semi-discretization, we integrate (1) over each cell, apply the divergence theorem to the convective and diffusive terms, and approximate the fluxes via
on each face . Here is the unit outward normal on face . Doing so, we get the spatial semi-discretization
where is the average of the solution over cell . We consider full discretizations of the form
is a time-averaged approximation of the combined flux . The specific form of depends on the time integration scheme. In Sections 3 and 4, we consider backward Euler and diagonally-implicit Runge-Kutta (DIRK) methods, respectively.
The basic idea used in this work was proposed almost 40 years ago in the hybrid schemes of Harten & Zwas harten1972self and in the flux-corrected transport (FCT) algorithm of Boris & Book boris1973flux . It has been employed in countless other methods proposed since then, and is explained neatly for instance in (leveque1992numerical, , Section 16.2) and kuzmin2012flux . The idea is to define two different numerical fluxes and , leading to two schemes, each of the form (5). The low-order flux is inaccurate but yields an update (5) that satisfies a desired bound. The high-order flux is more accurate but does not generally satisfy the bound. We then apply the scheme
where the limiters depend on and are chosen to maximize the accuracy while still enforcing the bound. We refer to the difference as the flux correction, since it does not approximate the flux itself but rather provides a high-order correction to it. Herein, we approximate the fluxes using the point value at the center of each face.
As described in Section 3, using backward Euler and the Lax-Friedrichs numerical flux yields low-order fluxes that are MPP with no restriction on the time step size. In Section 4, we use arbitrarily high-order methods based on WENO reconstruction and DIRK time integration to define . We present two approaches to choosing the values of the limiters; the first is based on Zalesak’s FCT limiters zalesak1979fully (but applied in time as well as space) while the second follows the recently-proposed monolithic convex limiting approach exGMCL .
3 Low-order scheme
In this section, we define low-order convective and diffusive fluxes, which we denote by and , respectively.
For the convective fluxes we use the local Lax-Friedrichs flux, also known as Rusanov’s flux, given by rusanov1961calculation
where denotes the solution average over cell and is an upper bound for the wave speed of the Riemann problem associated with face . In general, is a function of . For simplicity, we omit the dependence herein. For a uniform and structured grid, the diffusive flux can be taken as
Note that the scheme is mass-conservative since the right hand side above is antisymmetric with respect to the exchange of and . For the purely hyperbolic case () of (1), Guermond and Popov guermond2016invariant considered an equivalent representation of (9) in terms of upwinded averages
We remark that these states appear naturally in certain approximate Riemann solvers, where one assumes that the solution of the Riemann problem with data and consists of two traveling discontinuities, as shown in Figure 1. Then, is the intermediate state in the Riemann solution; see for instance leveque2002finite . These states satisfy guermond2016invariant ; kuzmin2020monolithic
We also define the arithmetic average states:
We now introduce the following quantities:
Importantly, the states are also bound preserving:
Finally, we can write the low-order spatial semi-discretization (9) as follows:
Since , and , the low-order spatial semi-discretization (13) is local extremum diminishing (LED) jameson1993computational . Namely, if is a local maximum, so can’t increase. Similarly, if is a local minimum, so can’t decrease. This leads to a semi-discretization that is MPP. By using the implicit Euler method, we achieve the MPP property in time also, with no restriction on the step size; see for example horvath1998positivity ; magiera2020constraint . The full low-order scheme is thus
where we use the superscript to refer to the low-order solution based on the backward Euler scheme with Lax-Friedrichs numerical fluxes. We can write (14) in terms of the low-order fluxes
[Time step restriction with Forward Euler] If we discretize (13) in time using the forward Euler method, we obtain
The solution is MPP provided
where is the mesh size.
4 High-order scheme
In scheme (6), are fluxes that improve the accuracy in space and time of the low-order fluxes . In this section, we define . Let us first define high-order convective and diffusive fluxes, which we denote by and , respectively.
For the high-order convective fluxes, we apply (7), after replacing the cell averages by high-order pointwise reconstructed values. Let denote a high-order approximation of in cell , based on weighted essentially non-oscillatory (WENO) reconstruction liu1994weighted ; jiang1996efficient . Then we set
For a uniform and structured mesh, the high-order diffusive fluxes can be given by
In principle, these fluxes should be integrated over each face , but the reconstruction required for this quadrature is very expensive. An economical alternative, which we use in Section 8, is to approximate the spatial integrand by the value at the midpoint of the face; this approach often reaps most of the benefits of the high-order WENO reconstruction at a reduced cost qiu2002construction . We correspondingly replace by in (17) and (18), where denotes the midpoint of face .
To integrate in time, we use high-order -stage diagonally implicit Runge-Kutta (DIRK) methods. Let , , and (with ) denote the Butcher coefficients of the DIRK method. The intermediate RK approximations to the cell averages are given by
The RK update is given by
are the high-order fluxes. Here we use the superscript to refer to the high-order solution based on DIRK schemes with WENO reconstruction. We use the high-order flux in scheme (6).
The high-order solution is not MPP due to violations introduced by the discretizations in space and time. Using scheme (6), with the flux limiters that we introduce in the next two sections, guarantees the RK update is MPP. However, the intermediate solutions might still violate the maximum principle. For some applications, preservation of the maximum principle is also needed for the intermediate solutions; we discuss ways to impose this in Section 7.
5 Flux corrected transport (FCT) limiting
In this section, we use the FCT method of boris1973flux ; zalesak1979fully to determine the flux limiters . Although both the low- and high-order methods are implicit, and hence require solving algebraic systems, the limiters are computed explicitly, as described below. The flux-limited update inherits the MPP properties of the low-order solution; that is, the solution is MPP with no time step restriction.
In the rest of this section, we follow kuzmin2012flux . We will determine flux limiters that guarantee
Using condition (23), this guarantees . The limiters are computed as follows:
Calculate the sum of positive and negative flux corrections:
Use the sums and the bounds , given by (23), to compute
Define the limiters by
Clearly, . The satisfaction of (23) is proven as follows:
and similarly for the lower bound . Since and , we have
which, by conservation of , implies the scheme (22) is mass conservative.
[Iterative FCT] In some of the numerical experiments from Section 8, we use the iterative FCT method to recover the high-order accuracy from the baseline scheme. The basic idea behind iterative FCT is to consider the quantity , which is the flux excluded by the limiters, and perform an extra limiting step given by
6 Global monolithic convex (GMC) limiting
In this section we use a different technique to determine limiters in (22) that will guarantee the MPP property (2). Namely, we follow the global monolithic convex (GMC) limiting approach from exGMCL . As in the previous section, and are the low- and high-order fluxes given by (15) and (21), respectively. Before defining the limiters , we need to rewrite scheme (22) in a form like that given in (exGMCL, , Section 3.2). Let us define the following quantities:
In terms of these, the low-order scheme (14) can be written as
Let us consider (6) where the high-order flux is computed via (21) at the beginning of the time step. The low-order flux is given by (15) and is treated implicitly; i.e., . Using (25), we rewrite scheme (6) as follows:
The MPP properties of (28) are guaranteed by the following theorem.
(Maximum principle) Let
Assume and that ’s are chosen to satisfy
Then given by (28) satisfies with no time step restriction.
7 Maximum principle preservation for intermediate stages
The procedures outlined in Sections 5 and 6 guarantee preservation of the maximum principle for the new solution , but not necessarily for the intermediate stages . For some applications, it may be important to guarantee the maximum principle for the intermediate stages, particularly if the system is not defined for values outside certain bounds. We consider two possible approaches:
Diagonally implicit Runge-Kutta methods
Apply the limiting procedure to each stage value of a given DIRK scheme. This can (at least formally) reduce the order of accuracy of the time integration scheme. In exGMCL , this approach was used with explicit Runge-Kutta methods. The authors did not observe loss of accuracy in their numerical experiments.
Methods with SSP stages
Start with a spatial semi-discretization that is MPP and let denote the time step under which it is MPP when discretized by the forward Euler method. Given a RK method, let denote the matrix of the Butcher coefficients , choose , and set
denote the column vector of lengthwith all entries equal to unity. Then if
7.1 Implicit Euler extrapolation methods
A particularly useful class of methods are those satisfying (32) for arbitrarily large values of . Such methods have stages that are unconditionally SSP (i.e., SSP under any step size) and can be constructed using extrapolation applied to the implicit Euler method (Hairer:ODEs2, , Section IV.9). These methods are nearly A-stable (specifically, they are -stable with close to 90 degrees) and can be constructed to have any order of accuracy. These methods are also highly parallelizable 2014_hork . The implicit Euler extrapolation algorithm for a method of order is given in Algorithm 1. For any fixed , this algorithm can be written as a Runge-Kutta method. As an example, we provide coefficients of the 4th-order method below. The coefficients are given in the standard Butcher form, although the implementation is done more efficiently using Algorithm 1.
For the second approach above, we need a spatial semi-discretization that is MPP, which can be obtained by applying the GMC limiters from Section 6 only to the spatial discretization. We refer to (exGMCL, , Section 2.2) for details. In Section 8.1.1, we test this methodology with a linear advection-diffusion problem in one-dimension. We recover the full accuracy of the underlying high-order scheme.
8 Numerical examples
In this section, we present a series of one- and two-dimensional numerical experiments to demonstrate the properties of the MPP algorithms we propose. For each problem, we consider the following four methods:
LLF-BE. Low-order (local Lax-Friedrichs) spatial discretization and backward Euler time integration; see Section 3 for details.
WENO-SDIRK. Fifth-order WENO spatial discretization and a fifth-order SDIRK time integration, whose Butcher tableau is given below; see Section 4 for details. This high-order scheme is used as the baseline high-order method for the following two MPP algorithms.
FCT-SDIRK. MPP algorithm presented in Section 5.
GMC-SDIRK. MPP algorithm presented in Section 6.
The fifth-order SDIRK method that we consider, which can be found in (kennedy2016diagonally, , Section 7.2.2) and references therein, has the following Butcher tableau:
In addition to the previous MPP algorithms, for the first one-dimensional test that we solve, we consider an algorithm that guarantees the intermediate solutions of the Runge-Kutta scheme are MPP. We do that only for one test to demonstrate that preserving the maximum principle for the intermediate stages does not destroy the accuracy properties of the underlying high-order scheme.
The spatial discretization is performed on uniform grids with elements. Let denote the -th element; then and for the one- and two-dimensional domains, respectively. The mesh spacing is denoted by and in the x- and y-direction, respectively. For the discretization in time, we use by default . To quantify the magnitude of the overshoots and undershoots, we report
Note that for any MPP solution. In practice, however, might be a small negative number on the order of machine precision, which indicates a small violation of the maximum principle. In practice it might be acceptable to clip these values, since the methods are conservative only up to machine precision. If the exact solution is available, we calculate and report the error
where is a fifth-order polynomial reconstruction of the numerical solution evaluated at . In addition, we report the corresponding experimental order of convergence (EOC).
8.1 Linear convection-diffusion
We start with the linear problem proposed in yang2016high . The problem is given by
with periodic boundary conditions. The coefficient controls the amount of dissipation and is the speed of advection. We take