# Local Characteristic Decomposition Based Central-Upwind Scheme

Alina Chertock, Shaoshuai Chu, Michael Herty, Alexander Kurganov,
and Mária Lukáčová-Medvid’ová Department of Mathematics, North Carolina State
University, Raleigh, NC 27695, USA<EMAIL_ADDRESS>of
Mathematics, Southern University of Science and Technology, Shenzhen, 518055,
China<EMAIL_ADDRESS>of Mathematics, RWTH Aachen
University, 52056 Aachen, Germany<EMAIL_ADDRESS>of
Mathematics, SUSTech 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, University of Mainz,
Germany<EMAIL_ADDRESS>
###### Abstract

We propose novel less diffusive schemes for conservative one- and two-
dimensional hyperbolic systems of nonlinear partial differential equations
(PDEs). The main challenges in the development of accurate and robust
numerical methods for the studied systems come from the complicated wave
structures, such as shocks, rarefactions and contact discontinuities, arising
even for smooth initial conditions. In order to reduce the diffusion in the
original central-upwind schemes, we use a local characteristic decomposition
procedure to develop a new class of central-upwind schemes. We apply the
developed schemes to the one- and two-dimensional Euler equations of gas
dynamics to illustrate the performance on a variety of examples. The obtained
numerical results clearly demonstrate that the proposed new schemes outperform
the original central-upwind schemes.

Key words: Local characteristic decomposition; central-upwind schemes;
hyperbolic systems of conservative laws; Euler equations of gas dynamics

AMS subject classification: 76M12, 65M08, 76N15, 35L65, 35L67.

##  1 Introduction

This paper focuses on numerical solutions of hyperbolic systems of
conservation laws, which in the two-dimensional (2-D) case, read as

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

where $x$ and $y$ are spatial variables, $t$ is the time,
$\bm{U}\in\mathbb{R}^{d}$ is a vector of unknowns,
$\bm{F}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ and
$\bm{G}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ are the $x$\- and $y$-flux functions,
respectively.

It is well-known that the solutions of (1.1) may develop complicated wave
structures including shocks, rarefactions, and contact discontinuities even
for infinitely smooth initial data, and thus developing highly accurate and
robust shock-capturing numerical methods for solving (1.1) is a challenging
task.

Since the pioneering works [7, 22, 9], a large number of various methods had
been introduced; see, e.g., the monographs and review papers [12, 23, 15, 37,
40, 2] and references therein. Here, we focus on finite-volume (FV) methods,
in which solutions are realized in terms of their cell averages and evolved in
time according to the following algorithm. First, a piecewise polynomial
interpolant is reconstructed out of the given cell averages, and then the
evolution step is performed using the integral form of (1.1). To this end, a
proper set of time-space control volumes has to be selected. Depending on this
selection, one may distinguish between two basic classes of finite-volume
methods: upwind and central schemes. In upwind schemes, the spatial part of
the control volumes coincides with the FV cells and therefore, one needs to
(approximately) solve (generalized) Riemann problems at every cell interface;
see, e.g., [9, 40, 2] and references therein. This helps the upwind schemes to
achieve very high resolution. At the same time, it might be quite complicated
and even impossible to solve the (generalized) Riemann problems for general
systems of conservation laws. Central schemes offer an attractive simple
alternative to the upwind ones. In staggered central schemes, first proposed
for one-dimensional (1-D) systems in [30] and then extended to higher order
[24, 29] and multiple number of dimensions [25, 13, 1], the control volumes
use the staggered spatial parts so that the FV cell interfaces remain inside
the control volumes. This allows to avoid even approximately solving any
Riemann problems, which makes staggered central schemes easy to implement for
a wide variety of hyperbolic systems. The major drawback of staggered central
schemes is, however, their relatively large numerical dissipation as they
basically average over the Riemann fans rather than resolving them.

In order to reduce the amount of excessive numerical dissipation present in
central schemes, a class of central-upwind (CU) schemes have been proposed in
[17, 19]. These schemes are based on nonuniform control volumes, whose spatial
size taken to be proportional to the local speeds of propagation. This allows
one to minimize the area over which the solution averaging occurs still
without (approximately) solving any (generalized) Riemann problems. The upwind
features of the CU schemes can be seen, for example, when they are applied to
simpler systems. For instance, the CU scheme from [17] reduces to the upwind
one when it is applied to a system whose Jacobian contains only positive (only
negative) eigenvalues. Another advantage of the CU schemes is related to the
fact that unlike the staggered central schemes, they admit a particularly
simple semi-discrete form. This observation is the basis of the modifications
of the CU schemes we propose in this paper.

Even though the CU schemes from [17, 19] are quite accurate, efficient and
robust tools for a wide variety of hyperbolic systems, higher resolution of
the numerical solutions can be achieved by further reducing numerical
dissipation. This can be done in a number of different ways, for example: (i)
by introducing a more accurate evolution procedure, which leads to a “built-
in” anti-diffusion term [18]; (ii) by implementing a numerical dissipation
switch to control the amount of numerical dissipation present in the CU
schemes [16]; (iii) by obtaining more accurate estimates for the one-sided
local speeds of propagation using the discrete Rankine-Hugoniot conditions
[8].

Another way to control the amount of numerical dissipation present in the CU
schemes or any other FV methods is by adjusting the nonlinear limiting
mechanism used in the piecewise polynomial reconstruction. It is well-known
that sharper reconstructions may lead to larger numerical oscillations and a
way to reduce these oscillations is to reconstruct the characteristic
variables rather than the conservative ones; see, e.g., [33]. This can be done
using the local characteristic decomposition (LCD); see, e.g., [3, 14, 28, 31,
33, 37, 42] and references therein.

In this paper, we modify the CU schemes from [17] by applying LCD to the
numerical diffusion part of the schemes. The obtained new LCD-based CU schemes
contain substantially smaller amount of numerical dissipation, which leads to
a significantly improved resolution of the computed solution compared with the
original CU schemes. As observed above, the key idea is that the new LCD-based
CU scheme reduces to the upwind scheme when applied to 1-D linear hyperbolic
systems

$\bm{U}_{t}+A\bm{U}_{x}=\bm{0},$ (1.2)

where $A$ is a constant matrix (disregarding the sign of the eigenvalues of
$A$). This feature suggests that the new CU schemes have more built-in
upwinding compared with their predecessors.

The paper is organized as follows. In §2, we briefly describe the 1-D second-
order FV CU scheme from [17]. In §3, we introduce the proposed new 1-D LCD-
based CU scheme and show that the developed 1-D scheme reduces to the upwind
scheme when applied to the linear system (1.2). In §4, we construct the 2-D
LCD-based CU scheme. Finally, in §5, we test the proposed schemes on a number
of 1-D and 2-D numerical examples for the Euler equations of gas dynamics. We
demonstrate high accuracy, efficiency, stability, and robustness of the new
LCD-based CU schemes, which clearly outperform the second-order CU scheme from
[17, 20].

##  2 1-D Central-Upwind Scheme: A Brief Overview

In this section, we consider the 1-D hyperbolic system of conservation laws

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

and describe the second-order semi-discrete CU scheme from [17]. To this end,
we assume that the computational domain is covered with the uniform cells
$C_{j}:=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ of size $\Delta x$ centered at
$x_{j}=(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}})/2$ and denote by
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}(t)$ cell averages of
$\bm{U}(\cdot,t)$ over the corresponding intervals $C_{j}$, that is,

$\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.$

We also assume that at certain time $t\geq 0$, the cell average values
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}$ are available and
from here on we suppress the time-dependence of all of the indexed quantities
for the sake of brevity.

According to [17], the computed cell averages (2.1) are evolved in time by
solving the following system of ordinary differential equations (ODEs):

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

where $\bm{\mathcal{F}}_{j+\frac{1}{2}}$ are the CU numerical fluxes given by

$\bm{\mathcal{F}}_{j+\frac{1}{2}}=\frac{{a_{j+\frac{1}{2}}^{+}}\bm{F}^{-}_{j+\frac{1}{2}}-{a_{j+\frac{1}{2}}^{-}}\bm{F}^{+}_{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)

Here,
$\bm{F}^{\pm}_{j+\frac{1}{2}}:=\bm{F}\big{(}\bm{U}^{\pm}_{j+\frac{1}{2}}\big{)}$
and $\bm{U}^{\pm}_{j+\frac{1}{2}}$ are the right/left-sided values
reconstructed out of the given set of cell averages $\\{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}\\}$. The one-sided local speeds of propagation
$a^{\pm}_{j+\frac{1}{2}}$ are estimated using the largest and the smallest
eigenvalues of the Jacobian
$A(\bm{U}):=\frac{\partial\bm{F}}{\partial\bm{U}}(\bm{U})$, which we denote by
$\lambda_{1}(A(\bm{U}))\leq\ldots\leq\lambda_{d}(A(\bm{U}))$. This can be
done, for example, by taking

$\displaystyle
a^{+}_{j+\frac{1}{2}}=\max\left\\{\lambda_{d}\big{(}A(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\lambda_{d}\big{(}A(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},0\right\\},$
(2.4) $\displaystyle
a^{-}_{j+\frac{1}{2}}=\min\left\\{\lambda_{1}\big{(}A(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\lambda_{1}\big{(}A(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},0\right\\}.$

The (formal) order of the scheme (2.2)–(2.4) is determined by the accuracy the
point values $\bm{U}^{\pm}_{j+\frac{1}{2}}$ are reconstructed with and the
order of the ODE solver used to integrate the ODE system (2.2) in time. In
this paper, we use the LCD-based second-order piecewise linear reconstruction
described in §2.1 and the three-stage third-order strong stability preserving
(SSP) Runge-Kutta solver; see, e.g., [10, 11].

###  2.1  LCD-Based Piecewise Linear Reconstruction

In order to ensure at least second-order of accuracy of the CU scheme
(2.2)–(2.4), one has to use a second-order piecewise linear reconstruction,
which will be non-oscillatory provided the numerical derivatives are computed
using a nonlinear limiter. A library of different limiters are available; see,
e.g., [26, 30, 39, 2, 12, 23, 40]. The limiters can be classified as
dissipative, compressive or overcompressive; see [26]. When applied to the
conservative variables $\bm{U}$, dissipative limiters may introduce excessive
numerical dissipation, which may results in oversmeared solution
discontinuities. On the other hand, compressive and overcompressive limiters
tend to lead to quite large oscillations near discontinuities and artificial
sharpening smooth parts of the solution (compressive limiters may lead to
kinks, that is, jumps in the first derivatives, while the overcompressive
limiters may lead to jumps in the solution itself); see, e.g., [26]. Part of
these difficulties can be overcome by reconstructing local characteristic
variables as it was demonstrated in [33] in the context of higher-order WENO
reconstructions.

In this paper, we use the minmod limiter, which is, according to [26], the
most compressive out of dissipative limiters. As this limiter typically leads
to quite large oscillations when applied to the conservative variables, we
implement it in the LCD framework. Specifically, we first introduce the matrix
$\widehat{A}_{j+\frac{1}{2}}=A(\widehat{\bm{U}}_{j+\frac{1}{2}})$, where
$A(\bm{U})=\frac{\partial\bm{F}}{\partial\bm{U}}(\bm{U})$ is the Jacobian and
$\widehat{\bm{U}}_{j+\frac{1}{2}}$ is either a simple average
$(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}+\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j+1})/2$ or another type of average of the $\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}$
and $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j+1}$ states.

As long as the system (1.1) is strictly hyperbolic, we compute the matrices
$R_{j+\frac{1}{2}}$ and $R^{-1}_{j+\frac{1}{2}}$ such that
$R^{-1}_{j+\frac{1}{2}}\widehat{A}_{j+\frac{1}{2}}R_{j+\frac{1}{2}}$ is a
diagonal matrix and introduce the local characteristic variables
$\bm{\Gamma}_{k}$ in the neighborhood of $x=x_{j+\frac{1}{2}}$:

$\bm{\Gamma}_{k}=R^{-1}_{j+\frac{1}{2}}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{k},\quad
k=j-1,\ldots,j+2.$

Equipped with the values $\bm{\Gamma}_{j-1}$, $\bm{\Gamma}_{j}$,
$\bm{\Gamma}_{j+1}$ and $\bm{\Gamma}_{j+2}$, we reconstruct $\bm{\Gamma}$ in
the cell $C_{j}$ by computing

$(\bm{\Gamma}_{x})_{j}={\rm
minmod}\left(2\,\frac{\bm{\Gamma}_{j}-\bm{\Gamma}_{j-1}}{\Delta
x},\,\frac{\bm{\Gamma}_{j+1}-\bm{\Gamma}_{j-1}}{2\Delta
x},\,2\,\frac{\bm{\Gamma}_{j+1}-\bm{\Gamma}_{j}}{\Delta x}\right),$ (2.5)

and

$(\bm{\Gamma}_{x})_{j+1}={\rm
minmod}\left(2\,\frac{\bm{\Gamma}_{j+1}-\bm{\Gamma}_{j}}{\Delta
x},\,\frac{\bm{\Gamma}_{j+2}-\bm{\Gamma}_{j}}{2\Delta
x},\,2\,\frac{\bm{\Gamma}_{j+2}-\bm{\Gamma}_{j+1}}{\Delta x}\right),$ (2.6)

where the minmod function, defined as

${\rm
minmod}(z_{1},z_{2},\ldots):=\begin{cases}\min_{j}\\{z_{j}\\}&\mbox{if}~{}z_{j}>0\quad\forall\,j,\\\
\max_{j}\\{z_{j}\\}&\mbox{if}~{}z_{j}<0\quad\forall\,j,\\\
0&\text{otherwise,}\end{cases}$ (2.7)

is applied in the component-wise manner. The slopes (2.5) and (2.6) allow to
evaluate

$\bm{\Gamma}^{-}_{j+\frac{1}{2}}=\bm{\Gamma}_{j}+\frac{\Delta
x}{2}(\bm{\Gamma}_{x})_{j}\quad\mbox{and}\quad\bm{\Gamma}^{+}_{j+\frac{1}{2}}=\bm{\Gamma}_{j+1}-\frac{\Delta
x}{2}(\bm{\Gamma}_{x})_{j+1},$

and obtain the corresponding point values of $\bm{U}$ by

$\bm{U}^{\pm}_{j+\frac{1}{2}}=R_{j+\frac{1}{2}}\bm{\Gamma}^{\pm}_{j+\frac{1}{2}}.$
(2.8)

###### Remark 2.1

In Appendix A, we provide a detailed explanation on how the average matrix
$\widehat{A}_{j+\frac{1}{2}}$ and the corresponding matrices
$R_{j+\frac{1}{2}}$ and $R^{-1}_{j+\frac{1}{2}}$ are computed in the case of
the Euler equation of gas dynamics.

##  3 1-D LCD-Based Central-Upwind Scheme

In this section, we introduce a new 1-D LCD-based CU scheme, in which the
amount of numerical dissipation is substantially reduced compared with the
original CU scheme (2.2)–(2.3). To this end, we first rewrite the numerical
flux (2.3) of the original CU scheme in the following form:

$\bm{\mathcal{F}}_{j+\frac{1}{2}}=\frac{\bm{F}_{j}+\bm{F}_{j+1}}{2}+\bm{D}_{j+\frac{1}{2}},$
(3.1)

where $\bm{F}_{j}:=\bm{F}(\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j})$ and
$\bm{D}_{j+\frac{1}{2}}$ is the numerical diffusion given by

$\displaystyle\bm{D}_{j+\frac{1}{2}}$
$\displaystyle=\frac{a^{+}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}}\left[\bm{F}^{-}_{j+\frac{1}{2}}-\frac{\bm{F}_{j}+\bm{F}_{j+1}}{2}\right]-\frac{a^{-}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}}\left[\bm{F}^{+}_{j+\frac{1}{2}}-\frac{\bm{F}_{j}+\bm{F}_{j+1}}{2}\right]$
$\displaystyle+\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),$

which, in turn, can be rewritten with the help of the matrix
$R_{j+\frac{1}{2}}$ introduced in §2.1 and using (2.8) as follows:

$\displaystyle\bm{D}_{j+\frac{1}{2}}$
$\displaystyle=R_{j+\frac{1}{2}}P_{j+\frac{1}{2}}R^{-1}_{j+\frac{1}{2}}\left[\bm{F}^{-}_{j+\frac{1}{2}}-\frac{\bm{F}_{j}+\bm{F}_{j+1}}{2}\right]+R_{j+\frac{1}{2}}M_{j+\frac{1}{2}}R^{-1}_{j+\frac{1}{2}}\left[\bm{F}^{+}_{j+\frac{1}{2}}-\frac{\bm{F}_{j}+\bm{F}_{j+1}}{2}\right]$
(3.2)
$\displaystyle+R_{j+\frac{1}{2}}Q_{j+\frac{1}{2}}\left(\bm{\Gamma}^{+}_{j+\frac{1}{2}}-\bm{\Gamma}^{-}_{j+\frac{1}{2}}\right).$

Here, $P_{j+\frac{1}{2}}$, $M_{j+\frac{1}{2}}$, and $Q_{j+\frac{1}{2}}$ are
the diagonal matrices

$\displaystyle P_{j+\frac{1}{2}}={\rm
diag}\big{(}(P_{1})_{j+\frac{1}{2}},\ldots,(P_{d})_{j+\frac{1}{2}}\big{)},\quad
M_{j+\frac{1}{2}}={\rm
diag}\big{(}(M_{1})_{j+\frac{1}{2}},\ldots,(M_{d})_{j+\frac{1}{2}}\big{)},$
(3.3) $\displaystyle Q_{j+\frac{1}{2}}={\rm
diag}\big{(}(Q_{1})_{j+\frac{1}{2}},\ldots,(Q_{d})_{j+\frac{1}{2}}\big{)}$

with

$\big{(}(P_{i})_{j+\frac{1}{2}},(M_{i})_{j+\frac{1}{2}},(Q_{i})_{j+\frac{1}{2}}\big{)}=\frac{1}{a^{+}_{{j+\frac{1}{2}}}-a^{-}_{j+\frac{1}{2}}}\big{(}a^{+}_{j+\frac{1}{2}},-a^{-}_{j+\frac{1}{2}},a^{+}_{j+\frac{1}{2}}a^{-}_{j+\frac{1}{2}}\big{)},\quad
i=1,\ldots,d.$ (3.4)

The main idea of the new the 1-D LCD-based CU scheme is to replace the
constant entries (3.4) in the diagonal matrices $P_{j+\frac{1}{2}}$,
$M_{j+\frac{1}{2}}$, and $Q_{j+\frac{1}{2}}$ with the corresponding
characteristic entries of the diagonal matrix
$R^{-1}_{j+\frac{1}{2}}\widehat{A}_{j+\frac{1}{2}}R_{j+\frac{1}{2}}$. This is
motivated by the fact that in the linear case with $\bm{F}(\bm{U})=A\bm{U}$,
the local characteristic speed might be different for each (diagonalized)
component $\bm{\Gamma}$. Hence, a diffusion that depends on the corresponding
local speed of the each component may lead to a sharper resolution of possible
discontinuities.

Proposed modifications lead to the semi-discrete scheme

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

where the numerical fluxes $\bm{\mathcal{F}}^{\rm LCD}_{j+\frac{1}{2}}$ are
the LCD modifications of (3.1):

$\bm{\mathcal{F}}^{\rm
LCD}_{j+\frac{1}{2}}=\frac{\bm{F}_{j}+\bm{F}_{j+1}}{2}+\bm{D}^{\rm
LCD}_{j+\frac{1}{2}},$ (3.6)

with the following LCD modification of the numerical diffusion term
(3.2)–(3.4):

$\displaystyle\bm{D}^{\rm LCD}_{j+\frac{1}{2}}$
$\displaystyle=R_{j+\frac{1}{2}}P^{\rm
LCD}_{j+\frac{1}{2}}R^{-1}_{j+\frac{1}{2}}\left[\bm{F}^{-}_{j+\frac{1}{2}}-\frac{\bm{F}_{j}+\bm{F}_{j+1}}{2}\right]+R_{j+\frac{1}{2}}M^{\rm
LCD}_{j+\frac{1}{2}}R^{-1}_{j+\frac{1}{2}}\left[\bm{F}^{+}_{j+\frac{1}{2}}-\frac{\bm{F}_{j}+\bm{F}_{j+1}}{2}\right]$
(3.7) $\displaystyle+R_{j+\frac{1}{2}}Q^{\rm
LCD}_{j+\frac{1}{2}}\left(\bm{\Gamma}^{+}_{j+\frac{1}{2}}-\bm{\Gamma}^{-}_{j+\frac{1}{2}}\right).$

Here,

$\displaystyle P^{\rm LCD}_{j+\frac{1}{2}}={\rm diag}\big{(}(P^{\rm
LCD}_{1})_{j+\frac{1}{2}},\ldots,(P^{\rm
LCD}_{d})_{j+\frac{1}{2}}\big{)},\quad M^{\rm LCD}_{j+\frac{1}{2}}={\rm
diag}\big{(}(M^{\rm LCD}_{1})_{j+\frac{1}{2}},\ldots,(M^{\rm
LCD}_{d})_{j+\frac{1}{2}}\big{)},$ $\displaystyle Q^{\rm
LCD}_{j+\frac{1}{2}}={\rm diag}\big{(}(Q^{\rm
LCD}_{1})_{j+\frac{1}{2}},\ldots,(Q^{\rm LCD}_{d})_{j+\frac{1}{2}}\big{)}$

with

$\displaystyle\big{(}(P^{\rm LCD}_{i})_{j+\frac{1}{2}},(M^{\rm
LCD}_{i})_{j+\frac{1}{2}},(Q^{\rm LCD}_{i})_{j+\frac{1}{2}}\big{)}$ (3.8)
$\displaystyle=\left\\{\begin{aligned}
&\frac{1}{(\lambda^{+}_{i})_{j+\frac{1}{2}}-(\lambda^{-}_{i})_{j+\frac{1}{2}}}\big{(}(\lambda^{+}_{i})_{j+\frac{1}{2}},-(\lambda^{-}_{i})_{j+\frac{1}{2}},(\lambda^{+}_{i})_{j+\frac{1}{2}}(\lambda^{-}_{i})_{j+\frac{1}{2}}\big{)}&&\mbox{if}~{}(\lambda^{+}_{i})_{j+\frac{1}{2}}-(\lambda^{-}_{i})_{j+\frac{1}{2}}>\varepsilon,\\\
&0&&\mbox{otherwise},\end{aligned}\right.$

and

$\displaystyle(\lambda^{+}_{i})_{j+\frac{1}{2}}$
$\displaystyle=\max\left\\{\lambda_{i}\big{(}A(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},\,\lambda_{i}\big{(}A(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\,0\right\\},$
$\displaystyle(\lambda^{-}_{i})_{j+\frac{1}{2}}$
$\displaystyle=\min\left\\{\lambda_{i}\big{(}A(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},\,\lambda_{i}\big{(}A(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\,0\right\\},$

for $i=1,\ldots,d$.

Finally, $\varepsilon$ in (3.8) is a very small desingularization constant,
taken $\varepsilon=10^{-10}$ in all of the numerical examples reported in §5.

It is crucial to note that the proposed 1-D LCD-based CU scheme reduces to the
second-order semi-discrete upwind scheme when applied to a linear hyperbolic
system (1.2) with constant coefficients. This is proven in the following
lemma.

###### Lemma 3.1

If

$\bm{F}(\bm{U})=A\bm{U},$ (3.9)

where $A$ is a constant matrix, then the scheme (3.5)–(3.7) becomes the
second-order semi-discrete upwind scheme.

Proof. We note that in the linear case with constant coefficients, the
corresponding matrix composed of the right eigenvectors of $A$ is also a
constant matrix, that is, $R_{j+\frac{1}{2}}\equiv R$, and matrices $P^{\rm
LCD}_{j+\frac{1}{2}}$ and $M^{\rm LCD}_{j+\frac{1}{2}}$ reduce to

$\displaystyle P^{\rm LCD}_{j+\frac{1}{2}}\equiv P={\rm diag}\left(\max\\{{\rm
sign}(\lambda_{1}),0\\},\ldots,\max\\{{\rm sign}(\lambda_{d}),0\\}\right),$
(3.10) $\displaystyle M^{\rm LCD}_{j+\frac{1}{2}}\equiv M={\rm
diag}\left(\min\\{{\rm sign}(\lambda_{1}),0\\},\ldots,\min\\{{\rm
sign}(\lambda_{d}),0\\}\right),$

while the matrix $Q^{\rm LCD}_{j+\frac{1}{2}}\equiv 0$ as
$\lambda^{+}_{i}\lambda^{-}_{i}=\max(\lambda_{i},0)\cdot\min(\lambda_{i},0)=0$.

Taking into account the above simplifications and substituting (3.9), (3.10)
into (3.7), yields the following expression for the numerical diffusion term:

$\displaystyle\bm{D}^{\rm
LCD}_{j+\frac{1}{2}}=RPR^{-1}A\left[\bm{U}^{-}_{j+\frac{1}{2}}-\frac{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}+\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j+1}}{2}\right]-RMR^{-1}A\left[\bm{U}^{+}_{j+\frac{1}{2}}-\frac{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}+\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j+1}}{2}\right]$ (3.11)
$\displaystyle=RPR^{-1}AR\left[R^{-1}\left(\bm{U}^{-}_{j+\frac{1}{2}}-\frac{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}+\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j+1}}{2}\right)\right]-RMR^{-1}AR\left[R^{-1}\left(\bm{U}^{+}_{j+\frac{1}{2}}-\frac{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}+\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j+1}}{2}\right)\right]$
$\displaystyle=R\Lambda^{+}\left(\bm{\Gamma}^{-}_{j+\frac{1}{2}}-\frac{\bm{\Gamma}_{j}+\bm{\Gamma}_{j+1}}{2}\right)+R\Lambda^{-}\left(\bm{\Gamma}^{+}_{j+\frac{1}{2}}-\frac{\bm{\Gamma}_{j}+\bm{\Gamma}_{j+1}}{2}\right),$

where

$\displaystyle\Lambda^{+}$ $\displaystyle=PR^{-1}AR=P\Lambda={\rm
diag}\left(\max\\{\lambda_{1},0\\},\ldots,\max\\{\lambda_{d},0\\}\right),$
$\displaystyle\Lambda^{-}$ $\displaystyle=MR^{-1}AR=M\Lambda={\rm
diag}\left(\min\\{\lambda_{1},0\\},\ldots,\min\\{\lambda_{d},0\\}\right).$

Finally, combining (3.9) and (3.11) together and including them into (3.6), we
obtain

$\displaystyle\bm{\mathcal{F}}^{\rm LCD}_{j+\frac{1}{2}}$
$\displaystyle={\frac{1}{2}}A\left(\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}+\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j+1}\right)+R\Lambda^{+}\left(\bm{\Gamma}^{-}_{j+\frac{1}{2}}-\frac{\bm{\Gamma}_{j}+\bm{\Gamma}_{j+1}}{2}\right)+R\Lambda^{-}\left(\bm{\Gamma}^{+}_{j+\frac{1}{2}}-\frac{\bm{\Gamma}_{j}+\bm{\Gamma}_{j+1}}{2}\right)$
(3.12) $\displaystyle={\frac{1}{2}}(A^{+}+A^{-})\left(\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}+\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j+1}\right)+A^{+}\left(\bm{U}^{-}_{j+\frac{1}{2}}-\frac{\bm{U}_{j}+\bm{U}_{j+1}}{2}\right)+A^{-}\left(\bm{U}^{+}_{j+\frac{1}{2}}-\frac{\bm{U}_{j}+\bm{U}_{j+1}}{2}\right)$
$\displaystyle=A^{+}\bm{U}^{-}_{j+\frac{1}{2}}+A^{-}\bm{U}^{+}_{j+\frac{1}{2}},$

where $A^{\pm}=R\Lambda^{\pm}R^{-1}$. This confirms that the scheme (3.5),
(3.12) is the second-order semi-discrete upwind scheme. $\hfill\blacksquare$

##  4 2-D LCD-Based Central-Upwind Scheme

In this section, we generalize the 1-D LCD-based CU scheme introduced in §3
for the 2-D hyperbolic system of conservation laws (1.1). We design the 2-D
LCD-based CU scheme in a “dimension-by-dimension” manner, so that it reads as

$\frac{{\rm d}\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}}{{\rm
d}t}=-\frac{\bm{\mathcal{F}}^{\rm
LCD}_{{j+\frac{1}{2}},k}-\bm{\mathcal{F}}^{\rm
LCD}_{{j-\frac{1}{2}},k}}{\Delta x}-\frac{\bm{\mathcal{G}}^{\rm
LCD}_{j,{k+\frac{1}{2}}}-\bm{\mathcal{G}}^{\rm
LCD}_{j,{k-\frac{1}{2}}}}{\Delta y},$

where

$\bm{\mathcal{F}}^{\rm
LCD}_{{j+\frac{1}{2}},k}=\frac{\bm{F}_{j,k}+\bm{F}_{j+1,k}}{2}+\bm{D}^{\rm
LCD}_{{j+\frac{1}{2}},k},\quad\bm{\mathcal{G}}^{\rm
LCD}_{j,{k+\frac{1}{2}}}=\frac{\bm{G}_{j,k}+\bm{G}_{j,k+1}}{2}+\bm{D}^{\rm
LCD}_{j,{k+\frac{1}{2}}}.$

Here, $\bm{F}_{j,k}:=\bm{F}(\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k})$,
$\bm{G}_{j,k}:=\bm{G}(\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k})$, and $\bm{D}^{\rm
LCD}_{{j+\frac{1}{2}},k}$ and $\bm{D}^{\rm LCD}_{j,{k+\frac{1}{2}}}$ are the
numerical diffusion terms defined by

$\displaystyle\bm{D}^{\rm LCD}_{{j+\frac{1}{2}},k}$
$\displaystyle=R_{{j+\frac{1}{2}},k}P^{\rm
LCD}_{{j+\frac{1}{2}},k}R^{-1}_{{j+\frac{1}{2}},k}\left[\bm{F}^{\rm
E}_{j,k}-\frac{\bm{F}_{j,k}+\bm{F}_{j+1,k}}{2}\right]$
$\displaystyle+R_{{j+\frac{1}{2}},k}M^{\rm
LCD}_{{j+\frac{1}{2}},k}R^{-1}_{{j+\frac{1}{2}},k}\left[\bm{F}^{\rm
W}_{j+1,k}-\frac{\bm{F}_{j,k}+\bm{F}_{j+1,k}}{2}\right]+R_{{j+\frac{1}{2}},k}Q^{\rm
LCD}_{{j+\frac{1}{2}},k}\left(\bm{\Gamma}^{\rm W}_{j+1,k}-\bm{\Gamma}^{\rm
E}_{j,k}\right),$ $\displaystyle\bm{D}^{\rm LCD}_{j,{k+\frac{1}{2}}}$
$\displaystyle=R_{j,{k+\frac{1}{2}}}P^{\rm
LCD}_{j,{k+\frac{1}{2}}}R^{-1}_{j,{k+\frac{1}{2}}}\left[\bm{G}^{\rm
N}_{j,k}-\frac{\bm{G}_{j,k}+\bm{G}_{j,k+1}}{2}\right]$ $\displaystyle-
R_{j,{k+\frac{1}{2}}}M^{\rm
LCD}_{j,{k+\frac{1}{2}}}R^{-1}_{j,{k+\frac{1}{2}}}\left[\bm{G}^{\rm
S}_{j,k+1}-\frac{\bm{G}_{j,k}+\bm{G}_{j,k+1}}{2}\right]+R_{j,{k+\frac{1}{2}}}Q^{\rm
LCD}_{j,{k+\frac{1}{2}}}\left[\bm{\Gamma}^{\rm S}_{j,k+1}-\bm{\Gamma}^{\rm
N}_{j,k}\right].$

The matrices $R_{{j+\frac{1}{2}},k}$, $R^{-1}_{{j+\frac{1}{2}},k}$ and
$R_{j,{k+\frac{1}{2}}}$, $R^{-1}_{j,{k+\frac{1}{2}}}$ are the matrices such
that
$R^{-1}_{{j+\frac{1}{2}},k}\widehat{A}_{{j+\frac{1}{2}},k}R_{{j+\frac{1}{2}},k}$
and
$R^{-1}_{j,{k+\frac{1}{2}}}\widehat{B}_{j,{k+\frac{1}{2}}}R_{j,{k+\frac{1}{2}}}$
are diagonal. Here,
$\widehat{A}_{{j+\frac{1}{2}},k}=A(\widehat{\bm{U}}_{{j+\frac{1}{2}},k})$,
$\widehat{B}_{j,{k+\frac{1}{2}}}=B(\widehat{\bm{U}}_{j,{k+\frac{1}{2}}})$ with
$A(\bm{U})=\frac{\partial\bm{F}(\bm{U})}{\partial\bm{U}}$,
$B(\bm{U})=\frac{\partial\bm{G}(\bm{U})}{\partial\bm{U}}$, and
$\widehat{\bm{U}}_{{j+\frac{1}{2}},k}$, $\widehat{\bm{U}}_{j,{k+\frac{1}{2}}}$
are either simple averages $(\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}+\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j+1,k})/2$, $(\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}+\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j,k+1})/2$ or another type of averages of $\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}$,
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j+1,k}$ and
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}$,
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k+1}$ states,
respectively. The numerical fluxes $\bm{F}^{\rm E,W}_{j,k}$ and $\bm{G}^{\rm
N,S}_{j,k}$ are defined by $\bm{F}^{\rm E,W}_{j,k}:=\bm{F}(\bm{U}^{\rm
E,W}_{j,k})$, $\bm{G}^{\rm N,S}_{j,k}:=\bm{G}(\bm{U}^{\rm N,S}_{j,k})$, and
the details of reconstructing the point values $\bm{U}^{\rm E,W}_{j,k}$ and
$\bm{U}^{\rm N,S}_{j,k}$, $\bm{\Gamma}^{\rm E,W}_{j,k}$ and $\bm{\Gamma}^{\rm
N,S}_{j,k}$ for the 2-D Euler equations of gas dynamics are provided in
Appendix B. Finally, the diagonal matrices $P^{\rm LCD}_{{j+\frac{1}{2}},k}$,
$M^{\rm LCD}_{{j+\frac{1}{2}},k}$, $Q^{\rm LCD}_{{j+\frac{1}{2}},k}$, $P^{\rm
LCD}_{j,{k+\frac{1}{2}}}$, $M^{\rm LCD}_{j,{k+\frac{1}{2}}}$, and $Q^{\rm
LCD}_{j,{k+\frac{1}{2}}}$ are defined by

$\displaystyle P^{\rm LCD}_{{j+\frac{1}{2}},k}={\rm diag}\left(\big{(}P^{\rm
LCD}_{1}\big{)}_{{j+\frac{1}{2}},k},\ldots,\big{(}P^{\rm
LCD}_{d}\big{)}_{{j+\frac{1}{2}},k}\right),$ $\displaystyle P^{\rm
LCD}_{j,{k+\frac{1}{2}}}={\rm diag}\left(\big{(}P^{\rm
LCD}_{1}\big{)}_{j,{k+\frac{1}{2}}},\ldots,\big{(}P^{\rm
LCD}_{d}\big{)}_{j,{k+\frac{1}{2}}}\right),$ $\displaystyle M^{\rm
LCD}_{{j+\frac{1}{2}},k}={\rm diag}\left(\big{(}M^{\rm
LCD}_{1}\big{)}_{{j+\frac{1}{2}},k},\ldots,\big{(}M^{\rm
LCD}_{d}\big{)}_{{j+\frac{1}{2}},k}\right),$ $\displaystyle M^{\rm
LCD}_{j,{k+\frac{1}{2}}}={\rm diag}\left(\big{(}M^{\rm
LCD}_{1}\big{)}_{j,{k+\frac{1}{2}}},\ldots,\big{(}M^{\rm
LCD}_{d}\big{)}_{j,{k+\frac{1}{2}}}\right),$ $\displaystyle Q^{\rm
LCD}_{{j+\frac{1}{2}},k}={\rm diag}\left(\big{(}Q^{\rm
LCD}_{1}\big{)}_{{j+\frac{1}{2}},k},\ldots,\big{(}Q^{\rm
LCD}_{d}\big{)}_{{j+\frac{1}{2}},k}\right),$ $\displaystyle Q^{\rm
LCD}_{j,{k+\frac{1}{2}}}={\rm diag}\left(\big{(}Q^{\rm
LCD}_{1}\big{)}_{j,{k+\frac{1}{2}}},\ldots,\big{(}Q^{\rm
LCD}_{d}\big{)}_{j,{k+\frac{1}{2}}}\right),$

where for all $i=1,\ldots,d$

$\displaystyle\left((P^{\rm LCD}_{i})_{{j+\frac{1}{2}},k},(M^{\rm
LCD}_{i})_{{j+\frac{1}{2}},k},(Q^{\rm LCD}_{i})_{{j+\frac{1}{2}},k}\right)$
$\displaystyle=\left\\{\begin{aligned}
&\frac{1}{\Delta(\lambda_{i})_{{j+\frac{1}{2}},k}}\left((\lambda^{+}_{i})_{{j+\frac{1}{2}},k},-(\lambda^{-}_{i})_{{j+\frac{1}{2}},k},(\lambda^{+}_{i})_{{j+\frac{1}{2}},k}(\lambda^{-}_{i})_{{j+\frac{1}{2}},k}\right)&&\mbox{if}~{}\Delta(\lambda_{i})_{{j+\frac{1}{2}},k}>\varepsilon,\\\
&0&&\mbox{otherwise},\end{aligned}\right.$ $\displaystyle\left((P^{\rm
LCD}_{i})_{j,{k+\frac{1}{2}}},(M^{\rm LCD}_{i})_{j,{k+\frac{1}{2}}},(Q^{\rm
LCD}_{i})_{j,{k+\frac{1}{2}}}\right)$ $\displaystyle=\left\\{\begin{aligned}
&\frac{1}{\Delta(\mu_{i})_{j,{k+\frac{1}{2}}}}\left((\mu^{+}_{i})_{j,{k+\frac{1}{2}}},-(\mu^{-}_{i})_{j,{k+\frac{1}{2}}},(\mu^{+}_{i})_{j,{k+\frac{1}{2}}}(\mu^{-}_{i})_{j,{k+\frac{1}{2}}}\right)&&\mbox{if}~{}\Delta(\mu_{i})_{j,{k+\frac{1}{2}}}>\varepsilon,\\\
&0&&\mbox{otherwise}.\end{aligned}\right.$

Here,
$\Delta(\lambda_{i})_{{j+\frac{1}{2}},k}:=(\lambda^{+}_{i})_{{j+\frac{1}{2}},k}-(\lambda^{-}_{i})_{{j+\frac{1}{2}},k}$,
$\Delta(\mu_{i})_{j,{k+\frac{1}{2}}}:=(\mu^{+}_{i})_{j,{k+\frac{1}{2}}}-(\mu^{-}_{i})_{j,{k+\frac{1}{2}}}$,
and

$\displaystyle(\lambda^{+}_{i})_{{j+\frac{1}{2}},k}=\max\left\\{\lambda_{i}\big{(}A(\bm{U}^{\rm
E}_{j,k})\big{)},\,\lambda_{i}\big{(}A(\bm{U}^{\rm
W}_{j+1,k})\big{)},\,0\right\\},$
$\displaystyle(\lambda^{-}_{i})_{{j+\frac{1}{2}},k}=\min\left\\{\lambda_{i}\big{(}A(\bm{U}^{\rm
E}_{j,k})\big{)},\,\lambda_{i}\big{(}A(\bm{U}^{\rm
W}_{j+1,k})\big{)},\,0\right\\},$
$\displaystyle(\mu^{+}_{i})_{j,{k+\frac{1}{2}}}=\max\left\\{\mu_{i}\big{(}B(\bm{U}^{\rm
N}_{j,k})\big{)},\,\mu_{i}\big{(}B(\bm{U}^{\rm
S}_{j,k+1})\big{)},\,0\right\\},$
$\displaystyle(\mu^{-}_{i})_{j,{k+\frac{1}{2}}}=\min\left\\{\mu_{i}\big{(}B(\bm{U}^{\rm
N}_{j,k})\big{)},\,\mu_{i}\big{(}B(\bm{U}^{\rm
S}_{j,k+1})\big{)},\,0\right\\},$

where $\lambda_{i}$ and $\mu_{i}$ are the eigenvalues of the Jacobians
$A(\bm{U})$ and $B(\bm{U})$: $\lambda_{1}(A)\leq\ldots\leq\lambda_{d}(A)$ and
$\mu_{1}(B)\leq\ldots\leq\mu_{d}(B)$, respectively.

##  5 Numerical Examples

In this section, we apply the proposed LCD-based CU schemes, which will be
referred to as the New CU schemes, to the 1-D and 2-D Euler equations of gas
dynamics described in Appendices A and B, respectively. We conduct several
numerical experiments and compare the performance of the New CU schemes with
that of the corresponding 1-D and 2-D second-order CU schemes from [17] and
[20], respectively, which will be referred to as the Old CU schemes.

In Examples 1–8, we take the specific heat ratio $\gamma=1.4$, while in
Example 9, we set $\gamma=5/3$. In all of the examples, we use the CFL number
0.4.

###  5.1 One-Dimensional Examples

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

In the first example, we consider the “shock-bubble” interaction problem taken
from [21]. The initial data, given by

$(\rho,u,p)(x,0)=\begin{cases}(13.1538,0,1)&\mbox{if}~{}|x|<0.25,\\\
(1.3333,-0.3535,1.5)&\mbox{if}~{}x>0.75,\\\
(1,0,1)&\mbox{otherwise},\end{cases}$

correspond to a left-moving shock, initially located at $x=0.75$, and a bubble
of radius 0.25, initially located at the origin. We implement solid wall
boundary conditions on the left and free boundary conditions on the right of
the computational domain $[-1,1]$. We compute the numerical solution until the
final time $t=3$ using both he New and Old CU schemes on the uniform mesh with
$\Delta x=1/100$. The obtained numerical results are presented in Figure 5.1
along with the reference solution computed by the Old CU scheme on a much
finer mesh with $\Delta x=1/2000$. As one can see, the results obtained by the
New CU scheme are sharper and less oscillatory compared to the Old CU scheme
results, which have lower resolution of the shock waves and contain an
oscillation at the contact discontinuity.

Figure 5.1: Example 1: Density ($\rho$) computed by the New and Old CU schemes
and zoom at $x\in[-0.65,-0.1]$ and $x\in[-0.5,-0.4]$.

#### Example 2—Shock-Entropy Wave Interaction Problem

In this example, we consider the shock-entropy problem from [38]. The initial
conditions,

$(\rho,u,p)(x,0)=\begin{cases}(1.51695,0.523346,1.805),&x<-4.5,\\\
(1+0.1\sin(20x),0,1),&x>-4.5,\end{cases}$

correspond to a forward-facing shock wave of Mach 1.1 interacting with high-
frequency density perturbations, that is, as the shock wave moves, the
perturbations spread ahead. We set free boundary condition at the both ends of
the computational domain $[-5,5]$. We apply both New and Old CU schemes and
compute the solutions until the final time $t=5$ on a uniform mesh with
$\Delta x=1/80$. The numerical results are shown in Figure 5.2 along with the
reference solution computed by the Old CU scheme on a much finer mesh with
$\Delta x=1/400$. One can observe that the New CU scheme produces
substantially more accurate results compared to those obtained by the Old CU
scheme. This can also be clearly seen on Figure 5.2 (right), where we zoom at
the area where the solution has smooth oscillatory structures.

Figure 5.2: Example 2: Density ($\rho$) computed by the New and Old CU schemes
and zoom at $x\in[0,3]$.

#### Example 3—Stationary Contact Wave, Traveling Shock and Rarefaction Wave

In the third example, we consider the strong shocks interaction problem
proposed in [43]. The initial conditions,

$(\rho,u,p)(x,0)=\begin{cases}(1,-19.59745,1000)&\mbox{if}~{}x<0.8,\\\
(1,-19.59745,0.01)&\mbox{otherwise},\\\ \end{cases}$

are prescribed in the computational domain $[-1,1]$, in which free boundary
conditions are implemented. We compute the numerical solutions until the final
time $t=0.03$ by both the New and Old CU schemes on a uniform mesh with
$\Delta x=1/100$, as well as the reference solution, which is obtained by the
Old CU scheme on a much finer mesh with $\Delta x=1/2000$. The numerical
results, plotted in Figure 5.3, show that both schemes produce non-oscillatory
numerical solutions, but the resolution of the contact wave achieved by the
New CU scheme is higher; see also Figure 5.3 (right), where we zoom at the
neighborhood of the contact wave.

Figure 5.3: Example 3: Density ($\rho$) computed by the New and Old CU schemes
and zoom at $x\in[0.75,0.95]$.

#### Example 4—Blast Wave Problem

In the final 1-D example, we solve the strong shocks interaction problem from
[43], which is considered on the interval $[0,1]$ with the solid wall boundary
conditions at both ends and subject to the following initial conditions:

$(\rho,u,p)(x,0)=\begin{cases}(1,0,1000),&x<0.1,\\\ (1,0,0.01),&0.1\leq x\leq
0.9,\\\ (1,0,100),&x>0.9.\end{cases}$

We compute the numerical solutions until the final time $t=0.038$ by both the
New and Old CU schemes on a uniform mesh with $\Delta x=1/400$ and, as before,
implement the Old CU scheme on a fine grid with $\Delta x=1/4000$ to obtain
the corresponding reference solution. The obtained results, presented in
Figure 5.4, demonstrate that the New CU scheme achieves slightly higher
resolution of the second density spike.

Figure 5.4: Example 4: Density ($\rho$) computed by the New and Old CU schemes
and zoom at $x\in[0.55,0.85]$.

###  5.2 Two-Dimensional Examples

#### Example 5—2-D Riemann Problem

In the first 2-D example, we consider Configuration 3 of the 2-D Riemann
problems taken from [20]; see also [34, 35, 44]. The initial conditions,

$(\rho(x,y,0),u(x,y,0),v(x,y,0),p(x,y,0))=\begin{cases}(1.5,0,0,1.5),&x>1,\,y>1,\\\
(0.5323,1.206,0,0.3),&x<1,\,y>1,\\\ (0.138,1.206,1.206,0.029),&x<1,\,y<1,\\\
(0.5323,0,1.206,0.3),&x>1,\,y<1,\end{cases}$

are prescribed in the computational domain $[0,1.2]\times[0,1.2]$, in which
the free boundary conditions are implemented on all of the four sides of the
domain. We compute the numerical solutions until the final time $t=1$ by both
the New and Old CU schemes on a uniform mesh with $\Delta x=\Delta y=3/2500$,
and plot the obtained results in Figure 5.5. As one can clearly see, the New
CU scheme outperforms the Old CU scheme in capturing a sideband instability of
the jet in the zones of strong along-jet velocity shear and the instability
along the jet’s neck.

Figure 5.5: Example 5: Density ($\rho$) computed by the Old (left) and New
(right) CU schemes.

#### Example 6—Explosion Problem

In this example, we consider the explosion problem studied in [8, 18, 16, 27].
This is a circularly symmetric problem with an initial circular region of
higher density and pressure with the following initial conditions,

$(\rho(x,y,0),u(x,y,0),v(x,y,0),p(x,y,0))=\begin{cases}(1,0,0,1),&x^{2}+y^{2}<0.16,\\\
(0.125,0,0,0.1),&\mbox{otherwise},\end{cases}$

prescribed in the computational domain $[0,1.5]\times[0,1.5]$. Solid wall
boundary conditions are imposed at $x=0$ and $y=0$, while free boundary
conditions are set at $x=1.5$ and $y=1.5$. It is well-known that the solution
of this initial-boundary value problem develops circular shock and contact
waves. While the shock wave is stable and a good numerical scheme should
contain a sufficient amount of numerical dissipation to capture the shock in a
stable, non-oscillatory manner, the contact wave is unstable and can only be
stabilized numerically by the numerical diffusion present in the scheme.
Therefore, this is a good benchmark to measure the amount of numerical
dissipation present in different schemes as one ideally wants to have as
little numerical dissipation as possible, but sufficient to stabilize the
shock wave.

We apply both New and Old CU schemes and compute the numerical solutions on a
uniform mesh with $\Delta x=\Delta y=3/800$ until the final time $t=3.2$. The
obtained results are presented in Figure 5.6. One can observe that compared
with the results obtained by the Old CU scheme, the contact curve captured by
the New CU scheme is substantially “curlier” and the mixing layer is slightly
wider (indicating a more severe instability), while the shock is still stable.

Figure 5.6: Example 6: Density ($\rho$) computed by the Old (left) and New
(right) CU schemes.

#### Example 7—Implosion Problem

In this example, we consider the implosion problem taken from [8, 18, 16, 27].
The initial conditions,

$(\rho(x,y,0),u(x,y,0),v(x,y,0),p(x,y,0))=\begin{cases}(0.125,0,0,0.14),&|x|+|y|<0.15,\\\
(1,0,0,1),&\mbox{otherwise},\end{cases}$

are prescribed in the computational domain $[0,0.3]\times[0,0.3]$ with solid
boundary conditions imposed at all of the four sides. We compute the numerical
solutions until the final time $t=2.5$ by both the New and Old CU schemes on a
uniform mesh with $\Delta x=\Delta y=1/2000$. The obtained results are
depicted in Figure 5.7 (top row). As one can clearly see, a jet generated by
the New CU scheme propagates much further in the direction of $y=x$ than the
jet produced by the Old CU scheme: This is attributed to a much smaller amount
of numerical dissipation present in the New CU scheme.

Figure 5.7: Example 7: Density ($\rho$) computed by the Old CU scheme using
$\Delta x=\Delta y=1/2000$ (top left), New CU schemes using $\Delta x=\Delta
y=1/2000$ (top right), Old CU scheme using $\Delta x=\Delta y=1/2500$ (bottom
left), and Old CU scheme using $\Delta x=\Delta y=3/10000$ (bottom right).

In this example, we also check the efficiency of the proposed LCD-based CU
scheme by performing a more thorough comparison between the New and Old CU
schemes: We now take into account additional computational cost of the New CU
scheme. To this end, we measure the CPU time consumed during the above
calculations by the New CU scheme and refine the mesh used by the Old CU
scheme to the level that exactly the same CPU time is consumed to compute both
of the numerical solutions. The corresponding grids are $\Delta x=\Delta
y=1/2000$ for the New CU scheme, and $\Delta x=\Delta y=1/2500$ for the Old CU
scheme (this solution is plotted in Figure 5.7 (bottom left). The obtained
numerical results indicate that the New CU scheme still achieves a much higher
resolution than the Old CU scheme.

It is instructive to check whether the solution computed by the Old CU scheme
on even finer mesh will be comparable with the New CU solution. To this end,
we refine the mesh to the level $\Delta x=\Delta y=3/10000$ and plot the
obtained Old CU scheme in Figure 5.7 (bottom right). As one can see, this
solution is similar to the New CU solution shown in Figure 5.7 (top right).
This shows that the New CU scheme is capable of resolving the same details on
a substantially coarser grid due to the smaller amount of numerical diffusion
present in the scheme.

#### Example 8—Kelvin-Helmholtz (KH) Instability

In this example, we study the KH instability, which develops in the test
problem taken from [5, 6, 8, 16, 32]. We take the following initial data:

$\displaystyle(\rho(x,y,0),u(x,y,0))=\begin{cases}(1,-0.5+0.5e^{(y+0.25)/L}),&y<-0.25,\\\
(2,0.5-0.5e^{(-y-0.25)/L}),&-0.25<y<0,\\\
(2,0.5-0.5e^{(y-0.25)/L}),&0<y<0.25,\\\
(1,-0.5+0.5e^{(0.25-y)/L}),&y>0.25,\end{cases}$ $\displaystyle
v(x,y,0)=0.01\sin(4\pi x),\qquad p(x,y,0)\equiv 1.5,$

where $L$ is a smoothing parameter (we take $L=0.00625$), which corresponds to
a thin shear interface with a perturbed vertical velocity field $v$ in the
conducted simulations. Periodic boundary conditions are imposed on all of the
four sides of the computational domain $[-0.5,0.5]\times[-0.5,0.5]$. We
compute the numerical solutions until the final time $t=4$ by both the New and
Old CU schemes on a uniform mesh with $\Delta x=\Delta y=1/1024$. The
numerical results at times $t=1$, 2.5, and 4 are presented in Figure 5.8. As
one can see, at the early time $t=1$, the vortex streets formed by the
adaptive scheme are more pronounced. These structures continue growing
exponentially in time showing much more complicated turbulent mixing captured
by the New CU scheme, especially at the final time $t=4$. Clearly, the New CU
scheme outperforms the Old one in achieving a higher resolution of the KH
instabilities.

Figure 5.8: Example 8: Time snapshots of density ($\rho$) computed by the Old
(left column) and New (right column) CU schemes at $t=1$ (top row), $t=2.5$
(middle row), and $t=4$ (bottom row).

For the KH instability problem, it is well-known that numerical solutions do
not converge when the mesh is refined, since the limiting solution is not a
weak solution but a so-called dissipative weak solution; see [4] for more
details. The latter can be approximated by the Cesàro averages computed from
the numerical solutions on meshes with different resolutions. As proved in
[4], the Cesàro averages converge strongly (in the $L^{q}$-norm with $1\leq
q<\infty$) to a dissipative weak solution. In order to approximate the
limiting solution, we consider the Cesàro averages of the densities computed
at the final time $t=4$ by the New and Old CU schemes. To this end, we
introduce the sequence of meshes with the cells of size ${1}/{2^{n}}$,
$n=5,\ldots,10$, and denote by $\rho(1/2^{n})$ the density computed on the
corresponding mesh. We then project the obtained coarser mesh solutions with
$n=5,\ldots,m-1$ onto the finer mesh with $n=m$ (the projection is carried out
using the minmod reconstruction [26, 30, 39] applied to the density field) and
denote the obtained results still by $\rho(1/2^{n})$, $n=5,\ldots,m$. After
this, we compute the Cesàro averages by

$\rho^{\rm C}(1/2^{m})=\frac{\rho(1/2^{5})+\cdots+\rho(1/2^{n})}{m-4},\quad
m=8,9,10.$

The obtained results, presented in Figure 5.9 for both the New and Old CU
schemes, indicate non-uniqueness of the limiting dissipative weak solution
that apparently depends on the choice of the numerical diffusion. This leads
to an interesting question on suitable selection criteria that can be studied
in the future.

Figure 5.9: Example 8: Cesàro averages of the density $\rho^{\rm C}$ computed
by the Old (left column) and New (right column) CU schemes at $t=4$.

#### Example 9—Rayleigh-Taylor (RT) Instability

In the last example taken from [8, 16, 36, 41], we investigate the RT
instability, which is a physical phenomenon occurring when a layer of heavier
fluid is placed on top of a layer of lighter fluid. The model is governed by
the 2-D Euler equations (1.1), (B.1)–(B.2) with added gravitational source
terms. In the studied setup, the gravitational force acts in the positive
$y$-direction and the modified system reads as

$\displaystyle\rho_{t}+(\rho u)_{x}+(\rho v)_{y}=0,$ $\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}=\rho,$ $\displaystyle
E_{t}+[u(E+p)]_{x}+[v(E+p)]_{y}=\rho v.$

This system is considered subject to the following initial conditions:

$(\rho(x,y,0),u(x,y,0),v(x,y,0),p(x,y,0))=\begin{cases}(2,0,-0.025\,c\cos(8\pi
x),2y+1),&y<0.5,\\\ (1,0,-0.025\,c\cos(8\pi
x),y+1.5),&\mbox{otherwise},\end{cases}$

where $c:=\sqrt{\gamma p/\rho}$ is the speed of sound, and the solid wall
boundary conditions at $x=0$ and $x=0.25$, and the following Dirichlet
boundary conditions at the top and bottom boundaries:

$(\rho,u,v,p)|_{y=1}=(1,0,0,2.5),\qquad(\rho,u,v,p)|_{y=0}=(2,0,0,1).$

We compute the numerical solutions in the computational domain
$[0,0.25]\times[0,1]$ discretized using a uniform mesh with $\Delta x=\Delta
y=1/1024$ until the final time $t=2.95$ by both the New and Old CU schemes.
The numerical results at times $t=1.95$ and 2.95, presented in Figure 5.10,
show a significant difference in performance of the New and Old CU schemes.
Indeed, the New CU scheme is capable of capturing the RT instability with a
much higher resolution of the complicated solution structure. Once again, this
indicates that the New CU scheme is substantially less dissipative than its
Old counterpart.

Figure 5.10: Example 9: Density ($\rho$) computed by the Old and New CU
schemes at different times.

### Acknowledgment

The work of A. Chertock was supported in part by NSF grant DMS-1818684. The
work of M. Herty was supported in part by the DFG (German Research Foundation)
through 20021702/GRK2326, 333849990/IRTG-2379, HE5386/18-1, 19-2, 22-1, 23-1
and under Germany’s Excellence Strategy EXC-2023 Internet of Production
390621612. The work of A. Kurganov was supported in part by NSFC grants
12171226 and 12111530004, and by the fund of the Guangdong Provincial Key
Laboratory of Computational Science and Material Design (No. 2019B030301001).
The work of M. Lukáčová-Medvid’ová has been funded by the DFG (German Research
Foundation) under the SFB/TRR 146 Multiscale Simulation Methods for Soft
Matter Systems. M. Lukáčová-Medvid’ová gratefully acknowledges support of the
Gutenberg Research College of University Mainz and the Mainz Institute of
Multiscale Modeling.

##  Appendix A LCD Matrices for the 1-D Euler Equations of Gas Dynamics

The 1-D Euler equations of gas dynamics read as the system (2.1) with

$\bm{U}=\begin{pmatrix}\rho\\\ \rho u\\\
E\end{pmatrix}\quad\mbox{and}\quad\bm{F}(\bm{U})=\begin{pmatrix}\rho u\\\ \rho
u^{2}+p\\\ u(E+p)\end{pmatrix},$

where $\rho$, $u$, $p$, and $E$ are the density, velocity, pressure and total
energy, respectively. The system is completed through the following equations
of state:

$E=\rho e+{\frac{1}{2}}\rho u^{2},\quad e=\frac{p}{(\gamma-1)\rho},$

where $e$ is the specific internal energy, and $\gamma$ is the parameter
representing the specific heat ratio.

We first compute the Jacobian

$\displaystyle
A(\bm{U})=\frac{\partial\bm{F}}{\partial\bm{U}}(\bm{U})=\begin{pmatrix}0&1&0\\\
\dfrac{\gamma-3}{2}u^{2}&(3-\gamma)u&\gamma-1\\\\[6.45831pt] -\dfrac{\gamma
u\hat{E}}{\rho}+(\gamma-1)u^{3}&H-(\gamma-1)u^{2}&\gamma u\end{pmatrix},$

where $H=\frac{E+p}{\rho}$ is the total specific enthalpy. We then introduce
the matrices

$\displaystyle\widehat{A}_{j+\frac{1}{2}}=\begin{pmatrix}0&1&0\\\
\dfrac{\gamma-3}{2}\hat{u}^{2}&(3-\gamma)\hat{u}&\gamma-1\\\\[6.45831pt]
-\dfrac{\gamma\hat{u}\hat{E}}{\hat{\rho}}+(\gamma-1)\hat{u}^{3}&\hat{H}-(\gamma-1)\hat{u}^{2}&\gamma\hat{u}\end{pmatrix},$

where $\hat{(\cdot)}$ stands for the following averages:

$\hat{\rho}=\frac{\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$\rho$\kern
0.0pt}}}_{j+1}}{2},~{}\hat{u}=\frac{u_{j}+u_{j+1}}{2},~{}\hat{p}=\frac{p_{j}+p_{j+1}}{2},~{}\hat{H}=\frac{\hat{E}+\hat{p}}{\hat{\rho}},~{}\hat{E}=\frac{\hat{p}}{\gamma-1}+{\frac{1}{2}}\hat{\rho}\hat{u}^{2},$

where

$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\mbox{and}\quad p_{j}=(\gamma-1)\left[\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}_{j}-{\frac{1}{2}}\rho_{j}u_{j}^{2}\right].$

Notice that all of the $\hat{(\cdot)}$ quantities have to have a subscript
index, that is, $\hat{(\cdot)}=\hat{(\cdot)}_{j+\frac{1}{2}}$, but we omit it
for the sake of brevity for all of the quantities except for $\widehat{A}$. We
then compute the matrix $R_{j+\frac{1}{2}}$ composed of the right eigenvectors
of $\widehat{A}_{j+\frac{1}{2}}$ and obtain

$R_{j+\frac{1}{2}}=\begin{pmatrix}1&1&1\\\
\hat{u}-\hat{c}&\hat{u}&\hat{u}+\hat{c}\\\
\hat{H}-\hat{u}\hat{c}&\dfrac{\hat{u}^{2}}{2}&\hat{H}+\hat{u}\hat{c}\end{pmatrix}\quad\mbox{and}\quad
R^{-1}_{j+\frac{1}{2}}=\dfrac{1}{\hat{\phi}}\begin{pmatrix}\dfrac{\hat{u}^{2}}{2}+\dfrac{\hat{u}\hat{\phi}}{2\hat{c}}&-\hat{u}-\dfrac{\hat{\phi}}{2\hat{c}}&1\\\\[5.16663pt]
2\hat{\phi}-{2\hat{H}}&{2\hat{u}}&-2\\\
\dfrac{\hat{u}^{2}}{2}-\dfrac{\hat{u}\hat{\phi}}{2\hat{c}}&-\hat{u}+\dfrac{\hat{\phi}}{2\hat{c}}&1\end{pmatrix},$

where $\hat{\phi}=2\hat{H}-\hat{u}^{2}$ and
$\hat{c}=\sqrt{\gamma\hat{p}/\hat{\rho}}$.

##  Appendix B LCD-Based Piecewise Linear Reconstruction for the 2-D Euler
Equations of Gas Dynamics

The 2-D Euler equations of gas dynamics read as (1.1) with

$\bm{U}=\begin{pmatrix}\rho\\\ \rho u\\\ \rho v\\\
E\end{pmatrix},\quad\bm{F}(\bm{U})=\begin{pmatrix}\rho u\\\ \rho u^{2}+p\\\
\rho uv\\\
u(E+p)\end{pmatrix},\quad\mbox{and}\quad\bm{G}(\bm{U})=\begin{pmatrix}\rho
v\\\ \rho uv\\\ \rho v^{2}+p\\\ v(E+p)\end{pmatrix},$ (B.1)

where $\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 is completed through
the following equations of state:

$E=\rho e+\frac{\rho}{2}(u^{2}+v^{2}),\quad e=\frac{p}{(\gamma-1)\rho}.$ (B.2)

where, as in the 1-D case, $e$ is the specific internal energy, and $\gamma$
is the parameter representing the specific heat ratio.

In what follows, we discuss the reconstruction of the point values
$\bm{U}^{\rm E,W}_{j,k}$ and $\bm{\Gamma}^{\rm E,W}_{j,k}$ (the point values
$\bm{U}^{\rm N,S}_{j,k}$ and $\bm{\Gamma}^{\rm N,S}_{j,k}$ can be computed in
a similar manner and we omit the details for the sake of brevity). To this
end, we first compute the Jacobian

$\displaystyle
A(\bm{U})=\frac{\partial\bm{F}}{\partial\bm{U}}(\bm{U})=\begin{pmatrix}0&1&0&0\\\
\dfrac{\gamma-3}{2}u^{2}+\dfrac{\gamma-1}{2}v^{2}&(3-\gamma)u&(1-\gamma)v&\gamma-1\\\
-uv&v&u&0\\\ -\dfrac{\gamma
uE}{\rho}+(\gamma-1)u(u^{2}+v^{2})&H-(\gamma-1)u^{2}&(1-\gamma)uv&\gamma
u\end{pmatrix},$

where, as in the 1-D case, $H=\frac{E+p}{\rho}$. We then introduce the
matrices

$\displaystyle\widehat{A}_{{j+\frac{1}{2}},k}=\begin{pmatrix}0&1&0&0\\\
\dfrac{\gamma-3}{2}\hat{u}^{2}+\dfrac{\gamma-1}{2}\hat{v}^{2}&(3-\gamma)\hat{u}&(1-\gamma)\hat{v}&\gamma-1\\\
-\hat{u}\hat{v}&\hat{v}&\hat{u}&0\\\
-\dfrac{\gamma\hat{u}\hat{E}}{\hat{\rho}}+(\gamma-1)\hat{u}(\hat{u}^{2}+\hat{v}^{2})&\hat{H}-(\gamma-1)\hat{u}^{2}&(1-\gamma)\hat{u}\hat{v}&\gamma\hat{u}\end{pmatrix},$

where $\hat{(\cdot)}$ stands for the following averages:

$\displaystyle\hat{\rho}=\frac{\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$\rho$\kern
0.0pt}}}_{j+1,k}}{2},\quad\hat{u}=\frac{u_{j,k}+u_{j+1,k}}{2},\quad\hat{v}=\frac{v_{j,k}+v_{j+1,k}}{2},\quad\hat{p}=\frac{p_{j,k}+p_{j+1,k}}{2},$
$\displaystyle\hat{H}=\frac{\hat{E}+\hat{p}}{\hat{\rho}},\quad\hat{E}=\frac{\hat{p}}{\gamma-1}+\frac{\hat{\rho}}{2}(\hat{u}^{2}+\hat{v}^{2}),$

where

$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\mbox{and}\quad
p_{j,k}=(\gamma-1)\left[\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}_{j,k}-\frac{\rho_{j,k}}{2}\left(u_{j,k}^{2}+v_{j,k}^{2}\right)\right].$

As in Appendix A, all of the $\hat{(\cdot)}$ quantities have to have a
subscript index, that is, $\hat{(\cdot)}=\hat{(\cdot)}_{{j+\frac{1}{2}},k}$,
which we omit for the sake of brevity for all of the quantities except for
$\widehat{A}$. We then compute the matrix $R_{{j+\frac{1}{2}},k}$ composed of
the right eigenvectors of $\widehat{A}_{{j+\frac{1}{2}},k}$ and obtain

$R_{{j+\frac{1}{2}},k}=\begin{pmatrix}1&1&0&1\\\
\hat{u}-\hat{c}&\hat{u}&0&\hat{u}+\hat{c}\\\ \hat{v}&\hat{v}&1&\hat{v}\\\
\hat{H}-\hat{u}\hat{c}&\dfrac{\hat{u}^{2}+\hat{v}^{2}}{2}&\hat{v}&\hat{H}+\hat{u}\hat{c}\end{pmatrix}$

and

$R^{-1}_{{j+\frac{1}{2}},k}=\dfrac{1}{\hat{\phi}}\begin{pmatrix}\dfrac{\hat{u}^{2}+\hat{v}^{2}}{2}+\dfrac{\hat{u}\hat{\phi}}{2\hat{c}}&-\hat{u}-\dfrac{\hat{\phi}}{2\hat{c}}&-\hat{v}&1\\\\[5.16663pt]
2\hat{\phi}-2\hat{H}&2\hat{u}&2\hat{v}&-2\\\
-\hat{v}\hat{\phi}&0&\hat{\phi}&0\\\
\dfrac{\hat{u}^{2}+\hat{v}^{2}}{2}-\dfrac{\hat{u}\hat{\phi}}{2\hat{c}}&-\hat{u}+\dfrac{\hat{\phi}}{2\hat{c}}&-\hat{v}&1\end{pmatrix},$

where $\hat{\phi}=2\hat{H}-\hat{u}^{2}-\hat{v}^{2}$ and
$\hat{c}=\sqrt{\gamma\hat{p}/\hat{\rho}}$.

We then introduce the local characteristic variables in the neighborhood of
$(x,y)=(x_{j+\frac{1}{2}},y_{k})$:

$\bm{\Gamma}_{\ell,k}=R^{-1}_{{j+\frac{1}{2}},k}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{\ell,k},\quad\ell=j-1,\ldots,j+2.$

Equipped with the values $\bm{\Gamma}_{j-1,k}$, $\bm{\Gamma}_{j,k}$,
$\bm{\Gamma}_{j+1,k}$, and $\bm{\Gamma}_{j+2,k}$, we compute

$(\bm{\Gamma}_{x})_{j,k}={\rm
minmod}\left(2\,\frac{\bm{\Gamma}_{j,k}-\bm{\Gamma}_{j-1,k}}{\Delta
x},\,\frac{\bm{\Gamma}_{j+1,k}-\bm{\Gamma}_{j-1,k}}{2\Delta
x},\,2\,\frac{\bm{\Gamma}_{j+1,k}-\bm{\Gamma}_{j,k}}{\Delta x}\right),$

and

$(\bm{\Gamma}_{x})_{j+1,\,k}={\rm
minmod}\left(2\,\frac{\bm{\Gamma}_{j+1,k}-\bm{\Gamma}_{j,k}}{\Delta
x},\,\frac{\bm{\Gamma}_{j+2,k}-\bm{\Gamma}_{j,k}}{2\Delta
x},\,2\,\frac{\bm{\Gamma}_{j+2,k}-\bm{\Gamma}_{j+1,k}}{\Delta x}\right),$

where the minmod function, defined in (2.7), is applied in the component-wise
manner. We then use these slopes to evaluate

$\bm{\Gamma}^{\rm E}_{j,k}=\,\bm{\Gamma}_{j,k}+\frac{\Delta
x}{2}(\bm{\Gamma}_{x})_{j,k}\quad\mbox{and}\quad\bm{\Gamma}^{\rm
W}_{j+1,k}=\bm{\Gamma}_{j+1,k}-\frac{\Delta x}{2}(\bm{\Gamma}_{x})_{j+1,k},$

and finally obtain the corresponding point values of $\bm{U}$ by

$\bm{U}^{\rm E,W}_{j+\frac{1}{2}}=R_{{j+\frac{1}{2}},k}\bm{\Gamma}^{\rm
E,W}_{{j+\frac{1}{2}},k}.$

## References

  * [1] P. Arminjon and M. C. Viallon, Généralisation du schéma de Nessyahu-Tadmor pour une équation hyperbolique à deux dimensions d’espace, C. R. Acad. Sci. Paris Sér. I Math., 320 (1995), pp. 85–88.
  * [2] M. Ben-Artzi and J. Falcovitz, Generalized Riemann problems in computational fluid dynamics, vol. 11 of Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, 2003.
  * [3] W. S. Don, D.-M. Li, Z. Gao, and B.-S. Wang, A characteristic-wise alternative WENO-Z finite difference scheme for solving the compressible multicomponent non-reactive flows in the overestimated quasi-conservative form, J. Sci. Comput., 82 (2020).  Paper No. 27, 24 pp.
  * [4] E. Feireisl, M. Lukáčová-Medviďová, H. Mizerová, and B. She, Numerical Analysis of Compressible Fluid Flows, vol. 20 of MS&A. Modeling, Simulation and Applications, Springer, Cham, 2021.
  * [5] E. Feireisl, M. Lukáčová-Medviďová, B. She, and Y. Wang, Computing oscillatory solutions of the Euler system via $\mathcal{K}$-convergence, Math. Models Methods Appl. Sci., 31 (2021), pp. 537–576.
  * [6] U. Fjordholm, S. Mishra, and E. Tadmor, On the computation of measure-valued solutions, Acta Numer., 25 (2016), pp. 567–679.
  * [7] K. O. Friedrichs, Symmetric hyperbolic linear differential equations, Comm. Pure Appl. Math., 7 (1954), pp. 345–392.
  * [8] N. K. Garg, A. Kurganov, and Y. Liu, Semi-discrete central-upwind Rankine-Hugoniot schemes for hyperbolic systems of conservation laws, J. Comput. Phys., 428 (2021), p. 110078.
  * [9] S. K. Godunov, A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics, Mat. Sb. (N.S.), 47 (1959), pp. 271–306.
  * [10] S. Gottlieb, D. Ketcheson, and C.-W. Shu, Strong stability preserving Runge-Kutta and multistep time discretizations, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2011.
  * [11] S. Gottlieb, C.-W. Shu, and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev., 43 (2001), pp. 89–112.
  * [12] J. S. Hesthaven, Numerical methods for conservation laws: From analysis to algorithms, Comput. Sci. Eng. 18, SIAM, Philadelphia, 2018.
  * [13] G. S. Jiang and E. Tadmor, Nonoscillatory central schemes for multidimensional hyperbolic conservation laws, SIAM J. Sci. Comput., 19 (1998), pp. 1892–1917.
  * [14] Y. Jiang, C.-W. Shu, and M. Zhang, An alternative formulation of finite difference weighted ENO schemes with Lax-Wendroff time discretization for conservation laws, SIAM J. Sci. Comput., 35 (2013), pp. A1137–A1160.
  * [15] D. I. Ketcheson, R. J. LeVeque, and M. J. del Razo, Riemann problems and Jupyter solutions, vol. 16 of Fundamentals of Algorithms, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2020.
  * [16] A. Kuganov, Y. Liu, and Y. Zeitlin, Numerical dissipation switch for two-dimensional central-upwind schemes, ESAIM: Math. Model. Num. Anal., 55 (2021), pp. 713–734.
  * [17] A. Kuganov, P. Noelle, and G. Petrova, Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations, SIAM J. Sci. Comput., 23 (2001), pp. 713–734.
  * [18] A. Kurganov and C.-T. Lin, On the reduction of numerical dissipation in central-upwind schemes, Commun. Comput. Phys., 2 (2007), pp. 141–163.
  * [19] A. Kurganov and E. Tadmor, New high-resolution semi-discrete central schemes for Hamilton-Jacobi equations, J. Comput. Phys., 160 (2000), pp. 720–742.
  * [20] , Solution of two-dimensional riemann problems for gas dynamics without Riemann problem solvers, Numer. Methods Partial Differential Equations, 18 (2002), pp. 584–608.
  * [21] A. Kurganov and R. Xin, Central-upwind schemes and contact discontinuities. Application to the Euler equations of gas dynamics.  In preparation.
  * [22] P. D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical computation, Comm. Pure Appl. Math., 7 (1954), pp. 159–193.
  * [23] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Appl. Math., Cambridge University Press, Cambridge, UK, 2002\.
  * [24] D. Levy, G. Puppo, and G. Russo, Central WENO schemes for hyperbolic systems of conservation laws, M2AN Math. Model. Numer. Anal., 33 (1999), pp. 547–571.
  * [25] K.-A. Lie and S. Noelle, An improved quadrature rule for the flux-computation in staggered central difference schemes in multidimensions, J. Sci. Comput., 18 (2003), pp. 69–81.
  * [26] , On the artificial compression method for second-order nonoscillatory central difference schemes for systems of conservation laws, SIAM J. Sci. Comput., 24 (2003), pp. 1157–1174.
  * [27] R. Liska and B. Wendrof, Comparison of several diference schemes on 1D and 2D test problems for the euler equations, SIAM J. Sci. Comput., 25 (2003), pp. 995–1017.
  * [28] H. Liu, A numerical study of the performance of alternative weighted ENO methods based on various numerical fluxes for conservation law, Appl. Math. Comput., 296 (2017), pp. 182–197.
  * [29] X. D. Liu and E. Tadmor, Third order nonoscillatory central scheme for hyperbolic conservation laws, Numer. Math., 79 (1998), pp. 397–425.
  * [30] H. Nessyahu and E. Tadmor, Nonoscillatory central differencing for hyperbolic conservation laws, J. Comput. Phys., 87 (1990), pp. 408–463.
  * [31] T. Nonomura and K. Fujii, Characteristic finite-difference WENO scheme for multicomponent compressible fluid analysis: overestimated quasi-conservative formulation maintaining equilibriums of velocity, pressure, and temperature, J. Comput. Phys., 340 (2017), pp. 358–388.
  * [32] J. Panuelos, J. Wadsley, and N. Kevlahan, Low shear difusion central schemes for particle methods, J. Comput. Phys., 414 (2020), p. 109454.
  * [33] J. Qiu and C.-W. Shu, On the construction, comparison, and local characteristic decomposition for high-order central WENO schemes, J. Comput. Phys., 183 (2002), pp. 187–209.
  * [34] C. W. Schulz-Rinne, Classification of the Riemann problem for two-dimensional gas dynamics, SIAM J. Math. Anal., 24 (1993), pp. 76–88.
  * [35] C. W. Schulz-Rinne, J. P. Collins, and H. M. Glaz, Numerical solution of the Riemann problem for two-dimensional gas dynamics, SIAM J. Sci. Comput., 14 (1993), pp. 1394–1394.
  * [36] J. Shi, Y.-T. Zhang, and C.-W. Shu, Resolution of high order WENO schemes for complicated flow structures, J. Comput. Phys., 186 (2003), pp. 690–696.
  * [37] C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes, Acta Numer., 5 (2020), pp. 701–762.
  * [38] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys., 77 (1988), pp. 439–471.
  * [39] P. K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal., 21 (1984), pp. 995–1011.
  * [40] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics: A practical introduction, Springer-Verlag, Berlin, Heidelberg, third ed., 2009.
  * [41] B.-S. Wang, W. S. Don, N. K. Garg, and A. Kurganov, Fifth-order A-WENO finite-difference schemes based on a new adaptive diffusion central numerical flux, SIAM J. Sci. Comput., 42 (2020), pp. A3932–A3956.
  * [42] B.-S. Wang, P. Li, Z. Gao, and W. S. Don, An improved fifth order alternative WENO-Z finite difference scheme for hyperbolic conservation laws, J. Comput. Phys., 374 (2018), pp. 469–477.
  * [43] P. Woodward and P. Colella, The numerical solution of two-dimensional fluid flow with strong shocks, J. Comput. Phys., 54 (1988), pp. 115–173.
  * [44] Y. Zheng, Systems of conservation laws. Two-dimensional Riemann problems, Progress in Nonlinear Differential Equations and their Applications, Birkhäuser Boston, Inc., Boston, MA, 2001.