Stabilizing Discontinuous Galerkin Methods Using Dafermos’ Entropy Rate Criterion: II – Systems of Conservation Laws and Entropy Inequality Predictors

A novel approach for the stabilization of the Discontinuous Galerkin method based on the Dafermos entropy rate crition is presented. First, estimates for the maximal possible entropy dissipation rate of a weak solution are derived. Second, families of conservative Hilbert-Schmidt operators are identified to dissipate entropy. Steering these operators us-ing the bounds on the entropy dissipation results in high-order accurate shock-capturing DG schemes for the Euler equations, satisfying the entropy rate criterion and an entropy inequality.


Introduction
Discontinuous Galerkin methods [5] are a popular tool to design numerical schemes for hyperbolic systems of conservation laws [8] ∂ f (u) An intriguing feature of DG methods is their ability to transfer the definition of a weak solution to a hyperbolic conservation law [15] ∀ϕ ∈ D : ( to the semidiscrete level [2,3,4].Using a method of lines approach this leads to the set of equations A cell of the subdivision T of the domain Ω

T
The left and right boundaries of cell T T l , T r The set of test functions D The space of ansatz functions for an cell T V T = span{ϕ 1 , ϕ 2 , . . ., ϕ N } Polynomials of degree of p in cell T V T,p L 2 projection of u onto V P V u.
Interpolation of u on V w.r.t. the collocation points (ξ k ) N k=1 Vector of nodal values in cell T at time t u T (t) ansatz function in cell T at position x and time t u T (x,t) The total entropy in cell T E u,T (t) The discrete total entropy in cell T E T (t) The inner product on T discretised using The inner product between u and v on cell T ⟨u, v⟩ T The p norm of u in cell T ∥u∥ T,p Exact solution to the initial condition u(x,t 0 ) after t − t 0 H(u(•,t 0 ),t − t 0 ).
Mean value of subcell k of N, u T as initial condition u T,N k The convex hull of a set A ch A Table 1: Notation used.As a general rule, quantities with only t as an argument are vectors of nodal values at a certain time.Values with x and t in their argument list are functions that were evaluated at these values.Objects with T added as exponent are approximations of the quantity in the cell T .
for every cell T ∈ T of a subdivision T of the domain into cells.The solution u(x,t) is approximated in every cell by u T (x,t) ∈ V T out of a finite dimensional space of ansatz functions V T .Using an approximation of the inner products as point evaluations results in the matrix vector form Sadly, the constructed schemes lack robustness and stability in the high order case and some stabilization measures and robustness enhancements are needed.Popular are overintegration, flux-differencing, modal filtering, sub-cells and (W)ENO recoveries [3,12,14,27,29,34,42].In this publication the procedure first presented in [21] will be refined, connections to some other stabilization techniques will be shown, and the technique will be tested on a catalog of problems for the Euler system of conservation laws.The method in [21] is based on the entropy rate admissibility criterion [6,11].An entropy [24] is a convex functional U : R m → R satisfying d f du dU du = dF du in conjunction with a entropy flux function F : R m → R. One can show that for this pair (F,U) holds ∂U(u(x,t)) ∂t in the sense of distributions [24].If the solution is smooth one can even show The entropy rate criterion states that the total entropy E u (t) = U(u(x,t)) dx of the selected weak solution u should reduce faster than the entropy of any other existing weak solution ũ ∀t > 0 : A numerical approximation of this total entropy can be defined as via a (positive) quadrature rule ω T k on each cell T ∈ T .The numerical enforcement of the criterion with respect to such a definition of the discrete entropy happened in [21] in three steps • Calculate the time derivative of the ansatz function d ũT dt using a DG scheme • Calculate an error prediction δ T for d ũT dt on T , i.e. d ũT dt − ∂ u ∂t T ≤ δ T .
• Correct the time derivative into the direction of the steepest entropy descent where h shall be the steepest descent direction that does not change the average value in cell T .
While the approach above is successful for scalar conservation laws [21] significant improvements can be made by introducing two refinements.The first one concerns the usage of an error indicator to estimate the entropy correction needed.We will instead show that it is possible to directly give bounds on how dissipative a weak solution can be.This will eliminate the need for the error indicator while allowing a faster convergence, because the derived bounds converge to zero significantly faster in the smooth case.A second refinement concerns the direction used for the entropy correction.DG methods can make use of modal filtering to remove unwanted high frequency modes from the solution [29].These filters can be sometimes expressed as viscosity, and we will devise correction directions that at the same time dissipate entropy and filter the solution from unwanted oscillations and thereby combine the dissipation and filtering.Our schemes will therefore follow the slightly different general layout of the dissipation mandated by the estimate.Here F * shall be a numerical entropy flux [35,36,37].
The procedure makes only use of the fact that in our cells there exist local ansatz functions and is therefore also applicable to similar schemes like the spectral volume (SV) method [40].The only difference would lie in the evaluation of a different scheme for the uncorrected derivative d ũ dt .One complication is brought in by the fact that entropy dissipation implies the non-smoothness of the solution, as otherwise the entropy equality applies.Therefore, dissipation can't happen in cells in the continuous setting, as polynomials are smooth.Instead, dissipation is a process taking place at the cell edges were our different ansatz functions transition.As we are not correcting the numerical fluxes used between cells dissipation will be centered in cells and not at cell edges, and we will show in section 3.1 how to work around this problem.

Entropy inequality predictors 2.1 Bounds for entropy and entropy dissipation
Our main tool to approximate the most dissipative weak solution using a DG method will be a bound on the derivative of the total entropy.We will derive a lower bound for the entropy dissipation Here θ ⊂ Ω shall be an arbitrary open subdomain of the complete domain.This value has to be smaller than zero for a solution that is admissible with respect to the classical entropy inequality (4).Further, we are interested in the entropy dissipation speed If this value is known one can estimate the total entropy's E u (t) derivative as To achieve our goal of estimating s θ we will view the problem in the setting of classical Finite-Volueme schemes [33] and go over to the limit ∆x → 0. In [7] it was shown that for scalar conservation laws the flux f of the solution to the Riemann problem u R (u l , u r ; x,t) is given by i.e. by entering the value of u into f that when entered into the flux yields the fastest entropy dissipation.In [22] it was shown that some approximate Riemann solvers, for example the local Lax-Friedrichs flux, can be also interpreted as approximate solutions to such variational descriptions of two-point fluxes.While the aforementioned results hold for semidiscrete schemes the new results below are new and aim at three point first order Finite-Difference/Finite-Volume schemes for systems of conservation laws.As one assumes piecewise constant functions in those first order methods any quadrature exact for constants will yield the same result in equation ( 5).As we only look at discrete time values in this part of the publication we will write E n u = E u (t n ) for the discrete total entropy at time level n.
Lemma 1.Let a system of hyperbolic conservation laws in conservation form and a strictly convex entropy pair (U, F) be given that is approximated by a Finite-Volume scheme with grid constant λ = ∆t ∆x .Then the original Lax-Friedrichs scheme has the fastest dissipation of the total entropy under all consistent and conservative three-point numerical schemes.
Proof.Assume f (u k , u k+1 ) is a consistent numerical two point flux minimizing the total entropy with maximal rate and let u l , u r be arbitrary in the domain of admissible values for the conserved variables.We apply a scheme using this flux to a Riemann problem, i.e. the initial data by analyzing the solution.As the flux is consistent it holds The scheme , as U is strictly convex.Entering this into the scheme's definition with u 0 0 = u l and u 0 1 = u r implies Rearranging for f (u l , u r ) shows and this is the classical Lax-Friedrichs flux and therefore uniquely determined by demanding maximal entropy rate.
This result shows that the classical LF scheme is the most direct realisation of a scheme satisfying Dafermos' entropy rate criterion and therefore justifies the use of the LF scheme in [20] as the most dissipative scheme possible for systems of conservation laws.Similar results are also known for scalar conservation laws.Tadmor showed in [35,36] that every monotonicity preserving scheme satisfying classical numerical entropy inequalities for a scalar conservation law has a viscosity coefficient less or equal to that of the LF scheme, and higher or equal than the viscosity coefficient of Godunov's scheme.Our result can be seen as a generalization of the LF part of this result to systems of conservation laws, as it states that the LF flux is the most dissipative flux for a selected time-step size.Using the scheme above one can derive estimates for the highest possible entropy dissipation in a time-step and using finite-differencing of this result, approximations for the lowest possible derivative of the total entropy with respect to time.corollary 1.The biggest possible entropy dissipation during a discrete time-step of a Finite-Volume scheme with grid constant λ = ∆t ∆x is given by the difference and an approximation to the total entropy's minimal derivative by The second estimate above degenerates for λ → 0 as the difference in entropy is in general finite between cells u k−1 , u k , u k+1 .A second, more refined, estimate is given by the following lemma based on the ideas from [18] and does not have these deficiencies.0 t x a l t = x a r t = x (a r , 1) 1.0 (a l , 1) (a) The case a l < 0 < a r 0 t x a l t = x a r t = x (a r , 1) (a l , 1) (b) The case 0 < a l < a r and vice-versa.
Figure 1: Layout of the integration areas in the proof (blue).As originally used in [18].To the left and right of the lines dx dt = a l and dx dt = a r the initial condition is unaltered.
Lemma 2. Given bounds on the fastest signal speed to the left a l and the highest signal speed to the right a r let M ≥ max(|a l | , |a r |).The maximum entropy dissipation of a Riemann problem solution on the interval θ = (−M, M) is bounded from below by The rate is bounded from below by The entropy dissipation is bounded by and its rate by Proof.The entropy of the initial condition in the interval [−M, M] is given by for any M > 0. Integrating over the triangle T = ch{(0, 0), (a l , 1), (a r , 1)} in spacetime and using the conservation law yields in conjunction with the Gauß divergence theorem, cf.figure 1.Here u lr shall denote the mean value of u(x, 1) on [a l , a r ] and is as apparent from the calculation above.Jensens inequality implies Therefore it follows for the entropy dissipation between t = 0 and t = 1 and using the invariance under transformations (x,t) → (µx, µt) for µ > 0 yields for the rate.To calculate the entropy dissipation s θ and its speed σ θ we just have to account for the entropy flowing in and out of the intervall θ using the entropy flux F. This is possible as u is constant to the left of (ta l ,t) and to the right of (ta r ,t).
The estimate above does not depend on any grid constant, and reduces to the previous one for −a l = c max = a r , λ c max = 1, and this is the CFL condition for the classical Lax-Friedrich scheme, i.e. both estimates are compatible.A Godunov type scheme using the HLL approximate Riemann solver is also compatible with the estimate above.The discrete total entropy after one time-step is still less or equal than the bound given above.Let λ c max ≤ 1 2 hold implying that the Riemann problems do not interact and u HLL (x,t n+1 ) be the picewise constant solution of the HLL solver as in figure 2, but not averaged over the cells, while u n+1 k shall be the corresponding cell averages.In this case the total discrete entropy at the next time-step is given by Therefore the discrete entropy of the approximate solution is lower than the entropy of any exact weak solution.The next subsection will move beyond first order schemes by generalizing this lower bound to one that also allows smooth solutions instead of piecewise constant ones.x t a l a r (b) Solutions of generalized Riemann problems lack scaling invariance.Still we assume that there exist bounds a l and a r that the waves from the interaction of u l (x) and u r (x) do not leave the cone [ta l ,ta r ] for small t.
x u(x,t) ta l ta r (c) The assumed solution used in the generalized entropy inequality predictor.Note that this solution follows the HLL idea of assuming a constant function in the wedge formed by ta l and ta r .
x u(x,t)

Asymptotic analysis based entropy inequality predictor
The entropy inequality predictor in this section will be based on an asymptotic analysis of the problem described in figure 3a, i.e. two smooth solutions splined together at an interface.An obstacle lies in the missing self-similarity.This is a difference to the previous part where the self-similarity of the initial condition and assumed self-similarity of the solution induced the existence of a self-similar, i.e. constant, speed of the entropy dissipation.We will therefore try to approximate for reasonably small |t 2 − t 1 | and the discontinuity at the interface in the interior of θ .The schemes in which we will use these entropy inequality predictors should converge with high orders for smooth solutions, necessitating a convergence of the predictor to zero with a high order for smooth solutions.This convergence is also dictated by the entropy equality for smooth solutions.If u l (x) and u r (x) are piecewise constant this problem is already solved by the methods described in the last subsection.We will therefore now reiterate through the proof of lemma 2 assuming that u l (x) and u r (x) are smooth functions.The missing selfsimilarity of Generalized Riemann problems [1], cf.figure 3b, defies the existence of the speed estimates a l and a r , and we therefore just assume that these speed estimates exist for small times.Further we assume that for small times the solutions left of (ta l ,t) and right of (ta r ,t) remain smooth, as no waves from the interaction arrive there and u l (x), u r (x) have bounded derivatives.
The average value u lr shall be determined by applying the conservation law to the triangle T = ch{(0, 0), (ta l ,t), (ta r ,t)} Dividing this equation by t and going over to the limit t → 0 results in using the continuity of the integrands and the mean value theorem of integration [26].Therefore it follows for vanishing t.Equation ( 7) stays also valid in the case of piecewise polynomial functions as initial conditions and for small t > 0. We can therefore conclude that a generalization of equation ( 8) holds in the form Accounting for the entropy flowing in and out of [−M, M] yields Applying the entropy equality to the subdomains that holds for small t > 0 because the solution stays smooth in the subdomains, allows us to restate this as Dividing by t and going over to the limit, using the limit of u lr and once more the mean value theorem, shows in this case also A significant problem of the derivation above lies in the fact that one can only estimate the entropy dissipation speed in the interval θ = (−M, M), but not in (−M, 0) as the true dissipation can be located anywhere in the cone [ta l ,ta r ].As the cells in our numerical tests will be layed out as in figure 3d We are therefore left with the problem of how to split this dissipation onto the two neighboring cells that have overlap with θ k+ 1 2 . This problem will be handled below in section 3.1.

Accounting for aliasing errors
In [21] one of the findings in the numerical tests section was that the entropy dissipation of the numerical solutions started already shortly prior a real entropy dissipating discontinuity formed.This was attributed to the fact that while the entropy of the exact solution is still constant as long as the solution is smooth this exact solution will in general not be representable in our ansatz space.It is therefore wise to dissipate entropy to arrive at a function that still lies in our space, and certainly better than selecting an ansatz function that has more entropy than the true solution.A similar issue could be the fact that in the L p norms, for p < ∞, near each piecewise continuous solution u T lies a C ∞ function that can be constructed via mollification.Therefore an infinitely small perturbation of u T in the usual norms leads to a vanishing entropy dissipation.Or, put differently, the dissipation bound as a functional is discontinuous in the L p spaces.While unsatisfactory let us remark that the functional is better behaved with respect to the BV semi norms.The discontinuity of the entropy dissipation bound is problematic with under-resolved solutions where a lucky, or in this case better to be considered unlucky, too smooth approximation of the solution in our piecewise polynomial spaces induces wrong, i.e. too conservative entropy dissipation predictions.
We are therefore interested in allowing our entropy inequality predictor to be also greedy, or one could say pessimistic, with respect to an under-resolved solution.The key to this strengthening is the following lemma.
Lemma 3 (Order of the entropy dissipation bound).The maximal entropy dissipation prediction (10) of a Riemann problem for a smooth flux function with smooth entropy-entropy flux pair vanishes quadratically with the jump of u at the interface Proof.As the entropy inequality holds it is clear that the entropy dissipation is non-positive in the sense of distributions.As we only allow entropy dissipative solutions the entropy dissipation on θ is a non-positive constant for a fixed jump.On the contrary (10) has to be zero for u l = u r and is smooth, implying that the line u l = u r consists of local maxima.Therefore, a power expansion of (10) around u l in u r has to be of the form with a negative semi-definite Hessian H ∈ R m×m .This proofs the claim.
We are therefore in the relaxing position that even if our approximations of u l , u r only satisfy ∥u l − u r ∥ ∈ O((∆x) p ) the corresponding estimate will converge significantly faster with σ θ (u l , u r ) ∈ O((∆x) 2p ).Our basic DG method predicts values for our solution u T in a Hilbertspace that is spanned by polynomials on every cell.In this case a suitable orthonormal basis is spanned by Legendre polynomials and the limits of these basis representations are L 2 functions.But as explained before our functional is not continuous on L 2 and our ansatz u T is only an approximation of a projection of the true solution onto our ansatz space.We can therefore try to exploit different projections of our ansatz, especially projections that assume less regularity of u T , and estimate our entropy dissipation with the strongest one encountered in all of these different approximations of u l and u r .A natural choice for projections on spaces assuming less regularity are projections on lower order polynomials.As the Legendre polynomials on each cell when truncated up to polynomial p are an orthogonal basis of the polynomials with degree less than or equal to p, is the ∂t onto V .Our approximation du T dt is not entropy dissipative in this example and should be corrected into the entropy dissipative half-space characterized by the normal − dU du .The direction ∇ 2 u that stems from a discretization of the heat kernel is suitable for this correction and has additional benefits compared to − dU du .While the diffusion also has a smoothing effect the addition of − dU du can even result in a sharpening effect.Higher even derivatives like ∇ 8 u will smooth the solution but will not result in a dissipation for all entropies.orthogonal projection onto these spaces given by discarding the higher order coefficients in the Legendre expansion of u T .We can truncate down to order p − 1 by discarding the highest coefficient and still achieve a convergence order of at least q = 2p − 2 > p of our entropy inequality predictor for p > 2. This can be summed up in the following procedure used above order p = 2.

• Assign
• Project the ansatz u T in every cell onto V T,p−1 using an orthonormal projection u T,p−1 = P V T,p−1 u(•,t).
as entropy inequality prediction.

Suitable dissipation directions and filtering
After deriving approximations for the entropy dissipation needed we will now determine how to correct the time derivative of the DG scheme to dissipate the amount of entropy needed.At the same time the resulting scheme is hopefully still high order accurate for entropy conservative solutions.In the scalar case the direction of the steepest descent of the entropy, corrected for conservation, was used for this purpose.This approach incurs several problems: • The direction of steepest entropy descent has in general no smoothing/filtering effect.
• Proving in the previous publication that a correction in the steepest descent direction with the length taken from the error indicator results in enough entropy dissipation was possible but resulted in highly technical arguments [21].• Dissipation stems from the viscous and parabolic history of hyperbolic conservation laws.A viscous flux x associated with a viscous regularization of a hyperbolic conservation law is proportional to the gradient of the solution for fixed viscosity and proportional to the gradient of the solution for constant viscosity.If a component of u is smooth with a low magnitude of the first and second derivative the viscous flux of this component will also only differ from the hyperbolic flux by a small margin.If our scheme is corrected with the steepest entropy descent direction one can ask if this correction can be expressed using some viscosity distribution ε(x) in the domain.This will be false in general.Even worse, the steepest gradient descend of the entropy can't be bounded using the first derivative of the respective component of the vector valued function u(x,t), incurring an infinitely large viscosity.All of the above reasons motivate us to devise alternative directions for the entropy correction.These alternative directions should have the following properties • The dissipation direction should have a filtering effect, i.e. when the direction only is used high order modes should be dissipated.• The direction should dissipate entropy.
• The dissipation should stem from a viscosity added to the hyperbolic flux.
Our new correction directions will be based on the construction of filters, i.e. operators that can regularize a solution u.A filter will in our case be a special Hilbert-Schmidt operator K [25].
Definition 1 (Filter).An operator K : L 2 (Ω) → L 2 (Ω) is said to be a filter if it is an integral operator whose pointwise evaluation results in a weighted average, i.e.
[Ku](x) = Ω k(x, y)u(y) dy, with ∀x ∈ Ω : is satisfied and the kernel k is of bounded Hilbert-Schmidt norm.
We are especially interested in conservative filters as they do not destroy the conservation of our basic schemes when they are applied on a per cell basis.

Lemma 4 (Conservative filter). A filter K
if it can be written as an integral operator with a kernel with mass one, i.e.
[Ku](x) = Ω u(y)k(x, y) dy, with ∀y ∈ Ω : Please note that the weighted average property is stated using the integration w.r.t. the second variable while the conservation results from the unit measure in the first variable.Obviously a convolution with a convolution kernel satisfying R k(y) dy = 1 satisfies both as k(x, y) = k(x − y) holds in this case, but not every operator satisfying these properties is a convolution.Especially when one is interested in bounded domains convolutions are not an option, but there still exist suitable smoothing operators.
Theorem 1 (Universally dissipative filters).A conservative filter K is dissipative for all convex entropies U, if it can be written as a conservative filter with a positive kernel, i.e.
Proof.Using Jensens inequality [30], the positivity and conservation of the filter allows us to show These theorem shows that the first and second bullet above can be satisfied by an integral operator with a suitable kernel.An example of a dissipation that can be identified with a positive conservative filter is the filtering by the time evolution of on the entire domain as the assorted filter has the heat kernel as kernel function [10], Further, this filtering obviously stems from viscosity and has therefore a direct physical interpretation.It is known that while a positive integral operator always dissipates entropy a high order finite-difference implementation will not dissipate all entropies [28] and similar theorems hold for higher even derivatives even in the analytic case.We will therefore outline how to construct a filter that is dissipative in the semidiscrete and fully discrete setting and can therefore be used as a descent direction.We begin by stating some discrete equivalents of the theorems above and will analyze if usual dissipations/filters satisfy this property.We will assume that ω k ≥ 0 is a positive quadrature rule on the cell T for the rest of the chapter and all notions of conservation for our filters will be centered around being conservative with respect to this quadrature rule.For a general DG method with dense mass matrix a quadrature can be calculated via ∑ l M lk = ω k , i.e. by entering the constant one into the discretised inner product, but positivity is not guaranteed in general.A general view of our plan could be to not discretise the second derivative, but its action as the generator of a Hilbert-Schmidt operator.We will therefore, when given a discrete filter, consider also its (discrete) generator.
Definition 2 (Conservative and positive filter generator).Let G ∈ R (p+1)×(p+1) be a square matrix.We call this matrix a filter generator if ∀k ∈ {1, . . ., p + 1} : Obviously, the definition of the conservative positive discrete filter mirrors the definition of such a filter in the continuous case using the quadrature rule.The definition of the averaging property on the other hand is not based on the quadrature rule, as this rule is not used when applying the filter pointwise Forward Euler steps connect the generators defined above with the filters, as we will see in the lemma below.
As the identity is conservative follows The positivity follows as for non-diagonal elements, k ̸ = l =⇒ (I + ∆tG) kl = ∆tG kl ≥ 0 is satisfied for any positive timestep size while the given restriction is needed to enforce It is clear that a discrete filter that is positive and conservative is also dissipative by reiterating through the arguments given above for the continuous case.Sadly, it is also true that while in the continuous case the filter which is generated by the second derivative, i.e. the heat kernel, is positive, the second derivative discretised in our DG method is not a positive generator and also does not generate a positive filter directly.We will therefore show how to design a generator generating an approximation of the heat kernel for forward Euler steps, thereby even allowing to prove the dissipativity of the entropy dissipation operator for finite time steps.The basis will be the heat equation with varying heat conductivity α(x) [19] on the (reference) element in conjunction with Neumann boundary conditions.The Neumann boundary conditions enforce the conservation of the resulting solution operator as any change of the cell mean values must happen through the numerical flux of the basic DG method.Discretising this problem [19] with the nodal basis of the basic DG method that is a continuous Galerkin method in this case because a single element is considered, yields As noted before, in general there exists no ∆t > 0 where I +∆t(−M −1 Q) is a positive operator because the negative elements in (−M −1 Q) prohibit it from being a positive generator.Yet the following theorem shows that the exact ODE solution to this problem for a t > 0 big enough is in fact eligible as a filter.
Theorem 2. If the quadrature ω is exact on V T , the solution of (11) for a positive initial condition u 0 ∈ R p+1 satisfies for all t > 0 Further, for a t > 0 big enough it follows ∀k ∈ {1, . . ., p + 1} : u k (t) ≥ 0.
As the quadrature is exact for the basis functions the same follows for the discretisation, and this shows the conservation.The matrix C used to describe the solution has the explicit form [25, sec. 34] u(0).
Multiplying this matrix with the vector v ∈ V representing the function 1 from the right reveals This already shows the second result as the nodal representation of C(t) must have unit row sum.The matrix −Q is negative semi-definite, the v vector is in its null space.If another linearly independent u ∈ V would be in its null space it would follow and this is a contradiction to ∂ u ∂ x ̸ = 0, as u was assumed non-constant.Therefore, there exists an orthonormal eigenvalue decomposition of the discretisation whose eigenvalues, apart from the constant eigenfunction ψ 1 = v with eigenvalue λ 1 = 0, are bounded away from zero, ∀k ∈ {1, . . ., p + 1} : We assume that the eigenvectors are sorted by increasing absolute value of the corresponding eigenvalues, The solution therefore converges to the average of u(0), as holds.Because a positive initial condition has a positive average the solution will converge to this positive average.
Using the theorem above we can construct filters ϒ simply by calculating the matrix C(t) = G used in the proof above.This matrix which maps an initial state onto the solution at time t is always a conservative filter, and when t is large enough also positive.In the implementation the suitable t was found using a bisection algorithm.Using ϒ = (C(t) − I)/t the corresponding generator can be found.We note in passing that numerous other possibilites exist to define a positive conservative filter as defined above, but that the method given above defines a filter than can be associated with viscosity.Lemma 6. Assume the null space of G consists only of constants.Then for a non-constant u and a strictly convex entropy U it holds If U is just convex only ≫≤≪ applies in the equation above.
Proof.The discrete dissipativity follows from the positive conservative filter property of G for λ ∆t > 0 small enough as in lemma 5 in conjunction with the strict convexity and Jensens inequality in the strict sense.Let now ∆t be fixed and small enough for all λ ∈ [0, 1], and denote by ε = E T (u + ∆tGu) − E T u < 0 the entropy dissipation for λ = 1.The convexity of U implies Entering this into the definition of the derivative of E T with respect to λ shows If U is not strictly convex the case ε = 0 is possible, reducing the result to ≫≤≪.
The last step consists of selecting a suitable viscosity distribution α, i.e. one that is zero at the endpoints.The standard mollifier is smooth and zero at the ends of the reference element.Further, even its derivatives vanish there.It was therefore selected.

Stable computation of the correction size required and timestep restrictions
After we have calculated the entropy dissipation needed and a suitable direction υ = Gu T one would guess we only have to calculate λ as in ( 6) via .
It turns out that this process is significantly more intricate than one would expect as this computation has to be stable with respect to roundoff errors.Further, our estimates on the entropy dissipation can only estimate the entropy dissipation that can take place at the interface between two adjacent cells, but are not able to give an estimate of how this dissipation is split between the two cells.Our method of calculating suitable values of λ T therefore consists of two steps.First, is calculated to enforce the per cell entropy dissipativity In a second step a correction to enforce an entropy rate high enough is determined for all θ ∈ Θ.Both corrections are then added together

ER
for all cells T ∈ T .Round-off errors tend to influence the calculation out of two reasons.The division by dU du , υ in equation ( 12) and ( 13) can approach a division by zero for a solution approaching a constant in the cell, as υ → 0 follows in this case.Further, we saw in lemma (1) that the entropy inequality predictor can vanish with a high order for smooth solutions, and an accurate DG scheme will also have a vanishing entropy error vanishing with a high order.The difference of these two values, i.e. the denominator of the fraction above, will in general not vanish that fast because round-off in the difference becomes important.Therefore λ will, for highly resolved smooth solutions, be to big because roundoff errors propagate into the calculation.Our solution to this problem is to calculate every time a λ is calculated by a division in the procedure above.Here, a shall be the nominator, b shall be the denominator and c shall be a suitable bound on the round-off error, a constant small with respect to a, b but large with respect to the machine precision.
In our implementation this is selected as c = √ 10 −16 , i.e. the square root of the machine precision for a solution scaled to be of unit magnitude.The addition of c can be seen as the one-dimensional version of Tikhonov regularization [23].Clipping the calculation of λ at 0 ensures that if a or b become negative from rounding errors λ will not become negative, i.e. λ υ will not be antidissipative.In a last step, the upper limit λ max is introduced for stability reasons as we want to enforce stability of If a Runge-Kutta time integration method can be written as convex combination of forward Euler steps, i.e. is Strong Stability Preserving (SSP) [16,31,32] and the time-steps satisfy ∆tλ ≤ 1 during every Euler step, the lemma 5 allows us to show that the solution is also entropy dissipative in the discrete case.If the time integration method used is just a conditionally stable Runge-Kutta method [9,41] we are interested in limiting the operator norm of ∆t ∥λ G∥ ≤ R in order to at least avoid a linear instability.The exact size depends on the time integration methods' stability region as we would like to fit the half-circle into the stability region of the method.
The tests below will focus on the cases p = 3 and p = 7 as the latter are popular in applications because they amount to 4 and 8 nodes, suitable for SIMD processor instructions.Results for values in between are essentially interpolatory to the ones reported for p = 3 and p = 7 and the source code is available to carry out tests for all values p > 0. Time integration will be carried out using the SSPRK(4, 3) method for most solutions, while the convergence analysis for p = 7 below will use the Hairer-Wanner DOPRI8 method, to achieve the needed convergence speed of the time integration.In all images below the ansatz functions of all cells are shown without any post-processing.

Shock tube tests
First, a series of shock tube tests was done to highlight the effectivity of the entropy correction in shock calculations as this is the primary aim of this publication.
The shock tube tests were always carried out for two different numbers of cells.First for N = N typ /(p + 1) cells, where N typ = 100 is the usual number of cells used in comparisons for Finite-Volume methods.This was done so that the same number of degrees of freedom has to be saved.The results look satisfactory and highlight the effectivity of the method in figures 5, 6, 7, 8.All shocks are sharp and concentrated to less than one cell width.Yet, only slight overshoots and oscillations are visible directly around the shocks.These distortions are confined to the cell directly next to the shock.Contact discontinuities are slightly smeared over one cell, but after they have been smeared to this width no additional smearing takes place.The computational complexity per timestep is still low as no recovery stencil selection has to be carried out and only 1/(p + 1) times the number of two-point fluxes need to be evaluated.Because some other publications use 100 cells also for DG methods we carried out the tests once more for N = 100 cells, amounting to 400 and 800 degrees of freedom for orders p = 3 and p = 7.

Numerical validation of the entropy rate criterion
To verify the entropy rate criterion the total entropy of the solution to the first shock tube above was compared to the solution calculated by a Lax-Friedrichs scheme with 3•10 4 cells.Similar comparisons were carried out in [20,21,22].Please note that the Godunov solver used previously was swapped for a LF scheme to evade the need for an exact Riemann solver.This is also supported by our finding in lemma 1 and corollary 1 as a Lax-Friedrichs solution therefore has to comply with the entropy rate criterion.A scheme should in these     comparisons have the same entropy dissipation rate (in the limit) as the Lax-Friedrichs scheme in the limit.Comparisons for orders p = 3 and p = 7 in figure 9 show that this seems to be the case.The DG scheme always has an entropy that lies below the entropy of the LF scheme.As the entropy inequality for vanishing viscosity solutions is also desirable it was also verified on a per-cell basis.We just note that the small positive violations in figure 9 are of the same magnitude as the precision achievable during the calculation of λ using our procedure with double precision floats.

Shu-Osher test
To showcase a combination of shocks and smooth areas the well established shock-sine interaction problem from [32, Problem 8] was tested.The initial conditions are given by The parameter ε was set to the canonical value of ε = 0.2.The results look satisfactory already when only N = 100 cells are used in the calculation.Yet, we note that this already corresponds to 400 and 800 degrees of freedom for the selected orders.When N = 200 cells are used the solution is nearly indistinguishable from the reference solution.

Convergence Analysis
While the main aim of our modification was to devise a new DG scheme usable for shockcapturing calculations the scheme also converges with high order of accuracy for smooth solutions in our experiments.As an example the solution of ρ 0 (x, 0) = 3.857153 + ε(x) sin(2x), v 0 (x, 0) = 2.0, p 0 (x, 0) = 10.33333, with ε(x) = e (x−3) 2 .
and periodic boundary conditions was calculated using our modified DG method.The analytical solution for this test problem is ρ(x,t) = 3.857153 + ε(x − 2t) sin(2x − 4t), v(x,t) = 2.0, p(x,t) = 10.33333, with suitable periodic boundary conditions.After the solution was calculated for N = {10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100} cells for p = 3 and with the same stepping up to 50 cells for p = 7 up to t = 5 the L 1 and L 2 errors were calculated.The convergence in figure 11 seems to take place with too high an order for the ansatz polynomials used.The reason for this could be that the accuracy of the basic scheme is significantly higher for these solutions than the accuracy of the corrected scheme, because the entropy dissipation estimate still falsely reports high amounts of entropy dissipation.When the grid is refined the entropy dissipation estimate converges with a higher speed than the basic scheme following lemma 3 and because the error introduced to enforce the dissipation dominates a higher convergence speed than expected is observed.

Timestep Analysis
An important result of any modification to a basic scheme can be an impact on the allowed timestep size.In the first part of this publication [21] this influence was tested by measuring the maximal timestep possible before a blow-up occurs.This was done once more.The maximal timestep possible for the first shock tube for orders p = 3 and p = 7 is shown in figure 12. Obviously this timestep is acceptable and when corrected for the larger maximal wave speed of the Riemann problem used for testing, larger than the timestep reported in the previous part, highlighting the superiority of the new dissipation direction.

Conclusion
The method described in [21] to enforce an entropy rate criterion for DG methods was improved.By using a direct indicator for the entropy dissipation the error indicator used before could be replaced, resulting in a lower dissipation in situations like contact discontinuities.For smooth solutions this new method to quantify the amount of dissipation needed converges significantly faster to zero than the error estimate used before, and therefore allows us to recover the convergence speed of the basic DG scheme that was reduced by one degree before.Further, the direct quantification of the entropy dissipation needed allowed us to consider different dissipation directions, especially combining smoothing and dissipation and therefore bridging into the field of modal filtering.The effectivity of the refined method was demonstrated for the Euler system of gas dynamics.The method is not only high order accurate but also able to handle shocks, contact discontinuities, and rarefactionwaves.The next logical steps can be the application to two-dimensional problems, the application of the designed entropy inequality predictors to other schemes like continuous Galerkin and Spectral Volume schemes, where several adjustments will be needed, and revisiting the splitting into a fully discrete scheme already explored in [21].The presented method to estimate the entropy dissipation needed could also be used with artificial viscosity shock-capturing as for example described in [13].

Competing Interests
The author has no relevant financial or non-financial interests to disclose.

Data Availability
The commented implementation of the schemes is available under https://github.com/simonius/dgdafermos.
8 Bibliography Entropy variables in cellT ∂U ∂ u (u T (x,t))Vector of nodal values of the entropy variables in cell T at t∂ U ∂ u T (t)Interpolation of the entropy variables in cellT on V ∂ U ∂ u T (x,t)The canonical inner product between a, b ∈ R n a • b or (a)•(b)

Figure 2 :
Figure 2: A set of noninteracting HLL approximate Riemann solutions.

Application of the entropy inequality predictor to an open overlap θ k+ 1 2 of
the domain -centered on the discontinuities.

Figure 3 :
Figure 3: Construction and application of the generalized HLL dissipation estimate.

Figure 4 :
Figure 4: Use of alternative dissipation directions.The direction − ∂ f ∂ x shall denote the L 2 projection of the exact solutions' derivative ∂ u∂t onto V .Our approximation du T dt is not entropy dissipative in this example and should be corrected into the entropy dissipative half-space characterized by the normal − dU du .The direction ∇ 2 u that stems from a discretization of the heat kernel is suitable for this correction and has additional benefits compared to − dU du .While the diffusion also has a smoothing effect the addition of − dU du can even result in a sharpening effect.Higher even derivatives like ∇ 8 u will smooth the solution but will not result in a dissipation for all entropies.

Lemma 5 (
Connecting generators and filters).It holds G conservative as generator =⇒ ϒ = I +∆tG conservative as filter.Let further ∆t max l |G ll | ≤ 1.Then it follows G positive as generator =⇒ ϒ positive as filter.Proof.We begin by showing the conservativity and filter property.It holds p+1 ∑ l=1

Figure 5 :
Figure 5: Shock tube 1 at t = 1.8 with 25 cells corresponding to 100 degrees of freedom and 100 cells corresponding to 400 degrees of freedom (p=3).

Figure 6 :
Figure 6: Shock tube 1 at t = 1.8 with 13 cells corresponding to 104 degrees of freedom and 100 cells corresponding to 800 degrees of freedom (p=7).

Figure 7 :
Figure 7: Shock tube 2 at t = 1.2 with 25 cells corresponding to 100 degrees of freedom and 100 cells corresponding to 400 degrees of freedom (p=3).

Figure 8 :
Figure 8: Shock tube 2 at t = 1.2 with 25 cells corresponding to 100 degrees of freedom and 100 cells corresponding to 400 degrees of freedom (p=7).

Figure 9 :
Figure 9: Entropy tests for the first shock tube.

Table 2 :
Used schemes in the numerical tests condition [32, 39, Problem I, Section 4.3.3 and Problem 6a]) is