Comparison between a priori and a posteriori slope limiters for high-order finite volume schemes

There have been significant advancements in the development of numerical schemes for hyperbolic conservation laws over the past few decades. With increasing demands for precision and reliability in simulations across fields like fluid dynamics, astrophysics, and climate modeling, there is a growing trend toward conducting simulations at higher resolutions.

More recently, there has been a shift in focus towards high-order numerical schemes, typically referring to those of order three or higher. These schemes offer several advantages over their lower-order counterparts, including improved resolution of fine details, reduced numerical diffusion over extended time scales, and the capability to capture sharper gradients and discontinuities. Moreover, in problems with smooth solutions, higher-order schemes exhibit exponentially lower errors as the order increases, meaning they can achieve the same accuracy as a given low-order scheme with a significantly lower resolution.

However, they also come with a significant drawback when dealing with solutions that contain discontinuities: Using higher-degree interpolation polynomials can create unphysical oscillations and severely overshoot/undershoot the true solution bounds. In problems involving the conservation of a scalar $u(t,\mathbf{x})$, it may be essential for the numerical solution to satisfy the maximum principle, which defines a range of values represented by $M=\max_{\mathbf{x}}u_{0}(\mathbf{x})$ and $m=\min_{\mathbf{x}}u_{0}(\mathbf{x})$, where $u_{0}(\mathbf{x})=u(0,\mathbf{x})$ is the initial condition. These spurious oscillations violate the maximum principle condition $u(t,\mathbf{x})\in[m,M]$, making it challenging to enforce physical properties such as the positivity of density or pressure. In fact, Godunov’s theorem asserts that any linear scheme higher than first-order will not strictly preserve the maximum principle of a scalar conservation law [4]. Maximum principle violations are catastrophic for solvers of many nonlinear conservation laws, rendering linear high-order schemes useless for such problems.

As far as linear schemes for conservation laws go, there are three popular spatial discretizations: finite difference (FD), finite volume (FV), and finite element (FE) methods. FD/FV methods can be made higher-order simply by increasing the size of the stencil used to interpolate spatial derivatives from nodes/cells and their neighbors. FV methods are particularly attractive for conservation laws because they can easily be defined to conserve solution quantities to numerical precision. Additionally, as demonstrated in the text, they remain stable with explicit Runge-Kutta methods even for large CFl factors.

In the realm of scientific computing, finite element (FE) methods, especially the discontinuous Galerkin (DG) method, have gained increasing popularity. In DG, the numerical solution is represented using a polynomial basis within a domain region known as an element. By using higher-degree polynomials across fewer elements, the number of discontinuities between solutions is relatively low compared to an equivalent FV scheme. This aspect is crucial, particularly when addressing these discontinuities involves costly Riemann solvers.

Many techniques have been developed to help numerical solutions preserve local monotonicity and avoid these violations. For instance, the artificial viscosity method adds a physical diffusion term to the underlying equation in order to smooth numerical solutions near discontinuities or elsewhere [5, 6, 7, 8, 9].

Another approach is to utilize a slope limiter, an operator designed to reduce the slopes of the piecewise reconstruction of a solution, thereby enforcing monotonicity. The key point is that this operator is nonlinear, allowing schemes to break free from the constraints of Godunov’s theorem. There are many examples of limited schemes in the literature. The MUSCL scheme introduced by Van Leer [10] is a second-order scheme that strictly preserves a given maximum principle. Later came the popular TVD limiter [11] and the PPM scheme [11]. The ENO [12] scheme and WENO scheme [13] were introduced as a means to preserve the high-order accuracy of FV schemes while reducing oscillations near discontinuities.

A strictly maximum principle-preserving (MPP) scheme working at arbitrarily high order was introduced by Zhang & Shu [1, 2]. Their method calculates the numerical solution as a convex combination of high-order and first-order updates. It is guaranteed to satisfy the maximum principle if a sufficiently small CFL factor is used along with first-order forward Euler integration or another Runge-Kutta method equivalent to a convex combination of Euler steps. These methods are known as Strong-Stability-Preserving (SSP) integrators. While the popular fourth-order explicit Runge-Kutta method is not classified as Strong Stability Preserving (SSP), it exhibits quasi-SSP behavior for certain problems [14]. Therefore, we will evaluate its compatibility with MPP schemes.

it behaves as a quasi-SSP integrator, enabling schemes to maintain a strict maximum principle.

While slope limiters are popular tools when solving discontinuous problems with FV schemes, and can even guarantee the avoidance of maximum principe violations, they have drawbacks for FE. These popular limiters often necessitate reducing the degree of interpolating polynomials to zero to prevent overshoots and undershoots. However, in FE methods where polynomial elements cover large spatial regions, reducing an element to first-order significantly compromises the accuracy of the overall solution.

An important advancement in addressing this issue of slope limiting for FE was made by considering a posteriori limiters. Unlike traditional a priori limiters, which use data at $t^{n}$ to compute a limited update at $t^{n+1}$ in explicit schemes, a posteriori limiters compute a candidate solution at $t^{n+1}$ and revise it until it satisfies the problem requirements. This concept, introduced by Krivodonova [15] and Clain et al. [3] as MOOD, involves reducing an element solution by one degree at a time, limiting the solution slope while only reverting to first-order when necessary.

Vilar & Abgrall [16] introduced a novel concept to further localize changes made to the high-order solution by the limiter. In their approach, a subgrid of FV cell averages is interpolated from any FE of the candidate solution that violate positivity or other prescribed condition. In FV, a posteriori limiting is executed on a cell-by-cell basis, identifying to called troubled cells whose fluxes need to be revised [3, 17, 18]. A robust fallback scheme computes the revision of this intermediate FV solution, after which a the high-degree polynomial of the element is reconstructed from these limited cell averages. This fallback scheme can be the FV MOOD scheme [3] or any of the other limited schemes discussed here.

The spectral difference (SD) method is another FE approach that can be shown to be equivalent to a quadrature free nodal DG scheme [19, 20]. Velasco et al. [18] took advantage of the fact that the SD method offers a natural interpretation as a nonuniform FV grid, avoiding the need to interpolate between two solution representations.

For FV schemes, the a priori and a posteriori limiting paradigms have fundamentally different implementations and trade-offs. Zhang & Shu’s limiter [1, 2] maintains a strict maximum principle for high-order schemes, but is computationally costly due to the reduced CFL condition and the many nodal reconstructions for the flux quadrature. The cost of the small CFL factor can be mostly mitigated by using an adaptive time-step size [21]. GPUs, renowned for their ability to execute numerous matrix multiplications as a single vectorized operation, might also potentially mitigate this issue. This is because each nodal reconstruction corresponds to a stencil that can be computed through matrix multiplication.

The cost of the nodal reconstructions is more severe in two dimensions or higher. While the two-dimensional finite volume (FV) flux integral can be reconstructed in various ways, Zhang & Shu’s limiter [1, 2] use a Gauss-Legendre quadrature. This quadrature method requires a node count per cell face that grows linearly with the polynomial degree. In contrast, transverse flux reconstruction [see e.g. 22, 23, 24] uses only one node per cell face, making it a significantly more economical option. We will investigate whether the Zhang & Shu [1, 2] a priori slope limiter can still be used in an MPP scheme when transverse flux reconstruction is exploited. Previous studies combining an a priori limiter with transverse flux reconstruction have resulted in significant violations of the maximum principle [23].

Another limitation of MPP, a priori limiters is that, at high-order and after long time-scales, their numerical solutions can be dominated by numerical artifacts and diffusion [25]. The high-order schemes in this case actually perform worse at long time-scales than, say, second-order MUSCL-Hancock. Perhaps the trade-off of a strictly MPP scheme is that is too stringent with the permitted slopes as a means to enforces monotonicty. It is unclear how the a posteriori limiting compares in this regard.

The a posteriori limiters presented here do not require the use of Gauss-Legendre quadrature in two dimensions, making them automatically more cost-effective. However, these limited schemes often violate the maximum principle of the problem at hand, as they lack a guarantee to the contrary. Vilar & Abgrall [16] discovered that extending the convex blending of revised and high-order fluxes to a region of neighbors around each troubled cell reduces the magnitude of these violations. Rueda-Ramírez et al [26] showed that a more sophisticated version of this blending can be formulated so as it guarantee the satisfaction of physical bounds on solution variables. Nevertheless, it remains unclear how well a high-order, MPP a posteriori limited scheme performs over long time scales and whether they suffer from the same dominance of numerical artifacts as the a priori limited schemes.

In this study, we develop a methodology to compare a priori and a posteriori limited explicit finite volume (FV) schemes at arbitrarily high orders. We chose FV as the base scheme over finite element (FE) methods, as modern FE schemes revert to FV when encountering violations, anyways. We solve the linear advection equation in one and two dimensions through a series of numerical experiments. The numerical solutions obtained from the two limiting paradigms are compared with respect to their maximum principle violations and their quality over long time-scales. Furthermore, we evaluate the computational costs of these methods on traditional CPU architecture as well as GPUs.

In Section 2, we give an overview the finite volume method and the various slope-limiting techniques compared in this work. Section 3 conveniently summarizes the a priori and a posteriori limited schemes studied in our numerical experiments. The results of those experiments are presented in Section 4, where we validate our high-order FV and MPP implementations and compare the high-order slope-limited schemes’ maximum principle violations, solution quality at long time-scales, and computational cost. Additional commentary regarding the experiments is provided in Section 5 and conclusions are drawn in Section 6.

We consider the scalar conservation law

for $\mathbf{f},\mathbf{x}\in\mathbb{R}^{d}$, where $u(\mathbf{x},t)$ is a conserved scalar, $\mathbf{f}$ is a flux function and $d$ is the number of dimensions. We begin with the one-dimensional case ($d=1$), and rewrite (1) into the initial value problem

In the finite volume formulation, the spatial domain is partitioned into finite intervals, or cells $I_{i}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$ and the cell volume average is defined at time $t^{n}$ as

where $h$ is the length of $I_{i}$, chosen to be uniform in this work. Combining (2) with (3) and applying the divergence theorem, the semi-discrete form is written

The nodal values $u_{i\pm\frac{1}{2}}=P_{i}(x_{i\pm\frac{1}{2}})$ are obtained by reconstructing $u(x,t^{n})$ as a polynomial $P_{i}(x)$ piecewise on each interval $I_{i}$ and evaluating it at the cell edges. $P_{i}(x)$ is a conservative interpolation polynomial such that $\int_{I_{i}}P_{i}(x)=\bar{u}_{i}^{n}$ and is found by taking the derivative of the Lagrange polynomial which interpolates the first moment of $u(x)$ in some kernel of cells containing $I_{i}$ [12, 27, 28]. The size of the kernel limits the accuracy the polynomial degree of $P_{i}(x)$ and, consequently, the order of accuracy of the finite volume scheme; $P_{i}(x)$ is of degree $p$ when it is reconstructed from a kernel of exactly $p+1$ neighboring cells, so a smooth $u(x)$ is approximated with $u(x)=P_{i}(x)+\mathcal{O}(h^{p+1})$ on $I_{i}$.

The evaluation of $P_{i}(x_{i\pm\frac{1}{2}})$ is summarized as stencil operations in Table 1 for various degrees $p$. We select only those stencils that are symmetric so as to simplify our implementation and minimize the number of cells that must be adjusted to enforce boundary conditions. For odd-degree polynomials (which require a kernel of even size), we take the average of the left- and right-biased stencils.

To resolve the Riemann problem presented by the nodal values $u_{i+\frac{1}{2}}^{-}$ and $u_{i+\frac{1}{2}}^{+}$ interpolated for cells $I_{i}$ and $I_{i+1}$, respectively, a numerical flux function is written as $\hat{f}_{i\pm\frac{1}{2}}=\mathcal{F}(u_{i+\frac{1}{2}}^{-},u_{i+\frac{1}{2}}^{+})$. $\mathcal{F}$ must be a Lipschitz continuous function of both arguments, physically consistent with $u$ ($F(u,u)=u$), and monotone, which is to say that it is non-decreasing in its first argument and non-increasing in its second argument [29]. We use here

which satisfies these requirements.

In the case of linear scalar advection with uniform advection velocity $a$, given by

the Riemann solver given by (5) is known as the upwind solution.

The semi-discrete scheme (4) is rewritten with numerical fluxes as

and the definition of the order $p+1$ finite volume dynamics, given by the high-order spatial discretization $\mathcal{L}_{p}$, is complete.