GP-MOOD: a positivity-preserving high-order finite volume method for hyperbolic conservation laws

Abstract We present an a posteriori shock-capturing finite volume method algorithm called GP-MOOD. The method solves a compressible hyperbolic conservative system at high-order solution accuracy in multiple spatial dimensions. The core design principle in GP-MOOD is to combine two recent numerical methods, the polynomial-free spatial reconstruction methods of GP (Gaussian Process) and the a posteriori detection algorithms of MOOD (Multidimensional Optimal Order Detection). We focus on extending GP’s flexible variability of spatial accuracy to an a posteriori detection formalism based on the MOOD approach. The resulting GP-MOOD method is a positivity-preserving method that delivers its solutions at high-order accuracy, selectable among three accuracy choices, including third-order, fifth-order, and seventh-order.


Introduction
The modern trend in designing numerical methods for astrophysical flow simulations bears critical design principles driven by the vital need for high-performance computing (HPC). In today's HPC, the hardware progression of the memory capacity per compute core has become gradually saturated. This hardware trend pushes the HPC community to put vast efforts to find more efficient ways to best exercise computing resources in pursuing computer simulations. From the mathematical perspective, one desirable way is developing highly efficient numerical algorithms, which has become an important subject in computational fluid dynamics (CFD) research fields as part of utilizing HPC resources as efficiently as possible.
In this paper, we delve into advancing high-arithmetic-intensity numerical approximation by using the Gaussian Process (GP) modeling. This work has been reported in a series of papers, which have appeared in our past studies by Reyes et al. (2018aReyes et al. ( , 2019; Reeves et al. (2021); Bourgeois and Lee (2022). Our approach is based on the theory of GP regression that furnishes a polynomial-free method, extended to both finite difference and finite volume methods as a conservative high-order solver. In our former studies, we showed that GP can deliver high-order accuracy at (2R + 1)th order, where R is an integer value that represents the GP stencil radius in the unit of grid-scale, e.g., Δx. We designed GP as an alternative to the conventional high-order polynomial-based reconstruction algorithms, featuring attractive algorithmic benefits such as GP's selectable order of accuracy within a single algorithmic framework.

D. Lee & R. Bourgeois
In designing high-order algorithms for astrophysical simulations, it is necessary to implement stable shock-capturing mechanisms to advance discrete solutions stably in the vicinity of sharp gradients. A majority of well-known shock-capturing methods rely on the so-called a priori shock-detection paradigm. This is the classical approach in most of the widely-used polynomial-based reconstruction schemes, in which each method detects the magnitude of local flow gradients in an a priori fashion before updating the solution, to evolve it stably while meeting underlying physical principles (e.g., positivity, conservation). The mathematical treatments relevant to a priori shock detection mechanisms are truly nonlinear, demanding to calculate cell-by-cell local switches, such as slope limiters, to prevent the evolution of unphysical oscillations near strong gradients. The use of nonlinear limiters accounts for an important portion of the computational cost of the simulation at the price of numerical stability, which is essential for non-smooth flow simulations. Putting this into a different perspective, one gains simulation speed-up by deactivating the nonlinear limiters of a shock-capturing algorithm in smooth flow simulations. For instance, our 2D smooth flow experiments show that one can obtain a factor of two or three gain, revealing the added expense of executing nonlinear switches in modern shock-capturing schemes. Another undesirable consequence of nonlinear limiters is the inherent numerical dissipation. It is the main operational mechanism of nonlinear limiters to gain stability at the cost of adding numerical dissipation by switching to a more stable low-order approximation at shocks and discontinuity. The study in Kent et al. (2014) shows that a large increase in numerical diffusion and dispersion errors are observed at the first time step in limited schemes, significantly reducing the effective grid resolution compared to the corresponding unlimited schemes.
The above discussion leads us to consider a completely different shock-capturing paradigm, called an a posteriori scheme. The first a posteriori shock-capturing method was proposed in Clain et al. (2011) and further investigated in Diot et al. ( , 2013; Diot (2012). Referred to as the MOOD (Multidimensional Optimal Order Detection) method, the new paradigm advances shock-dominant discontinuous solutions via the a posteriori fashion of "repeat-until-valid." The main focus of this paper is to combine GP's high-order linear spatial reconstruction methods for enhanced solution accuracy and the MOOD approach for an a posteriori shock-capturing strategy within GP. In what follows, we provide a brief description on how these two methods can be integrated to a new a posteriori GP-MOOD algorithm.

GP linear spatial reconstruction
By being a polynomial-free, radial kernel-based algorithm, GP is well-suited to reconstruct a pointwise fluid value at any arbitrary location using volume-averaged fluid quantities in a local stencil, called a GP stencil of radius R. A GP stencil is genuinely multi-dimensional. Its shape resembles a blocky-diamond centered at the (i, j) cell in the case of 2D, and a blocky-octahedron centered at the (i, j, k) cell in 3D. See Fig. 1 for two simple examples in 2D. The solution accuracy of the GP reconstruction on these GP stencils in Fig. 1 follows linearly as (2R + 1)th accurate, which is achieved easily by varying the size of the local GP stencil radius, R = 1, 2, . . . . The (2R + 1)th order GP reconstructor is primarily derived by formulating a GP kernel. In our former studies, we have used one of the popular GP kernels called the square exponential kernel (or SE), denoted by K and defined by where the length hyperparameter, , controls the characteristic length scale of the functions in the GP function space distributed with a prior mean function m(x) and a prior covariance function K(x, y). With this SE kernel, one applies the conditioning property of Bayes' theorem to the joint Gaussian distribution on the given grid cell data q ij to make GP's inference on f (x * ) given q ij . Denoted as x * = x i is a new point at which GP makes a probabilistic prediction of a function evaluation f (x * ). For example, x * = x i±1/2,j if one solves Riemann problems at cell face-centers. Here, the choice of functions f is agnostic and is only available probabilistically in GP in the sense of GP regression. Assuming a zero mean, the conditioning property furnishes a new pointwise posterior mean functioñ m(x * ) given bym  , x m ). The data-independent vector z T * = k T * K −1 is called the prediction vector, following the same convention in Reyes et al. (2018bReyes et al. ( , 2019; Reeves et al. (2021).
To use GP for finite volume methods where the input data q ij is given as volumeaveraged quantities (rather than pointwise quantities q ij in finite difference methods), we further integrate the kernels k * and K in Eq. (2.2) (see Reyes et al. (2018aReyes et al. ( , 2019; Reeves et al. (2021); Bourgeois and Lee (2022) for details). The output is a new GP finite volume reconstructor,m where z * and C are the integral versions of k * and K, respectively. We note here that Eq.
(2.3) is a linear reconstruction without any limited switches, unlike the conventional nonlinear limited polynomial-based reconstructions. The resulting GP reconstruction varies its order of accuracy depending on the size of the GP stencil (see Fig. 1) where the local data q ij is sampled from. The sizes of the GP kernels, k * and K, change accordingly to the size of the GP stencil, which are N × 1 and N × N , respectively. In the case of a uniform grid calculation, they need to be computed, transposed, inverted, and saved at the initial step only once when a computational grid is configured. They are then reused during simulations. For adaptively varying grid configurations such as adaptive mesh refinements (AMR), one can pre-compute them on each expected AMR grid resolution, save them once-and-for-all and reuse them during simulations. In both, it is shown that GP provides computational efficiency in computing (2R + 1)th accurate reconstructed pointwise Riemann states at any arbitrary locations, x * . This paper adopts multiple Gaussian quadrature points along each cell face.

The a posteriori MOOD shock-capturing method
The GP-MOOD method builds upon the existing MOOD methods by Clain et al. (2011);Diot et al. ( , 2013; Diot (2012). The primary difference between our GP-MOOD and the conventional MOOD approaches are two-fold, including (i) how the multidimensional spatial reconstruction is calculated (i.e., high-order GP vs. high-order polynomials) and (ii) the new addition of "Compressibility-Shock Detection" or CSD criterion to the MOOD iterative loop. See Fig. 2.
As displayed in Fig. 2, GP-MOOD bears the baseline "order-decrement" architecture of the original MOOD loop, which selectively re-calculate pointwise Riemann states U n ij,gm on "failed cells" according to the following three failed-cell detection criteria: (i) PAD (Physical Admissible Detection): positivity preservation on density and pressure variables, (ii) NAD (Numerical Admissibility Detection): numerical validity that monitors CAD (Computer science Admissibility Detection, i.e., NAN & Inf) and DMP (Discrete Maximum Principle) that monitors any excessive numerical oscillations, (iii) CSD (Compressibility-Shock Detection): compressibility and shock strengths, i.e., ∇p/p > σ p and ∇ · V < −σ v , where σ p > 0 and σ v > 0 are (heuristically) tunable threshold parameters. The GP-MOOD method uses these detection checks and drops GP's reconstruction order from the highest 7th-order GP-R3, to 5th-order GP-R2 and 3rdorder GP-R1, and down to the safest (also most diffusive and positivity-preserving) first-order. In general, the highest order solution in GP-MOOD can be any arbitrary high-order (2R + 1)th accurate GP scheme beyond the current choice of 7th-order, which needs to extend the GP stencil size further, e.g., the 41-point blocky-diamond stencil for the 9th-order GP-R4. Such an increase in the GP stencil size will result in higher arithmetic intensity, thereby slowing down the overall time-to-solution. Interested readers are encouraged to follow our recent study (Bourgeois and Lee 2022) for further details on GP-MOOD.

Results and Conclusion
We display two numerical results in Fig. 3. Shown on the left panel is a grid convergence result, tested on the 2D nonlinear isentropic vortex advection problem Shu (1998). The results of L 1 errors are reported on four different grid resolutions, N x = N y = 50, 100, 200, and 400. As shown, the convergence rates of the three GP-MOOD methods on the corresponding diamond GP stencil follow the analytical convergence rates (dotted lines) of (2R + 1), showing the expected 3rd-, 5th-, and 7th-order rates.
On the right, we present a 1D shock-tube test called the Shu-Osher problem (Shu and Osher 1989). Simulated results are displayed in a zoomed-in view of the entire profile to focus on the solution comparison over the high-frequency region. Overall, all results produce acceptable density profiles capturing the assumed high-frequency amplitudes reasonably well, conforming with the reference solution. The amplitude closest to the reference profile is achieved by GP-MOOD7, followed by GP-MOOD5, GP-MOOD3, and POL-MOOD3. It is pretty impressive to see how closely the solutions produced by GP-MOOD5 and GP-MOOD7 follow the reference profile, with only on the grid resolution 16 times lower than the reference solution.
Next, we test GP-MOOD on a highly compressible astrophysical problem involving strong shocks and discontinuities to demonstrate our method's robustness and shockcapturing capability. The test problem initializes two Mach 800 jets, one at the top of   the square domain [0, 1.5] × [0, 1.5], and the other at the bottom. There are two narrow slits, 0.7 ≤ x ≤ 0.8, at y = 0 and 1.5, through which the Mach 800 jets are injected into the domain via the inflow boundary condition fixed by the jet condition. See Bourgeois and Lee (2022) for more details. The two jets undergo a head-on collision, producing highly turbulent fluid motions that are progressively amplified as the two jets continue to make their ways in the opposite directions for t ≥ 0.003. The results in Fig. 4 illustrate that GP-MOOD can produce highly accurate turbulent flow dynamics of the two highly compressible jet evolution dynamics closely relevant to astrophysical applications. We remark that the algorithm's accuracy, stability, and positive-preserving are extremely crucial to successfully run this test problem.
We conclude this paper with a summary. We have developed a new GP-MOOD algorithm for a high-order hyperbolic algorithm. GP-MOOD combines (i) the high-order solution property of the GP linear reconstruction schemes (thereby furnishing affordable computational approximation cheaper than the conventional nonlinear limited reconstruction methods) and (ii) a new improved MOOD a posteriori strategy (thereby reducing numerical dissipation in the MOOD detection criteria). Our GP-MOOD is a strong positivity-preserving method by design, and monitors shocks and discontinuities by detecting a sequence of conditions in the MOOD loop in the a posteriori fashion.