On the order reduction of entropy stable DGSEM for the compressible Euler equations

Is the experimental order of convergence lower when using the entropy stable DGSEM-LGL variant? Recently, a debate on the question of the convergence behavior of the entropy stable nodal collocation discontinuous Galerkin spectral element method (DGSEM) with Legendre-Gauss-Lobatto nodes has emerged. Whereas it is well documented that the entropy conservative variant with no additional interface dissipation shows an odd-even behavior when testing its experimental convergence order, the results in the literature are less clear regarding the entropy stable version of the DGSEM-LGL, where explicit Riemann solver type dissipation is added at the element interfaces. We contribute to the ongoing discussion and present numerical experiments for the compressible Euler equations, where we investigate the effect of the choice of the numerical surface flux function. In our experiments, it turns out that the choice of the numerical surface flux has an impact on the convergence order. Penalty type numerical fluxes with high dissipation in all waves, such as the LLF and the HLL flux, appear to affect the convergence order negatively for odd polynomial degrees $N$, in contrast to the entropy conserving variant, where even polynomial degrees $N$ are negatively affected. This behavior is more pronounced in low Mach number settings. In contrast, for numerical surface fluxes with less dissipative behavior in the contact wave such as e.g. Roe's flux, the HLLC flux and the entropy conservative flux augmented with 5-wave matrix dissipation, optimal convergence rate of $N+1$ independent of the Mach number is observed.

In addition to the very robust dissipative entropy stable versions, it is also possible to construct virtually dissipation free variants by choosing appropriate element interface numerical fluxes.These entropy conserving variants all show an odd-even behavior when experimentally testing the order of convergence, e.g.[8,18], where the observed convergence order for even polynomial degrees N is N and for odd N is N + 1. Lately, a discussion emerged in the community, with interesting debates during the recent ICOSAHOM conference in London, where researchers reported non-optimal convergence behavior of the entropy stable DGSEM-LGL even with dissipative numerical surface fluxes, e.g.[5].
This paper contributes to this discussion and presents results of an experimental convergence order study for the compressible Euler equations with (i) the standard DGSEM with either Gauss and LGL nodes, (ii) the entropy stable DGSEM with LGL nodes.For these nodal schemes, we test the convergence order with different numerical surface fluxes and report the results depending on the Mach number of the test case.The remainder of the paper is organized as follows: in the next section we describe the numerical model for our numerical experiments, in Sec. 3 we present our observed experimental convergence orders for different configurations and draw our conclusion in Sec. 4.

Numerical Model
We consider the compressible Euler equations defined in the domain The state vector contains the conservative variables and the advective flux components are p, E are the mass density, fluid velocities, pressure and total energy.We close the system with the ideal gas assumption, which relates the total energy and pressure where γ denotes the adiabatic coefficient.
For our discretization, we subdivide the domain into non-overlapping hexahedral elements.For each element, we define a transfinite mapping to a unit reference space and use this mapping to transform the equations (1) from physical to reference space.A weak form is created by taking the inner product of the transformed equation with a test function.We use integration-by-parts for the flux term and approximate the resulting weak form as follows: the conservative variables are approximated by a polynomial in reference space with degree N, interpolated at the Gauss or LGL nodes.The volume fluxes are replaced by a standard interpolation of the non-linear flux function at the same Gauss/LGL nodes (standard DGSEM-Gauss or DGSEM-LGL), see e.g.[13].For the LGL variant, we are also able to introduce the split form volume integral based on entropy conserving and kinetic energy preserving numerical volume fluxes (Split-DGSEM), e.g.[12] and [17], resulting in either the entropy conserving or entropy stable DGSEM variants, depending on the choice of numerical surface flux.

Convergence results
In this section, we compare the convergence of the standard DGSEM and the entropy conservative and entropy stable discretization for different choices of the numerical flux and polynomial degrees N = 2, 3, 4, 5.
We choose the test case of a two-dimensional density wave, with a constant pressure and transported with a constant velocity, which was proposed for onedimensional convergence tests in [4].The density evolves as with a prescribed velocity (v 1 , v 2 ).The pressure is chosen as p = 1/γ with γ = 1.4,so that the sound speed ranges between c = 0.95 . . .
and Ma ≈ 3.5 with (v 1 , v 2 ) = (2.5, 2.4).The experimental order of convergence (EOC) is computed with the L 2 error of the density at t = 1.
The convergence study is performed with the open source, three-dimensional curvilinear split-form DG framework FLUXO (www.github.com/project-fluxo).As the test case is two-dimensional, we use fully periodic cartesian meshes of the domain [−1, 1] 3 with an equal number of elements in x-and y-directions and always 1 element in z-direction.Note that h 0 in the convergence tables refers to the coarsest mesh level, which is 4 2 elements for N = 2, 3 (h 0 = 1/2) and 2 2 elements for N = 4, 5 (h 0 = 1).
All simulation results are obtained with an explicit five stage, fourth order accurate low storage Runge-Kutta scheme [3], where a stable time step is computed according to the adjustable coefficient CFL∈ (0, 1], the local maximum wave speed, and the relative grid size, e.g.[9].We made sure that the time integrator did not influence the spatial convergence order, by adjusting the CFL number accordingly.

Volume disc.
Mach mesh Experimental order of convergence of L 2 error to the exact density (4), using the standard DGSEM-Gauss with HLL and Roe fluxes.Full order is marked with ( N + 1) and an order reduction with .
Experimental order of convergence of L 2 error to the exact density (4), using DGSEM-GL with HLL and Roe fluxes.Full order is marked with ( N + 1) and an order reduction with .

Standard DGSEM
The convergence of the standard DGSEM with Gauss-Legendre nodes (DGSEM-Gauss) and with Legendre-Gauss-Lobatto (DGSEM-LGL) is shown in Table 1 and Table 2, for the three Mach numbers and two choices of the numerical flux, namely the HLL and the Roe flux.The results of the local Lax-Friedrichs (LLF) and the HLLC fluxes are reported in the Appendix, as HLL is similar to LLF, and HLLC is the same as Roe, see Table 4 and Table 5.
For the HLL flux and the low Mach number Ma = 0.2, we observe an odd-even behavior with an order reduction for even polynomial degrees N = 2, 4. Also for Ma = 1.0, the convergence for even degrees is slightly affected, whereas for the high Mach number, all fluxes converge with full order.Comparing the L 2 errors of the finest mesh for HLL and Roe for the low Mach number, HLL is less accurate for N = 2, 4 and more accurate for N = 3, 5.
All numerical fluxes are approximate Riemann solvers, but the LLF and HLL only use the maximum wave speeds, whereas the HLLC and Roe also take the contact wave into account, and therefore keep the full order of the scheme for all Mach numbers for this test case.

Split-DGSEM +
Table 3 Experimental order of convergence of L 2 error to the exact density (4), using entropy conservative ECKEP flux and entropy stable HLL and ECKEP-Roe fluxes.Full order is marked with ( N + 1) and an order reduction with .

Entropy conservative and entropy stable DGSEM
Now, we investigate the order reduction of the entropy conservative and entropy stable discretizations.Here, the standard DGSEM volume integral is replaced by split-form formulation (Split-DGSEM) using a two-point entropy conservative and kinetic energy preserving flux ECKEP.If we choose the ECKEP flux at the surface, we get an entropy-conserving scheme.For entropy stability, we can use the LLF or HLL flux directly at the surface, or use the ECKEP flux and add a dissipation term, which must still satisfy the entropy inequality condition.In Winters et al. [17], such dissipation terms are carefully derived, using either only the maximum wave speed (LLF-type) or incorporating all waves (Roe-type), which we will refer to as ECKEP-LLF and ECKEP-Roe fluxes.
In Table 3, we summarize the convergence of the dissipation-free ECKEP flux, the HLL and ECKEP-Roe flux.The results for LLF and ECKEP-LF fluxes are found in the Appendix in Table 6, as they have the same convergence and error levels as the HLL flux.As expected, the dissipation-free surface flux (ECKEP) produces an order reduction for all Mach numbers for N = 3, 5, and for N = 2 full order not kept in the last refinement step.
If we simply use the HLL flux, we have an entropy stable scheme, but an order reduction for N = 2, 4 can be observed for the low Mach number flow, analogously to the standard DGSEM-LGL scheme.Interestingly, the odd-even behavior switches between entropy conserving and entropy stable fluxes.
The ECKEP-Roe entropy stable flux accounts for all waves of the Riemann problem and adjusts the dissipation for each wave accordingly, which gives full order convergence for all Mach numbers.

Conclusions
In this work, we report the convergence of standard DGSEM Gauss and Gauss-Lobatto schemes to entropy conservative (EC) and entropy stable (ES) DGSEM schemes for the Euler equations, as there have been findings of order reduction for EC and ES schemes.We choose a simple density transport test case on a periodic domain and investigate the influence of the Mach number of the transport velocity.
The EC scheme is dissipation free and an order reduction is observed by the convergence study presented here, confirming many similar observations found in literature.We also confirm that the ES scheme can have an order reduction for low Mach numbers, but only if the entropy stable numerical flux relies on simple approximate Riemann solvers such as local Lax-Friedrichs or HLL.If all waves are accounted for in the dissipation term of the entropy stable flux as presented in [17], the full order is observed for all Mach numbers.IN addition, we reproduce the same behavior for the standard DGSEM Gauss and Gauss-Lobatto schemes, where the LLF and HLL fluxes suffer from order reduction at low Mach number, and HLLC and Roe fluxes have full order for all Mach numbers.
We want to emphasize that the present convergence study should be seen merely as an observation, confirming that the numerical flux can have strong influence on the convergence order for both the standard DGSEM and the entropy stable DGSEM.Also, we stress that in our tests the order reduction is related to the form of the dissipation term in the numerical surface flux and is not related to the insufficient integration precision of the LGL-quadrature.
In the authors' experience, a convergence study using a manufactured solution technique can be misleading, as full convergence order is found independent of the choice of numerical flux.Hence, the introduction of a source term to balance the prescribed solution overcomes possible deficiencies of the surface fluxes, showing the limit of the manufactured solution technique in this context.In the Appendix, the convergence results of a manufactured solution are reported.

Volume disc.
Mach mesh Experimental order of convergence of L 2 error to the exact density (4), using DGSEM-Gauss with LLF and HLLC fluxes.Full order is marked with ( N + 1) and an order reduction with .
Note that the average Mach number in the domain is Ma = 0.8.Inserting (5) into the Euler equations, and using the fact that spatial and time derivatives are g = ∂ x 1 g = ∂ x 2 g = −∂ t g, we get an additional residual To solve the inhomogeneous problem, we subtract the residual from the approximate solution in each Runge-Kutta step.Moreover, we run the test case up to the final time t = 1.0.
In the convergence results for the standard DGSEM Gauss and Gauss-Lobatto, we see that the LLF flux still leads to an order reduction for N = 2, 4, whereas full order is found for the HLL, HLLC and Roe fluxes, see Table 7 and Table 8.

DGSEM-LGL +
Table 5 Experimental order of convergence of L 2 error to the exact density (4), using DGSEM-GL with LLF and HLLC fluxes.Full order is marked with ( N + 1) and an order reduction with .
In Table 9 the entropy conservative scheme shows again an order reduction for N = 3, 5, and the LLF-Type dissipation too, for N = 2, 4, and for this test case, all entropy stable schemes exhibit full order.
1.05.Thus, by changing the velocity, we change the Mach number of the flow Ma = | → v|/c.Three Mach numbers are chosen: Ma

Table 6
Experimental order of convergence of L 2 error to the exact density (4), using entropy stable LLF and ECKEP-LLF flux.Full order is marked with ( N + 1) and an order reduction with .

Table 7
Experimental order of convergence of L 2 error of density for the manufactured solution (5), using DGSEM-Gauss with LLF, HLL, HLLC and Roe fluxes.Full order is marked with ( N + 1) and an order reduction with .EOC L 2 ( ) EOC L 2 ( ) EOC L 2 ( ) EOC