Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting
Introduction
High-order numerical methods have been successfully applied to a wide range of applications [1], [2], [3]. These methods promise higher accuracy per degree of freedom when compared with traditional low-order methods, and have the potential to achieve high efficiency on modern computing architectures [4], [5], [6]. For example, in the field of computational fluid dynamics, such methods have seen particular success when applied to under-resolved turbulent flows, such as in the context of implicit large eddy simulation (ILES) [7], [8], [9]. However, a critical issue that must be addressed when applying high-order methods to convection-dominated problems is their robustness, especially in the context of nonlinear problems with discontinuous features such as shock waves [10], [11], [12].
In particular, discontinuous Galerkin (DG) methods have seen considerable success when applied to convection dominated problems [13]. These methods possess many desirable properties, such as arbitrary formal order of accuracy, and suitability for use with unstructured meshes. The robustness of DG methods is the subject of a large body of research [7], [14]. For scalar conservation laws and symmetric systems, the DG method satisfies a cell entropy inequality [15]. However, for general hyperbolic systems, such as the Euler equations, techniques such as flux differencing are required to ensure entropy stability [16], [17], [18], [19]. Furthermore, the use of high degree polynomials can introduce oscillations, and therefore limiters or artificial viscosity techniques are often used for bounds preservation, monotonicity, and shock capturing [20], [21], [22], [23].
An alternative to the above stabilization and limiting strategies is an approach developed by Guermond, Popov, and colleagues, based on invariant domain preserving (IDP) discretizations and convex limiting [24], [25], [26], [27]. A desirable property for numerical discretizations of hyperbolic conservation laws is invariant domain preservation: if the exact solution to the conservation law lies in a convex invariant set, then the numerical solution should as well [24]. This is a generalization of the concept of a discrete maximum principle, and will ensure that the discretization is bounds preserving, positivity preserving, and non-oscillatory. Suitable low-order invariant domain preserving (IDP) discretizations have been paired with high-order discretizations using convex limiting or algebraic flux correction strategies to obtain second-order accurate methods that preserve specific invariant domains [25], [26], [27].
In this work, we develop high-order discontinuous Galerkin methods that satisfy invariant domain preserving properties using a convex limiting strategy. The limiting strategy makes use of a novel sparse low-order IDP method whose stencil does not grow with the polynomial degree of the corresponding high-order method. Crucially, the accuracy of the low-order method does not degrade as the polynomial degree of the high-order method is increased, as is observed to occur with more naive graph viscosity approaches. Related strategies for sparsifying the convective operator for Bernstein basis finite element methods were previously developed by Kuzmin and colleagues [28]. The flux-corrected method is obtained by performing an efficient, dimension-by-dimension subcell flux correction procedure, blending the low-order IDP method with the high-order target DG method. The resulting method is conservative, and can satisfy any number of constraints on quasiconcave functionals specified by the user (cf. [25]). Since the accuracy of the low-order method does not decrease with the polynomial degree of the target method, we observe more accurate results using higher-order methods with a fixed number of degrees of freedom, even on problems with discontinuous solutions. This method can also be combined with a subcell resolution smoothness indicator to alleviate peak clipping effects near smooth extrema.
The structure of this paper is as follows. In Section 2, we formulate the high-order DG discretization, and state some key properties. In Section 3, we introduce a new low-order sparsified discretization that can be rendered invariant domain preserving using a graph viscosity approach. We develop a subcell flux correction strategy for blending the high-order (target) method and the low-order IDP method in Section 4. As specific examples, applications to the linear advection equation with variable velocity field and the Euler equations of gas dynamics are discussed. A number of numerical test cases demonstrating the effectiveness of the method on both scalar equations and systems of hyperbolic conservation laws are presented in Section 5. Finally, we end with some concluding remarks in Section 6.
Section snippets
Governing equations and discretization
Consider a system of hyperbolic conservation laws, with solution , . The flux function is given by . The spatial dimension is denoted , and the number of solution components is . The initial conditions are given by . Closely associated with the problem (1) is the following one dimensional Riemann problem, which will be important both for the definition of invariant domain preservation, and for the formulation of the discontinuous
Construction of the IDP low-order method
We now modify the discretization (13) with the goal of obtaining a method which is invariant domain preserving (IDP). In Guermond and Popov [24] this is done by adding a graph viscosity term that is based on the guaranteed maximum speed (GMS) of the hyperbolic system. One potential issue with this approach is that the amount of graph viscosity added to the discretization increases as the size of the discrete stencil increases. As a result, applying this technique to high-order methods with
Flux correction and limiting strategies
We now combine the low-order IDP discretization (33) with the high-order discretization (13), to obtain a bound-preserving scheme. To this end, we consider a forward Euler time discretization of both the low-order and high-order approximations (as mentioned above, the extension to high-order time integration is straightforward using SSP methods). The provisional high-order update for element is given by and the low-order IDP
Numerical examples
The method was implemented in the MFEM finite element framework [46], and is tested on a variety of benchmark problems, including scalar problems and hyperbolic systems, in 1D and 2D. For these test cases, unless stated otherwise, we integrate in time using the third-order SSP Runge–Kutta method [47]. The time step is chosen according to (38), with a CFL constant of , i.e. .
Conclusions
In this work, we have presented a discontinuous Galerkin spectral element method with convex limiting for hyperbolic conservation laws. This method preserves any specified set of invariant domain properties (e.g. local maximum principles, positivity of pressure and density, minimum principle for specific entropy, etc.). The method is based on an efficient dimension-by-dimension subcell blending of the target high-order (unlimited) DG-SEM method, and a low-order, invariant domain preserving
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
The author acknowledges H. Hajduk and D. Kuzmin for insightful conversations and comments on this work.
This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory, USA under Contract DE-AC52-07NA27344. LLNL-JRNL-808645. This document was prepared as an account of work sponsored by an agency of the United States government . Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees
References (58)
High order WENO and DG methods for time-dependent convection-dominated PDEs: a brief survey of several recent developments
J. Comput. Phys.
(2016)- et al.
On the eddy-resolving capability of high-order discontinuous Galerkin approaches to implicit LES/under-resolved DNS of Euler turbulence
J. Comput. Phys.
(2017) - et al.
An optimization-based approach for high-order accurate discretization of conservation laws with discontinuous solutions
J. Comput. Phys.
(2018) - et al.
High-order entropy stable finite difference schemes for nonlinear conservation laws: finite domains
J. Comput. Phys.
(2013) - et al.
Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws
J. Comput. Phys.
(2017) On discretely entropy conservative and entropy stable discontinuous Galerkin methods
J. Comput. Phys.
(2018)Limiters for high-order discontinuous Galerkin methods
J. Comput. Phys.
(2007)- et al.
Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws
Appl. Numer. Math.
(2004) - et al.
Invariant domain preserving discretization-independent schemes and convex limiting for hyperbolic systems
Comput. Methods Appl. Mech. Engrg.
(2019) Monolithic convex limiting for continuous finite element discretizations of hyperbolic conservation laws
Comput. Methods Appl. Mech. Engrg.
(2020)
Subcell flux limiting for high-order Bernstein finite element discretizations of scalar hyperbolic conservation laws
J. Comput. Phys.
Approximate tensor-product preconditioners for very high order discontinuous Galerkin methods
J. Comput. Phys.
Flux-corrected transport algorithms for continuous Galerkin methods based on high order Bernstein finite elements
J. Comput. Phys.
Finite element and finite difference methods for hyperbolic partial differential equations
Entropy-stable summation-by-parts discretization of the Euler equations on general curved elements
J. Comput. Phys.
Fully multidimensional flux-corrected transport algorithms for fluids
J. Comput. Phys.
Matrix-free subcell residual distribution for Bernstein finite elements: monolithic limiting
Comput. & Fluids
A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws
J. Comput. Phys.
Efficient implementation of essentially non-oscillatory shock-capturing schemes, II
J. Comput. Phys.
Entropy viscosity method for nonlinear conservation laws
J. Comput. Phys.
The numerical simulation of two-dimensional fluid flow with strong shocks
J. Comput. Phys.
CFD Vision 2030 Study: a Path to Revolutionary Computational AerosciencesTechnical Report NASA/CR-2014-218178
High-order CFD methods: current status and perspective
Internat. J. Numer. Methods Fluids
Efficiency of high order spectral element methods on petascale architectures
High Perform. Comput.
CEED ECP Milestone Report: Public Release of CEED 1.0Tech. Rep.
High-order matrix-free incompressible flow solvers with GPU acceleration and low-order refined preconditioners
High-order DNS and LES simulations using an implicit tensor-product discontinuous Galerkin method
Implicit LES of free and wall-bounded turbulent flows based on the discontinuous Galerkin/symmetric interior penalty method
Internat. J. Numer. Methods Fluids
An LES setting for DG-based implicit LES with insights on dissipation and robustness
Cited by (38)
High order entropy stable discontinuous Galerkin spectral element methods through subcell limiting
2024, Journal of Computational PhysicsFinite element-based invariant-domain preserving approximation of hyperbolic systems: Beyond second-order accuracy in space
2024, Computer Methods in Applied Mechanics and EngineeringOn the conservation property of positivity-preserving discontinuous Galerkin methods for stationary hyperbolic equations
2023, Journal of Computational PhysicsA positivity preserving strategy for entropy stable discontinuous Galerkin discretizations of the compressible Euler and Navier-Stokes equations
2023, Journal of Computational PhysicsCitation Excerpt :Secondly, we impose only one-sided positivity bounds for density and pressure. We do not perform any additional limiting, such as enforcing a local minimum entropy principle [22,34,35] or enforcing local bounds on the solution [36]. The motivation for this is that the minimum entropy principle is consistent with the compressible Euler equations but not the compressible Navier-Stokes equations.