Entropy Stable Discontinuous Galerkin Schemes on Moving Meshes for Hyperbolic Conservation Laws

This work is focused on the entropy analysis of a semi-discrete nodal discontinuous Galerkin spectral element method (DGSEM) on moving meshes for hyperbolic conservation laws. The DGSEM is constructed with a local tensor-product Lagrange-polynomial basis computed from Legendre–Gauss–Lobatto points. Furthermore, the collocation of interpolation and quadrature nodes is used in the spatial discretization. This approach leads to discrete derivative approximations in space that are summation-by-parts (SBP) operators. On a static mesh, the SBP property and suitable two-point flux functions, which satisfy the entropy condition from Tadmor, allow to mimic results from the continuous entropy analysis, if it is ensured that properties such as positivity preservation (of the water height, density or pressure) are satisfied on the discrete level. In this paper, Tadmor’s condition is extended to the moving mesh framework. We show that the volume terms in the semi-discrete moving mesh DGSEM do not contribute to the discrete entropy evolution when a two-point flux function that satisfies the moving mesh entropy condition is applied in the split form DG framework. The discrete entropy behavior then depends solely on the interface contributions and on the domain boundary contribution. The interface contributions are directly controlled by proper choice of the numerical element interface fluxes. If an entropy conserving two-point flux is chosen, the interface contributions vanish. To increase the robustness of the discretization we use so-called entropy stable two-point fluxes at the interfaces that are guaranteed entropy dissipative and thus give a bound on the interface contributions in the discrete entropy balance. The remaining boundary condition contributions depend on the type of the considered boundary condition. E.g. for periodic boundary conditions that are of entropy conserving type, our methodology with the entropy conserving interface fluxes is fully entropy conservative and with the entropy stable interface fluxes is guaranteed entropy stable. The presented proof does not require any exactness of quadrature in the spatial integrals of the variational forms. As it is the case for static meshes, these results rely on the assumption that additional properties like positivity preservation are satisfied on the discrete level. Besides the entropy stability, the time discretization of the moving mesh DGSEM will be investigated and it will be proven that the moving mesh DGSEM satisfies the free stream preservation property for an arbitrary s-stage Runge–Kutta method, when periodic boundary conditions are used. The theoretical properties of the moving mesh DGSEM will be validated by numerical experiments for the compressible Euler equations with periodic boundary conditions.

with summation-by-parts (SBP) operators and proved that two-point entropy conservative fluxes can be used to construct high-order schemes when the derivative approximations in space are SBP operators. A SBP operator provides a discrete analogue of the integrationby-parts formula [19,22,38]. It is worth to mention that the derivative matrix in the DGSEM provides a SBP operator, if the tensor-product Lagrange-polynomial basis is computed from Legendre-Gauss-Lobatto (LGL) points and interpolation and quadrature are collocated. Gassner et al. [23,24] showed that split forms of the partial differential equations can be discretely recovered when specific choices of numerical volume fluxes in the flux form volume integral of Fisher and Carpenter are chosen. Thus, an entropy stable DGSEM can be constructed by the following building blocks: (1) The derivative matrix satisfies the SBP property.
(2) There are two-point flux functions with Tadmor's discrete entropy condition that can be extended to high-order in a split form DG framework.
This methodology has been used in the construction of high-order entropy stable DGSEM on quadrilateral/hexahedral elements, e.g. [4,23,61], or on triangular/tetrahedral elements, e.g. [6,10,13]. All these methods are provably entropy stable and the semi-discrete entropy analysis for them is based merely on the properties of the SBP operators and the assumptions that the time integration is exact. Additionally, properties like positivity preservation (of the water height, density or pressure) must be satisfied on the discrete level. The exactness of quadrature in the spatial integrals of the variational form is not necessary. Available entropy stable moving mesh methods are for instance the low order continuous finite element method by Guermond et al. [27] and a spectral collocation based approach published during the extended review process of the current paper by Yamaleev et al. [66]. The remainder of the paper is organized as follows: The ALE transformation and continuous entropy analysis is presented in the Sects. 2.1, 2.2 and 2.3. The framework for the spectral element discretization with the SBP operator is given in the Sect. 2.4 and the DG split form framework is presented in the Sect. 2.5. The moving mesh DGSEM is finally presented in the Sect. 2.5. A discrete entropy analysis for the moving mesh DGSEM is given in the Sects. 2.6 and 2.7. Furthermore, in Sect. 2.8 it is proven that the moving mesh DGSEM satisfies the free stream preservation property. In Sect. 4, numerical examples with the compressible Euler equations are presented to validate our theoretical findings.

Entropy Stable DGSEM on Moving Meshes
The main goal of this work is the construction of an entropy stable moving mesh DGSEM. On static meshes, it is possible to construct high-order entropy stable DGSEM, if the derivative matrix is an SBP operator and entropy conservative two-point flux functions are available. This methodology has been used in the construction of high-order entropy stable DGSEM on quadrilateral/hexahedral elements, e.g. [4,23,61]. In this section, it will be shown that similar ideas can be used to construct high-order entropy stable moving mesh DGSEM. The construction of the entropy stable moving mesh DGSEM will be presented for an arbitrary symmetrizable and hyperbolic system of conservation laws The block vector nomenclature in [24] simplifies the analysis of the system (2.1) on curved elements. Thus, we translate the conservation law (2.1) in block vector notation. A block vector is highlighted by the double arrow The dot product of two block vectors is given by Thus, in particular, the spatial gradient of the conservative variables is defined by The gradient of a vector valued function g = [g 1 , g 2 , g 3 ] T is a second order tensor, written in matrix form as where ⊗ is the outer product of two vectors in a three dimensional space. The dot product (2.3) and the spatial gradient (2.6) are used to define the divergence of a block vector flux as Moreover, for a vector valued function g and the conservative variables, we have the product rule ∇ x · ( g u) = ∇ x · g u + g · ∇ x u (2.9) with respect to the dot products (2.3) and (2.4). These notations allow to write the conservation law (2.1) in the compact form (2.10)

Building Blocks of the ALE Transformation for Hexahedral Curved Meshes
In order to set up the moving mesh DGSEM in the Sect. 2.5, we make for all t ∈ [0, T ] the assumptions: (A1) For a fixed number K ∈ N the physical domain (t) can be subdivided into K time-dependent, non-overlapping and conforming hexahedral elements, e κ (t), κ = 1, . . . , K . These elements can have curved faces. (A2) The time-dependent elements e κ (t) are mapped into the spatial computational domain E = [−1, 1] 3 with a bijective isoparametric transfinite mapping Winters constructed in his PHD thesis [65] a mapping for this set up. Like in [65], it is assumed that the curved faces satisfy for all τ ∈ [0, T ] The location of the curved faces is sketched in Fig. 1. The curved faces of an element e κ (t) are approximated as interpolation polynomials up to degree N such that with the curved faces 1 ξ 1 , ξ 3 , τ , 2 ξ 1 , ξ 3 , τ , 3 ξ 1 , ξ 2 , τ , 4 ξ 2 , ξ 3 , τ , 5 ξ 1 , ξ 2 , τ , and The mapping x (t) = χ ξ, τ connects E and e κ (t) Mesh curving techniques are discussed by Hindenlang et al. [30] and methodologies to construct a moving mesh with the properties (A1)-(A3) are given in the literature e.g. the book of Huang and Russell [32,Chapter 6,Chapter 7]. In many situations the moving mesh methodology depends on the underlying problem, e.g. [45].
The mapping provides the grid velocity field It is desirable that the grid velocity is continuous, since the mesh should be conforming and watertight at each time level. The next statement provides conditions on the element boundaries to guarantee that the grid velocity becomes continuous.

Lemma 2.1
Let e 1 (t) and e 2 (t) be two neighboring elements which share one of the faces where l i , l = 1, 2, and i = 1, 2, 3, 4, 5, 6, are the faces of the element e l (t). Furthermore, suppose that the faces l i (·, ·, τ ) are continuously differentiable in the time interval [0, T ]. Then the grid velocity field is continuous in the points which belong to the face that the elements share.
In Appendix A the Lemma 2.1 is proven in two dimensions. The three dimensional proof can be done by the same argumentation.

Transformation of the Conservation Law onto a Reference Element
In the following, we show that the system (2.10) can be transformed from a time-dependent element e κ (t) on the reference element E. The mapping (A.1) provides the covariant basis vectors Hence, the contravariant block vector flux is given by Since the elements {e k (t)} K k=1 are time-dependent, the time evolution of the quantity J needs to be analyzed. Thus, we apply Jacobi's formula (cf. e.g. Bellman [3]) and obtain by (2.15), The chain rule formula and the identity (2.26) provide Next, we plug (2.30) into Eq. (2.31), apply the product rule (2.9) and rearrange. This provides the equation Finally, we combine the identities (2.27) and (2.32) to write the the conservation law (2.10) in the following form The formulation (2.33) is the representation of the system (2.10) on the time-independent reference element E for a time-dependent element e κ (t).

Remark 2.2
The metric identities (2.20) and the Eq. (2.30) provide the geometric conservation law (GCL) [18,28,41,42,46]. A numerical method to solve (2.10) on moving and deforming grids needs to satisfy both equations, otherwise the conservation properties of the conservation law (2.10) are not preserved. Farhat et al. [18,28,41] proved that the absence of these equations has a critical effect on the accuracy and stability of a moving mesh method. In particular, the preservation of constant states is no longer guaranteed, if the GCL is not satisfied on the discrete level.

Entropy Analysis in Three Dimensions
The system (2.1) is assumed to be symmetrizable. Thus, in particular, it is equipped with entropy/entropy flux pairs s, f s i , i = 1, 2, 3, (cf. Godunov [26] and Mock [49]). The strictly convex function s is the entropy function. The entropy function s provides the entropy variables w := ∂s ∂u , (2.35) and it follows by the chain rule The entropy flux functions and the flux functions in the conservation law are related and satisfy The identities (2.26) and (2.36) give Hence, we obtain with the identity (2.31) and the chain rule Therefore, the product rule provides the identity where we used the GCL (2.30) in the last step. Next, we apply the relation (2.37) for the entropy flux functions and obtain  Boundary conditions then need to be specified so that the bound on the entropy depends only on the boundary data. We will assume here that boundary data is given in a way that the right hand side in Eq. (2.44) is non-positive so that the entropy will not increase in time.
For discontinuous solutions the Eq. (2.44) is not satisfied, but under further assumptions it is possible to proof that a weak solution of (2.33) satisfies the inequality in the sense of distributions on E × (0, T ) (see Godlewski (2.46)

Building Blocks for the Spectral Element Approximation
A nodal approach is used for the spectral element approximation. The Lagrange basis functions are given by where the nodal points {ξ i } N i=0 are the LGL points. We note that ξ 0 = −1 and ξ N = 1. The Lagrange basis functions satisfy the cardinal property i ξ j = δ ji , (2.48) where δ ji is the Kronecker delta. On the reference element E = [−1, 1] 3 the solution and fluxes of the system (2.33) are approximated by tensor product Lagrange polynomials of degree N , e.g.
In the following, polynomial approximations are highlighted by capital letters, e.g. U is an approximation for the state vector u and F l , l = 1, 2, 3, are approximations for the fluxes f l , l = 1, 2, 3. The determinant J of the Jacobian matrix ∇ ξ χ is also approximated by tensor product Lagrange polynomials In particular, the interpolation operator for a function g is given by are sets of LGL points. Derivatives are approximated by exact differentiation of the polynomial interpolants. In general we have I N (g) = I N −1 g (cf. e.g. [7,36]), as differentiation and interpolation only commute if there are no discretization errors. However, the contravariant coordinate vectors need to be discretized in such a way that the metric identities (2.20) are satisfied on the discrete level, too. Kopriva [35] introduced the conservative curl form that computes From now on, the discrete contravariant coordinate vectors are denoted by Ja α β , when the curl form (2.52) has been used to compute these quantities.
Integrals are approximated by a tensor product extension of a 2N − 1 accurate LGL quadrature formula. Hence, interpolation and quadrature nodes are collocated. In one spatial dimension the LGL quadrature formula is given by where ω i , i = 0, . . . , N , are the quadrature weights and ξ i , i = 0, . . . , N , are the LGL quadrature points. The formula (2.54) motivates the definition of the discrete quantity for two functions f and g. We note that (2.55) satisfies wheren is the unit outward normal at the faces of the reference element E.
The spectral element approximation with LGL points for interpolation and quadrature provides a SBP operator Q = M D with the mass matrix M and the derivative matrix D. The mass matrix and the derivative matrix are given by (2.58) The important characteristic of this SBP operator is the property where B = diag (−1, 0, . . . , 0, 1). A SBP operator provides a discrete analogue of the integration-by-parts formula [19,22,38]. Finally, we note that in the LGL points

The Semi-discrete Discontinuous Galerkin Method
Now, we apply the notation introduced in Sect. 2.4 and construct a moving mesh DGSEM. We discretize the Eqs. (2.30) and (2.33) simultaneously. In this way, it is ensured that the Eq. (2.30) is satisfied on the discrete level [37,48,64]. First, we replace the solution u by (2.49), the Jacobian J by (2.50) and approximate the fluxes by the interpolation operator (2.51). Next, we multiply the GCL (2.30) by test functions ϕ ∈ P N (E), the Eq. (2.33) with ϕ ∈ P N (E, R p ), integrate the resulting equations and use integration-by-parts to separate boundary and volume contributions. The volume integrals in the variational form are approximated with the LGL quadrature. Then, we insert numerical surface fluxes ν * and ↔ G * at the spatial element interfaces. Afterwards, we use the SBP property (2.59) for the volume contribution to get the standard DGSEM in strong form: where we used the notation (2.55) and the notation (2.57) for the discrete surface integral.
The approximation of ν and the nonlinear flux ↔ g by the interpolation operator (2.51) causes aliasing errors in the standard strong form. The aliasing errors cannot be bounded and the errors are independent of the choice of the numerical surface flux. In Gassner [22] a detailed explanation and analysis of the aliasing problem is given. Furthermore, a specific reformulation of the volume integrals by using the skew-symmetry strategy has been developed to fix the aliasing problem. This approach has been enhanced and generalized by Gassner et al. in [23,24] with a technique developed for high-order FD schemes (LeFloch et al. [40] or Fisher and Carpenter [20]). The generalized approach is called split form DG framework. Here, we proceed similar as in [24] and replace the interpolation operators in the discrete volume integrals by derivative projection operators. The interpolation operator in the discrete equation for the GCL (2.30) is replaced by with the volume averages of the metric terms, e.g.
The derivative projection operator in the discrete equation for (2.33) is computed as in [24]. Thus, the operator is given by (2.64) We note that the discrete volume weighted contravariant vectors J a l , l = 1, 2, 3, in the derivative projection operator (2.62) and (2.64) are computed by the conservative curl form (2.52). The flux ↔ G EC is consistent and symmetric such that, e.g.
for i, j, k, m = 0, . . . , N . Furthermore, the flux functions G EC l , l = 1, 2, 3, satisfy for i, j, k, m = 0, . . . , N , the following discrete entropy conditions (2.67) The quantities and l are polynomial approximations which satisfy in the LGL points where W i jk , S i jk and F s l i jk are the nodal values of the polynomials Here, s represents an entropy for the system (2.1) with the corresponding entropy flux functions f s l , l = 1, 2, 3, and entropy variables w. Furthermore, the volume jumps in (2.67) are, e.g. (2.70) In Appendix B, flux functions with these properties are presented for the Euler equations.
Finally, for each element e κ (t) the semi-discrete moving mesh DGSEM can be written in the following form: The unit outward facing normal vector and surface element on the element side are constructed from the element metrics by To define the numerical surface fluxes in (2.71a) and (2.71b), we introduce notation for states at the LGL nodes along an interface between two spatial elements to be a primary "−" and complement the notation with a secondary "+" to denote the value at the LGL nodes on the opposite side. Then the orientated jump and the arithmetic mean at the interfaces are defined by When applied to vectors, the average and jump operators are evaluated separately for each vector component. Then the normal vector n is defined unique to point from the "−" to the "+" side. This notation allows to compute the contravariant surface numerical fluxes in We note that due to the assumptions made in Sect. 2.1, the mesh velocity is a continuous function and the averages reduce to the uniquely defined values on the surface. The contravariant surface numerical fluxes in (2.71b) are given bỹ where the Cartesian fluxes G EC l , l = 1, 2, 3, satisfy (2.65), (2.66), (2.67). We note that these fluxes are the baseline choices without interface dissipation, to get a baseline scheme that is entropy conservative.

Semi-discrete Entropy Conservation
The spatial integral of the entropy is bounded in time on the continuous level. Thus, it is desirable that a numerical method is stable in the sense that a discrete version of this integral is bounded in time, too. In the context of the moving mesh semi-discrete DGSEM (2.71), we are interested to find an upper bound for the quantitȳ where S = I N (s) is a polynomial approximation for the entropy s. Next, we prove the following statement for the semi-discrete moving mesh DGSEM.
Proof We proceed similar as in the continuous entropy analysis, use the polynomial approximation ϕ = I N (w) = W as test function in the Eq. (2.71b) and obtain (2.81) Next we multiply the Eq. (2.81) by ω i jk and compute the sum over all LGL nodes. This gives Since we assume time continuity for our semi-discrete analysis, we apply the product rule in time and obtain by (2.82) where we used in the last step the Eq. (2.71a) with the test function ϕ = . We note that the quantity is defined as a polynomial with the nodal values (2.68). In the Appendix C.1, the following equation is proven Here the polynomials F s l , l = 1, 2, 3, are given by (2.69). Moreover, we obtain by (2.73) and (2.74)   [14,31,53], is applied. Then the semi-discrete moving mesh DGSEM (2.71) satisfies the discrete entropy inequalityS (2.91)

Semi-discrete Entropy Stability
Entropy conservation can be merely expected when a reversible process is described by a system of PDEs. In general, conservation laws are describing irreversible processes with discontinuous solutions. Hence, it cannot be expected that the entropy conservative moving mesh DGSEM provides a physical meaningful discretization for the system (2.1). However, the entropy conservative flux at the element interfaces can be augmented by an artificial dissipation term. In the literature, there are different strategies to add dissipation to an entropy conservative flux. Here, dissipation is added via a matrix operator. This approach, for instance, has been used in the context of gas dynamics by Chandrashekar [8] or Winters et al. [63].
The conservative variables u can be written in dependence of the entropy variables w. Differentiation of the conservative variables u = u(w) provides the symmetric positive definite matrix ∂u ∂w , since the system (2.1) is assumed to be symmetrizable (cf. e.g. [29]). Thus, it follows by a Taylor expansion up to first order where the matrices R l , T l , depend on some averaged values of the states U − , U + and they are consistent with the right eigenvector matrix R l and the scaling matrix T l . The matrix | l | depends on the values λ l The matrix H l needs to be a symmetric positive definite matrix. Therefore, the matrix | l | has to be chosen carefully. In Appendix B.3, the matrices to construct the dissipation operator (2.96) for the compressible Euler equations are given. There it can be also seen which average values are used to evaluate the states U − , U + in the matrices.  77). We note that the dissipation operator (2.96) is not used to modify the entropy conservative fluxes in the derivative projection operator (2.64).
The numerical fluxesG EŜ n do not provide an entropy conservative scheme, but the result in Theorem 2.4 can be used to prove that the moving mesh DGSEM becomes entropy stable, such that the discrete mathematical entropy is bounded at any time by its initial data, when the numerical fluxesG EŜ n are used at the element interfaces and it is assumed that the time integration is exact and that properties like positivity preservation (of the water height, density or pressure) are satisfied on the discrete level.  The right hand side of the inequality vanishes, since the method is investigated with periodic boundary conditions. Hence, we obtain the the desired inequality (2.98) by integrating (2.102) over the temporal interval [0, T ].

Free Stream Preservation for the Moving Mesh DGSEM
In this section, we check the discretization of the geometric and metric terms in time. Since DG methods with the forward Euler discretization are unstable [9,16], we investigate directly the discretization by an explicit RK method with s ≥ 2 stages and the characteristic coefficients It is worth to mention that a Courant-Friedrichs-Lewy (CFL) restriction is necessary when an explicit s-stage RK method is used in the DG framework. In order to present the RK discretization of the semi-discrete DGSEM (2.71), it is beneficial to write the method in the equivalent nodal representation. This representation is for all i, j, k = 0, . . . , N , given by where the right hand sides are given by with the tensorial Lagrange polynomials i j k given by (2.47).
It should be noted that the solutions J i jk , i, j, k = 0, . . . , N , of the ordinary differential equations (ODEs) (2.103a) need to be positive. We note that the solutions J i jk , i, j, k = 0, . . . , N of the ODEs (2.103a) are not used to compute the volume weighted contravariant coordinate vectors J a l , l = 1, 2, 3, in the right hand sides (2.104) and (2.105). These vectors are computed from the mapping by the conservative curl form (2.52). Hence, the right hand sides (2.104) are are independent of J i jk , i, j, k = 0, . . . , N . Therefore, the solutions of the ODEs (2.103a) are positive, if the grid velocity does not cause to much distortion in the mesh which is ensured when the assumptions (A1)-(A3) are satisfied.
Next, the interval [0, T ] is divided in time points t n . The step size of the time discretization is t. The DGSEM solutions, the fluxes and the grid velocity field are approximated in the time points t n , e.g. U (t n ) ≈ U n . Then, the RK discretization of the semi-discrete DGSEM is given by for = 1, . . . , s : are sets of LGL points. Next, we prove that the fully-discrete split form RK-DGSEM (2.106) satisfies the free stream preservation property.

Theorem 2.8 Suppose the fully-discrete split form RK-DGSEM (2.106) is investigated with periodic boundary conditions and the solution of the scheme is given by
(2.107) Furthermore, since the metric terms are computed by the conservative curl form (2.52), we obtain (2.111) Hence, the solution of the RK stage (2.106b) is given by Since the parameter was arbitrary chosen, it follows U ( ) i jk for all = 1, . . . , s. In particular, it follows This completes the proof of Theorem 2.8.

Numerical Results
The numerical computations in this section are performed with the open source code FLEXI 1 and the three-dimensional high-order meshes for the simulations are generated with the open source tool HOPR. 2 We present tests on three dimensional moving hexahedral curved meshes for the compressible Euler equations. Based on these tests we will evaluate the theoretical findings of the previous sections. The three dimensional compressible Euler equations are given by The state vector and the components of the block vector flux, ↔ f , are given by where the conservative variables are the density ρ, the momentum ρ u = [ρu 1 , ρu 2 , ρu 3 ] T and the total energy E. In order to close the system, we assume an ideal gas such that the pressure is defined as where γ is the adiabatic exponent. We choose γ = 1.4 in the following experiments. The system (3.1) is investigated in the domain = [x min , x max ] 3 . At initial time t = 0 the domain is divided in K non-overlapping and conforming cartesian hexahedral elements e κ (0), κ = 1, . . . , K . For each element e κ (0), κ = 1, . . . , K , the temporal distribution of a grid point is given by (3.5) where L := x max − x min . In Fig. 2, we show a slice through a three dimensional mesh with K = 16 3 elements at initial time and at its maximal distortion. The mesh velocity is calculated by exact differentiation of Eq. (3.5). Note that the formula (3.5) is a common way to construct a deforming domain. For instance, similar formulas were used for the DG methods in [37,64] and the collocation method in [66]. Furthermore, the five stage fourth order low-storage two-register explicit RK method (RK4 (3) where h κ (t n ) is the minimum element size of e κ (t n ), C CFL ∈ (0, 1] and λ max is the largest advective wave speed at the current time level traveling in either the x 1 , x 2 , x 3 -direction.

Experimental Convergence Rates
In this section, we verify the high-order approximation of the moving DGSEM (2.71). For this purpose, we investigate the domain = [−1, 1] 3 and apply the method of manufactured solutions. Thus, we assume a solution of the form We plug solution (3.7) into the Euler system and compute the residual using a computer algebra system. This term is used as a source term in our convergence tests. We note that this term is handled and discretized as a solution independent part in the numerical computation.
We run the convergence test with periodic boundary conditions. Furthermore, the moving mesh DGSEM (2.71) is applied with the flux function in Appendix B.1 as volume and surface flux. In addition, the surface flux is stabilized by the dissipation operator in Appendix B.3. Besides using the grid point distribution given in (3.5), we also compute static reference solutions, by setting the grid velocity to zero. In this case, the moving mesh DGSEM (2.71) degenerates to the split form DGSEM for static meshes [23,24].
In Table 1, we list the experimental order of convergence (EOC) and L 2 errors for the conservative variables that we obtain for polynomials with odd degree N = 3 on a static mesh (top) and on a moving mesh (bottom). To calculate the L 2 norm, we interpolate the polynomial solution to a higher degree (at least twice the degree of the solution) and perform integration on that higher degree. The convergence rates on the moving mesh are not as good as on a static mesh, which can be justified by the high distortion in the mesh from the grid point distribution formula (3.5). However, with an increasing number of elements the same convergence rates as on a static mesh are almost reached. Moreover, the experimental order of convergence (EOC) and L 2 errors for the conservative variables that we obtain for polynomials with even degree N = 4 are listed in the Table 2. We observe a similar behavior as for the odd degree N = 3. This indicate the high-order approximation properties of the moving mesh DGSEM.

Entropy Analysis Validation
The three dimensional Euler equations (3.1) are equipped with the entropy/entropy flux pairs  The moving mesh DGSEM is used with N = 3 on a static mesh (top) and on a moving mesh (bottom) with the grid point distribution (3.5)

Table 2
Experimental order of convergence (EOC) and L 2 errors at time T = 5 for the Euler manufactured solution test (3.7) The moving mesh DGSEM is used with N = 4 on a static mesh (top) and on a moving mesh (bottom) with the grid point distribution (3.5) where ς = log pρ −γ . We are interested in the behavior of the discrete entropy conservation error (3.9) whereS (·) is computed by (2.78). Therefore, we investigate the inviscid Taylor-Green vortex (TGV) test case [56] in the domain = [0, 2π] 3 . The inviscid TGV can be a challenging test case regarding the robustness of a numerical scheme, partly because the dynamics produce arbitrarily small scales. The flow field is thus by design under-resolved, which makes it a suitable test case to investigate the entropy conservation properties of the scheme. The TGV evolves from the initial data (3.10) To render the simulation close to incompressible, the Mach number M 0 = 1 √ γ p 0 is set to 0.1 by adjusting the pressure correspondingly. We run the simulation with K = 16 3 elements and periodic boundary conditions. The final time is chosen to be T = 13. Furthermore, we apply the flux function in Appendix B.1 to compute the derivative projection operator (2.64). To calculate the discrete integral entropy, the SBP mass matrices are used directly. In Fig. 3 we present a log-log plot of the entropy conservation error for N = 3, 4. We note that the flux in Appendix B.1 is used as surface flux without a dissipation term in these computations, rendering the semi-discrete discretization fully entropy conserving. As expected, we observe the reduction of the remaining entropy conservation error according to the order of the RK method for decreasing CFL numbers. In Fig. 4 the temporal evolution of the entropy conservation errors S (T ) is given. The CFL number is set to C CFL = 0.125 and polynomial degrees N = 3 and N = 4 are used. We observe that the entropy conservation error S (T ) is constant in time (dashed line) when the flux in Appendix B.1 is used without a dissipation term. This indicates the entropy conservation in the TVG test case. On the other hand the entropy conservation error S (T ) is decreasing in time (solid line) when the surface flux is stabilized by the dissipation term in Appendix B.3. Thus, the moving mesh DGSEM is an entropy stable scheme in this test case. These observations agree with the results in Theorem 2.4 and Corollary 2.6.

Robustness Test
As has been stated in Sect. 3.2 and noted in literature [50,62], the inviscid TGV is a notoriously challenging test case for the stability of the numerical scheme. While for lower polynomial degrees calculations may be possible, high-order simulations are known to crash even if aliasing-reducing methods like polynomial dealiasing are used [50]. Thus, we use the TGV test case (3.10) to demonstrate the increased robustness of the entropy stable moving mesh DGSEM. To do so, we run the simulation up to T = 13 using a polynomial degree of N = 7 on three different meshes employing K 1 = 14 3 , K 2 = 19 3 and K 3 = 26 3 elements. These cases correspond to the most restrictive simulations from [50]. Again, the point distribution given in (3.5) is used. We use the flux function in Appendix B.1 as volume and surface flux and stabilize the surface flux by the dissipation operator in Appendix B.3.  Using the entropy stable moving DGSEM, we were able to run all simulations until final time. This shows that the consistent dissipation operators in combination with the entropy conservative volume fluxes can lead to superior stability properties. We note that simulations without dissipative surface fluxes crash before reaching the final time for higherorder simulations (N ≥ 3). This highlights the role that entropy conservation plays in the stabilization of challenging numerical problems.

Free Stream Preservation Validation
We consider the Euler equations (3.1) on the domain = [0, 2π] 3 with the initial data Table 3 Free stream preservation test for N = 3 (top) and N = 4 (bottom) The L ∞ errors measure the difference between the initial data (3.11) and the numerical solution at time T = 20 The entropy stable DGSEM is applied with the flux function in Appendix B.1 as volume and surface flux as well as the dissipation operator in Appendix B.3 to stabilize the surface flux. We apply K = 16 3 elements, the formula (3.5) to describe the displacement of the mesh points and periodic boundary conditions are used in the simulation. Furthermore, the final time is set to T = 20. In Table 3, we present the L ∞ errors between the initial data (3.11) and the numerical solution at time T = 20 for polynomials of degree N = 3 (top), N = 4 (bottom) and different CFL numbers C CFL . The errors are computed by super sampling the polynomial solution, at least doubling the amount of nodes per direction from the underlying solution. We observe that the errors are close to zero and vary slightly for the different CFL numbers. These results indicate the compliance of the free stream preservation property.

Conclusions
In this work a moving mesh DGSEM to solve non-linear conservation laws has been constructed and analyzed. The semi-discrete method is provably entropy stable and the free stream preservation property is satisfied for each explicit s-stage Runge-Kutta method. The moving mesh DGSEM has been presented for three dimensional conservation laws. The derivatives in space are approximated with high-order derivative matrices which are SBP operators. Furthermore, the split form DG framework [23,24] has been used to avoid aliasing in the discretization of the volume integrals. In addition, two-point flux functions with the generalized entropy condition (2.67) are used in the split form DG framework. These modules in the spatial discretization are the basis to prove that the moving mesh DGSEM is an entropy stable scheme, when periodic boundary conditions are used. Non-periodic boundary conditions require the construction of suitable dissipative boundary conditions to enforce the entropy stability. The discrete entropy analysis requires the assumptions that the time derivatives can be evaluated exactly and that properties like positivity preservation (of the water height, density or pressure) are satisfied on the discrete level. Operations like the integration-by-parts formula are mimicked by the SPB operators on the discrete level.
The three dimensional Euler equations have been considered to verify the proven properties of the moving mesh DGSEM in our numerical experiments. We presented convergence tests for smooth test problems to verify that the split form DG framework provides also on a moving mesh a high-order accurate approximation. Furthermore, the numerical robustness tests in the Sect. 3.3 emphasize the relevance of the entropy stable DGSEM, since the method was able to run the challenging inviscid TGV test case until final time.
since it holds the identity

B Entropy Stable Moving Mesh Euler Fluxes
We present entropy stable Cartesian fluxes G EC l , l = 1, 2, 3, for the compressible Euler equations (3.1) equipped with the entropy/entropy flux pairs (3.8). Then the entropy variables are given by and the entropy functionals are given by

B.1 Entropy Conservative Euler Flux Based on the Flux in [8]
In the following the logarithmic mean {{·}} log will be used. For two positive states a − and a + , the logarithmic mean is defined by The flux (B.7) is consistent with f 1 and symmetric. In particular, Chandrashekar proved that The state function (B.4) and the flux function (B.7) are used to construct the flux The flux (B.9) is consistent with g 1 = f 1 − ν 1 u, symmetric and it follows The fluxes (B.11), (B.12) are consistent with g 2 = f 2 − ν 2 u, g 3 = f 3 − ν 3 u, symmetric and satisfy Ranocha constructed in his PHD thesis [55] another KEPEC numerical flux for the compressible Euler equations. We proceed as in the Appendix B.1 and use the state (B.4) and Ranocha's flux to construct the following two-point flux functions The flux functions (B.15) are consistent with g 1 , g 2 , g 3 , symmetric and satisfy

B.3 Matrix Dissipation Term for the Euler Flux
We will use the Euler fluxes from the previous Appendices B.1 and B.2 with the dissipation operators form Winters et al. [63]. In the following, the matrices to construct the entropy stable dissipation operators (2.96) are listed. The average components of the dissipation term in the x 1 -direction are given by

C Proofs of Entropy Conservation for Advection Terms
In this section, we apply the following identities which result from the properties of the SBP operator Q N i, j=0 are generic nodal values. These identities can be proven in a similar way as the discrete split forms in Lemma 1 in [23]. Thus, we skip a proof in this paper.