# High-order variational Lagrangian schemes for compressible fluids

Guosheng Fu  Department of Applied and Computational Mathematics and
Statistics, University of Notre Dame, USA<EMAIL_ADDRESS>and  Chun Liu
Department of Applied Mathematics, Illinois Institute of Technology, USA.
<EMAIL_ADDRESS>
###### Abstract.

We present high-order variational Lagrangian finite element methods for
compressible fluids using a discrete energetic variational approach. Our
spatial discretization is mass/momentum/energy conserving and entropy stable.
Fully implicit time stepping is used for the temporal discretization, which
allows for a much larger time step size for stability compared to explicit
methods, especially for low-Mach number flows and/or on highly distorted
meshes. Ample numerical results are presented to showcase the good performance
of our proposed scheme.

###### Key words and phrases:

Discrete Energetic Variational Approach; Lagrangian Hydrodynamics; High-order
finite elements; Entropy stability

###### 1991 Mathematics Subject Classification:

65N30, 65N12, 76S05, 76D07

G. Fu’s research is partially supported by NSF grant DMS-2012031. C. Liu’s
research is partially supported by NSF grants DMS-1950868, DMS-2153029 and
DMS-2118181.

##  1\. Introduction

We are interested in the numerical simulation of multidimensional compressible
fluid flows. There are two typical choices: a Lagrangian framework in which
the computational mesh moves with the fluid velocity, and an Eulerian
framework in which the fluid flows through a fixed spatial mesh. Lagrangian
methods are widely used in many fields for multi-material flow simulations
such as astrophysics and inertial confinement fusion (ICF), due to their
distinguished advantage of being able to preserve the material interface
sharply. In this article, we are concerned exclusively with Lagrangian
methods.

To address problems in Lagrangian hydrodynamics, two fundamental approaches
are commonly used. The first one is the staggered grid hydrodynamics (SGH)
approach which employs a spatial discretization where the kinematic variables
such as position and velocity are continuous and defined at the mesh vertices,
and the thermodynamic variables including density, pressure, and internal
energy are discontinuous and defined at cell centers; see [35, 36]. Artificial
viscosity, as originally proposed in [35], is used to generate entropy
production and suppress numerical oscillations across shocks. The second
approach is the cell-centered hydrodynamics (CCH), which treats all
hydrodynamic variables as cell centered and use approximate Riemann solvers to
reconstruct a continuous velocity field; see [24, 9, 21, 6]. This process
naturally introduces a sufficient level of dissipation at shock boundaries
without the need of artificial viscosity.

While the original SGH scheme [35] does not preserve total energy
conservation, it can be recovered using the notion of “corner masses” and
“corner forces” as in the compatible hydro approach [5]. Energy-conservative
high-order finite element schemes were further developed in [10] which can be
viewed as a high-order generalization of the SGH approach since the key idea
of using a continuous kinematic approximation space and a discontinuous
thermodynamic approximation space were adopted. Advantages of the high-order
Lagrangian schemes in [10] over the low-order schemes in [5] include the
ability to more accurately capture the geometry of the flow and maintain
robustness with respect to mesh motion using high-order curved meshes, sharper
shock resolution and better symmetry preservation for symmetric flows.

Our spatial discretization follows the SGH approach using a (high-order)
continuous finite element space for the kinematic variables, and a (high-
order) discontinuous finite element space for the thermodynamic variables,
which turns out to be closely related to the high-order finite element scheme
[10]. However, we use a completely different derivation. While most of the
existing Lagrangian schemes are obtained by directly discretizing the
underling PDE system, the starting point of our spatial discretization is an
energy dissipation law. In particular, we present a class of variational
Lagrangian schemes for compressible flows using a discrete Energetic
variational approach, which is an analogue to the energetic variational
approach (EnVarA) [19, 14] in a semi-discrete level.

For a given energy-dissipation law and the kinematic (transport) relation, the
EnVarA [19, 14] provides a general framework to determine the dynamics of
system in a unique and well-defined way, through two distinct variational
processes: Least Action Principle (LAP) and Maximum Dissipation Principle
(MDP). This approach is originated from pioneering work of Onsager [26, 27]
and Rayleigh [33], and has been successfully applied to build up many
mathematical models [19, 34, 11, 14]. Most of existing EnVarA literature
focuses on the isothermal case where temperature variation is not allowed.
Here we adopt the approach used in [20] to model a thermodynamic system with
temperature variation using EnVarA. The obtained model is then discretized in
space using the discrete EnVarA, which leads to a high-order variational
Lagrangian finite element scheme. Main structures of the continuous model
including mass/momentum/energy conservation and entropy stability are
naturally preserved in the proposed spatial discretization. The resulting ODE
system is further discretized in time using high-order implicit time
integrators. Due to the special structure of the ODE system, a nonlinear
system only for the velocity degrees of freedom (DOFs) needs to be solved in
each time step. The allowed time step size for stability is drastically
improved over explicit time stepping, especially in the low Mach number
regime, at the expense of a nonlinear system solve. Ample numerical examples
are used to show the good performance of the proposed scheme.

We summarize the main features of our scheme:

  * •

Space-time high-order accuracy.

  * •

Mass/momentum/energy conservation and entropy stability for the spatial
discretization.

  * •

Implicit time stepping allows for large time step size especially for low Mach
number flows and/or on highly distorted meshes.

The rest of the paper is organized as follows. In Section 2, we present the
ideal gas model using EnVarA. In Section 3, a variational Lagrangian scheme is
constructed using discrete EnVarA. The implicit temporal discretization of the
resulting ODE system from Section 3 is then presented in Section 4. Numerical
examples are presented in Section 5, followed by a summary in Section 6. In
the Appendix, we briefly discuss our approach to the compressible isothermal
case where temperature is fixed.

##  2\. The energetic variational approach for ideal gas

The EnVarA [14] is a tool to derive the force balance equation starting from a
total energy $E^{\text{total}}$ and an energy dissipation rate functional
$\mathcal{D}$, where $E^{\text{total}}$ is the sum of the kinetic energy
$\mathcal{K}=\int_{\Omega}\frac{1}{2}\rho|\hbox{\boldmath$u$}|^{2}\mathrm{dx}$
and the Helmholtz free energy
$\mathcal{F}=\int_{\Omega}\psi(\rho,\theta)\mathrm{dx}$, and

$\displaystyle
2\mathcal{D}=\int_{\Omega}\eta|\nabla_{s}\hbox{\boldmath$u$}|^{2}+(\xi-\frac{2}{3}\eta)|\nabla\cdot\hbox{\boldmath$u$}|^{2}\mathrm{dx}$
(1)

is the rate of energy dissipation with dynamic viscosity $\eta$ and bulk
viscosity $\xi$, which may depend on both $\rho$ and $\theta$. Here
$\nabla_{s}$ denotes the symmetric gradient operator, $\rho$ is the density,
$\theta$ is the absolute temperature, and $u$ is the fluid velocity.

Thermodynamics of idea gas is well studied in the literature [1, 2, 3, 7, 8,
12, 22, 29]. In classical thermodynamics, the concept of free energy proves to
be useful [1, 29]. The internal energy and pressure of an ideal gas [22, 2]
have a linear relationship with temperature and the product of temperature and
density, respectively. Using this observation, Liu and Sulzbach [20] proposed
a definition for the free energy density of an ideal gas:

$\displaystyle\psi(\rho,\theta)=(c_{p}-c_{v})\theta\rho\log(\rho)-c_{v}\rho\theta\log(\theta),$
(2)

which we utilize in our current work. This definition includes the specific
heat at constant volume $c_{v}$ and specific heat at constant pressure
$c_{p}$. Associated with the free energy density (2), we define the three
thermodynamic variables, namely pressure $p$, internal energy (per unit mass)
$e$, and entropy (per unit mass) $s$:

$\displaystyle p:=$
$\displaystyle\;\psi_{\rho}\rho-\psi=(c_{p}-c_{v})\rho\theta,$ (3a)
$\displaystyle e:=$
$\displaystyle\;(\psi-\psi_{\theta}\theta)/\rho=c_{v}\theta,$ (3b)
$\displaystyle s:=$
$\displaystyle\;-\psi_{\theta}/\rho=\log(\theta^{c_{v}}/\rho^{c_{p}-c_{v}})+c_{v},$
(3c)

which are easily verified to satisfy the famous Gibbs equation (by chain rule)
[8]:

$\displaystyle\theta\partial s=\partial e+p\partial(\frac{1}{\rho}),$ (4)

where $\partial$ represents any differentiation.

In Lagrangian coordinates, we introduce the flow map:
$\hbox{\boldmath$x$}(\hbox{\boldmath$X$},t):\Omega^{0}\rightarrow\Omega^{t}$,
which satisfies the trajectory equation

$\displaystyle\frac{d}{dt}\hbox{\boldmath$x$}(\hbox{\boldmath$X$},t)=\hbox{\boldmath$u$}(\hbox{\boldmath$x$}(\hbox{\boldmath$X$},t),t),$
(5)

with initial condition
$\hbox{\boldmath$x$}(\hbox{\boldmath$X$},0)=\hbox{\boldmath$X$}$. Since we are
concerning with conventional ideal gas, we postulate the kinematics of the
temperature $\theta$ being transported along the trajectory:
$\frac{d}{dt}\theta=\theta_{t}+\hbox{\boldmath$u$}\cdot\nabla\theta$.

###  2.1. Kinematics: mass conservation

Within an Lagrangian control volume
$V^{t}:=\\{\hbox{\boldmath$x$}(\hbox{\boldmath$X$},t):\;\forall\hbox{\boldmath$X$}\in
V^{0}\\}$, mass does not change over time:

$\displaystyle\frac{d}{dt}\int_{V^{t}}\rho\,\mathrm{dx}=\int_{V^{0}}\frac{d}{dt}(\rho
J)\,\mathrm{dX}=0,$ (6)

where $J=\mathrm{Det}(\nabla_{X}\hbox{\boldmath$x$})$ is the Jacobian
determinant. Since equality (6) is valid for any control volume $V(t)$, there
must hold

$\displaystyle\frac{d}{dt}(\rho J)=0,\quad\text{ or
}\quad\rho(\hbox{\boldmath$x$}(\hbox{\boldmath$X$},t),t)=\rho_{0}(\hbox{\boldmath$X$})/J(\hbox{\boldmath$X$},t),$
(7)

where $\rho_{0}:\Omega^{0}\rightarrow\mathbb{R}^{+}$ is the initial density.
The equality (7) represents strong mass conservation. Writing mass
conservation (7) back to the Eulerian coordinates, we have
$\rho_{t}+\nabla\cdot(\rho\hbox{\boldmath$u$})=0$, which is often referred to
as the continuity equation. By abuse of notation, from it we can also get
$\delta\rho=-\nabla\cdot(\rho\delta\hbox{\boldmath$x$})$, where $\delta$
represents the variational derivative. Such relation will be used and made
clear in the derivations later in this paper.

###  2.2. Force balance: LAP and MDP

Here we combine the Least Action Principle (LAP) and Maximum Dissipation
Principle (MDP) [26, 27, 33] to derive the force balance equation. The action
functional for the system is

$\mathcal{A}=\int_{0}^{T}\left(\mathcal{K}-\mathcal{F}\right)\mathrm{dt}=\int_{0}^{T}\int_{\Omega^{t}}\left(\frac{1}{2}\rho|\hbox{\boldmath$u$}|^{2}-\psi(\rho,\theta)\right)\mathrm{dx}\mathrm{dt}.$

The LAP performs variation on the kinetic and free energies to derive the
inertial and conservative forces:

$\hbox{\boldmath$f$}_{\text{inertial}}=\frac{\delta\int_{0}^{T}\mathcal{K}\mathrm{dt}}{\delta\hbox{\boldmath$x$}},\quad\hbox{\boldmath$f$}_{\text{cons}}=\frac{\delta\int_{0}^{T}\mathcal{F}\mathrm{dt}}{\delta\hbox{\boldmath$x$}}.$

Taking variation on the kinetic energy and using mass conservation (7), we get

$\displaystyle\delta\int_{0}^{T}\mathcal{K}\mathrm{dt}=$
$\displaystyle\delta\int_{o}^{T}\int_{\Omega^{0}}\frac{1}{2}\rho_{0}(X)|\hbox{\boldmath$x$}_{t}|^{2}\mathrm{dXdt}=\int_{o}^{T}\int_{\Omega^{0}}\rho_{0}(X)\hbox{\boldmath$x$}_{t}\,\delta\hbox{\boldmath$x$}_{t}\,\mathrm{dXdt}$
$\displaystyle=$
$\displaystyle-\int_{o}^{T}\int_{\Omega^{0}}\rho_{0}(X)\hbox{\boldmath$x$}_{tt}\,\delta\hbox{\boldmath$x$}\,\mathrm{dXdt}$
$\displaystyle=$
$\displaystyle-\int_{o}^{T}\int_{\Omega^{0}}\rho_{0}(X)(\hbox{\boldmath$u$}_{t}+\hbox{\boldmath$u$}\cdot\nabla\hbox{\boldmath$u$})\,\delta\hbox{\boldmath$x$}\,\mathrm{dXdt}$
$\displaystyle=$
$\displaystyle-\int_{o}^{T}\int_{\Omega^{t}}\rho(\hbox{\boldmath$x$},t)\dot{\hbox{\boldmath$u$}}\,\delta\hbox{\boldmath$x$}\,\mathrm{dxdt}$
(8)

where we used the short-hand notation
$\dot{\hbox{\boldmath$u$}}:=\frac{d}{dt}\hbox{\boldmath$u$}(\hbox{\boldmath$x$}(\hbox{\boldmath$X$},t),t)$
for the material derivative. This implies that
$\hbox{\boldmath$f$}_{\text{inertial}}=-\rho\dot{\hbox{\boldmath$u$}}.$ Taking
variation on the Helmholtz free energy $\mathcal{F}$ (Hamilton’s principle of
virtual work) leads to

$\displaystyle\delta\int_{0}^{T}\mathcal{F}\mathrm{dt}=$
$\displaystyle\;\delta\int_{0}^{T}\int_{\Omega^{t}}\psi(\rho,\theta)\,\mathrm{dxdt}=\;\int_{0}^{T}\int_{\Omega^{t}}\psi_{\rho}\delta\rho+\psi_{\theta}\delta\theta\,\mathrm{dxdt}$
$\displaystyle=$
$\displaystyle\;\int_{0}^{T}\int_{\Omega^{t}}\psi_{\rho}(-\nabla\cdot(\rho\delta\hbox{\boldmath$x$}))+\psi_{\theta}(-\delta\hbox{\boldmath$x$}\cdot\nabla\theta)\,\mathrm{dxdt}$
$\displaystyle=$
$\displaystyle\;\int_{0}^{T}\int_{\Omega^{t}}(\rho\nabla\psi_{\rho}-\psi_{\theta}\nabla\theta)\cdot\delta\hbox{\boldmath$x$}\,\mathrm{dxdt}$
$\displaystyle=$ $\displaystyle\;\int_{0}^{T}\int_{\Omega^{t}}\nabla
p\cdot\delta\hbox{\boldmath$x$}\,\mathrm{dxdt}$ (9)

where we used mass conservation and the kinematic transport of the temperature
defined after formula (5) in the second row, and the definition of the
pressure (3a) in the last row:

$\nabla
p=\nabla(\psi_{\rho}\rho-\psi)=\rho\nabla\psi_{\rho}+\underbrace{\psi_{\rho}\nabla\rho-\psi_{\rho}\nabla\rho}_{=0}-\psi_{\theta}\nabla\theta.$

This gives the conservative force $\hbox{\boldmath$f$}_{\text{cons}}=\nabla
p.$

The MDP performs variation on the energy dissipation rate (1) to get the
dissipative force:

$\hbox{\boldmath$f$}_{\text{diss}}=\frac{\delta\mathcal{D}}{\delta\hbox{\boldmath$u$}},$

which implies

$\hbox{\boldmath$f$}_{\text{diss}}=-\nabla\cdot\left(\eta\nabla_{s}\hbox{\boldmath$u$}\cdot+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})\hbox{\boldmath$I$}\right),$

where $I$ is the identity matrix.

Combining these, we get the force balance equation [14]

$\displaystyle\frac{\delta\mathcal{A}}{\delta\hbox{\boldmath$x$}}=\frac{\delta\mathcal{D}}{\delta\hbox{\boldmath$u$}},$
(10)

which takes the following form

$\displaystyle\rho\dot{\hbox{\boldmath$u$}}+\nabla
p-\nabla\cdot\left(\eta\nabla_{s}\hbox{\boldmath$u$}+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})\hbox{\boldmath$I$}\right)=0,$
(11)

This is the usual momentum equation for the Navier-Stokes equations, with its
natural weak formulation

$\displaystyle\int_{\Omega^{t}}\left(\rho\dot{\hbox{\boldmath$u$}}\cdot\hbox{\boldmath$w$}-p\nabla\cdot\hbox{\boldmath$w$}+\eta\nabla_{s}\hbox{\boldmath$u$}\cdot\nabla_{s}\hbox{\boldmath$w$}+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})(\nabla\cdot\hbox{\boldmath$w$})\right)\mathrm{dx}=0,$
(12)

for any test function $w$ with homogeneous boundary conditions.

###  2.3. Internal energy and entropy equations

The Gibbs equation (4) naturally gives an update equation for the internal
energy:

$\displaystyle\rho\dot{e}=\rho\theta\dot{s}-p\rho\dot{(1/\rho)}=\rho\theta\dot{s}-p\nabla\cdot\hbox{\boldmath$u$}$
(13)

The rate of change of the entropy $s$ can be contributed by the entropy flux
$j$ and entropy production $\Delta$:

$\displaystyle\rho\dot{s}=\nabla\cdot\hbox{\boldmath$j$}+\Delta.$ (14)

The second law of thermodynamics states that entropy production is non-
negative: $\Delta\geq 0$. If we take the entropy flux $j$ being given by
Durhem relation

$\displaystyle\hbox{\boldmath$j$}\theta=q,$ (15)

and the heat flux

$\displaystyle q=\kappa\nabla\theta$ (16)

according to Fourier’s law in which $\kappa$ is the heat conductivity. From
here, the explicit expression of $\Delta$ can be derived via conservation of
total energy: there holds

$\displaystyle
0=\frac{d}{dt}\int_{\Omega}\left(\frac{1}{2}\rho|\hbox{\boldmath$u$}|^{2}+\rho
e\right)\mathrm{dx}=$
$\displaystyle\int_{\Omega}\rho\dot{\hbox{\boldmath$u$}}\cdot\hbox{\boldmath$u$}+\rho\dot{e}\mathrm{dx}.$
(17)

Using equations (12) with test function
$\hbox{\boldmath$w$}=\hbox{\boldmath$u$}$, (13), and (14), we get:

$\displaystyle 0=$
$\displaystyle\;\int_{\Omega}-(\eta|\nabla_{s}\hbox{\boldmath$u$}|^{2}+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})^{2})+\theta\left(\nabla\cdot(\frac{\kappa\nabla\theta}{\theta})+\Delta\right)\mathrm{dx}$
$\displaystyle=$
$\displaystyle\;\int_{\Omega}\left(-(\eta|\nabla_{s}\hbox{\boldmath$u$}|^{2}+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})^{2})-\frac{\kappa|\nabla\theta|^{2}}{\theta}+\theta\Delta\right)\mathrm{dx},$

where we applied the chain rule for the heat flux term and used the
homogeneous boundary condition $\nabla\theta\cdot\hbox{\boldmath$n$}=0$ on
$\partial\Omega$. This implies that we can take the entropy dissipation rate
as

$\displaystyle\Delta=\left(\eta|\nabla_{s}\hbox{\boldmath$u$}|^{2}+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})^{2}+\frac{\kappa|\nabla\theta|^{2}}{\theta}\right)/\theta,$
(18)

which satisfies the second law of thermodynamics as long as $\theta>0$. This
implies the following entropy equation:

$\displaystyle\rho\theta\dot{s}=\nabla\cdot(\kappa\nabla\theta)+\eta|\nabla_{s}\hbox{\boldmath$u$}|^{2}+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})^{2}.$
(19)

Plugging (14) and (18) back to the internal energy equation (13), we obtain:

$\displaystyle\rho\dot{e}=-p\nabla\cdot\hbox{\boldmath$u$}+\nabla\cdot(\kappa\nabla\theta)+\eta|\nabla_{s}\hbox{\boldmath$u$}|^{2}+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})^{2}.$
(20)

By (3b), equation (20) equivalently gives the dynamics of temperature
$\theta$, which has the heat equation as the leading term.

###  2.4. Summary

Combining the above results, we finally obtain the model equations:

$\displaystyle\dot{\hbox{\boldmath$x$}}=$
$\displaystyle\;\hbox{\boldmath$u$},$ (21a) $\displaystyle\dot{(\rho J)}=$
$\displaystyle\;0,$ (21b) $\displaystyle\rho\dot{\hbox{\boldmath$u$}}=$
$\displaystyle\;-\nabla
p+\nabla\cdot\left(\eta\nabla_{s}\hbox{\boldmath$u$}+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})\hbox{\boldmath$I$}\right),$
(21c) $\displaystyle\rho\dot{e}=$
$\displaystyle\;-p\nabla\cdot\hbox{\boldmath$u$}+\nabla\cdot(\kappa\nabla\theta)+(\eta|\nabla_{s}\hbox{\boldmath$u$}|^{2}+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})^{2}),$
(21d) where $e=c_{v}\theta,\quad p=(c_{p}-c_{v})\rho\theta.$ (21e)

This is nothing but the compressible Navier-Stokes equations of an ideal gas
[8]. Moreover, this system further satisfy the entropy equation (19). In the
next section we derive a variational Lagrangian scheme for this system, where
the spatial derivatives are evaluated by pulling back to the reference
configuration (Lagrange coordinates); see, e.g., [36, 6, 10].

##  3\. A discrete energetic variational approach: spatial discretization

In this section, we construct a variational Lagrangian scheme based on a
discrete EnVarA. For simplicity, we ignore heat conduction in the model, i.e.
we take $\kappa=0$ in this section.

###  3.1. Notation and the finite element spaces

Our grid-based scheme starts with a conforming triangulation
$\mathcal{T}_{h}^{0}=\\{T_{\ell}^{0}\\}_{\ell=1}^{N_{T}}$ of the initial
configuration $\Omega^{0}$ with $N_{T}$ elements, where we assume the element
$T_{\ell}^{0}:=\Phi_{T_{\ell}}^{0}(\widehat{T})$ is obtained from a polynomial
mapping $\Phi_{T_{\ell}}^{0}$ from the reference element $\widehat{T}$, which,
for simplicity, is a simplex or a hypercube.

We denote $\mathcal{P}^{k}(\widehat{T})$ as the polynomial space of degree no
greater than $k$ if $\widehat{T}$ is a reference simplex, or the tensor-
product polynomial space of degree no greater than $k$ in each direction if
$\widehat{T}$ is a reference hypercube, for $k\geq 1$. The mapped polynomial
space on a spatial physical element $T^{0}\in\mathcal{T}_{h}^{0}$ is denoted
as

$\mathcal{P}^{k}(T^{0}):=\\{\widehat{v}\circ(\Phi_{T}^{0})^{-1}:\;\forall\widehat{v}\in\mathcal{P}^{k}(\widehat{T}))\\}.$

We denote $\\{\widehat{\bm{\xi}}_{i}\\}_{i=1}^{N_{k}}$ as a set of quadrature
points with positive weights $\\{\widehat{\omega}_{i}\\}_{i=1}^{N_{k}}$ that
is accurate for polynomials of degree up to $2k+1$ on the reference element
$\widehat{T}$, i.e.,

$\displaystyle\int_{\widehat{T}}\widehat{f}\,\mathrm{dx}=\sum_{i=1}^{N_{k}}\widehat{\omega}_{i}\widehat{f}(\widehat{\bm{\xi}}_{i}),\quad\forall\widehat{f}\in\mathcal{P}^{2k+1}(\widehat{T}).$
(22)

Note that when $\widehat{T}$ is a reference square, we simply use the Gauss-
Legendre quadrature rule with $N_{k}=(k+1)^{2}$, which is optimal. On the
other hand, when $\widehat{T}$ is a reference simplex, the optimal choice of
quadrature rule is more complicated; see, e.g., [38, 37] and references cited
therein. Table 1 list the number $N_{k}$ for $0\leq k\leq 6$ of the symmetric
quadrature rules provided in [38].

| $k=0$ | $k=1$ | $k=2$ | $k=3$ | $k=4$ | $k=5$ | $k=6$
---|---|---|---|---|---|---|---
$N_{k}$ on Triangle | 1 | 6 | 7 | 15 | 19 | 28 | 37
$N_{k}$ on Tetrahedron | 1 | 8 | 14 | 36 | 61 | 109 | 171
Table 1. Number of quadrature points $N_{k}$ for the quadrature rule on a
simplex that is accurate up to degree $2k+1$ for $0\leq k\leq 6$.

The integration points and weights on a physical element $T_{\ell}^{0}$ are
simply obtained via mapping:
$\\{\bm{\xi}_{i}^{\ell}:=\Phi_{T_{\ell}^{0}}(\widehat{\bm{\xi}}_{i})\\}_{i=1}^{N_{k}}$,
and
$\\{\omega_{i}^{\ell}:=|\nabla\Phi_{T_{\ell}^{0}}(\widehat{\bm{\xi}}_{i})|\widehat{\omega}_{i}\\}_{i=1}^{N_{k}}$.
To simplify the notation, we denote the set of physical integration points and
weights

$\displaystyle\Xi_{h}^{k}:=$ $\displaystyle\;\\{\bm{\xi}_{i}^{\ell}:\;\;1\leq
i\leq N_{k},\,1\leq\ell\leq N_{T}\\},$ (23a) $\displaystyle\Omega_{h}^{k}:=$
$\displaystyle\;\\{\omega_{i}^{\ell}:\;\;1\leq i\leq N_{k},\,1\leq\ell\leq
N_{T}\\},$ (23b)

and denote $(\cdot,\cdot)_{h}$ as the discrete inner-product on the mesh
$\mathcal{T}_{h}^{0}$ using the quadrature points $\Xi_{h}^{k}$ and weights
$\Omega_{h}^{k}$:

$(\alpha,\beta)_{h}:=\sum_{\ell=1}^{N_{T}}\sum_{i=1}^{N_{k}}\alpha(\bm{\xi}_{i}^{\ell})\beta(\bm{\xi}_{i}^{\ell})\omega_{i}^{\ell}.$

We are now ready to present our continuous and discontinuous finite element
spaces:

$\displaystyle\hbox{\boldmath$V$}_{h}^{k}:=$
$\displaystyle\;\\{\hbox{\boldmath$v$}\in[H^{1}(\Omega^{0})]^{d}:\;\;\hbox{\boldmath$v$}|_{T_{\ell}^{0}}\in[\mathcal{P}^{k}(T_{\ell}^{0})]^{d},\;\;\forall
T_{\ell}^{0}\in\mathcal{T}_{h}^{0}\\},$ (24) $\displaystyle W_{h}^{k}:=$
$\displaystyle\;\\{w\in L^{2}(\Omega^{0}):\;\;w|_{T_{\ell}^{0}}\in
W^{k}(T_{\ell}^{0}),\;\;\forall T_{\ell}^{0}\in\mathcal{T}_{h}^{0}\\},$ (25)

where the local space

$W^{k}(T_{\ell}^{0}):=\mathcal{P}^{k}(T_{\ell}^{0})\oplus\delta
W_{k}(T_{\ell}^{0}),$

is associated with the integration rule in (22) such that $\dim
W^{k}(T_{\ell}^{0})=N_{k}$, and the nodal conditions

$\displaystyle\varphi_{i}^{\ell}({\bm{\xi}}_{j}^{\ell})=\delta_{ij},\quad\forall
1\leq j\leq N_{k},$ (26)

in which $\delta_{ij}$ is the Kronecker delta function determines a unique
solution $\varphi_{i}^{\ell}\in W^{k}(T_{\ell}^{0})$. This implies that
$\\{\varphi_{i}^{\ell}\\}_{i=1}^{N_{k}}$ is a set of nodal bases for the space
${W}^{k}(T_{\ell}^{0})$, i.e.,

$\displaystyle{W}^{k}(T_{\ell}^{0})=\mathrm{span}_{1\leq i\leq
N_{k}}\\{\varphi_{i}^{\ell}\\}.$ (27)

When $T_{\ell}^{0}$ is a mapped hypercube, we have $N_{k}=(k+1)^{2}$, hence
$\delta W_{k}(T_{\ell}^{0})=\emptyset$ and $W^{k}(T_{\ell})$ is simply the
(mapped) tensor product polynomial space $\mathcal{P}^{k}(T_{\ell}^{0})$. On
the other hand, when $T_{\ell}^{0}$ is a mapped simplex, we have $\dim\delta
W_{k}(T_{\ell}^{0})=N_{k}-\dim\mathcal{P}^{k}(T_{\ell}^{0})>0$ for $k\geq 1$.
However, we emphasize that the explicit expression of $\delta
W_{k}(T_{\ell}^{0})$ or the basis function $\phi_{i}^{\ell}$ does not matter
in our discretization, as only their nodal degrees of freedom (DOFs) on the
quadrature nodes will enter into the numerical integration. Any function
$\alpha(\hbox{\boldmath$X$},t)$ in $W_{h}^{k}$ (for fixed $t$) can be
expressed as

$\alpha(\hbox{\boldmath$X$},t)=\sum_{\ell=1}^{N_{T}}\sum_{i=1}^{N_{k}}{\sf
a}_{i}^{\ell}(t)\phi_{i}^{\ell}(\hbox{\boldmath$X$}),$

where $\\{{\sf a}_{i}^{\ell}(t)\\}$ are the unknown coefficients. We refer to
$W_{h}^{k}$ as the (discontinuous) integration rule space, which only
contains the $N_{T}\times N_{k}$ quadrature points and weights (23), and is
easy to implement in practice.

We further denote a set of basis functions for $\hbox{\boldmath$V$}_{h}^{k}$
as $\\{\hbox{\boldmath$\varphi$}_{i}\\}_{i=1}^{N_{V}^{k}}$, where $N_{V}^{k}$
is the dimension of $\hbox{\boldmath$V$}_{h}^{k}$. We use the continuous space
$\hbox{\boldmath$V$}_{h}^{k}$ to approximate the flow map
$\hbox{\boldmath$x$}_{h}$ and velocity $\hbox{\boldmath$u$}_{h}$, and the
discontinuous space $W_{h}^{k}$ to approximate the density $\rho_{h}$,
pressure $p_{h}$, internal energy $e_{h}$, temperature $\theta_{h}$, and
entropy $s_{h}$. More specifically, we have

$\displaystyle\hbox{\boldmath$x$}_{h}(\hbox{\boldmath$X$},t)=$
$\displaystyle\;\sum_{i=1}^{N_{V}^{k}}{\sf
x}_{i}(t)\hbox{\boldmath$\varphi$}_{i}(\hbox{\boldmath$X$}),\quad$
$\displaystyle\hbox{\boldmath$u$}_{h}(\hbox{\boldmath$X$},t)=\;\sum_{i=1}^{N_{V}^{k}}{\sf
u}_{i}(t)\hbox{\boldmath$\varphi$}_{i}(\hbox{\boldmath$X$}),$ (28a)
$\displaystyle\rho_{h}(\hbox{\boldmath$X$},t)=$
$\displaystyle\;\sum_{\ell=1}^{N_{T}}\sum_{i=1}^{N_{k}}{\rho}_{i}^{\ell}(t)\phi_{i}^{\ell}(\hbox{\boldmath$X$}),\quad$
$\displaystyle
p_{h}(\hbox{\boldmath$X$},t)=\;\sum_{\ell=1}^{N_{T}}\sum_{i=1}^{N_{k}}{\sf
p}_{i}^{\ell}(t)\phi_{i}^{\ell}(\hbox{\boldmath$X$}),$ (28b) $\displaystyle
e_{h}(\hbox{\boldmath$X$},t)=$
$\displaystyle\;\sum_{\ell=1}^{N_{T}}\sum_{i=1}^{N_{k}}{e}_{i}^{\ell}(t)\phi_{i}^{\ell}(\hbox{\boldmath$X$}),\quad$
$\displaystyle\theta_{h}(\hbox{\boldmath$X$},t)=\;\sum_{\ell=1}^{N_{T}}\sum_{i=1}^{N_{k}}{\theta}_{i}^{\ell}(t)\phi_{i}^{\ell}(\hbox{\boldmath$X$}),$
(28c) $\displaystyle s_{h}(\hbox{\boldmath$X$},t)=$
$\displaystyle\;\sum_{\ell=1}^{N_{T}}\sum_{i=1}^{N_{k}}{\sf
s}_{i}^{\ell}(t)\phi_{i}^{\ell}(\hbox{\boldmath$X$}),$ (28d)

where ${\sf X}_{h}=[{\sf x}_{1},\cdots,{\sf x}_{N_{V}^{k}}]^{T}$, ${\sf
U}_{h}=[{\sf u}_{1},\cdots,{\sf u}_{N_{V}^{k}}]^{T}$, ${\sf
R}_{h}=[{\sf\rho}_{1}^{1},\cdots,{\sf\rho}_{N_{k}}^{N_{T}}]^{T}$, ${\sf
P}_{h}=[{\sf p}_{1}^{1},\cdots,{\sf p}_{N_{k}}^{N_{T}}]^{T}$, ${\sf
E}_{h}=[{\sf e}_{1}^{1},\cdots,{\sf e}_{N_{k}}^{N_{T}}]^{T}$,
${\sf\Theta}_{h}=[{\sf\theta}_{1}^{1},\cdots,{\sf\theta}_{N_{k}}^{N_{T}}]^{T}$,
and ${\sf S}_{h}=[{\sf s}_{1}^{1},\cdots,{\sf s}_{N_{k}}^{N_{T}}]^{T}$ are the
time dependent coefficient vectors for $\hbox{\boldmath$x$}_{h}$,
$\hbox{\boldmath$u$}_{h}$, $\rho_{h}$, $p_{h}$, $e_{h}$, $\theta_{h}$, and
$s_{h}$, respectively. Note that by the thermodynamic relations (3), there are
only two independent thermodynamic variables. Here we take density $\rho$ and
temperature $\theta$ as the independent variables. The other variables will be
updated through the discrete formulas (31) below.

###  3.2. Trajectory equation and mass conservation

With the notation given in (28), the trajectory equation (5) simply implies
that

$\displaystyle{\sf X}_{h}^{\prime}(t)={\sf U}_{h}(t).$ (29)

We require mass conservation to be satisfied pointwise at the quadrature nodes
level, specifically, (7) implies that

$\displaystyle{\sf\rho}_{i}^{\ell}(t)=\rho_{0}(\hbox{\boldmath$\xi$}_{i}^{\ell})/J_{h}(\hbox{\boldmath$\xi$}_{i}^{\ell},t),\quad\forall
1\leq\ell\leq N_{T},\;1\leq i\leq N_{k},$ (30)

where the discrete Jacobian on the quadrature point
$\hbox{\boldmath$\xi$}_{i}^{\ell}$ is
$J_{h}(\hbox{\boldmath$\xi$}_{i}^{\ell},t):=\mathrm{Det}(\nabla_{X}\hbox{\boldmath$x$}_{h}(\hbox{\boldmath$\xi$}_{i}^{\ell},t))$.

###  3.3. Thermodynamic relations

We require the relations in (3) for the thermodynamic variables be satisfied
on the quadrature points level, which implies

$\displaystyle{\sf p}_{i}^{\ell}=$
$\displaystyle\;(c_{p}-c_{v})\rho_{i}^{\ell}\theta_{i}^{\ell},$ (31a)
$\displaystyle{\sf e}_{i}^{\ell}=$ $\displaystyle\;c_{v}\theta_{i}^{\ell},$
(31b) $\displaystyle{\sf s}_{i}^{\ell}=$
$\displaystyle\;c_{v}\log(\theta_{i}^{\ell})-(c_{p}-c_{v})\log(\rho_{i}^{\ell})+c_{v},$
(31c)

for all $1\leq i\leq N_{k}$ and $1\leq\ell\leq N_{T}$. It is easy to see that
the (pointwise) Gibbs equation (4) is satisfied:

$\displaystyle\rho_{i}^{\ell}\theta_{i}^{\ell}({\sf
s}_{i}^{\ell})^{\prime}=\rho_{i}^{\ell}({\sf
e}_{i}^{\ell})^{\prime}-\frac{{\sf
p}_{i}^{\ell}}{\rho_{i}^{\ell}}(\rho_{i}^{\ell})^{\prime}.$

By the density definition (30) and Jacobi’s formula, we have

$(\rho_{i}^{\ell})^{\prime}=-\rho_{i}^{\ell}\nabla\cdot\hbox{\boldmath$u$}_{h}(\hbox{\boldmath$\xi$}_{i}^{\ell}),$

which implies that

$\displaystyle\rho_{i}^{\ell}\theta_{i}^{\ell}({\sf
s}_{i}^{\ell})^{\prime}=\rho_{i}^{\ell}({\sf e}_{i}^{\ell})^{\prime}+{{\sf
p}_{i}^{\ell}}\nabla\cdot\hbox{\boldmath$u$}_{h}(\hbox{\boldmath$\xi$}_{i}^{\ell}).$
(32)

We note that the above pointwise relations also hold for the classical low-
order SGH schemes [35, 36, 5] where the thermodynamic variables are
approximated via piecewise constants, which, however, does not hold in general
for the high order scheme [10] due to the use of a different high-order
thermodynamic finite element space.

###  3.4. Discrete EnVarA and velocity equation

Instead of discretizing the force balance equation (21c), here we discretize
the energy law and use the EnVarA to derive the discrete force balance
equation directly. We denote the discrete action functional

$\displaystyle\mathcal{A}_{h}:=\int_{0}^{T}(\mathcal{K}_{h}-\mathcal{F}_{h})\mathrm{dt},$
(33)

where the discrete kinetic energy $\mathcal{K}_{h}$ and the discrete Helmholtz
free energy $\mathcal{F}_{h}$ are given as

$\mathcal{K}_{h}=\frac{1}{2}(\rho_{h}J_{h}\hbox{\boldmath$u$}_{h},\hbox{\boldmath$u$}_{h})_{h}=\sum_{\ell=1}^{N_{T}}\sum_{i=1}^{N_{k}}\frac{1}{2}\rho_{0}(\hbox{\boldmath$\xi$}_{i}^{\ell})|\hbox{\boldmath$u$}_{h}(\hbox{\boldmath$\xi$}_{i}^{\ell},t)|^{2}\omega_{i}^{\ell},$

and

$\mathcal{F}_{h}=(\psi(\rho_{h},\theta_{h})J_{h},1)_{h}=\sum_{\ell=1}^{N_{T}}\sum_{i=1}^{N_{k}}\psi\left(\rho_{h}(\hbox{\boldmath$\xi$}_{i}^{\ell},t),\theta_{h}(\hbox{\boldmath$\xi$}_{i}^{\ell},t)\right)J_{h}(\hbox{\boldmath$\xi$}_{i}^{\ell},t)\omega_{i}^{\ell}.$

Moreover, the discrete dissipation rate is given as

$\displaystyle\mathcal{D}_{h}:=\frac{1}{2}(\eta
J_{h}\nabla_{s}\hbox{\boldmath$u$}_{h},\nabla_{s}\hbox{\boldmath$u$}_{h})_{h}+\frac{1}{2}\left((\xi-\frac{2}{3}\eta)J_{h}\nabla\cdot\hbox{\boldmath$u$}_{h},\nabla\cdot\hbox{\boldmath$u$}_{h}\right)_{h}.$
(34)

The discrete force balance equation (10) is then

$\displaystyle\frac{\delta\mathcal{A}_{h}}{\delta{\sf
x}_{j}}=\frac{\delta\mathcal{D}_{h}}{\delta{\sf u}_{j}},\quad\forall 1\leq
j\leq N_{V}^{k}.$ (35)

Elementary calculation, using (29), (30) and Jacobi’s formula, yields that

$\displaystyle\frac{\delta\mathcal{A}_{h}}{\delta{\sf
x}_{j}}=-(\rho_{0}\dot{\hbox{\boldmath$u$}}_{h},\hbox{\boldmath$\varphi$}_{j})_{h}+(p_{h}J_{h},\nabla\cdot\hbox{\boldmath$\varphi$}_{j})_{h},$
(36)

where $\dot{\hbox{\boldmath$u$}}_{h}=\sum_{i=1}^{N_{V}^{k}}{\sf
u}_{i}(t)^{\prime}\hbox{\boldmath$\varphi$}_{i}$, and the pressure $p_{h}\in
W_{h}^{k}$ satisfies

$\displaystyle{\sf
p}_{i}^{\ell}=\psi_{\rho}(\rho_{i}^{\ell},\theta_{i}^{\ell})\rho_{i}^{\ell}-\psi(\rho_{i}^{\ell},\theta_{i}^{\ell})=(c_{p}-c_{v})\rho_{i}^{\ell}\theta_{i}^{\ell},$

according to (31a). We also have

$\displaystyle\frac{\delta\mathcal{D}_{h}}{\delta{\sf u}_{j}}=(\eta
J_{h}\nabla_{s}{\hbox{\boldmath$u$}}_{h},\nabla_{s}\hbox{\boldmath$\varphi$}_{j})_{h}+((\xi-\frac{2}{3}\eta)J_{h}\nabla\cdot{\hbox{\boldmath$u$}}_{h},\nabla\cdot\hbox{\boldmath$\varphi$}_{j})_{h},\quad\forall
1\leq j\leq N_{V}^{k}.$ (37)

Plugging (36) and (37) back to (35), and using the definition of the pressure,
we get the semi-discrete force balance equation:

$\displaystyle(\rho_{0}\dot{\hbox{\boldmath$u$}}_{h},\hbox{\boldmath$\varphi$}_{j})_{h}-\left((c_{p}-c_{v})\rho_{0}\theta_{h},\nabla\cdot\hbox{\boldmath$\varphi$}_{j}\right)_{h}+(\sigma_{h},\nabla\hbox{\boldmath$\varphi$}_{j})_{h}=0,\quad\forall
1\leq j\leq N_{V}^{k},$ (38)

where $\sigma_{h}$ is the viscous stress defined as

$\displaystyle\sigma_{h}:=\eta
J_{h}\nabla_{s}{\hbox{\boldmath$u$}}_{h}+(\xi-\frac{2}{3}\eta)J_{h}\nabla\cdot{\hbox{\boldmath$u$}}_{h}\hbox{\boldmath$I$}.$
(39)

Here the evaluation of spatial derivative terms shall be pulled back to the
initial configuration $\Omega^{0}$. In particular,

$\nabla\hbox{\boldmath$u$}_{h}=(F_{h})^{-1}\nabla_{X}\hbox{\boldmath$u$}_{h},$

where $F_{h}=\nabla_{X}\hbox{\boldmath$x$}_{h}$ is the deformation tensor.
Equation (38) provides an ODE system for the velocity coefficient vector ${\sf
U}_{h}(t)$. Taking test function $\hbox{\boldmath$\varphi$}_{j}$ as a constant
vector, we immediately obtain global momentum conservation:

$\frac{d}{dt}(\rho_{0}\hbox{\boldmath$u$}_{h},1)_{h}=(\rho_{0}\dot{\hbox{\boldmath$u$}}_{h},1)_{h}=0.$

###  3.5. Energy conservation and temperature equation

We use energy conservation to get an update equation for the temperature
coefficient vector ${\sf\Theta}_{h}(t)$. In the absence of heat conduction
($\kappa=0$), the spatial discretization of the internal energy equation (20)
leads to

$\displaystyle(\rho_{h}\dot{e}_{h},J_{h}\phi_{j}^{\ell})_{h}=-(p_{h}J_{h}\nabla\cdot\hbox{\boldmath$u$}_{h},\phi_{j}^{\ell})_{h}+(\sigma_{h}:\nabla\hbox{\boldmath$u$}_{h},\phi_{j}^{\ell})_{h},\quad\forall
1\leq\ell\leq N_{T},\,1\leq j\leq N_{k}.$ (40)

Equivalently, the equation (40) has the following pointwise form for the
coefficient vector ${\sf E}_{h}$:

$\displaystyle\rho_{0}(\hbox{\boldmath$\xi$}_{j}^{\ell})({\sf
e}_{j}^{\ell})^{\prime}=-{\sf
p}_{j}^{\ell}J_{h}(\hbox{\boldmath$\xi$}_{j}^{\ell})(\nabla\cdot\hbox{\boldmath$u$}_{h})(\hbox{\boldmath$\xi$}_{j}^{\ell})+\sigma_{h}(\hbox{\boldmath$\xi$}_{j}^{\ell}):\nabla\hbox{\boldmath$u$}_{h}(\hbox{\boldmath$\xi$}_{j}^{\ell}),$
(41)

for all $1\leq\ell\leq N_{T},\,1\leq j\leq N_{k}$. Plugging in the relations
(31a) and (31b) back to (41), we get the following ODE system for the
coefficent vector $\Theta_{h}(t)$:

$\displaystyle
c_{v}\rho_{0}(\hbox{\boldmath$\xi$}_{j}^{\ell})({\sf\theta}_{j}^{\ell})^{\prime}=-(c_{p}-c_{v})\rho_{0}(\hbox{\boldmath$\xi$}_{j}^{\ell})(\nabla\cdot\hbox{\boldmath$u$}_{h})(\hbox{\boldmath$\xi$}_{j}^{\ell}){\sf\theta}_{j}^{\ell}+\sigma_{h}(\hbox{\boldmath$\xi$}_{j}^{\ell}):\nabla\hbox{\boldmath$u$}_{h}(\hbox{\boldmath$\xi$}_{j}^{\ell}),$
(42)

One key observation is that (42) is a linear ODE system for ${\sf\Theta}_{h}$.
Moreover, taking $\hbox{\boldmath$\varphi$}_{j}=\hbox{\boldmath$u$}_{h}$ in
(38) and $\phi_{j}^{\ell}=1$ in (40) and adding, we obtain the total energy
conservation:

$\frac{d}{dt}\left(\frac{1}{2}\rho_{0}|\hbox{\boldmath$u$}_{h}|^{2}+\rho_{0}e_{h},1\right)_{h}=(\rho_{0}\dot{\hbox{\boldmath$u$}}_{h},\hbox{\boldmath$u$}_{h})_{h}+(\rho_{0}\dot{e}_{h},1)_{h}=0.$

###  3.6. The entropy equation

Combining the Gibbs equation (32) with the internal energy equation (41) and
simplifying, we obtain the ODE system satisfied by the entropy:

$\displaystyle
J_{h}(\hbox{\boldmath$\xi$}_{j}^{\ell})\rho_{i}^{\ell}\theta_{i}^{\ell}({\sf
s}_{i}^{\ell})^{\prime}=\sigma_{h}(\hbox{\boldmath$\xi$}_{j}^{\ell}):\nabla\hbox{\boldmath$u$}_{h}(\hbox{\boldmath$\xi$}_{j}^{\ell}),$
(43)

for all $1\leq\ell\leq N_{T}$, $1\leq j\leq N_{k}$. This implies that

$\displaystyle\frac{d}{dt}(\rho_{0}s_{h},1)_{h}=(\rho_{0}\dot{s}_{h},1)_{h}=(\frac{\sigma_{h}:\nabla\hbox{\boldmath$u$}_{h}}{\theta_{h}},1)_{h},$
(44)

where positivity of the right hand side, i.e., semi-discrete entropy
stability, is guaranteed as long as the temperature $\Theta_{h}>0$. We remark
that the entropy stability (44) is a direct consequence of our special choice
of (nodal) thermodynamic finite element space (25) and (27).

###  3.7. Artificial viscosity

To make the scheme robust even in the case of zero physical viscosities with
$\eta=\xi=0$, we add artificial viscosity [35] to the system so that shocks
can be dissipated. Specifically, we add to the stress term (39) an artificial
stress tensor $\sigma_{h}^{av}$ of the following form:

$\displaystyle\sigma_{h}\leftarrow\sigma_{h}+\sigma_{h}^{av},\quad\text{ where
}\sigma_{h}^{av}:=\mu_{av}J_{h}\nabla_{s}\hbox{\boldmath$u$}_{h},$ (45)

in which, following [10], the artificial viscosity coefficient $\mu_{av}$ is:

$\displaystyle\mu_{av}=\rho_{h}\left(q_{2}\ell_{s_{1}}^{2}|\Delta_{s_{1}}\hbox{\boldmath$u$}_{h}|+q_{1}\psi_{0}\psi_{1}\ell_{s_{1}}c_{s}\right)$
(46)

where $q_{1}$ and $q_{2}$ are linear and quadratic scaling coefficients,
$c_{s}=\sqrt{\gamma p_{h}/\rho_{h}}$ is the speed of sound with
$\gamma=c_{p}/c_{v}$ being the adiabatic constant,
$\Delta_{s_{1}}\hbox{\boldmath$u$}_{h}:={s_{1}\cdot\nabla\hbox{\boldmath$u$}_{h}\cdot
s_{1}}$ is the directional measure of compression and
$\ell_{s_{1}}=\ell_{0}|J_{h}s_{1}|$ is the directional length scale along the
direction $s_{1}$, and the two linear switches are
$\phi_{0}=\frac{|\nabla\cdot\hbox{\boldmath$u$}_{h}|}{\|\nabla\hbox{\boldmath$u$}_{h}\|}$,
and $\phi_{1}=\begin{cases}1,&\text{if
}\Delta_{s}\hbox{\boldmath$u$}_{h}<0,\\\ 0,&\text{if
}\Delta_{s}\hbox{\boldmath$u$}_{h}\geq 0.\\\ \end{cases}$ Here the direction
$s_{1}$ is the unit eigenvector of the symmetric tensor
$\nabla_{s}\hbox{\boldmath$u$}_{h}$ with the smallest eigenvalue
$\lambda_{1}$, i.e.,

$(\nabla_{s}\hbox{\boldmath$u$}_{h})s_{1}=\lambda_{1}s_{1},\quad|s_{1}|=1,\text{
and $\lambda_{1}$ is the smallest eigenvalue.}$

With this notation, we have
$\Delta_{s_{1}}\hbox{\boldmath$u$}_{h}=\lambda_{1}$. Moreover
$\ell_{0}=h_{0}/k$ is the mesh size of the initial domain divided by the
polynomial degree $k$. We refer interested reader to [10] and references cited
therein for more discussion on the choice of the artificial viscosity
coefficient.

###  3.8. Summary

The final form of the semi-discrete scheme is summarized in Algorithm 1 below.
This spatial discretization is high-order, mass/momentum/energy conserving,
and entropy stable.

Algorithm 1 Spatial discretization for model (21).

$\bullet$ Find
$\hbox{\boldmath$x$}_{h},\hbox{\boldmath$u$}_{h}\in\hbox{\boldmath$V$}_{h}^{k}$,
and $\theta_{h}\in W_{h}^{k}$ such that the ODE system  $\displaystyle{\sf
X}_{h}^{\prime}(t)=$ $\displaystyle\;{\sf U}_{h}(t),$
$\displaystyle(\rho_{0}\dot{\hbox{\boldmath$u$}}_{h},\hbox{\boldmath$\varphi$}_{h})_{h}-\left((c_{p}-c_{v})\rho_{0}\theta_{h},\nabla\cdot\hbox{\boldmath$\varphi$}_{h}\right)_{h}+(\sigma_{h},\nabla\hbox{\boldmath$\varphi$}_{h})_{h}=$
$\displaystyle\;0,\quad\forall\hbox{\boldmath$\varphi$}_{h}\in\hbox{\boldmath$V$}_{h}^{k},$
$\displaystyle({\sf\theta}_{j}^{\ell})^{\prime}+(\gamma-1)(\nabla\cdot\hbox{\boldmath$u$}_{h})(\hbox{\boldmath$\xi$}_{j}^{\ell}){\sf\theta}_{j}^{\ell}-\frac{\sigma_{h}(\hbox{\boldmath$\xi$}_{j}^{\ell})}{c_{v}\rho_{0}(\hbox{\boldmath$\xi$}_{j}^{\ell})}:\nabla\hbox{\boldmath$u$}_{h}(\hbox{\boldmath$\xi$}_{j}^{\ell})=$
$\displaystyle\;0,\quad\forall 1\leq\ell\leq N_{T},1\leq j\leq N_{k},$ holds
for the coefficient vectors ${\sf X}_{h}$, ${\sf U}_{h}$, and
${\sf\Theta}_{h}$, where the numerical stress
$\sigma_{h}=(\eta+\mu_{av})J_{h}\nabla_{s}{\hbox{\boldmath$u$}}_{h}+(\xi-\frac{2}{3}\eta)J_{h}\nabla\cdot{\hbox{\boldmath$u$}}_{h}\hbox{\boldmath$I$},$
in which the artificial viscosity $\mu_{av}$ is given in (46). Here the
notation (28) for the finite element approximations is used.

$\bullet$ The density approximation $\rho_{h}\in W_{h}^{k}$ satisfies mass
conservation (30), and the pressure, internal energy, and entropy
approximations $p_{h},e_{h},s_{h}\in W_{h}^{k}$ satisfy the thermodynamic
relations (31).

##  4\. Temporal discretization

In this section, we focus on the discretization of the ODE system in Algorithm
1. We use fully implicit time discretizations so that the fully discrete
scheme is robust for all mach numbers. We refer to [10, 30] for alternative
conservative explicit schemes.

###  4.1. First order energy dissipative scheme

Using implicit Euler for the time derivative terms, we arrive at the following
first order scheme: Given data
${\hbox{\boldmath$x$}}_{h}^{n-1},{\hbox{\boldmath$u$}}_{h}^{n-1}\in\hbox{\boldmath$V$}_{h}^{k}$
and ${\theta}_{h}^{n-1}\in W_{h}^{k}$ at time $t^{n-1}$, and time step size
$\delta t$, find solution
${\hbox{\boldmath$x$}}_{h}^{n},{\hbox{\boldmath$u$}}_{h}^{n}\in\hbox{\boldmath$V$}_{h}^{k}$
and ${\theta}_{h}^{n}\in W_{h}^{k}$ such that

$\displaystyle\frac{{\hbox{\boldmath$x$}}_{h}^{n}-{\hbox{\boldmath$x$}}_{h}^{n-1}}{\delta
t}=$ $\displaystyle\;{\hbox{\boldmath$u$}}_{h}^{n},$ (47a)
$\displaystyle(\rho_{0}\frac{{\hbox{\boldmath$u$}}_{h}^{n}-{\hbox{\boldmath$u$}}_{h}^{n-1}}{\delta
t},\hbox{\boldmath$\varphi$}_{h})_{h}-\left((c_{p}-c_{v})\rho_{0}\theta_{h}^{n},\nabla\cdot\hbox{\boldmath$\varphi$}_{h}\right)_{h}+(\sigma_{h}^{n},\nabla\hbox{\boldmath$\varphi$}_{h})_{h}=$
$\displaystyle\;0,\quad\forall\varphi_{h}\in\hbox{\boldmath$V$}_{h}^{k},$
(47b)
$\displaystyle\frac{{\sf\theta}_{j}^{\ell,n}-{\sf\theta}_{j}^{\ell,n-1}}{\delta
t}+(\gamma-1)\nabla\cdot\hbox{\boldmath$u$}_{h}^{n}(\hbox{\boldmath$\xi$}_{j}^{\ell}){\sf\theta}_{j}^{\ell,n}-\frac{\sigma_{h}^{n}(\hbox{\boldmath$\xi$}_{j}^{\ell}):\nabla\hbox{\boldmath$u$}_{h}^{n}(\hbox{\boldmath$\xi$}_{j}^{\ell})}{c_{v}\rho_{0}(\hbox{\boldmath$\xi$}_{j}^{\ell})}=$
$\displaystyle\;0,\quad\forall j,\ell,$ (47c)

where the stress

$\sigma_{h}^{n}=(\eta+\mu_{av}^{n-1})J_{h}^{n}\nabla_{s}{\hbox{\boldmath$u$}}_{h}^{n}+(\xi-\frac{2}{3}\eta)J_{h}^{n}\nabla\cdot{\hbox{\boldmath$u$}}_{h}^{n}\hbox{\boldmath$I$},$

in which the artificial viscosity coefficient $\mu_{av}^{n}$ is evaluated at
time $t^{n}$, and the Jacobian determinant
$J_{h}^{n}=|\nabla_{X}\hbox{\boldmath$x$}_{h}^{n}|$. The above system can be
solved by first expressing $\hbox{\boldmath$x$}_{h}^{n}$ and $\theta_{h}^{n}$
in terms of $\hbox{\boldmath$u$}_{h}^{n}$ using (47a) and (47c):

$\displaystyle\hbox{\boldmath$x$}_{h}^{n}=$
$\displaystyle\;\hbox{\boldmath$x$}_{h}^{n-1}+\delta
t\hbox{\boldmath$u$}_{h}^{n},$ (48a) $\displaystyle{\sf\theta}_{j}^{\ell,n}=$
$\displaystyle\;\frac{{\sf\theta}_{j}^{\ell,n-1}+\delta
t\frac{\sigma_{h}^{n}(\hbox{\boldmath$\xi$}_{i}^{\ell}):\nabla\hbox{\boldmath$u$}_{h}^{n}(\hbox{\boldmath$\xi$}_{i}^{\ell})}{c_{v}\rho_{0}(\hbox{\boldmath$\xi$}_{i}^{\ell})}}{1+\delta
t(\gamma-1)\nabla\cdot\hbox{\boldmath$u$}_{h}^{n}(\hbox{\boldmath$\xi$}_{i}^{\ell})},$
(48b)

and then solve the nonlinear system for $\hbox{\boldmath$u$}_{h}^{n}$ in (47b)
using (48). We use Newton’s method to solve this nonlinear system for the
velocity DOFs.

For the scheme (47), positivity of density on the quadrature points is
guaranteed as long as the Jacobian $J_{h}^{n}>0$ on these quadrature points.
And positivity of the temperature (hence positivity of pressure and internal
energy) is satisfied as long as the denominator of right hand side of (48b)
stays positive, i.e.,

$1+\delta
t(\gamma-1)\nabla\cdot\hbox{\boldmath$u$}_{h}^{n}(\hbox{\boldmath$\xi$}_{i}^{\ell})>0,\quad\forall
i,\ell.$

Moreover, strong mass conservation is satisfied due to (30), and global
momentum conservation is satisfied by taking $\hbox{\boldmath$\varphi$}_{j}$
to be a global constant in (47b). Finally, taking
$\hbox{\boldmath$\varphi$}_{j}=\hbox{\boldmath$u$}_{h}^{n}$ in (47b) and
combining with (47c), we get

$\displaystyle(\rho_{0}(\hbox{\boldmath$u$}_{h}^{n}-\hbox{\boldmath$u$}_{h}^{n-1}),\hbox{\boldmath$u$}_{h}^{n})_{h}+(\rho_{0}(e_{h}^{n}-e_{h}^{n-1}),1)_{h}=0,$

which implies that

$\displaystyle\left(\rho_{0}(\frac{1}{2}|\hbox{\boldmath$u$}_{h}^{n}|^{2}+e_{h}^{n}),1\right)_{h}-\left(\rho_{0}(\frac{1}{2}|\hbox{\boldmath$u$}_{h}^{n-1}|^{2}+e_{h}^{n-1}),1\right)_{h}=-\left(\rho_{0}\frac{1}{2}|\hbox{\boldmath$u$}_{h}^{n}-\hbox{\boldmath$u$}_{h}^{n-1}|^{2},1\right)_{h}\leq
0$

Hence, the total energy is dissipated over time for the scheme (47).

###  4.2. Second order energy-conservative scheme

Energy conservation can be recovered from the first order scheme (47) by
applying a time filter, which is the same as the midpoint rule; see [4]. This
time stepping algorithm has two steps and is recorded in Algorithm 2 for
reference.

Algorithm 2 The midpoint rule with Backward Euler – Forward Euler
implementation.

$\bullet$ Apply the backward Euler scheme (47) with half time step $\delta
t/2$ to get approximations at midpoint $t^{n-\frac{1}{2}}=t^{n-1}+\delta t/2$,
and denote the solutions as $\hbox{\boldmath$x$}_{h}^{n-\frac{1}{2}}$,
$\hbox{\boldmath$u$}_{h}^{n-\frac{1}{2}}$, and $\theta_{h}^{n-\frac{1}{2}}$.

$\bullet$ Apply a time filter (forward Euler) step to approximation the
solutions at time $t^{n}=t^{n-1}+\delta t$:
$\displaystyle\hbox{\boldmath$x$}_{h}^{n}=2\hbox{\boldmath$x$}_{h}^{n-\frac{1}{2}}-\hbox{\boldmath$x$}_{h}^{n-1},\;\;\hbox{\boldmath$u$}_{h}^{n}=2\hbox{\boldmath$u$}_{h}^{n-\frac{1}{2}}-\hbox{\boldmath$u$}_{h}^{n-1},\;\;\theta_{h}^{n}=2\theta_{h}^{n-\frac{1}{2}}-\theta_{h}^{n-1}.$
(49)

$\bullet$ The density, pressure, internal energy, and entropy approximations
$\rho_{h},p_{h},e_{h},s_{h}\in W_{h}^{k}$ are then recovered through (30) and
(31).

We notice that the BE step implies

$\displaystyle(\rho_{0}(\hbox{\boldmath$u$}_{h}^{n-\frac{1}{2}}-\hbox{\boldmath$u$}_{h}^{n-1}),\hbox{\boldmath$u$}_{h}^{n-\frac{1}{2}})_{h}+(\rho_{0}(e_{h}^{n-\frac{1}{2}}-e_{h}^{n-1}),1)_{h}=0,$

By the extrapolation relations (49), we have

$(\hbox{\boldmath$u$}_{h}^{n-\frac{1}{2}}-\hbox{\boldmath$u$}_{h}^{n-1})\cdot\hbox{\boldmath$u$}_{h}^{n-\frac{1}{2}}=\frac{1}{4}(|\hbox{\boldmath$u$}_{h}^{n}|^{2}-|\hbox{\boldmath$u$}_{h}^{n-1}|^{2}),\text{
and }e_{h}^{n-\frac{1}{2}}-e_{h}^{n-1}=\frac{1}{2}(e_{h}^{n-1}-e_{h}^{n-1}).$

Combining these equations, we get total energy conservation for this two-step
method:

$\left(\rho_{0}(\frac{1}{2}|\hbox{\boldmath$u$}_{h}^{n}|^{2}+e_{h}^{n}),1\right)_{h}=\left(\rho_{0}(\frac{1}{2}|\hbox{\boldmath$u$}_{h}^{n-1}|^{2}+e_{h}^{n-1}),1\right)_{h}.$

###  4.3. High-order BDF schemes

It is natural to extend the implicit Euler scheme (47) to higher order by
replacing the backward Euler time difference terms in (47) using higher-order
backward difference formulas (BDFs). For completeness, we record these high-
order schemes with uniform time stepping in Algoritm 3 below.

Algorithm 3 The BDF$[m]$ scheme.

$\bullet$ Given data
${\hbox{\boldmath$x$}}_{h}^{n-j},{\hbox{\boldmath$u$}}_{h}^{n-j}\in\hbox{\boldmath$V$}_{h}^{k}$
and ${\theta}_{h}^{n-j}\in W_{h}^{k}$ at time $t^{n-j}=(n-j)\delta t$ for
$j=1\cdots,m$, find solution
${\hbox{\boldmath$x$}}_{h}^{n},{\hbox{\boldmath$u$}}_{h}^{n}\in\hbox{\boldmath$V$}_{h}^{k}$
and ${\theta}_{h}^{n}\in W_{h}^{k}$ at time $t^{n}=n\delta t$ such that
$\displaystyle{D}_{\delta t}^{m}(\hbox{\boldmath$x$}_{h}^{n})=$
$\displaystyle\;{\hbox{\boldmath$u$}}_{h}^{n},$ (50a)
$\displaystyle(\rho_{0}{D_{\delta
t}^{m}({\hbox{\boldmath$u$}}_{h}^{n})},\hbox{\boldmath$\varphi$}_{h})_{h}-\left((c_{p}-c_{v})\rho_{0}\theta_{h}^{n},\nabla\cdot\hbox{\boldmath$\varphi$}_{h}\right)_{h}+(\sigma_{h}^{n},\nabla\hbox{\boldmath$\varphi$}_{h})_{h}=$
$\displaystyle\;0,\quad\forall\varphi_{h}\in\hbox{\boldmath$V$}_{h}^{k},$
(50b) $\displaystyle D_{\delta
t}^{m}({\sf\theta}_{j}^{\ell,n})+(\gamma-1)\nabla\cdot\hbox{\boldmath$u$}_{h}^{n}(\hbox{\boldmath$\xi$}_{j}^{\ell}){\sf\theta}_{j}^{\ell,n}-\frac{\sigma_{h}^{n}(\hbox{\boldmath$\xi$}_{j}^{\ell}):\nabla\hbox{\boldmath$u$}_{h}^{n}(\hbox{\boldmath$\xi$}_{j}^{\ell})}{c_{v}\rho_{0}(\hbox{\boldmath$\xi$}_{j}^{\ell})}=$
$\displaystyle\;0,\quad\forall j,\ell,$ (50c) where $D_{\delta t}^{m}(\alpha)$
is the approximation to the time derivative term $\alpha(t)^{\prime}$ using
BDF$[m]$, e.g., $D_{\delta
t}^{2}(\alpha^{n})=\frac{3\alpha^{n}-4\alpha^{n-1}+\alpha^{n-2}}{2\,\delta
t},\quad D_{\delta
t}^{3}(\alpha^{n})=\frac{11\alpha^{n}-18\alpha^{n-1}+9\alpha^{n-2}-2\alpha^{n-3}}{6\,\delta
t}.$

$\bullet$ The density, pressure, internal energy, and entropy approximations
$\rho_{h},p_{h},e_{h},s_{h}\in W_{h}^{k}$ are then recovered through (30) and
(31).

##  5\. Numerical results

We present numerical results in this section using the open-source finite-
element software NGSolve [31], https://ngsolve.org/. For all the simulation
results, we consider invisid models with zero physical viscosities.

We use a (variable time step size) BDF2 time stepping with high-order spatial
discretizations in Algorithm 3 for all the examples, except for Example 5.1
where higher order BDF time steppings are also used to verify the space/time
high-order accuracy of the proposed methods. For problems with shocks, we take
$q_{1}=0.5$ and $q_{2}=2$ as the default choice of the artificial viscosity
parameters in (46) unless otherwise stated.

We take the time step size as

$\displaystyle\delta
t=\min\\{\text{CFL}\frac{h_{min}}{|\hbox{\boldmath$u$}_{h}|+c_{s}}\\},$ (51)

where the length scale $h_{min}=h_{0}\alpha_{0}/k$ with $\alpha_{0}$ being the
minimal singular value of the Jacobian matrix
$\nabla_{X}\hbox{\boldmath$x$}_{h}$. Here the default choice of the CFL
constant is taken to be $\text{CFL}=1$ unless otherwise stated. The automatic
time-step control detailed in [10, 7.3] is also used for the examples with
shocks. The Newton’s method is used to solve the nonlinear system for velocity
DOFs in each time step, where the average iteration counts for all cases are
observed to be around 4-8. Most of the linearized systems in each Newton
iteration are solved using the sparse Cholesky factorization, with the only
exception of the low-Mach number cases in Example 5.6 where a direct pardiso
solver is used as the matrix failed to be positive definite therein.

###  5.1. Accuracy test: 2D Taylor-Green Vortex

We consider the invisid 2D Taylor-Green Vortex problem proposed in [10] to
check the high-order convergence of our proposed algorithm on deforming
domains. Following [10], we take $\gamma=5/3$ and add a source term

$e_{\sf src}(\hbox{\boldmath$x$},t)=\frac{3\pi}{8}(\cos(3\pi x)\cos(\pi
y)-\cos(\pi x)\cos(3\pi y))$

to the internal energy equation (20). The computational domain is a unit
square with wall boundary conditions, and the initial conditions are taken
such that the exact solutions are:

$\displaystyle\rho(\hbox{\boldmath$x$},t)=$ $\displaystyle\;1,$
$\displaystyle\hbox{\boldmath$u$}(\hbox{\boldmath$x$},t)=$
$\displaystyle\;(\sin(\pi x)\cos(\pi y),-\cos(\pi x)\sin(\pi y)),$
$\displaystyle p(\hbox{\boldmath$x$},t)=$
$\displaystyle\;\frac{1}{4}(\cos(2\pi x)+\cos(2\pi y))+1.$

We turn off artificial viscosity in the numerical simulations, and perform
mesh convergence studies for the $L^{2}$-errors of velocity and internal
energy at final time $t=0.5$ on a sequence of four consecutive uniform
rectangular meshes with size $2^{3+l}\times 2^{3+l}$ for $l=0,1,2,3$. We
consider the BDF$[m]$ scheme Algorithm 3 with polynomials of degree $m$ used
for the spatial discretization for $m=1,2,3,4$. The time step size is taken be
$\delta t=0.05/2^{l}$ for $l=0,1,2,3$. The midpoint rule Algorithm 2 (with
smaller time steps) is used to generate the starting values for high order BDF
schemes. History of convergence for the two $L^{2}$ errors at $t=0.5$ are
recorded in Table 2. We observe $m$-th order convergence for both variables
for all $1\leq m\leq 4$. Hence, spatial and temporal high order convergence is
achieved. This convergence rate is optimal for the BDF time stepping, but
suboptimal by one order for the spatial discretization. We find the
convergence behavior of the BDF$[m+1]$-$P^{m}$ scheme (not reported here for
simplicity) is similar to BDF$[m]$-$P^{m}$ for $1\leq m\leq 3$, which suggests
the one-order spatial convergence rate reduction is unavoidable for our
scheme.

We remark that our spatial discretization can be made slightly more efficient
by taking polynomial degree one order lower for the thermodynamic variables
than the flow map approximation, while still maintaining a similar convergence
behavior. However, since the major computational cost for our scheme is in the
nonlinear system solver in (47b). Such efficiency gain is not significant, and
will not be investigated further in this work.

$m$ | mesh | $\|\hbox{\boldmath$u$}_{h}-\hbox{\boldmath$u$}\|_{\Omega}$ | order | $\|e_{h}-e\|_{\Omega}$ | order
---|---|---|---|---|---
| $8\times 8$ | 1.052e-01 | – | 1.337e-01 | –
| $16\times 16$ | 4.131e-02 | 1.35 | 7.284e-02 | 0.88
1 | $32\times 32$ | 1.949e-02 | 1.08 | 3.697e-02 | 0.98
| $64\times 64$ | 9.710e-03 | 1.00 | 1.852e-02 | 1.00
| $8\times 8$ | 1.077e-02 | – | 1.263e-02 | –
| $16\times 16$ | 3.250e-03 | 1.73 | 3.578e-03 | 1.82
2 | $32\times 32$ | 7.933e-04 | 2.03 | 8.764e-04 | 2.03
| $64\times 64$ | 1.964e-04 | 2.01 | 2.192e-04 | 2.00
| $8\times 8$ | 5.590e-03 | – | 2.947e-03 | –
| $16\times 16$ | 5.809e-04 | 3.27 | 3.175e-04 | 3.21
3 | $32\times 32$ | 7.070e-05 | 3.04 | 4.181e-05 | 2.92
| $64\times 64$ | 8.766e-06 | 3.01 | 5.404e-06 | 2.95
| $8\times 8$ | 3.986e-03 | – | 2.510e-03 | –
| $16\times 16$ | 1.361e-04 | 4.87 | 7.305e-05 | 5.10
4 | $32\times 32$ | 5.690e-06 | 4.58 | 2.928e-06 | 4.64
| $64\times 64$ | 4.755e-07 | 3.58 | 3.304e-07 | 3.15
Table 2. Example 5.1: History of convergence for the BDF$[m]$-$P^{m}$ scheme
for $1\leq m\leq 4$.

###  5.2. 1D Shock Tube

We consider a simple 1D Riemann problem, the Sod shock tube on the domain
$\Omega=[-5,5]$ with initial condition

$(\rho,\hbox{\boldmath$u$},p)=(1,0,1)\text{ for
}x\in[-5,0],\quad(\rho,\hbox{\boldmath$u$},p)=(0.125,0,0.1)\text{ for
}x\in[0,5].$

Here $\gamma=1.4$. We apply the BDF2 scheme with polynomial degree $k=4$ on an
initial mesh with 20 uniform cells. The results at final time $t=2.0$ on all
$20\times(4+1)=100$ quadrature points are shown in Figure 1, along with the
deformed cells. The numerical approximation agrees with the exact Riemann
solution quite well even on this coarse mesh, where the shock is resolved
within 2 cells. As typical of Lagrangian schemes, the contact discontinuity is
captured without dissipation. We also observe the “wall heating” phenomenon in
the internal energy at the contact.

Figure 1. Example 5.2: Results at final time $t=0.2$ sampled on 100 quadrature
points. The gray lines are the deformed cell boundaries.

###  5.3. 2D Sedov explosion

The Sedov explosion [32] models the expanding wave by an intense explosion in
a perfect gas. It is a standard problem to test the ability of codes to
preserve the radial symmetry of shocks. The domain is a square
$\Omega=[0,1.2]\times[0,1.2]$. The initial condition is set to have unit
density and zero velocity, and also to have zero internal energy except at the
left bottom corner cell $T_{0}$, where its is a bilinear function whos value
is $\frac{0.2448\times 4}{\mathrm{area}(T_{0})}$ on the left/bottom corner
vertex and zero on the other three vertices. So the total initial internal
energy is $\int_{\Omega_{0}}\rho e\mathrm{dX}=0.2448.$ Symmetry boundary
conditions are imposed on the left and bottom boundaries, while free boundary
conditions are used on the top and right boundaries. The analytic solution at
time $t=1$ gives a shock at radius $r=1$ with a peak density of $6$.

We apply the BDF2 scheme with polynomial degree $k=4$ on initial uniform
rectangular meshes of size $N\times N$ with $N=16$ and $N=32$. The density
field on deformed meshes, and also the scattered plot of density versus radius
$r=\sqrt{x^{2}+y^{2}}$ on all quadrature points are shown in Figure 2. The
radial symmetry of the solution is preserved, and the numerical solution
agrees quite well with the analytic solution in Figure 2(c).

(a) $k=4$, $N=16$

(b) $k=4$, $N=32$

(c) den. v.s. rad.

Figure 2. Example 5.3: (a)-(b): Density contour on deformed domain at final
time $t=1.0$ using BDF2 time stepping with polynomial degree $k=4$ on
rectangular meshes with size $N\times N$. (c): Scattered plot of density v.s.
radius on all quadrature points at final time $t=1$.

###  5.4. 2D Noh explosion

The Noh explosion problem [25] consists of an ideal gas with $\gamma=5/3$,
initial density $\rho_{0}=1$, initial internal energy $e_{0}=0$, and initial
velocity
$\hbox{\boldmath$u$}_{0}=(-\frac{x}{\sqrt{x^{2}+y^{2}}},-\frac{y}{\sqrt{x^{2}+y^{2}}})$.
The computational domain is $\Omega=[0,1]\times[0,1]$. We use symmetry
boundary conditions on left and bottom boundaries, and free boundary
conditions on top and right boundaries. Similar to the previous case, this
problem has a radial symmetry, and the analytic solution at time $t=0.6$ gives
a shock at radius $r=0.2$ with a peak density of $16$.

We apply the BDF2 scheme with polynomial degree $k=4$ on initial uniform
rectangular meshes of size $N\times N$ with $N=16$ and $N=32$. For this
problem, we observe the default choice of artificial viscosity coefficients
with $q_{1}=0.5$ and $q_{2}=2$ leads to quite large post-shock oscillations.
So we increase these coefficients to $q_{1}=1$ and $q_{2}=4$ in our numerical
experiments reported here. The radial symmetry of the solution is preserved as
seen in (a) and (b) of Figure 3, and the numerical solution has a good
agreement with the analytic solution for $0.1<r<0.3$ in Figure 2(c), although
it is slightly oscillatory.

(a) $k=4$, $N=16$

(b) $k=4$, $N=32$

(c) den. v.s. rad.

Figure 3. Example 5.4: (a)-(b): Density contour on deformed domain at final
time $t=0.6$ using BDF2 time stepping with polynomial degree $k=4$ on
rectangular meshes with size $N\times N$. Here the gray square is the initial
domain. (c): Scattered plot of density v.s. radius on all quadrature points at
final time $t=0.6$.

###  5.5. Triple point problem

The triple point problem is a multimaterial test case proposed in [17]; see
also [13]. The initial data is shown in Fig. 4. The computational domain has a
rectangular shape with $7\times 3$ edge ratio. It includes three materials at
rest located in $\Omega_{1}=[0,1]\times[0,3]$,
$\Omega_{2}=[1,7]\times[0,1.5]$, and $\Omega_{3}=[1,7]\times[1.5,3]$,
initially forming a T-junction. The high-pressure material in $\Omega_{1}$
creates a shock wave moving to the right. Due to different material properties
in $\Omega_{2}$ and $\Omega_{3}$, the shock wave moves faster in $\Omega_{3}$,
which leads to vortex formation around the triple point. For pure Lagrangian
methods, there is a limit to how long this problem can be run due to vortex
generation. Similar to [10], we run simulation till time $t=3.3$.

Figure 4. Example 5.5: Initial data for the triple point problem.

We apply the BDF2 scheme with polynomial degree $k=4$ and on initial uniform
rectangular meshes of size $28\times 12$ and $56\times 24$. Here we reduce the
artificial viscosity coefficients to be $q_{1}=0.25$ and $q_{2}=1$, and take a
larger CFL number with $\text{CFL}=3$. A total of 238 time steps is used to
drive the solution to final time $t=3.3$ on the fine mesh with $56\times 28$
cells. We note that an explicit scheme would require about two order of
magnitude more time steps for a stable simulation. For example, the Laghos
code freely available in the github repository https://github.com/CEED/Laghos,
which implements the high-order Lagrangian finite element scheme in [10],
requires about $90,000$ time steps for $k=4$ on the fine mesh with RK4 time
stepping and CFL=1. We further note that the default choice of parameters also
work for this problem; we make the modifications to illustrate that the high-
order method is still robust with a larger time step size and a smaller
artificial viscosity.

The density plot in log-scale on the deformed meshes are shown in Figure 5. We
observe the material interfaces are sharply preserved as typical of Lagrangian
schemes, and the the shock locations are essentially the same, but the total
amount of ”roll-up” at the triple point increases as the mesh refines.

Figure 5. Example 5.5: Results at final time $t=3.3$.

###  5.6. Gresho vortex

The Gresho vortex problem [15] is an example of a stationary, incompressible
rotating flow around the origin in two spatial dimensions, where centrifugal
forces are exactly balanced by pressure gradients. It was first applied to the
compressible Euler equations in [18]. Here we use the low-Mach setup given in
[23]. The domain is a unit square $\Omega=(0,1)\times(0,1)$ with wall boundary
conditions, $\gamma=1.4$, and the initial conditions are

$\displaystyle\rho_{0}=1,\hbox{\boldmath$u$}_{0}=\;u_{\phi}\hbox{\boldmath$e$}_{\phi},p_{0}=\;\frac{1}{\gamma
M_{\max}^{2}}-\frac{1}{2}+$ $\displaystyle\begin{cases}12.5r^{2}&\text{if
}r<0.2,\\\ 4\mathrm{ln}(5r)+4-20r+12.5r^{2}&\text{if }0.2\geq r\geq 0.4,\\\
4\mathrm{ln}2-2&\text{if }r>0.4,\end{cases}$

where the radius $r=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}$, angular velocity
$u_{\phi}$ is

$u_{\phi}(r)=\begin{cases}5r&\text{if }r<0.2,\\\ 2-5r&\text{if }0.2\geq r\geq
0.4,\\\ 0&\text{if }r>0.4,\\\ \end{cases}$

the unit vector $\hbox{\boldmath$e$}_{\phi}=(-(y-0.5)/r,(x-0.5)/r)$, and
$M_{\max}$ is the parameter used to adjust the maximum Mach number of the
problem, where $M(r)=\frac{|\hbox{\boldmath$u$}|}{\sqrt{\gamma p/\rho}}$ is
the Mach number. The maximum number of the initial condition is achieved at
$r=0.2$ with $M(0.2)=M_{\max}$.

We apply the BDF2 scheme with polynomial degree $k=2$ on a $40\times 40$ mesh,
and $k=4$ on a $20\times 20$ mesh. The total number of velocity DOFs are the
same for the two cases. Since the problem is smooth, we turn off the
artificial viscosity. The period of one rotation for $r=0.2$ is $2\pi
r=0.4\pi$. We run simulation till final time $t=\frac{3}{4}\times{0.4\pi}$ so
that the internal flow has rotated $\frac{3}{4}\times 180=135$ degrees. The
default choice of time step size (51) is linear proportional to the Mach
number since $\max c_{s}\approx 1/M_{\max}$. We remove this Mach-number
dependency on time step size by multiplying the sound speed in (51) with the
maximum Mach number, i.e.,

$\delta
t=\min\\{\text{CFL}\frac{h_{min}}{|\hbox{\boldmath$u$}_{h}|+M_{\max}c_{s}}\\}.$

We use $M_{\max}=0.1,0.01,0.001$, and take the CFL number to be CFL=0.25 for
all cases. The relative Mach number $M/M_{\max}$ at final time
$t=\frac{3}{4}\times 0.4\pi$ for all cases are shown in Figure 6. The first
row of Figure 6 show results for the $k=2$ simulations, where we clearly
observe a locking phenomena as the Mach number decreases. On the other hand,
the higher order simulations with $k=4$ leads to almost identical results for
all three Mach numbers. This example illustrates the advantage of using a
higher order scheme over a low-order scheme in the low-Mach number regime.

Moreover, the total number of time steps for $k=4$ are between 900 to 1000 for
all three cases. If an explicit scheme (e.g., RK3) were to be used to solve
this problem, the total number of time steps will be about three orders of
magnitude larger when $M_{\max}=0.001$ due to sound speed based CFL
constraints. Finally, we note that the nonlinear system in each time step is
harder to solve as the Mach number $M_{\max}$ decreases. For example, Newton’s
method with sparse Cholesky direct solver works for $M_{\max}=0.1$ with CFL=1,
but it fails for $M_{\max}=0.001$, where we have to reduce CFL to be 0.25 and
replace the Cholesky solver by a pardiso solver, which indicates the linear
system for $M_{\max}=0.001$ in the Newton iteration is no longer positive
definite. The linear system solver issue in the low Mach number regime will be
further investigate in our future work.

(a) $k=2$, $N=40$, $M_{\max}=0.1$

(b) $k=2$, $N=40$, $M_{\max}=0.01$

(c) $k=2$, $N=40$, $M_{\max}=0.001$

(d) $k=4$, $N=20$, $M_{\max}=0.1$

(e) $k=4$, $N=20$, $M_{\max}=0.01$

(f) $k=4$, $N=20$, $M_{\max}=0.001$

Figure 6. Example 5.6: Density contour on deformed domain at final time
$t=1.0$ using BDF2 time stepping and Algorithm 1 with polynomial degree $k$ on
a rectangular mesh with size $N\times N$ for different choices of $k$ and $N$.
(parameters: CFL = 1.0, $q_{1}=0.5$, $q_{2}=2.0$)

###  5.7. Shock-bubble interaction

This test case corresponds to the interaction of shock wave with a cylindrical
Helium bubble surrounded by air at rest [28]. We use the same steup as in [13,
Section 8.4]. The initial domain is a rectangular box
$\Omega_{0}=(0,L)\times(-H/2,H/2)=(0,0.650)\times(-0.089,0.089)$, which
includes a circular bubble with center $(0.320,0)$ and radius $r_{b}=0.025$.
Initial data are shown in Figure 7(a). Wall boundary conditions at each
boundary is prescribed except at the right boundary, where we impose a piston-
like boundary condition with inward velocity $(-124.824,0)$. The left going
shock wave hits the bubble at time $t_{i}=668.153\times 10^{-6}$. The final
time of simulation is $t_{f}=t_{i}+674\times 10^{-6}=1342.153\times 10^{-6}$,
which corresponds to the time where experimental shadow-graph extracted from
[16] is displayed in [28].

We use two unstructured triangular meshes that is fitted to the bubble
boundary for this problem. The coarse mesh has 8324 triangular cells whose
mesh size is $h=H/32$, while the fine mesh has 33526 cells whose mesh size is
$h=H/64$; see Figure 7(b)-(c) for the zoom-in view of the two meshes around
the bubble. The BDF2 time stepping is used in combination with polynomial
degree $k=2$ and $k=4$ on these two meshes. The zoom-in view around the
deformed bubble at final time $t_{f}$ are shown in Figure 8. We observe that
the location and shape of the deformed bubble is similar for each simulation,
with a better resolution being obtained on a finer mesh with a higher order
polynomial degree. These shapes are also qualitatively similar to the
experimental Schlieren image in Figure 8(e) obtained from [16].

We display in Figure 9 the time evolution of the bubble at times $t=800\times
10^{-6},1100\times 10^{-6},1342.153\times 10^{-6}$ for $k=4$ on the two
meshes. We note that the results obtained with both meshes are quite similar.

(a) Geometry

(b) Coarse mesh

(c) Fine mesh

Figure 7. Example 5.7: Geometry setup and zoomed-in meshes around the bubble.

(a) $k$=2, coarse mesh

(b) $k$=2, fine mesh

(c) $k$=4, coarse mesh

(d) $k$=4, fine mesh

(e) Data from [16]

Figure 8. Example 5.7: Zoom on the deformed bubble at final time $t$=
1342.153e-6. Here (e) is the Schlieren image from experimental data [16].

(a) $k=4$, coarse mesh, $t$=800e-6

(b) $k=4$, fine mesh, $t$=800e-6

(c) $k=4$, coarse mesh, $t$=1100e-6

(d) $k=4$, fine mesh, $t$=1100e-6

(e) $k=4$, coarse mesh, $t$=1342.153e-6

(f) $k=4$, fine mesh, $t$=1342.153e-6

Figure 9. Example 5.7: Density contour in the vicinity of the bubble at
various times.

###  5.8. Multimaterial implosion in cylindrical geometry

The aim of this example is to access the capability of the implicit high-order
Lagrangian scheme to handle a multi-mode implosion in cylindrical geometry.
Here we consider a simple 1D multimaterial implosion problem on unstructured
2D triangular meshes. The problem consists of a low-density material with
$\rho_{1}=0.05$ in the radial range $r\in[0,1]$ surrounded by a shell of high
density material with $\rho_{2}=1.0$ in the radial range $r\in(1.0,1.2]$. Both
materials are initially at rest with pressure $p=0.1$ and adiabatic index
$\gamma=5/3$. This problem was originally proposed in [13] with a time
dependent pressure source on the outer radial surface $r=1.2$. Here we use the
modification in [10] that applies a constant radial velocity source of
$\hbox{\boldmath$u$}=-5(x,y)/1.2$ on the outer boundary, which drives a
cylindrical shock wave inwards. Due to the 1D setup (symmetry) , the material
interface should be a function of radius only for all time. The movement of
the interface radius against time is shown on the left panel of Figure 11. We
clearly observe the deceleration of the interface starting around $t=0.12$,
and the so-called stagnation phase is reached around $t=0.14$ where the radius
obtains its minimum value. The flow becomes Rayleigh-Taylor unstable after
this time as a small perturbation of the interface will grow exponentially as
function of time due to the fact that the light fluid inside is pushing the
heavy fluid outside after the stagnation phase. It is very challenging to
preserve the interface symmetry for a Lagrangian scheme on general
unstructured meshes, especially when time past the stagnation phase.

Due to symmetry of the problem, we take the computational domain to be a
quarter circle $\Omega=\\{(x,y):x\geq 0,y\geq 0,x^{2}+y^{2}\leq 1.2^{2}\\}$
with wall boundary conditions on the left and bottom boundaries. We apply the
BDF2 scheme with polynomial degree $k=2,4$ on two set of unstructured
triangular meshes. The coarse mesh has mesh size $h=0.05$ and is used for
$k=4$, while the fine mesh is a uniform refinement of the coarse mesh and is
used for $k=2$; see Figure 10(d) for the coarse mesh and Figure 10(a) for the
fine mesh. We take $q_{1}=1$ and $q_{2}=4$ as the artificial viscosity
coefficients and set CFL number equals 0.25. The two simulations have the same
number of velocity DOFs, and requires a similar amount of total time steps to
reach the final time $t=0.16$. Density plots in log-scale on the deformed
meshes are shown in Figure 10 at times $t=0.08$ and $t=0.16$, along with the
initial density at $t=0$. We observe while the results for $t=0.08$ are
similar for both cases, where radial symmetry of the material interface is
preserved quite well. The results for $k=4$ at $t=0.16$ (which past the
stagnation phase) is better than that for $k=2$ in terms of the interface
symmetry. In Figure 11 we plot the time evolution of the average radius of the
material interface, and the normalized standard deviation of this radius at
different times as an indication of the symmetry error over time. We observe
similar results for the average radius for both simulations, and a smaller
symmetry error for the higher order case.

(a) $k=2$, fine mesh, $t$=0

(b) $k=2$, fine mesh, $t$=0.08

(c) $k=2$, fine mesh, $t$=0.16

(d) $k=4$, coarse mesh, $t$=0

(e) $k=4$, coarse mesh, $t$=0.08

(f) $k=4$, coarse mesh, $t$=0.16

Figure 10. Example 5.7: Density contour in the vicinity of the bubble at
various times.

Figure 11. Example 5.7: Density contour in the vicinity of the bubble at
various times.

##  6\. Conclusion

We presented a class of high-order variational Lagrangian schemes for
compressible flow. The discrete EnVarA approach is used to derive the high-
order spatial finite element discretization. Features of our spatial
discretization include mass/momentum/energy conservation and entropy
stability. Fully implicit time stepping is then applied to the resulting ODE
system. Each time step requires a nonlinear system solve for the velocity DOFs
only. Ample numerical results are shown to support the good performance of the
proposed scheme. We plan to extend this pure Lagrangian scheme to the
arbitrary Eulerian Lagrangian framework, which has the potential to address
the mesh distortion issue.

## Appendix

Here we briefly discuss the isothermal case [14] where the temperature does
not change over time. In this case, the free energy (2) is a function of
density only: $\psi=\psi(\rho)$. Typical choices include
$\psi(\rho)=\alpha\rho^{\gamma}$ or $\psi(\rho)=\alpha\rho\log(\rho)$.

This model only has two thermodynamic variables: density $\rho$ and pressure

$p=\psi_{\rho}\rho-\psi.$

The EnVarA derivation leads to the following model equations; see also [14]:

$\displaystyle\dot{\hbox{\boldmath$x$}}=$
$\displaystyle\hbox{\boldmath$u$},\quad\quad\dot{(\rho J)}=\;0,$ (52a)
$\displaystyle\rho\dot{\hbox{\boldmath$u$}}=$ $\displaystyle\;-\nabla
p+\nabla\cdot\left(\eta\nabla_{s}\hbox{\boldmath$u$}+(\xi-\frac{2}{3}\eta)(\nabla\cdot\hbox{\boldmath$u$})\hbox{\boldmath$I$}\right),$
(52b)

where $p=\psi_{\rho}\rho-\psi$.

The spatial discretization for the model (52) is summarized below:

Algorithm 4 Variational Lagrangian scheme for the model (52).

Find
$\hbox{\boldmath$x$}_{h},\hbox{\boldmath$u$}_{h}\in\hbox{\boldmath$V$}_{h}^{k}$
such that  $\displaystyle\dot{\hbox{\boldmath$x$}}_{h}=$
$\displaystyle\;{\hbox{\boldmath$u$}}_{h},$ (53a)
$\displaystyle(\rho_{0}\dot{\hbox{\boldmath$u$}}_{h},\hbox{\boldmath$\varphi$}_{h})_{h}-\left(p_{h},\nabla\cdot\hbox{\boldmath$\varphi$}_{h}\right)_{h}+(\sigma_{h},\nabla\hbox{\boldmath$\varphi$}_{h})_{h}=$
$\displaystyle\;0,\quad\forall\varphi_{h}\in\hbox{\boldmath$V$}_{h}^{k},$
(53b) where the pressure $p_{h}\in W_{h}^{k}$ is a function of the density
approximation with
$p_{i}^{\ell}=\psi_{\rho}(\rho_{i}^{\ell})\rho_{i}^{\ell}-\psi(\rho_{i}^{\ell}),$
in which $\rho_{h}\in W_{h}^{k}$ satisfies mass conservation (30), and the
stress
$\sigma_{h}=(\eta+\mu_{av})J_{h}\nabla_{s}{\hbox{\boldmath$u$}}_{h}+(\xi-\frac{2}{3}\eta)J_{h}\nabla\cdot{\hbox{\boldmath$u$}}_{h}\hbox{\boldmath$I$},$
in which the artificial viscosity coefficient $\mu_{av}$ is given in (46).

It is clear that the spatial discretization in Algorithm 4 is mass and
momentum conservative. Next we prove Algorithm 4 is also entropy stable. Using
a similar derivation as in (9) and (36) and the definition of density and
pressure, we obtain

$\frac{d}{dt}(\psi(\rho_{h})J_{h},1)_{h}=-(p_{h}J_{h},\nabla\cdot\hbox{\boldmath$u$}_{h})_{h}.$

Taking test function $\hbox{\boldmath$\varphi$}_{h}=\hbox{\boldmath$u$}_{h}$
in (53b) and use the above relation, we get

$\frac{d}{dt}\underbrace{\left(\frac{1}{2}\rho_{0}|\hbox{\boldmath$u$}_{h}|^{2}+\psi(\rho_{h})J_{h},1\right)_{h}}_{\text{total
entropy}}=-(\sigma_{h},\nabla\hbox{\boldmath$u$}_{h})\leq 0.$

The ODE system in Algorithm 4 can be discretized using the implicit schemes
discussed in Section 4. Here we present a variational implicit time
discretization similar to the backward Euler scheme in Section 4.1. Given data
$\hbox{\boldmath$x$}^{n-1},\hbox{\boldmath$u$}_{h}^{n-1}\in V_{h}^{k}$ at time
$t^{n-1}$ and time step size $\delta t$, denote the following discrete energy
functional for the flow map $\hbox{\boldmath$x$}_{h}$:

$\displaystyle
E_{h}(\hbox{\boldmath$x$}_{h}):=\left(\frac{\rho_{0}|\hbox{\boldmath$x$}_{h}-\hbox{\boldmath$x$}_{h}^{n-1}-\delta
t\hbox{\boldmath$u$}_{h}^{n-1}|^{2}}{2\delta
t^{2}}+\psi\left(\rho_{0}/J_{h}\right)J_{h},1\right)_{h}+\frac{1}{2}(\Delta_{h},1)_{h},$
(54)

where the dissipation

$\Delta_{h}:=(\eta+\mu_{av}^{n-1})J_{h}|\nabla_{s}\frac{{\hbox{\boldmath$x$}}_{h}-{\hbox{\boldmath$x$}}_{h}^{n-1}}{\delta
t}|^{2}+(\xi-\frac{2}{3}\eta)J_{h}|\nabla\cdot\frac{{\hbox{\boldmath$x$}}_{h}-{\hbox{\boldmath$x$}}_{h}^{n-1}}{\delta
t}|^{2},$

with artificial viscosity $\mu_{av}^{n-1}$ explicitly evaluated at time level
$t^{n-1}$, and $J_{h}=|\nabla_{X}\hbox{\boldmath$x$}_{h}|$. The flow map at
next time level is obtained by solving the following minimization problem:

$\displaystyle\hbox{\boldmath$x$}_{h}^{n}:=\mathrm{argmin}_{\hbox{\boldmath$x$}_{h}\in
V_{h}^{k},J_{h}>0}E_{h}(\hbox{\boldmath$x$}_{h}).$ (55)

Assuming piecewise constant viscosity coefficients $\eta,\xi$, and
$\mu_{av}^{n-1}$, the Euler-Lagrangian equation for this minimization problem
is simply the Backward Euler scheme in (47) applied to the ODE system (53). We
note that such energy minimization interpretation is not available for the
non-isothermal case discussed in Section 4 due to the temperature equation
(47c).

## References

  * [1] R. Baierlein, Thermal Physics, Cambridge University Press, 1999.
  * [2] R. S. Berry, S. A. Rice, and J. Ross, Physical Chemistry, Oxford University Press, Cambridge, 2000.
  * [3] G. A. Bird, Molecular Gas Dynamics And The Direct Simulation Of Gas Flows, Clarendon Press, Oxford, 1994.
  * [4] J. Burkardt and C. Trenchea, Refactorization of the midpoint rule, Appl. Math. Lett., 107 (2020), pp. 106438, 7.
  * [5] E. J. Caramana, D. E. Burton, M. J. Shashkov, and P. P. Whalen, The construction of compatible hydrodynamics algorithms utilizing conservation of total energy, J. Comput. Phys., 146 (1998), pp. 227–262.
  * [6] J. Cheng and C.-W. Shu, A high order ENO conservative Lagrangian type scheme for the compressible Euler equations, J. Comput. Phys., 227 (2007), pp. 1567–1596.
  * [7] C. M. Dafermos, The second law of thermodynamics and stability, Arch. Rational Mech. Anal., 70 (1979), pp. 167–179.
  * [8] C. M. Dafermos, Hyperbolic conservation laws in continuum physics, vol. 325 of Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], Springer-Verlag, Berlin, fourth ed., 2016\.
  * [9] B. Després and C. Mazeran, Lagrangian gas dynamics in two dimensions and Lagrangian systems, Arch. Ration. Mech. Anal., 178 (2005), pp. 327–372.
  * [10] V. A. Dobrev, T. V. Kolev, and R. N. Rieben, High-order curvilinear finite element methods for Lagrangian hydrodynamics, SIAM J. Sci. Comput., 34 (2012), pp. B606–B641.
  * [11] B. Eisenberg, Y. Hyon, and C. Liu, Energy variational analysis of ions in water and channels: Field theory for primitive models of complex ionic fluids, The Journal of Chemical Physics, 133 (2010), p. 104104.
  * [12] J. L. Ericksen, Introduction to the thermodynamics of solids, vol. 131 of Applied Mathematical Sciences, Springer-Verlag, New York, revised ed., 1998.
  * [13] S. Galera, P.-H. Maire, and J. Breil, A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction, J. Comput. Phys., 229 (2010), pp. 5755–5787.
  * [14] M.-H. Giga, A. Kirshtein, and C. Liu, Variational modeling and complex fluids, in Handbook of Mathematical Analysis in Mechanics of Viscous Fluids, Y. Giga and A. Novotny, eds., Springer International Publishing, 2017, pp. 1–41.
  * [15] P. M. Gresho and S. T. Chan, On the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. I. Theory, Internat. J. Numer. Methods Fluids, 11 (1990), pp. 587–620.
  * [16] J. Haas and B. Sturtevant, Interaction of weak-shock waves, J. Fluid Mech., 181 (1987), pp. 41–76.
  * [17] M. Kucharik, R. V. Garimella, S. P. Schofield, and M. J. Shashkov, A comparative study of interface reconstruction methods for multi-material ALE simulations, J. Comput. Phys., 229 (2010), pp. 2432–2452.
  * [18] R. Liska and B. Wendroff, Comparison of several difference schemes on 1D and 2D test problems for the Euler equations, SIAM J. Sci. Comput., 25 (2003), pp. 995–1017.
  * [19] C. Liu, An introduction of elastic complex fluids: an energetic variational approach, in Multi-Scale Phenomena in Complex Fluids: Modeling, Analysis and Numerical Simulation, World Scientific, 2009, pp. 286–337.
  * [20] C. Liu and J.-E. Sulzbach, The Brinkman-Fourier system with ideal gas equilibrium, Discrete and Continuous Dynamical Systems, 42 (2022), pp. 425–462.
  * [21] P.-H. Maire, R. Abgrall, J. Breil, and J. Ovadia, A cell-centered Lagrangian scheme for two-dimensional compressible flow problems, SIAM J. Sci. Comput., 29 (2007), pp. 1781–1824.
  * [22] D. A. McQuarrie, Statistical Mechanics, Harper & Row, New York, 1976\.
  * [23] F. Miczek, F. Röpke, and P. Edelmann, New numerical solver for flows at various mach numbers, Astronomy & Astrophysics, 576 (2015), p. A50.
  * [24] C.-D. Munz, On Godunov-type schemes for Lagrangian gas dynamics, SIAM J. Numer. Anal., 31 (1994), pp. 17–42.
  * [25] W. Noh, Errors for calculations of strong shocks using an artificial viscosity and an artificial heat flux, J. Comput. Phys., 72 (1987), pp. 78–120.
  * [26] L. Onsager, Reciprocal relations in irreversible processes. i., Physical review, 37 (1931), p. 405.
  * [27] , Reciprocal relations in irreversible processes. ii., Physical review, 38 (1931), p. 2265.
  * [28] J. Quirk and S. Karni, On the dynamics of a shock-bubble interaction, J. Fluid Mech., 318 (1996), pp. 129–163.
  * [29] S. Salinas, Introduction to Statistical Physics, Springer, New York, 2001.
  * [30] A. Sandu, V. Tomov, L. Cervena, and T. Kolev, Conservative high-order time integration for Lagrangian hydrodynamics, SIAM J. Sci. Comput., 43 (2021), pp. A221–A241.
  * [31] J. Schöberl, C++11 Implementation of Finite Elements in NGSolve, 2014.  ASC Report 30/2014, Institute for Analysis and Scientific Computing, Vienna University of Technology.
  * [32] L. I. Sedov, Similarity and dimensional methods in mechanics, “Mir”, Moscow, 1982.  Translated from the Russian by V. I. Kisin.
  * [33] W. Strutt, J, Some general theorems relating to vibrations, Proceedings of the London Mathematical Society, 1 (1871), pp. 357–368.
  * [34] H. Sun and C. Liu, On energetic variational approaches in modeling the nematic liquid crystal flows, Discrete and Continuous Dynamical Systems, 23 (2009), pp. 455–475.
  * [35] J. Von Neumann and R. D. Richtmyer, A method for the numerical calculation of hydrodynamic shocks, J. Appl. Phys., 21 (1950), pp. 232–237.
  * [36] M. L. Wilkins, Use of artificial viscosity in multidimensional fluid dynamic calculations, J. Comput. Phys., 36 (1980), pp. 281–303.
  * [37] F. D. Witherden and P. E. Vincent, On the identification of symmetric quadrature rules for finite element methods, Comput. Math. Appl., 69 (2015), pp. 1232–1241.
  * [38] L. Zhang, T. Cui, and H. Liu, A set of symmetric quadrature rules on triangles and tetrahedra, J. Comput. Math., 27 (2009), pp. 89–96.