# A discontinuous Galerkin spectral element method for a nonconservative
compressible multicomponent flow model

Rémi Abgrall  Pratik Rai<EMAIL_ADDRESS>Florent Renac  Institute of
Mathematics, University of Zürich, Switzerland  CMAP, École Polytechnique,
Route de Saclay, 91128 Palaiseau Cedex, France  DAAA, ONERA, Université Paris
Saclay F-92322 Châtillon, France

###### Abstract

In this work, we propose an accurate, robust (the solution remains in the set
of states), and stable discretization of a nonconservative model for the
simulation of compressible multicomponent flows with shocks and material
interfaces. We consider the gamma-based model by Shyue [J. Comput. Phys., 142
(1998), 208–242] where each component follows a stiffened gas equation of
state (EOS). We here extend the framework proposed in Renac [J. Comput. Phys.,
382 (2019), 1–26] and Coquel et al. [J. Comput. Phys. 431 (2021), 110135] for
the discretization of hyperbolic systems, with both fluxes and nonconservative
products, to unstructured meshes with curved elements in multiple space
dimensions. The framework relies on a high-order discontinuous Galerkin
spectral element method (DGSEM) using collocation of quadrature and
interpolation points as proposed by Gassner [SIAM J. Sci. Comput., 35 (2013)]
in the case of hyperbolic conservation laws. We modify the integrals over
discretization elements where we replace the physical fluxes and
nonconservative products by two-point numerical fluctuations. The
contributions of this work are threefold. First, we analyze the semi-discrete
DGSEM discretization of general hyperbolic systems with conservative and
nonconservative terms and derive the conditions to obtain a scheme that is
high-order accurate, free-stream preserving, and entropy stable when excluding
material interfaces. Second, we design a three-point scheme with a HLLC solver
for the gamma-based model that does not require a root-finding algorithm for
the approximation of the nonconservative products. The scheme is proved to be
robust and entropy stable for convex entropies, to preserve uniform profiles
of pressure and velocity across material interfaces (material interface
preservation), and to satisfy a discrete minimum principle on the specific
entropy and maximum principles on the parameters of the EOS. Third, the HLLC
solver is applied at interfaces in the DGSEM, while we consider two kinds of
fluctuations in the integrals over discretization elements: the former is
entropy conservative (EC), while the latter preserves material interfaces
(CP). Time integration is performed using high-order strong-stability
preserving Runge-Kutta schemes. The fully discrete scheme is shown to preserve
material interfaces with CP fluctuations. Under a given condition on the time
step, both EC and CP fluctuations ensure that the cell-averaged solution
remains in the set of states; satisfy a minimum principle on any convex
entropy and maximum principles on the EOS parameters. These results have
allowed us to use existing limiters in order to restore positivity, and
discrete maximum principles of degrees-of-freedom within elements. Numerical
experiments in one and two space dimensions on flows with discontinuous
solutions support the conclusions of our analysis and highlight stability,
robustness and accuracy of the proposed DGSEM with either CP, or EC
fluctuations, while the scheme with CP fluctuations is shown to offer better
resolution capabilities.

Keywords.  Compressible multicomponent flows, Nonconservative hyperbolic
systems, Discontinuous Galerkin method, Summation-by-parts, Material interface
capturing, Entropy stability, High-order accuracy

AMS subject classifications.  65M12, 65M70, 76T10

##  1 Introduction

The discussion in this paper focuses on the approximation in multiple space
dimensions of a compressible multicomponent flow model in nonconservative
form. We consider a gamma-based model [58] for a mixture with a stiffened gas
equation of state (EOS) approximating components including both gas and
compressible liquids (hereafter referred to as the SG-gamma model). The model
is written in quasi-conservative form [1] to preserve velocity and pressure
profiles across material interfaces separating different components. The model
approximates mixture quantities and presents the main advantage of being
independent of the number of components. We are here interested in high-order,
robust (i.e., preserving the solution in the set of admissible states), and
entropy stable simulations of flows with shocks, material interfaces, and
complex interactions triggering small scale flow phenomena. Numerical
approximation of multicomponent and multiphase flows based on interface
capturing methods has been the subject of numerous works (see, e.g. [39, 56]
and references therein).

We discretize the SG-gamma model using the discontinuous Galerkin spectral
element method (DGSEM) with Gauss-Lobatto quadrature rules [42]. Such schemes
were first introduced in previous works [5, 20], and in [27] the DGSEM was
reinterpreted as a diagonal norm summation-by-parts (SBP) operator that falls
into the general framework of conservative elementwise flux differencing
schemes [23]. In this framework, using entropy conservative (EC) numerical
fluxes from Tadmor [60], semi-discrete EC finite-difference and spectral
collocation schemes have been derived in [23, 27]. The framework has been
extended to nonconservative hyperbolic systems on Cartesian meshes in [53, 15]
by using EC numerical fluctuations from [10]. In both frameworks, a semi-
discrete entropy inequality may be obtained by replacing physical fluxes and
nonconservative products with EC numerical fluxes or fluctuations within
discretization elements, while using entropy stable ones at interfaces between
elements. The design of the latter relies either on adding upwind-type
dissipation [36] to EC numerical fluxes and fluctuations [8, 24, 18, 28, 53,
15], or on designing approximate Riemann solvers [19, 54, 52, 51]. Note that
the SBP operators take into account the numerical quadrature when
approximating integrals compared to other schemes that require their exact
evaluation to achieve entropy stability [37, 33, 34].

Here, we extend the framework proposed in [53] to multidimensional
unstructured meshes with curved elements by using tensor multiplication of
quadrature rules and function basis [42] that satisfy geometric conservation
laws (the so-called metric identities [41]) at the discrete level. This
framework has been recently applied to the approximation with a well-balanced
DGSEM of balance laws with geometric source terms on multidimensional high-
order meshes in [68]. We here rather focus on specific properties of
discretizations of nonconservative multicomponent flows: preservation of
material interfaces, discrete conservation of physical fluxes, maximum
principles on purely transported quantities such as the EOS parameters, and
fully discrete entropy stability. The latter property presents some
difficulties due to the form of the entropy associated to SG-gamma model.
First, as a model for the mixture, the properties of the individual
components, such as mass and void fractions, are not known in general which
prevents the evaluation of the entropy variables necessary for the derivation
of EC fluctuations from the Castro et al. condition [10]. Then, the entropy is
not convex as is often the case in phase transition models [32], which will
restrict the entropy stability across shocks not interacting with material
fronts. We here circumvent these difficulties by considering a specific
entropy that allows the evaluation of the entropy variables. This entropy
satisfies a Gibbs relation and thus defines a complete EOS [46] and is concave
with respect to two thermodynamic intensive properties of the mixture, so
entropy stability can be ensured when excluding material interfaces. Let us
stress that material interfaces do not require the scheme to be entropy
stable, but rather a consistent approximation of the energy equation to ensure
pressure equilibrium [1]. We thus also design material interface preserving
(CP) fluctuations to preserve at the discrete level pressure and velocity
fields across such interfaces.

We then design a HLLC solver [63, 65] for the SG-gamma model that does not
require a root-finding algorithm to evaluate the nonconservative product in
contrast to other schemes [21, 12, 61]. We analyze the properties of a three-
point scheme using the HLLC solver and prove that the scheme is robust and
entropy stable for convex entropies defining a complete EOS, preserves uniform
profiles of pressure and velocity across material interfaces, and satisfies a
discrete minimum principle on the specific entropy and maximum principles on
the parameters of the EOS. We then apply the HLLC solver at mesh interfaces in
the numerical scheme and analyze the properties of the fully discrete scheme
with an explicit first-order Euler time integration. We derive conditions on
the time step so that the cell-averaged solution is a convex combination of
degrees-of-freedom (DOFs) and updates of three-point schemes, so the scheme
inherits the properties of the three-point scheme. In particular, the proposed
DGSEM satisfies a minimum principle on the entropy irrespective of the
fluctuations (EC or CP) that are used in the discretization elements. As a
consequence, the DGSEM with CP fluctuations within elements and the HLLC
solver at interfaces is able to handle shocks and to preserve material
interfaces. Time integration is performed using high-order strong-stability
preserving Runge-Kutta schemes that are convex combinations of explicit Euler
schemes, while linear scaling limiters [70, 71] are applied at the end of each
stage to impose positivity and maximum principles at all DOFs within
discretization elements.

This paper is organized as follows. The SG-gamma model and the entropy pair
are described in section 2. We then introduce the semi-discrete DGSEM on
multidimensional and high-order unstructured meshes in section 3. In section
4, we derive CP and EC numerical fluxes for the SG-gamma model, and in section
5, we propose the HLLC approximate Riemann solver and analyze its properties.
We recall the main properties of the fully-discrete DGSEM in section 6. A
posteriori limiters are described in section 7. The results are assessed by
numerical experiments in one and two space dimensions in section 8 and
concluding remarks about this work are given in section 9.

##  2 The SG-gamma model

###  2.1 Governing equations and thermodynamic model

In this work, we consider the gamma-based compressible multicomponent flow
model where each component is assumed to be a stiffened gas [58] and refer to
it as the SG-gamma model. The main interest in this model is that the number
of unknowns is independent of the number of components. We are here interested
in high-order approximations of the associated Cauchy problem in $d$ space
dimensions for flows with $n_{s}$ components:

$\displaystyle\partial_{t}{\textbf{u}}+\nabla\cdot\textbf{f}({\textbf{u}})+\textbf{c}({\textbf{u}})\nabla{\textbf{u}}$
$\displaystyle=0,\quad$ $\displaystyle{\bf x}\in\mathbb{R}^{d},t>0,$ (1a)
$\displaystyle{\textbf{u}}({\textbf{x}},0)$
$\displaystyle={\textbf{u}}_{0}({\textbf{x}}),\quad$ $\displaystyle{\bf
x}\in\mathbb{R}^{d},$ (1b)

where ${\textbf{x}}=(x_{1},\dots,x_{d})$ are the spatial coordinates,
$\textbf{f}({\textbf{u}})=\big{(}\textbf{f}_{1}({\textbf{u}}),\dots,\textbf{f}_{d}({\textbf{u}})\big{)}$,
$\textbf{c}({\textbf{u}})\nabla{\textbf{u}}=\sum_{i=1}^{d}\textbf{c}_{i}({\textbf{u}})\partial_{x_{i}}{\textbf{u}}$,
and

${\textbf{u}}=\begin{pmatrix}\rho\\\ \rho{\textbf{v}}\\\ \rho E\\\ \Gamma\\\
\Pi\end{pmatrix},\quad\textbf{f}({\textbf{u}})=\begin{pmatrix}\rho{\textbf{v}}^{\top}\\\
\rho{\textbf{v}}{\textbf{v}}^{\top}+\mathrm{p}{\textbf{I}}\\\ (\rho
E+\mathrm{p}){\textbf{v}}^{\top}\\\ 0\\\
0\end{pmatrix},\quad\textbf{c}({\textbf{u}})\nabla{\textbf{u}}=\begin{pmatrix}0\\\
0\\\ 0\\\ {\textbf{v}}\cdot\nabla\Gamma\\\
{\textbf{v}}\cdot\nabla\Pi\end{pmatrix},$ (2)

represent the vector of state variables, the physical fluxes and the
nonconservative products111The components of the third-order tensor in (eq.
1a) read ${\bf c}({\bf u})_{ijk}={\bf c}_{k}({\bf u})_{ij}$ with the ${\bf
c}_{k}({\bf u})$ in $\mathbb{R}^{n_{eq}\times n_{eq}}$, so the $i$th component
of ${\bf c}({\bf u})\nabla{\bf u}$ reads
$\sum_{j=1}^{n_{eq}}\sum_{k=1}^{d}{\bf c}({\bf u})_{ijk}\partial_{x_{k}}{\bf
u}_{j}$. Likewise, for ${\textbf{n}}=(n_{1},\dots,n_{d})^{\top}$, we have
$\textbf{c}({\textbf{u}}){\textbf{n}}=\sum_{i=1}^{d}n_{i}\textbf{c}_{i}({\textbf{u}})$.,
respectively. The mixture density, momentum, total energy, and internal energy
are defined as

$\rho=\sum_{i=1}^{n_{s}}\alpha_{i}\rho_{i},\quad\rho{\textbf{v}}=\sum_{i=1}^{n_{s}}\alpha_{i}\rho_{i}{\textbf{v}}_{i},\quad\rho
E=\rho e+\frac{1}{2}\rho{\textbf{v}}\cdot{\textbf{v}},\quad\rho
e=\sum_{i=1}^{n_{s}}\alpha_{i}\rho_{i}e_{i},$

where $\rho_{i}$, ${\textbf{v}}_{i}$, and $e_{i}$ represent the density, the
velocity vector and the specific internal energy of the $i$th component. The
model assumes immiscible phases and a saturation condition on the void
fractions $\alpha_{i}$:

$\sum_{i=1}^{n_{s}}\alpha_{i}=1.$ (3)

The partial pressures are related to partial densities and specific internal
energies through the stiffened gas EOS:

$\mathrm{p}_{i}(\rho_{i},e_{i})=(\gamma_{i}-1)\rho_{i}e_{i}-\gamma_{i}\mathrm{p}_{\infty_{i}},\quad
e_{i}=\mathrm{C_{v}}_{i}T_{i}+\frac{\mathrm{p}_{\infty_{i}}}{\rho_{i}},\quad
i=1,\dots,n_{s},$ (4)

where $\gamma_{i}=\mathrm{C_{p}}_{i}/\mathrm{C_{v}}_{i}>1$ is the ratio of
specific heats, $T_{i}$ is the temperature of the species, and
$\mathrm{p}_{\infty_{i}}\geqslant 0$ is a pressure-like constant. Observe that
when $\mathrm{p}_{\infty_{i}}=0$ in (eq. 4) we recover the polytropic EOS.
Assuming thermal equilibrium of the species,
$T_{i}(\rho_{i},e_{i})=T(\rho,e)$, $1\leqslant i\leqslant n_{s}$, the EOS for
the mixture is conveniently defined by [58]

$\frac{\mathrm{p}+\gamma\mathrm{p}_{\infty}}{(\gamma-1)}=\mathrm{p}\Gamma+\Pi=\rho
e,\quad\rho
e=\rho\mathrm{C_{v}}T+\sum_{i=1}^{n_{s}}\alpha_{i}\mathrm{p}_{\infty_{i}},$
(5)

where $\mathrm{C_{v}}=\sum_{i=1}^{n_{s}}Y_{i}\mathrm{C_{v}}_{i}$ denotes the
specific heat at constant volume of the mixture, the
$Y_{i}=\tfrac{\alpha_{i}\rho_{i}}{\rho}$ are the mass fractions of the
species, and the EOS parameters $\Gamma$ and $\Pi$ are defined by

$\Gamma=\frac{1}{\gamma-1}=\sum_{i=1}^{n_{s}}\frac{\alpha_{i}}{\gamma_{i}-1},\quad\Pi=\frac{\gamma\mathrm{p}_{\infty}}{\gamma-1}=\sum_{i=1}^{n_{s}}\frac{\alpha_{i}\gamma_{i}\mathrm{p}_{\infty_{i}}}{\gamma_{i}-1}.$
(6)

Hyperbolicity of the SG-gamma model requires that the solutions to (eq. 1)
belong to the set of states

$\Omega_{\mathrm{GM}}=\left\\{{\textbf{u}}\in\mathbb{R}^{n_{eq}}\mathrel{\mathop{\mathchar
58\relax}}\;\rho>0,\;{\textbf{v}}\in\mathbb{R}^{d},\;\rho
e>\mathrm{p}_{\infty},\;\Gamma>0,\;\Pi\geqslant 0\right\\},$ (7)

with $n_{eq}=d+4$. The matrix-valued function
$\sum_{i=1}^{d}n_{i}\big{(}\textbf{f}_{i}^{\prime}({\textbf{u}})+\textbf{c}_{i}({\textbf{u}})\big{)}$
in $\mathbb{R}^{n_{eq}\times n_{eq}}$ admits real eigenvalues

$\lambda_{1}({\textbf{u}})={\textbf{v}}\cdot{\textbf{n}}-c,\quad\lambda_{2}({\textbf{u}})=\dots=\lambda_{n_{eq}-1}({\textbf{u}})={\textbf{v}}\cdot{\textbf{n}},\quad\lambda_{n_{eq}}({\textbf{u}})={\textbf{v}}\cdot{\textbf{n}}+c,$
(8)

for all unit vector ${\textbf{n}}=(n_{1},\dots,n_{d})$, where
$c=\sqrt{\gamma(\gamma-1)(\rho e-\mathrm{p}_{\infty})/\rho}$ is the speed of
sound of the mixture. Here, $\\{\lambda_{i}\\}_{2\leqslant i\leqslant
n_{eq}-1}$ are associated to linearly degenerate (LD) fields, while
$\lambda_{1}$ and $\lambda_{n_{eq}}$ are associated to genuinely nonlinear
(GNL) fields. Observe that (eq. 1a) is not strictly hyperbolic as the
eigenvalues associated to the LD fields are not distinct.

###  2.2 Entropy pair

Solutions to (eq. 1) may develop discontinuities and (eq. 1a) has to be
understood in a weak sense. Weak solutions are not necessarily unique and (eq.
1) must be supplemented with further admissibility conditions to select the
physical solution. We here focus on entropy inequalities

$\partial_{t}\eta({\textbf{u}})+\nabla\cdot\textbf{q}({\textbf{u}})\leqslant
0,\qquad{\bf x}\in\mathbb{R}^{d},t>0,$ (9)

for some convex entropy – entropy flux pair
$(\eta({\textbf{u}}),\textbf{q}({\textbf{u}}))$. One common way to derive such
pair is by defining partial entropies of the species
$s_{i}(\rho_{i},\theta)=-\mathrm{C_{v}}_{i}\left(\ln\theta+(\gamma_{i}-1)\ln\rho_{i}\right)$,
where $\theta=1/T$ is the inverse of the temperature, that satisfies a Gibbs
relation

$T\mathrm{d}s_{i}=\mathrm{d}e_{i}-\frac{\mathrm{p}_{i}}{\rho^{2}_{i}}\mathrm{d}\rho_{i},$
(10)

so the mixture entropy reads

$\sum_{i=1}^{n_{s}}Y_{i}s_{i}=-\mathrm{C_{v}}\ln\theta-\sum_{i=1}^{n_{s}}Y_{i}(\gamma_{i}-1)\mathrm{C_{v}}_{i}\ln\rho_{i}.$

Unfortunately, we cannot use this entropy to derive EC fluxes in section 4 as
$Y_{i}$ is an unknown, that cannot be evaluated from u in (eq. 2). Hence, here
we consider an alternative pair

$\eta({\textbf{u}})=-\rho s,\quad\textbf{q}({\textbf{u}})=-\rho
s{\textbf{v}},$ (11)

where, upon introducing $\tau=\tfrac{1}{\rho}$, the covolume of the mixture,
the specific entropy reads

$s(\tau,e,\mathrm{C_{v}},\Gamma,\Pi)=\mathrm{C_{v}}\ln\left(\frac{\mathrm{p}+\mathrm{p}_{\infty}}{\rho^{\gamma}}\right)\overset{(\lx@cref{creftype~refnum}{eqn:
chap 5 Gamma and
Pi})}{=}\mathrm{C_{v}}\Big{(}\ln\Big{(}e-\frac{\Pi}{\Gamma+1}\tau\Big{)}+\frac{1}{\Gamma}\ln\tau-\ln\Gamma\Big{)}.$
(12)

The rationale for considering this entropy is as follows. First, we will see
in section 4.2 that the EC fluctuations can be explicitly computed without
knowledge of the individual mass fractions. Then, for smooth solutions of (eq.
1a), we have

$\partial_{t}\mathrm{C_{v}}+{\textbf{v}}\cdot\nabla\mathrm{C_{v}}=0,\quad\partial_{t}\tau+{\textbf{v}}\cdot\nabla\tau=\tau\nabla\cdot{\textbf{v}},\quad\partial_{t}e+{\textbf{v}}\cdot\nabla
e=-\tau\mathrm{p}\nabla\cdot{\textbf{v}},$

hence

$\partial_{t}s+{\textbf{v}}\cdot\nabla
s=\partial_{\tau}s\big{(}\partial_{t}\tau+{\textbf{v}}\cdot\nabla\tau\big{)}+\partial_{e}s\big{(}\partial_{t}e+{\textbf{v}}\cdot\nabla
e\big{)}=\tau\nabla\cdot{\textbf{v}}(\partial_{\tau}s-\mathrm{p}\partial_{e}s)=\mathrm{C_{v}}\nabla\cdot{\textbf{v}}\left(\frac{1}{\Gamma}-\frac{\mathrm{p}_{\infty}}{\rho
e-\mathrm{p}_{\infty}}-\frac{\mathrm{p}}{\rho
e-\mathrm{p}_{\infty}}\right)\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn:
chap 5 mixture EOS})}}}{{=}}0,$

so the mixture entropy is conserved, i.e., (eq. 9) is an equality. Moreover,
$(\tau,e)\mapsto s(\tau,e,\mathrm{C_{v}},\Gamma,\Pi)$ is obviously strictly
concave in $\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm GM$}}$ and from the
first equality in the above relation, we conclude that the inequality in (eq.
9) makes sense even if $\eta({\bf u})$ is not convex, which is also the case
of the physical entropy, $\sum_{i=1}^{n_{s}}Y_{i}s_{i}$, or in some phase
transition models [32]. Finally, differentiating (eq. 12) while fixing
$\mathrm{C_{v}}$, $\Gamma$ and $\Pi$, gives

$\mathrm{d}s=\mathrm{C_{v}}\left(\frac{\mathrm{d}\mathrm{p}}{\mathrm{p}+\mathrm{p}_{\infty}}-\gamma\frac{\mathrm{d}\rho}{\rho}\right)\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn:
chap 5 mixture
EOS})}}}{{=}}\mathrm{C_{v}}\left(\frac{(\gamma-1)\mathrm{d}(\rho
e)}{\mathrm{p}+\mathrm{p}_{\infty}}-\gamma\frac{\mathrm{d}\rho}{\rho}\right)=\frac{\rho\mathrm{C_{v}}}{\rho
e-\mathrm{p}_{\infty}}\left(\mathrm{d}e-\frac{\mathrm{p}}{\rho^{2}}\mathrm{d}\rho\right),$
(13)

so the entropy satisfies a Gibbs relation similar to (eq. 10) but with a
different temperature $\tilde{T}=\tfrac{1}{\rho\mathrm{C_{v}}}(\rho
e-\mathrm{p}_{\infty})$ instead of $T$ in (eq. 5). The entropy hence defines a
complete EOS with both pressure and temperature [46]. This latter observation
will be important in section 5 to prove entropy stability of the HLLC solver
through the existence of local minimum entropy principles (see proof of Lemma
5.1).

We end this section by deriving the entropy variables associated to the
entropy $\eta({\textbf{u}})$ in (eq. 11) when considering pure phases. For
pure phases, we have
$\mathrm{d}\Gamma=\mathrm{d}\Pi=\mathrm{d}\mathrm{C_{v}}=0$ so we obtain

$\displaystyle\mathrm{d}\eta=-\rho\mathrm{d}s-s\mathrm{d}\rho\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn:
chap 5 mixture physical
entropy})}}}{{=}}-\rho\mathrm{C_{v}}\frac{\mathrm{d}(\mathrm{p}+\mathrm{p}_{\infty})}{\mathrm{p}+\mathrm{p}_{\infty}}+\gamma\mathrm{C_{v}}\mathrm{d}\rho-s\mathrm{d}\rho$
$\displaystyle\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn:
chap 5 mixture EOS})}}}{{=}}-\rho\mathrm{C_{v}}(\gamma-1)\frac{\mathrm{d}\rho
e}{\mathrm{p}+\mathrm{p}_{\infty}}+\gamma\mathrm{C_{v}}\mathrm{d}\rho-s\mathrm{d}\rho$
$\displaystyle=-\rho\frac{(\gamma-1)\mathrm{C_{v}}}{\mathrm{p}+\mathrm{p}_{\infty}}\left(\mathrm{d}\rho
E-{\textbf{v}}\cdot\mathrm{d}\rho{\textbf{v}}+\frac{{\textbf{v}}\cdot{\textbf{v}}}{2}\mathrm{d}\rho\right)+\gamma\mathrm{C_{v}}\mathrm{d}\rho-s\mathrm{d}\rho$
$\displaystyle\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eq:def_zeta})}}}{{=}}-\zeta\mathrm{d}\rho
E+\zeta{\textbf{v}}\cdot\mathrm{d}\rho{\textbf{v}}+\left(\gamma\mathrm{C_{v}}-s-\zeta\frac{{\textbf{v}}\cdot{\textbf{v}}}{2}\right)\mathrm{d}\rho,$

where

$\zeta=\frac{(\gamma-1)\mathrm{C_{v}}\rho}{\mathrm{p}+\mathrm{p}_{\infty}},$
(14)

and we get the following expression of the entropy variables

$\bm{\vartheta}({\textbf{u}})=\frac{\partial}{\partial{\textbf{u}}}\eta({\textbf{u}})=\begin{pmatrix}\gamma\mathrm{C_{v}}-s-\zeta\frac{{\textbf{v}}\cdot{\textbf{v}}}{2}\\\
\zeta{\textbf{v}}\\\ -\zeta\\\ 0\\\ 0\end{pmatrix}.$ (15)

##  3 The discontinuous Galerkin spectral element method (DGSEM)

In this section, we recall the DGSEM framework [15, 42, 53] which is used to
discretize the Cauchy problem (eq. 1). Here, the space domain
$\Omega=\mathbb{R}^{d}$ is discretized using a mesh $\Omega_{h}$ consisting of
nonoverlapping and nonempty cells $\kappa$ (quadrangles for $d=2$ and
hexahedra for $d=3$) forming a partition of $\Omega$. By ${\cal E}_{h}$ we
denote the set of interfaces in $\Omega_{h}$. For the sake of clarity, we
introduce the DGSEM in two space dimensions $d=2$, as the extension to $d=3$
is straightforward while its derivation for $d=1$ can be found in [15, 53].

$\kappa=\kappa^{-}$$\kappa_{e}^{+}$$e$${\bf u}_{h}^{-}$${\bf u}_{h}^{+}$${\bf
n}_{e}^{k}$${\bf x}_{e}^{k}$${\bf
x}_{\kappa}^{ij}$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$
Figure 1: Notations for the mesh in two space dimensions ($d=2$): cell
$\kappa\in\Omega_{h}$ with quadrature points ${\textbf{x}}_{\kappa}^{ij}$
(bullets $\bullet$), edge $e\in\partial\kappa$ with quadrature point
${\textbf{x}}_{e}^{k}$, associated unit outward normal ${\textbf{n}}_{e}^{k}$,
traces of the approximate solution ${\textbf{u}}^{\pm}_{h}$ on $e$; and
adjacent cell $\kappa_{e}^{+}$ sharing edge $e$.

###  3.1 Numerical approximation and function space

Let us consider the reference element $\text{I}^{2}=[-1,1]^{2}$ with
coordinates $\bm{\xi}=(\xi,\eta)$ and the reference edge $\text{I}=[-1,1]$,
the functions ${\bf x}_{\kappa}(\bm{\xi})$ and ${\bf x}_{e}(\xi)$ map
reference to physical element and edge, respectively. The approximate solution
of (eq. 1) is sought in the function space of piecewise polynomials

$\mathcal{V}^{p}_{h}=\big{\\{}\phi\in
L^{2}(\Omega_{h})\mathrel{\mathop{\mathchar 58\relax}}\phi|_{\kappa}\circ{\bf
x}_{\kappa}(\bm{\xi})\in\mathcal{Q}^{p}(\text{I}^{2}),\forall\kappa\in\Omega_{h}\big{\\}},$

where $\mathcal{Q}^{p}(\text{I}^{2})$ denotes the space of polynomials over
the reference element $\text{I}^{2}$ formed by the tensor product of
polynomials of degree at most $p$ in each direction. The approximate solution
reads

${\textbf{u}}_{h}({\bf x},t)\mathrel{\mathop{\mathchar
58\relax}}=\sum^{p}_{i,j=0}\phi^{ij}_{\kappa}({\bf
x}){\textbf{U}}^{ij}_{\kappa}(t)\quad\forall{\bf
x}\in\kappa,\quad\kappa\in\Omega_{h},\quad t\geqslant 0,$ (16)

where $\\{\phi^{ij}_{\kappa}\\}_{0\leqslant i,j\leqslant p}$ constitutes a
basis of $\mathcal{V}^{p}_{h}$ restricted onto $\kappa$, with dimension
$(p+1)^{2}$, and $\\{{\textbf{U}}^{ij}_{\kappa}\\}_{0\leqslant i,j\leqslant
p}$ are the DOFs. Let $\ell_{0\leqslant k\leqslant p}$ denote the Lagrange
interpolation polynomials associated to the Gauss-Lobatto nodes over I:
$-1=\xi_{0}<\xi_{1}<\dots<\xi_{p}=1$. We define the basis functions as the
tensor products of these polynomials:

$\phi^{ij}_{\kappa}({\bf x})=\phi^{ij}({\bf
x}_{\kappa}(\bm{\xi}))=\ell_{i}(\xi)\ell_{j}(\eta),\quad 0\leqslant
i,j\leqslant p,$ (17)

which satisfy the following cardinality relation

$\phi^{ij}_{\kappa}({\bf
x}^{i^{\prime}j^{\prime}}_{\kappa})=\phi^{ij}_{\kappa}({\bf
x}_{\kappa}(\bm{\xi}_{i^{\prime}j^{\prime}}))=\ell_{i}(\xi_{i^{\prime}})\ell_{j}(\eta_{j^{\prime}})=\delta_{ii^{\prime}}\delta_{jj^{\prime}},\quad
0\leqslant i,i^{\prime},j,j^{\prime}\leqslant p,$ (18)

where $\delta_{ii^{\prime}}$ is the Kronecker delta. The DOFs are therefore
point values of the solution:
${\textbf{U}}^{ij}_{\kappa}(t)={\textbf{u}}_{h}({\bf x}^{ij}_{\kappa},t)$.
Likewise, the elements $\kappa$ are interpolated on the same grid of
quadrature points as the numerical solution, i.e., ${\bf
x}_{\kappa}(\bm{\xi})=\sum_{0\leqslant i,j\leqslant
p}\ell_{i}(\xi)\ell_{j}(\eta){\bf x}_{\kappa}^{ij}$ and
${\textbf{x}}_{e}(\xi)=\sum_{k=0}^{p}\ell_{i}(\xi){\textbf{x}}_{e}^{k}$ (see
Figure 1).

The integrals over elements and interfaces are approximated by using Gauss-
Lobatto quadrature rules where the quadrature and interpolation points are
collocated:

$\int_{\kappa}f({\bf
x})\mathrm{d}V\approx\sum^{p}_{i,j=0}\omega_{i}\omega_{j}J^{ij}_{\kappa}f({\bf
x}^{ij}_{\kappa}),\quad\int_{e}f({\bf
x})dS\approx\sum^{p}_{k=0}\omega_{k}J_{e}^{k}f({\bf x}^{k}_{e}),$ (19)

where $\omega_{i},\omega_{j}>0$ are the quadrature weights and
$J^{ij}_{\kappa}=J_{\kappa}({\textbf{x}}^{ij}_{\kappa})=|{\bf
x}_{\kappa}^{\prime}(\bm{\xi}_{ij})|$, and
$J_{e}^{k}=J_{e}({\textbf{x}}_{e}^{k})=|{\textbf{x}}_{e}^{\prime}(\xi_{k})|$.
The cell-averaged solution thus reads

$\langle{\textbf{u}}\rangle_{\kappa}(t)\mathrel{\mathop{\mathchar
58\relax}}=\sum^{p}_{i,j=0}\omega_{i}\omega_{j}\frac{J^{ij}_{\kappa}}{|\kappa|}{\textbf{U}}^{ij}_{\kappa}(t)\approx\frac{1}{|\kappa|}\int_{\kappa}{\textbf{u}}_{h}({\textbf{x}},t)\mathrm{d}V,$
(20)

where $|\kappa|=\sum^{p}_{i,j=0}\omega_{i}\omega_{j}J^{ij}_{\kappa}$ is the
volume of the cell $\kappa$.

We also introduce the discrete difference matrix with entries

$D_{kl}=\ell^{\prime}_{l}(\xi_{k}),\quad 0\leqslant k,l\leqslant p,$ (21)

where the property $\sum^{p}_{l=0}\ell_{l}\equiv 1$ implies

$\sum^{p}_{l=0}D_{kl}=0\quad\forall 0\leqslant k\leqslant p.$ (22)

The discrete difference matrix is known to satisfy the SBP property [42]:

$\omega_{k}D_{kl}+\omega_{l}D_{lk}=\delta_{kp}\delta_{lp}-\delta_{k0}\delta_{l0}\quad\forall
0\leqslant k,l\leqslant p.$ (23)

Finally, the discretization is assumed to satisfy the following metric
identities [41]

$\sum_{k=0}^{p}D_{ik}J_{\kappa}^{kj}\nabla\xi({\bf
x}_{\kappa}^{kj})+D_{jk}J_{\kappa}^{ik}\nabla\eta({\bf
x}_{\kappa}^{ik})=0\quad\forall 0\leqslant i,j\leqslant p,$ (24)

and volume and edge metric terms are related by the following relations when
evaluated at edges

$\displaystyle
J_{\kappa}^{pj}\nabla\xi({\textbf{x}}_{\kappa}^{pj})=J_{e}({\textbf{x}}_{\kappa}^{pj}){\bf
n}_{e}({\textbf{x}}_{\kappa}^{pj}),\quad
J_{\kappa}^{0j}\nabla\xi({\textbf{x}}_{\kappa}^{0j})=-J_{e}({\textbf{x}}_{\kappa}^{0j}){\bf
n}_{e}({\textbf{x}}_{\kappa}^{0j})\quad\forall 0\leqslant j\leqslant p,$ (25)
$\displaystyle
J_{\kappa}^{ip}\nabla\eta({\textbf{x}}_{\kappa}^{ip})=J_{e}({\textbf{x}}_{\kappa}^{ip}){\bf
n}_{e}({\textbf{x}}_{\kappa}^{ip}),\quad
J_{\kappa}^{i0}\nabla\eta({\textbf{x}}_{\kappa}^{i0})=-J_{e}({\textbf{x}}_{\kappa}^{i0}){\bf
n}_{e}({\textbf{x}}_{\kappa}^{i0})\quad\forall 0\leqslant i\leqslant p,$

where ${\bf n}_{e}$ denotes the unit normal to $e$ in $\partial\kappa$
pointing outward from $\kappa$ (see Figure 1).

###  3.2 Semi-discrete form

We start with the following semi-discrete weak form of (eq. 1a) (see, e.g.,
[25, 15, 53]): find ${\textbf{u}}_{h}$ in $({\cal V}_{h}^{p})^{n_{eq}}$ such
that

$\sum_{\kappa\in\Omega_{h}}\int_{\kappa}v_{h}\big{(}\partial_{t}{\textbf{u}}_{h}+\nabla\cdot\textbf{f}({\textbf{u}}_{h})+\textbf{c}({\textbf{u}}_{h})\nabla{\textbf{u}}_{h}\big{)}\mathrm{d}V+\sum_{e\in\mathcal{E}_{h}}\int_{e}v_{h}^{-}{\textbf{D}}^{-}({\textbf{u}}^{-}_{h},{\textbf{u}}^{+}_{h},{\textbf{n}}_{e})+v_{h}^{+}{\textbf{D}}^{+}({\textbf{u}}^{-}_{h},{\textbf{u}}^{+}_{h},{\textbf{n}}_{e})\mathrm{d}S=0\quad\forall
v_{h}\in{\cal V}_{h}^{p},$

where ${\textbf{u}}^{\pm}_{h}({\bf x},t)=\lim_{\varepsilon\downarrow
0}{\textbf{u}}_{h}({\bf x}\pm\varepsilon{\textbf{n}}_{e}({\textbf{x}}),t)$ are
the traces of ${\textbf{u}}_{h}$ at x on a given cell interface
$e\in\mathcal{E}_{h}$ (see Figure 1) and
${\textbf{D}}^{\pm}(\cdot,\cdot,\cdot)$ are the numerical fluctuations that
are applied at the interfaces and are considered under the form

$\displaystyle{\textbf{D}}^{-}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})$
$\displaystyle=\textbf{h}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-\textbf{f}({\textbf{u}}^{-})\cdot{\textbf{n}}+\textbf{d}^{-}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}}),$
(26a)
$\displaystyle{\textbf{D}}^{+}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})$
$\displaystyle=\textbf{f}({\textbf{u}}^{+})\cdot{\textbf{n}}-\textbf{h}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+\textbf{d}^{+}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}}),$
(26b)

to allow proper discretizations of each term in (eq. 1a). They satisfy the
consistency relations

$\textbf{h}({\textbf{u}},{\textbf{u}},{\textbf{n}})=\textbf{f}({\textbf{u}})\cdot{\textbf{n}},\quad\textbf{d}^{\pm}({\textbf{u}},{\textbf{u}},{\textbf{n}})=0\quad\forall{\textbf{u}}\in\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm
GM$}},$ (27)

and will be introduced in section 5.

The integrals over mesh elements should be modified where we replace the
physical fluxes and nonconservative products with numerical fluctuations of
the form [10, 11, 15, 53]

$\displaystyle{\textbf{D}}^{-}_{X}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})$
$\displaystyle=\textbf{h}_{X}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-\textbf{f}({\textbf{u}}^{-})\cdot{\textbf{n}}+\textbf{d}_{X}^{-}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}}),$
(28a)
$\displaystyle{\textbf{D}}^{+}_{X}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})$
$\displaystyle=\textbf{f}({\textbf{u}}^{+})\cdot{\textbf{n}}-\textbf{h}_{X}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+\textbf{d}_{X}^{+}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}}),$
(28b)

where the subscript X will correspond to either entropy conservative, X=ec, or
contact preserving, X=cp, numerical fluctuations, respectively, that will be
introduced in section 4.

We then substitute $v_{h}$ for the Lagrange interpolation polynomials (eq. 17)
and consider the quadrature rules (eq. 19). Using the discrete difference
matrix (eq. 21), the semi-discrete problem reads: find ${\textbf{u}}_{h}$ in
(eq. 16) such that

$\displaystyle\omega_{i}\omega_{j}J^{ij}_{\kappa}\frac{\mathrm{d}}{\mathrm{d}t}{\textbf{U}}^{ij}_{\kappa}$
$\displaystyle+\omega_{i}\omega_{j}\sum^{p}_{k=0}\left(D_{ik}\tilde{{\textbf{D}}}_{X}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})+D_{jk}\tilde{{\textbf{D}}}_{X}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{ik}_{\kappa},{\textbf{n}}_{i(j,k)})\right)$
(29)
$\displaystyle+\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\phi_{\kappa}^{ij}({\textbf{x}}_{e}^{k})\omega_{k}J_{e}^{k}{\textbf{D}}^{-}\big{(}{\textbf{U}}^{ij}_{\kappa},{\textbf{u}}_{h}^{+}({\textbf{x}}_{e}^{k},t),{\textbf{n}}_{e}^{k}\big{)}\quad\forall
t>0,\quad\kappa\in\Omega_{h},\quad 0\leqslant i,j\leqslant p,$

where by (eq. 18) $\phi_{\kappa}^{ij}({\textbf{x}}_{e}^{k})=1$ if
${\textbf{x}}_{\kappa}^{ij}={\textbf{x}}_{e}^{k}$ and
$\phi_{\kappa}^{ij}({\textbf{x}}_{e}^{k})=0$ else. The fluctuations in the
volume integrals read

$\displaystyle\tilde{{\textbf{D}}}_{X}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})$
$\displaystyle\mathrel{\mathop{\mathchar
58\relax}}={\textbf{D}}^{-}_{X}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-{\textbf{D}}^{+}_{X}({\textbf{u}}^{+},{\textbf{u}}^{-},{\textbf{n}}),$
(30a) $\displaystyle\underset{(\lx@cref{creftype~refnum}{eqn: chap 3
sum_Dkl_vanishes})}{\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn:
chap 5 EC fluxes abstract
form})}}}{{=}}}\textbf{h}_{X}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+\textbf{h}_{X}({\textbf{u}}^{+},{\textbf{u}}^{-},{\textbf{n}})+{\bf
d}_{X}^{-}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-{\bf
d}_{X}^{+}({\textbf{u}}^{+},{\textbf{u}}^{-},{\textbf{n}}),$ (30b)

and

${\textbf{n}}_{(i,k)j}=\frac{1}{2}\left(J^{ij}_{\kappa}\nabla\xi(\bm{\xi}_{ij})+J^{kj}_{\kappa}\nabla\xi(\bm{\xi}_{kj})\right),\quad{\textbf{n}}_{i(j,k)}=\frac{1}{2}\left(J^{ij}_{\kappa}\nabla\eta(\bm{\xi}_{ij})+J^{ik}_{\kappa}\nabla\eta(\bm{\xi}_{ik})\right)$
(31)

must be introduced to keep conservation of the physical fluxes [69] and
preserve uniform states. The numerical flux and fluctuations in (eq. 28)
satisfy the consistency conditions

$\textbf{h}_{X}({\textbf{u}},{\textbf{u}},{\textbf{n}})=\textbf{f}({\textbf{u}})\cdot{\textbf{n}},\quad\textbf{d}_{X}^{\pm}({\textbf{u}},{\textbf{u}},{\textbf{n}})=0\quad\forall{\textbf{u}}\in\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm
GM$}},$ (32)

and entropy conservative (EC) fluctuations
${\textbf{D}}_{ec}^{\pm}(\cdot,\cdot,\cdot)$ satisfy [10]

$\bm{\vartheta}({\textbf{u}}^{-})^{\top}{\textbf{D}}^{-}_{ec}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+\bm{\vartheta}({\textbf{u}}^{+})^{\top}{\textbf{D}}^{+}_{ec}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})=\llbracket{\textbf{q}}({\textbf{u}})\rrbracket\cdot{\textbf{n}}\quad\forall{\textbf{u}}^{\pm}\in\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm
GM$}},$ (33)

where $\llbracket v_{h}\rrbracket=v^{+}_{h}-v^{-}_{h}$ and
$\bm{\vartheta}({\textbf{u}})\mathrel{\mathop{\mathchar 58\relax}}=\nabla_{\bf
u}\eta({\textbf{u}})$ denotes the entropy variables. Additionally, the
interface fluctuations (eq. 26) in (eq. 29) are assumed to be entropy stable:

$\bm{\vartheta}({\textbf{u}}^{-})^{\top}{\textbf{D}}^{-}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+\bm{\vartheta}({\textbf{u}}^{+})^{\top}{\textbf{D}}^{+}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})\geqslant\llbracket{\textbf{q}}({\textbf{u}})\rrbracket\cdot{\textbf{n}}\quad\forall{\textbf{u}}^{\pm}\in\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm
GM$}}.$ (34)

Finally, the initial condition (eq. 1b) is projected onto the function space
by ${\textbf{U}}^{ij}_{\kappa}(0)={\textbf{u}}_{0}(x^{ij}_{\kappa})$ for all
$\kappa\in\Omega_{h}$ and $0\leqslant i,j\leqslant p$. The theorem below
summarizes the main properties of the semi-discrete scheme (eq. 29) for the
discretization of general systems of the form (eq. 1a). Points (ii) and (iv)
have been proved in [68] with slightly different volume fluctuations (see
Remark 3.1 for comments on fluctuations in (eq. 35)).

######  Theorem 3.1.

Let ${\textbf{D}}^{\pm}_{X}(\cdot,\cdot,\cdot)$ in (eq. 28) be consistent
fluctuations with ${\bf d}_{X}^{\pm}(\cdot,\cdot,\cdot)$ in (eq. 30)
satisfying

$\displaystyle{\bf
d}_{X}^{\pm}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})$
$\displaystyle=\mathcal{C}^{\pm}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})\llbracket{\textbf{u}}\rrbracket,$
(35a)
$\displaystyle\mathcal{C}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})$
$\displaystyle\mathrel{\mathop{\mathchar
58\relax}}=\mathcal{C}^{+}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+\mathcal{C}^{-}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}}),$
(35b)
$\displaystyle\mathcal{C}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+\mathcal{C}({\textbf{u}}^{+},{\textbf{u}}^{-},{\textbf{n}})$
$\displaystyle=\big{(}\textbf{c}({\textbf{u}}^{-})+\textbf{c}({\textbf{u}}^{+})\big{)}{\bf
n},$ (35c) $\displaystyle\mathcal{C}({\textbf{u}},{\textbf{u}},{\textbf{n}})$
$\displaystyle={\bf c}({\textbf{u}}){\bf n},$ (35d)

where ${\textbf{n}}=(n_{1},\dots,n_{d})^{\top}$, ${\bf c}({\textbf{u}}){\bf
n}=\sum_{i=1}^{d}n_{i}\textbf{c}_{i}({\textbf{u}})$, and
$\llbracket{\textbf{u}}\rrbracket={\textbf{u}}^{+}-{\textbf{u}}^{-}$, and let
${\textbf{D}}^{-}(\cdot,\cdot,\cdot)$ in (eq. 29) be consistent (eq. 27)
fluctuations. Then, the semi-discrete scheme (eq. 29) has the following
properties:

  1. (i)

it is a high-order accurate approximation of smooth enough solutions to (1);

  2. (ii)

it preserves uniform states (free-stream preservation);

  3. (iii)

the cell-averaged solution (eq. 20) satisfies the following cell-averaged
semi-discrete scheme

$\displaystyle|\kappa|\frac{\mathrm{d}}{\mathrm{d}t}\langle{\textbf{u}}_{h}\rangle_{\kappa}(t)$
$\displaystyle+\sum^{p}_{i,j,k=0}\omega_{i}\omega_{j}\left(D_{ik}{\bf
c}({\textbf{U}}_{\kappa}^{ij}){\textbf{n}}_{(i,k)j}{\textbf{U}}_{\kappa}^{kj}+D_{jk}{\bf
c}({\textbf{U}}_{\kappa}^{ij}){\textbf{n}}_{i(j,k)}{\textbf{U}}_{\kappa}^{ik}\right)$
(36)
$\displaystyle+\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\omega_{k}J_{e}^{k}\left({\bf
h}\big{(}{\textbf{u}}_{h}^{-}({\textbf{x}}_{e}^{k},t),{\textbf{u}}_{h}^{+}({\textbf{x}}_{e}^{k},t),{\textbf{n}}_{e}^{k}\big{)}+{\bf
d}^{-}\big{(}{\textbf{u}}_{h}^{-}({\textbf{x}}_{e}^{k},t),{\textbf{u}}_{h}^{+}({\textbf{x}}_{e}^{k},t),{\textbf{n}}_{e}^{k}\big{)}\right)\quad\forall
t>0,\quad\kappa\in\Omega_{h},$

which ensures a discretely conservative approximation of the physical fluxes
in (eq. 1a);

  4. (iv)

if the fluctuations ${\textbf{D}}_{X}(\cdot,\cdot,\cdot)$ in the volume
integral are further assumed to be EC (eq. 33) and the fluctuations at
interfaces ${\textbf{D}}(\cdot,\cdot,\cdot)$ are entropy stable (eq. 34) for a
convex entropy pair $(\eta({\textbf{u}}),{\textbf{q}}({\textbf{u}}))$, then
the following semi-discrete entropy inequality holds:

$|\kappa|\frac{\mathrm{d}}{\mathrm{d}t}\langle\eta({\textbf{u}}_{h})\rangle_{\kappa}+\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\omega_{k}J_{e}^{k}Q\big{(}{\textbf{u}}_{h}^{-}({\textbf{x}}_{e}^{k},t),{\textbf{u}}_{h}^{+}({\textbf{x}}_{e}^{k},t),{\textbf{n}}_{e}^{k}\big{)}\leqslant
0\quad\forall t>0,\quad\kappa\in\Omega_{h},$ (37)

with the consistent and conservative entropy flux

$Q({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})=\frac{1}{2}\left({\textbf{q}}({\textbf{u}}^{-})+{\textbf{q}}({\textbf{u}}^{+})\right)\cdot{\textbf{n}}+\frac{1}{2}\bm{\vartheta}({\textbf{u}}^{-})^{\top}{\textbf{D}}^{-}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-\frac{1}{2}\bm{\vartheta}({\textbf{u}}^{+})^{\top}{\textbf{D}}^{+}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}}).$
(38)

###### Proof.

Preliminary: By the metric identities (eq. 24) and (eq. 22), the metric terms
in (eq. 31) satisfy

$\sum^{p}_{k=0}D_{ik}{\textbf{n}}_{(i,k)j}+D_{jk}{\textbf{n}}_{i(j,k)}=\frac{1}{2}\sum^{p}_{k=0}D_{ik}\big{(}J_{\kappa}^{ij}\nabla\xi({\textbf{x}}_{\kappa}^{ij})+J_{\kappa}^{kj}\nabla\xi({\textbf{x}}_{\kappa}^{kj})\big{)}+D_{jk}\big{(}J_{\kappa}^{ij}\nabla\eta({\textbf{x}}_{\kappa}^{ij})+J_{\kappa}^{kj}\nabla\eta({\textbf{x}}_{\kappa}^{ik})\big{)}=0.$
(39)

High-order accuracy: It is sufficient to prove that the volume integral in
(eq. 29) is a high-order approximation of $\nabla\cdot{\bf
f}({\textbf{u}})+{\bf c}({\textbf{u}})\nabla{\textbf{u}}$ at
$({\textbf{x}}_{\kappa}^{ij},t)$ for smooth enough solutions. High-order
accuracy of the conservative term $\nabla\cdot{\bf f}({\textbf{u}})$ has been
proved in [13, 50] and the high-order accuracy of ${\bf
c}({\textbf{u}})\nabla{\textbf{u}}$ in one space dimension has been proved in
[53, Th. 3.2]. Since we are using tensor products of one-dimensional
operators, the proof of accuracy follows by considering each space dimension
independently.

Free-stream preservation: Let us assume that ${\textbf{U}}_{\kappa}^{ij}={\bf
u}$ in the space residuals for all $\kappa\in\Omega_{h}$ and $0\leqslant
i,j\leqslant p$, then by consistency: ${\textbf{D}}^{-}({\bf u},{\bf u},{\bf
n}_{e}^{k})=0$ and $\tilde{{\textbf{D}}}_{X}({\bf u},{\bf
u},{\textbf{n}})=2{\bf f}({\bf u})\cdot{\bf n}$. Further using (eq. 39), (eq.
29) becomes

$0=\omega_{i}\omega_{j}J^{ij}_{\kappa}\frac{\mathrm{d}}{\mathrm{d}t}{\textbf{U}}^{ij}_{\kappa}+\omega_{i}\omega_{j}\sum^{p}_{k=0}2{\bf
f}({\bf
u})(D_{ik}{\textbf{n}}_{(i,k)j}+D_{jk}{\textbf{n}}_{i(j,k)})=\omega_{i}\omega_{j}J^{ij}_{\kappa}\frac{\mathrm{d}}{\mathrm{d}t}{\textbf{U}}^{ij}_{\kappa}.$

Cell-averaged semi-discrete scheme: Summing up (eq. 29) over $0\leqslant
i,j\leqslant p$ gives

$|\kappa|\frac{\mathrm{d}}{\mathrm{d}t}\langle{\textbf{u}}_{h}\rangle_{\kappa}(t)+\sum^{p}_{j=0}\omega_{j}A_{j}+\sum^{p}_{i=0}\omega_{i}B_{i}+\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\omega_{k}J_{e}^{k}{\textbf{D}}^{-}\big{(}{\textbf{u}}_{h}^{-}({\textbf{x}}_{e}^{k},t),{\textbf{u}}_{h}^{+}({\textbf{x}}_{e}^{k},t),{\textbf{n}}_{e}^{k}\big{)}=0,$

where from (eq. 30), we have

$\displaystyle A_{j}$
$\displaystyle=\sum^{p}_{i,k=0}\omega_{i}D_{ik}\tilde{{\textbf{D}}}_{X}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})=\sum^{p}_{i,k=0}\omega_{i}D_{ik}\big{(}{\bf
h}_{X}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})+{\bf
h}_{X}({\textbf{U}}^{kj}_{\kappa},{\textbf{U}}^{ij}_{\kappa},{\textbf{n}}_{(i,k)j})$
$\displaystyle\hskip 170.71652pt+{\bf
d}^{-}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})-{\bf
d}^{+}({\textbf{U}}^{kj}_{\kappa},{\textbf{U}}^{ij}_{\kappa},{\textbf{n}}_{(i,k)j})\big{)}$
$\displaystyle B_{i}$
$\displaystyle=\sum^{p}_{j,k=0}\omega_{j}D_{jk}\tilde{{\textbf{D}}}_{X}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{ik}_{\kappa},{\textbf{n}}_{i(j,k)})=\sum^{p}_{j,k=0}\omega_{j}D_{jk}\big{(}{\bf
h}_{X}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{ik}_{\kappa},{\textbf{n}}_{i(j,k)})+{\bf
h}_{X}({\textbf{U}}^{ik}_{\kappa},{\textbf{U}}^{ij}_{\kappa},{\textbf{n}}_{i(j,k)})$
$\displaystyle\hskip 170.71652pt+{\bf
d}^{-}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{ik}_{\kappa},{\textbf{n}}_{i(j,k)})-{\bf
d}^{+}({\textbf{U}}^{ik}_{\kappa},{\textbf{U}}^{ij}_{\kappa},{\textbf{n}}_{i(j,k)})\big{)}$

Let us consider the first term. Using (eq. 25), we have

$\displaystyle A_{j}$ $\displaystyle\underset{(\lx@cref{creftype~refnum}{eqn:
numerical fluxes consistency
condition})}{\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn: chap
3 SBP})}}}{{=}}}J_{e}({\textbf{x}}_{\kappa}^{pj}){\bf
f}({\textbf{U}}_{\kappa}^{pj})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{pj})+J_{e}({\textbf{x}}_{\kappa}^{0j}){\bf
f}({\textbf{U}}_{\kappa}^{0j})\cdot{\bf n}_{e}({\textbf{x}}_{\kappa}^{0j})$
$\displaystyle+\sum^{p}_{i,k=0}\omega_{i}D_{ik}{\bf
h}_{X}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})-\omega_{k}D_{ki}{\bf
h}_{X}({\textbf{U}}^{kj}_{\kappa},{\textbf{U}}^{ij}_{\kappa},{\textbf{n}}_{(i,k)j})+\omega_{i}D_{ik}{\bf
d}^{-}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})+\omega_{k}D_{ki}{\bf
d}^{+}({\textbf{U}}^{kj}_{\kappa},{\textbf{U}}^{ij}_{\kappa},{\textbf{n}}_{(i,k)j})$
$\displaystyle\underset{(\lx@cref{creftype~refnum}{eqn: high-order
approximation proof}a,b)}{\overset{i\leftrightarrow
k}{=}}J_{e}({\textbf{x}}_{\kappa}^{pj}){\bf
f}({\textbf{U}}_{\kappa}^{pj})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{pj})+J_{e}({\textbf{x}}_{\kappa}^{0j}){\bf
f}({\textbf{U}}_{\kappa}^{0j})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{0j})+\sum^{p}_{i,k=0}\omega_{i}D_{ik}{\cal
C}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})({\textbf{U}}^{kj}_{\kappa}-{\textbf{U}}^{ij}_{\kappa})$
$\displaystyle\underset{(\lx@cref{creftype~refnum}{eqn: high-order
approximation proof
d})}{\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn: chap 3
SBP})}}}{{=}}}J_{e}({\textbf{x}}_{\kappa}^{pj})\big{(}{\bf
f}({\textbf{U}}_{\kappa}^{pj})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{pj})-{\bf c}({\textbf{U}}_{\kappa}^{pj}){\bf
n}_{e}({\textbf{x}}_{\kappa}^{pj}){\textbf{U}}_{\kappa}^{pj}\big{)}+J_{e}({\textbf{x}}_{\kappa}^{0j})\big{(}{\bf
f}({\textbf{U}}_{\kappa}^{0j})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{0j})-{\bf c}({\textbf{U}}_{\kappa}^{0j}){\bf
n}_{e}({\textbf{x}}_{\kappa}^{0j}){\textbf{U}}_{\kappa}^{0j}\big{)}$
$\displaystyle+\sum^{p}_{i,k=0}\omega_{i}D_{ik}{\cal
C}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j}){\textbf{U}}^{kj}_{\kappa}+\omega_{k}D_{ki}{\cal
C}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j}){\textbf{U}}^{ij}_{\kappa}$
$\displaystyle\underset{(\lx@cref{creftype~refnum}{eqn: high-order
approximation proof c})}{\overset{i\leftrightarrow
k}{=}}J_{e}({\textbf{x}}_{\kappa}^{pj})\big{(}{\bf
f}({\textbf{U}}_{\kappa}^{pj})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{pj})-{\bf c}({\textbf{U}}_{\kappa}^{pj}){\bf
n}_{e}({\textbf{x}}_{\kappa}^{pj}){\textbf{U}}_{\kappa}^{pj}\big{)}+J_{e}({\textbf{x}}_{\kappa}^{0j})\big{(}{\bf
f}({\textbf{U}}_{\kappa}^{0j})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{0j})-{\bf c}({\textbf{U}}_{\kappa}^{0j}){\bf
n}_{e}({\textbf{x}}_{\kappa}^{0j}){\textbf{U}}_{\kappa}^{0j}\big{)}$
$\displaystyle+\sum^{p}_{i,k=0}\omega_{i}D_{ik}\left({\bf
c}({\textbf{U}}_{\kappa}^{ij})+{\bf
c}({\textbf{U}}_{\kappa}^{kj})\right){\textbf{n}}_{(i,k)j}{\textbf{U}}_{\kappa}^{kj}$
$\displaystyle\underset{(\lx@cref{creftype~refnum}{eqn: high-order
approximation proof
d})}{\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn: chap 3
SBP})}}}{{=}}}J_{e}({\textbf{x}}_{\kappa}^{pj}){\bf
f}({\textbf{U}}_{\kappa}^{pj})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{pj})+J_{e}({\textbf{x}}_{\kappa}^{0j}){\bf
f}({\textbf{U}}_{\kappa}^{0j})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{0j})+\sum^{p}_{i,k=0}\omega_{i}D_{ik}{\bf
c}({\textbf{U}}_{\kappa}^{ij}){\textbf{n}}_{(i,k)j}({\textbf{U}}_{\kappa}^{kj}-{\textbf{U}}_{\kappa}^{ij}),$

where $i\leftrightarrow k$ indicates an inversion of indices $i$ and $k$ in
some of the terms. Likewise, we have

$B_{i}=J_{e}({\textbf{x}}_{\kappa}^{ip}){\bf
f}({\textbf{U}}_{\kappa}^{ip})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{ip})+J_{e}({\textbf{x}}_{\kappa}^{i0}){\bf
f}({\textbf{U}}_{\kappa}^{i0})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{i0})+\sum^{p}_{j,k=0}\omega_{j}D_{jk}{\bf
c}({\textbf{U}}_{\kappa}^{ij}){\textbf{n}}_{i(j,k)}({\textbf{U}}_{\kappa}^{ik}-{\textbf{U}}_{\kappa}^{ij}).$

From (eq. 39) we deduce

$\sum^{p}_{j=0}\omega_{j}A_{j}+\sum^{p}_{i=0}\omega_{i}B_{i}=\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\omega_{k}J_{e}^{k}{\bf
h}\big{(}{\textbf{u}}_{h}^{-}({\textbf{x}}_{e}^{k},t)\big{)}{\textbf{n}}_{e}^{k}+\sum^{p}_{i,j,k=0}\omega_{i}\omega_{j}\left(D_{ik}{\bf
c}({\textbf{U}}_{\kappa}^{ij}){\textbf{n}}_{(i,k)j}{\textbf{U}}_{\kappa}^{kj}+D_{jk}{\bf
c}({\textbf{U}}_{\kappa}^{ij}){\textbf{n}}_{i(j,k)}{\textbf{U}}_{\kappa}^{ik}\right),$

and we get (eq. 36).

Entropy stability: Let us now consider EC fluxes in the volume integral in the
semi-discrete DGSEM (eq. 29). We left multiply (29) by
$\bm{\vartheta}_{\kappa}^{ij}=\bm{\vartheta}({\textbf{u}}^{ij}_{\kappa})$ and
sum up over $0\leqslant i,j\leqslant p$ to get

$|\kappa|\frac{\mathrm{d}}{\mathrm{d}t}\langle\eta\rangle_{\kappa}(t)+\sum^{p}_{j=0}\omega_{j}C_{j}+\sum^{p}_{i=0}\omega_{i}E_{i}+\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\omega_{k}J_{e}^{k}\bm{\vartheta}\big{(}{\textbf{u}}_{h}^{-}({\textbf{x}}_{e}^{k},t)\big{)}\cdot{\textbf{D}}^{-}\big{(}{\textbf{u}}_{h}^{-}({\textbf{x}}_{e}^{k},t),{\textbf{u}}_{h}^{+}({\textbf{x}}_{e}^{k},t),{\textbf{n}}_{e}^{k}\big{)}=0,$

where

$C_{j}=\sum^{p}_{i,k=0}\omega_{i}D_{ik}\bm{\vartheta}_{\kappa}^{ij}\cdot\left({\bf
D}^{-}_{ec}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})-{\bf
D}^{+}_{ec}({\textbf{U}}^{kj}_{\kappa},{\textbf{U}}^{ij}_{\kappa},{\textbf{n}}_{(i,k)j})\right),$

and

$E_{i}=\sum^{p}_{j,k=0}\omega_{j}D_{jk}\bm{\vartheta}_{\kappa}^{ij}\cdot\left({\bf
D}^{-}_{ec}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{ik}_{\kappa},{\textbf{n}}_{i(j,k)})-{\bf
D}^{+}_{ec}({\textbf{U}}^{ik}_{\kappa},{\textbf{U}}^{ij}_{\kappa},{\textbf{n}}_{i(j,k)})\right).$

Using (eq. 33), we have

$\displaystyle C_{j}$
$\displaystyle=\sum^{p}_{i,k=0}\omega_{i}D_{ik}\left(\bm{\vartheta}_{\kappa}^{ij}\cdot{\bf
D}^{-}_{ec}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})+\bm{\vartheta}_{\kappa}^{kj}\cdot{\bf
D}^{-}_{ec}({\textbf{U}}^{kj}_{\kappa},{\textbf{U}}^{ij}_{\kappa},{\textbf{n}}_{(i,k)j})-\big{(}{\textbf{q}}({\textbf{U}}^{ij}_{\kappa})-{\textbf{q}}({\textbf{U}}^{kj}_{\kappa})\big{)}\cdot{\textbf{n}}_{(i,k)j}\right)$
$\displaystyle\underset{(\lx@cref{creftype~refnum}{eqn: numerical fluxes
consistency
condition})}{\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn: chap
3 SBP})}}}{{=}}}J_{e}({\textbf{x}}_{\kappa}^{pj}){\bf
q}({\textbf{U}}_{\kappa}^{pj})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{pj})+J_{e}({\textbf{x}}_{\kappa}^{0j}){\bf
q}({\textbf{U}}_{\kappa}^{0j})\cdot{\bf n}_{e}({\textbf{x}}_{\kappa}^{0j})$
$\displaystyle+\sum^{p}_{i,k=0}\omega_{i}D_{ik}\bm{\vartheta}_{\kappa}^{ij}\cdot{\bf
D}^{-}_{ec}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})-\omega_{k}D_{ki}\bm{\vartheta}_{\kappa}^{kj}\cdot{\bf
D}^{-}_{ec}({\textbf{U}}^{kj}_{\kappa},{\textbf{U}}^{ij}_{\kappa},{\textbf{n}}_{(i,k)j})-\big{(}\omega_{i}D_{ik}{\textbf{q}}({\textbf{U}}^{ij}_{\kappa})+\omega_{k}D_{ki}{\textbf{q}}({\textbf{U}}^{kj}_{\kappa})\big{)}\cdot{\textbf{n}}_{(i,k)j}$
$\displaystyle\underset{(\lx@cref{creftype~refnum}{eqn: chap 3
sum_Dkl_vanishes})}{\overset{i\leftrightarrow
k}{=}}J_{e}({\textbf{x}}_{\kappa}^{pj}){\bf
q}({\textbf{U}}_{\kappa}^{pj})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{pj})+J_{e}({\textbf{x}}_{\kappa}^{0j}){\bf
q}({\textbf{U}}_{\kappa}^{0j})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{0j})-\sum^{p}_{i,k=0}\omega_{i}D_{ik}J_{\kappa}^{kj}{\textbf{q}}({\textbf{U}}^{ij}_{\kappa})\nabla\xi({\textbf{x}}_{\kappa}^{kj}).$

Likewise

$E_{i}=J_{e}({\textbf{x}}_{\kappa}^{ip}){\bf
q}({\textbf{U}}_{\kappa}^{ip})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{ip})+J_{e}({\textbf{x}}_{\kappa}^{i0}){\bf
q}({\textbf{U}}_{\kappa}^{i0})\cdot{\bf
n}_{e}({\textbf{x}}_{\kappa}^{i0})-\sum^{p}_{j,k=0}\omega_{j}D_{jk}J_{\kappa}^{ik}{\textbf{q}}({\textbf{U}}^{ij}_{\kappa})\nabla\eta({\textbf{x}}_{\kappa}^{ik}),$

and again using the metric identities (eq. 24), we get

$|\kappa|\frac{\mathrm{d}}{\mathrm{d}t}\langle\eta\rangle_{\kappa}(t)+\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\omega_{k}J_{e}^{k}\left(\bm{\vartheta}({\textbf{x}}_{e}^{k})\cdot{\textbf{D}}^{-}\big{(}{\textbf{u}}_{h}^{-}({\textbf{x}}_{e}^{k},t),{\textbf{u}}_{h}^{+}({\textbf{x}}_{e}^{k},t),{\textbf{n}}_{e}^{k}\big{)}+{\textbf{q}}\big{(}{\textbf{u}}_{h}^{-}({\textbf{x}}_{e}^{k},t)\big{)}\cdot{\bf
n}_{e}^{k}\right)=0.$

Then, from (eq. 38) we have

$\bm{\vartheta}^{-}\cdot{\textbf{D}}^{-}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+{\textbf{q}}^{-}\cdot{\textbf{n}}=Q({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+\frac{1}{2}({\textbf{q}}^{-}-{\textbf{q}}^{+})\cdot{\textbf{n}}+\frac{1}{2}\bm{\vartheta}^{-}\cdot{\textbf{D}}^{-}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+\frac{1}{2}\bm{\vartheta}^{+}\cdot{\textbf{D}}^{+}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})$

and since the sum of the three last terms are non-negative by (eq. 34), we
obtain the desired entropy inequality (eq. 37). ∎

######  Remark 3.1.

Fluctuations with the following expressions fall into the category (eq. 35a):

$\mathcal{C}^{\pm}({\bf u}^{-},{\bf u}^{+},{\bf
n})=\frac{1}{2}\big{(}\alpha{\bf c}({\bf u}^{\pm})+(1-\alpha){\bf c}({\bf
u}^{\mp})\big{)}{\bf n},\quad 0\leqslant\alpha\leqslant 1,$

which belongs to the Volpert path family of schemes [66] and correspond to the
skew-symmetric splitting $\alpha{\bf c}\nabla{\bf
u}+(1-\alpha)(\nabla\cdot({\bf c}{\bf u})-(\nabla\cdot{\bf c}){\bf u})$ of the
nonconservative product. Relations (eq. 35b) and (eq. 35d) indeed correspond
to the consistency condition of the Volpert path family of schemes [66], while
(eq. 35c) is only necessary to get (eq. 36) which will be useful to prove
Theorem 6.1.

In the next two sections, we describe the fluctuations we use in the present
DGSEM (eq. 29) for the discretization of the SG-gamma model (eq. 1) and (eq.
2). In section 4, we first focus on designing CP and EC fluctuations that are
applied in the volume integral. We then design in section 5 a HLLC approximate
Riemann solver for (eq. 1a) that is applied at interfaces.

##  4 Numerical fluctuations for the volume integrals

In this section we focus on numerical fluctuations (eq. 28) for the scheme
(eq. 29) that will be applied in the volume integrals. We will make use of the
Leibniz identities, which we recall here: let
$a^{+},a^{-},b^{+},b^{-},c^{+},c^{-}$ in $\mathbb{R}$ and have finite values,
then we have

$\displaystyle\llbracket ab\rrbracket=\hskip 0.12135pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$a$}\hss}{\footnotesize{\hskip-0.12135pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\llbracket
b\rrbracket+\hskip 0.56876pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$b$}\hss}{\footnotesize{\hskip-0.56876pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\llbracket
a\rrbracket,\quad\llbracket abc\rrbracket=\hskip 0.12135pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$a$}\hss}{\footnotesize{\hskip-0.12135pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\big{(}\hskip
0.56876pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$b$}\hss}{\footnotesize{\hskip-0.56876pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\llbracket
c\rrbracket+\hskip 0.5526pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$c$}\hss}{\footnotesize{\hskip-0.5526pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\llbracket
b\rrbracket\big{)}+\hbox to0.0pt{\hskip
1.37863pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\footnotesize{\leavevmode\hbox{\set@color$bc$}}}\llbracket
a\rrbracket,$ (40)

where $\displaystyle\hskip 0.12135pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$a$}\hss}{\footnotesize{\hskip-0.12135pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}=\frac{a^{+}+a^{-}}{2}$
is the arithmetic mean and $\llbracket a\rrbracket=a^{+}-a^{-}$ the jump.

###  4.1 Contact preserving numerical fluxes

$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\Omega_{L}$$\Omega_{R}$$\Omega_{R}$
Figure 2: The above representation shows two cells, where the left cell
contains a material interface (line in color) which separates the states
$\Omega_{L}=\\{\alpha_{1L},\rho_{L},\Gamma_{L},\pi_{L}\\}$ and
$\Omega_{R}=\\{\alpha_{1R},\rho_{R},\Gamma_{R},\pi_{R}\\}$ with discontinuous
density and EOS parameters and uniform velocity v and pressure $\mathrm{p}$.
Note that the present scheme does not require the material discontinuity to be
aligned with mesh interfaces to preserve uniform pressure and velocity
profiles across it.

Here, we focus on deriving conditions that will ensure that the numerical
fluxes maintain uniform pressure and velocity profiles across an isolated
material interface. To this purpose, we introduce CP fluctuations in (eq. 28),
with X=cp, where
$\textbf{h}_{cp}=(h^{\rho}_{cp},\textbf{h}^{\rho{\textbf{v}}\top}_{cp},h^{\rho
E}_{cp},0,0)^{\top}$ is the vector of numerical fluxes for the conservative
equations of mass, momentum and energy, and
$\textbf{d}_{cp}^{\pm}=(0,0,0,d^{\pm}_{\Gamma},d^{\pm}_{\Pi})^{\top}$ is the
vector for the fluctuations for the nonconservative products in (eq. 1a).

For the sake of brevity, we here focus on volume fluctuations only, though the
same relations may be derived for the interface fluctuations by following the
same lines (see section 5.4.2). Instead of using the symmetrizer in (eq. 30),
we assume that the numerical fluxes are symmetric: ${\bf
h}_{cp}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})={\bf
h}_{cp}({\textbf{u}}^{+},{\textbf{u}}^{-},{\textbf{n}})$ without loss of
generality (they will be, see proposition 4.1). Let us represent the conserved
and nonconserved quantities in (eq. 1a) using
$A\in\\{\rho,\rho{\textbf{v}},\rho E\\}$ and $B\in\\{\Gamma,\Pi\\}$,
respectively. The discretization under the DGSEM for (eq. 1a) now reads

$\displaystyle\omega_{i}\omega_{j}J^{ij}_{\kappa}\frac{d}{dt}A^{ij}_{\kappa}$
$\displaystyle+2\omega_{i}\omega_{j}\sum^{p}_{k=0}\left(D_{ik}\textbf{h}^{A}_{cp}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})+D_{jk}\textbf{h}^{A}_{cp}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{ik}_{\kappa},{\textbf{n}}_{i(k,j)})\right)$
$\displaystyle+$
$\displaystyle\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\phi_{\kappa}^{ij}({\textbf{x}}_{e}^{k})\omega_{k}J_{e}^{k}{\textbf{D}}_{A}^{-}\big{(}{\textbf{U}}^{ij}_{\kappa},{\textbf{u}}_{h}^{+}({\textbf{x}}_{e}^{k},t),{\textbf{n}}_{e}^{k}\big{)}=0$
(41a)
$\displaystyle\omega_{i}\omega_{j}J^{ij}_{\kappa}\frac{d}{dt}B^{ij}_{\kappa}$
$\displaystyle+\omega_{i}\omega_{j}\sum^{p}_{k=0}\bigg{(}D_{ik}\big{(}\textbf{d}_{cp}^{B,-}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},\textbf{n}_{(i,k)j})-\textbf{d}_{cp}^{B,+}({\textbf{U}}^{kj}_{\kappa},{\textbf{U}}^{ij}_{\kappa},\textbf{n}_{(i,k)j})\big{)}$
$\displaystyle+$ $\displaystyle
D_{jk}\big{(}\textbf{d}_{cp}^{B,-}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{ik}_{\kappa},\textbf{n}_{i(k,j)})-\textbf{d}_{cp}^{B,+}({\textbf{U}}^{ik}_{\kappa},{\textbf{U}}^{ij}_{\kappa},\textbf{n}_{i(k,j)})\big{)}\bigg{)}+\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\phi_{\kappa}^{ij}({\textbf{x}}_{e}^{k})\omega_{k}J_{e}^{k}{\textbf{D}}_{B}^{-}\big{(}{\textbf{U}}^{ij}_{\kappa},{\textbf{u}}_{h}^{+}({\textbf{x}}_{e}^{k},t),{\textbf{n}}_{e}^{k}\big{)}=0$
(41b)

Now let us suppose that the initial condition consists of a material interface
with uniform velocity, v, and pressure, $\mathrm{p}$, and states $\rho_{L}$,
$\Gamma_{L}$ and $\pi_{L}$ in $\Omega_{L}$ and $\rho_{R}$, $\Gamma_{R}$ and
$\pi_{R}$ in $\Omega_{R}$ with
$\overline{\Omega_{L}}\cup\overline{\Omega_{R}}=\overline{\Omega}$, then so do
the DOFs, see Figure 2. We now derive conditions for the numerical fluxes (eq.
28) to preserve the uniform states in time.

We, first, focus on the velocity state and impose the semi-discrete scheme
(eq. 41) to satisfy a discrete counterpart to the differential relation
$\rho\mathrm{d}{\textbf{v}}=\mathrm{d}\rho{\textbf{v}}-{\textbf{v}}\mathrm{d}\rho=0$.
Ignoring the interface fluxes,
$\omega_{i}\omega_{j}J_{\kappa}^{ij}\rho_{\kappa}^{ij}\mathrm{d}_{t}{\textbf{v}}_{\kappa}^{ij}=0$
requires

$\omega_{i}\omega_{j}\sum^{p}_{k=0}\Big{(}\big{(}\textbf{h}^{\rho{\textbf{v}}}_{cp}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},\textbf{n}_{(i,k)j})-{\textbf{v}}h^{\rho}_{cp}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{kj}_{\kappa},\textbf{n}_{(i,k)j})\big{)}D_{ik}+\big{(}\textbf{h}^{\rho{\textbf{v}}}_{cp}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{ik}_{\kappa},\textbf{n}_{i(k,j)})-{\textbf{v}}h^{\rho}_{cp}({\textbf{U}}^{ij}_{\kappa},{\textbf{U}}^{ik}_{\kappa},\textbf{n}_{i(k,j)})\big{)}D_{jk}\Big{)}=0,$

and a sufficient condition reads

$\textbf{h}^{\rho{\textbf{v}}}_{cp}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})=\tilde{{\textbf{v}}}({\textbf{u}}^{-},{\textbf{u}}^{+})h^{\rho}_{cp}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+\tilde{\mathrm{p}}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}}),$
(42)

where $\tilde{{\textbf{v}}}$ and $\tilde{\mathrm{p}}$ are any consistent
discretizations of the velocity vector v and pressure $\mathrm{p}{\bf n}$.
Similarly, a semi-discrete equation for the pressure (eq. 5) can be obtained
by using $\Gamma\mathrm{d}\mathrm{p}=\mathrm{d}\rho
E-(\tfrac{1}{2}{\textbf{v}}\cdot{\textbf{v}})\mathrm{d}\rho-\mathrm{p}\mathrm{d}\Gamma-\mathrm{d}\Pi$
from (eq. 6), and again ignoring surface contributions in (eq. 41),
$\omega_{i}\omega_{j}J_{\kappa}^{ij}\Gamma_{\kappa}^{ij}\mathrm{d}_{t}\mathrm{p}({\textbf{U}}_{\kappa}^{ij})=0$
requires

$\displaystyle\omega_{i}\omega_{j}\sum^{p}_{k=0}$ $\displaystyle
D_{ik}\bigg{(}2h^{\rho
E}_{cp}(\textbf{U}^{ij}_{\kappa},\textbf{U}^{kj}_{\kappa},\textbf{n}_{(i,k)j})-{\textbf{v}}\cdot{\textbf{v}}h^{\rho}_{cp}(\textbf{U}^{ij}_{\kappa},\textbf{U}^{kj}_{\kappa},{\textbf{n}}_{(i,k)j})-\mathrm{p}\left(d^{-}_{\Gamma}(\textbf{U}^{ij}_{\kappa},\textbf{U}^{kj}_{\kappa},\textbf{n}_{(i,k)j})-d^{+}_{\Gamma}(\textbf{U}^{kj}_{\kappa},\textbf{U}^{ij}_{\kappa},\textbf{n}_{(i,k)j})\right)$
$\displaystyle-d^{-}_{\Pi}(\textbf{U}^{ij}_{\kappa},\textbf{U}^{kj}_{\kappa},\textbf{n}_{(i,k)j})+d^{+}_{\Pi}(\textbf{U}^{kj}_{\kappa},\textbf{U}^{ij}_{\kappa},\textbf{n}_{(i,k)j})\bigg{)}$
$\displaystyle+$ $\displaystyle D_{jk}\bigg{(}2h^{\rho
E}_{cp}(\textbf{U}^{ij}_{\kappa},\textbf{U}^{ik}_{\kappa},\textbf{n}_{i(k,j)})-{\textbf{v}}\cdot{\textbf{v}}h^{\rho}_{cp}(\textbf{U}^{ij}_{\kappa},\textbf{U}^{ik}_{\kappa},{\textbf{n}}_{i(j,k)})-\mathrm{p}\left(d^{-}_{\Gamma}(\textbf{U}^{ij}_{\kappa},\textbf{U}^{ik}_{\kappa},\textbf{n}_{i(k,j)})-d^{+}_{\Gamma}(\textbf{U}^{ik}_{\kappa},\textbf{U}^{ij}_{\kappa},\textbf{n}_{i(k,j)})\right)$
$\displaystyle-d^{-}_{\Pi}(\textbf{U}^{ij}_{\kappa},\textbf{U}^{ik}_{\kappa},\textbf{n}_{i(k,j)})+d^{+}_{\Pi}(\textbf{U}^{ik}_{\kappa},\textbf{U}^{ij}_{\kappa},\textbf{n}_{i(k,j)})\bigg{)}=0,$

and subtracting the trivial quantity
$2\omega_{i}\omega_{j}\sum^{p}_{k=0}\big{(}{\bf f}^{\rho
E}(\textbf{U}^{ij}_{\kappa})-\tfrac{{\textbf{v}}\cdot{\textbf{v}}}{2}{\bf
f}^{\rho}(\textbf{U}^{ij}_{\kappa})\big{)}\cdot(D_{ik}{\textbf{n}}_{(i,k)j}+D_{jk}{\textbf{n}}_{i(j,k)})=0$,
from (eq. 39), a sufficient condition reads

$h^{\rho E}_{cp}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-(\rho
E^{-}+\mathrm{p}){\textbf{v}}\cdot{\textbf{n}}=\frac{{\textbf{v}}\cdot{\textbf{v}}}{2}\big{(}h^{\rho}_{cp}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-\rho^{-}{\textbf{v}}\cdot{\textbf{n}}\big{)}\\\
+\frac{\mathrm{p}}{2}\left(d^{-}_{\Gamma}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-d^{+}_{\Gamma}({\textbf{u}}^{+},{\textbf{u}}^{-},{\textbf{n}})\right)+\frac{1}{2}\big{(}d^{-}_{\Pi}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-d^{+}_{\Pi}({\textbf{u}}^{+},{\textbf{u}}^{-},{\textbf{n}})\big{)}.$
(43)

We can now propose CP fluxes for the volume integral.

######  Proposition 4.1.

Numerical fluxes of the form (eq. 28) where

${\textbf{h}}_{cp}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})=\begin{pmatrix}\hskip
0.17343pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\rho$}\hss}{\footnotesize{\hskip-0.17343pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\,\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}}\\\
\hskip 0.17343pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\rho$}\hss}{\footnotesize{\hskip-0.17343pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\,\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\,\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}}+\hskip
0.0pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\mathrm{p}$}\hss}{\footnotesize{\hskip
0.0pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}{\textbf{n}}\\\ (\hbox
to0.0pt{\hskip
3.4078pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\footnotesize{\leavevmode\hbox{\set@color$\rho
E$}}}+\hskip 0.0pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\mathrm{p}$}\hss}{\footnotesize{\hskip
0.0pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}})\hskip 0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}}\\\
0\\\
0\end{pmatrix},\quad\textbf{d}_{cp}^{\pm}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})=\frac{1}{2}{\textbf{v}}^{\pm}\cdot{\textbf{n}}\begin{pmatrix}0\\\
0\\\ 0\\\ \llbracket\Gamma\rrbracket\\\ \llbracket\Pi\rrbracket\end{pmatrix},$
(44)

preserve the uniform pressure and velocity fields across contact
discontinuities and material interfaces for the SG-gamma model (eq. 1a) and
(eq. 2), with the mixture EOS (eq. 5).

###### Proof.

Checking that condition (eq. 42) holds for (eq. 44) is direct, wile using the
mixture EOS (eq. 5), we get

$\frac{\mathrm{p}}{2}\left(d^{-}_{\Gamma}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-d^{+}_{\Gamma}({\textbf{u}}^{+},{\textbf{u}}^{-},{\textbf{n}})\right)+\frac{1}{2}\big{(}d^{-}_{\Pi}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})-d^{+}_{\Pi}({\textbf{u}}^{+},{\textbf{u}}^{-},{\textbf{n}})\big{)}=\frac{{\textbf{v}}\cdot{\textbf{n}}}{2}\big{(}\mathrm{p}\llbracket\Gamma\rrbracket+\llbracket\Pi\rrbracket\big{)}\overset{(\lx@cref{creftype~refnum}{eqn:
chap 5 Leibniz
rule})}{=}\frac{{\textbf{v}}\cdot{\textbf{n}}}{2}\left(\llbracket\rho
E\rrbracket-\frac{{\textbf{v}}\cdot{\textbf{v}}}{2}\llbracket\rho\rrbracket\right),$

so (eq. 43) holds as well. ∎

######  Remark 4.1.

The CP numerical fluxes (eq. 44) are similar to the one proposed in [40]. Here
we have modified the contributions towards the energy equation to satisfy (eq.
43).

###  4.2 Entropy conservative numerical fluxes

We, now, propose EC fluxes for the SG-gamma model that are applied to the
modified volume integral in (eq. 29) and, according to Theorem 3.1, these
numerical fluxes will contribute to the entropy stability of the numerical
scheme.

######  Proposition 4.2.

Consider the entropy pair (eq. 11) with (eq. 12), then fluctuations of the
form (eq. 28) with

${\textbf{h}}_{ec}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})=\begin{pmatrix}\hat{\rho}\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}}\\\
{\bf h}^{\rho{\textbf{v}}}_{ec}\\\ h^{\rho E}_{ec}\\\ 0\\\
0\end{pmatrix},\quad\textbf{d}_{ec}^{\pm}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})=\frac{1}{2}{\textbf{v}}^{\pm}\cdot{\textbf{n}}\begin{pmatrix}0\\\
0\\\ 0\\\ \llbracket\Gamma\rrbracket\\\ \llbracket\Pi\rrbracket\end{pmatrix},$
(45)

where

${\bf
h}^{\rho{\textbf{v}}}_{ec}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})=\hat{\rho}(\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}})\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}+\left(\hskip
0.17343pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\rho$}\hss}{\footnotesize{\hskip-0.17343pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\hbox
to0.0pt{\hskip
9.07979pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\footnotesize{\leavevmode\hbox{\set@color$(\gamma-1)$}}}\hbox
to0.0pt{\hskip
5.58pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\footnotesize{\leavevmode\hbox{\set@color$\left(\displaystyle\frac{\mathrm{C_{v}}}{\zeta}\right)$}}}-\hbox
to0.0pt{\hskip
2.52pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\footnotesize{\leavevmode\hbox{\set@color$\mathrm{p}_{\infty}$}}}\right){\textbf{n}},\quad
h^{\rho
E}_{ec}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})=\left(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\left(\displaystyle\frac{\mathrm{C_{v}}}{\zeta}\right)")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{3.87498pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{\left(\displaystyle\frac{\mathrm{C_{v}}}{\zeta}\right)}{\tmpbox}+\frac{{\textbf{v}}^{-}\cdot{\textbf{v}}^{+}}{2}\right)\hat{\rho}\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}}+\left(\hskip
0.17343pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\rho$}\hss}{\footnotesize{\hskip-0.17343pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\hbox
to0.0pt{\hskip
9.07979pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\footnotesize{\leavevmode\hbox{\set@color$(\gamma-1)$}}}\hbox
to0.0pt{\hskip
5.58pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\footnotesize{\leavevmode\hbox{\set@color$\left(\displaystyle\frac{\mathrm{C_{v}}}{\zeta}\right)$}}}\right)\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}},$

and
$\tfrac{\mathrm{C_{v}}}{\zeta}=\tfrac{\Gamma}{\rho}(\mathrm{p}+\tfrac{\Pi}{\Gamma+1})$
from (eq. 14), are EC in the sense (eq. 33) when excluding material
interfaces, i.e.,
$\llbracket\Gamma\rrbracket\equiv\llbracket\Pi\rrbracket\equiv 0$.

###### Proof.

As we consider pure phases, the system (eq. 1a) is conservative and (eq. 33)
reduces to the Tadmor condition [60]

$\textbf{h}_{ec}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})\cdot\llbracket\bm{\vartheta}\rrbracket-\llbracket\bm{\psi}({\textbf{u}})\rrbracket\cdot{\textbf{n}}=0\quad\forall{\textbf{u}}^{\pm}\in\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm
GM$}},$

where the entropy variables $\bm{\vartheta}$ are defined in (eq. 15), and
$\bm{\psi}\equiv\textbf{f}^{\top}\bm{\vartheta}-\textbf{q}$ is the entropy
potential and reads

$\bm{\psi}({\textbf{u}})\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn:
chap 5 entropy variables})}}}{{=}}-\zeta(\rho
E+\mathrm{p}){\textbf{v}}+\zeta{\textbf{v}}(\rho{\textbf{v}}\cdot{\textbf{v}}+\mathrm{p})+\left(\gamma\mathrm{C_{v}}-s-\zeta\frac{{\textbf{v}}\cdot{\textbf{v}}}{2}\right)\rho{\textbf{v}}+\rho
s{\textbf{v}}=(\gamma\mathrm{C_{v}}-\zeta
e)\rho{\textbf{v}}=\left((\gamma-1)\mathrm{C_{v}}\rho-\mathrm{p}_{\infty}\zeta\right){\textbf{v}},$

so

$\llbracket\bm{\psi}({\textbf{u}})\rrbracket\cdot{\textbf{n}}\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn:
chap 5 Leibniz rule})}}}{{=}}\big{(}\hbox to0.0pt{\hskip
13.65979pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\footnotesize{\leavevmode\hbox{\set@color$(\gamma-1)\mathrm{C_{v}}$}}}\llbracket\rho{\textbf{v}}\rrbracket-\hbox
to0.0pt{\hskip
2.52pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\footnotesize{\leavevmode\hbox{\set@color$\mathrm{p}_{\infty}$}}}\llbracket\zeta{\textbf{v}}\rrbracket\big{)}\cdot{\textbf{n}}.$
(46)

Using the definition of $\zeta$ in (eq. 14), the entropy of the mixture (eq.
12) may be reformulated as

$s=\mathrm{C_{v}}\ln\left(\frac{\mathrm{p}+\mathrm{p}_{\infty}}{\rho^{\gamma}}\right)=-\mathrm{C_{v}}\ln\zeta-(\gamma-1)\mathrm{C_{v}}\ln\rho+\mathrm{C_{v}}\ln\big{(}(\gamma-1)\mathrm{C_{v}}\big{)},$
(47)

then we have

$\displaystyle\textbf{h}^{\top}_{ec}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})\llbracket\bm{\vartheta}({\textbf{u}})\rrbracket$
$\displaystyle=\llbracket\gamma\mathrm{C_{v}}-s-\zeta\frac{{\textbf{v}}\cdot{\textbf{v}}}{2}\rrbracket\hat{\rho}\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}}+\left(\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\hat{\rho}+\frac{(\gamma-1)\mathrm{C_{v}}\hskip
0.17343pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\rho$}\hss}{\footnotesize{\hskip-0.17343pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}}{\hskip
0.53125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\zeta$}\hss}{\footnotesize{\hskip-0.53125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}}-\mathrm{p}_{\infty}\right)\llbracket\zeta{\textbf{v}}\rrbracket\cdot{\textbf{n}}$
$\displaystyle\hskip
113.81102pt-\llbracket\zeta\rrbracket\left(\left(\frac{\mathrm{C_{v}}}{\hat{\zeta}}+\frac{{\textbf{v}}^{-}\cdot{\textbf{v}}^{+}}{2}\right)\hat{\rho}+\frac{(\gamma-1)\mathrm{C_{v}}}{\hskip
0.53125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\zeta$}\hss}{\footnotesize{\hskip-0.53125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}}\hskip
0.17343pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\rho$}\hss}{\footnotesize{\hskip-0.17343pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\right)\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}}$
$\displaystyle\underset{(\lx@cref{creftype~refnum}{eqn: chap 5 Leibniz
rule})}{\stackrel{{\scriptstyle\text{(\lx@cref{creftype~refnum}{eqn: chap 5
physical entropy of
conservative})}}}{{=}}}\left(\frac{\llbracket\zeta\rrbracket}{\hat{\zeta}}+(\gamma-1)\frac{\llbracket\rho\rrbracket}{\hat{\rho}}\right)\mathrm{C_{v}}\hat{\rho}\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}}+\llbracket\zeta{\textbf{v}}\rrbracket\left(\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\hat{\rho}+\frac{(\gamma-1)\mathrm{C_{v}}\hskip
0.17343pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\rho$}\hss}{\footnotesize{\hskip-0.17343pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}}{\hskip
0.53125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\zeta$}\hss}{\footnotesize{\hskip-0.53125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}}-\mathrm{p}_{\infty}\right)\cdot{\textbf{n}}$
$\displaystyle\hskip
113.81102pt-\llbracket\zeta\rrbracket\left(\left(\frac{\mathrm{C_{v}}}{\hat{\zeta}}+\frac{{\textbf{v}}^{-}\cdot{\textbf{v}}^{+}}{2}\right)\hat{\rho}+\frac{(\gamma-1)\mathrm{C_{v}}}{\hskip
0.53125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\zeta$}\hss}{\footnotesize{\hskip-0.53125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}}\hskip
0.17343pt\hbox
to0.0pt{\leavevmode\hbox{\set@color$\rho$}\hss}{\footnotesize{\hskip-0.17343pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\right)\hskip
0.125pt\hbox
to0.0pt{\leavevmode\hbox{\set@color${\textbf{v}}$}\hss}{\footnotesize{\hskip-0.125pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}}}\cdot{\textbf{n}}$
$\displaystyle=\big{(}(\gamma-1)\mathrm{C_{v}}\llbracket\rho{\textbf{v}}\rrbracket-\mathrm{p}_{\infty}\llbracket\zeta{\textbf{v}}\rrbracket\big{)}\cdot{\textbf{n}},$

which cancels out with (eq. 46), so the proof is complete. ∎

##  5 An HLLC Riemann solver for the cell interfaces

We now look for two-point numerical fluctuations at interfaces and design an
HLLC approximate Riemann solver for the SG-gamma model (eq. 1a) and analyze
its properties. For the sake of generality, we assume that the entropy
$\eta({\textbf{u}})$ in (eq. 9) is convex, which excludes material interfaces.

###  5.1 One-dimensional Riemann problem

We are interested in approximating solutions to the following Riemann problem
in a given unit direction ${\bf n}$ in $\mathbb{R}^{d}$:

$\partial_{t}{\textbf{u}}+\partial_{x}\textbf{f}({\textbf{u}})\cdot{\textbf{n}}+\textbf{c}_{{\textbf{n}}}({\textbf{u}})\partial_{x}{\textbf{u}}=0,$
(48a)

where $x={\bf x}\cdot{\textbf{n}}$ and
$\textbf{c}_{{\textbf{n}}}({\textbf{u}})=\textbf{c}({\textbf{u}}){\textbf{n}}=\sum^{d}_{i=1}n_{i}\textbf{c}_{i}({\textbf{u}})$,
together with initial data

${\textbf{u}}_{0}(x)=\left\\{\begin{array}[]{ll}{\textbf{u}}_{L},&x<0,\\\
{\textbf{u}}_{R},&x>0.\end{array}\right.$ (48b)

By
$\bm{\mathcal{W}}(\tfrac{x}{t},{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})$
we denote the exact entropy weak solution to (eq. 48) for $t>0$. Following
[31, 26], we integrate (eq. 48a) over the control volume
$[-\tfrac{h}{2},\tfrac{h}{2}]\times[0,\Delta t]$ with $h>0$ and $\Delta t>0$
the space and time steps, respectively. Using (eq. 48b), we obtain

$\int^{\frac{h}{2}}_{-\frac{h}{2}}\bm{\mathcal{W}}\left(\frac{x}{\Delta
t};{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}}\right)\mathrm{d}x-\frac{h}{2}({\textbf{u}}_{L}+{\textbf{u}}_{R})+\Delta
t\left(\textbf{f}({\textbf{u}}_{R})-\textbf{f}({\textbf{u}}_{L})\right)\cdot{\textbf{n}}+\int^{\Delta
t}_{0}\int^{\frac{h}{2}}_{-\frac{h}{2}}\textbf{c}_{{\textbf{n}}}({\textbf{u}})\partial_{x}{\textbf{u}}\mathrm{d}x\mathrm{d}t=0.$
(49)

Note that both $\Gamma$ and $\Pi$ are continuous across shocks and
discontinuous across the intermediate contact wave. Let us introduce the two
last components $\textbf{e}_{\Gamma}=(0,0,0,1,0)^{\top}$ and
$\textbf{e}_{\Pi}=(0,0,0,0,1)^{\top}$ of the canonical basis in
$\mathbb{R}^{n_{eq}}$, and define $u\mathrel{\mathop{\mathchar
58\relax}}={\textbf{v}}\cdot{\textbf{n}}$. By (eq. 49) and (eq. 2) we have

$\displaystyle\int^{\Delta
t}_{0}\int^{\frac{h}{2}}_{-\frac{h}{2}}\textbf{c}_{{\textbf{n}}}({\textbf{u}})\partial_{x}{\textbf{u}}\mathrm{d}x\mathrm{d}t$
$\displaystyle=\frac{h}{2}({\textbf{u}}_{L}+{\textbf{u}}_{R})-\Delta
t\left(\textbf{f}({\textbf{u}}_{R})-\textbf{f}({\textbf{u}}_{L})\right)\cdot{\textbf{n}}-\int^{\frac{h}{2}}_{-\frac{h}{2}}\bm{\mathcal{W}}\left(\frac{x}{\Delta
t};{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}}\right)\mathrm{d}x$
$\displaystyle=\frac{h}{2}(\Gamma_{L}+\Gamma_{R})\textbf{e}_{\Gamma}+\frac{h}{2}(\Pi_{L}+\Pi_{R})\textbf{e}_{\Pi}-\left(\frac{h}{2}-u^{\star}\Delta
t\right)(\Gamma_{R}\textbf{e}_{\Gamma}+\Pi_{R}\textbf{e}_{\Pi})-\left(u^{\star}\Delta
t+\frac{h}{2}\right)(\Gamma_{L}\textbf{e}_{\Gamma}+\Pi_{L}\textbf{e}_{\Pi})$
$\displaystyle=\Delta
tu^{\star}\big{(}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi}\big{)},$

where $u_{L}$, $u^{\star}$, and $u_{R}$ are the normal velocity components in
the left state, the star region, and the right states, respectively. The
integral form for (eq. 1a) thus reads

$\frac{1}{\Delta
t}\int^{\frac{h}{2}}_{-\frac{h}{2}}\bm{\mathcal{W}}\left(\frac{x}{\Delta
t};{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}}\right)\mathrm{d}x-\frac{h}{2\Delta
t}({\textbf{u}}_{L}+{\textbf{u}}_{R})+\big{(}\textbf{f}({\textbf{u}}_{R})-\textbf{f}({\textbf{u}}_{L})\big{)}\cdot{\textbf{n}}+u^{\star}\big{(}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi}\big{)}=0.$
(50)

Likewise integrating (eq. 9) in the direction n over the control volume
$[-\tfrac{h}{2},\tfrac{h}{2}]\times[0,\Delta t]$ gives

$\frac{1}{\Delta
t}\int^{\frac{h}{2}}_{-\frac{h}{2}}\eta\left(\bm{\mathcal{W}}\left(\frac{x}{\Delta
t};{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}}\right)\right)\mathrm{d}x\mathrm{d}t-\frac{h}{2\Delta
t}\big{(}\eta({\textbf{u}}_{L})+\eta({\textbf{u}}_{R})\big{)}+{\textbf{q}}({\textbf{u}}_{R})\cdot{\textbf{n}}-{\textbf{q}}({\textbf{u}}_{L})\cdot{\textbf{n}}\leqslant
0.$ (51)

###  5.2 Three-point schemes and the Godunov method

It will be convenient for the analysis of the HLLC solver to consider one-
dimensional three-point numerical schemes in fluctuation form [47]

${\bf U}_{j}^{n+1}-{\bf U}_{j}^{n}+\frac{\Delta t}{h}\big{(}{\bf D}^{-}({\bf
U}_{j}^{n},{\bf U}_{j+1}^{n},{\bf n})+{\bf D}^{+}({\bf U}_{j-1}^{n},{\bf
U}_{j}^{n},{\bf n})\big{)}=0,$ (52)

where ${\bf D}^{\pm}(\cdot,\cdot,\cdot)$ are assumed to be consistent (eq.
27). Here, ${\bf U}_{j}^{n}$ approximates the cell-averaged solution in the
$j$th cell of size $h$ at time $t^{(n)}=n\Delta t$. In the Godunov method, the
fluctuations are defined by solving exact Riemann problems (eq. 48) centered
at every interface $j+\tfrac{1}{2}$ of coordinate $x_{j+\frac{1}{2}}$, between
cells $j$ and $j+1$ (see Figure 3), with ${\bf U}_{j}^{n}$ and ${\bf
U}_{j+1}^{n}$ as the left and right initial data, respectively. The solution
at time $t^{n+1}$ in the $j$th cell is defined as the cell-average of the
exact solution at $t^{(n+1)}=t^{(n)}+\Delta t$:

$\displaystyle{\textbf{U}}^{n+1}_{j}$
$\displaystyle=\frac{1}{h}\left(\int^{x_{j}}_{x_{j-\frac{1}{2}}}\bm{\mathcal{W}}\left(\frac{x}{\Delta
t};{\textbf{U}}^{n}_{j-1},{\textbf{U}}^{n}_{j},{\textbf{n}}\right)\mathrm{d}x+\int^{x_{j+\frac{1}{2}}}_{x_{j}}\bm{\mathcal{W}}\left(\frac{x}{\Delta
t};{\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\textbf{n}}\right)\mathrm{d}x\right),$
(53) $\displaystyle={\textbf{U}}^{n}_{j}-\frac{\Delta
t}{h}\left(\frac{h}{2\Delta t}{\textbf{U}}^{n}_{j}-\frac{1}{\Delta
t}\int^{x_{j}}_{x_{j-\frac{1}{2}}}\bm{\mathcal{W}}\left(\frac{x}{\Delta
t};{\textbf{U}}^{n}_{j-1},{\textbf{U}}^{n}_{j},{\textbf{n}}\right)\mathrm{d}x+\frac{h}{2\Delta
t}{\textbf{U}}^{n}_{j}-\frac{1}{\Delta
t}\int^{x_{j+\frac{1}{2}}}_{x_{j}}\bm{\mathcal{W}}\left(\frac{x}{\Delta
t};{\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\textbf{n}}\right)\mathrm{d}x\right),$

with $x_{j}=x_{j-\frac{1}{2}}+\tfrac{h}{2}=x_{j+\frac{1}{2}}-\tfrac{h}{2}$
(see Figure 3), and takes the form (eq. 52) with the fluctuations defined as

$\displaystyle{\textbf{D}}^{-}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})$
$\displaystyle=\frac{h}{2\Delta t}{\textbf{u}}_{L}-\frac{1}{\Delta
t}\int^{0}_{-\frac{h}{2}}\bm{\mathcal{W}}\left(\frac{x}{\Delta
t};{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}}\right)\mathrm{d}x$ (54a)
$\displaystyle=\textbf{f}\big{(}\bm{\mathcal{W}}(0;{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})\big{)}\cdot{\textbf{n}}-\textbf{f}({\textbf{u}}_{L})\cdot{\textbf{n}}+\min(u^{\star},0)\big{(}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi}\big{)},$
$\displaystyle{\textbf{D}}^{+}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})\mathrel{\mathop{\mathchar
58\relax}}$ $\displaystyle=\frac{h}{2\Delta t}{\textbf{u}}_{R}-\frac{1}{\Delta
t}\int^{\frac{h}{2}}_{0}\bm{\mathcal{W}}\left(\frac{x}{\Delta
t};{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}}\right)\mathrm{d}x$ (54b)
$\displaystyle=\textbf{f}({\textbf{u}}_{R})\cdot{\textbf{n}}-\textbf{f}\big{(}\bm{\mathcal{W}}(0;{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})\big{)}\cdot{\textbf{n}}+\max(u^{\star},0)\big{(}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi}\big{)}.$

$h$$x$$x_{j-\frac{1}{2}}$$x_{j}$$x_{j+\frac{1}{2}}$${j-1}$${j}$${j+1}$ Figure
3: Notations for the mesh used for the three-point scheme (eq. 52).

Note that ${\textbf{D}}^{\pm}(\cdot,\cdot,\cdot)$ in (eq. 54) satisfy the
following path-conservation property [47]

${\textbf{D}}^{-}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})+{\textbf{D}}^{+}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})=\textbf{f}({\textbf{u}}_{R})\cdot{\textbf{n}}-\textbf{f}({\textbf{u}}_{L})\cdot{\textbf{n}}+u^{\star}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+u^{\star}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Pi},$
(55)

for a path $\bm{\phi}\mathrel{\mathop{\mathchar
58\relax}}[0,1]\times\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm
GM$}}\times\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm
GM$}}\rightarrow\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm GM$}}$ such that

$u^{\star}(\Gamma_{R}-\Gamma_{L})=\int_{0}^{1}u\big{(}\bm{\phi}(s;{\textbf{u}}_{L},{\textbf{u}}_{R})\big{)}\partial_{s}\bm{\phi}(s;{\textbf{u}}_{L},{\textbf{u}}_{R})\cdot\textbf{e}_{\Gamma}ds,\quad
u^{\star}(\Pi_{R}-\Pi_{L})=\int_{0}^{1}u\big{(}\bm{\phi}(s;{\textbf{u}}_{L},{\textbf{u}}_{R})\big{)}\partial_{s}\bm{\phi}(s;{\textbf{u}}_{L},{\textbf{u}}_{R})\cdot\textbf{e}_{\Pi}ds.$

The interface Riemann problems are assumed to be noninteracting through the
imposition of a half CFL condition:

$\frac{\Delta
t}{h}\max_{j\in\mathbb{Z}}|\lambda|_{max}({\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\bf
n})\leqslant\frac{1}{2},$ (56)

where $|\lambda|_{max}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})$ is an
upper bound of the absolute value of the signal speeds in the Riemann problem
(eq. 48).

Finally, invoking a Jensen’s inequality in (eq. 53) for any convex entropy
function in (eq. 11), we obtain the following discrete entropy inequality [31]
consistent with (eq. 9):

$\eta({\textbf{U}}^{n+1}_{j})\leqslant\eta({\textbf{U}}_{j})-\frac{\Delta
t}{h}\left(Q({\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\textbf{n}})-Q({\textbf{U}}^{n}_{j-1},{\textbf{U}}^{n}_{j},{\textbf{n}})\right),$
(57)

with the consistent entropy flux $Q({\textbf{u}}_{L},{\textbf{u}}_{R},{\bf
n})={\textbf{q}}\left(\bm{\mathcal{W}}(0;{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})\right)\cdot{\textbf{n}}$.

###  5.3 HLLC Riemann solver

The HLLC solver [65, 2] is a simple solver [9] with four uniform states
separated by simple discontinuities, see Figure 4:

$\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\left(\frac{x}{t};{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}}\right)=\left\\{\begin{array}[]{ll}{\textbf{u}}_{L},&\frac{x}{t}<s_{L},\\\
{\textbf{u}}^{\star}_{L},&s_{L}<\frac{x}{t}<s^{\star},\\\
{\textbf{u}}^{\star}_{R},&s^{\star}<\frac{x}{t}<s_{R},\\\
{\textbf{u}}_{R},&s_{R}<\frac{x}{t},\end{array}\right.$ (58)

where the wave $s^{\star}$ approximates the speed of the intermediate contact
wave. In contrast to [31], we here require an approximate consistency of the
HLLC solver (eq. 58) with the integral form (eq. 50). Integrating (eq. 58)
over $[-\tfrac{h}{2},\tfrac{h}{2}]$ at time $\Delta t$ gives

$\int^{\frac{h}{2}}_{-\frac{h}{2}}\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\left(\frac{x}{\Delta
t};{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}}\right)\mathrm{d}x=\left(\frac{h}{2}-s_{R}\Delta
t\right){\textbf{u}}_{R}+\Delta
t(s_{R}-s^{\star}){\textbf{u}}^{\star}_{R}+\Delta
t(s^{\star}-s_{L}){\textbf{u}}^{\star}_{L}+\left(s_{L}\Delta
t+\frac{h}{2}\right){\textbf{u}}_{L},$ (59)

and using (eq. 50) with
$\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm HLLC$}}$ and
$s^{\star}$ in place of $\bm{\mathcal{W}}$ and $u^{\star}$, we obtain

$s_{R}({\textbf{u}}_{R}^{\star}-{\textbf{u}}_{R})+s^{\star}({\textbf{u}}^{\star}_{L}-{\textbf{u}}^{\star}_{R})+s_{L}({\textbf{u}}_{L}-{\textbf{u}}_{L}^{\star})+{\bf
f}({\textbf{u}}_{R})\cdot{\textbf{n}}-{\bf
f}({\textbf{u}}_{L})\cdot{\textbf{n}}+s^{\star}\big{(}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi}\big{)}=0.$
(60)

$-\frac{h}{2}$$\frac{h}{2}$$s_{L}$$s_{R}$$s^{*}$$\mathbf{u}_{L}$$\mathbf{u}_{R}$$0$$\Delta
t$$t$$x$$\Gamma_{L},\Pi_{L}$$\Gamma_{R},\Pi_{R}$$s_{L}\Delta t$$s_{R}\Delta
t$$\mathbf{u}_{L}^{*}$$\mathbf{u}_{R}^{*}$ Figure 4: Wave pattern of the HLLC
solver (eq. 58).

The component on mass conservation in the above relation gives

$\rho_{R}(s_{R}-u_{R})+\rho_{L}(u_{L}-s_{L})=\rho_{R}^{\star}(s_{R}-s^{\star})+\rho_{L}^{\star}(s^{\star}-s_{L}),$

which is satisfied by further requiring the half-consistency conditions [3]
which give

$Q_{L}=\rho_{L}(u_{L}-s_{L})=\rho_{L}^{\star}(s^{\star}-s_{L})>0,\quad
Q_{R}=\rho_{R}(s_{R}-u_{R})=\rho_{R}^{\star}(s_{R}-s^{\star})>0,$ (61)

and will be referred to as the mass fluxes [64]. Note that (eq. 61) will be
shown to be important for the proof of entropy stability in section 5.4.1 and
are also usually invoked through the satisfaction of some jump relations
across the $s_{L}$ and $s_{R}$ waves [65, 2] (see below).

From the approximate global consistency relation (eq. 60) we also deduce

$\displaystyle Q_{R}(u_{R}^{\star}-u_{R})+Q_{L}(u_{L}^{\star}-u_{L})$
$\displaystyle=\mathrm{p}_{L}-\mathrm{p}_{R},$ $\displaystyle
Q_{R}({\textbf{v}}_{R}^{\perp}-{\textbf{v}}_{R}^{\star\perp})+Q_{L}({\textbf{v}}_{L}^{\perp}-{\textbf{v}}_{L}^{\star\perp})$
$\displaystyle=0,$ $\displaystyle
Q_{R}(E_{R}^{\star}-E_{R})+Q_{L}(E_{L}^{\star}-E_{L})$
$\displaystyle=\mathrm{p}_{L}u_{L}-\mathrm{p}_{R}u_{R},$ $\displaystyle
s_{R}(\Gamma_{R}-\Gamma_{R}^{\star})+s_{L}(\Gamma_{L}^{\star}-\Gamma_{L})$
$\displaystyle=s^{\star}(\Gamma_{R}-\Gamma_{R}^{\star})+s^{\star}(\Gamma_{L}^{\star}-\Gamma_{L}),$
$\displaystyle s_{R}(\Pi_{R}-\Pi_{R}^{\star})+s_{L}(\Pi_{L}^{\star}-\Pi_{L})$
$\displaystyle=s^{\star}(\Pi_{R}-\Pi_{R}^{\star})+s^{\star}(\Pi_{L}^{\star}-\Pi_{L}),$

where ${\textbf{v}}^{\perp}={\textbf{v}}-u{\textbf{n}}$ denotes the velocity
component perpendicular to n. The second and two latter conditions impose

${\textbf{v}}_{L}^{\star\perp}={\textbf{v}}_{L}^{\perp},\quad{\textbf{v}}_{R}^{\star\perp}={\textbf{v}}_{R}^{\perp},\quad\Gamma^{\star}_{L}=\Gamma_{L},\quad\Gamma^{\star}_{R}=\Gamma_{R},\quad\Pi^{\star}_{L}=\Pi_{L},\quad\Pi^{\star}_{R}=\Pi_{R},$
(62)

meaning that ${\textbf{v}}^{\perp}$, $\Gamma$, and $\Pi$ are continuous across
shocks and may be discontinuous across $s^{\star}$ in agreement with the fact
that ${\textbf{v}}^{\perp}$, $\Gamma$, and $\Pi$ are associated to LD fields.
The remaining unknowns can be computed by imposing the Rankine-Hugoniot
relations across $s_{L}$ and $s_{R}$

$\textbf{f}^{\star}_{L}-\textbf{f}({\textbf{u}}_{L})\cdot{\textbf{n}}=s_{L}({\textbf{u}}^{\star}_{L}-{\textbf{u}}_{L}),\quad\textbf{f}^{\star}_{R}-\textbf{f}({\textbf{u}}_{R})\cdot{\textbf{n}}=s_{R}({\textbf{u}}^{\star}_{R}-{\textbf{u}}_{R}),$
(63)

where ${\bf f}_{X=L,R}^{\star}$ means that we are considering
$\mathrm{p}^{\star}_{X}$ as an unknown instead of evaluating the pressure via
the EOS (eq. 5) and ${\textbf{u}}_{X}^{\star}$. This leads to the following
definition of the velocity in the star regions:

$u^{\star}_{L}=u_{L}-\frac{\mathrm{p}^{\star}_{L}-\mathrm{p}_{L}}{Q_{L}},\quad
u^{\star}_{R}=u_{R}+\frac{\mathrm{p}^{\star}_{R}-\mathrm{p}_{R}}{Q_{R}},$ (64)

and to the following jump relation across the $s^{\star}$ wave:

${\bf f}_{R}^{\star}-{\bf f}_{L}^{\star}=s^{\star}\big{(}{\bf
u}_{R}^{\star}-{\bf
u}_{L}^{\star}-(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}-(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi}\big{)}.$

Further imposing continuity of the velocity and the pressure across the
intermediate wave, $u^{\star}_{L}=u^{\star}_{R}=s^{\star}$ and
$p^{\star}_{L}=p^{\star}_{R}=p^{\star}$, gives the expression for the pressure
and velocity in the star region:

$\mathrm{p}^{\star}=\frac{Q_{L}\mathrm{p}_{R}+Q_{R}\mathrm{p}_{L}+Q_{R}Q_{L}(u_{L}-u_{R})}{Q_{L}+Q_{R}},\quad
s^{\star}=\frac{Q_{L}u_{L}+Q_{R}u_{R}+\mathrm{p}_{L}-\mathrm{p}_{R}}{Q_{L}+Q_{R}}.$
(65)

The other states are then directly obtained:

$\displaystyle\rho^{\star}_{L}=\frac{u_{L}-s_{L}}{s^{\star}-s_{L}}\rho_{L},\;e^{\star}_{L}=e_{L}+(s^{\star}-u_{L})\left(\frac{s^{\star}-u_{L}}{2}-\frac{\mathrm{p}_{L}}{Q_{L}}\right),\;E^{\star}_{L}=e^{\star}_{L}+\frac{1}{2}{\textbf{v}}_{L}^{\star}\cdot{\textbf{v}}_{L}^{\star}=E_{L}+(s^{\star}-u_{L})\left(s^{\star}-\frac{\mathrm{p}_{L}}{Q_{L}}\right),$
(66a)
$\displaystyle\rho^{\star}_{R}=\frac{s_{R}-u_{R}}{s_{R}-s^{\star}}\rho_{R},\;e^{\star}_{R}=e_{R}+(s^{\star}-u_{R})\left(\frac{s^{\star}-u_{R}}{2}+\frac{\mathrm{p}_{R}}{Q_{R}}\right),\;E^{\star}_{R}=e^{\star}_{R}+\frac{1}{2}{\textbf{v}}_{R}^{\star}\cdot{\textbf{v}}_{R}^{\star}=E_{R}+(s^{\star}-u_{R})\left(s^{\star}+\frac{\mathrm{p}_{R}}{Q_{R}}\right).$
(66b)

As for the Godunov method in (eq. 54), we now define the fluctuations through

$\displaystyle{\bf D}^{-}({\bf u}_{L},{\bf u}_{R},{\bf n})$
$\displaystyle=\frac{h}{2\Delta t}{\bf u}_{L}-\frac{1}{\Delta
t}\int_{-\frac{h}{2}}^{0}\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\Big{(}\frac{x}{\Delta t};{\bf u}_{L},{\bf u}_{R},{\bf n}\Big{)}dx,$
(67a) $\displaystyle{\bf D}^{+}({\bf u}_{L},{\bf u}_{R},{\bf n})$
$\displaystyle=\frac{h}{2\Delta t}{\bf u}_{R}-\frac{1}{\Delta
t}\int_{0}^{\frac{h}{2}}\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\Big{(}\frac{x}{\Delta t};{\bf u}_{L},{\bf u}_{R},{\bf n}\Big{)}dx,$
(67b)

and plugging (eq. 58) into (eq. 67) gives

$\displaystyle{\textbf{D}}^{-}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})=\left\\{\begin{array}[]{ll}0,&0<s_{L},\\\
s_{L}({\textbf{u}}^{\star}_{L}-{\textbf{u}}_{L}),&s_{L}<0<s^{\star},\\\
\textbf{f}_{R}^{\star}-\textbf{f}({\textbf{u}}_{L})\cdot{\textbf{n}}+s^{\star}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+s^{\star}(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi},&s^{\star}<0<s_{R},\\\
\textbf{f}({\textbf{u}}_{R})\cdot{\textbf{n}}-\textbf{f}({\textbf{u}}_{L})\cdot{\textbf{n}}+s^{\star}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+s^{\star}(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi},&s_{R}<0,\end{array}\right.$
(68e)
$\displaystyle{\textbf{D}}^{+}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})=\left\\{\begin{array}[]{ll}\textbf{f}({\textbf{u}}_{R})\cdot{\textbf{n}}-\textbf{f}({\textbf{u}}_{L})\cdot{\textbf{n}}+s^{\star}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+s^{\star}(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi},&0<s_{L},\\\
\textbf{f}({\textbf{u}}_{R})\cdot{\textbf{n}}-\textbf{f}_{R}^{\star}+s^{\star}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+s^{\star}(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi},&s_{L}<0<s^{\star},\\\
s_{R}({\textbf{u}}_{R}-{\textbf{u}}_{R}^{\star}),&s^{\star}<0<s_{R},\\\
0,&s_{R}<0.\end{array}\right.$ (68j)

Note that by construction the numerical fluxes (eq. 68) satisfy the relation

${\textbf{D}}^{-}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})+{\textbf{D}}^{+}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})=\textbf{f}({\textbf{u}}_{R})\cdot{\textbf{n}}-\textbf{f}({\textbf{u}}_{L})\cdot{\textbf{n}}+s^{\star}(\Gamma_{R}-\Gamma_{L})\textbf{e}_{\Gamma}+s^{*}(\Pi_{R}-\Pi_{L})\textbf{e}_{\Pi}$

similar to the path-conservation property (eq. 55) [47]. Note that for (eq.
55) to hold, one needs $s^{\star}=u^{\star}$ which would require a root-
finding algorithm to evaluate $u^{\star}$. We here follow another strategy
where we approximate $u^{\star}$ by $s^{\star}$ in (eq. 65), this is justified
by the fact that the nonconservative product is here associated to a LD field.

Finally, we define the updated cell-averaged solution as

${\bf
U}_{j}^{(n+1)}=\frac{1}{h}\int_{x_{j-\frac{1}{2}}}^{x_{j}}\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\Big{(}\frac{x}{\Delta t};{\bf U}_{j-1}^{n},{\bf
U}_{j}^{n},{\textbf{n}}\Big{)}dx+\frac{1}{h}\int_{x_{j}}^{x_{j+\frac{1}{2}}}\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\Big{(}\frac{x}{\Delta t};{\bf U}_{j}^{n},{\bf
U}_{j+1}^{n},{\textbf{n}}\Big{)}dx,$ (69)

so the HLLC solver may be recast as a three-point scheme (eq. 52) with (eq.
68).

###  5.4 Properties of the HLLC solver

In this section, we analyse the properties of the numerical scheme (eq. 52)
using fluxes (eq. 68), where the time step $\Delta t>0$ is assumed to satisfy
the CFL condition

$\frac{\Delta
t}{h}\max_{j\in\mathbb{Z}}\left(\left|s_{R}({\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\bf
n})\right|,\left|s_{L}({\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\bf
n})\right|\right)\leqslant\frac{1}{2},$ (70)

where the wave speeds $s_{L}$ and $s_{R}$ will be defined in section 5.4.4.

####  5.4.1 Discrete entropy inequality

We are here interested in the nonlinear stability of the scheme (eq. 52) and
follow [3] to use the local entropy minimum principles and first prove in
Theorem 5.1 the entropy inequality in integral form (eq. 57) for the HLLC
solution (eq. 58). We then prove in Lemma 5.1 that the local entropy minimum
principles hold for the intermediate states.

######  Theorem 5.1.

Suppose that condition (eq. 70) on the time step holds and that the
intermediate states in the HLLC solver (eq. 58), satisfy
${\textbf{u}}^{\star}_{L},{\textbf{u}}^{\star}_{R}\in\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm
GM$}}$ together with the following local minimum entropy principles

$s({\textbf{u}}^{\star}_{L})\geqslant s({\textbf{u}}_{L}),\quad
s({\textbf{u}}^{\star}_{R})\geqslant s({\textbf{u}}_{R}),$ (71)

for the specific entropy (eq. 12). Then, the three-point scheme (eq. 52)
satisfies an entropy inequality (eq. 57) with the consistent numerical flux

$Q({\textbf{U}}^{n}_{j-1},{\textbf{U}}^{n}_{j},{\textbf{n}})={\textbf{q}}({\textbf{U}}^{n}_{j})\cdot{\textbf{n}}+\frac{1}{\Delta
t}\int^{x_{j}}_{x_{j-\frac{1}{2}}}\eta\left(\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\left(\frac{x}{\Delta
t};{\textbf{U}}^{n}_{j-1},{\textbf{U}}^{n}_{j},{\textbf{n}}\right)\right)\mathrm{d}x-\frac{h}{2\Delta
t}\eta\left({\textbf{U}}^{n}_{j}\right).$ (72)

###### Proof.

We first prove the entropy inequality in integral form (eq. 51) using (eq.
59), therefore, we have

$\displaystyle\int^{\frac{h}{2}}_{-\frac{h}{2}}\eta\left(\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\left(\frac{x}{\Delta
t};{\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}}\right)\right)\mathrm{d}x=$
$\displaystyle\left(s_{L}\Delta
t+\frac{h}{2}\right)\eta({\textbf{u}}_{L})+(s^{\star}-s_{L})\Delta
t\eta({\textbf{u}}^{\star}_{L})+(s_{R}-s^{\star})\Delta
t\eta({\textbf{u}}^{\star}_{R})+\left(\frac{h}{2}-s_{R}\Delta
t\right)\eta({\textbf{u}}_{R}),$
$\displaystyle\overset{(\lx@cref{creftype~refnum}{eqn: chap 5 entropy
pair})}{=}$
$\displaystyle\frac{h}{2}\big{(}\eta({\textbf{u}}_{L})+\eta({\textbf{u}}_{R})\big{)}-\Delta
t\big{(}s_{L}\rho_{L}s({\textbf{u}}_{L})+(s^{\star}-s_{L})\rho^{\star}_{L}s({\textbf{u}}^{\star}_{L})\big{)}$
$\displaystyle-\Delta
t\big{(}(s_{R}-s^{\star})\rho^{\star}_{R}s({\textbf{u}}^{\star}_{R})-s_{R}\rho_{R}s({\textbf{u}}_{R})\big{)},$
$\displaystyle\overset{(\lx@cref{creftype~refnum}{eqn: chap 5 intermediate
entropy condition})}{\leqslant}$
$\displaystyle\frac{h}{2}\big{(}\eta({\textbf{u}}_{L})+\eta({\textbf{u}}_{R})\big{)}-\Delta
ts({\textbf{u}}_{L})\left(s_{L}\rho_{L}+(s^{\star}-s_{L})\rho^{\star}_{L}\right)-\Delta
t\big{(}(s_{R}-s^{\star})\rho^{\star}_{R}-s_{R}\rho_{R}\big{)}s({\textbf{u}}_{R}),$
$\displaystyle=$
$\displaystyle\frac{h}{2}\big{(}\eta({\textbf{u}}_{L})+\eta({\textbf{u}}_{R})\big{)}-\Delta
t\big{(}\rho_{R}{\textbf{v}}_{R}s({\textbf{u}}_{R})-\rho_{L}{\textbf{v}}_{L}s({\textbf{u}}_{L})\big{)}\cdot{\textbf{n}},$

where we have used the half-consistency conditions (eq. 61). As a consequence,
setting ${\textbf{u}}_{L}={\textbf{U}}^{n}_{j}$ and
${\textbf{u}}_{R}={\textbf{U}}^{n}_{j+1}$, the numerical flux (eq. 72)
satisfies

$Q({\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\textbf{n}})\leqslant{\textbf{q}}({\textbf{U}}^{n}_{j})\cdot{\textbf{n}}-\frac{1}{h}\int^{x_{j+\frac{1}{2}}}_{x_{j}}\eta\left(\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\left(\frac{x}{\Delta
t};{\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\textbf{n}}\right)\right)\mathrm{d}x+\frac{h}{2\Delta
t}\eta({\textbf{U}}^{n}_{j}),$ (73)

and using (eq. 69) we have the following relation through Jensen’s inequality
for the convex entropy function (eq. 11)

$\displaystyle\eta({\textbf{U}}^{n+1}_{j})$
$\displaystyle\leqslant\frac{1}{h}\int^{x_{j}}_{x_{j-\frac{1}{2}}}\eta\left(\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\left(\frac{x}{\Delta
t};{\textbf{U}}^{n}_{j-1},{\textbf{U}}^{n}_{j},{\textbf{n}}\right)\right)\mathrm{d}x+\frac{1}{h}\int^{x_{j+\frac{1}{2}}}_{x_{j}}\eta\left(\bm{\mathcal{W}}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}\left(\frac{x}{\Delta
t};{\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\textbf{n}}\right)\right)\mathrm{d}x,$
$\displaystyle\underset{(\lx@cref{creftype~refnum}{eqn: chap 5 HLLC discrete
entropy flux inequality})}{\overset{(\lx@cref{creftype~refnum}{eqn: chap 5
entropy flux for entropy ineq})}{\leqslant}}\frac{\Delta
t}{h}\left(Q({\textbf{U}}^{n}_{j-1},{\textbf{U}}^{n}_{j},{\textbf{n}})-{\textbf{q}}({\textbf{U}}^{n}_{j})\cdot{\textbf{n}}+\frac{h}{2\Delta
t}\eta({\textbf{U}}^{n}_{j}))\right)+\frac{\Delta
t}{h}\left(-Q({\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\textbf{n}})+{\textbf{q}}({\textbf{U}}^{n}_{j})\cdot{\textbf{n}}+\frac{h}{2\Delta
t}\eta({\textbf{U}}^{n}_{j})\right),$
$\displaystyle=\eta({\textbf{U}}^{n}_{j})-\frac{\Delta
t}{h}\left(Q({\textbf{U}}^{n}_{j},{\textbf{U}}^{n}_{j+1},{\textbf{n}})-Q({\textbf{U}}^{n}_{j-1},{\textbf{U}}^{n}_{j},{\textbf{n}})\right).$

∎

######  Lemma 5.1.

There exist wave speed estimates $s_{L}$ and $s_{R}$ large enough that: (i)
bound the minimum and maximum wave speeds in the exact entropy weak solution
of the Riemann problem (eq. 48), (ii) satisfy the interlacing condition
$s_{L}<s^{\star}<s_{R}$, (iii) ensure that the local minimum entropy
principles (eq. 71) hold.

###### Proof.

These results are consequences of [3, Prop. 3.2] and we now show that the
required assumptions on the half domains separated by $\tfrac{x}{t}=s^{\star}$
hold, see [3, after Prop. 3.2]. This is justified here because the EOS (eq. 5)
does not change accross the extreme waves $s_{L}$ and $s_{R}$, hence
$\tfrac{x}{t}=s^{\star}$ separates domains with one unique equivalent pure
phase and its associated entropy (eq. 12) which is strictly convex and
satisfies the Gibbs principle (eq. 13). The required conditions across the
extreme waves are: first, the half-consistency relations hold through (eq.
61); then, the following quantities must be invariant across the $s_{L}$ and
$s_{R}$ waves [3, § 4.2]:

$\mathrm{p}^{\star}_{X}+\frac{Q^{2}_{X}}{\rho^{\star}_{X}}=\mathrm{p}_{X}+\frac{Q^{2}_{X}}{\rho_{X}},\quad
e^{\star}_{X}-\frac{(\mathrm{p}^{\star}_{X})^{2}}{2Q^{2}_{X}}=e_{X}-\frac{\mathrm{p}^{2}_{X}}{2Q^{2}_{X}},\quad
X=L,R.$

The first relation is a direct consequence of (eq. 64), with
$u_{L}^{\star}=u_{R}^{\star}=s^{\star}$:

$\mathrm{p}_{L}^{\star}=\mathrm{p}_{L}+Q_{L}(u_{L}-s^{\star}),\quad\mathrm{p}_{R}^{\star}=\mathrm{p}_{R}+Q_{R}(s^{\star}-u_{R}),$

and (eq. 61). For the second relation, we inject the above relations in the
expression of $e_{X=L,R}^{\star}$ in (eq. 66) to get

$e_{X}^{\star}-\mathrm{p}_{X}=\frac{\mathrm{p}_{X}-\mathrm{p}_{X}^{\star}}{Q_{X}}\left(\frac{\mathrm{p}_{X}-\mathrm{p}_{X}^{\star}}{2Q_{X}}-\frac{\mathrm{p}_{X}-}{Q_{X}}\right)=-\frac{(\mathrm{p}_{X}^{\star})^{2}-\mathrm{p}_{X}^{2}}{2Q_{X}},\quad
X=L,R,$

which concludes the proof. ∎

We finally link the discrete entropy inequality (eq. 51) to the entropy stable
character of the numerical fluctuations [10], so that the HLLC fluxes can be
used at the interfaces in the DG scheme to prove the semi-discrete entropy
inequality established in Theorem 3.1.

######  Corollary 5.1.

Let a three-point scheme of the form (eq. 52) to discretize (eq. 1a). Then,
the entropy inequality (eq. 57) implies the entropy stability of the numerical
fluxes in the sense of (eq. 34).

###### Proof.

The proof relies on similar arguments as the ones used in [6, Lemma 2.8] in
the conservative setting. Let
${\textbf{U}}^{n}_{j-1}={\textbf{U}}^{n}_{j}={\textbf{u}}_{L}$ and
${\textbf{U}}^{n}_{j+1}={\textbf{u}}_{R}$, then from (eq. 52), we obtain
${\textbf{U}}^{n+1}_{j}={\textbf{u}}_{L}-\tfrac{\Delta
t}{h}\textbf{D}^{-}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})$ and (eq.
57) gives

$\eta\left({\textbf{u}}_{L}-\frac{\Delta
t}{h}\textbf{D}^{-}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})\right)\leqslant\eta({\textbf{u}}_{L})-\frac{\Delta
t}{h}\left(Q({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})-{\textbf{q}}({\textbf{u}}_{L})\cdot{\textbf{n}}\right).$

Likewise, using ${\textbf{U}}^{n}_{j-1}={\textbf{u}}_{L}$ and
${\textbf{U}}^{n}_{j}={\textbf{U}}^{n}_{j+1}={\textbf{u}}_{R}$, we get
${\textbf{U}}^{n+1}_{j}={\textbf{u}}_{R}-\tfrac{\Delta
t}{h}\textbf{D}^{+}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})$ and

$\eta\left({\textbf{u}}_{R}-\frac{\Delta
t}{h}\textbf{D}^{+}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})\right)\leqslant\eta({\textbf{u}}_{R})-\frac{\Delta
t}{h}\left({\textbf{q}}({\textbf{u}}_{R})\cdot{\textbf{n}}-Q({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})\right).$

Summing both equations and letting $\Delta t\to 0^{+}$ with fixed $h$, we have
up to ${\cal O}(\Delta t)$:

$\eta({\textbf{u}}_{L})-\frac{\Delta
t}{h}\bm{\vartheta}({\textbf{u}}_{L})\cdot\textbf{D}^{-}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})+\eta({\textbf{u}}_{R})-\frac{\Delta
t}{h}\bm{\vartheta}({\textbf{u}}_{R})\cdot\textbf{D}^{+}({\textbf{u}}_{L},{\textbf{u}}_{R},{\textbf{n}})\leqslant\eta({\textbf{u}}_{L})+\eta({\textbf{u}}_{R})-\frac{\Delta
t}{h}\left({\textbf{q}}({\textbf{u}}_{R})-{\textbf{q}}({\textbf{u}}_{L})\right)\cdot{\textbf{n}},$

and simplifying terms gives (eq. 34). ∎

Finally, as an immediate consequence of the local minimum entropy principles
(eq. 71) and the definition of the updated solution (eq. 69), the HLLC solver
also satisfies a discrete minimum principle on the specific physical entropy
$s({\bf u})$:

$s({\bf U}_{j}^{n+1})\geqslant\min\big{(}s({\bf U}_{j-1}^{n}),s({\bf
U}_{j}^{n}),s({\bf U}_{j+1}^{n})\big{)}.$ (74)

####  5.4.2 Preservation of material interfaces and pure phases

We, first, prove that the three-point scheme (eq. 52), with numerical fluxes
(eq. 68), preserves material interfaces [1] by following the same analysis as
in section 4.1. Let us assume that the left and right states satisfy
${\textbf{v}}_{L}={\textbf{v}}_{R}={\textbf{v}}$ and
$\mathrm{p}_{L}=\mathrm{p}_{R}=\mathrm{p}$, then (eq. 65) gives
$\mathrm{p}^{\star}=\mathrm{p}$ and $s^{\star}=u$, respectively. As a result,
the discrete requirements
$\mathrm{d}(\rho{\textbf{v}})={\textbf{v}}\mathrm{d}\rho$ and $\mathrm{d}\rho
E=(\tfrac{1}{2}{\textbf{v}}\cdot{\textbf{v}})\mathrm{d}\rho+\mathrm{p}\mathrm{d}\Gamma+\mathrm{d}\Pi$
applied to the three-point scheme (eq. 52) impose

$D^{\pm}_{\rho u}=uD^{\pm}_{\rho},\quad D^{\pm}_{\rho
E}=\frac{{\textbf{v}}\cdot{\textbf{v}}}{2}D^{\pm}_{\rho}+\mathrm{p}D^{\pm}_{\Gamma}+D^{\pm}_{\Pi},$
(75)

and are obviously satisfied. Now assume that $\Gamma_{L}=\Gamma_{R}$ and
$\Pi_{L}=\Pi_{R}$, then $D^{\pm}_{\Gamma}=0$ and $D^{\pm}_{\Pi}=0$ and the
fluxes reduce to the conservative HLLC solver for the Euler equations so pure
phases are preserved.

####  5.4.3 Positivity of the solution

Since the updated solution (eq. 69) is the cell-average of the superposition
of approximate Riemann solutions from the HLLC solver (eq. 58), and assuming
that the left and right states are positive, it is sufficient to prove
positivity of the intermediate states in the Riemann solution [22]. Note that
the intermediate states should be evaluated from the intermediate fluxes [2],
${\bf f}_{X=L,R}^{\star}$ in (eq. 63) since we are imposing the pressure via
(eq. 65), and not evaluating it from the EOS (eq. 5), see section 5.3.

The proof for positivity of density in the star region can be directly stated
following its definition in (eq. 66) since $s_{L}<u_{L}$, $u_{R}<s_{R}$ from
the wave estimates in section 5.4.4, and the interlacing property
$s_{L}<s^{\star}<s_{R}$ in Lemma 5.1. According to (eq. 7) hyperbolicity of
(eq. 1a) and positivity of the solution in the star region require satisfying
$\rho^{\star}_{X}e^{\star}_{X}>\mathrm{p}_{\infty_{X}}$, $X=L,R$, which gives
for the left intermediate state:

$\displaystyle\rho^{\star}_{L}\left(E^{\star}_{L}-\frac{(s^{\star})^{2}}{2}\right)>\mathrm{p}_{\infty_{L}}\overset{(\lx@cref{creftype~refnum}{eqn:
chap 5 HLLC star variables})}{\underset{(\lx@cref{creftype~refnum}{eqn: mass
flux HLLC})}{\Leftrightarrow}}$
$\displaystyle\rho^{\star}_{L}\left(E_{L}+s^{\star}(s^{\star}-u_{L})-(s^{\star}-u_{L})\frac{\mathrm{p}_{L}}{\rho_{L}(s_{L}-u_{L})}-\frac{(s^{\star})^{2}}{2}\right)>\mathrm{p}_{\infty_{L}}$
$\displaystyle\Leftrightarrow$
$\displaystyle\rho^{\star}_{L}\left(e_{L}+\frac{(u_{L}-s^{\star})^{2}}{2}-(s^{\star}-u_{L})\frac{\mathrm{p}_{L}}{\rho_{L}(u_{L}-s_{L})}\right)>\mathrm{p}_{\infty_{L}}$
$\displaystyle\overset{(\lx@cref{creftype~refnum}{eqn: chap 5 mixture
EOS})}{\underset{(\lx@cref{creftype~refnum}{eqn: mass flux
HLLC})}{\Leftrightarrow}}$
$\displaystyle\left(\frac{u_{L}-s_{L}}{s^{\star}-s_{L}}\right)\left(\frac{\mathrm{p}_{L}+\gamma_{L}\mathrm{p}_{\infty_{L}}}{\gamma_{L}-1}\right)+\left(\frac{u_{L}-s_{L}}{s^{\star}-s_{L}}\right)\frac{\rho_{L}(u_{L}-s^{\star})^{2}}{2}-\left(\frac{s^{\star}-u_{L}}{s^{\star}-s_{L}}\right)\mathrm{p}_{L}-\mathrm{p}_{\infty_{L}}>0$
$\displaystyle\overset{s^{\star}>s_{L}}{\Leftrightarrow}$
$\displaystyle(u_{L}-s_{L})\left(\frac{\mathrm{p}_{L}+\gamma_{L}\mathrm{p}_{\infty_{L}}}{\gamma_{L}-1}\right)+(u_{L}-s_{L})\frac{\rho_{L}}{2}\sigma^{2}+\sigma\mathrm{p}_{L}-(s^{\star}-s_{L})\mathrm{p}_{\infty_{L}}>0$
$\displaystyle\Leftrightarrow$
$\displaystyle(u_{L}-s_{L})\frac{\rho_{L}}{2}\sigma^{2}+(\mathrm{p}_{L}+\mathrm{p}_{\infty_{L}})\sigma+(u_{L}-s_{L})\left(\frac{\mathrm{p}_{L}+\mathrm{p}_{\infty_{L}}}{\gamma_{L}-1}\right)>0,$

where $\sigma=u_{L}-s^{\star}$. The sign for this inequality holds for all
$\sigma\in\mathbb{R}$ if the discriminant $\mathcal{D}$ of the above quadratic
equation is negative:

$\mathcal{D}=(\mathrm{p}_{L}+\mathrm{p}_{\infty_{L}})^{2}-2\rho_{L}(u_{L}-s_{L})^{2}\left(\frac{\mathrm{p}_{L}+\mathrm{p}_{\infty_{L}}}{\gamma_{L}-1}\right)$
(76)

Using $\gamma_{L}(\mathrm{p}_{L}+\mathrm{p}_{\infty_{L}})=\rho_{L}c_{L}^{2}$,
$\mathcal{D}<0$ implies $s_{L}<u_{L}-\sqrt{(\gamma_{L}-1)/2\gamma_{L}}c_{L}$
which is satisfied by the wave speed estimates in section 5.4.4 since
$\tfrac{\gamma_{L}-1}{2\gamma_{L}}<1$. A similar result holds for the right
intermediate state. Note that (76) is the same as in [2] with
$\mathrm{p}_{X}+\mathrm{p}_{\infty_{X}}$ instead of $\mathrm{p}_{X}$.

Finally, let us prove positivity of $\Gamma$ and $\Pi$ through a discrete
maximum principle. The discrete equation for $X\in\\{\Gamma,\Pi\\}$ in (52)
with fluctuations (eq. 68) reads

$\displaystyle X^{n+1}_{j}$ $\displaystyle=X^{n}_{j}-\frac{\Delta
t}{h}\left(\min(s^{*}_{j+\frac{1}{2}},0)\left(X^{n}_{j+1}-X^{n}_{j}\right)+\max(s^{*}_{j-\frac{1}{2}},0)\left(X^{n}_{j}-X^{n}_{j-1}\right)\right)$
$\displaystyle=\left(1-\frac{\Delta
t}{h}\left(\max(s^{*}_{j-\frac{1}{2}},0)-\min(s^{*}_{j+\frac{1}{2}},0)\right)\right)X^{n}_{j}-\frac{\Delta
t}{h}\min(s^{*}_{j+\frac{1}{2}},0)X^{n}_{j+1}+\frac{\Delta
t}{h}\max(s^{*}_{j-\frac{1}{2}},0)X^{n}_{j-1}$ (77)

which shows that $X_{j}^{n+1}$ is a convex combination of the $X_{i=j\pm
1,j}^{n}$ under the CFL condition (eq. 70).

####  5.4.4 Wave speed estimates

The fan of waves of the HLLC solver must contain the fan of waves of the exact
Riemann problem (eq. 48). This is in particular required to ensure the local
entropy minimum principles (eq. 71), see [3, Prop. 3.2]. Direct wave speed
estimates have been proposed that do not require one to solve the exact
Riemann problem [62, 6]. Note that the usual estimate
$S_{R}=-S_{L}=\max(|u_{L}|+c_{L},|u_{R}|+c_{R})$ may be wrong due to the Lax
entropy condition [43] across a shock. On the other hand, the time steps can
be affected by overestimating wave speed estimates through (eq. 70). We
propose the following wave speeds estimates

$s_{L}=u_{L}-\tilde{c}_{L},\quad s_{R}=u_{R}+\tilde{c}_{R},$

where $u_{X}={\textbf{v}}_{X}\cdot{\textbf{n}}$, $X=L,R$, and

$\text{if}\;\mathrm{p}_{R}\geqslant\mathrm{p}_{L}\mathrel{\mathop{\mathchar
58\relax}}\;\left\\{\begin{array}[]{l}\tilde{c}_{L}=c_{L}+\frac{\gamma+1}{2}\max\left(\frac{\mathrm{p}_{R}-\mathrm{p}_{L}}{\rho_{R}c_{R}}+u_{L}-u_{R},0\right),\\\
\tilde{c}_{R}=c_{R}+\frac{\gamma+1}{2}\max\left(\frac{\mathrm{p}_{L}-\mathrm{p}_{R}}{\rho_{L}\tilde{c}_{L}}+u_{L}-u_{R},0\right),\end{array}\right.\quad\text{else}\mathrel{\mathop{\mathchar
58\relax}}\;\left\\{\begin{array}[]{l}\tilde{c}_{R}=c_{R}+\frac{\gamma+1}{2}\max\left(\frac{\mathrm{p}_{L}-\mathrm{p}_{R}}{\rho_{L}c_{L}}+u_{L}-u_{R},0\right),\\\
\tilde{c}_{L}=c_{L}+\frac{\gamma+1}{2}\max\left(\frac{\mathrm{p}_{R}-\mathrm{p}_{L}}{\rho_{R}\tilde{c}_{R}}+u_{L}-u_{R},0\right),\end{array}\right.$

and $\gamma=\max(\gamma_{L},\gamma_{R})$,
$c_{X}=\sqrt{\gamma(\mathrm{p}_{X}+\mathrm{p}_{\infty_{X}})/\rho_{X}}$ for
$X=L,R$. These estimates will bound the wave speeds in the exact Riemann
solution in the case of pure phases [6]. Moreover, the above definition of
$\gamma$ enables one to bound the signal speeds in the case of polytropic
gases, $\mathrm{p}_{\infty_{L}}=\mathrm{p}_{\infty_{R}}=0$ [45, 53]. Although,
it is difficult to guarantee such properties in the general case when
$\Gamma_{L}\neq\Gamma_{R}$ and $\Pi_{L}\neq\Pi_{r}$, these estimates proved to
be robust in the present numerical experiments.

##  6 Properties of the proposed DGSEM

We recall here the main properties of the DGSEM proposed in this work for the
discretization of the SG-gamma model (eq. 1) with a stiffened gas EOS (eq. 5).

###  6.1 Semi-discrete scheme

The semi-discrete scheme (eq. 29) with the EC fluxes (eq. 45) in the volume
integral and the HLLC flux (eq. 68) at interfaces satisfies a semi-discrete
entropy inequality, see Theorem 3.1. In contrast, the scheme with the CP
fluxes (eq. 44) in the volume integral, along with the HLLC solver at the
interfaces, preserves uniform pressure and velocity profiles across material
interfaces. Both are high-order accurate and preserve uniform states.

###  6.2 Fully discrete scheme

We again consider the HLLC solver (eq. 68) at interfaces and both CP and EC
fluctuations in the volume integrals. We restrict ourselves to the use of a
one-step first-order explicit time discretization. High-order time integration
will be achieved by using a strong-stability preserving explicit Runge-Kutta
method from [57] that is a convex combination of forward Euler steps and thus
keeps the properties of the first-order in time scheme under some condition on
the time step. The fully discrete DGSEM reads

$\omega_{i}\omega_{j}J_{\kappa}^{ij}\frac{{\bf U}_{\kappa}^{ij,n+1}-{\bf
U}_{\kappa}^{ij,n}}{\Delta t^{(n)}}+{\bf R}_{\kappa}^{ij}({\bf
u}_{h}^{(n)})=0\quad\forall\kappa\in\Omega_{h},\;0\leqslant i,j\leqslant
p,\;n\geqslant 0,$ (78)

where $\Delta t^{(n)}=t^{(n+1)}-t^{(n)}>0$ is the time step, ${\bf
U}_{\kappa}^{ij,n}={\bf U}_{\kappa}^{ij}(t^{(n)})$, ${\bf u}_{h}^{(n)}={\bf
u}_{h}(\cdot,t^{(n)})$, and the vector of space residuals ${\bf
R}_{\kappa}^{ij}(\cdot)$ is defined by (eq. 29). Theorem 6.1 below summarizes
the properties of the above scheme with fluctuations (eq. 28) and (eq. 26).
These properties are independent of the volume fluctuations and therefore hold
for both EC and CP fluctuations for the SG-gamma model (eq. 2), though the
minimum entropy principle holds when excluding material interfaces. Likewise,
any other interface fluctuation that satisfy properties (i) to (iv) below can
be used in place of the HLLC solver. The derivation of the CFL condition will
rely on the work in [7, Lemma 3.4], that proves that there exist pseudo-
equilibrium states ${\textbf{u}}_{\kappa}^{\star,n}$ in
$\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm GM$}}$ and finite wave speed
estimates $\lambda_{\kappa}^{\star,n}>0$ such that

$\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\omega_{k}J_{e}^{k}{\bf
h}^{\scalebox{0.8}{$\scriptscriptstyle\rm Rus$}}\big{(}{\bf
u}_{\kappa}^{\star,n},{\bf u}_{h}^{-}({\bf x}_{e}^{k},t^{(n)}),{\bf
n}_{e}^{k}\big{)}=0,\quad\lambda_{\kappa}^{\star,n}\geqslant\max_{0\leqslant
k\leqslant p,\;e\in\partial\kappa}\big{(}|\lambda|_{max}\Big{(}{\bf
u}_{\kappa}^{\star,n},{\bf u}_{h}^{-}({\bf
x}_{e}^{k},t^{(n)}),{\textbf{n}}_{e}^{k}\big{)}\Big{)},$ (79)

where ${\bf h}^{\scalebox{0.8}{$\scriptscriptstyle\rm Rus$}}({\bf u}^{-},{\bf
u}^{+},{\textbf{n}})=\tfrac{1}{2}\big{(}{\bf f}({\bf u}^{-})+{\bf f}({\bf
u}^{+})\big{)}\cdot{\textbf{n}}-\tfrac{\lambda_{\kappa}^{\star,n}}{2}({\bf
u}^{+}-{\bf u}^{-})$ is the Rusanov flux and $\lambda_{\kappa}^{\star,n}$
bounds the fans of waves in the exact Riemann problems (eq. 48) with
${\textbf{n}}={\textbf{n}}_{e}^{k}$ and left and right states ${\bf
u}_{\kappa}^{\star,n}$ and ${\bf u}_{h}^{-}({\bf x}_{e}^{k},t^{(n)})$, see
(eq. 80).

######  Theorem 6.1.

Let us consider the numerical scheme (eq. 78) with consistent fluctuations
(eq. 28) of the form (eq. 35) in the volume integrals and consistent interface
fluctuations (eq. 26) such that the associated three-point scheme (eq. 52):
(i) is robust: ${\bf U}_{j\in\mathbb{Z}}^{n\geqslant
0}\in\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm GM$}}$; (ii) satisfies a
discrete entropy inequality (eq. 57) with consistent numerical flux; (iii)
satisfies a discrete minimum principle (eq. 74) on the specific physical
entropy $s({\textbf{u}})$; (iv) satisfies discrete maximum principles (eq. 77)
on $\Gamma$ and $\Pi$. Then, the updated cell-averaged solution $\langle{\bf
u}_{h}^{(n+1)}\rangle_{\kappa}$ of the numerical scheme is a convex
combination of DOFs at time $t^{(n)}$ and updates of three-point schemes under
the following conditions on the time step:

$\displaystyle\Delta t^{(n)}\max_{e\in{\cal E}_{h}}\max_{0\leqslant k\leqslant
p}\frac{\omega_{k}J_{e}^{k}}{\min(\tilde{\omega}_{\kappa^{\pm}}({\textbf{x}}_{e}^{k})J_{\kappa^{\pm}}({\textbf{x}}_{e}^{k}))}\max\big{(}\big{|}s^{n}_{L_{j+\frac{1}{2}}}\big{|},\big{|}s^{n}_{R_{j-\frac{1}{2}}}\big{|},\lambda_{\kappa^{-}}^{\star,n},\lambda_{\kappa^{+}}^{\star,n}\big{)}$
$\displaystyle\leqslant\frac{1}{2},$ (80a) $\displaystyle\Delta
t^{(n)}\max_{\kappa\in\Omega_{h}}\max_{0\leqslant i,j\leqslant
p}\sum_{k=0}^{p}\frac{\omega_{k}\Big{(}\omega_{j}D_{ki}{\textbf{v}}_{\kappa}^{kj,n}\cdot{\textbf{n}}_{(i,k)j}+\omega_{i}D_{kj}{\textbf{v}}_{\kappa}^{ik,n}\cdot{\textbf{n}}_{i(j,k)}-\sum\limits_{e\in\partial\kappa}\phi_{\kappa}^{ij}({\textbf{x}}_{e}^{k})J_{e}^{k}\min\big{(}s^{\star}({\textbf{x}}_{e}^{k},t^{(n)}),0\big{)}\Big{)}}{\omega_{i}\omega_{j}J_{\kappa}^{ij}}$
$\displaystyle\leqslant 1,$ (80b)

where $\lambda_{\kappa}^{\star,n}$ is defined in (eq. 79),
$\tilde{\omega}_{\kappa}({\textbf{x}}_{e}^{k})=\tfrac{\omega_{i}\omega_{j}}{d}$
if ${\textbf{x}}_{e}^{k}={\textbf{x}}_{\kappa}^{ij}$ is a vertex of $\kappa$,
$\tilde{\omega}_{\kappa}({\textbf{x}}_{e}^{k})=\tfrac{\omega_{i}\omega_{j}}{d-1}$
if ${\textbf{x}}_{e}^{k}={\textbf{x}}_{\kappa}^{ij}$ is on an edge, and
$\tilde{\omega}_{\kappa}({\textbf{x}}_{e}^{k})=\omega_{i}\omega_{j}$ if
${\textbf{x}}_{e}^{k}={\textbf{x}}_{\kappa}^{ij}$ is on a face. Assuming that
the DOFs ${\textbf{U}}^{ij,n}_{\kappa}$ are in
$\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm GM$}}$ for all
$\kappa\in\Omega_{h}$ and $0\leqslant i,j\leqslant p$, the scheme (eq. 78)
guarantees positivity of the cell-averaged solution:

$\langle{\textbf{u}}^{n+1}_{h}\rangle_{\kappa}\in\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm
GM$}}\quad\forall\kappa\in\Omega_{h},$

together with a minimum principle on the specific entropy when excluding
material interfaces

$s(\langle{\bf u}_{h}^{(n+1)}\rangle_{\kappa})\geqslant\min\big{\\{}s({\bf
U}_{\kappa}^{ij,n})\mathrel{\mathop{\mathchar 58\relax}}\;0\leqslant
i,j\leqslant p\big{\\}}\cup\big{\\{}s\big{(}{\bf u}_{h}^{+}({\bf
x}_{e}^{k},t^{(n)})\big{)}\mathrel{\mathop{\mathchar
58\relax}}\;e\in\partial\kappa,0\leqslant k\leqslant
p\big{\\}}\quad\forall\kappa\in\Omega_{h},$

and maximum principles on the EOS parameters $Y$ in $\\{\Gamma,\Pi\\}$:

$\min{\cal S}_{\kappa}(Y_{h}^{(n)})\leqslant\langle
Y_{h}^{(n+1)}\rangle_{\kappa}\leqslant\max{\cal
S}_{\kappa}(Y_{h}^{(n)})\;\forall\kappa\in\Omega_{h},\quad{\cal
S}_{\kappa}(Y_{h}^{(n)})\mathrel{\mathop{\mathchar
58\relax}}=\\{Y_{\kappa}^{ij}\mathrel{\mathop{\mathchar 58\relax}}0\leqslant
i,j\leqslant p\\}\cup\big{\\{}Y_{h}^{+}({\bf
x}_{e}^{k},t)\mathrel{\mathop{\mathchar
58\relax}}e\in\partial\kappa,0\leqslant k\leqslant p\big{\\}}.$

###### Proof.

From (eq. 36) and (eq. 78), the cell-averaged discrete scheme for the $d+2$
first components of (eq. 1a), $Y$ in $\\{\rho,\rho{\textbf{v}},\rho E\\}$,
reads

$\langle Y_{h}^{(n+1)}\rangle_{\kappa}=\langle
Y_{h}^{(n)}\rangle_{\kappa}-\frac{\Delta
t^{(n)}}{|\kappa|}\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\omega_{k}J_{e}^{k}{\bf
h}_{Y}^{\scalebox{0.8}{$\scriptscriptstyle\rm HLLC$}}\big{(}{\bf
u}_{h}^{-}({\bf x}_{e}^{k},t^{(n)}),{\bf u}_{h}^{+}({\bf
x}_{e}^{k},t^{(n)}),{\bf n}_{e}^{k}\big{)},$

where ${\bf h}_{Y}^{\scalebox{0.8}{$\scriptscriptstyle\rm
HLLC$}}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})={\textbf{D}}_{Y}^{-}({\textbf{u}}^{-},{\textbf{u}}^{+},{\textbf{n}})+{\bf
f}_{Y}({\textbf{u}}^{-})\cdot{\textbf{n}}$ is the numerical flux in
conservation form associated to the HLLC fluctuations (eq. 68). We first
decompose $\langle Y_{h}^{(n)}\rangle_{\kappa}$ in (eq. 20) as the convex
combination

$\displaystyle\langle
Y_{h}^{(n)}\rangle_{\kappa}=\sum_{i,j=1}^{p-1}\omega_{i}\omega_{j}\frac{J_{\kappa}^{ij}}{|\kappa|}Y_{\kappa}^{ij,n}+\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\tilde{\omega}_{\kappa}({\textbf{x}}_{e}^{k})\frac{J_{\kappa}({\textbf{x}}_{e}^{k})}{|\kappa|}Y_{h}(x_{e}^{k},t^{(n)}),$

then following [7, Sec. 4.2], we add $\tfrac{\Delta t^{(n)}}{|\kappa|}$ times
the flux balance in (eq. 79) to the cell-averaged scheme to get

$\displaystyle\langle Y_{h}^{(n+1)}\rangle_{\kappa}=$
$\displaystyle\sum_{i,j=1}^{p-1}\omega_{i}\omega_{j}\frac{J_{\kappa}^{ij}}{|\kappa|}Y_{\kappa}^{ij,n}+\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\tilde{\omega}_{\kappa}({\textbf{x}}_{e}^{k})\frac{J_{\kappa}({\textbf{x}}_{e}^{k})}{|\kappa|}\bigg{(}Y_{h}^{-}({\bf
x}_{e}^{k},t^{(n)})$ $\displaystyle-\frac{\Delta
t^{(n)}\omega_{k}J_{e}^{k}}{\tilde{\omega}_{\kappa}({\textbf{x}}_{e}^{k})J_{\kappa}({\textbf{x}}_{e}^{k})}\Big{(}{\bf
h}_{Y}^{\scalebox{0.8}{$\scriptscriptstyle\rm HLLC$}}\big{(}{\bf
u}_{h}^{-}({\bf x}_{e}^{k},t^{(n)}),{\bf u}_{h}^{+}({\bf
x}_{e}^{k},t^{(n)}),{\bf n}_{e}^{k}\big{)}-{\bf
h}_{Y}^{\scalebox{0.8}{$\scriptscriptstyle\rm Rus$}}\big{(}{\bf
u}_{\kappa}^{\star,n},{\bf u}_{h}^{-}({\bf x}_{e}^{k},t^{(n)}),{\bf
n}_{e}^{k}\big{)}\Big{)}\bigg{)}.$

As a consequence, the first $d+2$ components of
$\langle{\textbf{u}}_{h}^{(n+1)}\rangle_{\kappa}$ are convex combinations of
quantities in $\Omega_{\scalebox{0.8}{$\scriptscriptstyle\rm GM$}}$ under (eq.
80a) which confirms positivity of $\langle\rho_{h}^{(n+1)}\rangle_{\kappa}$
and $\rho e(\langle{\textbf{u}}_{h}^{(n+1)}\rangle_{\kappa})$ by concavity of
$\rho e({\textbf{u}})=\rho
E-\tfrac{\rho{\textbf{v}}\cdot\rho{\textbf{v}}}{2\rho}$. Excluding material
interfaces, we also get the minimum entropy principle.

Likewise, the two last components $Y$ in $\\{\Gamma,\Pi\\}$ satisfy

$\displaystyle\langle Y_{h}^{(n+1)}\rangle_{\kappa}$ $\displaystyle=\langle
Y_{h}^{(n)}\rangle_{\kappa}-\frac{\Delta
t^{(n)}}{|\kappa|}\Bigg{(}\sum^{p}_{i,j,k=0}\omega_{i}\omega_{j}\left(D_{ik}{\textbf{v}}^{ij,n}_{\kappa}\cdot{\textbf{n}}_{(i,k)j}Y_{\kappa}^{kj,n}+D_{jk}{\textbf{v}}^{ij,n}_{\kappa}\cdot{\textbf{n}}_{i(j,k)}Y_{\kappa}^{ik,n}\right)+\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\omega_{k}J_{e}^{k}{\textbf{D}}_{Y}^{-}\big{(}{\bf
u}_{h}^{-}({\bf x}_{e}^{k},t^{(n)}),{\bf u}_{h}^{+}({\bf
x}_{e}^{k},t^{(n)}),{\bf n}_{e}^{k}\big{)}\Bigg{)}.$

Inverting indices $i$ and $k$, then $j$ and $k$ in the volume integral, and
using the form of the interface fluctuations $D_{Y}^{-}({\bf u}^{-},{\bf
u}^{+},{\bf n})=\min(s^{\star},0)(Y^{+}-Y^{-})$, and observing that by (eq.
18) we have $\sum_{0\leq i,j\leq
p}\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\phi_{\kappa}^{ij}({\textbf{x}}_{e}^{k})f({\textbf{x}}_{e}^{k})=\sum_{e\in\partial\kappa}\sum_{k=0}^{p}f({\textbf{x}}_{e}^{k})$,
we get

$\displaystyle\langle Y_{h}^{(n+1)}\rangle_{\kappa}$
$\displaystyle=\sum_{i,j=0}^{p}\left(\omega_{i}\omega_{j}\frac{J_{\kappa}^{ij}}{|\kappa|}-\frac{\Delta
t^{(n)}}{|\kappa|}\sum_{k=0}^{p}\Big{(}\omega_{k}\omega_{j}D_{ki}{\bf
v}_{\kappa}^{kj,n}\cdot{\textbf{n}}_{(i,k)j}+\omega_{i}\omega_{k}D_{kj}{\bf
v}_{\kappa}^{ik,n}\cdot{\textbf{n}}_{i(j,k)}-\sum_{e\in\partial\kappa}\omega_{k}J_{e}^{k}\min\big{(}s^{\star}({\textbf{x}}_{e}^{k},t^{(n)}),0\big{)}\Big{)}\right)Y_{\kappa}^{ij,n}$
$\displaystyle-\frac{\Delta
t^{(n)}}{|\kappa|}\sum_{e\in\partial\kappa}\sum_{k=0}^{p}\omega_{k}J_{e}^{k}\min\big{(}s^{\star}({\textbf{x}}_{e}^{k},t^{(n)}),0\big{)}Y_{h}^{+}({\textbf{x}}_{e}^{k},t^{(n)}),$

and invoking the metric identities in (eq. 39), $\langle
Y_{h}^{(n+1)}\rangle_{\kappa}$ is indeed a convex combination of the DOFs
associated to $Y_{h}^{(n)}$ under (eq. 80). Positivity and the minimum and
maximum principles follow directly. ∎

Note that (eq. 80) correspond to sharp estimates of the time step. Assuming
hypercube elements of size $h$ and estimating the wave speeds by
$\lambda_{\infty}$, (eq. 80a) leads to $\Delta
t^{(n)}\tfrac{\lambda_{\infty}}{h}\leqslant\tfrac{1}{2dp(p+1)}$, which results
in a time step reduced by a factor $2$ when compared with [70, Th. 3.1] which
however assumes Cartesian meshes. Then, assuming a uniform velocity field,
${\textbf{v}}_{\kappa}^{ij,n}={\textbf{v}}$, (eq. 80b) leads to $\Delta
t^{(n)}\tfrac{\|{\bf v}\|_{\infty}}{h}\leqslant\tfrac{1}{dp(p+1)}$, thus
indicating that this condition will be usually satisfied by (eq. 80a).

##  7 A posteriori limiters

Properties of the discrete DGSEM in Theorem 6.1 hold for the cell-averaged
solution and a posteriori limiters are applied at the end of each Runge-Kutta
stage to extend these properties to all DOFs within elements. Here we describe
these limiters for (eq. 78) that are similar to the ones proposed in [71, 70,
15, 67]. As per their basic principle the DOFs are limited at time $t^{(n+1)}$
through a linear scaling around the cell-average (eq. 20):

$\tilde{{\textbf{U}}}^{ij,n+1}_{\kappa}=\theta_{\kappa}\big{(}{\textbf{U}}^{ij,n+1}_{\kappa}-\langle{\textbf{u}}_{h}^{(n+1)}\rangle_{\kappa}\big{)}+\langle{\textbf{u}}_{h}^{(n+1)}\rangle_{\kappa}\quad\forall
0\leqslant i,j\leqslant p,\quad\kappa\in\Omega_{h},$

where $0\leqslant\theta_{\kappa}\leqslant 1$ is the limiter coefficient. We
here apply successive limiters $(\theta_{\kappa}^{\rho},\theta_{\kappa}^{\rho
e},\theta_{\kappa}^{\Gamma},\theta_{\kappa}^{\Pi})$ on:

  * •

mixture density:

$\tilde{\rho}^{ij,n+1}_{\kappa}=\theta^{\rho}_{\kappa}\big{(}\rho^{ij,n+1}_{\kappa}-\langle\rho_{h}^{(n+1)}\rangle_{\kappa}\big{)}+\langle\rho_{h}^{(n+1)}\rangle_{\kappa},\quad\theta^{\rho}_{\kappa}=\min\left(\frac{\langle\rho_{h}^{(n+1)}\rangle_{\kappa}-\epsilon}{\langle\rho_{h}^{(n+1)}\rangle_{\kappa}-\rho^{\min}_{\kappa}},1\right),\quad\rho^{\min}_{\kappa}=\min_{0\leqslant
i,j\leqslant p}\rho^{ij,n+1}_{\kappa};$ (81)

  * •

EOS parameters $\Gamma$ and $\Pi$:

$\tilde{Y}^{ij,n+1}_{\kappa}=\theta^{Y}_{\kappa}\left(Y^{ij,n+1}_{\kappa}-\langle
Y_{h}^{(n+1)}\rangle_{\kappa}\right)+\langle
Y_{h}^{(n+1)}\rangle_{\kappa},\quad\theta^{Y}_{\kappa}=\min\left(\frac{\langle
Y_{h}^{(n+1)}\rangle_{\kappa}-m_{Y}}{\langle
Y_{h}^{(n+1)}\rangle_{\kappa}-Y^{\min}_{\kappa}},\frac{M_{Y}-\langle
Y_{h}^{(n+1)}\rangle_{\kappa}}{Y^{\max}_{\kappa}-\langle
Y_{h}^{(n+1)}\rangle_{\kappa}},1\right),\quad Y\in\\{\Gamma,\Pi\\},$ (82)

where $Y^{\min}_{\kappa}=\min_{0\leqslant i,j\leqslant p}Y^{ij,n+1}_{\kappa}$,
$Y^{\max}_{\kappa}=\max_{0\leqslant i,j\leqslant p}Y^{ij,n+1}_{\kappa}$, and

$m_{\Gamma}=\min_{1\leqslant i\leqslant n_{s}}\frac{1}{\gamma_{i}-1},\quad
M_{\Gamma}=\max_{1\leqslant i\leqslant n_{s}}\frac{1}{\gamma_{i}-1},\quad
m_{\Pi}=\min_{1\leqslant i\leqslant
n_{s}}\frac{\gamma_{i}\mathrm{p}_{\infty_{i}}}{\gamma_{i}-1},\quad
M_{\Pi}=\max_{1\leqslant i\leqslant
n_{s}}\frac{\gamma_{i}\mathrm{p}_{\infty_{i}}}{\gamma_{i}-1};$

  * •

mixture total internal energy:

$\tilde{Y}^{ij,n+1}_{\kappa}=\theta^{\rho
e}_{\kappa}\big{(}Y^{ij,n+1}_{\kappa}-\langle
Y_{h}^{(n+1)}\rangle_{\kappa}\big{)}+\langle
Y_{h}^{(n+1)}\rangle_{\kappa},\quad Y\in\\{\rho,\rho{\textbf{v}},\rho E\\},$
(83)

where $\rho e({\textbf{u}})=\rho
E-\tfrac{\rho{\textbf{v}}\cdot\rho{\textbf{v}}}{2\rho}$ and

$\theta^{\rho e}_{\kappa}=\min\left(\min_{0\leqslant i,j\leqslant
p}\left(\frac{\rho
e(\langle{\textbf{u}}_{h}\rangle^{n+1}_{\kappa})-\tilde{\mathrm{p}}_{\infty_{\kappa}}^{ij,n+1}-\epsilon}{\rho
e(\langle{\textbf{u}}_{h}\rangle^{n+1}_{\kappa})-\rho
e^{ij,n+1}_{\kappa}}\right),1\right).$ (84)

Note that in (eq. 81) and (eq. 84), $0<\epsilon\ll 1$ is a small parameter,
which we set as $\epsilon=10^{-8}$ in our numerical tests. The limiters (eq.
81) and (eq. 83) thus guarantee that $\tilde{\rho}^{0\leqslant ij\leqslant
p,n+1}_{\kappa}>0$ and $\tilde{\rho e}^{0\leqslant ij\leqslant
p,n+1}_{\kappa}>\tilde{\mathrm{p}}_{\infty_{\kappa}}^{0\leqslant ij\leqslant
p,n+1}$. Recalling (eq. 6), we have $m_{Y}\leqslant Y\leqslant M_{Y}$, with