# Low-Dissipation Central-Upwind Schemes for Compressible Multifluids

Shaoshuai Chu, Alexander Kurganov, and Ruixiao Xin Department of Mathematics
and Shenzhen International Center for Mathematics, Southern University of
Science and Technology, Shenzhen, 518055, China;
<EMAIL_ADDRESS>of Mathematics, Shenzhen International
Center for Mathematics, and Guangdong Provincial Key Laboratory of
Computational Science and Material Design, Southern University of Science and
Technology, Shenzhen, 518055, China<EMAIL_ADDRESS>of
Mathematics, Southern University of Science and Technology, Shenzhen, 518055,
China<EMAIL_ADDRESS>
###### Abstract

We introduce second-order low-dissipation (LD) path-conservative central-
upwind (PCCU) schemes for the one- (1-D) and two-dimensional (2-D) multifluid
systems, whose components are assumed to be immiscible and separated by
material interfaces. The proposed LD PCCU schemes are derived within the flux
globalization based PCCU framework and they employ the LD central-upwind
(LDCU) numerical fluxes. These fluxes have been recently proposed in [A.
Kurganov and R. Xin, J. Sci. Comput., 96 (2023), Paper No. 56] for the single-
fluid compressible Euler equations and we rigorously develop their multifluid
extensions. In order to achieve higher resolution near the material
interfaces, we track their locations and use an overcompressive SBM limiter in
their neighborhoods, while utilizing a dissipative generalized minmod limiter
in the rest of the computational domain. We first develop a second-order
finite-volume LD PCCU scheme and then extend it to the fifth order of accuracy
via the finite-difference alternative weighted essentially non-oscillatory
(A-WENO) framework. We apply the developed schemes to a number of 1-D and 2-D
numerical examples to demonstrate the performance of the new schemes.

Key words: Low-dissipation central-upwind schemes, path-conservative central-
upwind schemes, flux globalization, affine-invariant WENO-Z interpolation,
compressible multifluids.

AMS subject classification: 76M12, 65M08, 76M20, 65M20, 76N30.

##  1 Introduction

In this paper, we focus on the development of highly accurate and conservative
finite-volume methods for compressible multifluids, which are assumed to be
immiscible. The studied two-dimensional (2-D) multifluid system reads as

$\displaystyle\rho_{t}+(\rho u)_{x}+(\rho v)_{y}=0,$ (1.1) $\displaystyle(\rho
u)_{t}+(\rho u^{2}+p)_{x}+(\rho uv)_{y}=0,$ $\displaystyle(\rho v)_{t}+(\rho
uv)_{x}+(\rho v^{2}+p)_{y}=0,$ $\displaystyle
E_{t}+[u(E+p)]_{x}+[v(E+p)]_{y}=0.$

Here, $x$ and $y$ are spatial variables, $t$ is the time, $\rho$ is the
density, $u$ and $v$ are the $x$\- and $y$-velocities, $p$ is the pressure,
and $E$ is the total energy. The system (1.1) is closed through the following
equation of state (EOS) for each of the fluid components:

$p=(\gamma-1)\left[E-\frac{\rho}{2}(u^{2}+v^{2})\right]-\gamma\pi_{\infty},$
(1.2)

where the parameters $\gamma$ and $\pi_{\infty}$ represent the specific heat
ratio and stiffness parameter, respectively. When $\pi_{\infty}\equiv 0$, the
system (1.1)–(1.2) reduces to the ideal gas multicomponent case.

The fluid components can be identified by the variable $\phi$, such as the
specific heat ratio $\gamma$ (or a certain function of $\gamma$), the mass
fraction of the fluid component in the fluid mixture, or a level-set function
designed to track the interfaces between the fluid components; see, e.g., [2,
3, 15, 36, 41, 33, 51, 50, 10] and references therein. The state variable
$\phi$ propagates with the fluid velocity and thus satisfies the following
advection equation:

$\phi_{t}+u\phi_{x}+v\phi_{y}=0.$ (1.3)

The system (1.1)–(1.3) is a nonlinear hyperbolic system of PDEs and thus its
solutions may develop complicated wave structures including shocks,
rarefactions, and contact discontinuities. In the single-fluid regime, that
is, when $\gamma\equiv{\rm Const}$ and $\pi_{\infty}\equiv{\rm Const}$, the
system (1.1)–(1.3) reduces to the Euler equations of gas dynamics, which can
be numerically solved by finite-volume (FV) methods; see, e.g., the monographs
[20, 30, 45] and references therein. However, a straightforward application of
single-fluid FV methods to the multifluid system (1.1)–(1.3) may generate
spurious pressure and velocity oscillations, which typically originate near
the material interface and then spread all over the computational domain; see,
e.g., the review paper [2] and references therein.

In recent years, a variety of FV methods capable of capturing material
interfaces in a non-oscillatory manner have been proposed. A fully
conservative approach was first developed in [42], where the pressure and
velocity remained constant across the material interface. This approach is
robust but may suffer obvious drawbacks when strong shocks pass through the
fluid interface. The quasi-conservative approach was first introduced in [1],
where pressure and velocity non-disturbing condition at an isolated material
interface was introduced to analyze and derive the spatial discretization. The
resulting schemes reduced the numerical oscillations effectively with the help
of a quasi-conservative discretization. There are also many locally
nonconservative approaches designed to prevent pressure/velocity oscillations
by sacrificing the conservation property near material interfaces. The
conservation errors in these approaches are typically small and decay after
the mesh is refined. The pressure-based hybrid algorithms [24, 9] are obtained
by switching from the conservative energy equation to the nonconservative
pressure one near the interfaces. The ghost-cell methods based on the single-
fluid interpolations leading to two different single-fluid numerical fluxes at
the material interfaces (placed at the cell interfaces at each time step) were
introduced in [3, 15]. The interface tracking method [10] is based on the
interpolation between the single-fluid data from both sides of the interface
and ignoring the “mixed” cell data. We note that both the ghost fluid and
interface tracking approaches are very robust in the 1-D case, but their
multidimensional extensions are rather cumbersome. For several high-order WENO
schemes for compressible multifluids, we refer the reader to [13, 19, 22, 38,
39].

In this paper, our objective is to develop highly accurate and non-oscillatory
numerical schemes for the so-called $\gamma$-based multifluid systems studied
in [43, 8], which in the 2-D case read as (1.1)–(1.2) together with the
equations (1.3) for the state variables $\Gamma:=1/(\gamma-1)$ and
$\Pi:=\gamma\pi_{\infty}/(\gamma-1)$, which we recast as follows:

$\Gamma_{t}+(u\Gamma)_{x}+(v\Gamma)_{y}=\Gamma(u_{x}+v_{y}),\quad\Pi_{t}+(u\Pi)_{x}+(v\Pi)_{y}=\Pi(u_{x}+v_{y}).$
(1.4)

The resulting system is nonconservative (in fact, it can be rewritten in the
conservative form, but as it was shown in [43], a nonconservative form is
preferable for designing an accurate numerical method) and the nonconservative
terms on the right-hand side require a special treatment.

We numerically solve the system (1.1)–(1.2), (1.4) and its one-dimensional
(1-D) version by the Riemann-problem-solver-free central-upwind (CU) schemes,
which were originally introduced in [25, 27, 28] for general multidimensional
hyperbolic systems of conservative laws, and then extended to nonconservative
hyperbolic systems in [6], where path-conservative CU (PCCU) schemes were
introduced. The PCCU schemes were extended to the flux globalization framework
allowing to treat a wider variety of nonconservative systems in [4, 5, 7, 26].

The aforementioned PCCU schemes are based on the CU numerical fluxes from [27,
25], which have relatively large numerical dissipation preventing high
resolution of contact waves/material interfaces. In the recent work [29], we
have introduced a new way of reducing the numerical dissipation present in the
CU schemes and introduced the low-dissipation CU (LDCU) schemes. In these
schemes, the dissipation is reduced at the projection step, performed after
the numerical solution is evolved to the new time level. The novel projection
is based on a subcell resolution technique, which introduces several degrees
of freedom that can be utilized to better approximate contact waves and
material interfaces. In [29], we have designed the LDCU schemes for both the
1-D and 2-D single-fluid compressible Euler equations and in this paper, we
extend the LDCU schemes to the $\gamma$-based multifluid models. The extension
is carried out in the flux globalization PCCU framework and results in new
flux globalization based LD PCCU schemes.

We also extend the proposed LD PCCU schemes to the fifth order of accuracy
using the framework of the finite-difference alternative WENO (A-WENO) schemes
developed in [21, 34, 35, 47, 49, 48, 11]. Our new fifth-order schemes are
based on the LD PCCU numerical fluxes, a new, more efficient way to
approximate the high-order A-WENO correction terms (see [12]), and the
recently proposed fifth-order affine-invariant WENO-Z (Ai-WENO-Z)
interpolation [14, 31, 46] applied to the local characteristic variables with
the help of the local characteristic decomposition (LCD).

This paper is organized as follows. §2 is devoted to the 1-D LD PCCU scheme.
In §2.1, we give an overview of the flux globalization based PCCU schemes and
develop such scheme for the $\gamma$-based multifluid model. The LD PCCU
scheme is derived in §2.2, where we prove that the new scheme preserves
constant velocity and pressure across isolated material interfaces (this
ensures lack of pressure/velocity-based oscillations). In §2.3, the second-
order LD PCCU scheme is extended to the fifth-order flux globalization based
LD Ai-WENO PCCU scheme. In §3, we present the 2-D extensions of the new LD
PCCU schemes. In §4, we test the proposed LD PCCU schemes together with their
fifth-order versions on a number of 1-D and 2-D numerical examples. Finally,
in §5, we give some concluding remarks and comments.

##  2 One-Dimensional Algorithms

In this section, we present the new 1-D flux globalization based LD PCCU
schemes for the 1-D compressible multifluids.

###  2.1 Flux Globalization Based Path-Conservative Central-Upwind Schemes

We begin with a brief overview of the flux globalization based PCCU scheme,
which was introduced in [26] for the general nonconservative system

$\bm{U}_{t}+\bm{F}(\bm{U})_{x}=B(\bm{U})\bm{U}_{x},$

which can be rewritten in the following quasi-conservative form:

$\bm{U}_{t}+\bm{K}(\bm{U})_{x}=0,\quad\bm{K}(\bm{U})=\bm{F}(\bm{U})-\bm{R}(\bm{U}),$
(2.1)

where $\bm{U}(\bm{x},t)\in\mathbb{R}^{d}$ is the vector of unknowns,
$\bm{F}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is a flux, $B\in\mathbb{R}^{d\times
d}$,

$\bm{R}(\bm{U}):=\int\limits^{x}_{\widehat{x}}B(\bm{U})\bm{U}_{\xi}(\xi,t)\,{\rm
d}\xi,$

and $\widehat{x}$ is an arbitrary number.

We first introduce a uniform mesh consisting of the finite-volume cells
$C_{j}:=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ of size
$x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}\equiv\Delta x$ centered at
$x_{j}=(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}})/2$, $j=1,\ldots,N$ and set
$\widehat{x}=x_{\frac{1}{2}}$. We assume that at a certain level $t$, an
approximate solution, realized in terms of the cell averages

$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}(t):\approx\frac{1}{\Delta x}\int\limits_{C_{j}}\bm{U}(x,t)\,{\rm
d}x,$

is available. The cell averages $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}$ are then evolved in
time by solving the following system of ODEs (see [26]):

$\frac{{\rm d}}{{\rm d}t}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}=-\frac{\bm{{\cal
K}}_{j+\frac{1}{2}}-\bm{{\cal K}}_{j-\frac{1}{2}}}{\Delta x},$ (2.2)

where $\bm{{\cal K}}_{j+\frac{1}{2}}$ are the CU numerical fluxes

$\bm{{\cal
K}}_{j+\frac{1}{2}}=\frac{a^{+}_{j+\frac{1}{2}}\bm{K}^{-}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\bm{K}^{+}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}}+\frac{a^{+}_{j+\frac{1}{2}}a^{-}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}}\left(\bm{U}^{+}_{j+\frac{1}{2}}-\bm{U}^{-}_{j+\frac{1}{2}}\right).$
(2.3)

Notice that all of the indexed quantities in (2.2) and (2.3) as well as other
indexed quantities introduced below depend on $t$, but from now on we will
omit this dependence for the sake of brevity. In (2.3),
$\bm{U}^{\pm}_{j+\frac{1}{2}}$ are the left/right-sided point values of
$\bm{U}$ at the cell interfaces $x_{j+\frac{1}{2}}$, which are computed using
a piecewise linear reconstruction, which will be discussed in §2.1.1, and
$a^{\pm}_{j+\frac{1}{2}}$ are the one-sided local speeds of propagation, which
can be estimated using the largest and smallest eigenvalues of the matrix
$\frac{\partial\bm{F}}{\partial\bm{U}}(\bm{U})-B(\bm{U})$. The global fluxes
$\bm{K}^{\pm}_{j+\frac{1}{2}}$ are obtained using the relation in (2.1),
namely, by

$\bm{K}^{\pm}_{j+\frac{1}{2}}=\bm{F}^{\pm}_{j+\frac{1}{2}}-\bm{R}^{\pm}_{j+\frac{1}{2}},$
(2.4)

where $\bm{F}^{\pm}_{j+\frac{1}{2}}:=\bm{F}(\bm{U}^{\pm}_{j+\frac{1}{2}})$ and
the point values of the global variable $\bm{R}$ are computed as follows.
First, we set $\bm{R}^{-}_{\frac{1}{2}}:=\bm{0}$ and then evaluate

$\bm{R}^{+}_{\frac{1}{2}}=\bm{B}_{\bm{\Psi},{\frac{1}{2}}},$ (2.5)

and recursively

$\bm{R}^{-}_{j+\frac{1}{2}}=\bm{R}^{+}_{j-\frac{1}{2}}+\bm{B}_{j},\quad\bm{R}^{+}_{j+\frac{1}{2}}=\bm{R}^{-}_{j+\frac{1}{2}}+\bm{B}_{\bm{\Psi},{j+\frac{1}{2}}},\quad
j=1,\ldots,N.$ (2.6)

In (2.5) and (2.6), $\bm{B}_{j}$ and $\bm{B}_{\bm{\Psi},{j+\frac{1}{2}}}$ are
obtained using a proper quadrature for the integrals in

$\bm{B}_{j}\approx\int\limits_{C_{j}}B(\bm{U})\bm{U}_{x}\,{\rm
d}x\quad\mbox{and}\quad\bm{B}_{\bm{\Psi},{j+\frac{1}{2}}}=\int\limits^{1}_{0}B\big{(}\bm{\Psi}_{j+\frac{1}{2}}(s)\big{)}\bm{\Psi}^{\prime}_{j+\frac{1}{2}}(s)\,{\rm
d}s,$ (2.7)

where,
$\bm{\Psi}_{j+\frac{1}{2}}(s):=\bm{\Psi}\big{(}s;\bm{U}^{-}_{j+\frac{1}{2}},\bm{U}^{+}_{j+\frac{1}{2}}\big{)}$
is a sufficiently smooth path connecting the states
$\bm{U}^{-}_{j+\frac{1}{2}}$ and $\bm{U}^{+}_{j+\frac{1}{2}}$, that is,

$\bm{\Psi}:[0,1]\times\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{d},\quad\bm{\Psi}(0;\bm{U}^{-}_{j+\frac{1}{2}},\bm{U}^{+}_{j+\frac{1}{2}})=\bm{U}^{-}_{j+\frac{1}{2}},\quad\bm{\Psi}(1;\bm{U}^{-}_{j+\frac{1}{2}},\bm{U}^{+}_{j+\frac{1}{2}})=\bm{U}^{+}_{j+\frac{1}{2}}.$

####  2.1.1 Application to the Compressible Multifluid System

We apply the flux globalization based PCCU scheme to the 1-D $\gamma$-based
multifluid system

$\displaystyle\rho_{t}+(\rho u)_{x}=0,$ (2.8) $\displaystyle(\rho u)_{t}+(\rho
u^{2}+p)_{x}=0,$ $\displaystyle E_{t}+[u(E+p)]_{x}=0,$
$\displaystyle\Gamma_{t}+(u\Gamma)_{x}=\Gamma u_{x},$
$\displaystyle\Pi_{t}+(u\Pi)_{x}=\Pi u_{x},$

completed with the following EOS:

$p=(\gamma-1)\left[E-\frac{\rho}{2}u^{2}\right]-\gamma\pi_{\infty}.$ (2.9)

The system (2.8) can be rewritten in the equivalent quasi-conservative form
(2.1) with

$\displaystyle\bm{U}:=(\rho,\rho
u,E,\Gamma,\Pi)^{\top},\quad\bm{F}(\bm{U})=(\rho u,\rho
u^{2}+p,u(E+p),u\Gamma,u\Pi)^{\top},$ (2.10)
$\displaystyle\mbox{and}\quad\bm{R}(\bm{U})=\bigg{(}0,0,0,\int\limits^{x}_{\widehat{x}}\Gamma
u_{\xi}\,{\rm d}\xi,\int\limits^{x}_{\widehat{x}}\Pi u_{\xi}\,{\rm
d}\xi\bigg{)}^{\top}.$

We first discuss the reconstruction procedure for recovering the point values
$\bm{U}^{\pm}_{j+\frac{1}{2}}$ out of the cell averages $\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}$.
Since the variables $u$ and $p$ are continuous across material interfaces
(contact waves), we reconstruct the primitive variables
$\bm{V}:=(\rho,u,p,\Gamma,\Pi)^{\top}$ instead of the conservative ones. To
this end, we compute the cell centered values of $u$ and $p$,

$u_{j}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})_{j}}{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}_{j}},\quad p_{j}=\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j}}\left[\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}_{j}-\frac{((\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern
0.0pt}}})_{j})^{2}}{2\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j}}-\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern
0.0pt}}}_{j}\right],$ (2.11)

and then construct the linear pieces

$\widetilde{\bm{V}}_{j}(x)=\bm{V}_{j}+(\bm{V}_{x})_{j}(x-x_{j}),\quad x\in
C_{j},$ (2.12)

where $\bm{V}_{j}:=(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}_{j},u_{j},p_{j},\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}_{j},\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern
0.0pt}}}_{j})^{\top}$ and the slopes $(\bm{V}_{x})_{j}$ are supposed to be
computed with the help of a nonlinear limiter to ensure a non-oscillatory
nature of (2.12). In the numerical experiments reported in §4, we implement a
simple adaptive limiting strategy and use different limiters near and away
from the material interfaces. To this end, we need to detect the location of
the interfaces. In the two-fluid case, this can be done as follows. We first
introduce $\widehat{\Gamma}:=(\Gamma_{\rm I}+\Gamma_{\rm II})/2$, where
$\Gamma_{\rm I}$ and $\Gamma_{\rm II}$ are the values of $\Gamma$ for the
first and second fluid, respectively. We then assume that the interface is
located either in cell $C_{j}$ or $C_{j+1}$ if

$(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j}-\widehat{\Gamma})(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j+1}-\widehat{\Gamma})<0.$ (2.13)

In these two cells as well as in the neighboring cells $C_{j-1}$ and
$C_{j+2}$, we use the overcompressive SBM limiter [32]:

$(\bm{V}_{x})_{j}=\phi^{\rm SBM}_{\theta,\tau}\left(\frac{\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j+1}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern 0.0pt}}}_{j}}{\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j-1}}\right)\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern 0.0pt}}}_{j+1}-\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j}}{\Delta x},$ (2.14)

where the two-parameter function

$\phi^{\rm SBM}_{\theta,\tau}(r):=\left\\{\begin{aligned}
&0&&\mbox{if}~{}r<0,\\\ &\min\\{r\theta,1+\tau(r-1)\\}&&\mbox{if}~{}0<r\leq
1,\\\ &r\phi^{\rm
SBM}_{\theta,\tau}\Big{(}\frac{1}{r}\Big{)}&&\mbox{otherwise},\end{aligned}\right.$
(2.15)

is applied in a componentwise manner with $\tau=-0.5$, which belongs to the
overcompressive range of values of $\tau$; see [32]. Away from the material
interfaces, we use a dissipative generalized minmod limiter which is given by
the same formulae (2.14)–(2.15), but with $\tau=0.5$. In both areas, we use
$\theta=1.3$.

###### Remark 2.1

Notice that the generalized minmod limiter can be written in a simpler form
[32, 37, 44] by

$(\bm{V}_{x})_{j}={\rm minmod}\left(\theta\,\frac{\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j+1}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern 0.0pt}}}_{j}}{\Delta
x},\,\frac{\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j+1}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern 0.0pt}}}_{j-1}}{2\Delta
x},\,\theta\,\frac{\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern 0.0pt}}}_{j}-\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j-1}}{\Delta x}\right),$

with the minmod function defined by

${\rm minmod}(c_{1},c_{2},\ldots)=\left\\{\begin{aligned}
&\min(c_{1},c_{2},\ldots)&&{\rm if}~{}c_{i}>0,~{}\forall i,\\\
&\max(c_{1},c_{2},\ldots)&&{\rm if}~{}c_{i}<0,~{}\forall i,\\\
&0&&\mbox{otherwise}.\end{aligned}\right.$

###### Remark 2.2

If the number of fluid components is more than two, detecting interface cells
becomes a more complicated task. A reasonable extension of the strategy used
in (2.13) should be developed for the problem at hand.

Equipped with $(\bm{V}_{x})_{j}$, we then use (2.12) to obtain

$\bm{V}^{-}_{j+\frac{1}{2}}=\lim_{x\to
x^{-}_{j+\frac{1}{2}}}\widetilde{\bm{V}}(x)=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j}+\frac{\Delta
x}{2}(\bm{V}_{x})_{j},\quad\bm{V}^{+}_{j+\frac{1}{2}}=\lim_{x\to
x^{+}_{j+\frac{1}{2}}}\widetilde{\bm{V}}(x)=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j+1}-\frac{\Delta x}{2}(\bm{V}_{x})_{j+1},$ (2.16)

and the corresponding point values of the conservative variables $\bm{U}$:

${\bm{U}}^{\pm}_{j+\frac{1}{2}}=\left(\rho^{\pm}_{j+\frac{1}{2}},\rho^{\pm}_{j+\frac{1}{2}}u^{\pm}_{j+\frac{1}{2}},E^{\pm}_{j+\frac{1}{2}},\Gamma^{\pm}_{j+\frac{1}{2}},\Pi^{\pm}_{j+\frac{1}{2}}\right)^{\top},\quad
E^{\pm}_{j+\frac{1}{2}}=\Gamma^{\pm}_{j+\frac{1}{2}}p^{\pm}_{j+\frac{1}{2}}+\frac{\rho^{\pm}_{j+\frac{1}{2}}}{2}\big{(}u^{\pm}_{j+\frac{1}{2}}\big{)}^{2}+\Pi^{\pm}_{j+\frac{1}{2}}.$
(2.17)

We then compute the point values $\bm{K}^{\pm}_{j+\frac{1}{2}}$ following
(2.4)–(2.6). First, we use (2.10) and (2.17) to obtain

$\bm{F}^{\pm}_{j+\frac{1}{2}}=\left(\rho^{\pm}_{j+\frac{1}{2}}u^{\pm}_{j+\frac{1}{2}},\rho^{\pm}_{j+\frac{1}{2}}\big{(}u^{\pm}_{j+\frac{1}{2}}\big{)}^{2}+p^{\pm}_{j+\frac{1}{2}},u^{\pm}_{j+\frac{1}{2}}\big{(}E^{\pm}_{j+\frac{1}{2}}+p^{\pm}_{j+\frac{1}{2}}\big{)},\Gamma^{\pm}_{j+\frac{1}{2}}u^{\pm}_{j+\frac{1}{2}},\Pi^{\pm}_{j+\frac{1}{2}}u^{\pm}_{j+\frac{1}{2}}\right)^{\top},$
(2.18)

and then evaluate $\bm{B}_{j}$ in (2.7) by substituting there the piecewise
linear reconstructions (2.12) of $u$, $\Gamma$, and $\Pi$, which results in

$\bm{B}_{j}=\bigg{(}0,0,0,\frac{\Gamma^{-}_{j+\frac{1}{2}}+\Gamma^{+}_{j-\frac{1}{2}}}{2}\big{(}u^{-}_{j+\frac{1}{2}}-u^{+}_{j-\frac{1}{2}}\big{)},\frac{\Pi^{-}_{j+\frac{1}{2}}+\Pi^{+}_{j-\frac{1}{2}}}{2}\big{(}u^{-}_{j+\frac{1}{2}}-u^{+}_{j-\frac{1}{2}}\big{)}\bigg{)}^{\top}.$
(2.19)

Next, in order to obtain $\bm{B}_{\bm{\Psi},{j+\frac{1}{2}}}$ in (2.7), a
proper path connecting the states
$(u^{-}_{j+\frac{1}{2}},\Gamma^{-}_{j+\frac{1}{2}},\Pi^{-}_{j+\frac{1}{2}})$
and
$(u^{+}_{j+\frac{1}{2}},\Gamma^{+}_{j+\frac{1}{2}},\Pi^{+}_{j+\frac{1}{2}})$
needs to be used, for instance, a simple linear path:

$\displaystyle\bm{\Psi}^{u}_{j+\frac{1}{2}}(s)=u^{-}_{j+\frac{1}{2}}+s\big{(}u^{+}_{j+\frac{1}{2}}-u^{-}_{j+\frac{1}{2}}\big{)},~{}~{}\bm{\Psi}^{\Gamma}_{j+\frac{1}{2}}(s)=\Gamma^{-}_{j+\frac{1}{2}}+s\big{(}\Gamma^{+}_{j+\frac{1}{2}}-\Gamma^{-}_{j+\frac{1}{2}}\big{)},$
(2.20)
$\displaystyle\bm{\Psi}^{\Pi}_{j+\frac{1}{2}}(s)=\Pi^{-}_{j+\frac{1}{2}}+s\big{(}\Pi^{+}_{j+\frac{1}{2}}-\Pi^{-}_{j+\frac{1}{2}}\big{)}.$

Substituting (2.20) into (2.7) then results in

$\bm{B}_{\bm{\Psi},{j+\frac{1}{2}}}=\bigg{(}0,0,0,\frac{\Gamma^{+}_{j+\frac{1}{2}}+\Gamma^{-}_{j+\frac{1}{2}}}{2}\big{(}u^{+}_{j+\frac{1}{2}}-u^{-}_{j+\frac{1}{2}}\big{)},\frac{\Pi^{+}_{j+\frac{1}{2}}+\Pi^{-}_{j+\frac{1}{2}}}{2}\big{(}u^{+}_{j+\frac{1}{2}}-u^{-}_{j+\frac{1}{2}}\big{)}\bigg{)}^{\top}.$
(2.21)

Finally, the one-sided local-speeds of propagation $a^{\pm}_{j+\frac{1}{2}}$
can be estimated by

$a^{+}_{j+\frac{1}{2}}=\max\left\\{u^{-}_{j+\frac{1}{2}}+c^{-}_{j+\frac{1}{2}},u^{+}_{j+\frac{1}{2}}+c^{+}_{j+\frac{1}{2}},0\right\\},\quad
a^{-}_{j+\frac{1}{2}}=\min\left\\{u^{-}_{j+\frac{1}{2}}-c^{-}_{j+\frac{1}{2}},u^{+}_{j+\frac{1}{2}}-c^{+}_{j+\frac{1}{2}},0\right\\},$

where $c:=\sqrt{\left[(1+\Gamma)p+\Pi\right]/(\Gamma\rho)}$.

###### Remark 2.3

The numerical fluxes (2.3) are slightly different from those present in [26]
as one of the goals in [26] was to make the resulting scheme well-balanced and
thus different values of $\bm{U}^{\pm}_{j+\frac{1}{2}}$ were used in (2.3) and
in the computation of $\bm{K}^{\pm}_{j+\frac{1}{2}}$.

###  2.2 Flux Globalization Based Low-Dissipation PCCU Schemes

In this section, we derived a modified, LD version of the flux globalization
based PCCU scheme presented in §2.1. We follow the idea used in [29], where
the LDCU scheme for the single-fluid compressible Euler equations has been
introduced. In order to extend the LDCU scheme to the multifluid case, we now
go through all of the derivation steps and begin with the development of the
fully discrete LD PCCU scheme.

####  2.2.1 Fully Discrete Scheme

We assume that the computed cell averages $\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,n}_{j}:\approx\frac{1}{\Delta x}\int_{C_{j}}\bm{U}(x,t^{n})\,{\rm
d}x$ are available at a certain time level $t=t^{n}$ and use them to
reconstruct a second-order piecewise linear interpolant consisting of the
linear pieces $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,n}_{j}+(\bm{U}_{x})^{n}_{j}(x-x_{j})$, $x\in C_{j}$, where the
slopes $(\bm{U}_{x})^{n}_{j}$ are obtained using a certain nonlinear limiter.
We then estimate the local speeds of propagation $a^{\pm}_{j+\frac{1}{2}}$,
introduce the corresponding points
$x^{n}_{{j+\frac{1}{2}},\ell}:=x_{j+\frac{1}{2}}+a^{-}_{j+\frac{1}{2}}\Delta
t^{n}$ and
$x^{n}_{{j+\frac{1}{2}},r}:=x_{j+\frac{1}{2}}+a^{+}_{j+\frac{1}{2}}\Delta
t^{n}$, and integrate the system (2.1) over the space-time control volumes,
which consist of the “smooth”,
$[x_{{j-\frac{1}{2}},r},x_{{j+\frac{1}{2}},\ell}]\times[t^{n},t^{n+1}]$, and
“nonsmooth”,
$[x_{{j+\frac{1}{2}},\ell},x_{{j+\frac{1}{2}},r}]\times[t^{n},t^{n+1}]$, ones,
where $t^{n+1}:=t^{n}+\Delta t^{n}$. This way the solution is evolved in time
and upon the completion of the evolution step, we obtain the intermediate cell
averages

$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\rm\,int}_{j+\frac{1}{2}}=\frac{1}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}}\Big{\\{}$
$\displaystyle\bm{U}^{n}_{{j+\frac{1}{2}},r}a^{+}_{j+\frac{1}{2}}-\frac{(\bm{U}_{x})^{n}_{j+1}}{2}\big{(}a^{+}_{j+\frac{1}{2}}\big{)}^{2}\Delta
t^{n}-\bm{U}^{n}_{{j+\frac{1}{2}},\ell}a^{-}_{j+\frac{1}{2}}+\frac{(\bm{U}_{x})^{n}_{j}}{2}\big{(}a^{-}_{j+\frac{1}{2}}\big{)}^{2}\Delta
t^{n}$ (2.22)
$\displaystyle-\left[\bm{K}\big{(}\bm{U}^{n+{\frac{1}{2}}}_{{j+\frac{1}{2}},r}\big{)}-\bm{K}\big{(}\bm{U}^{n+{\frac{1}{2}}}_{{j+\frac{1}{2}},\ell}\big{)}\right]\Big{\\}}$

and

$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\rm\,int}_{j}$
$\displaystyle=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,n}_{j}+\frac{(\bm{U}_{x})^{n}_{j}}{2}\big{(}a^{+}_{j-\frac{1}{2}}+a^{-}_{j+\frac{1}{2}}\big{)}\Delta
t^{n}$ (2.23) $\displaystyle-\frac{\Delta t^{n}}{\Delta
x-\big{(}a^{+}_{j-\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\big{)}\Delta
t^{n}}\left[\bm{K}\big{(}\bm{U}^{n+{\frac{1}{2}}}_{{j+\frac{1}{2}},\ell}\big{)}-\bm{K}\big{(}\bm{U}^{n+{\frac{1}{2}}}_{{j-\frac{1}{2}},r}\big{)}\right];$

see [29] for details. In (2.22)–(2.23), the point values of $\bm{U}$ at
$\big{(}x^{n}_{{j+\frac{1}{2}},\ell},t^{n}\big{)}$ and
$\big{(}x^{n}_{{j+\frac{1}{2}},r},t^{n}\big{)}$ are computed using the
piecewise linear reconstruction of $\bm{U}$, namely,

$\displaystyle\bm{U}^{n}_{{j+\frac{1}{2}},\ell}:=\widetilde{\bm{U}}(x^{n}_{{j+\frac{1}{2}},\ell},t^{n})=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,n}_{j}+(\bm{U}_{x})^{n}_{j}\Big{(}\frac{\Delta
x}{2}+a^{-}_{j+\frac{1}{2}}\Delta t^{n}\Big{)},$
$\displaystyle\bm{U}^{n}_{{j+\frac{1}{2}},r}:=\widetilde{\bm{U}}(x^{n}_{{j+\frac{1}{2}},r},t^{n})=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,n}_{j+1}-(\bm{U}_{x})^{n}_{j+1}\Big{(}\frac{\Delta
x}{2}-a^{+}_{j+\frac{1}{2}}\Delta t^{n}\Big{)},$

and the point values of $\bm{U}$ at
$\big{(}x^{n}_{{j+\frac{1}{2}},\ell},t^{n+{\frac{1}{2}}}\big{)}$ and
$\big{(}x^{n}_{{j+\frac{1}{2}},r},t^{n+{\frac{1}{2}}})$, are obtained using
the Taylor expansions about $\big{(}x_{{j+\frac{1}{2}},r},t^{n}\big{)}$ and
$\big{(}x_{{j+\frac{1}{2}},\ell},t^{n}\big{)}$, respectively, and the fact
that $\bm{U}_{t}=-\bm{K}(\bm{U})_{x}$:

$\bm{U}^{n+{\frac{1}{2}}}_{{j+\frac{1}{2}},\ell}=\bm{U}^{n}_{{j+\frac{1}{2}},\ell}-\frac{\Delta
t^{n}}{2}\bm{K}\big{(}\bm{U}^{n}_{{j+\frac{1}{2}},\ell}\big{)}_{x},\quad\bm{U}^{n+{\frac{1}{2}}}_{{j+\frac{1}{2}},r}=\bm{U}^{n}_{{j+\frac{1}{2}},r}-\frac{\Delta
t^{n}}{2}\bm{K}\big{(}\bm{U}^{n}_{{j+\frac{1}{2}},r}\big{)}_{x}.$ (2.24)

Here, the slopes $\bm{K}\big{(}\bm{U}^{n}_{{j+\frac{1}{2}},\ell}\big{)}_{x}$
and $\bm{K}\big{(}\bm{U}^{n}_{{j+\frac{1}{2}},r}\big{)}_{x}$ can be computed
with the help of a certain nonlinear limiter; see [25] for details.

Next, the intermediate solution, realized in terms of $\\{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j}\\}$ and $\\{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}\\}$, is projected onto the original grid. To this end, we
need to construct the interpolant

$\widetilde{\bm{U}}^{\,\rm int}(x)=\sum_{j}\left\\{\widetilde{\bm{U}}^{\,\rm
int}_{j+\frac{1}{2}}(x){\cal
X}_{[x_{{j+\frac{1}{2}},\ell},x_{{j+\frac{1}{2}},r}]}+\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j}{\cal X}_{[x_{{j-\frac{1}{2}},r},x_{{j+\frac{1}{2}},\ell}]}\right\\},$
(2.25)

where $\cal X$ denotes the characteristic function of the corresponding
intervals. We set

$\widetilde{\bm{U}}^{\,\rm int}_{j+\frac{1}{2}}(x)=\left\\{\begin{aligned}
&\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}},&&x<x_{j+\frac{1}{2}},\\\ &\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}},&&x>x_{j+\frac{1}{2}},\end{aligned}\right.$ (2.26)

where the values $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$ and $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$ are determined in several steps. First, according to
the local conservation requirement, the conditions

$\displaystyle a^{+}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}=(a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}})\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}$ (2.27)

have to be satisfied, and the rest of the relations on $\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$ and $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$ are to be established for the problem at hand. In
fact, the conservation of $\Gamma$ and $\Pi$ components of $\bm{U}$ is not
physically essential, but the relation (2.27) for $\Gamma$ and $\Pi$ are
crucial for proving physically relevant properties of the resulting flux
globalization based LD PCCU scheme; see §2.2.3.

For the $\gamma$-based multifluid model (2.8)–(2.9), we follow the single-
fluid approach from [29] and make the projection step sharp and accurate for
the contact waves, which are linearly degenerate and thus affected by the
excessive numerical dissipation much more than nonlinear shock waves. In order
to design such projection step, we consider an isolated contact wave
consisting of the jump discontinuities in $\rho$, $\gamma$, and $\Pi$
propagating in the region with constant $u$ and $p$, and make both $u$ and $p$
to be constant across the cell interface, namely, we set

$\displaystyle\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,L}_{j+\frac{1}{2}}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,R}_{j+\frac{1}{2}}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}},$ (2.28) $\displaystyle\frac{1}{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}}\bigg{(}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\frac{\big{(}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,L}_{j+\frac{1}{2}}\big{)}^{2}}{2\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}\bigg{)}=\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}}\bigg{(}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}-\frac{\big{(}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,R}_{j+\frac{1}{2}}\big{)}^{2}}{2\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}\bigg{)},$

where we have used the EOS (2.9). Next, we solve (2.27) and (2.28) for
$(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho
u$\kern 0.0pt}}})^{\rm int,L}_{j+\frac{1}{2}}$, $(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,R}_{j+\frac{1}{2}}$, $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$, and $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$, and express these quantities in terms of
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$, $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$, $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$, $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$, $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$, and $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$:

$\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,L}_{j+\frac{1}{2}}=\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j+\frac{1}{2}},\quad\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}=\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm int}_{j+\frac{1}{2}}$
(2.29) $\displaystyle\hskip
36.98866pt+\frac{a^{+}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}\big{)}}{2\big{(}a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\big{)}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}\big{)}^{2}}\big{(}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j+\frac{1}{2}}\big{)}^{2}+\frac{a^{+}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}\big{)}}{\big{(}a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\big{)}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}},$ $\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,R}_{j+\frac{1}{2}}=\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j+\frac{1}{2}},\quad\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}=\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm int}_{j+\frac{1}{2}}$
$\displaystyle\hskip
36.98866pt+\frac{a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}\big{)}}{2\big{(}a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\big{)}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}\big{)}^{2}}\big{(}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j+\frac{1}{2}}\big{)}^{2}+\frac{a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}\big{)}}{\big{(}a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\big{)}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}.$

Notice that after enforcing (2.27) and (2.28), we are left with three degrees
of freedom, which we use to make the profiles of $\rho$, $\Gamma$, and $\Pi$
across the cell interface as sharp as possible and, at the same time, non-
oscillatory. This is achieved in the same way as in [29], namely, by setting

$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+\frac{\delta^{\rho}_{j+\frac{1}{2}}}{a^{-}_{j+\frac{1}{2}}},\quad\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+\frac{\delta^{\Gamma}_{j+\frac{1}{2}}}{a^{-}_{j+\frac{1}{2}}},\quad\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+\frac{\delta^{\Pi}_{j+\frac{1}{2}}}{a^{-}_{j+\frac{1}{2}}},$
(2.30) $\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+\frac{\delta^{\rho}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}},\quad\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+\frac{\delta^{\Gamma}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}},\quad\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+\frac{\delta^{\Pi}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}},$

where

$\displaystyle\delta^{\rho}_{j+\frac{1}{2}}={\rm
minmod}\left(-a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\rho^{\rm
int}_{{j+\frac{1}{2}},\ell}\big{)},a^{+}_{j+\frac{1}{2}}\big{(}\rho^{\rm
int}_{{j+\frac{1}{2}},r}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}\big{)}\right),$ (2.31)
$\displaystyle\delta^{\Gamma}_{j+\frac{1}{2}}={\rm
minmod}\left(-a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\Gamma^{\rm
int}_{{j+\frac{1}{2}},\ell}\big{)},a^{+}_{j+\frac{1}{2}}\big{(}\Gamma^{\rm
int}_{{j+\frac{1}{2}},r}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}\big{)}\right),$
$\displaystyle\delta^{\Pi}_{j+\frac{1}{2}}={\rm
minmod}\left(-a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\Pi^{\rm
int}_{{j+\frac{1}{2}},\ell}\big{)},a^{+}_{j+\frac{1}{2}}\big{(}\Pi^{\rm
int}_{{j+\frac{1}{2}},r}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}\big{)}\right),$

and $\bm{U}^{\,\rm
int}_{{j+\frac{1}{2}},\ell}\approx\bm{U}(x_{{j+\frac{1}{2}},\ell},t^{n+1})$
and $\bm{U}^{\,\rm
int}_{{j+\frac{1}{2}},r}\approx\bm{U}(x_{{j+\frac{1}{2}},r},t^{n+1})$ are
evaluated similarly to (2.24):

${\bm{U}}^{\,\rm
int}_{{j+\frac{1}{2}},\ell}=\bm{U}^{n}_{{j+\frac{1}{2}},\ell}-\Delta
t^{n}\bm{K}\big{(}\bm{U}^{n}_{{j+\frac{1}{2}},\ell}\big{)}_{x},\quad{\bm{U}}^{\,\rm
int}_{{j+\frac{1}{2}},r}=\bm{U}^{n}_{{j+\frac{1}{2}},r}-\Delta
t^{n}\bm{K}\big{(}\bm{U}^{n}_{{j+\frac{1}{2}},r}\big{)}_{x}.$ (2.32)

We then complete the construction of $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$ and $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$ by substituting (2.30) into (2.29), which results in

$\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,L}_{j+\frac{1}{2}}=(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j+\frac{1}{2}}+\frac{\delta^{\rho}_{j+\frac{1}{2}}}{a^{-}_{j+\frac{1}{2}}}\,u^{\rm
int}_{j+\frac{1}{2}},\quad(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,R}_{j+\frac{1}{2}}=(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j+\frac{1}{2}}+\frac{\delta^{\rho}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}}\,u^{\rm
int}_{j+\frac{1}{2}},$ (2.33) $\displaystyle\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}=\bigg{(}1+\frac{\delta^{\Gamma}_{j+\frac{1}{2}}}{a^{-}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}\bigg{)}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+\frac{\delta^{\rho}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\delta^{\Gamma}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}{2a^{-}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}\big{(}u^{\rm
int}_{j+\frac{1}{2}}\big{)}^{2}+\frac{\delta^{\Pi}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\delta^{\Gamma}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}{a^{-}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}},$ $\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}=\bigg{(}1+\frac{\delta^{\Gamma}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}\bigg{)}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+\frac{\delta^{\rho}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\delta^{\Gamma}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}{2a^{+}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}\big{(}u^{\rm
int}_{j+\frac{1}{2}}\big{)}^{2}+\frac{\delta^{\Pi}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\delta^{\Gamma}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}},$

where $u^{\rm int}_{j+\frac{1}{2}}:=(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j+\frac{1}{2}}/\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}$.

Equipped with (2.30)–(2.33), we finalize the projection step by integrating
the piecewise constant interpolant (2.25)–(2.26) over the cell $C_{j}$. This
leads to the following cell averages at the time level $t=t^{n+1}$:

$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,n+1}_{j}$
$\displaystyle=\frac{1}{\Delta x}\int\limits_{C_{j}}\widetilde{\rho}^{\,\rm
int}(x)\,{\rm d}x=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\rm\,int}_{j}+\frac{\Delta
t^{n}}{\Delta x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j-\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}\Big{]}$ (2.34)
$\displaystyle\stackrel{{\scriptstyle(\ref{2.2.10})}}{{=}}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}^{\rm\,int}_{j}+\frac{\Delta t^{n}}{\Delta
x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}+\delta^{\rho}_{j-\frac{1}{2}}-\delta^{\rho}_{j+\frac{1}{2}}\Big{]},$
$\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{n+1}_{j}$
$\displaystyle=\frac{1}{\Delta x}\int\limits_{C_{j}}(\widetilde{\rho u})^{\rm
int}(x)\,{\rm d}x$ (2.35) $\displaystyle=(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j}+\frac{\Delta t^{n}}{\Delta
x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,R}_{j-\frac{1}{2}}-(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int,L}_{j+\frac{1}{2}}-(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j}\big{)}\Big{]}$
$\displaystyle\stackrel{{\scriptstyle(\ref{2.2.13})}}{{=}}(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j}+\frac{\Delta t^{n}}{\Delta
x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j-\frac{1}{2}}-(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{j+\frac{1}{2}}-(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm int}_{j}\big{)}$
$\displaystyle\hskip 99.58464pt+\delta^{\rho}_{j-\frac{1}{2}}u^{\rm
int}_{j-\frac{1}{2}}-\delta^{\rho}_{j+\frac{1}{2}}u^{\rm
int}_{j+\frac{1}{2}}\Big{]},$ $\displaystyle\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,n+1}_{j}$
$\displaystyle=\frac{1}{\Delta x}\int\limits_{C_{j}}\widetilde{E}^{\,\rm
int}(x)\,{\rm d}x=\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\rm int}_{j}+\frac{\Delta
t^{n}}{\Delta x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,R}_{j-\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm int}_{j}\big{)}\Big{]}$
(2.36)
$\displaystyle\stackrel{{\scriptstyle(\ref{2.2.13})}}{{=}}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}^{\rm\,int}_{j}+\frac{\Delta t^{n}}{\Delta
x}\bigg{[}a^{+}_{j-\frac{1}{2}}\big{(}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm int}_{j}\big{)}$
$\displaystyle\hskip
52.63777pt+\frac{\delta^{\rho}_{j-\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\delta^{\Gamma}_{j-\frac{1}{2}}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}}{2\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}}\big{(}u^{\rm
int}_{j-\frac{1}{2}}\big{)}^{2}-\frac{\delta^{\rho}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\delta^{\Gamma}_{j+\frac{1}{2}}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}{2\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}\big{(}u^{\rm int}_{j+\frac{1}{2}}\big{)}^{2}$
$\displaystyle+\frac{\delta^{\Pi}_{j-\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}-\delta^{\Gamma}_{j-\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}}-\frac{\delta^{\Pi}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\delta^{\Gamma}_{j+\frac{1}{2}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}+\frac{\delta^{\Gamma}_{j-\frac{1}{2}}}{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}-\frac{\delta^{\Gamma}_{j+\frac{1}{2}}}{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}\bigg{]},$ $\displaystyle\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}^{\,n+1}_{j}$ $\displaystyle=\frac{1}{\Delta
x}\int\limits_{C_{j}}\widetilde{\Gamma}^{\,\rm int}(x)\,{\rm
d}x=\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}^{\rm\,int}_{j}+\frac{\Delta t^{n}}{\Delta
x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{j-\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}\Big{]}$ (2.37)
$\displaystyle\stackrel{{\scriptstyle(\ref{2.2.10})}}{{=}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}^{\rm\,int}_{j}+\frac{\Delta t^{n}}{\Delta
x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}+\delta^{\Gamma}_{j-\frac{1}{2}}-\delta^{\Gamma}_{j+\frac{1}{2}}\Big{]},$
$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,n+1}_{j}$
$\displaystyle=\frac{1}{\Delta x}\int\limits_{C_{j}}\widetilde{\Pi}^{\,\rm
int}(x)\,{\rm d}x=\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\rm\,int}_{j}+\frac{\Delta
t^{n}}{\Delta x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,R}_{j-\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}\Big{]}$ (2.38)
$\displaystyle\stackrel{{\scriptstyle(\ref{2.2.10})}}{{=}}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern
0.0pt}}}^{\rm\,int}_{j}+\frac{\Delta t^{n}}{\Delta
x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}+\delta^{\Pi}_{j-\frac{1}{2}}-\delta^{\Pi}_{j+\frac{1}{2}}\Big{]}.$

####  2.2.2 Semi-Discrete Scheme

Finally, we pass to the semi-discrete limit $\Delta t^{n}\to 0$ in
(2.34)–(2.38) and proceed as in [29, §2.2] to end up with the semi-
discretization (2.2) with the modified (compared with (2.3)) numerical fluxes

$\bm{{\cal
K}}_{j+\frac{1}{2}}=\frac{a^{+}_{j+\frac{1}{2}}\bm{K}^{-}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\bm{K}^{+}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}}+\frac{a^{+}_{j+\frac{1}{2}}a^{-}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}}\left(\bm{U}^{+}_{j+\frac{1}{2}}-\bm{U}^{-}_{j+\frac{1}{2}}\right)+\bm{q}_{j+\frac{1}{2}},$
(2.39)

where

$\displaystyle\bm{q}_{j+\frac{1}{2}}$
$\displaystyle=q^{\rho}_{j+\frac{1}{2}}\Big{(}1,u^{*}_{j+\frac{1}{2}},\frac{\big{(}u^{*}_{j+\frac{1}{2}}\big{)}^{2}}{2},0,0\Big{)}^{\top}$
(2.40)
$\displaystyle+q^{\Gamma}_{j+\frac{1}{2}}\bigg{(}0,0,\frac{1}{\Gamma^{*}_{j+\frac{1}{2}}}\Big{[}E^{*}_{j+\frac{1}{2}}-\frac{\big{(}(\rho
u)^{*}_{j+\frac{1}{2}}\big{)}^{2}}{2\rho^{*}_{j+\frac{1}{2}}}-\Pi^{*}_{j+\frac{1}{2}}\Big{]},1,0\bigg{)}^{\top}+q^{\Pi}_{j+\frac{1}{2}}(0,0,1,0,1)^{\top}$

is a built-in “anti-diffusion” term. In (2.40),

$\displaystyle\bm{U}^{*}_{j+\frac{1}{2}}=\frac{a^{+}_{j+\frac{1}{2}}\bm{U}^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\bm{U}^{-}_{j+\frac{1}{2}}-\big{(}\bm{K}^{+}_{j+\frac{1}{2}}-\bm{K}^{-}_{j+\frac{1}{2}}\big{)}}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}},\quad
u^{*}_{j+\frac{1}{2}}=\frac{(\rho
u)^{*}_{j+\frac{1}{2}}}{\rho^{*}_{j+\frac{1}{2}}},$ (2.41) $\displaystyle
q^{\rho}_{j+\frac{1}{2}}={\rm
minmod}\left(-a^{-}_{j+\frac{1}{2}}\big{(}\rho^{*}_{j+\frac{1}{2}}-\rho^{-}_{j+\frac{1}{2}}\big{)},a^{+}_{j+\frac{1}{2}}\big{(}\rho^{+}_{j+\frac{1}{2}}-\rho^{*}_{j+\frac{1}{2}}\big{)}\right),$
$\displaystyle q^{\Gamma}_{j+\frac{1}{2}}={\rm
minmod}\left(-a^{-}_{j+\frac{1}{2}}\big{(}\Gamma^{*}_{j+\frac{1}{2}}-\Gamma^{-}_{j+\frac{1}{2}}\big{)},a^{+}_{j+\frac{1}{2}}\big{(}\Gamma^{+}_{j+\frac{1}{2}}-\Gamma^{*}_{j+\frac{1}{2}}\big{)}\right),$
$\displaystyle q^{\Pi}_{j+\frac{1}{2}}={\rm
minmod}\left(-a^{-}_{j+\frac{1}{2}}\big{(}\Pi^{*}_{j+\frac{1}{2}}-\Pi^{-}_{j+\frac{1}{2}}\big{)},a^{+}_{j+\frac{1}{2}}\big{(}\Pi^{+}_{j+\frac{1}{2}}-\Pi^{*}_{j+\frac{1}{2}}\big{)}\right).$

###### Remark 2.4

If we replace (2.26) with the following limited linear piece:

$\widetilde{\bm{U}}^{\,\rm int}_{j+\frac{1}{2}}(x)=\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+(\bm{U}_{x})^{\rm
int}_{j+\frac{1}{2}}\left(x-\frac{x_{{j+\frac{1}{2}},r}+x_{{j+\frac{1}{2}},\ell}}{2}\right),$

where

$(\bm{U}_{x})^{\rm int}_{j+\frac{1}{2}}={\rm
minmod}\Bigg{(}\,\frac{\bm{U}^{\rm
int}_{{j+\frac{1}{2}},r}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}}{\delta},\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\bm{U}^{\rm
int}_{{j+\frac{1}{2}},\ell}}{\delta}\Bigg{)},\quad\delta:=\frac{\Delta
t^{n}}{2}\big{(}a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\big{)},$

we will end up with an alternative built-in “anti-diffusion” term

$\bm{q}_{j+\frac{1}{2}}=-\frac{a^{+}_{j+\frac{1}{2}}a^{-}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}}\,{\rm
minmod}\left(\bm{U}^{*}_{j+\frac{1}{2}}-\bm{U}^{-}_{j+\frac{1}{2}},\bm{U}^{+}_{j+\frac{1}{2}}-\bm{U}^{*}_{j+\frac{1}{2}}\right),$

which was introduced in [25]. This will result in another flux globalization
based PCCU scheme for the $\gamma$-based multifluid system. In the numerical
results reported in §4, we will compare the behavior of this PCCU scheme with
the proposed LD PCCU one.

####  2.2.3 Properties of the Semi-Discrete Scheme

In this section, we establish two important properties the designed semi-
discrete scheme satisfies:

1\. When initially $\gamma\equiv{\rm Const}$ and $\pi_{\infty}\equiv{\rm
Const}$, that is, if initially the system contains a single fluid, then
$\gamma$ and $\pi_{\infty}$ will stay constant for all $t$;

2\. At the isolated material interface at which initially $u\equiv{\rm Const}$
and $p\equiv{\rm Const}$, both $u$ and $p$ will stay constant for all $t$.

To this end, we prove the following theorem.

###### Theorem 2.5

1\. If $\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j}\equiv\widehat{\Gamma}={\rm Const}$ and $\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern
0.0pt}}}_{j}\equiv\widehat{\Pi}={\rm Const}$ for all $j$ at a certain time
level $t$, then

$\frac{{\rm d}}{{\rm d}t}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}_{j}=\frac{{\rm d}}{{\rm
d}t}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}_{j}\equiv 0,\quad\forall j.$

2\. If at a certain time level $t=t^{n}$, $u^{n}_{j}\equiv\widehat{u}={\rm
Const}$ and $p^{n}_{j}\equiv\widehat{p}={\rm Const}$ for all $j$, then at the
next time level $t=t^{n+1}$,

$u^{n+1}_{j}\equiv\widehat{u}\quad{\rm and}\quad
p^{n+1}_{j}\equiv\widehat{p},\quad\forall j,$ (2.42)

provided the ODE system (2.2), (2.39)–(2.41) is discretized using the forward
Euler method.

Proof: 1\. We only show that $\frac{{\rm d}}{{\rm d}t}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j}\equiv 0$ as $\frac{{\rm d}}{{\rm d}t}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}_{j}\equiv
0$ can be proved in a similar way. Since the point values of $\Gamma$ at the
cell interfaces are obtained using the piecewise linear reconstruction (2.16),
we have
$\Gamma^{+}_{j+\frac{1}{2}}=\Gamma^{-}_{j+\frac{1}{2}}=\Gamma^{+}_{j-\frac{1}{2}}\equiv\widehat{\Gamma}$,
which we substitute into (2.18), (2.19), and (2.21) to obtain

$\big{(}F^{(4)}\big{)}^{\pm}_{j+\frac{1}{2}}=\widehat{\Gamma}u^{\pm}_{j+\frac{1}{2}},\quad
B^{(4)}_{j}=\widehat{\Gamma}\big{(}u^{-}_{j+\frac{1}{2}}-u^{+}_{j-\frac{1}{2}}\big{)},\quad
B^{(4)}_{\bm{\Psi},{j+\frac{1}{2}}}=\widehat{\Gamma}\big{(}u^{+}_{j+\frac{1}{2}}-u^{-}_{j+\frac{1}{2}}\big{)}.$
(2.43)

We then use (2.4)–(2.6) to compute the flux differences

$\displaystyle\big{(}K^{(4)}\big{)}^{+}_{j+\frac{1}{2}}-\big{(}K^{(4)}\big{)}^{-}_{j+\frac{1}{2}}=\big{(}F^{(4)}\big{)}^{+}_{j+\frac{1}{2}}-\big{(}F^{(4)}\big{)}^{-}_{j+\frac{1}{2}}-B^{(4)}_{\bm{\Psi},{j+\frac{1}{2}}}\stackrel{{\scriptstyle(\ref{pr1.1})}}{{=}}0,$
(2.44)
$\displaystyle\big{(}K^{(4)}\big{)}^{-}_{j+\frac{1}{2}}-\big{(}K^{(4)}\big{)}^{+}_{j-\frac{1}{2}}=\big{(}F^{(4)}\big{)}^{-}_{j+\frac{1}{2}}-\big{(}F^{(4)}\big{)}^{+}_{j-\frac{1}{2}}-B^{(4)}_{j}\stackrel{{\scriptstyle(\ref{pr1.1})}}{{=}}0,$

which we substitute into (2.41) to verify that

$q^{\Gamma}_{j+\frac{1}{2}}=0.$ (2.45)

Finally, we substitute (2.44)–(2.45) into (2.2), (2.39)–(2.40) to end up with
$\frac{{\rm d}}{{\rm d}t}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}_{j}\equiv 0$. This
completes the proof of the first part of the theorem.

2\. Since the point values of $u$ and $p$ at the cell interfaces are computed
using the piecewise linear reconstruction (2.16), we have
$u^{+}_{j+\frac{1}{2}}=u^{-}_{j+\frac{1}{2}}=u^{+}_{j-\frac{1}{2}}\equiv\widehat{u}$
and
$p^{+}_{j+\frac{1}{2}}=p^{-}_{j+\frac{1}{2}}=p^{+}_{j-\frac{1}{2}}\equiv\widehat{p}$
at the time level $t=t^{n}$. This results in

$\displaystyle\bm{U}^{\pm}_{j+\frac{1}{2}}=\big{(}\rho^{\pm}_{j+\frac{1}{2}},\rho^{\pm}_{j+\frac{1}{2}}\widehat{u},E^{\pm}_{j+\frac{1}{2}},\Gamma^{\pm}_{j+\frac{1}{2}},\Pi^{\pm}_{j+\frac{1}{2}}\big{)}^{\top},\quad
E^{\pm}_{j+\frac{1}{2}}=\widehat{p}\,\Gamma^{\pm}_{j+\frac{1}{2}}+\frac{\widehat{u}^{\,2}}{2}\rho^{\pm}_{j+\frac{1}{2}}+\Pi^{\pm}_{j+\frac{1}{2}},$
(2.46)
$\displaystyle\bm{F}^{\pm}_{j+\frac{1}{2}}=\big{(}\rho^{\pm}_{j+\frac{1}{2}}\widehat{u},\rho^{\pm}_{j+\frac{1}{2}}\widehat{u}^{\,2}+\widehat{p},\widehat{u}\,(E^{\pm}_{j+\frac{1}{2}}+\widehat{p}\,),\Gamma^{\pm}_{j+\frac{1}{2}}\widehat{u},\Pi^{\pm}_{j+\frac{1}{2}}\widehat{u}\big{)}^{\top},$

and hence, after substituting $u^{\pm}_{j+\frac{1}{2}}$ into (2.19) and
(2.21), we obtain

$\bm{B}_{j}=\bm{B}_{\bm{\Psi},{j+\frac{1}{2}}}\equiv\bm{0}.$

The latter equality implies
$\bm{K}^{\pm}_{j+\frac{1}{2}}=\bm{F}^{\pm}_{j+\frac{1}{2}}$, which together
with (2.46) results in

$\bm{F}^{+}_{j+\frac{1}{2}}-\bm{F}^{-}_{j+\frac{1}{2}}=\widehat{u}\,(\bm{U}^{+}_{j+\frac{1}{2}}-\bm{U}^{+}_{j+\frac{1}{2}}),$

so that the first line in (2.41) can be rewritten as

$\bm{U}^{*}_{j+\frac{1}{2}}=\frac{a^{+}_{j+\frac{1}{2}}\bm{U}^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\bm{U}^{-}_{j+\frac{1}{2}}-\widehat{u}\,\big{(}\bm{U}^{+}_{j+\frac{1}{2}}-\bm{U}^{-}_{j+\frac{1}{2}}\big{)}}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}},\quad
u^{*}_{j+\frac{1}{2}}=\widehat{u}.$ (2.47)

We then use (2.40), (2.41), (2.46), and (2.47) to compute the “anti-diffusion”
term, which reduces to

$\bm{q}_{j+\frac{1}{2}}=\Big{(}q^{\rho}_{j+\frac{1}{2}},\widehat{u}\,q^{\rho}_{j+\frac{1}{2}},\widehat{p}\,q^{\Gamma}_{j+\frac{1}{2}}+\frac{\widehat{u}^{\,2}}{2}\,q^{\rho}_{j+\frac{1}{2}}+q^{\Pi}_{j+\frac{1}{2}},q^{\Gamma}_{j+\frac{1}{2}},q^{\Pi}_{j+\frac{1}{2}}\Big{)}^{\top},$
(2.48)

and then substitute (2.46) and (2.48) into (2.39) to obtain the numerical
fluxes

$\bm{{\cal K}}_{j+\frac{1}{2}}=\Big{(}{\cal
K}^{(1)}_{j+\frac{1}{2}},\widehat{u}\,{\cal
K}^{(1)}_{j+\frac{1}{2}}+\widehat{p},\widehat{p}\,{\cal
K}^{(4)}_{j+\frac{1}{2}}+\frac{\widehat{u}^{\,2}}{2}{\cal
K}^{(1)}_{j+\frac{1}{2}}+{\cal
K}^{(5)}_{j+\frac{1}{2}}+\widehat{u}\widehat{p},{\cal
K}^{(4)}_{j+\frac{1}{2}},{\cal K}^{(5)}_{j+\frac{1}{2}}\Big{)}^{\top},$

which, in turn, are substituted into (2.2) to obtain the following semi-
discrete relations:

$\frac{{\rm d}}{{\rm d}t}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern
0.0pt}}})^{n}_{j}=\widehat{u}\,\frac{{\rm d}}{{\rm d}t}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}^{\,n}_{j},\quad\frac{{\rm d}}{{\rm d}t}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}^{\,n}_{j}=\widehat{p}\,\frac{{\rm d}}{{\rm d}t}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}^{\,n}_{j}+\frac{\widehat{u}^{\,2}}{2}\frac{{\rm d}}{{\rm
d}t}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,n}_{j}+\frac{{\rm
d}}{{\rm d}t}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,n}_{j}.$

These relations are discretized using the forward Euler method, which results
in

$\displaystyle\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern
0.0pt}}})^{n+1}_{j}-(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{n}_{j}}{\Delta
t}=\widehat{u}\,\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}^{\,n+1}_{j}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,n}_{j}}{\Delta t},$
(2.49) $\displaystyle\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,n+1}_{j}-\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}^{\,n}_{j}}{\Delta t}=\widehat{p}\,\frac{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}^{\,n+1}_{j}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,n}_{j}}{\Delta
t}+\frac{\widehat{u}^{\,2}}{2}\cdot\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}^{\,n+1}_{j}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,n}_{j}}{\Delta
t}+\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern
0.0pt}}}^{\,n+1}_{j}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,n}_{j}}{\Delta t}.$

Finally, we substitute (2.11) expressed at both time levels $t=t^{n}$ and
$t=t^{n+1}$ into (2.49) to end up with (2.42). This completes the proof of the
second part of the theorem.$\hfill\blacksquare$

###### Remark 2.6

The second part of Theorem 2.5 is still true if the forward Euler time
discretization is replaced with another strong stability preserving (SSP) ODE
solver; see, e.g.,[16, 17].

###### Remark 2.7

As in [29], the computation of numerical fluxes in (2.39) should be
desingularized to avoid division by zero or very small numbers. If
$a^{+}_{j+\frac{1}{2}}<\varepsilon_{0}$ and
$a^{-}_{j+\frac{1}{2}}>-\varepsilon_{0}$ for a small positive
$\varepsilon_{0}$, we replace the fluxes $\bm{{\cal K}}_{j+\frac{1}{2}}$ with

$\bm{{\cal
K}}_{j+\frac{1}{2}}=\frac{\bm{K}\big{(}\bm{U}^{-}_{j+\frac{1}{2}}\big{)}+\bm{K}\big{(}\bm{U}^{+}_{j+\frac{1}{2}}\big{)}}{2}.$

In all of the numerical examples reported in §4, we have taken
$\varepsilon_{0}=10^{-12}$.

###  2.3 Flux Globalization Based LD Ai-WENO PCCU Scheme

In this section, we extend the second-order flux globalization based LD PCCU
schemes from §2.2 to the fifth order of accuracy within the A-WENO framework.

The semi-discrete fifth-order LD Ai-WENO PCCU scheme for the 1-D quasi-
conservative system (2.1) reads as

$\frac{{\rm d}}{{\rm
d}t}\bm{U}_{j}=-\frac{\bm{\mathcal{H}}_{j+\frac{1}{2}}-\bm{\mathcal{H}}_{j-\frac{1}{2}}}{\Delta
x},$ (2.50)

where $\bm{U}_{j}:\approx\bm{U}(x_{j},t)$ and the fifth-order numerical fluxes
$\bm{\mathcal{H}}_{j+\frac{1}{2}}$ are defined by

$\bm{\mathcal{H}}_{j+\frac{1}{2}}=\bm{{\cal K}}_{j+\frac{1}{2}}-\frac{\Delta
x}{24}(\bm{K}_{xx})_{j+\frac{1}{2}}+\frac{7(\Delta
x)^{3}}{5760}(\bm{K}_{xxxx})_{j+\frac{1}{2}}.$ (2.51)

Here, $\bm{{\cal K}}_{j+\frac{1}{2}}$ are the finite-volume fluxes
(2.39)–(2.41), and $(\bm{K}_{xx})_{j+\frac{1}{2}}$ and
$({\bm{K}_{xxxx}})_{j+\frac{1}{2}}$ are approximations of the second- and
fourth-order spatial derivatives of $\bm{K}$ at $x=x_{j+\frac{1}{2}}$, which
we compute using the finite-difference approximations recently proposed in
[12]:

$\displaystyle(\bm{K}_{xx})_{j+\frac{1}{2}}=\frac{1}{12(\Delta
x)^{2}}\left[-\bm{{\cal K}}_{j-\frac{3}{2}}+16\bm{{\cal
K}}_{j-\frac{1}{2}}-30\bm{{\cal K}}_{j+\frac{1}{2}}+16\bm{{\cal
K}}_{j+\frac{3}{2}}-\bm{{\cal K}}_{j+\frac{5}{2}}\right],$ (2.52)
$\displaystyle(\bm{K}_{xxxx})_{j+\frac{1}{2}}=\frac{1}{(\Delta
x)^{4}}\left[\bm{{\cal K}}_{j-\frac{3}{2}}-4\bm{{\cal
K}}_{j-\frac{1}{2}}+6\bm{{\cal K}}_{j+\frac{1}{2}}-4\bm{{\cal
K}}_{j+\frac{3}{2}}+\bm{{\cal K}}_{j+\frac{5}{2}}\right].$

The resulting semi-discrete scheme (2.50)–(2.52) will be fifth-order accurate
provided the point values $\bm{U}^{\pm}_{j+\frac{1}{2}}$ are calculated using
a fifth-order interpolation. To this end, we apply the recently proposed
fifth-order Ai-WENO-Z interpolation [14, 31, 46], which we briefly describe in
Appendix A.

##  3 Two-Dimensional Algorithms

In this section, we extend the proposed 1-D flux globalization based LD A-WENO
PCCU schemes to the 2-D $\gamma$-based multifluid system (1.1), (1.2), (1.4).
This system can be written in the vector form

$\bm{U}_{t}+\bm{F}(\bm{U})_{x}+\bm{G}(\bm{U})_{y}=B(\bm{U})\bm{U}_{x}+C(\bm{U})\bm{U}_{y},$

or, equivalently, in the quasi-conservative form

$\bm{U}_{t}+\bm{K}(\bm{U})_{x}+\bm{L}(\bm{U})_{y}=\bm{0}$ (3.1)

with

$\displaystyle\bm{K}(\bm{U})=\bm{F}(\bm{U})-\bm{R}(\bm{U}),$
$\displaystyle\bm{L}(\bm{U})=\bm{G}(\bm{U})-\bm{S}(\bm{U}),$ (3.2)
$\displaystyle\bm{R}(\bm{U})=\int\limits^{x}_{\widehat{x}}B(\bm{U})\bm{U}_{\xi}\,{\rm
d}\xi,$
$\displaystyle\bm{S}(\bm{U})=\int\limits^{y}_{\widehat{y}}C(\bm{U})\bm{U}_{\eta}\,{\rm
d}\eta.$

Here,

$\displaystyle\bm{U}:=(\rho,\rho u,\rho v,E,\Gamma,\Pi)^{\top},$ (3.3)
$\displaystyle\bm{F}(\bm{U})=(\rho u,\rho u^{2}+p,\rho
uv,u(E+p),u\Gamma,u\Pi)^{\top},$ $\displaystyle
B(\bm{U})\bm{U}_{x}=\big{(}0,0,0,0,\Gamma u_{x},\Pi u_{x}\big{)}^{\top},$
$\displaystyle\bm{G}(\bm{U})=(\rho v,\rho uv,\rho
v^{2}+p,v(E+p),v\Gamma,v\Pi)^{\top},$ $\displaystyle
C(\bm{U})\bm{U}_{y}=\big{(}0,0,0,0,\Gamma v_{y},\Pi v_{y}\big{)}^{\top}.$

We first introduce a uniform mesh consisting of the finite-volume cells
$C_{j,k}:=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]\times[y_{k-\frac{1}{2}},y_{k+\frac{1}{2}}]$
of the uniform size $\Delta x\Delta y$ with
$x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}\equiv\Delta x$ and
$y_{k+\frac{1}{2}}-y_{k-\frac{1}{2}}\equiv\Delta y$ centered at
$(x_{j},y_{k})$ with $x_{j}=(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}})/2$ and
$(y_{k-\frac{1}{2}}+y_{k+\frac{1}{2}})/2$, $j=1,\ldots,N_{x}$,
$k=1,\ldots,N_{y}$.

We assume that at certain time level $t$, an approximate solution, realized in
terms of the cell averages $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j,k}:\approx\frac{1}{\Delta x\Delta
y}\iint_{C_{j,k}}\bm{U}(x,y,t)\,{\rm d}x\,{\rm d}y$, is available. These cell
averages are then evolved in time by solving the following system of ODEs:

$\frac{{\rm d}}{{\rm d}t}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}=-\frac{\bm{{\cal
K}}_{{j+\frac{1}{2}},k}-\bm{{\cal K}}_{{j-\frac{1}{2}},k}}{\Delta
x}-\frac{\bm{{\cal L}}_{j,{k+\frac{1}{2}}}-\bm{{\cal
L}}_{j,{k-\frac{1}{2}}}}{\Delta y},$ (3.4)

where the $x$\- and $y$-numerical fluxes are

$\displaystyle\bm{{\cal
K}}_{{j+\frac{1}{2}},k}=\frac{a^{+}_{{j+\frac{1}{2}},k}\bm{K}^{-}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}\bm{K}^{+}_{{j+\frac{1}{2}},k}}{a^{+}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}}+\frac{a^{+}_{{j+\frac{1}{2}},k}a^{-}_{{j+\frac{1}{2}},k}}{a^{+}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}}\left(\bm{U}^{+}_{{j+\frac{1}{2}},k}-\bm{U}^{-}_{{j+\frac{1}{2}},k}\right)+\bm{q}_{{j+\frac{1}{2}},k},$
(3.5) $\displaystyle\bm{{\cal
L}}_{j,{k+\frac{1}{2}}}=\frac{b^{+}_{j,{k+\frac{1}{2}}}\bm{L}^{-}_{j,{k+\frac{1}{2}}}-b^{-}_{j,{k+\frac{1}{2}}}\bm{L}^{+}_{j,{k+\frac{1}{2}}}}{b^{+}_{j,{k+\frac{1}{2}}}-b^{-}_{j,{k+\frac{1}{2}}}}+\frac{b^{+}_{j,{k+\frac{1}{2}}}b^{-}_{j,{k+\frac{1}{2}}}}{b^{+}_{j,{k+\frac{1}{2}}}-b^{-}_{j,{k+\frac{1}{2}}}}\left(\bm{U}^{+}_{j,{k+\frac{1}{2}}}-\bm{U}^{-}_{j,{k+\frac{1}{2}}}\right)+\bm{q}_{j,{k+\frac{1}{2}}}.$
(3.6)

The one-sided point values $\bm{U}^{\pm}_{{j+\frac{1}{2}},k}$ and
$\bm{U}^{\pm}_{j,{k+\frac{1}{2}}}$ at the cell interfaces
$(x_{j+\frac{1}{2}},y_{k})$ and $(x_{j},y_{k+\frac{1}{2}})$, respectively, are
obtained as follows. We first use the cell averages $\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}$
to compute the point values of $u$, $v$, and $p$ at the cell centers:

$u_{j,k}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern
0.0pt}}})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k}},\quad
v_{j,k}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern
0.0pt}}})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k}},\quad
p_{j,k}=\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k}}\left[\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}_{j,k}-\frac{\big{(}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern
0.0pt}}})_{j,k}\big{)}^{2}+\big{(}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern
0.0pt}}})_{j,k}\big{)}^{2}}{2\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k}}-\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern
0.0pt}}}_{j,k}\right],$

and then construct the linear pieces to approximate the primitive variables
${\bm{V}}=(\rho,u,v,p,\Gamma,\Pi)^{\top}$:

$\widetilde{\bm{V}}_{j,k}(x,y)=\bm{V}_{j,k}+(\bm{V}_{x})_{j,k}(x-x_{j})+(\bm{V}_{y})_{j,k}(y-y_{k}),\quad(x,y)\in
C_{j,k},$ (3.7)

where $\bm{V}_{j,k}:=(\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}_{j,k},u_{j,k},v_{j,k},p_{j,k},\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}_{j,k},\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern
0.0pt}}}_{j,k})^{\top}$ and $(\bm{V}_{x})_{j,k}$ and $(\bm{V}_{y})_{j,k}$ are
the slopes, which are supposed to be computed using a nonlinear limiter.

As in the 1-D case, we use different limiters near and away from the material
interfaces, which need to be detected. In the two-fluid case, we check whether

$(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j.k}-\widehat{\Gamma})(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j+1,k}-\widehat{\Gamma})<0,$ (3.8)

where, as before, $\widehat{\Gamma}=(\Gamma_{\rm I}+\Gamma_{\rm II})/2$. If
(3.8) is satisfied, we then use the overcompressive SBM limiter,

$(\bm{V}_{x})_{\ell,k}=\phi^{\rm
SBM}_{\theta,\tau}\left(\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{\ell+1,k}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{\ell,k}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{\ell,k}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{\ell-1,k}}\right)\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{\ell+1,k}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern 0.0pt}}}_{\ell,k}}{\Delta x},$
(3.9)

for $\ell=j-1$, $j$, $j+1$, and $j+2$. In (3.9), the function $\phi^{\rm
SBM}_{\theta,\tau}(r)$, given by (2.15), is applied in a componentwise manner
with $\tau=-0.5$ and $\theta=1.3$. Otherwise, that is, away from the material
interface, we use a dissipative generalized minmod limiter which is also given
by the same formulae (3.9), (2.15), but with $\tau=0.5$ and $\theta=1.3$. We
proceed similarly in the $y$-direction: we use

$(\bm{V}_{y})_{j,m}=\phi^{\rm SBM}_{\theta,\tau}\left(\frac{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j,m+1}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern 0.0pt}}}_{j,m}}{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j,m}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j,m-1}}\right)\frac{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j,m+1}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern 0.0pt}}}_{j,m}}{\Delta y},$ (3.10)

with $\tau=-0.5$ and $\theta=1.3$ for $m=k-1$, $k$, $k+1$, and $k+2$ if

$(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k}-\widehat{\Gamma})(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k+1}-\widehat{\Gamma})<0,$ (3.11)

is satisfied, and with $\tau=0.5$ and $\theta=1.3$ otherwise.

Equipped with $(\bm{V}_{x})_{j,k}$ and $(\bm{V}_{y})_{j,k}$, we use (3.7) to
obtain

$\displaystyle\bm{V}^{-}_{{j+\frac{1}{2}},k}=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j,k}+\frac{\Delta
x}{2}(\bm{V}_{x})_{j,k},\quad\bm{V}^{+}_{{j+\frac{1}{2}},k}=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j+1,k}-\frac{\Delta x}{2}(\bm{V}_{x})_{j+1,k},$
$\displaystyle\bm{V}^{-}_{j,{k+\frac{1}{2}}}=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j,k}+\frac{\Delta
y}{2}(\bm{V}_{y})_{j,k},\quad\bm{V}^{+}_{j,{k+\frac{1}{2}}}=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{V}$\kern
0.0pt}}}_{j,k+1}-\frac{\Delta y}{2}(\bm{V}_{y})_{j,k+1},$

and then the corresponding point values of the conservative variables $\bm{U}$
are

${\bm{U}}^{\pm}_{\ell,m}=\left(\rho^{\pm}_{\ell,m},\rho^{\pm}_{\ell,m}u^{\pm}_{\ell,m},\rho^{\pm}_{\ell,m}v^{\pm}_{\ell,m},E^{\pm}_{\ell,m},\Gamma^{\pm}_{\ell,m},\Pi^{\pm}_{\ell,m}\right)^{\top},$

for $(\ell,m)=({j+\frac{1}{2}},k)$ and $(j,{k+\frac{1}{2}})$. Here,
$E^{\pm}_{\ell,m}=\Gamma^{\pm}_{\ell,m}p^{\pm}_{\ell,m}+\rho^{\pm}_{\ell,m}((u^{\pm}_{\ell,m})^{2}+(v^{\pm}_{\ell,m})^{2})/2+\Pi^{\pm}_{\ell,m}$.

The global fluxes $\bm{K}^{\pm}_{{j+\frac{1}{2}},k}$ and $\bm{{\cal
L}}_{j,{k+\frac{1}{2}}}$ in (3.5)–(3.6) are obtained using the relations in
(3.2), namely, by

$\bm{K}^{\pm}_{{j+\frac{1}{2}},k}=\bm{F}^{\pm}_{{j+\frac{1}{2}},k}-\bm{R}^{\pm}_{{j+\frac{1}{2}},k},\quad\bm{L}^{\pm}_{j,{k+\frac{1}{2}}}=\bm{G}^{\pm}_{j,{k+\frac{1}{2}}}-\bm{S}^{\pm}_{j,{k+\frac{1}{2}}},$

where
$\bm{F}^{\pm}_{{j+\frac{1}{2}},k}:=\bm{F}(\bm{U}^{\pm}_{{j+\frac{1}{2}},k})$,
$\bm{G}^{\pm}_{j,{k+\frac{1}{2}}}:=\bm{G}(\bm{U}^{\pm}_{j,{k+\frac{1}{2}}})$,
and the point values of the global variables $\bm{R}$ and $\bm{S}$ are
computed as follows. First, we set $\widehat{x}=x_{\frac{1}{2}}$ and
$\widehat{y}=y_{\frac{1}{2}}$ so that $\bm{R}^{-}_{{\frac{1}{2}},k}:=\bm{0}$
and $\bm{S}^{-}_{j,{\frac{1}{2}}}:=\bm{0}$. We then evaluate
$\bm{R}^{+}_{{\frac{1}{2}},k}=\bm{B}_{\bm{\Psi},{\frac{1}{2}},k}$ and
$\bm{S}^{+}_{j,{\frac{1}{2}}}=\bm{B}_{\bm{\Psi},j,{\frac{1}{2}}}$ and
recursively compute the rest of the required point values:

$\displaystyle\bm{R}^{-}_{{j+\frac{1}{2}},k}=\bm{R}^{+}_{{j-\frac{1}{2}},k}+\bm{B}^{x}_{j,k},\quad\bm{R}^{+}_{{j+\frac{1}{2}},k}=\bm{R}^{-}_{{j+\frac{1}{2}},k}+\bm{B}_{\bm{\Psi},{j+\frac{1}{2}},k},\quad
j=1,\ldots,N_{x},$
$\displaystyle\bm{S}^{-}_{j,{k+\frac{1}{2}}}=\bm{S}^{+}_{j,{k-\frac{1}{2}}}+\bm{B}^{y}_{j,k},\quad\bm{S}^{+}_{j,{k+\frac{1}{2}}}=\bm{S}^{-}_{j,{k+\frac{1}{2}}}+\bm{B}_{\bm{\Psi},j,{k+\frac{1}{2}}},\quad
k=1,\ldots,N_{y}.$

Here, $\bm{B}^{x}_{j,k}$, $\bm{B}_{\bm{\Psi},{j+\frac{1}{2}},k}$,
$\bm{B}^{y}_{j,k}$, and $\bm{B}_{\bm{\Psi},j,{k+\frac{1}{2}},k}$ are evaluated
precisely the same way as in §2.1.1:

$\displaystyle\bm{B}^{x}_{j,k}=\Big{(}0,0,0,0,\frac{\Gamma^{-}_{{j+\frac{1}{2}},k}+\Gamma^{+}_{{j-\frac{1}{2}},k}}{2}\big{(}u^{-}_{{j+\frac{1}{2}},k}-u^{+}_{{j-\frac{1}{2}},k}\big{)},\frac{\Pi^{-}_{{j+\frac{1}{2}},k}+\Pi^{+}_{{j-\frac{1}{2}},k}}{2}\big{(}u^{-}_{{j+\frac{1}{2}},k}-u^{+}_{{j-\frac{1}{2}},k}\big{)}\Big{)}^{\top},$
$\displaystyle\bm{B}_{\bm{\Psi},{j+\frac{1}{2}},k}=\Big{(}0,0,0,0,\frac{\Gamma^{+}_{{j+\frac{1}{2}},k}+\Gamma^{-}_{{j+\frac{1}{2}},k}}{2}\big{(}u^{+}_{{j+\frac{1}{2}},k}-u^{-}_{{j+\frac{1}{2}},k}\big{)},\frac{\Pi^{+}_{{j+\frac{1}{2}},k}+\Pi^{-}_{{j+\frac{1}{2}},k}}{2}\big{(}u^{+}_{{j+\frac{1}{2}},k}-u^{-}_{{j+\frac{1}{2}},k}\big{)}\Big{)}^{\top},$
$\displaystyle\bm{B}^{y}_{j,k}=\Big{(}0,0,0,0,\frac{\Gamma^{-}_{j,{k+\frac{1}{2}}}+\Gamma^{+}_{j,{k-\frac{1}{2}}}}{2}\big{(}v^{-}_{j,{k+\frac{1}{2}}}-v^{+}_{j,{k-\frac{1}{2}}}\big{)},\frac{\Pi^{-}_{j,{k+\frac{1}{2}}}+\Pi^{+}_{j,{k-\frac{1}{2}}}}{2}\big{(}v^{-}_{j,{k+\frac{1}{2}}}-v^{+}_{j,{k-\frac{1}{2}}}\big{)}\Big{)}^{\top},$
$\displaystyle\bm{B}_{\bm{\Psi},j,{k+\frac{1}{2}}}=\Big{(}0,0,0,0,\frac{\Gamma^{+}_{j,{k+\frac{1}{2}}}+\Gamma^{-}_{j,{k+\frac{1}{2}}}}{2}\big{(}v^{+}_{j,{k+\frac{1}{2}}}-v^{-}_{j,{k+\frac{1}{2}}}\big{)},\frac{\Pi^{+}_{j,{k+\frac{1}{2}}}+\Pi^{-}_{j,{k+\frac{1}{2}}}}{2}\big{(}v^{+}_{j,{k+\frac{1}{2}}}-v^{-}_{j,{k+\frac{1}{2}}}\big{)}\Big{)}^{\top}.$

Finally, the one-sided local-speeds of propagation
$a^{\pm}_{{j+\frac{1}{2}},k}$ and $b^{\pm}_{j,{k+\frac{1}{2}}}$ can be
estimated by

$\displaystyle
a^{+}_{{j+\frac{1}{2}},k}=\max\left\\{u^{-}_{{j+\frac{1}{2}},k}+c^{-}_{{j+\frac{1}{2}},k},u^{+}_{{j+\frac{1}{2}},k}+c^{+}_{{j+\frac{1}{2}},k},0\right\\},$
$\displaystyle
a^{-}_{{j+\frac{1}{2}},k}=\min\left\\{u^{-}_{{j+\frac{1}{2}},k}-c^{-}_{{j+\frac{1}{2}},k},u^{+}_{{j+\frac{1}{2}},k}-c^{+}_{{j+\frac{1}{2}},k},0\right\\},$
$\displaystyle
b^{+}_{j,{k+\frac{1}{2}}}=\max\left\\{v^{-}_{j,{k+\frac{1}{2}}}+c^{-}_{j,{k+\frac{1}{2}}},v^{+}_{j,{k+\frac{1}{2}}}+c^{+}_{j,{k+\frac{1}{2}}},0\right\\},$
$\displaystyle
b^{-}_{j,{k+\frac{1}{2}}}=\min\left\\{v^{-}_{j,{k+\frac{1}{2}}}-c^{-}_{j,{k+\frac{1}{2}}},v^{+}_{j,{k+\frac{1}{2}}}-c^{+}_{j,{k+\frac{1}{2}}},0\right\\}.$

###  3.1 “Built-in” Anti-Diffusion

In this section, we discuss the derivation of the “built-in” anti-diffusion
terms $\bm{q}_{{j+\frac{1}{2}},k}$ and $\bm{q}_{j,{k+\frac{1}{2}}}$ in
(3.5)–(3.6) in a “dimension-by-dimension” manner following the idea introduced
in [29].

In order to derive the formula for $\bm{q}_{{j+\frac{1}{2}},k}$, we consider
the 1-D restriction of the system (3.1)–(3.3) along the lines $y=y_{k}$:

$\bm{U}_{t}(x,y_{k},t)+\bm{K}\big{(}\bm{U}(x,y_{k},t)\big{)}_{x}=\bm{0},\quad
k=1,\ldots,N_{y}.$ (3.12)

We then can go through all of the steps in the derivation of the 1-D fully
discrete scheme for the systems in (3.12) following §2.2.1 up to (2.26), which
now reads as

$\widetilde{\bm{U}}^{\,\rm
int}_{{j+\frac{1}{2}},k}(x,y_{k})=\left\\{\begin{aligned} \hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{{j+\frac{1}{2}},k},\quad x<x_{j+\frac{1}{2}},\\\ \hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{{j+\frac{1}{2}},k},\quad x>x_{j+\frac{1}{2}},\end{aligned}\right.$

and the corresponding local conservation requirements (2.27) become

$a^{+}_{{j+\frac{1}{2}},k}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{{j+\frac{1}{2}},k}=\big{(}a^{+}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}\big{)}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}.$ (3.13)

In addition to the six conservation constraints given by (3.13), we have six
degrees of freedom, which we use as in the 1-D case to enforce a sharp
approximation of quasi 1-D isolated contact waves propagating in the
$x$-direction. To this end, we enforce the continuity of $u$ and $p$ across
the cell interfaces $x=x_{j+\frac{1}{2}}$ by setting

$\displaystyle\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\,\rm
int,L}_{{j+\frac{1}{2}},k}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{{j+\frac{1}{2}},k}}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\,\rm
int,R}_{{j+\frac{1}{2}},k}}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{{j+\frac{1}{2}},k}},$ $\displaystyle\frac{1}{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{{j+\frac{1}{2}},k}}\bigg{(}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,L}_{{j+\frac{1}{2}},k}-\frac{\big{(}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\,\rm
int,L}_{{j+\frac{1}{2}},k}\big{)}^{2}+\big{(}(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern 0.0pt}}})^{\,\rm
int,L}_{{j+\frac{1}{2}},k}\big{)}^{2}}{2\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{{j+\frac{1}{2}},k}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,L}_{{j+\frac{1}{2}},k}\bigg{)}$ $\displaystyle=$
$\displaystyle\,\frac{1}{\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{{j+\frac{1}{2}},k}}\bigg{(}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,R}_{{j+\frac{1}{2}},k}-\frac{\big{(}(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\,\rm
int,R}_{{j+\frac{1}{2}},k}\big{)}^{2}+\big{(}(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern 0.0pt}}})^{\,\rm
int,R}_{{j+\frac{1}{2}},k}\big{)}^{2}}{2\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{{j+\frac{1}{2}},k}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,R}_{{j+\frac{1}{2}},k}\bigg{)},$

and then proceed as in §2.2.1 to enforce sharp (yet, non-oscillatory) jumps of
the $\rho$-, $\rho v$-, $\Gamma$-, and $\Pi$-components. This leads to the
following formulae analogous to (2.30)–(2.31):

$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{{j+\frac{1}{2}},k}$ $\displaystyle=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}+\frac{\delta^{\rho}_{{j+\frac{1}{2}},k}}{a^{-}_{{j+\frac{1}{2}},k}},\quad$
$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{{j+\frac{1}{2}},k}$ $\displaystyle=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}+\frac{\delta^{\rho}_{{j+\frac{1}{2}},k}}{a^{+}_{{j+\frac{1}{2}},k}},$
$\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern 0.0pt}}})^{\,\rm
int,L}_{{j+\frac{1}{2}},k}$ $\displaystyle=(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern 0.0pt}}})^{\,\rm
int}_{{j+\frac{1}{2}},k}+\frac{\delta^{\rho
v}_{{j+\frac{1}{2}},k}}{a^{-}_{{j+\frac{1}{2}},k}},\quad$
$\displaystyle(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern 0.0pt}}})^{\,\rm
int,R}_{{j+\frac{1}{2}},k}$ $\displaystyle=(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern 0.0pt}}})^{\,\rm
int}_{{j+\frac{1}{2}},k}+\frac{\delta^{\rho
v}_{{j+\frac{1}{2}},k}}{a^{+}_{{j+\frac{1}{2}},k}},$
$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,L}_{{j+\frac{1}{2}},k}$ $\displaystyle=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}+\frac{\delta^{\Gamma}_{{j+\frac{1}{2}},k}}{a^{-}_{{j+\frac{1}{2}},k}},\quad$
$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int,R}_{{j+\frac{1}{2}},k}$ $\displaystyle=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}+\frac{\delta^{\Gamma}_{{j+\frac{1}{2}},k}}{a^{+}_{{j+\frac{1}{2}},k}},$
$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,L}_{{j+\frac{1}{2}},k}$ $\displaystyle=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}+\frac{\delta^{\Pi}_{{j+\frac{1}{2}},k}}{a^{-}_{{j+\frac{1}{2}},k}},\quad$
$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int,R}_{{j+\frac{1}{2}},k}$ $\displaystyle=\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}+\frac{\delta^{\Pi}_{{j+\frac{1}{2}},k}}{a^{+}_{{j+\frac{1}{2}},k}},$

where

$\displaystyle\delta^{\rho}_{{j+\frac{1}{2}},k}={\rm
minmod}\left(-a^{-}_{{j+\frac{1}{2}},k}\big{[}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}-\big{(}\rho^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{\ell}\big{]},\,a^{+}_{{j+\frac{1}{2}},k}\big{[}\big{(}\rho^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{r}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}\big{]}\right),$ $\displaystyle\delta^{\rho
v}_{{j+\frac{1}{2}},k}={\rm
minmod}\left(-a^{-}_{{j+\frac{1}{2}},k}\big{[}(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern 0.0pt}}})^{\,\rm
int}_{{j+\frac{1}{2}},k}-\big{(}(\rho v)^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{\ell}\big{]},\,a^{+}_{{j+\frac{1}{2}},k}\big{[}\big{(}(\rho
v)^{\rm int}_{{j+\frac{1}{2}},k}\big{)}_{r}-(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern 0.0pt}}})^{\,\rm
int}_{{j+\frac{1}{2}},k}\big{]}\right),$
$\displaystyle\delta^{\Gamma}_{{j+\frac{1}{2}},k}={\rm
minmod}\left(-a^{-}_{{j+\frac{1}{2}},k}\big{[}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}-\big{(}\Gamma^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{\ell}\big{]},\,a^{+}_{{j+\frac{1}{2}},k}\big{[}\big{(}\Gamma^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{r}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}\big{]}\right),$
$\displaystyle\delta^{\Pi}_{{j+\frac{1}{2}},k}={\rm
minmod}\left(-a^{-}_{{j+\frac{1}{2}},k}\big{[}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}-\big{(}\Pi^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{\ell}\big{]},\,a^{+}_{{j+\frac{1}{2}},k}\big{[}\big{(}\Pi^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{r}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Pi$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}\big{]}\right).$

Here, the values $\big{(}\rho^{\rm int}_{{j+\frac{1}{2}},k}\big{)}_{\ell}$,
$\big{(}\rho^{\rm int}_{{j+\frac{1}{2}},k}\big{)}_{r}$, $\big{(}(\rho v)^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{\ell}$, $\big{(}(\rho v)^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{r}$, $\big{(}\Gamma^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{\ell}$, $\big{(}\Gamma^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{r}$, $\big{(}\Pi^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{\ell}$, and $\big{(}\Pi^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{r}$, are obtained using the Taylor expansions
as it was done in (2.32).

We then proceed as in the remaining part of §2.2.1 and complete the derivation
of the fully discrete scheme (not shown here for the sake of brevity), and
after this, we pass to the semi-discrete limit and end up with the LD PCCU
flux (3.5) with the following “built-in” anti-diffusion term:

$\bm{q}_{{j+\frac{1}{2}},k}=\left(q^{\rho}_{{j+\frac{1}{2}},k},u^{*}_{{j+\frac{1}{2}},k}q^{\rho}_{{j+\frac{1}{2}},k},q^{\rho
v}_{{j+\frac{1}{2}},k},q^{E}_{{j+\frac{1}{2}},k},q^{\Gamma}_{{j+\frac{1}{2}},k},q^{\Pi}_{{j+\frac{1}{2}},k}\right)^{\top}.$
(3.14)

Here,

$\displaystyle\bm{U}^{*}_{{j+\frac{1}{2}},k}=\frac{a^{+}_{{j+\frac{1}{2}},k}\bm{U}^{+}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}\bm{U}^{-}_{{j+\frac{1}{2}},k}-\big{(}\bm{K}^{+}_{{j+\frac{1}{2}},k}-\bm{K}^{-}_{{j+\frac{1}{2}},k}\big{)}}{a^{+}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}},\quad
u^{*}_{{j+\frac{1}{2}},k}=\frac{(\rho
u)^{*}_{{j+\frac{1}{2}},k}}{\rho^{*}_{{j+\frac{1}{2}},k}},$ $\displaystyle
q^{\rho}_{{j+\frac{1}{2}},k}={\rm
minmod}\left(-a^{-}_{{j+\frac{1}{2}},k}\big{(}\rho^{*}_{{j+\frac{1}{2}},k}-\rho^{-}_{{j+\frac{1}{2}},k}\big{)},a^{+}_{{j+\frac{1}{2}},k}\big{(}\rho^{+}_{{j+\frac{1}{2}},k}-\rho^{*}_{{j+\frac{1}{2}},k}\big{)}\right),$
$\displaystyle q^{\rho v}_{{j+\frac{1}{2}},k}={\rm
minmod}\left(-a^{-}_{{j+\frac{1}{2}},k}\big{(}(\rho
v)^{*}_{{j+\frac{1}{2}},k}-(\rho
v)^{-}_{{j+\frac{1}{2}},k}\big{)},a^{+}_{{j+\frac{1}{2}},k}\big{(}(\rho
v)^{+}_{{j+\frac{1}{2}},k}-(\rho v)^{*}_{{j+\frac{1}{2}},k}\big{)}\right),$
$\displaystyle q^{\Gamma}_{{j+\frac{1}{2}},k}={\rm
minmod}\left(-a^{-}_{{j+\frac{1}{2}},k}\big{(}\Gamma^{*}_{{j+\frac{1}{2}},k}-\Gamma^{-}_{{j+\frac{1}{2}},k}\big{)},a^{+}_{{j+\frac{1}{2}},k}\big{(}\Gamma^{+}_{{j+\frac{1}{2}},k}-\Gamma^{*}_{{j+\frac{1}{2}},k}\big{)}\right),$
$\displaystyle q^{\Pi}_{{j+\frac{1}{2}},k}={\rm
minmod}\left(-a^{-}_{{j+\frac{1}{2}},k}\big{(}\Pi^{*}_{{j+\frac{1}{2}},k}-\Pi^{-}_{{j+\frac{1}{2}},k}\big{)},a^{+}_{{j+\frac{1}{2}},k}\big{(}\Pi^{+}_{{j+\frac{1}{2}},k}-\Pi^{*}_{{j+\frac{1}{2}},k}\big{)}\right),$
$\displaystyle
q^{E}_{{j+\frac{1}{2}},k}=\frac{a^{+}_{{j+\frac{1}{2}},k}a^{-}_{{j+\frac{1}{2}},k}}{a^{+}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}}\left\\{\frac{(d^{\Gamma})^{-}_{{j+\frac{1}{2}},k}\Bigg{(}(\rho
v)^{*}_{{j+\frac{1}{2}},k}+\dfrac{q^{\rho
v}_{{j+\frac{1}{2}},k}}{a^{+}_{{j+\frac{1}{2}},k}}\Bigg{)}^{2}}{2\Bigg{(}\rho^{*}_{{j+\frac{1}{2}},k}+\dfrac{q^{\rho}_{{j+\frac{1}{2}},k}}{a^{+}_{{j+\frac{1}{2}},k}}\Bigg{)}}-\frac{(d^{\Gamma})^{+}_{{j+\frac{1}{2}},k}\Bigg{(}(\rho
v)^{*}_{{j+\frac{1}{2}},k}+\dfrac{q^{\rho
v}_{{j+\frac{1}{2}},k}}{a^{-}_{{j+\frac{1}{2}},k}}\Bigg{)}^{2}}{2\Bigg{(}\rho^{*}_{{j+\frac{1}{2}},k}+\dfrac{q^{\rho}_{{j+\frac{1}{2}},k}}{a^{-}_{{j+\frac{1}{2}},k}}\Bigg{)}}\right\\}$
$\displaystyle\hskip
29.87547pt+\frac{\big{(}u^{*}_{{j+\frac{1}{2}},k}\big{)}^{2}}{2}q^{\rho}_{{j+\frac{1}{2}},k}+\frac{q^{\Gamma}_{{j+\frac{1}{2}},k}}{\Gamma^{*}_{{j+\frac{1}{2}},k}}\bigg{[}E^{*}_{{j+\frac{1}{2}},k}-\frac{\big{(}(\rho
u)^{*}_{{j+\frac{1}{2}},k}\big{)}^{2}}{2\rho^{*}_{{j+\frac{1}{2}},k}}-\Pi^{*}_{{j+\frac{1}{2}},k}\bigg{]}+q^{\Pi}_{{j+\frac{1}{2}},k},$

where

$(d^{\Gamma})^{\pm}_{{j+\frac{1}{2}},k}=1-\frac{q^{\Gamma}_{{j+\frac{1}{2}},k}}{a^{\pm}_{{j+\frac{1}{2}},k}\Gamma^{*}_{{j+\frac{1}{2}},k}}.$

Similarly, the “built-in” anti-diffusion term in the LD PCCU flux (3.6) is

$\bm{q}_{j,{k+\frac{1}{2}}}=\left(q^{\rho}_{j,{k+\frac{1}{2}}},q^{\rho
u}_{j,{k+\frac{1}{2}}},v^{*}_{j,{k+\frac{1}{2}}}q^{\rho}_{j,{k+\frac{1}{2}}},q^{E}_{j,{k+\frac{1}{2}}},q^{\Gamma}_{j,{k+\frac{1}{2}}},q^{\Pi}_{j,{k+\frac{1}{2}}}\right)^{\top}.$
(3.15)

with

$\displaystyle\bm{U}^{*}_{j,{k+\frac{1}{2}}}=\frac{b^{+}_{j,{k+\frac{1}{2}}}\bm{U}^{+}_{j,{k+\frac{1}{2}}}-b^{-}_{j,{k+\frac{1}{2}}}\bm{U}^{-}_{j,{k+\frac{1}{2}}}-\big{(}\bm{L}^{+}_{j,{k+\frac{1}{2}}}-\bm{L}^{-}_{j,{k+\frac{1}{2}}}\big{)}}{b^{+}_{j,{k+\frac{1}{2}}}-b^{-}_{j,{k+\frac{1}{2}}}},\quad
v^{*}_{j,{k+\frac{1}{2}}}=\frac{(\rho
v)^{*}_{j,{k+\frac{1}{2}}}}{\rho^{*}_{j,{k+\frac{1}{2}}}},$ $\displaystyle
q^{\rho}_{j,{k+\frac{1}{2}}}={\rm
minmod}\left(-b^{-}_{j,{k+\frac{1}{2}}}\big{(}\rho^{*}_{j,{k+\frac{1}{2}}}-\rho^{-}_{j,{k+\frac{1}{2}}}\big{)},b^{+}_{j,{k+\frac{1}{2}}}\big{(}\rho^{+}_{j,{k+\frac{1}{2}}}-\rho^{*}_{j,{k+\frac{1}{2}}}\big{)}\right),$
$\displaystyle q^{\rho u}_{j,{k+\frac{1}{2}}}={\rm
minmod}\left(-b^{-}_{j,{k+\frac{1}{2}}}\big{(}(\rho
u)^{*}_{j,{k+\frac{1}{2}}}-(\rho
u)^{-}_{j,{k+\frac{1}{2}}}\big{)},b^{+}_{j,{k+\frac{1}{2}}}\big{(}(\rho
u)^{+}_{j,{k+\frac{1}{2}}}-(\rho u)^{*}_{j,{k+\frac{1}{2}}}\big{)}\right),$
$\displaystyle q^{\Gamma}_{j,{k+\frac{1}{2}}}={\rm
minmod}\left(-b^{-}_{j,{k+\frac{1}{2}}}\big{(}\Gamma^{*}_{j,{k+\frac{1}{2}}}-\Gamma^{-}_{j,{k+\frac{1}{2}}}\big{)},b^{+}_{j,{k+\frac{1}{2}}}\big{(}\Gamma^{+}_{j,{k+\frac{1}{2}}}-\Gamma^{*}_{j,{k+\frac{1}{2}}}\big{)}\right),$
$\displaystyle q^{\Pi}_{j,{k+\frac{1}{2}}}={\rm
minmod}\left(-b^{-}_{j,{k+\frac{1}{2}}}\big{(}\Pi^{*}_{j,{k+\frac{1}{2}}}-\Pi^{-}_{j,{k+\frac{1}{2}}}\big{)},b^{+}_{j,{k+\frac{1}{2}}}\big{(}\Pi^{+}_{j,{k+\frac{1}{2}}}-\Pi^{*}_{j,{k+\frac{1}{2}}}\big{)}\right),$
$\displaystyle
q^{E}_{j,{k+\frac{1}{2}}}=\frac{b^{+}_{j,{k+\frac{1}{2}}}b^{-}_{j,{k+\frac{1}{2}}}}{b^{+}_{j,{k+\frac{1}{2}}}-b^{-}_{j,{k+\frac{1}{2}}}}\left\\{\frac{(d^{\Gamma})^{-}_{j,{k+\frac{1}{2}}}\Bigg{(}(\rho
u)^{*}_{j,{k+\frac{1}{2}}}+\dfrac{q^{\rho
u}_{j,{k+\frac{1}{2}}}}{b^{+}_{j,{k+\frac{1}{2}}}}\Bigg{)}^{2}}{2\Bigg{(}\rho^{*}_{j,{k+\frac{1}{2}}}+\dfrac{q^{\rho}_{j,{k+\frac{1}{2}}}}{b^{+}_{j,{k+\frac{1}{2}}}}\Bigg{)}}-\frac{(d^{\Gamma})^{+}_{j,{k+\frac{1}{2}}}\Bigg{(}(\rho
u)^{*}_{j,{k+\frac{1}{2}}}+\dfrac{q^{\rho
u}_{j,{k+\frac{1}{2}}}}{b^{-}_{j,{k+\frac{1}{2}}}}\Bigg{)}^{2}}{2\Bigg{(}\rho^{*}_{j,{k+\frac{1}{2}}}+\dfrac{q^{\rho}_{j,{k+\frac{1}{2}}}}{b^{-}_{j,{k+\frac{1}{2}}}}\Bigg{)}}\right\\}$
$\displaystyle\hskip
29.87547pt+\frac{\big{(}v^{*}_{j,{k+\frac{1}{2}}}\big{)}^{2}}{2}q^{\rho}_{j,{k+\frac{1}{2}}}+\frac{q^{\Gamma}_{j,{k+\frac{1}{2}}}}{\Gamma^{*}_{j,{k+\frac{1}{2}}}}\bigg{[}E^{*}_{j,{k+\frac{1}{2}}}-\frac{\big{(}(\rho
v)^{*}_{j,{k+\frac{1}{2}}}\big{)}^{2}}{2\rho^{*}_{j,{k+\frac{1}{2}}}}-\Pi^{*}_{j,{k+\frac{1}{2}}}\bigg{]}+q^{\Pi}_{j,{k+\frac{1}{2}}},$

where

$(d^{\Gamma})^{\pm}_{j,{k+\frac{1}{2}}}=1-\frac{q^{\Gamma}_{j,{k+\frac{1}{2}}}}{b^{\pm}_{j,{k+\frac{1}{2}}}\Gamma^{*}_{j,{k+\frac{1}{2}}}}.$

###### Remark 3.1

As in [29], the computation of numerical fluxes in (3.5) should be
desingularized to avoid division by zero or very small numbers. First, if
$a^{+}_{{j+\frac{1}{2}},k}<\varepsilon_{0}$ and
$a^{-}_{{j+\frac{1}{2}},k}>-\varepsilon_{0}$ for a small positive
$\varepsilon_{0}$, we replace the flux $\bm{{\cal K}}_{{j+\frac{1}{2}},k}$
with

$\bm{{\cal
K}}_{{j+\frac{1}{2}},k}=\frac{\bm{K}\big{(}\bm{U}^{-}_{{j+\frac{1}{2}},k}\big{)}+\bm{K}\big{(}\bm{U}^{+}_{{j+\frac{1}{2}},k}\big{)}}{2}.$

Similarly, if $b^{+}_{j,{k+\frac{1}{2}}}<\varepsilon_{0}$ and
$b^{-}_{j,{k+\frac{1}{2}}}>-\varepsilon_{0}$, we replace the flux $\bm{{\cal
L}}_{j,{k+\frac{1}{2}}}$ with

$\bm{{\cal
L}}_{j,{k+\frac{1}{2}}}=\frac{\bm{L}\big{(}\bm{U}^{-}_{j,{k+\frac{1}{2}}}\big{)}+\bm{L}\big{(}\bm{U}^{+}_{j,{k+\frac{1}{2}}}\big{)}}{2}.$

In addition, the computation of the energy numerical fluxes have to be
desingularized even in the case when only one of the local speeds is very
small. In particular,

$\displaystyle\mbox{if}~{}a^{+}_{{j+\frac{1}{2}},k}<\varepsilon_{0}~{}{\rm
but}~{}a^{-}_{{j+\frac{1}{2}},k}<-\varepsilon_{0},$ $\displaystyle\mbox{we
take}~{}~{}{\cal
K}^{(4)}_{{j+\frac{1}{2}},k}=u^{-}_{{j+\frac{1}{2}},k}\big{(}E^{-}_{{j+\frac{1}{2}},k}+p^{-}_{{j+\frac{1}{2}},k}\big{)};$
$\displaystyle\mbox{if}~{}a^{-}_{{j+\frac{1}{2}},k}>-\varepsilon_{0}~{}{\rm
but}~{}a^{+}_{{j+\frac{1}{2}},k}>\varepsilon_{0},$ $\displaystyle\mbox{we
take}~{}~{}{\cal
K}^{(4)}_{{j+\frac{1}{2}},k}=u^{+}_{{j+\frac{1}{2}},k}\big{(}E^{+}_{{j+\frac{1}{2}},k}+p^{+}_{{j+\frac{1}{2}},k}\big{)};$
$\displaystyle\mbox{if}~{}b^{+}_{j,{k+\frac{1}{2}}}<\varepsilon_{0}~{}{\rm
but}~{}b^{-}_{j,{k+\frac{1}{2}}}<-\varepsilon_{0},$ $\displaystyle\mbox{we
take}~{}~{}{\cal
L}^{(4)}_{j,{k+\frac{1}{2}}}=v^{-}_{j,{k+\frac{1}{2}}}\big{(}E^{-}_{j,{k+\frac{1}{2}}}+p^{-}_{j,{k+\frac{1}{2}}}\big{)};$
$\displaystyle\mbox{if}~{}b^{-}_{j,{k+\frac{1}{2}}}>-\varepsilon_{0}~{}{\rm
but}~{}b^{+}_{j,{k+\frac{1}{2}}}>\varepsilon_{0},$ $\displaystyle\mbox{we
take}~{}~{}{\cal
L}^{(4)}_{j,{k+\frac{1}{2}}}=v^{+}_{j,{k+\frac{1}{2}}}\big{(}E^{+}_{j,{k+\frac{1}{2}}}+p^{+}_{j,{k+\frac{1}{2}}}\big{)}.$

As in the 1-D case, we take $\varepsilon_{0}=10^{-12}$ in all of the numerical
examples.

###  3.2 Flux Globalization Based LD Ai-WENO PCCU Scheme

In this section, we extend the 2-D second-order flux globalization based LD
PCCU schemes from §3.1 to the fifth order of accuracy within the A-WENO
framework.

The semi-discrete fifth-order LD Ai-WENO PCCU scheme for the 2-D quasi-
conservative system (3.1)–(3.2) reads as

$\frac{{\rm d}}{{\rm
d}t}\bm{U}_{j,k}=-\frac{\bm{\mathcal{H}}^{x}_{{j+\frac{1}{2}},k}-\bm{\mathcal{H}}^{x}_{{j-\frac{1}{2}},k}}{\Delta
x}-\frac{\bm{\mathcal{H}}^{y}_{j,{k+\frac{1}{2}}}-\bm{\mathcal{H}}^{y}_{j,{k-\frac{1}{2}}}}{\Delta
y},$ (3.16)

where the fifth-order numerical fluxes
$\bm{\mathcal{H}}^{x}_{{j+\frac{1}{2}},k}$ and
$\bm{\mathcal{H}}^{y}_{j,{k+\frac{1}{2}}}$ are defined by

$\displaystyle\bm{\mathcal{H}}^{x}_{{j+\frac{1}{2}},k}=\bm{{\cal
K}}_{{j+\frac{1}{2}},k}-\frac{\Delta
x}{24}(\bm{K}_{xx})_{{j+\frac{1}{2}},k}+\frac{7(\Delta
x)^{3}}{5760}(\bm{K}_{xxxx})_{{j+\frac{1}{2}},k},$ (3.17)
$\displaystyle\bm{\mathcal{H}}^{y}_{j,{k+\frac{1}{2}}}=\bm{{\cal
L}}_{j,{k+\frac{1}{2}}}-\frac{\Delta
y}{24}(\bm{L}_{yy})_{j,{k+\frac{1}{2}}}+\frac{7(\Delta
y)^{3}}{5760}(\bm{L}_{yyyy})_{j,{k+\frac{1}{2}}}.$

Here, $\bm{{\cal K}}_{{j+\frac{1}{2}},k}$ and $\bm{{\cal
L}}_{j,{k+\frac{1}{2}}}$ are the finite-volume fluxes (3.5), and
$(\bm{K}_{xx})_{{j+\frac{1}{2}},k}$, $({\bm{K}_{xxxx}})_{{j+\frac{1}{2}},k}$,
$(\bm{L}_{yy})_{j,{k+\frac{1}{2}}}$, and
$({\bm{L}_{yyyy}})_{j,{k+\frac{1}{2}}}$ are approximations of the second- and
fourth-order spatial derivatives of $\bm{K}$ at
$(x,y)=(x_{j+\frac{1}{2}},y_{k})$ and $\bm{L}$ at
$(x,y)=(x_{j},y_{k+\frac{1}{2}})$, respectively. We compute these quantities
using the finite-difference approximations analogous to (2.52):

$\displaystyle(\bm{K}_{xx})_{{j+\frac{1}{2}},k}=\frac{1}{12(\Delta
x)^{2}}\left[-\bm{{\cal K}}_{j-\frac{3}{2},k}+16\bm{{\cal
K}}_{{j-\frac{1}{2}},k}-30\bm{{\cal K}}_{{j+\frac{1}{2}},k}+16\bm{{\cal
K}}_{j+\frac{3}{2},k}-\bm{{\cal K}}_{j+\frac{5}{2},k}\right],$
$\displaystyle(\bm{K}_{xxxx})_{{j+\frac{1}{2}},k}=\frac{1}{(\Delta
x)^{4}}\left[\bm{{\cal K}}_{j-\frac{3}{2},k}-4\bm{{\cal
K}}_{{j-\frac{1}{2}},k}+6\bm{{\cal K}}_{{j+\frac{1}{2}},k}-4\bm{{\cal
K}}_{j+\frac{3}{2},k}+\bm{{\cal K}}_{j+\frac{5}{2},k}\right],$
$\displaystyle(\bm{L}_{yy})_{j,{k+\frac{1}{2}}}=\frac{1}{12(\Delta
y)^{2}}\left[-\bm{{\cal L}}_{j,k-\frac{3}{2}}+16\bm{{\cal
L}}_{j,{k-\frac{1}{2}}}-30\bm{{\cal L}}_{j,{k+\frac{1}{2}}}+16\bm{{\cal
L}}_{j,k+\frac{3}{2}}-\bm{{\cal L}}_{j,k+\frac{5}{2}}\right],$
$\displaystyle(\bm{L}_{yyyy})_{j,{k+\frac{1}{2}}}=\frac{1}{(\Delta
y)^{4}}\left[\bm{{\cal L}}_{j,k-\frac{3}{2}}-4\bm{{\cal
L}}_{j,{k-\frac{1}{2}}}+6\bm{{\cal L}}_{j,{k+\frac{1}{2}}}-4\bm{{\cal
L}}_{j,k+\frac{3}{2}}+\bm{{\cal L}}_{j,k+\frac{5}{2}}\right].$

As in the 1-D case, the resulting semi-discrete scheme (3.16)–(3.17) is fifth-
order accurate provided the point values $\bm{U}^{\pm}_{{j+\frac{1}{2}},k}$
and $\bm{U}^{\pm}_{j,{k+\frac{1}{2}}}$ are calculated using a fifth-order
interpolation. To this end, we apply the recently proposed fifth-order Ai-
WENO-Z interpolation [14, 31, 46], which can be performed in the $x$\- and
$y$-directions similarly to the 1-D case discussed in Appendix A; we omit the
details for the sake of brevity.

###### Remark 3.2

As in the 1-D case, one needs to apply the Ai-WENO-Z interpolation procedure
in the local characteristic variables to reduce the magnitude of the numerical
oscillations. In Appendix B, we provide a detailed explanation on how to apply
the LCD to the 2-D system (3.1)–(3.3).

##  4 Numerical Examples

In this section, we apply the developed schemes to several 1-D and 2-D
numerical examples and compare the performance of the second-order flux
globalization based PCCU, the second-order flux globalization based LD PCCU,
and the fifth-order flux globalization based LD Ai-WENO PCCU schemes, which
will be referred to as the PCCU, LD PCCU, and Ai-WENO schemes, respectively.

In all of the numerical examples, we have solved the ODE systems (2.2),
(2.50), (3.4), and (3.16) using the three-stage third-order strong stability
preserving (SSP) Runge-Kutta method (see, e.g.,[16, 17]) and used the CFL
number 0.45.

###  4.1 One-Dimensional Examples

#### Example 1—“Shock-Bubble” Interaction

In the first example, we consider the “shock-bubble” interaction problem,
which is a two-fluid modification of a single-fluid example from [29]. The
initial data are given by

$(\rho,u,p;\gamma,\pi_{\infty})(x,0)=\begin{cases}(13.1538,0,1;5/3,0),&|x|<0.25,\\\
(1.3333,-0.3535,1.5;1.4,0),&x>0.75,\\\
(1,0,1;1.4,0),&\mbox{otherwise,}\end{cases}$

which correspond to a left-moving shock, initially located at $x=0.75$, and a
resting “bubble” with a radius of 0.25, initially located at the origin. These
inital data are considered in the computational domain $[-1,2]$ subject to the
solid wall boundary conditions imposed at both $x=-1$ and $x=2$.

We apply the studied PCCU, LD PCCU, and Ai-WENO schemes to this initial-
boundary value problem and compute its numerical solutions until the final
time $t=3$ on a uniform mesh with $\Delta x=1/100$. The obtained densities and
velocities are presented in Figures 4.1 and 4.2 together with the reference
solution computed by the PCCU scheme on a much finer mesh with $\Delta
x=1/2000$. As one can see, the LD PCCU scheme achieves sharper resolution than
the PCCU one, and the use of a fifth-order Ai-WENO scheme further enhances the
resolution.

Figure 4.1: Example 1: Density $\rho$ (left) and zoom at $x\in[-0.55,-0.4]$
(right) .

Figure 4.2: Example 1: Velocity $u$ (left) and zoom at $x\in[-0.2,0.9]$
(right).

#### Example 2—Water-Air “Shock-Bubble” Interaction

In the second 1-D example, which is a 1-D modification of an example from [8],
we consider a gas-liquid multifluid system, where the liquid component is
modeled by the EOS (2.9) with $\pi_{\infty}\gg 1$. The initial conditions,

$(\rho,u,p;\gamma,\pi_{\infty})=\begin{cases}(0.05,0,1;1.4,0),&|x-6|<3,\\\
(1.325,-68.525,19153;4.4,6000),&x>11.4,\\\
(1,0,1;4.4,6000),&\mbox{otherwise,}\end{cases}$

correspond to the left-moving shock, initially located at $x=11.4$, and a
resting air “bubble” with a radius 3, initially located at $x=6$. The initial
conditions are prescribed in the computational domain $[0,18]$ subject to the
free boundary conditions.

We compute the numerical solution until the final time $t=0.045$ on a uniform
mesh with $\Delta x=1/10$ by the PCCU, LD PCCU, and Ai-WENO schemes and plot
the obtained density in Figure 4.3 together with the reference solution
computed by the PCCU scheme on a much finer mesh with $\Delta x=1/400$. As one
can observe, the LD PCCU solution achieves sharper resolution (especially of
the contact wave located at about $x=3$) compared with its PCCU counterpart,
but it produces certain oscillations. The Ai-WENO scheme, on the other hand,
achieves even higher resolution and the obtained fifth-order results are
oscillation-free. We believe that this is thanks to the fact that the Ai-
WENO-Z interpolation is performed using the LCD.

Figure 4.3: Example 2: Density $\rho$ (left) and zoom at $x\in[2.5,4.8]$
(right).

#### Example 3—Water-Air Model With a Very Stiff Equation of State

In the last 1-D example taken from [2, 10], we consider another gas-liquid
multifluid system with the water component modeled using even stiffer EOS than
the one used in Example 2. The initial conditions that correspond to a severe
water-air shock tube problem,

$(\rho,u,p;\gamma,\pi_{\infty})=\begin{cases}(1000,0,10^{9};4.4,6\cdot
10^{8}),&x<0.7,\\\ (50,0,10^{5};1.4,0),&x>0.7,\end{cases}$

are prescribed in the computational domain $[0,1]$ subject to the free
boundary conditions.

We compute the numerical solutions by the studied PCCU, LD PCCU, and Ai-WENO
schemes until the final time $t=0.00025$ on a uniform mesh with $\Delta
x=1/400$. The obtained densities are shown in Figure 4.4 along with the
reference solution computed by the PCCU scheme on a much finer mesh with
$\Delta x=1/6400$. One can observe that all of the studied schemes produce
non-oscillatory numerical solutions, and the Ai-WENO solution is slightly
sharper than the solutions computed by the PCCU and LD PCCU schemes.

Figure 4.4: Example 3: Density $\rho$ (left) and zoom at $[0.79,0.87]$
(right).

###  4.2 Two-Dimensional Examples

In this section, we present four 2-D numerical examples. In all of them, we
plot the Schlieren images of the magnitude of the density gradient field,
$|\nabla\rho|$. To this end, we have used the following shading function:

$\exp\bigg{(}-\frac{80|\nabla\rho|}{\max(|\nabla\rho|)}\bigg{)},$

where the numerical derivatives of the density are computed using standard
central differencing.

#### Example 4—Shock-Helium Bubble Interaction

In the first 2-D example taken from [40, 10], a shock wave in the air hits the
light resting bubble which contains helium. We take the following initial
conditions:

$(\rho,u,v,p;\gamma,p_{\infty})=\begin{cases}(4/29,0,0,1;5/3,0),&\mbox{in
region A},\\\ (1,0,0,1;1.4,0),&\mbox{in region B},\\\
(4/3,-0.3535,0,1.5;1.4,0),&\mbox{in region C},\end{cases}$

where regions A, B, and C are outlined in Figure 4.5, and the computational
domain is $[-3,1]\times[-0.5,0.5]$. We impose the solid wall boundary
conditions on the top and bottom and the free boundary conditions on the left
and right edges of the computational domain.

Figure 4.5: Initial setting for the 2-D numerical examples.

We compute the numerical solutions until the final time $t=3$ on a uniform
mesh with $\Delta x=\Delta y=1/500$. In Figures 4.6 and 4.7, we present
different stages of the shock-bubble interaction computed by the PCCU, LD
PCCU, and Ai-WENO schemes. Notice that the bubble changes its shape and
propagates to the left, but in order to focus on the details of the bubble
structure, we only zoom at $[\sigma,\sigma+1]\times[-0.5,0.5]$ square area
containing the bubble ($\sigma$ is decreasing in time from -0.5 to -1.6). As
one can observe from the numerical results, the bubble interface develops very
complex structures after the bubble is hit by the shock, especially at large
times. The obtained results are in a good agreement with the experimental
findings presented in [18] and the numerical results reported in [40, 10, 9].
At the same time, one can see that at a small time $t=0.5$, the resolution of
the bubble interface is significantly improved by the use of the LD PCCU and
especially the Ai-WENO schemes. At larger times, the interface develops
instabilities which are smeared by a more dissipative PCCU scheme. The
differences in the achieved resolution of the small solution structures become
even more pronounced at the larger times $t=2$, 2.5, and especially at $t=3$.
This clearly indicates that the LD PCCU scheme outperforms the PCCU one, and
the further improvement in the Ai-WENO results is much more obvious than in
the 1-D examples.

Figure 4.6: Example 4: Shock-helium bubble interaction by the PCCU (left
column), LD PCCU (middle column), and Ai-WENO (right column) schemes at times
$t=0.5$, 1, and 1.5.

Figure 4.7: Same as in Figure 4.6, but at larger times $t=2$, 2.5, and 3.

#### Example 5—Shock-R22 Bubble Interaction

In the second 2-D example also taken from [40, 10], a shock wave in the air
hits the heavy resting bubble which contains R22. The initial conditions are

$(\rho,u,v,p;\gamma,\pi_{\infty})=\begin{cases}(3.1538,0,0,1;1.249,0),&\mbox{in
region A},\\\ (1,0,0,1;1.4,0),&\mbox{in region B},\\\
(4/3,-0.3535,0,1.5;1.4,0),&\mbox{in region C}.\end{cases}$

The regions A, B, and C are the same as in Example 4 and they are specified in
Figure 4.7. In this example, we impose the same boundary conditions and use
the same computational domain as in Example 4.

We compute the numerical solutions until the final time $t=3$ on a uniform
mesh with $\Delta x=\Delta y=1/500$. In Figures 4.8 and 4.9, we present
different stages of the shock-bubble interaction computed by the PCCU, LD
PCCU, and Ai-WENO schemes. As one can see, the bubble changes its shape and
propagates to the left, and in order to focus on the details of the bubble
structure, we only zoom at $[\sigma,\sigma+1]\times[-0.5,0.5]$ square area
containing the bubble ($\sigma$ is decreasing in time from -0.5 to -1.15).
Compared with Example 4, the bubble moves to the left a little slower and
develops totally different structures as the R22 is heavier than the Helium.
The obtained results are in a good agreement with the numerical results
reported in [40, 10, 9]. Similar to Example 4, at a small time $t=0.5$, the
resolution of the bubble interface is significantly improved by the use of
either the LD PCCU or Ai-WENO schemes, and the improvement in this example is
even more pronounced. By the time $t=1$, the interface develop instabilities
which are smeared by a more dissipative PCCU scheme. As time progresses, the
solutions develop very complex small structures, which are better resolved by
the schemes containing smaller amount of numerical dissipation, namely, by the
LD PCCU and Ai-WENO schemes.

Figure 4.8: Example 5: Shock-R22 bubble interaction by the PCCU (left column),
LD PCCU (middle column), and Ai-WENO (right column) schemes at times $t=0.5$
and 1.

Figure 4.9: Same as in Figure 4.8, but at larger times $t=1.5$, 2, 2.5, and 3.

#### Example 6—A Cylindrical Explosion Problem

In this 2-D example, which is a modification of an example from [52], we
consider the case where a cylindrical explosive source is located between an
air-water interface and an impermeable wall. The initial conditions are given
by

$(\rho,u,v,p;\gamma,\pi_{\infty})=\begin{cases}(1.27,0,0,8290;2,0),&\mbox{$(x-5)^{2}+(y-2)^{2}<1$},\\\
(0.02,0,0,1;1.4,0),&\mbox{$y>4$},\\\
(1,0,0,1;7.15,3309),&\mbox{otherwise},\end{cases}$

the solid wall boundary conditions are imposed at the bottom, and the free
boundary conditions are prescribed on the other sides of the computational
domain $[0,10]\times[0,6]$.

Notice that the initial conditions contain three—not two—different fluids, and
therefore we need to modify the way the fluid interfaces are detected as the
criteria (3.8) and (3.11) are applicable in the two-fluid case only. We first
set $\Gamma_{\rm I}=1/(\gamma_{\rm I}-1)$, $\Gamma_{\rm II}={1}/{(\gamma_{\rm
II}-1)}$, and $\Gamma_{III}=1/(\gamma_{\rm III}-1)$, where $\gamma_{\rm I}=2$,
$\gamma_{\rm II}=1.4$, and $\gamma_{\rm III}=7.15$ are the specific heat
ratios for the three fluids. We then introduce
$\widehat{\Gamma}_{1}:=(\Gamma_{\rm I}+\Gamma_{\rm III})/2$,
$\widehat{\Gamma}_{2}:=(\Gamma_{\rm II}+\Gamma_{\rm III})/2$, and replace the
conditions (3.8) and (3.11) with

$(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k}-\widehat{\Gamma}_{1})(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j+1,k}-\widehat{\Gamma}_{1})<0\quad{\rm or}\quad(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k}-\widehat{\Gamma}_{2})(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j+1,k}-\widehat{\Gamma}_{2})<0$

and

$(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k}-\widehat{\Gamma}_{1})(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k+1}-\widehat{\Gamma}_{1})<0\quad{\rm or}\quad(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k}-\widehat{\Gamma}_{2})(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\Gamma$\kern
0.0pt}}}_{j,k+1}-\widehat{\Gamma}_{2})<0$

respectively.

We compute the numerical solutions until the final time $t=0.02$ on a uniform
mesh with $\Delta x=\Delta y=1/80$ by the studied PCCU, LD PCCU, and Ai-WENO
schemes. In Figure 4.10, we present time snapshots of the obtained results,
which qualitatively look very similar to the numerical results reported in
[52]. As one can see, the LD PCCU scheme captures both the material interfaces
and many of the developed wave structures substantially sharper than the PCCU
scheme and the use of the Ai-WENO scheme further enhances the resolution. In
particular, one can observe more pronounced small structures in the Ai-WENO
solution (especially inside the bubble) at the large time $t=0.02$; see Figure
4.11, where we zoom at the bubble area. This example, once again, clearly
indicates that the LD PCCU and Ai-WENO schemes outperform the PCCU one.

Figure 4.10: Example 6: Solutions computed by the PCCU (left column), LD PCCU
(middle column), and Ai-WENO (right column) schemes at times $t=0.008$, 0.014,
and $0.02$.

Figure 4.11: Example 6: Solutions computed by the PCCU (left), LD PCCU
(middle), and Ai-WENO (right) schemes at time $t=0.02$; zoom at the bubble
area.

#### Example 7—Water-Air Shock-Bubble Interaction

In the last 2-D example, which is taken from [8], we consider the interaction
of a shock in water with a gas bubble. The initial conditions

$(\rho,u,v,p;\gamma,\pi_{\infty})=\begin{cases}(0.0012,0,0,1;1.4,0),&\mbox{$(x-6)^{2}+(y-6)^{2}<9$},\\\
(1.325,-68.525,0,19153;4.4,6000),&\mbox{$x>11.4$},\\\
(1,0,0,1;4.4,6000),&\mbox{otherwise},\end{cases}$

correspond to a cylindrical air bubble impacted by a Mach 1.72 shock initiated
in water. In this example, the initial data are prescribed in the
computational domain $[0,12]\times[0,12]$ with the solid wall boundary
conditions imposed on the top and bottom and the free boundary conditions on
the left and right edges of the computational domain.

We compute the numerical solutions by the studied PCCU, LD PCCU, and Ai-WENO
schemes until the final time $t=0.045$ on a uniform mesh with $\Delta x=\Delta
y=3/200$. In Figures 4.12 and 4.13, we present different stages of the
interaction process. As one can see, the bubble containing air will be
compressed by the water, propagates to the left, and changes its shape until
losing its integrity and breaking up. The obtained results are in a good
qualitative agreement with the numerical results reported in [8]. As in the
previous examples, the resolution of the bubble interface is significantly
improved by the use of the LD PCCU and Ai-WENO schemes, especially for the
small times $t=0.0204$, 0.0305, and 0.0368; see Figure 4.12. At the same time,
the differences near the bubble interfaces between the LD PCCU and Ai-WENO
solutions are minor.

Figure 4.12: Example 7: Shock-air bubble interaction computed by the PCCU
(left column), LD PCCU (middle column), and Ai-WENO (right column) schemes at
times $t=0.0204$, 0.0305, and 0.0368.

Figure 4.13: Example 7: Same as in Figure 4.12, but at larger times $t=0.0405$
and 0.045.

###### Remark 4.1

In order to preserve the positivity of the pressure computed by the Ai-WENO
scheme, we have used the following algorithm. We first detect the material
interfaces by (3.8) and (3.11) and then switch to the second-order LD PCCU
scheme (with the one-sided point values computed by the SBM limiter (2.15),
(3.7), (3.9), (3.10) in the overcompressive regime with $\theta=1.3$ and
$\tau=-0.5$) in the neighboring four cells at each side of the interface. This
results in a hybrid “mixed-order” scheme, which still achieves higher
resolution compared with the second-order LD PCCU scheme.

##  5 Conclusion

In this paper, we have developed flux globalization based low-dissipation (LD)
path-conservative central-upwind (PCCU) schemes for one- (1-D) and two-
dimensional (2-D) compressible multifluids. The LD PCCU schemes are based on
the flux globalization based PCCU schemes and employ the recently proposed
LDCU fluxes to reduce the numerical dissipation present in the original PCCU
schemes. In order to further enhance the resolution of material interfaces, we
track their locations and use the overcompressive SBM limiter in their
neighborhoods, while in the rest of the computational domain, a dissipative
generalized minmod limiter is utilized. The new second-order finite-volume
method is then extended to the fifth order of accuracy via the finite-
difference A-WENO framework. We have applied the developed schemes to a number
of 1-D and 2-D examples and the obtained numerical results clearly demonstrate
that both of the LD PCCU and LD A-WENO schemes outperform the flux
globalization based PCCU scheme that employs the original central-upwind
numerical flux from [25]. At the same time, these examples show that the
fifth-order LD Ai-WENO scheme enhances the resolution achieved by the second-
order LD PCCU scheme.

### Acknowledgment

The work of A. Kurganov was supported in part by NSFC grant 12171226 and by
the fund of the Guangdong Provincial Key Laboratory of Computational Science
and Material Design (No. 2019B030301001).

##  Appendix A 1-D Fifth-Order Ai-WENO-Z Interpolation

In this appendix, we briefly describe the fifth-order Ai-WENO-Z interpolation
recently introduced in [14, 31, 46].

Assume that the point values $W_{j+\ell}$ of a certain quantity $W$ at the