A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations of the Euler Equations
AA Subcell Finite Volume Positivity-Preserving Limiterfor DGSEM Discretizations of the Euler Equations
Andr´es M. Rueda-Ram´ırez *† Gregor J. Gassner ∗ Abstract
In this paper, we present a positivity-preserving limiter for nodal Discontinuous Galerkindisctretizations of the compressible Euler equations. We use a Legendre-Gauss-Lobatto(LGL) Discontinuous Galerkin Spectral Element Method (DGSEM) and blend it locallywith a consistent LGL-subcell Finite Volume (FV) discretization using a hybrid FV/DGSEMscheme that was recently proposed for entropy stable shock capturing. We show that ourstrategy is able to ensure robust simulations with positive density and pressure when us-ing the standard and the split-form DGSEM. Furthermore, we show the applicability ofour FV positivity limiter in extremely under-resolved vortex dominated simulations and inproblems with shocks.
Keywords : Positivity Preservation, Discontinuous Galerkin Methods, Compressible Flow
Discontinuous Galerkin (DG) methods provide high-order accuracy and local conservation onarbitrary grids when dealing with advection-dominated problems. The combination of high-order DG methods with high-order explicit Runge-Kutta schemes for the time integration leadsto a highly accurate and easily parallelizable family of methods, e.g., [3, 8]. However, it hasbeen observed that in extremely under-resolved fluid flow simulations, or in the presence ofshocks, the oscillations of the high-order polynomials can lead to non-physical negative valuesof density or pressure and cause the scheme to crash.A particular class of DG methods, known as Discontinuous Galerkin Spectral ElementMethods (DGSEM), can be combined with a split-form representation of the governing equa-tions [6] or entropy-stable fluxes [1] to obtain non-linear stability and hence increased robust-ness. Nevertheless, these schemes still suffer from positivity issues in the presence of shocks[7, 12] and under-resolved turbulence.An effective and rather commonly used strategy to preserve physical (positive) values ofdensity and pressure in high-order DG methods was proposed by Zhang and Shu [14]. Thismethod shrinks the high-order numerical solution in each element around its mean value untila valid solution is obtained. Although this method is often very effective, we observed that itsometimes causes large jumps in the numerical solution, which may trigger instabilities andlead to restrictions in the time-step size. * Department of Mathematics and Computer Science, University of Cologne, Weyertal 86-90, 50931 Cologne,Germany. † Corresponding author: Andr´es M. Rueda-Ram´ırez ([email protected]). a r X i v : . [ m a t h . NA ] F e b n this paper, we propose a novel positivity-preserving limiter for DGSEM discretizationsthat is based on a subcell FV discretization. We make use of the hybrid FV/DGSEM scheme thatwas recently developed by Hennemann et al. [7] for shock capturing. The hybrid FV/DGSEMscheme blends a Legendre-Gauss-Lobatto (LGL) DGSEM discretization with a compatibleLGL-subcell Finite Volume (FV) scheme in an element-wise manner. In the present work,we present a simple strategy to compute the blending coefficient of each element at the end ofevery Runge-Kutta stage, such that the positivity condition is never violated. The Euler equations can be written as a system of conservation laws ∂ u ∂t + ∇ · ↔ f ( u ) = , (1)where the state quantities are density, momentum and total energy, u = [ ρ, ρ(cid:126)v, ρE ] T , and theflux reads ↔ f ( u ) = ρ(cid:126)vρ(cid:126)v ⊗ (cid:126)v + Ip(cid:126)v ( ρE + p ) , (2)where I is the identity matrix, and the pressure is computed with the calorically perfect gasassumption, p = ( γ − ρe, γ is the heat capacity ratio, and e = E − (cid:107) (cid:126)v (cid:107) / is the internalenergy. We make use of the blended scheme by Hennemann et al. [7] that combines a high-orderDGSEM discretization with a first-order FV method on the corresponding LGL-subcell gridto obtain an accurate scheme with robust shock-capturing properties. As both schemes areconstructed on the same subcell LGL distribution, they are directly compatible and thus thesemi-discrete version of (1) is given by ˙ u = (1 − α ) ˙ u DG + α ˙ u FV , (3)where ˙ u DG is the high-order DG discretization of the PDE, ˙ u FV is the first-order FV discretiza-tion, and α ∈ [0 , is an element-local blending coefficient that depends on a shock indicator.In the present work, we are not necessarily interested in using this blending approach for shockcapturing, but, as we will show, it can be used for positivity preservation.To obtain the DGSEM-discretization of (1), the simulation domain is tessellated into K non-overlapping elements and all variables are approximated within each element by piece-wise Lagrange interpolating polynomials of degree N on LGL nodes. These polynomials arecontinuous in each element and may be discontinuous at the element interfaces. Furthermore,(1) is multiplied by an arbitrary polynomial (test function) of degree N and numerically inte-grated by parts inside each element of the mesh with an LGL quadrature rule of N + 1 pointson a reference element, ξ ∈ [ − , , to obtain ˙ u DG j = F Vol j + δ jN J ω j (cid:16) f N − ˆ f ( N,R ) (cid:17) − δ j J ω j (cid:16) f − ˆ f ( L, (cid:17) , j = 0 , ..., N. (4)2n (4), ω j is the reference-space quadrature weight, δ ij is the Kronecker delta, J is the geom-etry mapping Jacobian from reference space to physical space, ˆ f ( i,j ) = ˆ f ( u i , u j ) is the surfacenumerical flux, which we evaluate with the left and right outer states, L and R , to account forthe jump of the solution across cell interfaces, and F Vol is the volume integral term.From (4), we can recover the standard and the split-form (flux differencing) DGSEM dis-cretizations by choosing the volume integral term as F Vol , Std j = − J N (cid:88) i =0 D ji f i , or F Vol , Split j = − J N (cid:88) i =0 D ji f ∗ ( j,i ) , (5)respectively, where D ji = (cid:96) (cid:48) i ( ξ j ) is the so-called derivative matrix, defined in terms of theLagrange interpolating polynomials, { (cid:96) i } Ni =0 , and f ∗ ( j,i ) = f ∗ ( u j , u i ) is a volume numerical two-point numerical flux. We can construct a provably entropy-stable (ES) DGSEM scheme bychoosing an entropy conserving flux for f ∗ , e.g. [2], and an entropy stable flux for ˆ f , e.g. [6].The compatible LGL co-located first-order FV discretization proposed by Hennemann et al.[7], which interprets the nodal values as the subcell mean values, reads ˙ u FV j = 1 J ω j (cid:16) ˆ f ( j,j − − ˆ f ( j,j +1) (cid:17) , (6)where the numerical fluxes on the element boundaries match with the DGSEM boundary val-ues and are evaluated with the left and right outer states, ˆ f (0 , − := ˆ f ( u , u L ) , ˆ f ( N,N +1) :=ˆ f ( u N , u R ) . We obtain the evolution of the state variables over time by applying a high-order explicit strongstability-preserving Runge-Kutta (SSP-RK) method to (3).In the SSP-RK methods of Spiteri and Ruuth [13], the solution is updated in every RK stage, s , with u s +1 = a ss u s + ∆ t s ˙ u s + s − (cid:88) i =1 (cid:0) a si u i + ∆ tb si ˙ u i (cid:1) , (7)where a si and b si are the RK coefficients, ∆ t s = b ss ∆ t is the so-called time-step size of the s RK stage.
First-order FV schemes are positivity preserving in density and pressure for the Euler equationsif a forward Euler time-discretization is used and an appropriate Riemann solver is selected[11]. Therefore, if α = 1 , the solution that is obtained is safe (positive in density and pressure), u s +1safe = a ss u s + ∆ t s ˙ u s, FV + s − (cid:88) i =1 (cid:0) a si u i + ∆ tb si ˙ u i (cid:1) . (8)As a result, it is possible to modify α locally, such that the discrete solution remains positivein density and pressure throughout the simulation. We propose a goal condition based on thecomputed safe solution ρ ≥ βρ safe , p ≥ βp safe , (9)3ith β ∈ (0 , . In other words, we allow only a certain deviation (controlled by the value β )below the safe solution, which is a stricter requirement than positivity. We do not consider upperbounds for density or pressure, as we are focusing on positivity. Note that upper bounds areequally straightforward to impose and that upper bounds might be necessary to avoid oscillatorysolutions. Let us obtain the corrected blending coefficient for density first. This is easy, because the densityis one of the conserved quantities, and as such, it depends linearly on the blending coefficient.Using (7) and (8), we get for every degree of freedom of an element, ρ safe − ρ = ∆ t s ( ˙ ρ s, FV − ˙ ρ s ) . (10)Summing and subtracting βρ safe we obtain ρ safe − ρ new + a ρ = ∆ t s ( ˙ ρ s, FV − ˙ ρ s ) , (11)where ρ new = βρ safe and a ρ = βρ safe − ρ .If a ρ > , requirement (9) is not fulfilled and α must be corrected such that ρ = ρ new . Werewrite (11) as ρ safe − ρ new = ∆ t s ( ˙ ρ s, FV − ˙ ρ s, new ) , (12)where ˙ ρ s, new = ˙ ρ s + a ρ ∆ t s , such that we can compute a new blending coefficient a posteriori toensure the fulfillment of the density criterion, α new = α + a ρ ∆ t s ( ˙ ρ s, FV − ˙ ρ s, DG ) , (13)where α is the blending coefficient that was used to obtain ˙ u .The corrected solution and its time derivative can be then computed with ∆ α = α new − α as u s +1new = u s +1 + ∆ α ∆ t s ( ˙ u FV − ˙ u DG ) , ˙ u s +1new = ˙ u s +1 + ∆ α ( ˙ u FV − ˙ u DG ) . (14)Note that, since we use element-local blending coefficients, α new must be taken as the max-imum of the corrected blending coefficients computed for all the degrees of freedom of anelement. After correction of density, we check the resulting approximation for pressure. If condition (9)is not fulfilled for the pressure, one must solve the following nonlinear equation for α new , g ( α new ) = p ( u new ( α new )) − βp safe = 0 , (15)which can be done using an iterative Newton’s method, α n +1new = α n − g ( α n ) ∂p ( α n ) /∂α . (16)where ∂p/∂α can be easily obtained using the chain rule, ∂p∂α = ∂p∂ u ∂ u ∂α . (17)4he first term in the right-hand side is obtained from (3) and (7), and the second term canbe derived from the equation of state, ∂ u ∂α = ∆ t s ( ˙ u FV − ˙ u DG ) , ∂p∂ u = ( γ − (cid:18) (cid:107) (cid:126)v (cid:107) , − (cid:126)v, (cid:19) . (18)We finally compute the corrected solution and its time derivative using (14) and the newelement-wise blending coefficient, α new , which is again taken as the maximum of the correctedblending coefficients for all the degrees of freedom of an element.In our proposed a posteriori positivity-preserving limiter, we first check density values andcorrect invalid states, and then check pressure values and correct if needed. The reason is that, inmany cases, performing the density correction removes invalid values of pressure automatically.Furthermore, the density correction is computationally cheaper than the pressure correctionsince the former can be computed with a linear, single-step process, while the latter needs aniterative method. For the examples of this section, we use the 2-register SSP-RK of fourth order and five stagesof Spiteri and Ruuth [13], where the time step is computed with the typical CFL condition forthe DGSEM [9] with CFL = 0 . . To control the positivity conditions (9), we choose the value β = 0 . . To test the capabilities of the FV positivity limiter, we first apply it to simulate an invis-cid Kelvin-Helmholtz instability (KHI). This setup is quite challenging for a high-order DGmethod, as it is severely under-resolved with an effective Reynolds number Re = ∞ . Theinitial condition is given by ρ ( t = 0) = 12 + 34 B, p ( t = 0) = 1 ,v ( t = 0) = 12 ( B − , v ( t = 0) = 110 sin(2 πx ) , (19)with B = tanh (15 y + 7 . − tanh(15 y − . . The simulation domain,
Ω = [ − , , is complemented with periodic boundary conditions.For this example, we tessellate Ω using K = 64 quadrilateral elements, represent the solutionwith polynomials of degree N = 3 , and run the simulation until t = 25 .The initial condition, (19), provides a Mach number Ma ≤ . , which makes compressibilityeffects relevant, but does not cause shocks to develop in the simulation domain. Thus, we donot apply shock capturing by connecting α to a troubled cell indicator, e.g. [7], but instead justuse α = 0 as our baseline high-order discretization.The Euler equations are discretized using the standard DGSEM and the entropy-stableDGSEM (ES-DGSEM). In both frameworks, we employ a traditional Rusanov scheme for thesurface numerical fluxes, ˆ f , and for the ES-DGSEM we use the entropy-conserving and kineticenergy preserving flux of Chandrashekar [2] for f ∗ .Figure 1 shows the evolution of the total mathematical entropy over time in the entire sim-ulation domain, S Ω ( t ) = − (cid:90) Ω ρsγ − x, (20)5here s = ln ( pρ − γ ) is the thermodynamic entropy. Note that the conventional use of theminus sign in (20) implies that S Ω decreases in the simulation domain if the second law ofthermodynamics is fulfilled.As can be observed in Figure 1, the ES-DGSEM fulfills the second law of thermodynamics,but the standard DGSEM does not. Note, however, that entropy stability is only guaranteedfor positive solutions and that the ES-DGSEM has no inbuilt positivity preservation. It can beobserved (close up plot (b) in the figure) that both schemes indeed crash at the beginning of thesimulation due to positivity problems. However, with the proposed positivity limiter strategy,both schemes are able to run until the end of the simulation. t S St. DGSEM + posit.ES-DGSEM + posit.St. DGSEMES-DGSEM (a) Entire simulation. t S St. DGSEM + posit.ES-DGSEM + posit.St. DGSEMES-DGSEM (b) Beginning of the simulation.
Figure 1: Total entropy evolution over time (a) and detailed zoom at the beginning of the simu-lation (b) for the Kelvin-Helmholtz instability simulation.Figure 2 illustrates the evolution of the blending coefficient over time in the whole domainusing a sample time ∆ τ = 0 . . The quantities max( α ) and ¯ α are defined as max( α )( t ) = max τ ∈ [ t − ∆ τ,t ] (cid:16) K max k =1 α k ( τ ) (cid:17) , ¯ α ( t ) = 1 n s n s (cid:88) s =1 (cid:32) K K (cid:88) k =1 α sk (cid:33) , (21)where k denotes the element index, α sk is the blending coefficient of element k at the RK stage s ,and n s is the number of RK stages taken from t − ∆ τ to t . A value ¯ α = 1 means that the entiredomain uses a first-order FV method, whereas ¯ α = 0 means that it uses only the high-order DGmethod.In Figure 2(a) it is evident that only a small amount of first-order dissipation is needed tostabilize the simulation, as the maximum blending coefficient for any element of the domain isalways below α ≤ . for the standard DGSEM, and α ≤ . for the ES-DGSEM.The FV stabilization of the high-order schemes is added very locally. As can be observedin Figure 2(b), the amount of the domain that is discretized with a first-order FV method is lessthan . during the entire simulation when the base-line scheme is the standard DGSEM,and less than . when the base-line scheme is the ES-DGSEM. Moreover, in the latter case,the entire domain is discretized with a pure ES-DGSEM method during most of the simulation.Finally, Figure 3 shows the density and the blending coefficient of both stabilized high-orderschemes for a snapshot at t = 3 . . Since the blending occurs sporadically and in small amounts,it is very difficult to capture a blending coefficient α > at a given snapshot. Therefore, weplot the maximum blending coefficient of each element of the domain over a sampling time ∆ τ = 0 . . It can be observed in Figure 3 that the positivity limiter preserves the symmetriesof the problem, and that the standard DGSEM needs more stabilization and is more distortedthan the ES-DGSEM. 6 t m a x () St. DGSEM + posit.ES-DGSEM + posit. (a) Maximum blending coefficient. t (b) Amount of low order method used in the do-main. Figure 2: Evolution of the blending coefficient over time for the Kelvin-Helmholtz instabilitysimulation.
The Sedov blast problem describes the evolution of a radially symmetrical blast wave that ex-pands from an initial concentration of density and pressure into a homogeneous medium. Forthe initial condition, we assume a gas in rest, v ( t = 0) = v ( t = 0) = 0 , with a Gaussiandistribution of density and pressure, ρ ( t = 0) = ρ + 14 πσ ρ exp (cid:18) − r σ ρ (cid:19) p ( t = 0) = p + γ − πσ p exp (cid:18) − r σ p (cid:19) , (22)where we choose σ ρ = 0 . and σ p = 0 . . Furthermore, the ambient density is set to ρ = 1 and the ambient pressure to p = 10 − .We complement the simulation domain, Ω = [ − . , . , with periodic boundary con-ditions. Furthermore, we tessellate Ω using K = 64 quadrilateral elements, represent thesolution with polynomials of degree N = 7 , and run the simulation until t = 20 .The smooth initial condition defined by (22) quickly evolves into an expanding shock frontthat bounces back into the domain because of the periodic boundary conditions. This generatesa complex shock pattern with multiple shock-shock interactions, which causes the transition toturbulence (Re = ∞ ).To deal with the strong shocks of this simulation, we use the indicator proposed by Henne-mann et al. [7], which is a modification of [10]. We use the gas pressure for the indicator, suchthat α reads α = 11 + exp (cid:0) − s T ( E − T ) (cid:1) , (23)where s = 9 . is the so-called sharpness, T ( N ) = 0 . · − . N +1) . is the threshold, and E is the total energy of the two highest modes of the pressure solution. Furthermore, we performone spatial propagation sweep of the blending coefficient to avoid large jumps in the blendingcoefficient from element to element. The blending coefficient of element k is adjusted with α k = max E { α, . α E } , where α E denotes the blending coefficient of any neighbor element to k . In this example, we use the HLLE scheme [4] for the surface numerical fluxes of bothframeworks, ˆ f , which has been shown to be positivity preserving [5]. For the ES-DGSEM weuse the entropy-conserving and kinetic energy preserving flux of Chandrashekar [2] for f ∗ .7 a) St. DGSEM - Density. (b) St. DGSEM - Blending coefficient.(c) ES-DGSEM - Density. (d) ES-DGSEM - Blending coefficient. Figure 3: Snapshot of the KHI simulations at t = 3 . for the standard DGSEM (top) and theES-DGSEM (bottom). The blending coefficient is taken as the maximum over all the SSP-RKstages of a sampling interval ∆ τ = 0 . .Figure 4 shows the evolution of the total mathematical entropy over time in the entire sim-ulation domain. It can be observed that, in spite of the use of the shock indicator, both thestandard and the ES-DGSEM need positivity-preserving stabilization. The standard DGSEMwithout FV stabilization crashes due to positivity problems when the shocks start forming at t ≈ . . The ES-DGSEM simulation without stabilization is able to run longer, but finallycrashes at t ≈ . . However, both the standard and the ES-DGSEM are able to run until theend of the simulation with the proposed positivity limiter strategy.Figure 5 shows the evolution of the correction to the blending coefficient over time forthe entire simulation using a sample time ∆ τ = 0 . . Note that the illustrated quantities arecomputed with (21) using ∆ α = α new − α instead of α , such that we can appreciate the actionof the positivity limiter.As can be observed in Figure 5(a), the ES-DGSEM only needs a small amount of extrafirst-order dissipation to stabilize the last part of the simulation ( t > . ), where the blendingcoefficient correction is always below ∆ α ≤ . . On the other hand, the standard DGSEMneeds very high corrections of the blending coefficient at the beginning of the simulation, where ∆ α oscillates between . and . It is clear that the shock indicator, (23), was tuned for the ES-DGSEM. 8 .0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 t S St. DGSEM + posit.ES-DGSEM + posit.St. DGSEMES-DGSEM (a) Entire simulation. t S St. DGSEM + posit.ES-DGSEM + posit.St. DGSEMES-DGSEM (b) Beginning of the simulation.
Figure 4: Total entropy evolution over time (a) and detailed zoom at the beginning of the simu-lation (b) for the Sedov blast simulation. t m a x () St. DGSEM + posit.ES-DGSEM + posit. (a) Maximum blending coefficient. t (b) Amount of low order method used in the do-main. Figure 5: Evolution of the correction to the blending coefficient over time for the Sedov blastsimulation.Even though the standard DGSEM needs high amounts of stabilization in some elements,Figure 5(b) shows that the amount of the domain that is discretized with a first-order FV correc-tion is less than . at the beginning of the simulation, and less than . at the end of thesimulation when the base-line scheme is the standard DGSEM. Moreover, when the baselinescheme is the ES-DGSEM, the amount of first-order FV correction is less than . .In Figure 6, we plot log( ρ ) for different snapshots, as this quantity allows to appreciatethe different scales that appear in the problem. Note that the standard DGSEM solutions aremore distorted than the ones obtained with the ES-DGSEM, as was also observed in the KHIsimulation (Figure 3).Vortices develop in the domain due to the appearance of Kelvin-Helmholtz instabilities inthe density layer, which are in turn triggered by the shock-shock interactions. In fact, it is theshock-driven vortical dominated flow part that causes the ES-DGSEM without positivity limiterto crash (see Figures 4(a), 6(b) and 6(d)). On the other hand, Figures 6(a) and 6(c) confirm thatthe standard DGSEM without positivity limiter crashes when shocks are forming, a long timebefore the shock front bounces on the boundaries.Finally, Figure 7 shows the distribution of α for snapshots of the simulations at t = 2 and t = 12 . . A comparison of the blending coefficient contours with Figure 6 reveals that the9 a) St. DGSEM, t = 2 . (b) St. DGSEM, t = 12 . .(c) ES-DGSEM, t = 2 . (d) ES-DGSEM, t = 12 . . Figure 6: Snapshots of the density logarithm for the Sedov blast simulations at t = 2 (left) and t = 12 . (right) for the standard DGSEM (top) and the ES-DGSEM (bottom).indicator, (23), correctly detects the shocks and is not triggered on turbulent areas of the domain.The action of the positivity limiter is evident in Figure 7(a): the blending coefficient correctionaffects elements in the post-shock region of the standard DGSEM, where the spatial propagationof the blending coefficient had already acted. This behavior reveals that an appropriate shockindicator for the standard DGSEM should detect the post-shock regions and explains the jumpsof ∆ α between . and that we observed in Figure 5(a). We presented a positivity-preserving limiter for DGSEM discretizations of the compressibleEuler equations that is based on a subcell Finite Volume method. We show that our strategy isable to ensure robust simulations with positive density and pressure when using the standardand split-form DGSEM if a positivity preserving Riemann solver is used. We tested our schemewith the Rusanov and HLLE Riemann solvers in extremely under-resolved vortex dominatedsimulations and in problems with shocks.
Acknowledgments
This work has received funding from the European Research Councilthrough the ERC Starting Grant EXTREME, with the ERC grant agreement no. 714487.10 a) St. DGSEM, t = 2 . (b) St. DGSEM, t = 12 . .(c) ES-DGSEM, t = 2 . (d) ES-DGSEM, t = 12 . . Figure 7: Snapshots of the blending coefficient for the Sedov blast simulations at t = 2 (left)and t = 12 . (right) for the standard DGSEM (top) and the ES-DGSEM (bottom). References [1] M. H. Carpenter, T. C. Fisher, E. J. Nielsen, and S. H. Frankel. Entropy stable spectralcollocation schemes for the Navier-Stokes Equations: Discontinuous interfaces.
SIAMJournal on Scientific Computing , 36(5):B835–B867, 2014.[2] P. Chandrashekar. Kinetic energy preserving and entropy stable finite volume schemesfor compressible euler and Navier-Stokes equations.
Communications in ComputationalPhysics , 14(5):1252–1286, 2013.[3] B. Cockburn and C. W. Shu. Runge-Kutta Discontinuous Galerkin methods forconvection-dominated problems.
Journal of Scientific Computing , 16(3):173–261, 2001.[4] B. Einfeldt. On Godunov-type methods for gas dynamics.
SIAM Journal on NumericalAnalysis , 25(2):294–318, 1988.[5] B. Einfeldt, C. D. Munz, P. L. Roe, and B. Sj ¨ogreen. On Godunov-type methods near lowdensities.
Journal of Computational Physics , 92(2):273–295, 1991.[6] G. J. Gassner, A. R. Winters, and D. A. Kopriva. Split form nodal discontinuous Galerkinschemes with summation-by-parts property for the compressible Euler equations.
Journalof Computational Physics , 327:39–66, 2016.117] S. Hennemann, A. M. Rueda-Ram´ırez, F. J. Hindenlang, and G. J. Gassner. A provablyentropy stable subcell shock capturing approach for high order split form DG for the com-pressible Euler equations.
Journal of Computational Physics , page 109935, 2020.[8] F. Hindenlang, G. J. Gassner, C. Altmann, A. Beck, M. Staudenmaier, and C. D. Munz.Explicit discontinuous Galerkin methods for unsteady problems.
Computers and Fluids ,61:86–93, 2012.[9] N. Krais, A. Beck, T. Bolemann, H. Frank, D. Flad, G. Gassner, F. Hindenlang, M. Hoff-mann, T. Kuhn, M. Sonntag, and C. D. Munz. FLEXI: A high order discontinuousGalerkin framework for hyperbolic–parabolic conservation laws.
Computers and Mathe-matics with Applications , 81:186–219, oct 2021.[10] P.-O. Persson and J. Peraire. Sub-Cell Shock Capturing for Discontinuous Galerkin Meth-ods. , pages 1–13, 2006.[11] B. Perthame and C. W. Shu. On positivity preserving finite volume schemes for Eulerequations.
Numerische Mathematik , 73(1):119–130, 1996.[12] A. M. Rueda-Ram´ırez, S. Hennemann, F. J. Hindenlang, A. R. Winters, and G. Gassner.
AnEntropy Stable Nodal Discontinuous Galerkin Method for the resistive MHD Equations.Part II: Subcell Finite Volume Shock Capturing . 2020.[13] R. J. Spiteri and S. J. Ruuth. A new class of optimal high-order strong-stability-preservingtime discretization methods.
SIAM Journal on Numerical Analysis , 40(2):469–491, 2002.[14] X. Zhang and C. W. Shu. On positivity-preserving high order discontinuous Galerkinschemes for compressible Euler equations on rectangular meshes.