A New Family of Thermodynamically Compatible Discontinuous Galerkin Methods for Continuum Mechanics and Turbulent Shallow Water Flows

In this work we propose a new family of high order accurate semi-discrete discontinuous Galerkin (DG) finite element schemes for the thermodynamically compatible discretization of overdetermined first order hyperbolic systems. In particular, we consider a first order hyperbolic model of turbulent shallow water flows, as well as the unified first order hyperbolic Godunov–Peshkov–Romenski (GPR) model of continuum mechanics, which is able to describe at the same time viscous fluids and nonlinear elastic solids at large deformations. Both PDE systems treated in this paper belong to the class of hyperbolic and thermodynamically compatible systems, since both satisfy an entropy inequality and the total energy conservation can be obtained as a direct consequence of all other governing equations via suitable linear combination with the corresponding thermodynamic dual variables. In this paper, we mimic this process for the first time also at the semi-discrete level at the aid of high order discontinuous Galerkin finite element schemes. For the GPR model we directly discretize the entropy inequality and obtain total energy conservation as a consequence of a suitable discretization of all other evolution equations. For turbulent shallow water flows we directly discretize the nonconservative evolution equations related to the Reynolds stress tensor and obtain total energy conservation again as a result of the thermodynamically compatible discretization. As a consequence, for continuum mechanics the new DG schemes satisfy a cell entropy inequality directly by construction and thanks to the discrete thermodynamic compatibility they are provably nonlinearly stable in the energy norm for both systems under consideration.


Introduction
Hyperbolic systems find a wide range of applications in science and engineering, like the Euler equations of compressible gasdynamics in aerospace engineering, the magnetohydrodynamics (MHD) equations in plasma and astrophysics, the Maxwell equations in computational electromagnetism, the equations of nonlinear hyperelasticity which describe the dynamics of nonlinear elasto-plastic solids under large deformations, or the shallow water equations in geophysics and environmental engineering. A common aspect of these equations is that they can all be written under the form of a symmetrizable and thermodynamically compatible hyperbolic (SHTC) system.
In his groundbreaking work An interesting class of quasilinear systems, Godunov discovered the link between symmetric hyperbolicity in the sense of Friedrichs [44] and thermodynamic compatibility, see [52]. The idea of Godunov was rediscovered independently 10 years later by Friedrichs and Lax [45]. Relevant contributions to the field can also be found in the works of Boillat [13] and Ruggeri [88]. In a series of subsequent papers Godunov and Romenski and collaborators were able to show that the class of SHTC systems covers a rather wide range of known mathematical models, including magnetohydrodynamics [53], nonlinear hyperelasticity [54], compressible multi-phase flows [84,86] and also fluid and solid mechanics in special and general relativity, see [56,85]. More complex continuum mechanics with torsion was included in the class of SHTC systems in [78] and very recently a link between the class of SHTC systems and Hamiltonian mechanics was discovered, see [76]. For a general presentation of the SHTC framework the reader is referred to [83] and [55]. All SHTC systems can be derived from an underlying variational principle and the total energy potential appearing in the Lagrangian has therefore a privileged role in the SHTC framework. This special role of the total energy potential is also justified by the connection with Hamiltonian mechanics. A common point of view adopted in the framework of most SHTC systems is therefore that the entropy density is seen as a primary evolution variable, while the conservation law for the total energy density is seen as the consequence, since total energy conservation can be obtained as an extra conservation law after a suitable linear combination of all the other evolution equations with the so-called thermodynamic dual variables or main field variables [88], which are the partial derivatives of the total energy density with respect to the conservative variables.
Since the seminal work of Tadmor [92] on entropy compatible schemes, many schemes have been developed in order to obtain thermodynamic compatibility on the discrete level according to the ideas of Friedrichs and Lax [45], i.e. discretizing the energy conservation law directly and obtaining the entropy inequality as a consequence of the other equations. For recent works on high order entropy-compatible finite volume and discontinuous Galerkin finite element schemes for systems of conservation laws the reader is referred to [27,28,30,38,43,46,48,60,70,81,82,87,89] and references therein. Entropy-compatible schemes for non-conservative hyperbolic equations were introduced, for example, in [2,42]. Most of the above-mentioned schemes rely on the framework introduced by Tadmor in [92], while a completely different and very simple and general framework for the construction of com-patible numerical schemes that satisfy extra conservation laws at the discrete level was very recently forwarded in [1,2,[4][5][6]. For the sake of completeness, we also would like to mention recent developments in the context of thermodynamically compatible schemes for Lagrangian hydrodynamics, see [8,25,72].
Up to now, there are only very few numerical schemes which try to construct a dual algorithm to the framework of Tadmor and in line with the philosophy of SHTC systems, namely numerical schemes that discretize the entropy inequality directly and which obtain the total energy equation as a natural consequence of the thermodynamically compatible discretization of all other equations. First attempts have been documented recently in [3,[20][21][22], where thermodynamically compatible finite volume schemes were introduced that obtain the total energy conservation as a consequence of all the other equations, thus creating the missing dual algorithms to the entropy-compatible schemes of Tadmor. For the construction of a thermodynamically compatible flux for the underlying Euler and shallow water subsystems our method makes use of the Godunov parametrization in terms of a generating potential.
It is the aim of the present paper to extend the above mentioned thermodynamically compatible finite volume schemes to the discontinuous Galerkin finite element framework, i.e. we want to construct new thermodynamically compatible DG schemes that directly discretize the entropy inequality and which obtain discrete total energy conservation as a consequence. Therefore, our new DG schemes will satisfy a cell entropy inequality by construction and can be proven to be nonlinearly stable in the energy norm. We want to stress that is is not the aim of this paper to develop better numerical schemes compared to existing ones, which directly discretize the total energy conservation equation and which obtain the entropy inequality as a consequence. In this paper, instead, we want to promote a radically new concept in the field of entropy-compatible DG schemes. For an introduction to high order DG schemes, the reader is referred to the seminal papers of Cockburn and Shu [32][33][34][35][36][37], while the first proof of nonlinear L 2 stability and verification of a discrete cell entropy inequality for DG schemes was presented for nonlinear scalar conservation laws in [63].
We also would like to clearly point out the two main limitations of the DG method presented in this paper: first, all mathematical proofs rely on the assumption of an exact calculation of all integrals, as in [63]. While at least in principle arbitrary accuracy can be achieved via suitable adaptive numerical quadrature, this assumption is rather restrictive for practical purposes. In the actual implementation we have always used N + 1 Gauss-Legendre quadrature points to approximate the integrals for DG schemes of polynomial approximation degree N . Future work will include a generalization of the HTC DG schemes presented in this paper to account also for numerical quadrature errors. Second, we have not yet developed suitable thermodynamically compatible limiters, but we rather rely on a simple artificial viscosity approach to stabilize our scheme in the presence of shock waves or steep gradients. Future research concerning the two previous points will be needed, but is out of scope of the present work.
The rest of this article is organized as follows. In Sect. 2 we present the two governing PDE systems that are discretized in this paper, namely the unified first order hyperbolic GPR model of continuum mechanics [22,40,77,83] and the first order hyperbolic model of unsteady turbulent shallow water flows introduced and studied by Gavrilyuk et al. in [49,50,62]. In Sect. 3 we first present the new semi-discrete DG schemes in one space dimension, for pedagogical reasons and to facilitate the reading, and subsequently we extend it also to multiple space dimensions. A set of numerical results obtained with the new HTC DG schemes is presented for both governing PDE systems in Sect. 4. The paper is rounded-off by some concluding remarks and an outlook to future work in Sect. 5.

Godunov-Peshkov-Romenski (GPR) Model of Continuum Mechanics
We consider the following first order hyperbolic model of continuum mechanics regularized with vanishing viscosity terms and which goes back to the work of Godunov [52], Godunov and Romenski [54,55,83] and Peshkov and Romenski, see [40,77]: is a vanishing viscosity and the nonnegative entropy production term due to the viscous terms is given by The positivity of this term comes from > 0 and since we assume a positive temperature T > 0 and that the Hessian of the total energy potential is at least positive semi-definite, H i j := ∂ 2 q i q j E ≥ 0. Throughout this paper, we use the notations ∂ p = ∂/∂ p and ∂ 2 pq = ∂ 2 /(∂ p∂q) for the first and second partial derivatives w.r.t. generic coordinates or quantities p and q, which may also be vectors or components of a vector. Furthermore, we make use of the Einstein summation convention over repeated indices. Last but not least, in some occasions we also use bold face symbols in order to denote vectors and matrices, e.g. q = {q i } and A = {A ik }, and so on. In the above model the four contributions to the total energy density are with the components of the metric tensor G and its trace-free partG given by G ik = A ji A jk andG ik = G ik − 1 3 G mm δ ik . The total energy density associated with the Euler subsystem is E 12 = E 1 + E 2 . The vector of thermodynamic dual variables reads The pressure is defined as p = ρ ∂ ρ E +ρv i ∂ ρv i E +ρ S ∂ ρ S E −E = ρ 2 ∂ ρ E, the stress tensors due to shear stress and thermal stress are, respectively, while the heat flux vector is given by The total energy flux F k can be decomposed into a contribution from the Euler subsystem (black terms) and the remaining terms (red terms) as F k = F 12 k + F 34 k , with Note that for our convenience, we use the opposite sign in the definition of the stress tensor compared to the generally accepted notation. Furthermore, θ 1 (τ 1 ) > 0 and θ 2 (τ 2 ) > 0 are two algebraic functions of the state vector q and the positive relaxation times τ 1 > 0 and τ 2 > 0: with ρ 0 and T 0 being some reference density and temperature. It is easy to check that (1f) is a consequence of (1a)-(1e), i.e.
In [40] a formal asymptotic analysis of the model (1a)-(1f) was carried out, revealing that in the stiff limit the stress tensor σ ik and the heat flux h k tend to i.e. when the relaxation times τ 1 , τ 2 → 0, the Navier-Stokes-Fourier equations are retrieved, with effective shear viscosity μ = 1 6 ρ 0 c 2 s τ 1 and heat conductivity κ = ρ 0 T 0 c 2 h τ 2 .

Turbulent Shallow Water Flows (TSW)
The second model considered is the following overdetermined hyperbolic model for turbulent shear shallow water flows (TSW) in multiple space dimensions, which has been recently proposed by Gavrilyuk and Ivanova et al. in [50] and which was also applied and studied in [12,29,62]. In this paper, we employ the reformulation proposed in [21] in terms of a new evolution variable Q that allows to rewrite the Reynolds stress tensor P as P = QQ T . Note that the original Reynolds stress tensor P is symmetric and positive definite, while Q is not symmetric. There are three reasons for such a decomposition: first, the interesting analogy with the GPR model of continuum mechanics, where the symmetric positive definite metric tensor G = A T A is written in terms of the non-symmetric distortion field A (the inverse deformation gradient for purely elastic materials), see [21] and eqn. (16) below; second, there is a purely numerical advantage, since the evolution of Q instead of P allows the discrete trace of P to be naturally non-negative, while this is non-trivial for high order DG schemes when evolving P directly, see also [21]; third, the total energy potential is a quadratic function in terms of Q, while it is not in terms of P. The quadratic dependence on Q is important for the construction of thermodynamically compatible schemes, as those presented in [21] and in this paper. For more details concerning the physical interpretation of the PDE for Q and the various motivations why to choose this evolution variable rather than the original Reynolds stress tensor P, see [21]. With P ik = Q im Q km written in terms of the new evolution variable Q and the notation ∂ m = ∂/∂ x m the turbulent shallow water system can be rewritten as an overdetermined PDE system as follows: with the evolution variables h = h(x, t), hv = hv(x, t), Q = Q(x, t), the Reynolds stress tensor P = P(x, t), and the total energy Here, E 12 = E 1 + E 2 is the total energy potential of the shallow water subsystem in (11a)-(11b) (black terms) and E 3 is the total energy associated with the new object Q ik (red terms).
In what follows, we will again decompose the inviscid part of the total energy flux in (11d) as with the abbreviations that will be used later. Here, F 12 corresponds to the energy flux related to the shallow water subsystem (black terms) and F 34 to the energy flux related to the work of the Reynolds stress tensor (red terms). The production term, ik , which is needed to achieve the consistency of (11a)-(11c) with the total energy conservation law (11d), reads In (14) the vector q = q i = (h, hv i , Q ik ) indicates the vector of primary state variables and ∂ 2 q i q j E is the Hessian matrix of the total energy potential with respect to these state variables. One can show that the Hessian matrix is positive definite for sufficiently small turbulent kinetic energy Q i j Q i j compared to gh, namely Q i j Q i j < gh, which is a rather mild assumption, see [21] for details. In [21] it was also shown that the energy conservation law (11d) is a consequence of Eqs. (11a)-(11c). Concerning the entropy inequality associated with the turbulent shallow water equations, in [21,62] it was shown that the following entropy inequality is a consequence of system (11): with δ ik the usual Kronecker symbol and |Q| = det(Q) denoting the determinant of Q. However, obtaining a discrete analogue of (15) is clearly out of scope of the present paper, since it would require a provably compatible discretization of the Jacobi identity of the time derivative of the determinant of a matrix. Note that the system (11) contains non-conservative terms that do not vanish across shock waves, see [62]. This poses many challenges from a theoretical as well as from a numerical point of view. As already pointed out in [21,62] the key ingredient for a successful numerical discretization of (11) is the discrete consistency of the numerical scheme with the total energy conservation law. In [62] this discrete consistency was imposed via a suitable wave splitting approach that made use of total energy conservation in each sub step of the split scheme. In [21] total energy conservation was instead achieved either via a provably thermodynamically compatible discretization, as used in this paper, or via a suitable rescaling of Q at the end of each time step, such as to enforce discrete energy conservation. For a more detailed discussion on this topic, see [21]. Note that the governing PDE for Q (11c) is formally very similar to the governing PDE for the distortion field A ik in nonlinear hyperelasticity, (1d), which after application of the product rule reads As one can easily see, the order of the matrix-product in the third term on the left hand side of (11c) and (16) is exchanged.

The Godunov Form of Hyperbolic Systems of Conservation Laws
The Godunov form, [52], of the inviscid Euler subsystem in (1) and of the shallow water subsystem in (11) (black terms) reads The so-called generating potential L 12 , which is the Legendre transform of the total energy density related to the Euler/shallow water subsystem E 12 = E 1 + E 2 , is defined as It is easy to check that with the above parametrization of the hyperbolic system according to Godunov, the extra conservation law for the total energy density takes the form with the related total energy flux F 12 k .

General Formulation
Both PDE systems presented previously, i.e. the GPR model (1) and the hyperbolic model for turbulent shallow water flows (11) can be cast into the following general form with the extra conservation law for the total energy density Here, the flux tensor f k (q) is related to the terms that directly fall into the formalism of Godunov (black terms), the flux tensor h k (q) and the nonconservative product B k (q)∂ k q include the terms that are not directly contained in the Godunov formalism (red terms), the dissipation is ∂ m ( ∂ m q) with the associated entropy production P (blue terms) and the relaxation source terms are denoted by S(q) (green terms). In the extra conservation law (22) the total energy flux is denoted by F k = F k (q).
Since p = ∂ q E and therefore p · ∂ t q = ∂ t E, for thermodynamic compatibility of the full system (21) with (22) the following identity must hold: The above condition can be divided into two separate equalities, which mean compatibility of the flux terms and of the non-conservative products with the corresponding contribution to the total energy flux. Moreover, we also have i.e. compatibility of the dissipation terms with the production term and is required in order to have compatibility of the algebraic source terms.

Review of Thermodynamically Compatible Finite Volume Schemes
In this section we first briefly recall the new class of thermodynamically compatible HTC finite volume schemes forwarded in [21,22] for the GPR model and for the hyperbolic model of turbulent shallow water flows and in [20] for the MHD equations. All these schemes obtain the total energy conservation law as a consequence of all other equations. Recalling these schemes here is useful since their thermodynamically compatible numerical flux and the appropriate discretization of the fluctuations will be also the basis of the high order DG schemes proposed in this paper. Given a spatial control volume denoted by with circumcenter x , and denoting one of its neighbors by r and the common edge by ∂ r , n r = (n r 1 , n r 2 ) T being the outward pointing unit normal vector to the face ∂ r , with n r = −n r , and denoting by N the set of neighbors of cell , then the thermodynamically compatible finite volume schemes [20][21][22] read as follows: The thermodynamically compatible numerical flux in normal direction F q , q r · n r is based on a discrete Godunov formalism, making use of the generating potential L 12 and a path integral along a straight line segment path in the main field (thermodynamic dual) variables, which leads to the so-called p-scheme proposed in [21], with the path given by The path integral is calculated numerically at the aid of a Gauss-Legendre quadrature of sufficient accuracy. It is easy to check, see [20][21][22], that with the Godunov parametrization of the physical flux f k = ∂ p (v k L 12 ) in terms of the generating potential L 12 = p · q − E 12 the numerical flux (29) satisfies the compatibility property and since F 12 , after subtracting p r · f r k − p · f k n r k from both sides of the above equation, we obtain the property in terms of the fluctuations F r − f k n r k and f r k n r k − F r related to the numerical flux in normal direction. Moreover, the flux containing the numerical viscosity reads with an associate non-negative entropy production term that will be specified later. If present in the mathematical model, the algebraic source term must obey a pointwise compatibility condition as follows while the fluctuations related to the remaining part of the energy potential (red terms) satisfy with n r = −n r . The discrete fluctuation terms, R q , q r · n r , are specific for each model and are detailed below. In this paper we use the terminology fluctuation for both, flux differences, as well as for the jump terms that arise at the element boundaries due to the non-conservative terms in the PDE system (21), see [27,75]. The fluctuations are also related to the wave propagation form of the scheme, see [68,69]. GPR model. Following [22], the thermodynamically compatible discretization of the fluctuations R q , q r · n r for the GPR model (1) reads: where and P q , q r = 0, 0, r , 0, 0 T , is the non-negative entropy production term due to the numerical viscosity with Turbulent shallow water flows. According to [21] the thermodynamically compatible discretization of the fluctuations R q , q r · n r for the turbulent shallow water model (11) reads: where and P q , q r = 0, 0, r ik T , is the non-negative entropy production term due to the numerical viscosity with In both cases the Roe matrix of the Hessian ∂ 2 qqẼ r must satisfy the Roe property It can be easily obtained via the path integral where this time a segment path in q variables is used: Throughout this paper all path integrals are computed numerically at the aid of a Gauss-Legendre quadrature formula of sufficient accuracy. In practice, we found that 3 quadrature points are enough, see [22] for a detailed analysis on the influence of the accuracy of the numerical quadrature on total energy conservation.

Thermodynamically Compatible DG Schemes in One Space Dimension
In one space dimension, the computational domain, , is covered by non-overlapping elements of the form and the discrete solution is expanded in terms of nodal basis functions ϕ m (x) as withq m (t) the time dependent degrees of freedom and N the maximum polynomial degree of the basis functions. In this paper, we use Lagrange interpolation polynomials passing through the Gauss-Legendre quadrature nodes. Consequently, we have the discrete main field and the discrete total energy density To ease notation, we will make use of the following abbreviations: Multiplication of (21) by a test function ϕ k , integration over a cell , substitution of (48), integration by parts and addition of the jump terms ∂ x ϕ ± 1 2 k V ± 1 2 with positive sign on the right hand side lead to the DG scheme where the different terms are defined as described below.
According to [20][21][22] and as shown above for the finite volume case, the numerical flux is obtained by a path integral following the ideas on path conservative methods presented in [26,75]: with n = (1, 0, 0), in 1D, and the path given by The point values p are given by the point-wise evaluation p the boundary-extrapolated values of the discrete numerical solution q h . The numerical flux, by construction, verifies the following Roe-type property, see [20][21][22]: The viscous numerical flux is inspired by the work of Gassner, Lörcher and Munz [47] and reads We stress that in our DG scheme above the jump terms V + 1 2 have the opposite sign compared to the ones derived in [47] via the usual integration by parts back and forth. This choice of the sign is necessary in order to prove semi-discrete total energy conservation and thus nonlinear stability in the energy norm, see the additional comments in the proof of Theorem 2. They read For the GPR model, the discrete entropy production term related to the viscous terms is while for the turbulent shallow water equations we have P k = (0, 0, TSW For both models, the fluctuations are defined according to (36) and (41), respectively, which in one space dimension reduces to with n = (1, 0, 0). They verify the discrete compatibility relation

Thermodynamically Compatible DG Schemes in Multiple Space Dimensions
In multiple space dimensions, the computational domain is discretized via non-overlapping ]. The discrete solution for the vector of conservative variables reads as usual As in the one-dimensional case we employ nodal basis functions, defined as the tensor products of the one-dimensional basis functions described in the previous section. The maximum polynomial approximation degree per dimension is again denoted by N and N in (60) denotes the number of degrees of freedom per element, i.e. (N + 1) d for a tensor-product basis in d space dimensions. In multiple space dimensions, the thermodynamically compatible DG scheme for the hyperbolic system (21), denoted by HTC DG scheme in the following, can be derived via multiplication of (21) with a test function ϕ k (x), integration by parts of the flux terms and via the use of suitable thermodynamically compatible fluxes and fluctuations defined later: Here, q h , q r h are the boundary extrapolated values of the DG solution computed at the left and right sides of the element boundary, respectively. To compute the thermodynamically compatible numerical fluxes F(q h , q r h ) · n and the fluctuations R(q h , q r h ) · n we can make use of those already introduced for the finite volume scheme (29), (36), (41), which satisfy the compatibility relations (31) and (35). The viscous fluxes read while the viscous jump terms are given by We emphasize again that the term in (61) including the jump V q , q r has the opposite sign compared to the one given in [47] in order to prove thermodynamic compatibility, see the proof of Theorem 2. The discrete entropy production term, while for the TSW model it reads P k = (0, 0, TSW Theorem 1 (Cell entropy inequality) For the GPR model the HTC DG scheme (61) satisfies the cell entropy inequality Proof To obtain a cell entropy inequality for the GPR model, we employ the test function ϕ k = 1 in (61) that, together with the notation and (64), yields Substituting (36), we get At the aid of Gauss' theorem and thanks to π ≥ 0 and ∂ 2 qq E ≥ 0, we obtain the sought cell entropy inequality: for vanishing boundary fluxes.
Proof Using the Godunov parametrization of the flux f m = ∂ p (v m L) with the shorthand notation L := L 12 and summing up all equations in (61) after multiplication byp i m yields The compatibility condition (25) and the point-wise compatibility of the source terms, (27), Using (24), the notation F r = F(q h , q r h ) · n and adding and subtracting 1 2 Next, we add and subtract 1 2 ∂ i (v m L) r · n m d S and use the compatibility condition of the fluctuations, (35), Rearranging terms and using property (31) of the thermodynamically compatible flux, we get We now address the numerical dissipation terms. By adding and subtracting 1 where the second equality follows from the definition of the dissipative terms (62)- (63) and the last one is obtained introducing the production term and making use of the Roe property of the Hessian ∂ 2 ppL r . The two terms indicated by the red dashed lines cancel only due to the particular choice of the opposite sign of the jump terms related to V in (61) compared to the sign given in [47]. Since the main objective of this paper is the construction of provably thermodynamically compatible DG schemes that leads to semi-discrete total energy conservation, the sign is dictated by the present proof. Substitution of the former relation into (69) leads to where the sum of the second, third and fourth terms is a consistent approximation of the total energy flux. Integrating over the domain , we obtain nonlinear stability in the energy norm since the sum of all internal fluxes cancels and the boundary fluxes vanish.

Numerical Results
In this section, we present several numerical test cases aiming at assessing the proposed semidiscrete HTC DG schemes for both the GPR model and the hyperbolic turbulent shallow water system. As time integrator the classical fourth order Runge-Kutta method is used for all test problems shown below. Besides, the time step is set according to the CFL-type condition with h = min( x, y) the characteristic mesh spacing, |λ max | the maximum absolute value of the eigenvalues in the domain and CFL < d, where d is the number of space dimensions.

GPR Model of Continuum Mechanics
In all numerical tests carried out in the fluid limit of the GPR model the relaxation time τ 1 is obtained from the relation μ = 1 6 ρ 0 c 2 s τ 1 for given shear sound speed c s , reference density ρ 0 and dynamic viscosity μ. If not stated otherwise, the artificial viscosity is by default set to = 0, the reference density is set to ρ 0 = 1 and the specific heat at constant volume is set to c v = 1. We furthermore set γ = 1.4 for all tests.

Numerical Convergence Study
The order of accuracy of the new HTC DG scheme is verified experimentally for the Euler subsystem, i.e. for the black terms in (1) at the aid of the well-known isentropic vortex problem, see [61]. The parameters for the GPR model are set to c s = 0, c h = 0 and the artificial viscosity is set to = 0. The problem is solved until a final time of t = 0.25 in a periodic domain = [0, 10] 2 . Since the flow is isentropic, the entropy is constant and therefore the corresponding velocity, temperature, density and pressure profiles are with r 2 = (x − 5) 2 + (y − 5) 2 and the vortex strength ε = 5. The above vortex is a steady solution of the Euler equations, hence the initial condition is equal to the exact solution of the problem also for all later times. The numerical convergence study of the HTC DG schemes is carried out with different polynomial approximation degrees N on a sequence of successively refined Cartesian meshes composed of N x × N y elements. The L 2 error norms computed at the final time are reported in Table 1, together with the numerically observed convergence rates of the scheme for the density ρ, the momentum density ρv 1 and the entropy density ρ S.
The numerical convergence rates are optimal for the density (N +1), while they are not for the momentum density. Instead, the entropy density reaches orders between 2N and 2N + 2. The explanation for this interesting observation is the following: the present test problem is isentropic and for = 0 the only mechanism that generates entropy in the HTC DG scheme is the jump term 1 2 s max (q + h − q − h ) in the numerical viscosity flux G. But since the jumps tend to zero with order between N + 1 2 to N + 1 and the production term k in the entropy inequality is quadratic in the jump, for the proposed HTC DG scheme we actually expect twice the convergence order for the entropy density for all isentropic flows. This seems to be indeed a very interesting feature of our new HTC DG scheme, which is directly based on the discretization of the entropy inequality (1c) rather than on the total energy conservation law (1f), unlike standard DG schemes.

Simple Shear Motion in Solids and Fluids
Next, we simulate the time evolution of a simple isolated shear layer in one space dimension. In the fluid limit of the model, i.e. for τ 1 1, we set κ = μ and the reference solution for v 2 is simply given by the exact solution of the first problem of Stokes for the incompressible Navier-Stokes equations, see e.g. [11,15,24,40], and which reads v 2 (x, t) = v 0 erf 1 2 x √ νt (72) with ν = μ/ρ 0 . In the solid limit of the GPR model a reference solution can be obtained by solving the system (1) on a very fine mesh of 10000 cells using a classical second order accurate MUSCL-Hancock type TVD finite volume scheme, see [96] for details. In all cases we set the artificial viscosity to = 10 −6 . The numerical results obtained with the new HTC DG scheme are shown in Fig. 1 and in all cases an excellent agreement between numerical and reference solutions can be observed.

Riemann Problems
In this section, we apply the new HTC DG scheme to five Riemann problems in the domain = [−0.5, +0.5] with left and right initial states given in Table 2 and with initial discontinuity located in x c = 0, if not stated otherwise. We consider test cases for both, the full GPR model (1) (RP4 and RP5), as well as the Euler subsystem (RP1, RP2, RP3), i.e. just the black terms in (1). The exact Riemann solver for the Euler equations can be found in the wellknown textbook [96]. For the full GPR model we generate a numerical reference solution by solving (1) on a very fine mesh of 128000 control volumes using a classical second order TVD finite volume scheme of the MUSCL-Hancock type, see [96]. For the computation of the reference solution, we solve the total energy conservation law (1f) instead of the entropy inequality (1c), while in the HTC DG scheme we solve the entropy inequality (1c) instead of Table 2 Initial states left (L) and right (R) for density ρ, velocity v = (v 1 , v 2 , 0) and pressure p for Riemann problems RP1-RP3 (Euler subsystem) and for RP4-RP5 (full GPR model) the energy Eq. (1f). As such, the reference solution is really obtained in a completely different manner compared to the numerical scheme proposed in this paper. For the last two Riemann problems, RP4 and RP5, we define the initial conditions for the distortion field A and for the specific thermal impulse J to A = I and J = 0. Furthermore, we set c s = c h = 1. For RP4 we furthermore choose the relaxation times so that μ = λ = 10 −5 and for RP5 we simply set τ 1 = τ 2 = 10 20 . The artificial viscosity is set to = 10 −5 in all cases. In Figs. 2, 3, 4 and 5 we compare the reference solutions with the numerical results obtained with the new HTC DG scheme proposed in this paper. The employed mesh resolution is provided for each test case in the corresponding figure caption and the polynomial approximation degree for the HTC DG scheme is set to N = 9. For all Riemann problems we can observe an excellent agreement between numerical solution and reference solution.

Viscous Shock Wave
Here we apply our new HTC DG scheme to a stationary viscous shock wave with a shock Mach number M s = 2. The Prandtl number in the fluid is set to Pr= 0.75, hence an exact solution of the compressible Navier-Stokes equations exists, see e.g. [9,14] and [40] for a detailed description of the test problem and for details concerning the computation of the exact solution. The problem is solved in the domain = [−0.5, +0.5] and the shock wave

Solid Rotor Problem
As fourth test, we study the solid rotor problem proposed in [15,22], by solving the GPR The reference solution is provided by a second order MUSCL-Hancock scheme on 512 × 512 control volumes, see [15,22]. In Fig. 7 the numerical results obtained with the HTC DG scheme are compared with the reference solution. We observe a very good agreement, also with the results published previously in [15,22].

Lid-Driven Cavity
We now apply our new thermodynamically compatible HTC DG scheme to the well-known lid-driven cavity problem, see [51]. This test can be used to validate compressible flow solvers in the low Mach number regime, see e.g. [11,24,94] and was already successfully solved with the GPR model, see [15,22,40]. The two-dimensional computational domain is = [0, 1] 2 and the initial condition is given by ρ = 1, v = 0, p = 10 2 , A = I and J = 0. The lid velocity on the upper boundary is v = (1, 0, 0), while solid no-slip walls, v = 0, are imposed everywhere else. As model parameters we use c s = 8, c h = 2, τ 2 = 10 −2 and μ = 10 −2 , i.e. the associated Reynolds number is Re = 100, while the characteristic Mach number of the flow based on the lid velocity is about M = 0.08. The new HTC DG scheme is run until t = 10 using 256 × 256 elements with polynomial approximation degree N = 3 and = 10 −3 . The obtained computational results are depicted in Fig. 8, where also a comparison with the reference solution of Ghia et al. [51] is shown and which is based on the solution of the incompressible Navier-Stokes equations. One can observe an excellent agreement between the numerical results obtained with the new HTC DG scheme applied to the GPR model and the reference solution. We stress that the reference solution has been obtained with a different numerical scheme that was applied to a different PDE system.
We would like to clarify that the lid-driven cavity generates pressure peaks in the upper corners of the domain due to the discontinuous velocity field in the boundary conditions, which requires limiting of high order DG schemes. This is the reason why the simulation was carried out with artificial viscosity, see also [22] for the same artificial viscosity parameter used in thermodynamically compatible finite volume schemes.

Double Shear Layer
In this last test case concerning the GPR model, we solve the double shear layer problem of [10]. This benchmark was also used in [11,15,22,24,40,93,94] to assess the behaviour of compressible flow solvers in the weakly compressible regime, including applications to the GPR model. The computational domain is = [0, 1] 2 , the boundary conditions are periodic everywhere and the test is run until a final time of t = 1.8. The initial condition reads We run our new HTC DG scheme on a mesh of 1024 × 1024 elements, using a polynomial approximation degree of N = 3 and = 1 × 10 −6 . In Fig. 9 the temporal evolution of the distortion field component A 12 is shown. One can observe how the initial shear layers develop into several vortex structures. A more detailed analysis of the flow has been provided in [10,11,15,24,40]. Our computational results agree very well with those obtained with the thermodynamically compatible finite volume scheme [22] run on 4000 × 4000 control volumes.

Numerical Convergence Study
In order to verify the convergence of the new HTC DG scheme for the turbulent shallow water model experimentally we consider a manufactured solution test proposed in [21]. The test is run on the computational domain = [0, 2π] 2 and with periodic boundary conditions everywhere. The solution of system (11) is defined as follows: Q(x, t) = q 0 sin(x) cos(y) cos(t) − sin(x) cos(y) cos(t) − cos(x) sin(y) cos(t) cos(x) sin(y) cos(t) (73) with the total energy This choice of Q(x, t) yields a Reynolds stress tensor of the form with c 1 = sin(x) cos(y) and c 2 = sin(y) cos(x). To complete the definition of the problem we set h 0 = 1 and q 0 = 0.5. Let us remark that to get the sought solution, (73)-(74), a set of source terms must be added to the right hand side of (11). The expressions of the source terms can be simply calculated by substitution of (73)- (74) in (11). The simulation is run until a final time of t = 0.1 using the new HTC DG schemes proposed in this paper and using polynomial degrees N ∈ {2, 3, 4, 5}. The artificial viscosity coefficient is set to = 0 inside each element, while for this test case the penalty parameter inside the Rusanov-type flux G is set to η = s max . The errors in L 2 norm obtained for h, hv 1 and Q 11 are reported in Table 3. From the obtained results we can conclude that, overall, the obtained experimental order of accuracy of the scheme is N + 1 for the water depth h, while for the other quantities we observe N + 1 for odd and N for even polynomial approximation degrees N .

Riemann Problems
In this section, we address a set of three Riemann problems proposed and solved in [21] for the turbulent shallow water model (11). An exact Riemann solver for the original model [50,62] based on the evolution of the Reynolds stress tensor P instead of its decomposition Q was recently forwarded in [74]. For alternative numerical schemes applied to the hyperbolic shear shallow water model the reader is referred to [12,21,29,74]. An important difficulty present in the model (11) and also in the original model proposed in [50] is its non-conservative formulation, with non-conservative products active across genuinely nonlinear fields. This makes a proper numerical discretization in the presence of shock waves extremely challenging. In [21] it was shown that thermodynamic compatibility, in particular compatibility with the total energy conservation law, is a necessary key ingredient for a correct discretization of the equations. Thermodynamic compatibility can be either achieved via the original wave splitting proposed in [50], or by the unsplit schemes introduced in [21], which either enforce energy conservation by explicitly solving the energy equation and an appropriate rescaling The simulations were run on Cartesian meshes of N x × N x elements up to time t = 0.1 of the object Q, or via a thermodynamically compatible discretization that is consistent with (11) at a discrete level, including the entropy production terms ik , or via a direct numerical simulation (DNS) of the viscous equations in the vanishing viscosity limit. In the following we solve three Riemann problems in the one-dimensional domain = [0, 1] with initial data given in Table 4, where also the final times are reported for each test case. The state variables which are not explicitly indicated in Table 4 are set to zero, i.e. v 1 = 0 and Q 21 = 0. All simulations are run with the thermodynamically compatible DG scheme proposed in this paper, using 4096 elements and a polynomial approximation degree of N = 9. The artificial viscosity is set to = 10 −6 in RP1 and RP2 and to = 2 × 10 −6 for RP3. The obtained computational results are shown in Figs. 10, 11 and 12, where also a comparison with the scheme of Ivanova et al. [62] is provided, together with a comparison The numerical results obtained with the HTC DG scheme are in excellent agreement with both reference solutions, even for RP3 with its complex wave pattern. This highlights the importance of the discrete thermodynamic compatibility when discretizing non-conservative hyperbolic equations like system (11).

Roll Waves
The last test case concerns the simulation of one-dimensional roll waves and was originally proposed and numerically solved in [50]. The obtained numerical results can be compared with the experimental data provided by Brock in [16,17]. The one-dimensional domain = [0, L] is periodic with length L = 1.3. According to [50] the initial condition reads h = h 0 (1 + a sin(2π x/L)) with a = 0.05, v 1 = gh 0 tan θ/C f , v 2 = 0 and Q = 1 2 ϕh 2 I. The remaining parameters of this test problem are θ = 0.05011, ∂ x b = tan(θ ), h 0 = 0.00798, C f = 0.0036, C r = 0.00035 and ϕ = 22.76, see [50]. In order to simulate this test case, the following dissipative source terms must be added to the right hand side of (11): in the momentum equation (11b) we add the term −C f √ v m v m v i to the right hand side and in the PDE for Q ik , Eqn. (11c), we add −α/h Q ik , with α = max 0, C r trP−ϕh 2 trP 2 v 3 to the right hand side, see [50] and [21] for details. It is obvious that in the presence of the dissipative source terms, total energy in the turbulent flow is no longer conserved.
We run the new thermodynamically compatible DG scheme proposed in this paper using 512 elements, a polynomial approximation degree of N = 5 and an artificial viscosity of = 4 × 10 −5 . Simulations are run until a final time of t = 12.5. As already suggested in [21], the bottom slope is simply implemented as an algebraic source. In Fig. 13 the numerical results obtained with the new HTC DG scheme are compared with the experimental profile of Brock documented in [16,17,21,50]. In order to perform this comparison, the spatial coordinate has been normalized by the reference length L and the wather depth has been normalized by h 0 , i.e. we plot h/h 0 over x/L. The numerical results have been shifted in x direction so that the numerical shock location coincides with the experimental one. In total two periods of the simulation are shown. Overall, one can note a good agreement between the numerical results and the experimental reference data.

Conclusion
In this paper we have introduced new high order accurate and thermodynamically compatible discontinuous Galerkin finite element schemes for the Godunov-Peshkov-Romenski (GPR) model of continuum mechanics and for the first order hyperbolic model of turbulent shallow water flows of Gavrilyuk et al. [50,62]. The key feature of our new methods is that they discretize the entropy inequality directly as a primary evolution equation, while total energy conservation is achieved as a mere consequence of the thermodynamically compatible discretization of all other equations. The new DG methods satisfy a cell entropy inequality by construction, and can be proven to be nonlinearly stable in the energy norm due to the thermodynamic compatibility. The new DG schemes proposed in this paper can therefore be seen as the missing dual algorithms to known entropy consistent DG schemes, which usually discretize the total energy conservation law directly and obtain the entropy inequality as a consequence, see e.g. [30,38,48,60,70,87,89].
To the very best of our knowledge, the DG schemes presented in this paper are the first provably thermodynamically compatible DG schemes for the unified model of continuum mechanics of Godunov, Peshkov and Romenski [54,77] and for the turbulent shallow water model of Gavrilyuk et al. [50,62]. It is also the first time that a DG scheme has been developed which does not directly discretize the total energy conservation law, but which obtains total energy conservation as a consequence of a compatible discretization of all the other PDE. Compared to existing HTC finite volume schemes [20][21][22] the mathematical proofs for the DG framework are more complex, in particular concerning the design of a proper thermodynamically compatible numerical viscosity, for which the jump terms V must have opposite sign compared to classical symmetric interior penalty Galerkin schemes [47].
We have applied our new DG schemes to a wide range of test problems in one and two space dimensions, obtaining overall a very good agreement with available exact, numerical and even experimental reference solutions. As in previous work on thermodynamically compatible finite volume schemes [20][21][22] the high order DG schemes presented in this paper rely on a path-integral and a Godunov parametrization of the physical flux f k = ∂ p (v k L) in terms of a generating potential L in order to construct a thermodynamically compatible flux for the inviscid Euler and shallow water subsystems.
We clearly emphasize again that for all proofs provided in this paper we assume all integrals to be exact. Further work on provably thermodynamically compatible DG schemes including numerical quadrature errors will be necessary in the future, also considering the recent ideas of Abgrall et al. [1,3].
In the future, we plan to extend our thermodynamically compatible DG schemes also to the magnetohydrodynamics (MHD) equations and to the conservative SHTC system of compressible two-fluid flows proposed and studied in [71,84,86,95]. The thermodynamically compatible DG schemes presented in this paper have been analyzed only in the semi discrete setting. The fully discrete case is clearly out of scope of the present paper. Further research will also consider fully-discrete thermodynamically compatible time discretizations, either making use of the line integral approach of Brugnano and Iavernaro [18,19], or by extending the exactly conservative fully-discrete FV scheme forwarded in [22] for the Euler subsystem to the full GPR model and other more complex SHTC systems. However, the approach presented in [22] requires the calculation of q = q(p) which is very complex in case of the GPR model. For alternative extensions to the fully discrete case, at least for the Euler subsystem, see e.g. [22,73,79,80].
Moreover, we will investigate whether it is possible to extend thermodynamically compatible DG schemes in such a way as to preserve curl and divergence involution constraints exactly, see e.g. [15,39] for involution-preserving semi-implicit discretizations on staggered meshes. We will also consider an extension of the framework of thermodynamically compatible schemes to the class of semi-implicit hybrid finite-volume / finite-element methods [11,23,24] on staggered unstructured meshes.
The development of provably thermodynamically compatible limiters for the DG scheme presented in this paper was out of scope of this work and will be left to future studies. Instead of sophisticated limiters, a simple artificial viscosity approach has been used in this paper. In the future we will consider also slope and moment limiters [64,65], provably positivity preserving limiters [97], as well as the use of cell-centered thermodynamically compatible finite volume schemes [3,[20][21][22] as a posteriori subcell limiters, similar to the ideas on subcell limiting for DG schemes presented in [41,60,90,91].
A fundamental mathematical property of any reasonable numerical scheme for nonlinear systems of hyperbolic conservation laws that was not investigated in this paper at all was the invariant domain preserving property (IDP), such as the positivity of density and entropy. Future work will consider provably IDP extensions of our methods, making use of the mathematical techniques presented in the seminal work of Guermond and Popov et al. concerning provably invariant domain preserving schemes [31,[57][58][59] and in the work of Kuzmin et al. [7,[65][66][67] concerning bound-preserving algebraic flux limiters and slope limiters for high order continuous and discontinuous Galerkin finite element schemes.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.