High Order Semi-implicit WENO Schemes for All Mach Full Euler System of Gas
Dynamics

Sebastiano Boscarino 111Department of Mathematics and Computer Science,
University of Catania, Catania, 95127. E-mail<EMAIL_ADDRESS>S.
Boscarino and G. Russo are members of the INdAM Research group GNCS. They
would like to thank the Italian Ministry of Instruction, University and
Research (MIUR) to support this research with funds coming from PRIN Project
2017 (No. 2017KKJP4X, entitled “Innovative numerical methods for evolutionary
partial differential equations and applications”). G. Russo has also been
supported by ITN-ETN Horizon 2020 Project ModCompShock, Modeling and
Computation on Shocks and Interfaces, Project Reference 642768 and S.
Boscarino has been supported by the University of Catania (“Piano della
Ricerca 2016/2018, Linea di intervento 2”)., Jing-Mei Qiu222Department of
Mathematics, University of Delaware, Newark, DE 19716. E-mail:
<EMAIL_ADDRESS>Research is supported by NSF grant NSF-DMS-1818924, Air
Force Office of Scientific Research FA9550-18-1-0257. , Giovanni Russo
333Department of Mathematics and Computer Science, University of Catania,
Catania, 95125. E-mail<EMAIL_ADDRESS>Tao Xiong 444Corresponding
author. School of Mathematical Sciences, Fujian Provincial Key Laboratory of
Mathematical Modeling and High-Performance Scientific Computing, Xiamen
University, Xiamen, Fujian, P.R. China, 361005. Email<EMAIL_ADDRESS>Research is supported by NSFC No. 11971025, NSF of Fujian Province No.
2019J06002, the Strategic Priority Research Program of Chinese Academy of
Sciences Grant No. XDA25010401, the Science Challenge Project No. TZ2016002.

Abstract. In this paper, we propose high order semi-implicit schemes for the
all Mach full Euler equations of gas dynamics. Material waves are treated
explicitly, while acoustic waves are treated implicitly, thus avoiding severe
CFL restrictions for low Mach flows. High order accuracy in time is obtained
by semi-implicit temporal integrator based on the IMEX Runge-Kutta (IMEX-RK)
framework. High order in space is achieved by finite difference WENO schemes
with characteristic-wise reconstructions adapted to the semi-implicit IMEX-RK
time discretization. The schemes are proven to be asymptotic preserving and
asymptotically accurate as the Mach number vanishes. Besides, they can well
capture discontinuous solutions in the compressible regime, especially for two
dimensional Riemann problems. Numerical tests in one and two space dimensions
will illustrate the effectiveness of the proposed schemes.

Keywords. All Mach number, full Euler Equations, Asymptotic Preserving,
Asymptotically Accurate, finite difference WENO, characteristic-wise
reconstruction.

AMS subject classification. 35L65, 35B40, 35C20, 76N15, 76M20, 76M45, 76B47,
65M06.

##  1 Introduction

Computational fluid dynamics (CFD) has been a very active research field in
the past decades. Numerical methods developed in this area generally can be
divided into two categories, which are classified by the dimensionless Mach
number. For moderate to high Mach number compressible effects have to be taken
into account, while for low Mach number the flow can be considered
incompressible or weakly compressible. For compressible flows, most numerical
solutions are obtained by Godunov type shock capturing schemes for
compressible Euler equations, which have the structure of a hyperbolic system
of conservation laws [34, 44, 21, 40, 13], while for the incompressible flows,
preserving incompressibility and resolving vortex dynamics are among the main
purposes [11, 43, 22].

There are, however, circumstances in which flows with a wide range of Mach
number appear, making it desirable to develop numerical methods which can be
applied for fluid flows at any speed, as already shown in the pioneering work
of Harlow and Amsden [26, 27]. However, due to the different physical
mechanisms and mathematical characteristics for the governing equations at
different speeds, the development of efficient and effective numerical methods
to capture flows with different compressibility is challenging [32, 18], and a
lot or progress has been made only recently. For hyperbolic systems, waves
propagate at finite speeds. Numerical methods have to resolve all the space
and time scales that characterize these waves. Most shock capturing schemes
devoted to such systems are obtained by explicit time discretization, and the
time step has to satisfy a stability restriction, known as the CFL condition:
it is limited by the size of the spatial mesh divided by the fastest wave
speed. For compressible flows with Mach number greater than, say, $0.1$, such
a restriction is not a problem: indeed, if one is interested in resolving all
the waves, accuracy and stability restrictions on space and time
discretization are of similar nature. However, for low Mach flows, acoustic
waves usually carry a negligible amount of energy. If one is not interested in
resolving them, then the system becomes stiff: stability limitations on the
time step are much stricter than the restrictions imposed by accuracy [45,
46]. In such cases, one may resort to implicit time discretization to avoid
the acoustic CFL restriction. However, shock capturing schemes are highly non-
linear, and a naive implicit version of them risks to be very inefficient.
Furthermore, numerical viscosities for Godunov-type schemes are inversely
proportional to the Mach number, introducing excessive numerical dissipations
on the slow waves [17]. Preconditioning techniques are adopted to cure the
large numerical diffusion as discussed in [45, 47, 35], but such techniques
are effectively applicable only if Mach numbers are not too small.

On the other hand, as the Mach number vanishes, the flow converges to the
incompressible limit. For the full Euler equations, at the incompressible
limit, the density remains constant along the fluid particle trajectories, and
the pressure waves propagate with infinite speed. Mass conservation equation
reduces to the incompressibilty condition $\nabla\cdot{\bf u}=0$ on the
velocity field [12], so that the pressure and the density are decoupled. The
pressure turns out to act as a Lagrange multiplier to enforce
incompressibility of the flow [32]. A rigorous proof for the compressible flow
converging to the incompressible one as the Mach number goes to zero is given
in [30]. An effective approach to deal with low Mach flows is given by
pressure-based algorithms, such as, for example, the one by Casulli and
Greenspan [10], in which a semi-implicit treatment of the pressure is
incorporated in a scheme for compressible flow. The authors use an upwind
discretization on the material wave, and an implicit equation for the
pressure, which is solved by a SOR-type method. Several authors have
subsequently worked on the development of semi-implicit methods [32, 36] based
on low-Mach asymptotic [30]. However, many of such schemes are specifically
designed to deal with low Mach flows. When the fluid flow is compressible at
large speed, shock discontinuities may form and propagate. In these cases, it
is necessary to resort to conservative schemes (density-based schemes) which
correctly capture possible shocks.

Recently several papers have been written along these lines, see for example
[16, 23, 41, 20, 8, 7] for isentropic Euler and Navier-Stokes equations, or
[14, 42, 19, 48, 18] for full Euler and Navier-Stokes equations. However, most
finite volume and finite difference schemes for full Euler equations are
second order accurate in space and time, while existing high order schemes
developed for all Mach flows in isentropic Euler equations are not robust
enough to be directly extended to the full Euler equations.

The aim of the present paper is to design a high order finite difference shock
capturing scheme for the full compressible Euler equations. Finite difference
weighted essentially non-oscillatory (WENO) schemes are used for spatial
discretization, while high order semi-implicit IMEX Runge-Kutta (SI-IMEX-RK)
methods are adopted for time discretization. The stability and accuracy of the
method is independent of the Mach number $\varepsilon$. The scheme is able to
capture shocks and discontinuities in an essentially non-oscillatory fashion
in the compressible regime, meanwhile the scheme is a good incompressible
solver in the limit of vanishing $\varepsilon$. A key feature of the scheme is
the implicit treatment of acoustic waves, while material waves are treated
explicitly by WENO reconstructions of numerical fluxes. In particular, a
suitable local Lax-Friedrich flux with characteristic-wise WENO reconstruction
has been adopted for the explicit convective terms, while component-wise WENO
reconstructions with zero numerical diffusion is used for implicit acoustic
terms. In order to avoid the nonlinearity from the equation of state (EOS), a
semi-implicit treatment similar to the one adopted in [8] is used, leading to
a linearized elliptic equation for the pressure, as described in Section 3.2.
Another essential ingredient of the scheme design is to split the pressure
into a thermodynamic pressure and a hydrodynamic one, using a similar idea
adopted in [14] but with a fixed splitting parameter $\alpha$. The
thermodynamic pressure is used for the characteristic reconstructions, while
the hydrodynamic pressure is obtained by solving an elliptic system. We show
that the resulting scheme is asymptotic preserving (AP) and asymptotic
accurate [28, 29], i.e., it is a consistent and high order discretization of
the compressible Euler equations and, in the limit as $\varepsilon\to 0$, with
$\Delta x$ and $\Delta t$ fixed, if becomes a consistent and high order
discretization of the incompressible Euler equations.

The rest of the paper is organised as follows. We recall the low Mach limit
for the full compressible Euler equations in Section 2. We start Section 3 by
introducing a first order semi-implicit scheme in time, then we describe the
extension to high order time discretization using IMEX methods in Section 3.2,
and we close the section with a description of high order spatial
discretization obtained by characteristic-wise and component-wise WENO
strategies. The asymptotic preserving (AP) and asymptotic accuracy (AA)
properties of the scheme are given in Section 4. Numerical tests are performed
in Section 5, with the conclusion drew in the last section.

##  2 Low Mach limit for the full Euler equations

We consider the compressible Euler equations for an ideal gas in the non-
dimensional form [37, 17, 14]:

$\left\\{\begin{array}[]{ll}\rho_{t}+\nabla\cdot(\rho\mathbf{u})\,=\,0,\\\\[8.53581pt]
(\rho\mathbf{u})_{t}+\nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u})+\frac{1}{{\varepsilon^{2}}}{\nabla
p}\,=\,0,\\\\[8.53581pt]
E_{t}+\nabla\cdot[(E+p)\mathbf{u}]\,=\,0,\end{array}\right.$ (2.1)

with the EOS for a polytropic gas satisfying

$E=\frac{p}{\gamma-1}+\frac{\varepsilon^{2}}{2}\rho|\textbf{u}|^{2},$ (2.2)

where $\gamma>1$ being the ratio of specific heats. The parameter
$\varepsilon$ represents a global Mach number characterizing the non
dimensionalization. System (2.1) is _hyperbolic_ and the eigenvalues along the
direction n are: $\lambda_{1}=\textbf{u}\cdot\textbf{n}-c_{s}/\varepsilon$,
$\lambda_{2}=\textbf{u}\cdot\textbf{n}$,
$\lambda_{3}=\textbf{u}\cdot\textbf{n}-c_{s}/\varepsilon$ with $\displaystyle
c_{s}=\sqrt{\gamma p/\rho}$.

When the reference Mach number is of order one, namely
$\varepsilon=\mathcal{O}(1)$, modern shock capturing methods are able to
compute the formation and evolution of shocks and other complex structures
with high resolutions at a reasonable cost. On the other hand, when the flows
are slow compared to the speed of sound, i.e. $\varepsilon\ll 1$, we are near
the incompressible regime. In such a situation, pressure waves become very
fast compared to material waves. Standard explicit shock-capturing methods
require a CFL time restriction dictated by the sound speed $c_{s}/\varepsilon$
to integrate the system. This leads to the stiffness in time, see e.g. [23,
16, 14], where the time discretization is constrained by a stability condition
given by $\Delta t<\Delta x/\lambda_{\max}=\mathcal{O}(\varepsilon\Delta x),$
here $\Delta t$ is the time step size, $\Delta x$ is the mesh size and
$\lambda_{\max}=\max_{\Omega}(|\textbf{u}|+c_{s}/\varepsilon)$ on the
computational domain $\Omega$. This restriction results in an increasingly
large computational time for low Mach fluid flows. Moreover, excessive
numerical viscosity (scales as $\varepsilon^{-1}$) in standard upwind schemes,
leads to highly inaccurate solutions [45, 46]. Thus, it is of challenge and
great importance to design numerical schemes, not only for shock-capturing,
but also with consideration on stability and consistency in the incompressible
limit, i.e. with _asymptotic preserving_ (AP) property. In fact, in the low
Mach limit, one is not interested in resolving the pressure waves; instead the
fluid pressure serves as a Lagrangian multiplier in preserving the
incompressibility of the velocity field. For the theoretical analysis of
convergence from compressible flow to incompressible equations, such as
$\varepsilon\to 0$, we refer to Klainerman and Majda [30, 31] for a rigorous
study in this low Mach limit.

Here we recall the classical formal derivation of the incompressible Euler
equations from the rescaled compressible Euler equations for an ideal gas
(2.1) with the EOS (2.2). We consider an asymptotic expansion _ansatz_ for the
following two main variables:

$\begin{array}[]{l}p(\textbf{x},t)=p_{0}(\textbf{x},t)+\varepsilon^{2}\,p_{2}(\textbf{x},t)+\cdots,\\\\[8.53581pt]
u(\textbf{x},t)={\bf u}_{0}(\textbf{x},t)+\varepsilon\,{\bf
u}_{1}(\textbf{x},t)+\cdots,\end{array}$ (2.3)

and insert them into the full Euler equations (2.1). First, for the leading
order $\mathcal{O}(\varepsilon^{-2})$, we formally find: $\nabla p_{0}({\bf
x},t)=0$, i.e. pressure $p_{0}$ is constant in space, up to fluctuations of
$\mathcal{O}(\varepsilon^{2})$ and by (2.2), assuming for $E$ an expansion as
$p$ in (2.3), we get $E_{0}=p_{0}/(\gamma-1)$. Then we formally find to
$\mathcal{O}(1)$,

$\left\\{\begin{array}[]{l}\displaystyle\partial_{t}\rho_{0}+\nabla\cdot(\rho_{0}{\bf
u}_{0})\,=\,0,\\\\[8.53581pt] \displaystyle\partial_{t}(\rho_{0}{\bf
u}_{0})+\nabla\cdot\left(\rho_{0}{\bf u}_{0}\otimes{\bf u}_{0}\right)+\nabla
p_{2}\,=\,0,\\\\[8.53581pt]
\displaystyle\partial_{t}E_{0}+\nabla\cdot\left[(E_{0}+p_{0}){\bf
u}_{0}\right]\,=\,0.\end{array}\right.$ (2.4)

From the energy equation in (2.4), with the EOS (2.2) and denoting
$d/dt=\partial/\partial t+{\bf u}_{0}\cdot\nabla$, we have

$\displaystyle\nabla\cdot{\bf u}_{0}=-\frac{1}{p_{0}\gamma}\frac{dp_{0}}{dt}.$
(2.5)

Now integrating the equation (2.5) over a spatial domain $\Omega$ with no-slip
or periodic boundary conditions, we get $\int_{\Omega}\nabla\cdot{\bf
u}_{0}\,d{\bf x}=\int_{\partial\Omega}{\bf u}_{0}\cdot{\bf n}\,ds=0$, where
${\bf n}$ is the unit outward normal vector along $\partial\Omega$, and this
implies $p_{0}$ is constant both in space and time, i.e. $p_{0}=$ Const. Then,
back into (2.5), one finds the divergence constraint $\nabla\cdot{\bf
u}_{0}=0$ in the zero Mach number limit. In summary, we have:

$\left\\{\begin{array}[]{l}p_{0}=\textrm{Const.},\quad\nabla\cdot{\bf{\bf
u}}_{0}=0,\\\\[8.53581pt]
\displaystyle\partial_{t}\rho_{0}+\nabla\cdot(\rho_{0}{\bf
u}_{0})=0,\\\\[8.53581pt] \displaystyle\partial_{t}(\rho_{0}{\bf
u}_{0})+\nabla\cdot\left(\rho_{0}{\bf u}_{0}\otimes{\bf u}_{0}\right)+\nabla
p_{2}=0.\end{array}\right.$ (2.6)

As a direct consequence of $\nabla\cdot{\bf u}_{0}=0$, we obtain from the
mass-continuity equation $d\rho/dt=0$, namely the material derivative of the
density is zero. This means that the density is constant along the particle
trajectories. In particular, if the initial density is constant in space, the
density of the fluid is constant in space and time. Note that

$p_{2}=\lim_{\varepsilon\to 0}\frac{1}{\varepsilon^{2}}(p-p_{0})$ (2.7)

is implicitly defined by the constraint $\nabla\cdot\textbf{u}_{0}=0$, which
satisfies the following elliptic equation:

$\displaystyle-\nabla\cdot\left(\frac{1}{\rho_{0}}\nabla
p_{2}\right)=\nabla\cdot(({\bf u}_{0}\cdot\nabla){\bf u}_{0}).$ (2.8)

Finally, we assume the initial condition is well-prepared [30, 31, 17, 32],
that is, the initial condition for (2.3) is compatible with the equations at
various orders of $\varepsilon$:

$\left\\{\begin{array}[]{l}p(t=0,\textbf{x})=p_{0}+\varepsilon^{2}\,p_{2}(0,\textbf{x})+\cdots\\\\[8.53581pt]
{\bf u}(t=0,\textbf{x})={{\bf
u}}_{0}(\textbf{x})+\mathcal{O}(\varepsilon),\end{array}\right.$ (2.9)

with $p_{0}=\textrm{Const.}$ and $\nabla\cdot{{\bf u}}_{0}=0$ and we impose
$\rho(0,{\bf x})=\rho_{0}({\bf x})$, with $\rho_{0}({\bf x})$ being a strictly
positive function. Note that well-prepared initial conditions are required if
we want that the solution to the $\varepsilon$-dependent problem smoothly
converges to the solution of the limit incompressible problem. Furthermore,
well-prepared initial condition is an important requirement to design AP
schemes. It is crucial to preserve the constant state for leading order terms
of $p$ and $E$, as well as that the divergence free constraint on the leading
order term of ${\bf u}$. For an arbitrary initial condition, an initial layer
will appear, which require a numerical resolution at the $\varepsilon$-scale.

##  3 Numerical scheme

In this section, we aim to construct and analyze a class of high order finite
difference schemes with the AP property for unsteady compressible flows, when
the Mach number $\varepsilon$ spans several orders of magnitude. The features
of our scheme are the following: we design a semi-implicit IMEX (SI-IMEX) time
discretization strategy, so that the scheme is stable with a time stepping
constraint independent of the Mach number $\varepsilon$, the AP property is
preserved in the zero Mach number limit and the scheme can be implemented in a
semi-implicit manner [3, 2, 1, 4] to enable effective and efficient numerical
implementations. Our scheme preserves the incompressible velocity field in the
zero Mach number limit by involving an elliptic solver for the hydrostatic
pressure. In this section, the strategy of numerical discretizations is
different from the traditional method-of-lines approach, since we first
perform time discretization to ensure AP property, after which we apply a
suitable space discretization. In particular we adopt high order WENO
strategies with characteristic reconstructions tailored to IMEX-type methods
in time. The final scheme can successfully capture shocks in the compressible
regime, and efficiently solves the equations in the low Mach regime, with CFL
condition depending only on fluid velocity.

###  3.1 Semi-implicit treatment

We introduce our proposed strategy of implicit and explicit time
discretizations, which is similar in spirit to the first order scheme in [14]
with a slight modification. We emphasize the special treatment to avoid
solving nonlinear equations required by a fully implicit scheme. We rewrite
(2.1) as

$\frac{dU}{dt}=-\nabla\cdot\mathcal{F}_{E}-\nabla\cdot\mathcal{F}_{SI}$ (3.10)

where $U=(\rho,\rho{\bf u},E)^{T}$ and

$\mathcal{F}_{E}\doteq\left(\begin{array}[]{c}\displaystyle{\bf
q}_{E}\\\\[8.53581pt] \displaystyle\left(\frac{{\bf q}_{E}\otimes{\bf
q}_{E}}{\rho_{E}}\right)+{\alpha}\,p_{E}\mathbb{I}\\\\[8.53581pt]
\displaystyle
0\end{array}\right),\qquad\mathcal{F}_{SI}\doteq\left(\begin{array}[]{c}\displaystyle
0\\\\[8.53581pt]
\displaystyle\frac{1-{\alpha}\varepsilon^{2}}{\varepsilon^{2}}p_{I}\mathbb{I}\\\\[8.53581pt]
\displaystyle\frac{E_{E}+p_{E}}{\rho_{I}}{\bf q}_{I}\end{array}\right).$
(3.11)

Subscripts $E$ and $SI$ of $\mathcal{F}$ indicate the explicit and semi-
implicit treatment of the first and the second term respectively. Several
remarks are in order:

  1. 1.

The parameter $\alpha$ determines the splitting between the explicit and
implicit contribution of the pressure: the former ensures that the explicit
part is still hyperbolic, with a much smaller sound speed than the physical
one, when $\varepsilon\ll 1$, while the latter will ensure much milder
stability restrictions. As $\varepsilon$ increases, the explicit contribution
becomes more and more relevant. This form of splitting is similar to the one
in [14], but differs in two main aspects: one is the implicit treatment of
$\rho$ in the energy equation in (3.11), which gives a better asymptotic
preserving property as will be elaborated in Section 4; the other is the
choice of the parameter $\alpha$ in splitting the pressure. In [14], $\alpha$
is chosen depending on the Mach number $\varepsilon$. In our proposed scheme,
however, $\alpha$ is chosen to be equal to 1 for all $\varepsilon<1$, and
$\alpha=1/\varepsilon^{2}$ for $\varepsilon\geq 1$.

  2. 2.

We use the following EOS for (2.2) to avoid the nonlinearity in the semi-
implicit scheme:

$\displaystyle E_{E}=\frac{1}{\gamma-1}p_{E}+\varepsilon^{2}\frac{|{\bf
q}_{E}|^{2}}{2\rho_{E}},\quad\displaystyle
E_{I}=\frac{1}{\gamma-1}p_{I}+\varepsilon^{2}\frac{|{\bf
q}_{E}|^{2}}{2\rho_{E}}.$ (3.12)

The subscripts $E$ and $I$ of ${\bf q}$, $\rho$, $E$ and $p$ indicate the
explicit and implicit treatments of the corresponding variables, respectively.

  3. 3.

For the implicit term $p_{I}$ in (3.11), it is convenient to introduce a
pressure perturbation $p_{I,2}$, corresponding to the hydrodynamic pressure in
the incompressible limit, defined as

$p_{I,2}\doteq\frac{p_{I}-\bar{p}_{E}}{{\varepsilon^{2}}},$ (3.13)

where $\bar{p}_{E}$ denotes the spatial average of the pressure $p_{E}$
computed from the EOS (3.12). In this way the term $p_{I,2}$ will remain
finite even as $\varepsilon\to 0$. Then we have

$\frac{1}{\varepsilon^{2}}\nabla p_{I}=\nabla p_{I,2}.$ (3.14)

As an example, we present the scheme as well as the flow chart to update the
numerical solution $U^{n+1}=(\rho^{n+1},{\bf q}^{n+1},E^{n+1})^{T}$ for the
first order semi-implicit scheme solving system (2.1). We focus on the time
discretization while keeping the space continuous, whose discretizations will
be discussed in detail in Section 3.3:

$\displaystyle\frac{\rho^{n+1}-\rho^{n}}{{\Delta t}}+\nabla\cdot{\bf
q}^{n}=0,$ (3.15a) $\displaystyle\frac{{\bf q}^{n+1}-{\bf q}^{n}}{{\Delta
t}}+\nabla\cdot\left(\frac{{\bf q}^{n}\otimes{\bf
q}^{n}}{\rho^{n}}+p^{n}\,{\mathbb{I}}\right)+\frac{1-{\varepsilon^{2}}}{\varepsilon^{2}}\,\nabla
p^{n+1}=0,$ (3.15b) $\displaystyle\frac{E^{n+1}-E^{n}}{{\Delta
t}}+\nabla\cdot\left(\frac{E^{n}+p^{n}}{\rho^{n+1}}{\bf q}^{n+1}\right)=0,$
(3.15c)

with

$\mathcal{F}_{E}(U^{n})\doteq\left(\begin{array}[]{c}\displaystyle{\bf
q}^{n}\\\\[8.53581pt] \displaystyle\left(\frac{{\bf q}^{n}\otimes{\bf
q}^{n}}{\rho^{n}}\right)+p^{n}\,\mathbb{I}\\\\[8.53581pt] \displaystyle
0\end{array}\right),\quad\mathcal{F}_{SI}(U^{n},U^{n+1})\doteq\left(\begin{array}[]{c}\displaystyle
0\\\\[8.53581pt]
\displaystyle\frac{1-\varepsilon^{2}}{\varepsilon^{2}}\,p^{n+1}\,\mathbb{I}\\\\[8.53581pt]
\displaystyle\frac{E^{n}+p^{n}}{\rho^{n+1}}{\bf q}^{n+1}\end{array}\right).$
(3.16)

The flow chart based on the semi-implicit scheme is the following:

  1. 1.

update $\rho^{n+1}$ from (3.15a).

  2. 2.

we rewrite (3.15b) as

${\bf q}^{n+1}={{\bf q}^{*}}-{{\Delta
t}}\frac{1-{\varepsilon^{2}}}{\varepsilon^{2}}\nabla p^{n+1},$ (3.17)

with ${{\bf q}^{*}}={{\bf q}^{n}}-{{\Delta t}}\nabla\cdot\left(\frac{{\bf
q}^{n}\otimes{\bf q}^{n}}{\rho^{n}}+p^{n}{\mathbb{I}}\right)$. We substitute
${\bf q}^{n+1}$ into (3.15c) to get

$E^{n+1}=E^{*}+{{\Delta
t}}^{2}\frac{1-{\varepsilon^{2}}}{\varepsilon^{2}}\nabla\cdot\left(H^{n}\nabla
p^{n+1}\right),$ (3.18)

where $H^{n}=(E^{n}+p^{n})/\rho^{n+1}$ and $E^{*}=E^{n}-{{\Delta
t}}\nabla\cdot(H^{n}{{\bf q}^{*}})$.

  3. 3.

Now we replace $E^{n+1}$ by $p^{n+1}/(\gamma-1)+\varepsilon^{2}{({\bf
q}^{n})^{2}}/{(2\rho^{n})}$ in (3.18) using (3.12), and with the introduction
of (3.13) i.e. $p^{n+1}_{2}\doteq({p^{n+1}-\bar{p}^{n}})/{{\varepsilon^{2}}}$
where $p_{I}=p^{n+1}$ and $\bar{p}^{n}=\bar{p}_{E}$, we rewrite (3.18) as

$\frac{\varepsilon^{2}}{\gamma-1}p^{n+1}_{2}-{{\Delta
t}}^{2}(1-\varepsilon^{2})\nabla\cdot\left(H^{n}\nabla
p_{2}^{n+1}\right)=E^{**},$ (3.19)

where $E^{**}=E^{*}-\bar{p}^{n}/({\gamma-1})-\varepsilon^{2}{({\bf
q}^{n})^{2}}/({2\rho^{n}})$ is explicitly computed. We obtain an elliptic
equation (3.19) for $p^{n+1}_{2}$.

  4. 4.

Finally, we update ${\bf q}^{n+1}$ from (3.17), and then $E^{n+1}$ from
(3.15c).

Note that if $\varepsilon\geq 1$ the implicit pressure contribution in (3.15b)
vanishes, so the momentum ${\bf q}^{n+1}$ is evaluated explicitly. With
updated $\rho^{n+1}$ and ${\bf q}^{n+1}$, $E^{n+1}$ can also be updated in an
explicit way from (3.15c).

###  3.2 High order semi-implicit (SI) temporal discretization using IMEX

####  3.2.1 High order SI-IMEX R-K scheme for all Mach number full Euler
system

We generalize the first order scheme (3.15) to high order in the framework of
IMEX R-K methods. In order to do so, we follow the idea introduced in [3]. We
consider an autonomous system of the form

$U^{\prime}={\mathcal{H}}(U,U),\quad U(t_{0})=U_{0},$

where $\mathcal{H}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$
is a sufficiently regular mapping. We assume an explicit treatment of the
first argument of ${\mathcal{H}}$ (using subscript ”E”), and an implicit
treatment to the second argument (using subscript ”I” ), i.e.

$\left\\{\begin{array}[]{l}U^{\prime}_{E}=\mathcal{H}(U_{E},U_{I}),\\\\[8.53581pt]
U^{\prime}_{I}=\mathcal{H}(U_{E},U_{I}),\end{array}\right.$ (3.20)

with initial conditions

$U_{E}(t_{0})=U_{0},\quad U_{I}(t_{0})=U_{0}.$ (3.21)

Then system (3.10) with (3.11) can be rewritten in the form (3.20) where
$U\doteq(U_{E},U_{I})^{T}$ with $U_{E}=(\rho_{E},{\bf q}_{E},E_{E})$, and
$U_{I}=(\rho_{I},{\bf q}_{I},E_{I})$, and
$\mathcal{H}(U_{E},U_{I})=-\nabla\cdot\mathcal{F}_{E}-\nabla\cdot\mathcal{F}_{SI}$.

System (3.20) is a particular case of the _partitioned system_ , [25]. One can
apply an IMEX R-K scheme to (3.20), using the corresponding pair of Butcher
$tableau$ [9],

$\begin{array}[]{c|c}\tilde{c}&\tilde{A}\\\ \hline\cr\vspace{-0.25cm}\hfil\\\
&\tilde{b^{T}}\end{array},\ \ \ \ \ \begin{array}[]{c|c}{c}&{A}\\\
\hline\cr\vspace{-0.25cm}\hfil\\\ &{b^{T}}\end{array},$ (3.22)

where $\tilde{A}=(\tilde{a}_{ij})$ is an $s\times s$ matrix for an explicit
scheme, with $\tilde{a}_{ij}=0$ for $j\geq i$ and $A=({a}_{ij})$ is an
$s\times s$ matrix for an implicit scheme. For the implicit part of the
methods, we use a diagonally implicit scheme, i.e. $a_{ij}=0$, for $j>i$, in
order to guarantee simplicity and efficiency in solving the algebraic
equations corresponding to the implicit part of the discretization. The
vectors $\tilde{c}=(\tilde{c}_{1},...,\tilde{c}_{s})^{T}$,
$\tilde{b}=(\tilde{b}_{1},...,\tilde{b}_{s})^{T}$, and
$c=(c_{1},...,c_{s})^{T}$, $b=(b_{1},...,b_{s})^{T}$ complete the
characterization of the scheme. The coefficients $\tilde{c}$ and $c$ are given
by the usual relation

$\displaystyle\tilde{c}_{i}=\sum_{j=1}^{i-1}\tilde{a}_{ij},\ \ \
c_{i}=\sum_{j=1}^{i}a_{ij}.$ (3.23)

Later, for the consideration of AP property, we consider stiffly accurate (SA)
implicit schemes, i.e., the implicit part of the Butcher table satisfies the
condition $b^{T}=e^{T}_{s}A$, with $e_{s}=(0,...,0,1)$ and $c_{s}=1$. We will
see that SA guarantees that the numerical solution is identical to the last
internal stage value of the scheme.

Now an SI-IMEX R-K scheme applied to (3.20) reads

$U^{(i)}_{E}=U^{n}_{E}+\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\mathcal{H}(U^{(j)}_{E},U^{(j)}_{I}),\quad
U^{(i)}_{I}=U^{n}_{I}+\Delta
t\sum_{j=1}^{i}{a}_{ij}\mathcal{H}(U^{(j)}_{E},U^{(j)}_{I}),$ (3.24a)
$U^{n+1}_{E}=U^{n}_{E}+\Delta
t\sum_{i=1}^{s}\tilde{b}_{i}\mathcal{H}(U^{(i)}_{E},U^{(i)}_{I}),\quad
U^{n+1}_{I}=U^{n}_{I}+\Delta
t\sum_{i=1}^{s}{b}_{i}\mathcal{H}(U^{(i)}_{E},U^{(i)}_{I}).$ (3.24b)

Apparently, by writing in the form (3.20), we increased the computational cost
since we double the number of variables. However, the cost is mainly due to
the number of function evaluations, which essentially depends on the semi-
implicit scheme, because in (3.24) the identical term
$\mathcal{H}(U^{j}_{E},U^{j}_{I})$ appears in both explicit and implicit part.
Furthermore, if $\tilde{b}_{i}=b_{i}$ for all $i$, then $U^{n}_{E}=U^{n}_{I}$
for all $n>0$ (provided $U^{0}_{E}=U^{0}_{I}$) and therefore the duplication
of variables is not necessary. Such a property is of particular relevance for
designing time discretization schemes which are asymptotic preserving.

We rewrite the scheme (3.24) in the following new form for the convenience of
further discussion

$U^{(i)}_{E}=U^{n}+\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\mathcal{H}(U^{(j)}_{E},U^{(j)}_{I}),\quad{U}_{I}^{(i)}=\tilde{U}_{I}^{(i)}+\Delta
t\,{a}_{ii}\mathcal{H}(U^{(i)}_{E},U^{(i)}_{I}),$ (3.25)

where $\tilde{U}_{I}^{(i)}=U^{n}+\Delta
t\sum_{j=1}^{i-1}{a}_{ij}\mathcal{H}(U^{(j)}_{E},U^{(j)}_{I})$. Denoting
$U^{n}=U^{n}_{E}=U^{n}_{I}$ for all $n>0$ (due to $\tilde{b}_{i}=b_{i}$ for
all $i$), finally we have

$U^{n+1}=U^{n}+\Delta
t\sum_{i=1}^{s}{b}_{i}\mathcal{H}(U^{(i)}_{E},U^{(i)}_{I}).$ (3.26)

Similar to the first order scheme, we split the pressure into explicit
($p_{E}$) and implicit ($p_{I}$) terms. In order to avoid the nonlinearity in
the semi-implicit step we use the following equations of state (3.12), i.e.,
for the generic stages $i=1,\ldots,s$

$\displaystyle
E^{(i)}_{E}=\frac{1}{\gamma-1}p^{(i)}_{E}+\varepsilon^{2}\frac{|{\bf
q}^{(i)}_{E}|^{2}}{2\rho^{(i)}_{E}},$ (3.27a) $\displaystyle
E^{(i)}_{I}=\frac{1}{\gamma-1}p^{(i)}_{I}+\varepsilon^{2}\frac{|{\bf
q}^{(i)}_{E}|^{2}}{2\rho^{(i)}_{E}}.$ (3.27b)

######  Remark 3.1.

The first order scheme (3.15) is the same as applying the following Butcher
table to (3.20)

$\mbox{\bf Type A I:}\quad\begin{array}[]{c|c}0&0\\\
\hline\cr&1\end{array}\qquad\begin{array}[]{c|c}1&1\\\
\hline\cr&1\end{array}.$ (3.28)

Formally applying the above tableau (3.28) to the partitioned system (3.20),
it reads

$\begin{array}[]{ll}U^{(1)}_{E}&=U^{n},\\\\[5.69054pt]
U^{(1)}_{I}&=U^{n}+\Delta
t\mathcal{H}(U^{(1)}_{E},U^{(1)}_{I}),\\\\[5.69054pt] U^{n+1}&=U^{n}+\Delta
t\mathcal{H}(U^{(1)}_{E},U^{(1)}_{I}),\end{array}$

which is the first order scheme (3.15), considering $U^{(1)}_{E}=U^{n}$ and
$U^{(1)}_{I}=U^{n+1}$.

####  3.2.2 Construction of an IMEX R-K scheme

Next, we construct an IMEX R-K Butcher tableau for semi-implicit
discretization of all-Mach full Euler system. The construction is based on the
following considerations for accuracy and for handling non-well prepared
initial data.

  1. 1.

We require that the implicit part of the IMEX R-K scheme is SA. With such an
assumption, one can derive the AP and AA properties of the scheme as discussed
in Section 4. Note that if the implicit part of the scheme is $A$-stable, SA
is a sufficient condition to make it $L$-stable, see [25].

  2. 2.

The invertibility of the implicit matrix $A$ of the SI-IMEX R-K scheme is
important to handle non-well prepared initial conditions and for proper
initialization of the hydrodynamic pressure. For further related discussions
and analysis, we refer the reader to the following papers [38, 6, 5].

  3. 3.

To synchronize $U_{E}$ and $U_{I}$, the weights for the final stage of the
double tableaus should be the same, i.e. $\tilde{b}_{i}=b_{i}$, $i=1\cdots s$.
That is, we keep only one set of numerical solution in the process of updating
[3]. Alternatively, we can select a different vector of weights for the
$U_{E}(t)$, say $\tilde{b}_{i}\neq b_{i}$, which will provide a lower/higher
order approximation of the solution for $U_{E}(t)$; this can be used to
implement a procedure of automatic time step control [24].

  4. 4.

In order to simplify order conditions, given the equations of state used in
the semi-implicit step (3.27), we impose that $c_{i}=\tilde{c}_{i},\,i=2\cdots
s$. Note that if matrix $A$ is invertible we have $c_{1}\neq\tilde{c}_{1}=0$
[38].

Based on the above considerations, we design a high-order IMEX scheme with
matrix $A$ invertible, which is SA and satisfies $c_{i}=\tilde{c}_{i}$ for
$i=2,\cdots,s$. We impose the conditions required for a third order scheme
[38] and obtain the following double Butcher tableau

Explicit : $\displaystyle\begin{array}[]{c|cccc}0&0&0&0&0\\\
\gamma&\gamma&0&0&0\\\ 0.717933260754&0.435866521508&0.282066739245&0&0\\\
1&-0.733534082748750&2.150527381100&-0.416993298352&0\\\ \hline\cr
0&0&1.208496649176&-0.644363170684&\gamma\end{array},$ (3.34) Implicit :
(3.35) $\displaystyle\begin{array}[]{c|cccc}\gamma&\gamma&0&0&0\\\
\gamma&0&\gamma&0&0\\\ 0.717933260754&0&0.282066739245&\gamma&0\\\
1&0&1.208496649176&-0.644363170684&\gamma\\\
\hline\cr&0&1.208496649176&-0.644363170684&\gamma\end{array},$ (3.41)

with $\gamma=0.435866521508$. We call it SI-IMEX(4,4,3), the triplet
$(s,\sigma,p)$ characterizing the number of stages of the implicit scheme
($s=4$), the number of stages of the explicit scheme ($\sigma=4$ ) and the
order of the scheme ($p=3$).

###  3.3 High order spatial discretization for all-Mach fluid flows

Below we describe our spatial discretization strategies that incorporate WENO
mechanism to capture shocks in the compressible regime and produce a high
order incompressible solver for the flow in the zero Mach limit. One major
difficulty, hence the new ingredient in the scheme design, in extending the
high order AP scheme for the isentropic Euler system [7] to the full Euler
system is about the pressure. In the isentropic case, the pressure is an
explicit function of $\rho$, whereas in the full Euler case the pressure comes
from the EOS involving all conserved variables. As such, we split the pressure
into the part involving $p_{E}$ obtained from the EOS (3.27a), and the part
$p_{I,2}$ obtained by formulating an elliptic equation from a semi-implicit
solver of the system. In this process, the spatial discretization in the
scheme formulation becomes critical for its robustness. There are several
ingredients in our spatial discretization that differentiate our approach from
the existing ones in the literature. Firstly, we carefully apply the WENO
procedure for spatial reconstructions of fluxes for the full Euler system
which is a nontrivial generalization from the isentropic one [7]. In
particular, we propose to apply the fifth order characteristic-wise WENO
procedure to $\mathcal{F}_{E}$ (in (3.11)) of the full Euler system so that in
the compressible regime oscillations could be best controlled; and to apply a
component-wise WENO procedure for the flux functions of $\mathcal{F}_{SI}$ (in
(3.11)). We found that such a choice is optimal in compromising the need from
both compressible and incompressible fluid solvers. The characteristic-wise
WENO turns out to be important for compressible Euler system, as shown in the
numerical section. Secondly, following the strategy in our previous work [7],
we apply a compact high order spatial discretization to second order
derivatives terms in (3.47). This is discussed in Section 3.3.2.

####  3.3.1 WENO spatial discretization

Without loss of generality, we describe our algorithms in a 2D setting. Let

$U=(\rho,\rho u,\rho v,E)^{T},$
$\nabla\cdot\mathcal{F}_{E}=\partial_{x}\mathcal{F}^{x}_{E}+\partial_{y}\mathcal{F}^{y}_{E},\quad\nabla\cdot\mathcal{F}_{SI}=\partial_{x}\mathcal{F}^{x}_{SI}+\partial_{y}\mathcal{F}^{y}_{SI}$

with flux functions in $x$ and $y$ directions

$\mathcal{F}^{x}_{E}=(\rho u,\rho u^{2}+p,\rho
uv,0)^{T},\quad\mathcal{F}^{y}_{E}=(\rho v,\rho uv,\rho v^{2}+p,0)^{T},$
$\mathcal{F}^{x}_{SI}=\left(0,\frac{1-\varepsilon^{2}}{\varepsilon^{2}}p,0,\frac{E+p}{\rho}(\rho
u)\right)^{T},\quad\mathcal{F}^{y}_{SI}=\left(0,0,\frac{1-\varepsilon^{2}}{\varepsilon^{2}}p,\frac{E+p}{\rho}(\rho
v)\right)^{T}.$

We consider a rectangular domain discretized by a uniform Cartesian grid, with
meshsize $\Delta x=\Delta y$, and grid points $(x_{i},y_{j})$, $i=1,\cdots
N_{x}$, $j=1\cdots N_{y}$, located at cell centers. We use the subscript
$\cdot_{i,j}$ to denote the solution point values at $(x_{i},y_{j})$, and
$\cdot_{i\pm 1/2,j}$ and $\cdot_{i,j\pm 1/2}$ for the reconstructed numerical
fluxes in approximating $x$\- an $y$\- derivatives respectively.

Below we first present two types of WENO spatial discretizations, i.e, a
characteristic-wise WENO and a component-wise WENO, for explicit and semi-
implicit parts as specified in (3.11), respectively.

  1. 1.

Characteristic-wise WENO $\nabla_{CW}$ with Lax-Friedrichs splitting for
discretizing $\nabla\cdot\mathcal{F}_{E}$. For compressible hyperbolic systems
with shocks, WENO reconstructions in the component-wise fashion may still lead
to oscillations and a local characteristic decomposition will be needed [39].
We let $\mathcal{F}^{x}_{E,i,j}$ be the flux function at each grid point
$(x_{i},y_{j})$. To approximate the
$\partial_{x}\mathcal{F}^{x}_{E}|_{x_{i},y_{j}}$, we fix a $j$-index and
perform a global Lax-Friedrichs splitting for the flux terms, i.e.

$\mathcal{F}^{x}_{E,i,j}=\mathcal{F}^{x,+}_{E,i,j}+\mathcal{F}^{x,-}_{E,i,j},\;\forall
i$

with

$\mathcal{F}^{x,\pm}_{E,i,j}=\frac{1}{2}(\mathcal{F}^{x}_{E,i,j}\pm\Lambda\,U_{E,i,j}),\quad\Lambda=\max_{u,v}\\{|u|+|v|+\min\\{1,\frac{1}{\varepsilon}\\}\,c_{s}\\}$
(3.42)

with $c_{s}=\sqrt{\gamma p/\rho}$ being the sound speed multipled by
$\varepsilon$, and the $\max$ is taken over appropriate range of $u$, $v$.
Notice that such a choice of $\Lambda$ is the same as in our previous work
[7]. We then project $\\{\mathcal{F}^{x,\pm}_{E,i,j}\\}_{i=1}^{N_{x}}$ to
local characteristics directions for the full Euler system in the compressible
regime (corresponding to $\varepsilon=1$), perform WENO reconstruction of
fluxes there and then project the reconstructed fluxes back, to obtain the
flux terms as

$\hat{\mathcal{F}}^{x}_{E,i+\frac{1}{2},j}=\hat{\mathcal{F}}^{x,+}_{E,i+\frac{1}{2},j}+\hat{\mathcal{F}}^{x,-}_{E,i+\frac{1}{2},j}.$

Then

$\partial_{x}\mathcal{F}^{x}_{E}|_{x_{i},y_{j}}\approx\frac{1}{\Delta
x}(\hat{\mathcal{F}}^{x}_{E,i+\frac{1}{2},j}-\hat{\mathcal{F}}^{x}_{E,i-\frac{1}{2},j}).$

Similar procedure could be performed to approximate spatial derivatives in the
$y$-direction. We introduce the notation of $\nabla_{CW}$ for the
characteristic-wise WENO spatial discretization of $\mathcal{F}_{E}$.

  2. 2.

A component-wise WENO $\nabla_{W}$ with Lax-Friedrichs flux splitting for
discretizing $\nabla\cdot\mathcal{F}_{SI}$. We apply a component-wise WENO
reconstruction of fluxes for discretizing $\nabla\cdot\mathcal{F}_{SI}$
denoted as $\nabla_{W}\cdot\mathcal{F}_{SI}$. In particular, a Lax-Friedrichs
flux splitting, followed by a component-wise WENO reconstruction is performed
to obtain the numerical fluxes.

####  3.3.2 Flowchart and compact discretization of second order derivatives

With the introduction of characteristic-wise and component-wise WENO
$\nabla_{CW}$ and $\nabla_{W}$, we summarize the flow chart of the SI-IMEX R-K
time discretization (3.25)-(3.26), coupled with WENO spatial discretization,
for solving all-Mach full Euler system. In this process, a compact high order
discretization for second order spatial derivatives in (3.47) is applied.

  1. 1.

Start from $U^{n}$ at time $t^{n}$, we first compute $U_{E}^{(i)}$ from (3.25)

$\displaystyle\rho^{(i)}_{E}=\rho^{n}-\Delta
t\,\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla_{CW}\cdot{\bf q}_{E}^{(j)},\\\ $
(3.43a) $\displaystyle{{\bf q}}^{(i)}_{E}={\bf q}^{n}-{\Delta
t}\,\sum_{j=1}^{i-1}\tilde{a}_{ij}\left(\nabla_{CW}\cdot\left(\frac{{\bf
q}^{(j)}_{E}\otimes{\bf
q}^{(j)}_{E}}{\rho^{(j)}_{E}}+p^{(j)}_{E}\,{\mathbb{I}}\right)+(1-\varepsilon^{2})\frac{\nabla_{W}p_{I}^{(j)}}{\varepsilon^{2}}\right),\\\
$ (3.43b) $\displaystyle{E}^{(i)}_{E}=E^{n}-\Delta
t\,\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla_{W}\cdot\left(\bar{H}^{(j)}{{\bf
q}}^{(j)}_{I}\right).$ (3.43c)

  2. 2.

We compute $\tilde{U}^{(i)}$ in (3.25).

$\displaystyle\tilde{\rho}^{(i)}=\rho^{n}-\Delta
t\,\sum_{j=1}^{i-1}{a}_{ij}\nabla_{CW}\cdot{\bf q}_{E}^{(j)},$ (3.44a)
$\displaystyle\tilde{{\bf q}}^{(i)}={\bf q}^{n}-{\Delta
t}\,\sum_{j=1}^{i-1}{a}_{ij}\left(\nabla_{CW}\cdot\left(\frac{{\bf
q}^{(j)}_{E}\otimes{\bf
q}^{(j)}_{E}}{\rho^{(j)}_{E}}+p^{(j)}_{E}\,{\mathbb{I}}\right)+(1-\varepsilon^{2})\frac{\nabla_{W}p_{I}^{(j)}}{\varepsilon^{2}}\right),$
(3.44b) $\displaystyle\tilde{E}^{(i)}=E^{n}-\Delta
t\,\sum_{j=1}^{i-1}a_{ij}\nabla_{W}\cdot\left(\bar{H}^{(j)}{{\bf
q}}^{(j)}_{I}\right),$ (3.44c)

with $\displaystyle\bar{H}^{(j)}=(E^{(j)}_{E}+p^{(j)}_{E})/\rho^{(j)}_{I}$.

  3. 3.

Solve $U_{I}^{(i)}$ from (3.25).

    1. (a)

In components, $U_{I}^{(i)}$ satisfies

$\displaystyle{\rho}^{(i)}_{I}=\displaystyle\tilde{\rho}^{(i)}-\Delta
t\,a_{ii}\nabla_{CW}\cdot{\bf q}^{(i)}_{E},$ (3.45a) $\displaystyle{{\bf
q}}^{(i)}_{I}=\displaystyle\tilde{\tilde{{\bf q}}}^{(i)}-\Delta
t\,a_{ii}\left(\frac{1-\varepsilon^{2}}{\varepsilon^{2}}\nabla
p_{I}^{(i)}\right),$ (3.45b) $\displaystyle E_{I}^{(i)}=\tilde{E}^{(i)}-\Delta
ta_{ii}\nabla\cdot(\bar{H}_{i}{\bf q}_{I}^{(i)}),$ (3.45c)

where

$\tilde{\tilde{{\bf q}}}^{(i)}=\displaystyle\tilde{{\bf q}}^{(i)}-\Delta
t\,a_{ii}\left(\nabla_{CW}\cdot\left(\frac{{\bf q}^{(i)}_{E}\otimes{\bf
q}^{(i)}_{E}}{\rho^{(i)}_{E}}+p^{(i)}_{E}\,{\mathbb{I}}\right)\right).$ (3.46)

    2. (b)

To solve the system (3.45), we substitute ${{\bf q}}^{(i)}_{I}$ of (3.45b)
into (3.45c) and obtain

$\displaystyle E_{I}^{(i)}=E^{**}_{I}+(1-\varepsilon^{2})\Delta
t^{2}\,a_{ii}^{2}\nabla\left(\bar{H}^{(i)}\left(\frac{\nabla
p_{I}^{(i)}}{\varepsilon^{2}}\right)\right),$ (3.47)

with $E_{I}^{(i)}$ following the EOS (3.27b) and
$E_{I}^{**}=\tilde{E}^{(i)}-\Delta
t\,a_{ii}\nabla_{W}\cdot\left(\bar{H}^{(i)}\tilde{\tilde{{\bf
q}}}^{(i)}\right)$. (3.47) is an implicit equation about $p^{(i)}_{I}$. By the
introduction of pressure perturbation $p^{(i)}_{I,2}$ in (3.13), we solve
$p^{(i)}_{I,2}$ from

$\displaystyle\frac{\varepsilon^{2}}{\gamma-1}\,p_{I,2}^{(i)}=E^{***}_{I}+(1-\varepsilon^{2})\Delta
t^{2}\,a_{ii}^{2}\nabla\left(\bar{H}^{(i)}\left(\nabla
p_{I,2}^{(i)}\right)\right),$ (3.48)

where $E^{***}_{I}=E^{**}_{I}-\bar{p}_{E}/(\gamma-1)-\varepsilon^{2}|{\bf
q}^{(i)}_{E}|^{2}/(2\rho^{(i)}_{E})$. Notice that in the process of
substitution to obtain (3.47) or (3.48), the gradient and the divergence
operators are kept to be continuous, obtaining the second order operator
$\nabla\left(\bar{H}^{(i)}\nabla p_{I,2}^{(i)}\right)$. The second order
spatial derivative is then discretized by a compact discretization as proposed
in [7].

    3. (c)

With $p^{(i)}_{I,2}$ solved from (3.48), we update ${\bf q}^{(i)}_{I}$ from

$\displaystyle{{\bf q}}^{(i)}_{I}=\displaystyle\tilde{\tilde{{\bf
q}}}^{(i)}-\Delta
t\,a_{ii}\left(1-\varepsilon^{2}\right)\nabla_{W}p_{I,2}^{(i)},$ (3.49)

and successively update $E^{(i)}_{I}$ from

$\displaystyle E_{I}^{(i)}=\tilde{E}^{(i)}-\Delta
ta_{ii}\nabla_{W}\cdot(\bar{H}_{i}{\bf q}_{I}^{(i)}).$ (3.50)

  4. 4.

Finally, update the numerical solution $U^{n+1}=U^{(s)}_{I}$ with the
assumption on the SA property of IMEX R-K schemes.

This completes the description of the high order SI-IMEX R-K time
discretization to the all-Mach full Euler equations.

##  4  Asymptotic preserving (AP) and Asymptotically Accurate (AA) property

###  4.1 AP property

In this section we prove the AP property of scheme (3.15). In particular, we
prove that its limiting scheme is consistent with the continuous limit model
(2.6) at $\varepsilon=0$. We focus on the AP analysis on time discretizations,
while keeping the space continuous. We assume that the data at time $t^{n}$
are well-prepared in the sense of (2.3), i.e. $p^{n}({\bf x}):=p(t^{n},{\bf
x})$ and ${\bf u}^{n}({\bf x}):={\bf u}(t^{n},{\bf x})$ admit the
decomposition:

$p^{n}({\bf x})=p^{n}_{0}+\varepsilon^{2}\,p_{2}^{n}({\bf x}),\quad{\bf
u}^{n}({\bf x})={\bf u}_{0}^{n}({\bf x})+\mathcal{O}(\varepsilon),$ (4.51)

where $p^{n}_{0}=(\gamma-1)E^{n}_{0}$ is a constant and $\nabla\cdot{\bf
u}^{n}_{0}({\bf x})=0$.

Then we consider an expansion in powers of $\varepsilon$ of the form (2.3) and
we plugin it into the semi-discrete scheme (3.15). Equating to zero the
$\mathcal{O}(\varepsilon^{-2})$ term we have

$\nabla p_{0}^{n+1}=0.$ (4.52)

Equating to zero the $\mathcal{O}(\varepsilon^{0})$ terms, we have

$\displaystyle\frac{\rho_{0}^{n+1}-\rho_{0}^{n}}{\Delta
t}+\nabla\cdot(\rho^{n}_{0}{\bf u}_{0}^{n})=0,$ (4.53a)
$\displaystyle\frac{\rho^{n+1}_{0}{\bf u}_{0}^{n+1}-\rho_{0}^{n}{\bf
u}_{0}^{n}}{\Delta t}+\nabla\cdot(\rho^{n}_{0}{\bf u}_{0}^{n}\otimes{\bf
u}_{0}^{n})+\nabla p_{2}^{n+1}=0,$ (4.53b)
$E_{0}^{n+1}=p_{0}^{n+1}/(\gamma-1),$ (4.53c)
$\displaystyle\frac{E_{0}^{n+1}-E_{0}^{n}}{\Delta
t}+\nabla\cdot\left(\bar{H}^{n}_{0}(\rho_{0}^{n+1}{\bf
u}_{0}^{n+1})\right)=0,$ (4.53d)

with

$\bar{H}^{n}_{0}=\frac{p_{0}^{n}+E_{0}^{n}}{\rho_{0}^{n+1}}=\frac{\gamma}{\gamma-1}\frac{p^{n}_{0}}{\rho^{n+1}_{0}}.$
(4.54)

By (4.53c) and (4.54), (4.53d) is equivalent to

$\displaystyle\frac{p_{0}^{n+1}-p_{0}^{n}}{\Delta
t}+\gamma\,p^{n}_{0}\,\nabla\cdot{\bf u}_{0}^{n+1}=0.$ (4.55)

Note that (4.55) is the time discretization of the limiting equation for the
pressure (2.5).

From (4.52), it follows that $p^{n+1}_{0}$ is independent of space. Now, as in
the continuous case, integrating (4.55) over a spatial domain $\Omega$ and
assuming some boundary conditions (for example, no-slip or periodic), we get
$p^{n+1}_{0}=p^{n}_{0}$, i.e. $p^{n+1}_{0}$ is also independent of time. Using
this in the equation (4.55), we obtain the divergence-free condition for the
velocity $\nabla\cdot{\bf u}_{0}^{n+1}=0$.

This yields the following theorem:

######  Theorem 4.1.

The time-discrete scheme (4.53) is asymptotic preserving in the sense that at
the leading order asymptotic expansion the equations (4.53a), (4.53b) with
$p^{n+1}_{0}$ constant and $\nabla\cdot{\bf u}_{0}^{n+1}=0$, are a consistent
approximation of the incompressible Euler equations (2.6).

Another noticeable feature of the time-discrete scheme (4.53) is that we can
obtain a time-discrete version for the elliptic equation (2.8). We get it by
applying the divergence operator to the momentum equation (4.53b), after some
algebraic manipulations, and making use of the density equation (4.53a). The
resulting relation is

$\nabla\cdot\left(\frac{\rho^{n}_{0}}{\rho_{0}^{n+1}}({\bf
u}_{0}^{n}\cdot\nabla){\bf
u}_{0}^{n}\right)=-\nabla\cdot\left(\frac{1}{\rho_{0}^{n+1}}\nabla
p_{2}^{n+1}\right),$

which is a consistent discretization of Eq. (2.8) because
$\rho_{0}^{n}/\rho_{0}^{n+1}=1+\mathcal{O}(\Delta t)$.

Note that once it is established that $\nabla\cdot{\bf u}_{0}^{n+1}=0$, from
the first equation of (2.4) it follows that if the initial density
$\rho_{0}^{0}$ is constant, then it remains constant for any later time.

######  Remark 4.2.

Note that in [14] the authors proposed a scheme similar to (3.15), but with a
slight difference in the explicit treatment of the density in the energy
equation, i.e.

$\frac{E^{n+1}-E^{n}}{{\Delta
t}}+\nabla\cdot\left(\frac{E^{n}+p^{n}}{\rho^{n}}{\bf q}^{n+1}\right)=0.$

Their scheme is AP, i.e. it is consistent with the limiting equations (2.4) in
the incompressible regime. However, as $\varepsilon\to 0$, they prove that the
divergence-free condition on the velocity is explicitly satisfied up to the
order of the approximation with the CFL condition independent of the Mach-
number $\varepsilon$, i.e. $\nabla\cdot{\bf u}_{0}^{n+1}=\mathcal{O}(\Delta
t$) (see the Proposition in Section 4.2 in [14]). In our case, instead, scheme
(3.15) has the correct discrete divergence-free condition for the leading
order velocity ${\bf u}^{n+1}_{0}$, i.e. $\nabla\cdot{\bf u}_{0}^{n+1}=0$.

###  4.2 AA property

The AP property guarantees only the consistency of the scheme, but in general
the AP property does not guarantee the high order accuracy of IMEX schemes in
the limit for $\varepsilon\ll 1$, i.e. as $\varepsilon\to 0$ the order of
accuracy may degrade. In what follows, we first formally state the definition
of the AA property, and then we recognize that the SA condition is crucial to
guarantee the AA property of our SI-IMEX R-K scheme. Similarly as the AP
property, here we focus on the AA analysis on time discretizations, while
keeping space continuous.

######  Definition 4.3.

A semi-implicit temporal discretization for the Euler system (2.1) is said to
be asymptotically accurate (AA), if it maintains its order of temporal
accuracy for the limiting system (2.4) when $\varepsilon\to 0$.

######  Proposition 4.4.

Consider an SI-IMEX R-K scheme (3.25)-(3.26) of order $r$ applied to system
(2.1) in a bounded domain $\Omega\subset\mathbb{R}$ with zero Neumann
condition. Assume that the implicit part of the IMEX R-K scheme is SA and that
the initial conditions $(\rho^{0}({\bf x}),\rho^{0}({\bf x}){\bf u}^{0}({\bf
x}),p^{0})^{T}$ are well-prepared in the form of (2.9). Let us denote by
$(\rho^{1}({\bf x};\varepsilon)$, $\rho^{1}({\bf x};\varepsilon){\bf
u}^{1}({\bf x};\varepsilon)$, $p^{1}({\bf x};\varepsilon))^{T}$ the numerical
solution after one time step. Then we have:

$\lim_{\varepsilon\to 0}p^{1}({\bf
x};\varepsilon)=p_{*},\quad\lim_{\varepsilon\to 0}\nabla\cdot{\bf u}^{1}({\bf
x};\varepsilon)=0,$ (4.56)

with $p_{*}$ a constant.

Furthermore, let ${\bf V}_{inc}({\bf x},t)=(\rho_{inc}({\bf
x},t),{\rho_{inc}({\bf x},t){\bf u}_{inc}({\bf x},t)},p_{inc}({\bf x},t))^{T}$
be the exact solution of the incompressible Euler equations (2.6) with the
same initial data. Then one has the following one-step error estimate

$\lim_{\varepsilon\to 0}{\bf V}^{1}({\bf x};\varepsilon)={\bf V}_{inc}({\bf
x},\Delta t)+\mathcal{O}(\Delta t^{r+1}),$ (4.57)

i.e., the scheme is AA.

Proof. We consider the first step from $t^{0}=0$ to $t^{1}=\Delta t$ for the
SI-IMEX R-K scheme (3.25)-(3.26) of order $r$ applied to system (2.1) with
well-prepared initial data (4.51):

$\rho^{0}({\bf x})=\rho_{inc}^{0}+\varepsilon^{2}\,\rho_{2}^{0}({\bf x}),\quad
p^{0}({\bf x})=p_{*}+\varepsilon^{2}\,p_{2}^{0}({\bf x}),\quad{\bf u}^{0}({\bf
x})={\bf u}_{inc}^{0}+\mathcal{O}(\varepsilon),$

where $\rho^{0}_{inc}:=\rho_{inc}({\bf x},0)$, ${\bf u}^{0}_{inc}:={\bf
u}_{inc}({\bf x},0)$ and by well prepared assumption we have: $p_{inc}({\bf
x},0):=p_{*}$ constant independent of time and space, and $\nabla\cdot{\bf
u}^{0}_{inc}=0$.

Now we consider a formal $\varepsilon$-expansion of the quantities
$U_{I}^{(i)}=(\rho^{(i)}_{I},{\bf q}^{(i)}_{I},E^{(i)}_{I})^{T}$, and
$U_{E}^{(i)}=(\rho^{(i)}_{E},{\bf q}^{(i)}_{E},E^{(i)}_{E})^{T}$ with ${\bf
q}^{(i)}_{I}=\rho^{(i)}_{I}{\bf u}^{(i)}_{I}$ and ${\bf
q}^{(i)}_{E}=\rho^{(i)}_{E}{\bf u}^{(i)}_{E}$, as an example for the density
and pressure:

$\begin{array}[]{lll}&\rho^{(i)}_{I}=\rho^{(i)}_{0,I}+\varepsilon^{2}\,\rho^{(i)}_{2,I}+...,&\rho^{(i)}_{E}=\rho^{(i)}_{0,E}+\varepsilon^{2}\,\rho^{(i)}_{2,E}+...,\\\\[8.53581pt]
&p^{(i)}_{I}=p^{(i)}_{0,I}+\varepsilon^{2}\,p^{(i)}_{2,I}+...,&p^{(i)}_{E}=p^{(i)}_{0,E}+\varepsilon^{2}\,p^{(i)}_{2,E}+....\end{array}$
(4.58)

In order to prove the theorem, we use the mathematical induction.

  * •

_Asymptotic accuracy for the internal stages $i=1,\cdots,s$_.

Case $i=1$ leads to the same AP analysis for the scheme (3.15) with $\Delta t$
replaced by $a_{11}\Delta t$. To prove the result for $i>1$ onwards, we make
use of the the induction hypothesis, assuming the property holds for $j\leq
i-1$ and prove that it holds for $j=i$. Then for $j=1,\cdots,i-1$ we have

$\displaystyle p_{0,E}^{(j)}=p_{*},\quad
E_{0,E}^{(j)}=\frac{p_{*}}{\gamma-1},\quad\nabla\cdot{\bf u}_{0,I}^{(j)}=0.$
(4.59)

Now we insert the expansions (4.58) into the explicit step in (3.25), and we
get for the energy equation:

$E^{(i)}_{0,E}=E_{inc}^{0}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla\cdot\left(\bar{H}^{(j)}_{0}{\bf
q}_{0,I}^{(j)}\right),$ (4.60)

with $E_{inc}^{0}=p_{*}/(\gamma-1)$ and for $j=1,\cdots,i-1$,

$\displaystyle\bar{H}_{0}^{(j)}=\frac{E^{(j)}_{0,E}+p^{(j)}_{0,E}}{\rho^{(j)}_{0,I}}=\frac{\gamma}{\gamma-1}\frac{p_{*}}{\rho^{(j)}_{0,I}}.$
(4.61)

Now by (4.59), (4.61) and $E^{(i)}_{0,E}=p^{(i)}_{0,E}/(\gamma-1)$, from
(4.60) we obtain

$p^{(i)}_{0,E}=p_{*}-\Delta
t\,\gamma\,p_{*}\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla\cdot{\bf
u}^{(j)}_{0,I}=p_{*},$ (4.62)

that is, $p^{(i)}_{0,E}=p_{*}$ is constant. Then
$E^{(i)}_{0,E}=p_{*}/(\gamma-1)$ is also constant for the stage $i$.

From (3.25), to $\mathcal{O}(1)$ we obtain for the density and momentum
equations

$\rho_{0,E}^{(i)}=\rho^{0}_{inc}\displaystyle-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\nabla\cdot{\bf q}^{(j)}_{0,E},$ (4.63)

and

${\bf q}_{0,E}^{(i)}={\bf q}_{inc}^{0}-\Delta
t\sum_{j=1}^{i-1}\tilde{a}_{ij}\left(\nabla\cdot\left({\rho^{(j)}_{0,E}}\,{\bf
u}^{(j)}_{0,E}\otimes{\bf u}^{(j)}_{0,E}\right)+\nabla p_{2,I}^{(j)}\right),$
(4.64)

with ${\bf q}_{inc}^{0}=(\rho{\bf u})_{inc}^{0}$ and $\nabla p^{(j)}_{0,E}=0$
for $j=1,\cdots,i-1$.

Similarly, inserting expansions (4.58) into (3.44b), up to $\mathcal{O}(1)$ we
get for the intermediate explicit step $\tilde{U}_{0}^{(i)}$ in (3.25):

$\tilde{{\bf q}}_{0}^{(i)}={\bf q}^{0}_{inc}\displaystyle-\Delta
t\sum_{j=1}^{i-1}a_{ij}\left(\nabla\cdot\left({\rho^{(j)}_{0,E}}{\bf
u}^{(j)}_{0,E}\otimes{\bf u}^{(j)}_{0,E}\right)+\nabla p_{2,I}^{(j)}\right),$
(4.65)

where from (3.13) and (4.59), it follows $\nabla\bar{p}^{(j)}_{E}=0$ for
$j=1,\cdots,i-1$. Furthermore, from (3.44a) we have

$\tilde{\rho}_{0}^{(i)}=\rho^{0}_{inc}\displaystyle-\Delta
t\sum_{j=1}^{i-1}a_{ij}\nabla\cdot{\bf q}^{(j)}_{0,E},$ (4.66)

and using (4.61), we get

$\tilde{E}_{0}^{(i)}=E^{0}_{inc}\displaystyle-\Delta
t\sum_{j=1}^{i-1}a_{ij}\nabla\cdot(\bar{H}^{(j)}_{0}{\bf
q}_{0,I}^{(j)})=\frac{p_{*}}{\gamma-1}-\frac{\gamma\,p_{*}}{\gamma-1}\Delta
t\sum_{j=1}^{i-1}a_{ij}\nabla\cdot{\bf u}_{0,I}^{(j)}.$ (4.67)

Thus, by (4.59) we get from (4.67)

$\tilde{E}_{0}^{(i)}=\frac{p_{*}}{\gamma-1}.$ (4.68)

Now from (3.24b) and (4.68), it follows for the energy equation

$E_{0,I}^{(i)}=\tilde{E}^{(i)}_{0}\displaystyle-\Delta
t\,a_{ii}\nabla\cdot(\bar{H}^{(i)}_{0}{\bf
q}_{0,I}^{(i)})=\frac{p_{*}}{\gamma-1}\displaystyle-\Delta
t\,a_{ii}\frac{\gamma\,p_{*}}{\gamma-1}\nabla\cdot{\bf u}_{0,I}^{(i)}.$

Considering the EOS (3.27b) to zeroth order in $\varepsilon$, we get
$E_{0,I}^{(i)}=p_{0,I}^{(i)}/(\gamma-1)$, and we obtain for the pressure

$p_{0,I}^{(i)}=p_{*}+\Delta t\,\gamma\,p_{*}\,a_{ii}\nabla\cdot{\bf
u}_{0,I}^{(i)}.$

Integrating it over spatial bounded domain $\Omega$, and assuming some
boundary conditions (for example, no-slip or periodic,) we first obtain
$p_{0,I}^{(i)}=p_{*}$ and by this we get $\nabla\cdot{\bf u}_{0,I}^{(i)}=0$ at
the stage $i$.

Finally, from (3.24b), considering (4.65) and (4.66), we get for the density
and momentum,

$\rho_{0,I}^{(i)}=\rho^{0}_{inc}\displaystyle-\Delta
t\sum_{j=1}^{i}a_{ij}\nabla\cdot{\bf q}^{(j)}_{0,E},$ (4.69a) ${\bf
q}_{0,I}^{(i)}={\bf q}_{inc}^{0}-\Delta
t\sum_{j=1}^{i}a_{ij}\left(\nabla\cdot\left({\rho^{(j)}_{0,E}}\,{\bf
u}^{(j)}_{0,E}\otimes{\bf u}^{(j)}_{0,E}\right)+\nabla p_{2,I}^{(j)}\right),$
(4.69b)

where from (3.13) and $p^{(i)}_{0,E}=p_{*}$ it follows in the equation of the
momentum $\nabla\bar{p}^{(i)}_{E}=0$ for $i$.

Then equations (4.63), (4.64), 4.69a), (4.69b), with constant limiting
pressure, i.e. $p_{0,E}^{(i)}=p_{0,I}^{(i)}=p_{*}$ and, the divergence free
leading order velocity, i.e., $\nabla\cdot{\bf u}_{0,I}^{(i)}=0$, provide the
discretization of system (2.6) for the internal stage $i$ of SI-IMEX R-K
scheme. This shows that in the limit $\varepsilon\to 0$, the scheme becomes
the same SI-IMEX RK time-discrete scheme for the incompressible Euler
equations (2.6).

  * •

_Asymtotic accuracy for the numerical solution_.

Assuming that SI-IMEX R-K scheme (3.22) is SA, then the numerical solution
coincides with the last internal stage $s$, and then by setting $i=s$, we get

$p^{1}_{0}=p^{(s)}_{0,I}=p_{*},\quad\nabla\cdot{\bf u}_{0}^{1}=\nabla\cdot{\bf
u}_{0,I}^{(s)}=0,$ (4.70)

i.e., we have (4.56).

Now if we denote by ${\bf V}_{inc}({\bf x},t)=(\rho_{inc}({\bf
x},t),\rho_{inc}({\bf x},t){\bf u}_{inc}({\bf x},t),p_{inc}({\bf x},t))^{T}$
the exact solutions of (2.6), with initial data ${\bf V}_{inc}({\bf
x},0)=(\rho^{0}({\bf x}),\rho^{0}({\bf x}){\bf u}^{0}({\bf x}),p^{0})^{T}$,
from equations (4.63), (4.64), (4.69b) and (4.69a) with (4.70), one gets in
the limit case $\varepsilon=0$, a SI-IMEX R-K scheme of order $r$ for the
numerical solutions of equations (2.6), that is, the SI-IMEX R-K
scheme(3.25)-(3.26) of order $r$ is asymptotically accurate (AA), and the
conclusion (4.57) is obtained.

##  5 Numerical Tests

Through an extensive set of 1D and 2D numerical tests, we will show that our
scheme is uniformly stable, effective and can capture the correct asymptotic
limit. We use the 3rd order SI-IMEX(4,4,3) scheme in time (3.2.2), 4th order
compact central difference discretization for 2nd order derivatives in the
elliptic equation (3.47), and 5th order finite difference WENO reconstruction
in space for both $\nabla_{W}$ and $\nabla_{CW}$ in Section 3.3. Overall the
scheme is 4th order in space and 3rd order in time, denoted as “S4T3”. On the
other hand, an “S2T3” scheme refers to using a 2nd order TVD reconstruction
instead of a 5th order WENO reconstruction in the S4T3 scheme. Reference
solutions are computed by a 5th order finite difference WENO scheme with 3rd
order explicit R-K time discretization [39], denoted as ”WENO5RK3”. For
simplicity, we all take $\gamma=1.4$ with an ideal EOS. The time step is
$\Delta t=\text{CFL}\,\Delta x/\Lambda,$ where
$\Lambda=\max_{x}\left(|u|+\min(1/\varepsilon,1)\,c_{s}\right)$ in 1D, and
$\Lambda=\max_{x,y}\left(|u|+|v|+\min(1/\varepsilon,1)\,c_{s}\right)$ in 2D
respectively. $c_{s}=\sqrt{\gamma{p}/{\rho}}$ is the scaled sound speed. We
take $\text{CFL}=0.25$ for all tests.

######  Example 5.1.

(Two colliding acoustic pulses [32, 37].) This problem is defined on the
domain $-L\leq x\leq L=2/\varepsilon$ with periodic boundary condition and
initial data

$\left\\{\begin{array}[]{ll}\rho(x,0)=\rho_{0}+\frac{1}{2}\varepsilon\rho_{1}\left(1-\cos\left(2\pi
x/L\right)\right),&\rho_{0}=0.955,\quad\rho_{1}=2.0;\\\\[2.84526pt]
u(x,0)=\frac{1}{2}u_{0}\,\text{sign}(x)\left(1-\cos\left(2\pi
x/L\right)\right),&u_{0}=2\sqrt{\gamma};\\\\[2.84526pt]
p(x,0)=p_{0}+\frac{1}{2}\varepsilon p_{1}\left(1-\cos\left(2\pi
x/L\right)\right),&p_{0}=1.0,\quad p_{1}=2\gamma.\end{array}\right.$ (5.1)

We first test the order of accuracy for our scheme. Due to $u(x,0)$ in (5.1)
is not smooth enough, in order to observe more than 2nd order accuracy, we
modify it as

$u(x,0)=u_{0}\sin\left(2\pi x/L\right)\left(1-\cos\left(2\pi
x/L\right)\right)/2.$ (5.2)

Others are the same as in (5.1). We take $\varepsilon=10/11$ to avoid order
reduction from the high order IMEX time discretization. We take mesh sizes
with $N_{k}=2^{k}\cdot N_{0}$, where $k=0,1,2,3$ and $N_{0}=40$. Reference
solutions are computed with $N=2560$. The errors are computed by comparing the
numerical solutions for the pressure $p_{E}$ to its reference solution. In
Table 5.1, we show the errors and orders at time $t=0.1$. The order of
convergence is between $4$ and $5$, due to the dominance of spatial errors.

Table 5.1: Example 5.1. Convergence test for the two colliding acoustic pulses problem with initial condition (5.1), and the velocity is replaced by (5.2). $t=0.1$. $\varepsilon=10/11$.  $N$ | 40 | 80 | 160 | 320
---|---|---|---|---
$L^{1}$ error | 1.62E-02 | 9.97E-04 | 3.54E-05 | 1.34E-06
order | – | 4.02 | 4.82 | 4.72

Then we take $\varepsilon=1/11$ with initial condition (5.1) and compute the
solution up to $T=1.63$ by both “S2T3” and “S4T3” schemes. A reference
solution is computed with $N=2200$. In Fig. 5.1, we compare the reference
solution with the numerical ones obtained with $N=22$ grid points at time
$t=1.63$. We can see that the results of the “S4T3” scheme on the very coarse
mesh match the reference solutions better than the “S2T3” scheme.

Figure 5.1: Example 5.1. Two colliding acoustic pulses with initial condition
(5.1). $t=1.63$.

######  Example 5.2.

(1D shock tube problem.) In this example, we consider two 1D shock tube
problems in the compressible regime when the Mach number is of
$\mathcal{O}(1)$. We take the initial data, one is the Sod problem, where

$(\rho,u,p)=(1,0,1),\quad\text{ if }x<0.5;\quad(0.125,0,0.1),\quad\text{
otherwise};$ (5.3)

another is the Lax problem, where

$(\rho,u,p)=(0.445,0.698,3.528),\quad\text{ if
}x<0.5;\quad(0.5,0,0.571),\quad\text{ otherwise},$ (5.4)

both with $\varepsilon=1$ on the domain $x\in[0,1]$. Reflective boundary
conditions are considered and we take $N=50$. We compute the numerical
solution by both “S2T3” and “S4T3” schemes. The results are shown in Fig. 5.2,
and compared to the exact solutions. For both the Sod and Lax shock tube
problems, the results match the exact solutions very well, which show that our
scheme in the moderate Mach regime ($\varepsilon$ of order 1) can capture
strong discontinuities without any observable numerical oscillations.

Figure 5.2: Example 5.2. 1D shock tube problem for $\varepsilon=1$. Top: the
Sod problem at $t=0.2$. Bottom: the Lax problem at $t=0.16$. Mesh grid:
$N=50$. The solid lines are the exact solutions.

######  Example 5.3.

(2D Convergence test.) We take the initial data

$\left\\{\begin{array}[]{rcl}\rho(x,y,0)&=&1+\varepsilon^{2}\,\sin^{2}(2\pi(x+y)),\\\\[2.84526pt]
\rho(x,y,0)\,u(x,y,0)&=&\sin(2\pi(x-y))+\varepsilon^{2}\,\sin(2\pi(x+y)),\\\\[2.84526pt]
\rho(x,y,0)\,v(x,y,0)&=&\sin(2\pi(x-y))+\varepsilon^{2}\,\cos(2\pi(x+y)),\end{array}\right.$
(5.5)

on the domain $\Omega=[0,1]^{2}$ with periodic boundary conditions. Initially
we take $p=\rho^{\gamma}$. We choose $N=N_{x}=N_{y}$ and refine the mesh size
by $N_{k}=2^{k}\cdot N_{0}$, for $k=0,1,2$, with $N_{0}=32$. The numerical
errors are computed by comparing the numerical solutions of momentum $q_{2}$
to the reference solution, which is computed by the “S4T3” scheme on the mesh
$N_{x}=N_{y}=256$. In Table 5.2, we show the errors and orders at time
$t=0.02$. Around 4th order for $\varepsilon=1$ and 5th order for
$\varepsilon=10^{-6}$ are observed. For the intermediate value of
$\varepsilon=10^{-2}$, order reduction is observed, which is a typical
behaviour of high order IMEX schemes for these multi-scale problems.

Table 5.2: Example 5.3. Convergence test for 2D full Euler equations with initial condition (5.5) at $t=0.02$. Errors are computed by comparison with a reference solution.  | $\varepsilon=1$ | $\varepsilon=10^{-2}$ | $\varepsilon=10^{-6}$
---|---|---|---
$N$ |  $L^{1}$ error | order |  $L^{1}$ error | order |  $L^{1}$ error | order
32 | 4.64e-03 | – | 2.45e-03 | – | 7.26e-05 | –
64 | 2.82e-05 | 4.05 | 2.68e-03 | – | 1.79e-06 | 5.34
128 | 1.34e-06 | 4.39 | 1.43e-03 | 0.91 | 4.81e-08 | 5.22

######  Example 5.4.

(2D Riemann problem.) This is a 2D Riemann problem in the high Mach regime. We
take $\varepsilon=1$ and $\Omega=[-1,1]^{2}$. The initial data are defined in
four quadrants, which are modified from Configuration 3 in [33] with four
shocks

$(\rho,u,v,p)=\left\\{\begin{array}[]{lll}(1.5,0,0,1.5),&x\geq 0.8,&y\geq
0.8;\\\\[2.84526pt] (0.5323,1.206,0,0.3)),&x<0.8,&y\geq 0.8;\\\\[2.84526pt]
(0.138,1.206,1.206,0.029),&x<0.8,&y<0.8;\\\\[2.84526pt]
(0.5323,0,1.206,0.3),&x\geq 0.8,&y<0.8.\end{array}\right.$ (5.6)

The second one is the Configuration 5 in [33] with four contact
discontinuities, where the initial data are

$(\rho,u,v,p)=\left\\{\begin{array}[]{lll}(1,-0.75,-0.5,1),&x\geq 0.5,&y\geq
0.5;\\\\[2.84526pt] (2,-0.75,0.5,1),&x<0.5,&y\geq 0.5;\\\\[2.84526pt]
(1,0.75,0.5,1),&x<0.5,&y<0.5;\\\\[2.84526pt] (3,0.75,-0.5,1),&x\geq
0.5,&y<0.5.\end{array}\right.$ (5.7)

We take mesh size $N_{x}\times N_{y}=400\times 400$. For this example, three
different approaches are compared to the reference solution. The first one is
our main approach, which is splitting by taking $\alpha=1$ in (3.11), using
both spatial discretizations $\nabla_{CW}$ and $\nabla_{W}$ described in
Section 3.3 (“A1”); another one uses the same splitting, but replaces
$\nabla_{CW}$ all by $\nabla_{W}$ (“A2”); the third one is no splitting,
namely, we take $\alpha=0$ in (3.11) and only $\nabla_{W}$ is used (“A3”). In
Fig. 5.3, we show the surface plots of the density at $T=0.8$ for the initial
data (5.6). We can see that the solution of A1 is the closest to the reference
solution and no obvious oscillations are observed. For the other two
approaches A2 and A3, the solutions do not perform well, numerical
oscillations can be clearly seen in the middle region. We also show the
cutting plots along two different lines. A1 approach is observed to perform
better than the A2 and A3 approaches. In Fig. 5.4, we present same results for
the 2D Riemann initial data (5.7) at $T=0.23$. Similar observations can be
made.

Figure 5.3: Example 5.4. Surface and cut plots of the density for 2D Riemann
problem (5.6) at $T=0.8$. Mesh grid: $400\times 400$. Top left: reference
solution; top middle: the approach A1, splitting with characteristic-wise
reconstruction $\nabla_{CW}$; bottom left: the approach A2, splitting with
only component-wise reconstruction $\nabla_{W}$; bottom middle: the approach
A3, no splitting. Top right and bottom right are 1D cuts of solutions.

Figure 5.4: Example 5.4. Surface plots of the density for 2D Riemann problem
(5.7) at $T=0.23$. Mesh grid: $400\times 400$. Top left: reference solution;
top middle: the approach A1, splitting with characteristic-wise reconstruction
$\nabla_{CW}$; bottom left: the approach A2, splitting with only component-
wise reconstruction $\nabla_{W}$; bottom middle: the approach A3, no
splitting. Top right and bottom right are 1D cuts of solutions.

######  Example 5.5.

(Gresho vortex [37, 8, 7].) This is the time-dependent rotational Gresho
vortex problem for the full Euler system. Initially a vortex is centered at
$(x_{0},y_{0})=(0.5,0.5)$ with radius $R=0.4$ in the domain $[0,1]^{2}$. The
initial background state is set as:
$\rho_{\infty}=1,\,\,{\bf{u}}_{\infty}=(u_{\infty},0),\,\,p_{\infty}=1,\,\,c_{\infty}=\sqrt{\gamma
p_{\infty}/\rho_{\infty}}=\sqrt{\gamma}.$ The transverse velocity for the
vortex is given by

$\displaystyle u_{\theta}(r)=\left\\{\begin{array}[]{lll}2r/R,&\text{ if
}&0\leq r<R/2,\\\\[2.84526pt] 2(1-r/R),&\text{ if }&R/2\leq
r<R,\\\\[2.84526pt] 0,&\text{ if }&r\geq R,\end{array}\right.$ (5.11)

and the corresponding velocity components are

$u(x,y,0)=u_{\infty}-\frac{y-y_{0}}{r}u_{\theta}(r),\,\,v(x,y,0)=\frac{x-x_{0}}{r}u_{\theta}(r).$

The centrifugal force is balanced by the pressure gradient, so the (scaled)
pressure is given by

$\displaystyle
p(r)=p_{\infty}+\varepsilon^{2}\left\\{\begin{array}[]{lll}2(r/R)^{2}+2-\log
16,&\text{ if }&0\leq r<R/2,\\\\[2.84526pt]
2(r/R)^{2}-4(2r/R+\log(r/R))+6,&\text{ if }&R/2\leq r<R,\\\\[2.84526pt]
0,&\text{ if }&r\geq R,\end{array}\right.$ (5.15)

where $r=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}$. Periodic boundary conditions on both
directions are used and the mesh is $N_{x}\times N_{y}=100\times 100$. The
background velocity is $u_{\infty}=0.1$ moving in the $x$-direction, and we
take $\varepsilon=10^{-1},10^{-2},10^{-6}$. The rotation period is
$2\pi/\omega=R\pi$ if $r\leq R/2$. We take $T=R\pi$ as one rotating period
[7].

We define the ratio between the local and the maximum Mach number as

$M_{\text{ratio}}=\sqrt{[(u-u_{\infty})^{2}+v^{2}]/(\gamma p/\rho)}.$ (5.16)

In Fig. 5.5 (top left), we show the surface plot of $M_{\text{ratio}}$ for the
initial condition with $\varepsilon=10^{-2}$, which is very similar for other
$\epsilon$’s. As time evolves, the solution will rotate while moving in the
$x$-direction, but its shape will be kept. However, numerically the shape will
be damped due to numerical viscosity. This is a standard example to check
whether the numerical viscosity greatly depends on the parameter $\varepsilon$
or not. For better illustration, we show the cutting plots along $x=0.5$ at
three different times: $1/4$ period, $1/2$ period and $1$ full period, and
compare them to the corresponding initial shapes, for
$\varepsilon=10^{-1},10^{-2},10^{-6}$ respectively. In these plots, the vortex
are shifted to the center by $-u_{\infty}t$ periodically. We can see that the
shapes are preserved relatively well for our ”S4T3” scheme. In Fig. 5.5, we
also show the time evolution of the averaged kinetic energy, which is defined
as

$\text{kinetic energy
}=\sum_{i,j}\left[(u(x_{i},y_{j},t)-u_{\infty})^{2}+v(x_{i},y_{j},t)^{2}\right]/N_{x}N_{y}.$
(5.17)

We can see the conservation of kinetic energy is well maintained, and refining
the mesh can greatly improve the conservation of the kinetic energy.

Figure 5.5: Example 5.5, Gresho vortex problem. Top left: $M_{\text{ratio}}$
in (5.16) at $T=R\pi$ on the mesh grid $100\times 100$. Top middle: time
evolution of the kinetic energy (5.17) relative to its initial value on the
mesh grid $100\times 100$. Top right: time evolution of the kinetic energy
(5.17) relative to its initial value on the mesh grid $200\times 200$. Bottom:
cuts along $x=0.5$ at four different times, initial ($t=0$), $1/4$ period
($t=1/4p$), $1/2$ period ($t=1/2p$) and $1$ period ($t=1p$). The results are
shifted to the center by $-u_{\infty}t$ periodically. From left to right:
$\varepsilon=10^{-1}$, $\varepsilon=10^{-2}$, $\varepsilon=10^{-6}$, on the
mesh grid $100\times 100$.

######  Example 5.6.

(Incompressible flow.) Finally we consider two problems in the incompressible
flow regime by taking $\varepsilon=10^{-6}$. One is the shear flow problem on
$[0,2\pi]^{2}$ with

$v(x,y,0)=\delta\cos(x),\quad
u(x,y,0)=\left\\{\begin{array}[]{ll}\tanh((y-\frac{\pi}{2})/\rho),&\rm{if}\,\,y\leq\pi,\\\\[5.69054pt]
\tanh((\frac{3\pi}{2}-y)/\rho),&\rm{if}\,\,y>\pi.\end{array}\right.$ (5.18)

The other is the Kelvin-Helmholtz instability problem [15] on
$[0,4\pi]\times[0,2\pi]$ with

$u(x,y,0)=\cos(y),\quad v(x,y,0)=0.03\sin(0.5x).$ (5.19)

The initial density and pressure for both cases are taken to be $1$. For the
shear flow problem, we run the solution up to $T=6$ on a mesh grid $256\times
256$, while $T=40$ for the Kelvin-Helmholtz instability problem on the mesh
grid $N_{x}\times N_{y}=256\times 128$. The vorticity $\omega=v_{x}-u_{y}$ for
both cases are shown in Fig. 5.6, where $v_{x}$ and $u_{y}$ are discretized by
the 4th order central difference. We observe that it is comparable to the
results for the isentropic case as in [7].

We also show the time evolution of the divergence error $u_{x}+v_{y}$ for the
velocity in Fig. 5.6. The divergence is computed by a 5th order finite
difference WENO reconstruction with zero viscosity, which mimics what we have
done in the numerical scheme. We also compare it with the linear 4th order
central difference discretization. For both cases, the divergence error is
increasing with time. When very fine structures are no longer supported by the
mesh, the divergence error suffers from a sudden increase. We would remark
that when the flow is incompressible or weakly incompressible, without
discontinuities in the initial condition, high order linear reconstructions
would perform better than WENO reconstructions in preserving the divergence
and in resolving solution structures.

Figure 5.6: Example 5.6, the vorticity $\omega=v_{x}-u_{y}$ and time history
of $L^{\infty}$ norm for the divergence $u_{x}+v_{y}$. Left: the shear flow
problem (5.18) at $T=6$ with mesh grid $256\times 256$. Right: the Kelvin-
Helmholtz instability problem (5.19) at $T=40$ with mesh grid $256\times 128$.

##  6 Conclusion

In this paper we present a high-order semi-implicit IMEX RK WENO scheme for
the full compressible Euler equations in the case of all-Mach flows. We
combine the semi-implicit IMEX RK discretization in time with high order
finite difference WENO space discretizations. The EOS is treated in a semi-
implicit manner, therefore requiring the solution of linearized elliptic
systems at each time step. Characteristic reconstructions are proposed in the
semi-implicit framework with a fixed splitting for the pressure. The scheme is
proven to be asymptotic preserving and asymptotically accurate in the
incompressible limit. Numerical tests have demonstrated the efficiency and
effectiveness of our proposed approach.

## References

  * [1] S. Boscarino, R. Bürger, P. Mulet, G. Russo, and L. M. Villada, Linearly implicit IMEX Runge–Kutta methods for a class of degenerate convection-diffusion problems, SIAM Journal on Scientific Computing, 37 (2015), pp. B305–B331.
  * [2] S. Boscarino, R. Bürger, P. Mulet, G. Russo, and L. M. Villada, On linearly implicit IMEX Runge-Kutta methods for degenerate convection-diffusion problems modeling polydisperse sedimentation, Bulletin of the Brazilian Mathematical Society, New Series, 47 (2016), pp. 171–185.
  * [3] S. Boscarino, F. Filbet, and G. Russo, High order semi-implicit schemes for time dependent partial differential equations, Journal of Scientific Computing, 68 (2016), pp. 975–1001.
  * [4] S. Boscarino, P. G. LeFloch, and G. Russo, High-order asymptotic-preserving methods for fully nonlinear relaxation problems, SIAM Journal on Scientific Computing, 36 (2014), pp. A377–A395.
  * [5] S. Boscarino and L. Pareschi, On the asymptotic properties of IMEX Runge–Kutta schemes for hyperbolic balance laws, Journal of Computational and Applied Mathematics, 316 (2017), pp. 60–73.
  * [6] S. Boscarino, L. Pareschi, and G. Russo, Implicit-explicit Runge–Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit, SIAM Journal on Scientific Computing, 35 (2013), pp. A22–A51.
  * [7] S. Boscarino, J.-M. Qiu, G. Russo, and T. Xiong, A high order semi-implicit IMEX WENO scheme for the all-Mach isentropic Euler system, Journal of Computational Physics, 392 (2019), pp. 594–618.
  * [8] S. Boscarino, G. Russo, and L. Scandurra, All mach number second order semi-implicit scheme for the Euler equations of gas dynamics, Journal of Scientific Computing, (2017), pp. 1–35.
  * [9] J. C. Butcher, Numerical Methods for Ordinary Differential Equations, Third Edition, John Wiley & Sons Ltd, 2016.
  * [10] V. Casulli and D. Greenspan, Pressure method for the numerical solution of transient, compressible fluid flows, International Journal for Numerical Methods in Fluids, 4 (1984), pp. 1001–1012.
  * [11] A. J. Chorin, Numerical solution of the Navier-Stokes equations, Mathematics of computation, 22 (1968), pp. 745–762.
  * [12] A. J. Chorin and J. Marsden, A Mathematical Introduction to Fluid Mechanics, Springer Verlag, 1993.
  * [13] B. Cockburn, C. Johnson, C.-W. Shu, and E. Tadmor, Advanced numerical approximation of nonlinear hyperbolic equations: lectures given at the 2nd session of the Centro Internazionale Matematico Estivo (CIME) held in Cetraro, Italy, June 23-28, 1997, Springer, 2006.
  * [14] F. Cordier, P. Degond, and A. Kumbaro, An asymptotic-preserving all-speed scheme for the Euler and Navier-Stokes equations, Journal of Computational Physics, 231 (2012), pp. 5685–5704.
  * [15] N. Crouseilles, M. Mehrenberger, and E. Sonnendrücker, Conservative semi-lagrangian schemes for vlasov equations, Journal of Computational Physics, 229 (2010), pp. 1927–1953.
  * [16] P. Degond and M. Tang, All speed scheme for the low Mach number limit of the isentropic Euler equations, Communications in Computational Physics, 10 (2011), pp. 1–31.
  * [17] S. Dellacherie, Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number, Journal of Computational Physics, 229 (2010), pp. 978–1016.
  * [18] F. Denner, F. Evrard, and B. van Wachem, Conservative finite-volume framework and pressure-based algorithm for flows of incompressible, ideal-gas and real-gas fluids at all speeds, Journal of Computational Physics, 409 (2020), p. 109348.
  * [19] F. Denner, C.-N. Xiao, and B. van Wachem, Pressure-based algorithm for compressible interfacial flows with acoustically-conservative interface discretization, Journal of Computational Physics, 367 (2018), pp. 192–234.
  * [20] G. Dimarco, R. Loubère, and M.-H. Vignal, Study of a new asymptotic preserving scheme for the Euler system in the low Mach number limit, SIAM Journal on Scientific Computing, 39 (2017), pp. A2099–A2128.
  * [21] E. Godlewski and P.-A. Raviart, Numerical Approximation of Hyperbolic Systems of Conservation Laws, Springer, 2014.
  * [22] J.-L. Guermond, P. Minev, and J. Shen, An overview of projection methods for incompressible flows, Computer Methods in Applied Mechanics and Engineering, 195 (2006), pp. 6011–6045.
  * [23] J. Haack, S. Jin, and J.-G. Liu, An all-speed asymptotic-preserving method for the isentropic Euler and Navier-Stokes equations, Communications in Computational Physics, 12 (2012), pp. 955–980.
  * [24] E. Hairer, S. Nørsett, and G. Wanner, Solving ordinary differential equations: Nonstiff problems, vol. 1, Springer Verlag, 1993.
  * [25] E. Hairer and G. Wanner, Solving ordinary differential equations II: stiff and differential algebraic problems, vol. 2, Springer Verlag, 1993.
  * [26] F. Harlow and A. Amdsen, Numerical calculation of almost incompressible flow, Journal of Computational Physics, 3 (1968), pp. 80–93.
  * [27] F. Harlow and A. Amdsen, A numerical fluid dynamics calculation method for all flow speeds, Journal of Computational Physics, 8 (1971), pp. 197–213.
  * [28] S. Jin, Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations, SIAM Journal on Scientific Computing, 21 (1999), pp. 441–454.
  * [29] S. Jin, Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review, Lecture Notes for Summer School on Methods and Models of Kinetic Theory (M&MKT), Porto Ercole (Grosseto, Italy), (2010), pp. 177–216.
  * [30] S. Klainerman and A. Majda, Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids, Communications on pure and applied Mathematics, 34 (1981), pp. 481–524.
  * [31] S. Klainerman and A. Majda, Compressible and incompressible fluids, Communications on Pure and Applied Mathematics, 35 (1982), pp. 629–651.
  * [32] R. Klein, Semi-implicit extension of a Godunov-type scheme based on low Mach number asymptotics. I: One-dimensional flow, J. Comput. Phys., Vol. 121 (1995), pp. pp. 213–237.
  * [33] P. D. Lax and X.-D. Liu, Solution of two-dimensional Riemann problems of gas dynamics by positive schemes, SIAM Journal on Scientific Computing, 19 (1998), pp. 319–340.
  * [34] R. J. LeVeque, Finite volume methods for hyperbolic problems, vol. 31, Cambridge university press, 2002.
  * [35] F. Miczek, F. Röpke, and P. Edelmann, A new numerical solver for flows at various Mach numbers, Astronomy & Astrophysics, Vol. 576 (2015), p. A50.
  * [36] C.-D. Munz, S. Roller, R. Klein, and K. J. Geratz, The extension of incompressible flow solvers to the weakly compressible regime, Computers & Fluids, 32 (2003), pp. 173–196.
  * [37] S. Noelle, G. Bispen, K. R. Arun, M. Lukáčová-Medvidová, and C.-D. Munz, A weakly asymptotic preserving low Mach number scheme for the Euler equations of gas dynamics, SIAM Journal on Scientific Computing, 36 (2014), pp. B989–B1024.
  * [38] L. Pareschi and G. Russo, Implicit-explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxation, Journal of Scientific computing, 25 (2005), pp. 129–155.
  * [39] C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, (1998), pp. 325–432.
  * [40] C.-W. Shu, High order weighted essentially nonoscillatory schemes for convection dominated problems, SIAM review, 51 (2009), pp. 82–126.
  * [41] M. Tang, Second order method for isentropic euler equation in the low mach number regime, Kinetic and Related Models, 5 (2012), pp. 155–184.
  * [42] M. Tavelli and M. Dumbser, A pressure-based semi-implicit space–time discontinuous galerkin method on staggered unstructured meshes for the solution of the compressible navier–stokes equations at all mach numbers, Journal of Computational Physics, 341 (2017), pp. 341–376.
  * [43] R. Temam, Navier-Stokes Equations: Theory and Numerical Analysis, AMS Chelsea Publishing, 1984.
  * [44] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics: a practical introduction, Springer, 2009.
  * [45] E. Turkel, Preconditioned methods for solving the incompressible and low speed compressible equations, Journal of Computational Physics, 72 (1987), pp. 277–298.
  * [46] E. Turkel, A. Fiterman, and B. van Leer, Preconditioning and the limit to the incompressible flow equations, NASA CR-191500, (1993).
  * [47] C. Viozat, Implicit upwind schemes for low Mach number compressible flows, PhD thesis, Inria, 1997.
  * [48] J. Zeifang, J. Schütz, K. Kaiser, A. Beck, M. Lukáčová-Medvidová, and S. Noelle, A novel full-Euler low Mach number IMEX splitting, Communication in Computational Physics, 27 (2020), pp. 292–320.