Monolithic Convex Limiting for Legendre-Gauss-Lobatto Discontinuous Galerkin Spectral Element Methods

We extend the monolithic convex limiting (MCL) methodology to nodal discontinuous Galerkin spectral element methods (DGSEM). The use of Legendre-Gauss-Lobatto (LGL) quadrature endows collocated DGSEM space discretizations of nonlinear hyperbolic problems with properties that greatly simplify the design of invariant domain preserving high-resolution schemes. Compared to many other continuous and discontinuous Galerkin method variants, a particular advantage of the LGL spectral operator is the availability of a natural decomposition into a compatible subcell flux discretization. Representing a high-order spatial semi-discretization in terms of intermediate states, we perform flux limiting in a manner that keeps these states and the results of Runge-Kutta stages in convex invariant domains. Additionally, local bounds may be imposed on scalar quantities of interest. In contrast to limiting approaches based on predictor-corrector algorithms, our MCL procedure for LGL-DGSEM yields nonlinear flux approximations that are independent of the time-step size and can be further modified to enforce entropy stability. To demonstrate the robustness of MCL/DGSEM schemes for the compressible Euler equations, we run simulations for challenging setups featuring strong shocks, steep density gradients and vortex dominated flows.

The algebraic flux correction (AFC) schemes that we review and modify in the present paper are based on the methodology that is currently known as convex limiting [6,8,22].The underlying design philosophy traces its origins to localized flux-corrected transport (FCT) algorithms for scalar conservation laws [15,27,28].The first extension to nonlinear hyperbolic systems was proposed by Guermond et al. [8].In contrast to Zalesak's multidimensional FCT limiter [29] and its edge-based generalizations to continuous finite element methods for the Euler equations [5,15,[30][31][32], convex limiting approaches enforce preservation of local and global bounds by constraining individual fluxes rather than sums of fluxes.In the explicit case, the local extremum diminishing (LED) and/or invariant domain preserving (IDP) properties of flux-limited approximations are shown using representations in terms of intermediate states that stay in convex admissible sets [8,13].
All of the aforementioned FCT algorithms belong to the family of AFC schemes in which the computation of a property-preserving low-order predictor is followed by an anti-diffusive correction stage.The monolithic convex limiting (MCL) methodology developed in [6] differs from such fractional-step approaches in that limited anti-diffusive fluxes are incorporated into the residual of the semi-discrete scheme.The resulting nonlinear system of ordinary differential equations has a well-defined steady state, and the use of implicit time integrators is an option.The IDP property of the explicit version can be shown following the analysis of the low-order (local Lax-Friedrichs) method in [33].Moreover, the validity of (semi-)discrete entropy inequalities can be enforced using limiter-based or dissipation-based fixes [34,35].
The first successful extensions of FCT and MCL to high-order finite elements [14,15,22,24,36] used Bernstein polynomials as local basis functions.In this context, a key to achieving optimal accuracy lies in the use of sparse discrete gradient/Laplacian operators and subcell flux limiting techniques.The discontinuous Galerkin spectral element methods (DGSEM) proposed by Pazner [37], Lin et al. [38], and Rueda-Ramírez et al. [39] extend subcell convex limiting of FCT type to Legendre-Gauss-Lobatto (LGL) bases.The underlying low-order method has the structure of the subcell finite volume scheme employed in [23].The high-order DGSEM discretization also admits a natural sparse representation in terms of subcell fluxes between neighbor nodes.Hence, there is no need for artificial flux reconstructions or decompositions.Moreover, the mass matrices of collocated LGL-DGSEM approximations are diagonal and the discrete gradient/divergence operators possess summation-by-parts (SBP) properties, which are needed to achieve entropy stability [40].
As an alternative to the LGL versions [37-39, 41, 42] of high-order FCT algorithms and sophisticated limiters for Bernstein finite elements [22,35,43], we introduce a tailor-made LGL-DGSEM counterpart of Hajduk's [22] subcell MCL scheme for conservation laws.In fact, the proposed methodology is also applicable to any other spatial semi-discretization that produces sparse discrete gradient operators with SBP properties, such as Gauss-DGSEM discretizations [44] or general SBP discretizations of nonconservative systems of balance laws [45,46].The flux constraints of the MCL procedure and steadystate solutions are independent of the time step.In the context of subcell flux limiting for the Euler equations of gas dynamics, the density, momentum, and total energy fluxes are limited sequentially to enforce local bounds for the density, individual velocity components, and specific total energy.If the pressure becomes negative, a simple scaling limiter is applied.No local bounds are imposed on the physical entropy because the limiter-based fixes proposed in [34,35,43] guarantee entropy stability under less restrictive constraints.The above limiting strategy enables us to achieve high resolution without sacrificing any important properties or using impractically small time steps.
The remainder of this paper is organized as follows.In Section 2, we briefly present the LGL-DGSEM, derive the LGL-DGSEM subcell MCL method, and discuss some of its properties.In Section 3, we use the LGL-DGSEM/MCL method to perform challenging simulations of the compressible Euler equations, and present some comparisons with FCT/IDP strategies.Finally, we draw our conclusions in Section 4.

Numerical Methods
In this work, we deal with hyperbolic systems of conservation laws of the form where Ω ⊆ R D is a computational domain.The number of space dimensions is D ∈ {1, 2, 3}.The vector u( x, t) ∈ R neq of conserved quantities depends on the space location x and the time instant t.The flux function 1) is equipped with an initial condition, u(•, 0) = u 0 , and suitable boundary conditions on ∂Ω.
For brevity and better readability, we introduce the methods under investigation in the simple context of a one-dimensional (D = 1) conservation law or system.All algorithms to be discussed admit straightforward tensor-product extensions to two and three space dimensions, and to curvilinear grids.

The Discontinuous Galerkin Spectral Element Method
Let T = {Ω 1 , . . ., Ω K } be a tessellation of the domain Ω into K nonoverlapping elements.Within each element, we approximate the solution u by a polynomial of degree N .A piecewise-polynomial DG approximation u DG ≈ u may be discontinuous at the element interfaces.We seek u DG in the space φ| Ω e ∈ P N (Ω e ) ∀Ω e ∈ T }.
Restricting our attention to a single element Ω e , we multiply (1) by an arbitrary polynomial test function φ ∈ (P N (Ω e )) neq , integrate the weighted residual over Ω e , and perform integration by parts to obtain the weak form of the local conservation law.Since u DG is generally not uniquely defined at the element interfaces, we calculate f ≈ ↔ f • n| ∂Ω e using an approximate Riemann solver that receives two one-sided limits and returns a numerical flux.
The Legendre-Gauss-Lobatto (LGL) discontinuous Galerkin spectral element method (DGSEM) is a so-called nodal collocation variant of the DG method.It produces discrete gradient/divergence operators that possess summation-by-parts (SBP) properties [47].The restriction of u DG to Ω e is represented using Lagrange basis functions that are associated with (N +1) D LGL interpolation points.The quadrature rule for numerical integration on Ω e uses the LGL collocation nodes on the reference element Ω = [−1, 1] D .A mapping F e : Ω → Ω e is used for transformations from the reference space to the physical space ( ξ → x for ξ ∈ Ω and x = F e (ξ) ∈ Ω e ).After some manipulations, the evolution equation for the ith local degree of freedom of a one-dimensional LGL-DGSEM discretization of (1) on Ω e can be written as [39,48] where J denotes the constant determinant of the Jacobian of the mapping from the reference element, ω i denotes the reference-space quadrature weight, and δ ij is the Kronecker delta of the node indices i and j.The numerical fluxes f(0,L) and f(N,R) are calculated using the inner and outer limits of u DG on the boundaries of the element Ω e = (x e 0 , x e N ) containing the LGL nodal point x e i .The strong form derivative matrix S = ( Sik ) N i,k=0 admits the representation where B := diag(−1, 0, . . ., 0, 1) is the so-called boundary evaluation matrix.
The entries are defined using the derivatives of the Lagrange basis polynomials whose entries we denote by S ik , the discretized volume integral can be expressed in terms of two-point numerical fluxes f * (i,k) [49].The semi-discrete scheme is equivalent to (3) if the standard average f * (i,k) = (f i + f k )/2 is used.However, additional robustness can be achieved with other choices of the volumetric numerical flux f * (i,k) .For instance, some two-point approximations to fluxes of the Euler equations guarantee kinetic energy preservation [50], entropy conservation/dissipation [51,52], pressure equilibrium preservation [53], or all of these properties together [54,55].
All diagonal-norm SBP discretizations of conservation laws (and hence also the LGL-DGSEM considered here) can be written in the so-called flux-differencing form [49] where the indices −1 and N + 1 refer to the outer states.The symmetric and consistent fluxes f DG (i,j) = f DG (j,i) are defined by [39,49] f Note that the flux f DG (i,j) is multiplied by the one-dimensional unit normal n (i,j) ∈ {−1, 1} in (5).The normal fluxes n (i,j) f DG (i,j) are anti-symmetric, that is, n (i,j) f DG (i,j) = −n (j,i) f DG (j,i) .Hence, (4) has local (subcell-level) conservation properties, as required by the Lax-Wendroff theorem [56].
Remark 1 Let ∆x i = Jω i .Then (5) corresponds to the subcell finite volume scheme We adopt the two-subscript notation because it is better suited for flux-based finite element discretizations.
Remark 2 Since the LGL-DGSEM is a diagonal-norm SBP operator, its representation in the flux-differencing form ( 5)-( 8) is readily available.Other DG approximations with dense mass matrices need the application of a sparsification operator to recover the flux-differencing form.Examples of decompositions into subcell fluxes can be found, e.g., in [22,26,35].For a DG method using Bernstein polynomials of degree N > 1 as local basis functions, neq sparse linear systems of size (N + 1) × (N + 1) need to be solved for each element in each Runge-Kutta stage [22,35].Vilar [25] showed that it is possible to obtain a flux-differencing formula for any (modal or nodal) representation of the DG solution if one expresses the test function as a combination of so-called subresolution basis functions and exploits existing relationships to the histopolation theory.An adaptation to unstructured triangular grids was proposed by Vilar and Abgrall [26], who parametrized u DG in terms of subcell averages that satisfy a two-dimensional version of (5).The calculation of subcell fluxes involves solving small linear systems with sparse graph Laplacians again.
Remark 3 Mateo-Gabín et al. [44] showed that (both standard and split-form versions of) the Legendre-Gauss DGSEM scheme can also be written in the flux-differencing form with explicit staggered fluxes and a diagonal mass matrix.As a result, most of the algorithms to be presented in this paper are applicable to the Legendre-Gauss DGSEM.However, the treatment of inter-element fluxes and projection operators requires additional analysis and, possibly, appropriate modifications.

Monolithic Convex Limiting
The monolithic convex limiting (MCL) methodology [6,22,35] is a subcell flux correction procedure that combines a high-order baseline discretization with a compatible and invariant domain preserving low-order scheme.The validity of physical and numerical admissibility conditions is enforced using a representation in terms of intermediate states (similarly to the predictor-corrector approaches proposed in [8,13,[37][38][39]).Flux limiters for semi-discrete MCL schemes can be designed to enforce entropy stability conditions in addition to local and/or global maximum principles [43,57].To minimize the levels of loworder numerical dissipation, localized subcell limiting procedures are used for high-order finite elements [22,35,43].Moreover, sequential MCL algorithms for systems support the possibility of using individually chosen correction factors for different conserved or derived quantities [6,22].

Low-Order Invariant Domain Preserving Scheme
As explained in [23], one can obtain a low-order finite volume scheme that is compatible with the LGL-DGSEM discretization by interpreting the nodal values of the DGSEM scheme as mean values of the subcells.Let where f FV (i,j) is a low order numerical approximation to the flux between nodes i and j.Such a subcell FV scheme exhibits the same structure as (5).Therefore, the two schemes are compatible and can be hybridized.
It has been shown that ( 9) is invariant domain preserving (IDP) for the first-order Rusanov (also known as local Lax-Friedrichs, LLF) fluxes [37,39] f FV (i,j) = where λ max (i,j) = λ max (j,i) > 0 is an upper bound for the maximum wave speed of the Riemann problem with the initial states u i and u j .Estimation of this speed is addressed, e.g., in [58].In our notation, the presence of (i − j) in the dissipative part of (10) ensures that the flux f FV (i,j) = f FV (j,i) is symmetric.The multiplication by the unit normal n (i,j) = −n (j,i) makes it anti-symmetric.
Inserting the Rusanov fluxes (10) into (9) and using the forward Euler method for time integration yields a fully discrete version of the low-order scheme.Following Guermond and Popov [33], we write it in the form using the auxiliary bar states If the time-step size ∆t satisfies the CFL condition then the result u FV,n+1 i of the explicit update ( 11) is a property-preserving convex combination of the states u n i and u n (i,i±1) .We discuss the time-step restriction in Section 2.2.4.
The most important property of the so-defined LLF bar states is that they preserve all convex invariants of initial value problems for hyperbolic systems, as shown in [33] in the context of a continuous (multi-)linear finite element discretization.In fact, u (i,j) defined by ( 12) is the intermediate state of the HLL approximate Riemann solver [59].Positivity preservation and the validity of entropy conditions can be deduced from this interpretation.
Remark 4 The bar states (12) of the (semi-discrete or fully discrete) low-order LLF scheme are symmetric in the sense that u (i,j) = u (j,i) for any pair of adjacent nodes with local indices i ∈ {0, . . ., N } and j ∈ {i − 1, i + 1}.
Remark 5 To obtain compatible low-order IDP schemes for general high-order DG methods, it is necessary to first replace the discrete gradient/divergence operators with sparse approximations and then apply low-order dissipation (e.g., using a sparse graph Laplacian operator as in [14,22]).For our LGL-DGSEM scheme and, in general, for all diagonal norm SBP operators, the low-order IDP scheme ( 9) is readily available and compatible with the flux-differencing form (5) of (4).

Limiting Procedure
To enforce relevant inequality constraints, we replace ( 5) with (cf.[22,35]) In the simplest case, the hybrid subcell fluxes f(i,j) ∈ R neq are given by where • denotes the Hadamard (or component-wise) product.The scalarvalued components of α (i,j) ∈ R neq are weights that attain values between 0 and 1.The high-order DG method (5) and the low-order FV scheme (9) can be recovered using α (i,j) = 1 and α (i,j) = 0, respectively.As detailed in the next section, the computation of α (i,j) might not be numerically well posed.Therefore, it should be avoided in practical implementations if there is a direct way to calculate the fluxes f(i,j) .
, then the DG and FV schemes use the same two-point flux approximation on the boundaries of Ω e = (x e 0 , x e N ).In this case, formula (15) will produce f(i,j) = f DG (i,j) = f FV (i,j) for any choice of α (i,j) .This desirable property of boundary fluxes is specific to the LGL-DGSEM discretization because it includes the boundary nodes.It leads to a very local implementation of the limiting procedure.
Since the low-order component of f(i,j) is provably IDP, the purpose of subcell limiting is to constrain the anti-diffusive components of the fluxes f DG With this notation, the high-order update can be written as Using the representation of the flux difference f ) in terms of the bar states defined by (12), we find that We can now define the bar state of the high-order method as and cast (18) into the bar state form which has the same structure as (11).Following the derivation of MCL schemes for Lagrange and Bernstein finite elements [6,22,35], we replace the DG bar states defined by (19) with ) is a limited approximation to ∆ f(i,j) .In contrast to the low-order component u (i,j) = u (j,i) , the limited bar state (21) is generally not symmetric due to the skew-symmetry of ∆ f Lim (i,j) = −∆ f Lim (j,i) .The MCL bar states u Lim (i,j) should satisfy the same inequality constraints as u (i,j) and stay as close as possible to the high-order target u DG (i,j) .The forward Euler time discretization should be replaced with a high-order Runge-Kutta method.For SSP-RK schemes with forward Euler stages, the IDP property can be shown in the same way as for the low-order scheme [6,22].A general Runge-Kutta method may require flux limiting in time [60,61].
The convex limiting techniques employed in [8,13,[37][38][39] differ from MCL in that they split the computation of u DG,n+1 i into a low-order IDP update and an anti-diffusive correction stage.This predictor-corrector strategy is also used in older FCT-type algorithms for finite element discretizations of hyperbolic systems [5,15,[30][31][32].In contrast to MCL, the resulting schemes have no semidiscrete counterparts.Moreover, the bounds of the limiting constraints depend on the time step.Depending on the application, this peculiarity of FCT/IDP approaches may be an advantage or a disadvantage.
We will now describe the computation of the limited anti-diffusive fluxes ∆ f Lim (i,j) for the MCL version with generic bounds.Appropriate definitions of the bounds are discussed in Section 2.2.3.

Limiter for conservative quantities
The simplest limiting strategy for systems is to treat each equation as a scalar conservation law and to limit the anti-diffusive fluxes of each conserved variable individually.Let ρ be a scalar component of u.We denote by ρ (i,j) and ∆ f ρ,Lim (i,j) the corresponding components of u (i,j) and ∆ f Lim (i,j) , respectively.To keep ρ Lim (i,j) we impose the inequality constraints A positive/negative anti-diffusive flux ∆ f ρ,Lim may violate the upper/lower bound for ρ Lim (i,j) or the lower/upper bound for ρ Lim (j,i) .Introducing we define [6,22] ∆ f ρ,Lim It is easy to verify that conditions (22) are met for this choice of ∆ f ρ,Lim (i,j) .Moreover, there exists is a convex combination of the FV and DG fluxes.Kuzmin [6] noticed that the computation of α ρ (i,j) = ∆ f ρ,Lim (i,j) /∆ f ρ (i,j) is unnecessary and numerically ill posed in the case of a small nonvanishing denominator.The direct computation of the limited anti-diffusive flux ( 23) is therefore preferable in practice.

Sequential limiter for "primitive" quantities
In some situations, we are interested in imposing bounds on scalar quantities that are not included in the state vector u.If the quantity of interest represents the ratio of two conservative variables, we can use the sequential limiting approach proposed in [6,9].For instance, components of the velocity field, v = (ρ v)/ρ, and the total specific energy, E = (ρE)/ρ, of the Euler equations of gas dynamics (see Appendix A) might belong to the set of control variables.
Let ρ and ρφ be generic conservative variables.To ensure that the first stage of a sequential MCL algorithm [6,22] limits ∆ f ρ (i,j) using (23).The second stage limits ∆ f ρφ (i,j) using a discrete version of the product rule (ρφ) = ρ φ + φ ρ.The bar states of the φ variable are symmetric in the LGL-DGSEM version.This property is a further advantage compared to Bernstein-basis DG methods [22].
The inequality constraints to be enforced in the second stage are given by where ∆ĝ φ,Lim It is easy to verify that conditions (26) are equivalent to where ∆ f ρφ,Lim We use the bounding fluxes This definition, which is similar to (23), guarantees the validity of ( 28) and of the corresponding constraints for the flux-corrected bar state (ρφ) Lim (j,i) .The limited anti-diffusive flux ∆ f ρφ,Lim (i,j) is calculated using formula (29).

Pressure limiter
When solving the compressible Euler equations of gas dynamics (Appendix A), we require the pressure and internal energy to be non-negative at all times.Positivity preservation is guaranteed if the limited bar states satisfy To enforce (31), we apply a synchronized limiting factor α p (i,j) ∈ [0, 1] to all components of ∆ f Lim (i,j) .The limited bar states become and the prelimited anti-diffusive fluxes ∆ f Lim (i,j) are replaced with α p (i,j) ∆ f Lim (i,j) .Dropping the superscript p for better readability and introducing the scaled bar states w (i,j) := λ max (i,j) u (i,j) , we translate (31) into the quadratic inequalities where Following Kuzmin [6], we notice that α 2 ≤ α for α ∈ [0, 1].Therefore, (33) holds under the linear sufficient condition P (i,j) α ≤ Q (i,j) , where We conclude that the pressure fix can be performed using This definition exploits the property that the bar states of low-order LGL-DGSEM are symmetric.The general formula for α (i,j) is more involved [6].
To ensure continuous dependence of the limited fluxes α p (i,j) ∆ f(i,j) on the data, one may replace P (i,j) with the upper bound [6] We explore this possibility in the present paper.In the descriptions of our numerical experiments, we call the pressure limiter (35) that uses (34) "sharp".The one that uses (36) instead of (34) is referred to as "cautious".As we show in the Numerical Results section, the cautious pressure fix can add much more numerical dissipation to the scheme than its sharp counterpart.

Semi-discrete entropy limiter
In this work, we also use the semi-discrete entropy limiter developed in [57] for MCL schemes.Semi-discrete entropy stability of a DG or FV method is guaranteed if the numerical fluxes satisfy Tadmor's shuffle condition [62][63][64] v where Ψ (i,j) := Ψ j −Ψ i denotes the jump operator, Ψ is the so-called entropyflux potential, and v is the vector of entropy variables.See Appendix A for the definition of these quantities for the Euler equations.Let ∆ f Lim (i,j) be a limited anti-diffusive flux that is constrained to preserve local and/or global bounds for all scalar quantities of interest.Define using a correction factor α s (i,j) ∈ [0, 1] such that Tadmor's condition ( 37) is fulfilled.The substitution of ( 38) into (37) yields a linear inequality constraint for α s (i,j) ∈ [0, 1].We enforce this constraint using (cf.[57]) where is a small positive number and is the rate of entropy production before the application of α s (i,j) .
Remark 9 All correction tools described in this section lead to closed-form expressions for limited fluxes or correction factors.This distinguishes our approach from FCT/IDP alternatives that require solving nonlinear equations [8,37,39].

Definition of Bounds
We distinguish between global and local bounds.Global bounds enforce physical admissibility conditions, such as positivity of the density ρ and pressure p in the case of the compressible Euler equations of gas dynamics.Preservation of these bounds is a prerequisite for running challenging simulations without crashing.If a lower bound ρ min i ≥ 0 is used in the limiter for ρ and a (sharp or cautious) pressure limiter is applied in the final stage of the sequential limiting procedure, then the MCL scheme is positivity preserving in this sense.
The imposition of local bounds, on the other hand, makes it possible to avoid spurious oscillations within the global bounds and to improve the shock-capturing capabilities of the method.The corresponding numerical admissibility conditions are frequently formulated as local maximum or minimum principles.The inequality constraints of our MCL method are feasible if they are satisfied by the low-order bar states.Therefore, these states must be built into the definition of the upper and lower bounds.For instance, the value of a conservative or primitive quantity φ at node i may be treated as numerically admissible if it is bounded by φ (i,j) , (40) where N * (i) = {i−1, i+1} and φ i is the solution at node i of the previous time step.In the multidimensional case, the integer set N * (i) contains the indices of all nodes j = i such that the flux f(i,j) appears in the evolution equation for the nodal state u i .
We can rewrite (40) as where N (i) = {i − 1, i, i + 1} is the integer set containing i and the indices of all neighboring nodes of i, as φ (i,i) = φ i by definition (12).If the MCL bounds (41) are too tight, the flux-corrected scheme may fail to achieve the optimal order of accuracy in smooth regions.Wider bounds can be constructed by including the values of φ FV,n+1 at node i and its neighbors belonging to cells that physically contain the nodal point [14,15].The inclusion of extrapolated states makes it possible to guarantee linearity preservation on general meshes [6].Alternatively, local bounds can be relaxed using smoothness indicators [43] to avoid a potential loss of accuracy due to unnecessary limiting.
The use of the low-order bar states to define the bounds is also common in FCT/IDP methods [5,8,30,33].Similarly to the MCL version, further nodal states φ FV,n+1 j can be incorporated into the definition of φ min i and φ max i .The use of smoothness indicators is also an option [8,37].In FCT/IDP methods, it is also possible to define the bounds with the low-order solution at the next time step instead of using the low-order bar states, e.g.[37,39,65], This leads to more dissipative schemes, but the time-step restriction might be relaxed in some situations (see next section for details).A peculiarity of FCT/IDP schemes is the fact that the bounds for the limited anti-diffusive fluxes are inversely proportional to the time step.As a consequence, the quality of flux-limited approximations often exhibits strong dependence on the CFL number.
In this paper, we use the tight local bounds (41) for both MCL and FCT methods without any relaxation.

Explicit Time-Step Restrictions
The monolithic convex limiting method enforces nodal bounds by limiting the interface fluxes, such that all high-order bar states satisfy the bounds.This strategy enforces the bounds of the discrete solution if ( 20) is a convex sum.As a result, we obtain the following CFL-like time-step restriction in 1D: In general, the time-step restriction scales as where D is the number of spatial dimensions of the problem, m i is the diagonal mass matrix entry of node i, and λ max,d i is the maximum wave speed at node i in the coordinate direction d.For instance, with D = 1 we have λ max,1 i = max{λ max (i,i−1) , λ max (i,i+1) }.Since the CFL condition ( 44) is derived from the bar-states representation, MCL and FCT/IDP methods that use bar-state bounds (41) need to fulfill it to be able to keep the solution within bounds.Unfortunately, ( 44) is more restrictive than the typical CFL stability condition of the low-order method, especially for multiple space dimensions [66]: where ∆x d i is the size of the subcell in direction d.To the authors' knowledge, the only way to circumvent the strict CFL condition (44) and still keep the solution within prescribed bounds is to use FCT/IDP methods and avoid the bar-state bounds altogether.For instance, one can use bounds computed from the (robust) low-order solution at the next time step (42).Although this strategy can reduce the computing time in some situations, it might lead to non-physical solutions since the low order method is only provably positivity preserving (see, e.g., [67]) when the strict CFL condition (44) is fulfilled.
Since FCT/IDP methods apply limiting to enforce a fully discrete fulfillment of the bounds (i.e. after the time update is done), their anti-diffusion correction (and hence the spatial discretization) depends on the time-step size.On the other hand, MCL applies the limiting at the semi-discrete (spatial) level without considering the temporal discretization.As a result, the amount of dissipation depends on the CFL number for FCT/IDP methods, but not for the MCL approach.In fact, we have included a numerical investigation of the CFL dependence of the schemes and demonstrate that there is no obvious time convergence observable for FCT/IDP.

Numerical Results
To test the convergence and robustness properties of the MCL/DGSEM schemes, we run simulations with the compressible Euler equations of gas dynamics (see Appendix A).
In all cases, we use the Rusanov (LLF) numerical flux for the compatible robust low-order subcell scheme as well as for the surface fluxes in the highorder DG method.The CFL condition is given by (44).

Convergence Test
We simulate the advection of a density wave with initial condition (46) in the computational domain Ω = [−1, 1] 2 , and quantify the L 2 error of the solution as the mesh is refined.
We use the entropy-conserving and kinetic energy preserving flux of Ranocha [54] for the volume numerical flux of the split-form DGSEM method, choose a heat capacity ratio of γ = 1.4,and use CFL=0.9.
We first study the effect of imposing global bounds on the solution with the MCL approach.To do so, we impose strict positivity of density and pressure for all limited bar states, The positivity of density is enforced with a one-sided MCL limiter for conservative variables, and the positivity of pressure with the sharp pressure positivity limiter.Tables 1 and 2 show the L 2 error of all solution quantities and the experimental order of convergence (EOC) for the MCL/DGSEM scheme that imposes global bounds (positivity) for density and pressure.The EOC is computed for four different Cartesian meshes with N e elements per spatial direction.We observe an EOC ≈ N + 1 for approximations with polynomial degree N .
We now study the effect of imposing local bounds with the MCL approach.To do so, we impose local minima and maxima on the density, velocity and specific total energy of all bar states using the sequential MCL limiter, min    Tables 3 and 4 show the L 2 error of all solution quantities and the EOC for the MCL/DGSEM scheme that imposes local bounds (sequential limiting).In this case, the experimental order of convergence is at most second order, independent of the polynomial degree, which indicates that the local bounds are too strict to achieve an EOC > 2. Hence, what is done typically, one needs to relax the strict bounds to restore high-order accuracy, for instance by combining the local bounds of the MCL limiter with a smoothness sensor/indicator to avoid limiting smooth extrema of the approximate solution.

Kelvin-Helmholtz Instability
We consider the inviscid two-dimensional Kelvin-Helmholtz instability (KHI) setup, e.g., presented in [39,41].Due to its high density contrast and compressibility effects, the test case is challenging for nodal high-order methods when the (under-resolved) vortical structures of the KHI develop and evolve.In fact, the standard LGL-DGSEM method requires limiting to ensure robustness (positivity in this case).The initial condition is given by with B = tanh (15y + 7.5) − tanh(15y − 7.5).
We tessellate the simulation domain, Ω = [−1, 1] 2 , using 64 × 64 quadrilateral elements, use periodic boundary conditions, represent the solution with polynomials of degree N = 7, and run the simulation until the final time t = 10.Moreover, we discretize the Euler equations using the split-form DGSEM and the entropy-conserving and kinetic energy preserving flux of Ranocha [54] for the volume numerical fluxes, and select CFL=0.9 We first study the effect of imposing global bounds on the solution with the MCL and FCT/IDP approaches.With the MCL method, it is possible to impose strict positivity of density and pressure for all limited bar states (47) using a one-sided MCL limiter for conservative quantities and the sharp pressure positivity limiter.However, with the FCT/IDP approach presented in [39] a positive threshold greater than zero is necessary, as otherwise some nodes might get an invalid vacuum state.Hence, for the FCT/IDP variant, we consider the heuristic positivity-preserving method of Rueda-Ramírez [41] in a subcell-wise manner, i.e., we impose lower bounds for density and pressure that depend on the FV solution, with β = 0.1.We note, that this is a somewhat stricter requirement than strict positivity.
Figure 1 shows the density contours at different stages of the KHI simulation using the MCL and FCT/IDP limiters with global bounds -both approaches run stably until the final time.Even though this simulation setup is extremely sensitive to the discretization scheme [39] regarding the shape and form of the vortex roll-ups, the two approaches to impose global bounds produce remarkably similar looking results.
We now study the effect of imposing local bounds on the solution with the MCL and FCT/IDP approaches.With the MCL method, we use the standard sequential limiting to impose local minima and maxima on the density, velocity and total energy (48).With the FCT/IDP method, we impose local minima and maxima on the density, and local minima on the specific entropy, min  where the condition on the modified specific entropy, η = eρ 1−γ , guarantees the fulfillment of a discrete entropy inequality [13].Moreover, η is an efficient choice since it is computationally cheaper to evaluate than the specific entropy s = ln(p/ρ γ ) and it improves the convergence of the Newton method that is used in FCT/IDP methods to solve the non-linear equation to obtain the limiting factor [13,68].Note that condition ( 48) is typical for MCL and condition ( 51) is standard for FCT/IDP.Figure 2 shows the density contours at different stages of the KHI simulation using the MCL and FCT/IDP limiters with local bounds.Again, the solutions are very comparable between the two approaches, even though the limiting techniques and bounds are different.When comparing Figures 1 and 2, it is evident that the methods that impose local bounds add more numerical dissipation than the methods with global bounds, as expected.There is a smaller range of scales apparent with local bounds, especially at larger times.
Finally, we compute the mean limiting factor as where e ∈ [1, K] denotes the element index, K is the number of elements of the domain, i, j ∈ [0, N ] are the node indexes, N is the polynomial degree, α e ij (t) is the limiting factor of node ij of element e at time t, and V is the total area of the domain.Since in MCL methods the limiting is done for each interface and each equation without using a limiting factor, we first compute the effective limiting factor for each interface and each equation, and then compute a nodal α e ij using the average over all interfaces of node ij, We present a plot of the evolution of the mean limiting factors for the KHI simulations in Figure 3.To obtain limiting factors between 0 and 1, the factors for the momentum and energy equations of MCL are the effective scaling of the auxiliary flux, i.e., where is a very small number.We plot the quantity (1 − α) for FCT/IDP methods since the limiting factor of our FCT/IDP methods [39] is defined inversely as for MCL methods.A mean limiting factor ᾱ = 1 means that the discretization uses the unlimited high-order scheme everywhere, whereas a mean limiting factor ᾱ = 0 means that the anti-diffusion fluxes, ∆ f for FCT/IDP, or ∆ f ρ or ∆ĝ φ for MCL, are set to zero everywhere.

Inviscid Bow Shock Upstream of a Blunt Body
We consider the supersonic flow over a 2D blunt body that produces a detached bow shock to test the performance of the MCL/DGSEM method on curvilinear grids.This problem setup was proposed as an advanced test case for the Fifth International Workshop on High-Order CFD Methods [69].
The left boundary of the domain is a circular segment with origin (3.85, 0) and radius 5.9, the blunt body has a flat front of length 1 connected with two quarter circles of radius 0.5, and the right boundary is located at x = 0.The heat capacity ratio is set to γ = 1.4 and the initial condition is the constant state ρ(x, y) = 1.4,p(x, y) = 1, which corresponds to a Mach number Ma = 1.4.
For the blunt body we use a reflecting wall boundary condition, while for the other boundaries we use characteristics-based inflow/outflow boundaries, on which the external state is selected depending on the flow conditions normal to the boundary.
We use the split-form DGSEM with the entropy-conserving and kinetic energy preserving flux of Ranocha [54], polynomial degree N = 5, an isoparametric mapping of the geometry, MCL limiting with global (47) and local (48) bounds, and a conforming mesh with 36 elements distributed regularly on the inflow and wall boundaries and 24 elements distributed regularly on the outflow boundaries.To impose positivity of pressure, we use again the sharp pressure positivity limiter.
Figure 4 shows the pressure contours for the bow shock simulations at time t = 10 using MCL limiters with global and local bounds.Figure 4a shows that global positivity bounds are enough to keep the simulation running until the end time, but spurious oscillations appear near the shock.The use of local bounds removes the spurious oscillations near the shock, as can be observed in Figure 4b.

Sedov Blast Explosion
The Sedov blast problem is a very challenging simulation setup with a strong circular shock that describes the evolution of a symmetrical blast wave expanding from an initial concentration of density and pressure into a gas at rest.For the initial condition, we use the standard setup from the FLASH astrophysical code [70].The gas is initially at rest, v 1 (t = 0) = v 2 (t = 0) = 0, the density is constant ρ(t = 0) = 1, the atmospheric pressure is p 0 = 10 −5 , and we insert a quantity of dimensionless energy e = 1 into a small region of radius r 0 = 0.21875 at the center of the grid, with r = x 2 + y 2 .We tessellate the simulation domain, Ω = [−2, 2] 2 , with 64 × 64 quadrilateral elements, use periodic boundary conditions (however the final time is small enough, such that this does not matter), run the simulations with the split-form DGSEM and the entropy-conserving and kinetic energy preserving We test different variants of the MCL limiter with local bounds to identify which variants of the MCL limiter can handle the strong shocks of this test, and visualize the nodal limiting factors computed with (53).
The first variant of the MCL limiter that we test (from now on referred to as MCL limiter A) uses the standard MCL sequential limiter (48) first, then the sharp positivity limiter with global bounds (47), and the semi-discrete entropy limiter (38) at last.
Figure 5 shows the density contours and limiting factors for the Sedov blast simulation using the MCL limiter A at the final time t = 3.The text on the limiting factor plots indicates to which flux the factors are applied.Since both the pressure positivity limiter and the semi-discrete entropy limiter act on all components of the anti-diffusive flux, ∆ f , we indicate in brackets if the limiting factor corresponds to the pressure (p) or the entropy (ds/dt) limiters.Some artifacts at the shock front of the blast wave can be observed with this standard version of the MCL limiter.
The second variant of the MCL limiter that we test (from now on referred to as MCL limiter B) uses the MCL limiter for conservative quantities with local bar-state bounds for the density, min then computes the effective limiting factor for density, where is a very small number, applies α ρ (i,j) to all other conservative quantities (∆ f ρv1 , ∆ f ρv2 , and ∆ f ρE ), then uses the sharp positivity limiter with global bounds (47), and the semi-discrete entropy limiter (38) at last.
Figure 6 shows the density contours and limiting factors for the Sedov blast simulation using the MCL limiter B at the final time t = 3.In this case, the density, pressure, and the semi-discrete entropy limiter act on all components of the anti-diffusive flux.Note that the pressure positivity limiter needs to act less, but the semi-discrete entropy limiter needs to act more than in the MCL limiter A. Moreover, the MCL limiter B removes the Carbuncle-like artefacts at the shock fronts of the blast.
The last variant of the MCL limiter that we test (from now on referred to as MCL limiter C) is a combination of MCL limiters A and B. It uses the MCL limiter for conservative quantities with local bar-state bounds for the density (57), computes the effective limiting coefficient using (58) and applies it to the other conservative quantities, then uses the sequential limiter to impose the bounds on "primitive" quantities (48), then imposes global bounds on density and pressure (47) with the sharp positivity limiter, and then applies the semi-discrete entropy limiter (38) at last. of the anti-diffusive flux.Note that the sequential limiter for the velocities and the total specific energy needs to apply less limiting than in MCL limiter A due to the action of the density limiting coefficient on the momentum and total energy fluxes.The pressure positivity limiter needs to apply less limiting than in MCL limiters A and B, and the semi-discrete entropy limiter does not need to act at all for the snapshot at t = 3.As with MCL limiter B, the Carbunclelike artefacts are no longer present, but the resulting scheme is clearly more dissipative.

High-Mach Astrophysical Jet
To test the robustness of the MCL techniques, we simulate a setup inspired by an astrophysical jet application with Mach number Ma ≈ 2000, which was originally proposed by Ha et al. [71].This extreme benchmark case has been used to stress-test shock-capturing techniques for high-order methods [11,39,72].The computational domain, Ω = [−0.5,0.5] 2 , is filled with a mono-atomic gas (γ = 5/3) at rest with ρ(x, y) = 0.5, p(x, y) = 0.4127,  We solve this problem using 256 × 256 quadrilateral elements of degree N = 3, use periodic boundary conditions for the top and bottom boundaries and characteristics-based inflow/outflow boundary conditions for the left and right boundaries, the entropy-conserving and kinetic energy preserving flux of Ranocha [54] for the volume numerical flux of the DGSEM method, and different MCL and FCT/IDP limiters and CFL numbers.
We first compare the MCL limiter variant C from the previous section with the FCT/IDP method with local bar-state bounds for the density and specific entropy (51) at different CFL numbers.Figure 8 shows that the amount of vortical structures in the density contours at t = 10 −3 obtained with the FCT/IDP method is highly dependent on the CFL number, whereas the simulations that use the MCL limiter show a weaker dependence on the CFL number.Table 5 further shows that the total number of time steps needed to reach t = 10 −3 depends linearly on the CFL number for MCL methods, but not for FCT/IDP methods.
The dependence of the spatial discretization on the time-step size for FCT/IDP methods causes the number of vortical structures to be highly dependent on the CFL number and the total number of time steps to be not inversely proportional to the CFL number.In fact, the amount of dissipation is reduced for small CFL numbers, which leads to lower minimum densities and higher maximum pressures in FCT/IDP, as can be seen in Table 5. Something like a feedback effect occurs when the lower densities and higher pressures increase the speed of sound in the medium, which in turn reduces the time-step size even more, which again reduces dissipation (indirectly, due to lowering the bounds of the a-posteriori limiting approach in FCT/DIP).While having reduced dissipation is in general of course desirable, here it is more subtle as one buys the low dissipation with a strongly increased number of time steps, i.e., with a strongly increased CPU time.
Finally, we compare the difference between the sharp pressure positivity limiter ((35) with (34)) and the cautious pressure positivity limiter ((35) with (36)).Figure 9 shows the density contours and pressure limiter limiting factors at t = 10 −3 obtained with the MCL limiters A and C, and the cautious and sharp pressure positivity limiters.For this particular case, it is clear that the cautious pressure positivity limiter adds significant dissipation in the shear layer of the yet, which suppresses the appearance of vortical structures for both MCL limiters.On the other hand, the sharp pressure positivity limiter adds just enough dissipation to maintain the pressure of all bar states non-negative, and hence allows the development of turbulence.

Conclusion
In this paper, we have extended the monolithic convex limiting method to nodal discontinuous Galerkin methods (DGSEM) that use Legendre-Gauss-Lobatto (LGL) points.We have shown that the collocated nature of the LGL-DGSEM approximation greatly simplifies the design of MCL limiters and the derivation of compatible low-order invariant domain preserving schemes, which are needed to apply the MCL strategy.We have demonstrated the versatility of LGL-DGSEM/MCL methods to solve challenging simulation setups featuring supersonic, hypersonic and turbulent flow regimes.We have compared the performance of MCL methods and predictor-corrector-type flux corrected transport (FCT) subcell limiting methods.Unlike FCT methods, the amount of dissipation (and hence the spatial discretization) obtained with MCL techniques does not depend on the time-step size.As a result, time convergence can be expected for MCL but not for FCT schemes in problems that need stabilization.

Supplementary information. .
The methods used in this paper were implemented in the opensource high-order DG code Trixi.jl[73][74][75].We refer the interested reader to our reproducibility repository (https://github.com/amrueda/paper2023 MCL LGL-DGSEM), where we provide detailed instructions of how to duce the numerical results that we present here.

Fig. 1 :
Fig. 1: Density contours for the Kelvin-Helmholtz simulations with MCL and FCT/IDP (positivity limiting for density and pressure).DGSEM results with polynomial degree N = 7 and 64 × 64 elements.
(a) Global bounds (b) Local bounds

Fig. 4 :
Fig. 4: Pressure contours for the bow shock simulations with MCL at time t = 10.DGSEM results with polynomial degree N = 5 and 24 × 36 elements.Left plot shows the simulation with global bounds only (positivity) and right plot shows the result with local bounds (sequential limiting).

Fig. 5 :
Fig. 5: Contours of the density and limiting factors for the DGSEM simulation with MCL limiter A (sequential limiting, pressure and entropy limiters) of the Sedov blast test at t = 3.The text on the limiting factor plots indicates to which flux the factors are applied.

Fig. 6 :
Fig. 6: Contours of the density and limiting factors for the DGSEM simulation with MCL limiter B (density factor for all equations, pressure and entropy limiters) of the Sedov blast test at t = 3.The text on the limiting factor plots indicates to which flux the factors are applied.

Figure 7 Fig. 7 :
Fig. 7: Contours of the density and limiting factors for the DGSEM simulation with MCL limiter C (density factor for all equations, sequential limiting, pressure and entropy limiters) of the Sedov blast test at t = 3.The text on the limiting factor plots indicates to which flux the factors are applied.
(x, y) = 0, and on the left boundary there is a hypersonic inflow with ρ(−0.5, y B ) = 5, p(−0.5, y B ) = 0.4127, v 1 (−0.5, y B ) = 800, v 2 (−0.5, y B ) = 0, for y B ∈ [−0.05, 0.05], which corresponds to a Mach number of Ma = 2156.91with respect to the speed of sound of the jet gas, and Ma = 682.08 with respect to the speed of sound of the ambient gas.

Fig. 8 :
Fig. 8: Density contours for the astrophysical jet simulations with MCL limiter C (density factor for all equations, sequential limiting, pressure and entropy limiters) and FCT/IDP (density and specific entropy) strategies.DGSEM results with polynomial degree N = 3 and 256 × 256 elements.
(a) MCL C cautious (b) MCL C sharp (c) MCL A cautious (d) MCL A sharp

Fig. 9 :
Fig. 9: Density contours for the astrophysical jet simulations with MCL limiters A and C using the cautious and sharp pressure limiters.DGSEM results with polynomial degree N = 3 and 256 × 256 elements.CFL=0.9.

Table 1 :
L 2 errors and EOC for the convergence test with positivity limiting (global bounds), N = 3.

Table 2 :
L 2 errors and EOC for the convergence test with positivity limiting (global bounds), N = 4.

Table 3 :
L 2 errors and EOC for the convergence test with sequential limiting (local bounds), N = 3.

Table 4 :
L 2 errors and EOC for the convergence test with sequential limiting (local bounds), N = 4.

Table 5 :
Total number of time steps, density and pressure range at t = 10 −3 for the MCL and FCT/IDP simulations of the astrophysical jet as a function of the CFL number