# New Adaptive Low-Dissipation Central-Upwind Schemes

Shaoshuai Chu and Alexander Kurganov 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>
###### Abstract

We introduce new second-order adaptive low-dissipation central-upwind (LDCU)
schemes for the one- and two-dimensional hyperbolic systems of conservation
laws. The new adaptive LDCU schemes employ the LDCU numerical fluxes (recently
proposed in [A. Kurganov and R. Xin, J. Sci. Comput., 96 (2023), Paper No.
56]) computed using the point values reconstructed with the help of adaptively
selected nonlinear limiters. To this end, we use a smoothness indicator to
detect “rough” parts of the computed solution, where the piecewise linear
reconstruction is performed using an overcompressive limiter, which leads to
extremely sharp resolution of shock and contact waves. In the “smooth” areas,
we use a more dissipative limiter to prevent appearance of artificial kinks
and staircase-like structures there. In order to avoid oscillations, we
perform the reconstruction in the local characteristic variables obtained
using the local characteristic decomposition. We test two different smoothness
indicators and apply the developed schemes to the one- and two-dimensional
Euler equations of gas dynamics. The obtained numerical results clearly
demonstrate that the new adaptive LDCU schemes outperform the original ones.

Key words: Low-dissipation central-upwind schemes, minmod-based smoothness
indicator, weak local residual, overcompressive limiters, dissipative
limiters, Euler equations of gas dynamics.

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

##  1 Introduction

This paper focuses on developing new adaptive numerical methods for the
hyperbolic systems of conservation laws, which in the one- (1-D) and two-
dimensional (2-D) cases, read as

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

and

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

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

It is well-known that even when the initial data are smooth, solutions of
(1.1) and (1.2) can produce extremely complex nonsmooth wave patterns
including shocks, rarefactions, and contact discontinuities. This makes it
quite challenging to develop accurate and reliable shock-capturing numerical
methods for (1.1) and (1.2).

A library of numerical methods for the studied systems have been introduced
since the pioneering works of Friedrichs [9], Lax [28], and Godunov [13]. We
refer the reader to the monographs and review papers [22, 43, 29, 18, 2, 39]
and references therein, where one can find a description of many existing
numerical methods. In this paper, we restrict our consideration to semi-
discrete finite-volume (FV) methods, where the solution, represented in terms
of its cell averages, is evolved in time with the help of the numerical
fluxes, computed, in turn, using the reconstructed point values of $\bm{U}$ at
the boundaries of the FV cells. Many of such schemes are upwind in the sense
that their numerical fluxes are based on either exact or approximate solution
of the (generalized) Riemann problems arising at each cell interface. We,
however, focus on the Riemann-solver-free central-upwind (CU) schemes, which
provide one with accurate, efficient and robust tools for a wide variety of
hyperbolic systems. The CU schemes belong to the class of non-oscillatory
central schemes, but they have a certain upwind nature as they rely on the
local one-sided speeds of propagation, which can be estimated using the
largest and smallest eigenvalues of the corresponding Jacobians. The original
CU schemes from [24, 26] contain relatively large amount of numerical
dissipation, which was reduced in [25] and recently in [27], where built-in
“anti-diffusion” terms were introduced. The amount of numerical dissipation
can be also reduced by applying the local characteristic decomposition (LCD)
technique to the numerical diffusion of the CU fluxes; see [3].

In this paper, we use the low-dissipation CU (LDCU) numerical fluxes from
[27], and further enhance the resolution of the “rough” parts of the computed
solution by applying a new scheme adaption approach: The point values used to
evaluate the LDCU fluxes are reconstructed with the help of adaptively
selected nonlinear limiters, which are, in general, required to make the
reconstructed point values non-oscillatory. A variety of limiters are
available; see, e.g., [30, 32, 42, 2, 18, 29, 43] and references therein. Many
of the limiters can be classified as dissipative, compressive, or
overcompressive as it was done in [30]. The use of compressive and
overcompressive limiters leads to very sharp resolution of discontinuous parts
of the approximated solution, while dissipative limiters may smear the jumps.
At the same time, applying compressive and overcompressive limiters in the
smooth areas typically results in the artificial sharpening of the smooth
solution profiles, that is, in the appearance of kinks or staircase-like
structures, or even non-physical jumps.

We therefore switch between different limiters. To this end, we need to
automatically detect “rough” (nonsmooth) parts of the computed solution with
the help of a smoothness indicator (SI). Many different SIs are readily
available; see, e.g., [7, 4, 1, 12, 11, 16, 33, 35, 44, 10] and references
therein. In this paper, we test two different SIs: a slightly modified minmod
(MM)-based shock indicator from [47] (see also [17, 41]) and a SI based on
weak local residuals (WLR) from [23] (see also [20, 4, 21]). In the areas
identified as being “rough”, we use the overcompressive SBM limiters from
[30], while switching to the dissipative Minmod2 limiter elsewhere. It is
well-known that the use of any of these two limiters may lead to numerical
oscillations in the vicinities of shock and contact discontinuities. In order
to reduce these oscillations, we perform the reconstruction in the local
characteristic variables rather than in the conservative or primitive ones
(this strategy was advocated in, e.g., [34]). We switch to the characteristic
variables using the LCD, which is often used in the context of high-order
schemes, but can also be implemented to enhance the resolution of second-order
schemes; see, e.g., [3, 19, 34, 39] and references therein.

The paper is organized as follows. In §2, we review the recently proposed 1-D
LDCU scheme from [27]. We then introduce the adaptive schemes that employ
either the MM- or WLR-based SI to detect the “rough” areas. In §3, we extend
the proposed adaptive LDCU scheme to the 2-D case. In §4, we apply the
developed schemes to a number of 1-D and 2-D numerical examples for the Euler
equations of gas dynamics. We demonstrate that the adaptive LDCU schemes
contain substantially smaller amount of numerical dissipation and achieve much
higher resolution compared with the LDCU schemes based on the Minmod2 limiters
applied throughout the entire computational domain. Finally, we give some
concluding remarks in §5.

##  2 One-Dimensional Scheme Adaption Algorithm

In this section, we consider the 1-D conservation laws (1.1) and describe the
1-D adaptive algorithm.

###  2.1 1-D Low-Dissipation Central-Upwind (LDCU) Schemes

Assume that the computational domain is covered with the uniform cells
$C_{j}:=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ with
$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$ 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 suppose that at a certain time $t\geq 0$, the point values of the computed
solution $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}(t)$ are available.
Note that all of the indexed quantities are time-dependent, but from here on,
we will suppress the time-dependence of all of the indexed quantities for the
sake of brevity.

According to the semi-discrete LDCU scheme from [27], the computed cell
averages are evolved in time by numerically 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{{\cal F}}_{j+\frac{1}{2}}-\bm{{\cal
F}}_{j-\frac{1}{2}}}{\Delta x},$ (2.1)

where $\bm{{\cal F}}_{j+\frac{1}{2}}$ are the LDCU numerical fluxes defined by

$\bm{{\cal
F}}_{j+\frac{1}{2}}\big{(}\bm{U}_{j+\frac{1}{2}}^{-},\bm{U}_{j+\frac{1}{2}}^{+}\big{)}=\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)+\bm{q}_{j+\frac{1}{2}}.$

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 point values of
$\bm{U}$ at the cell interface $x=x_{j+\frac{1}{2}}$. The point values
$\bm{U}^{\pm}_{j+\frac{1}{2}}$ are 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}\\}$ using a proper
nonlinear limiter; see §2.1.1. 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})$,
$\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\big{\\{}\lambda_{d}\big{(}A(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\lambda_{d}\big{(}A(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},0\big{\\}},$
$\displaystyle
a^{-}_{j+\frac{1}{2}}=\min\big{\\{}\lambda_{1}\big{(}A(\bm{U}^{+}_{j+\frac{1}{2}})\big{)},\lambda_{1}\big{(}A(\bm{U}^{-}_{j+\frac{1}{2}})\big{)},0\big{\\}}.$

Finally, $\bm{q}_{j+\frac{1}{2}}$ is a built-in “anti-diffusion” term, which
can be derived for a particular system at hand. For instance, we consider the
1-D Euler equations of gas dynamics, which read as (1.1) with

$\bm{U}=\big{(}\rho,\rho u,E\big{)}^{\top}\quad{\rm
and}\quad\bm{F}=\big{(}\rho u,\rho u^{2}+p,u(E+p)\big{)}^{\top}.$ (2.2)

Here, $\rho$, $u$, $p$, and $E$ are the density, velocity, pressure, and total
energy, respectively, and the system is completed through the following
equations of state (EOS) for ideal gases:

$p=(\gamma-1)\Big{[}E-{\frac{1}{2}}\rho u^{2}\Big{]},$ (2.3)

where the parameter $\gamma$ represents the specific heat ratio. For the Euler
system (1.1), (2.2), (2.3), the “anti-diffusion” term $\bm{q}_{j+\frac{1}{2}}$
has been rigorously derived in [27] and it is given by

$\bm{q}_{j+\frac{1}{2}}={\rm
minmod}\big{(}-a^{-}_{j+\frac{1}{2}}(\rho^{*}_{j+\frac{1}{2}}-\rho^{-}_{j+\frac{1}{2}}),a^{+}_{j+\frac{1}{2}}(\rho^{+}_{j+\frac{1}{2}}-\rho^{*}_{j+\frac{1}{2}})\big{)}\begin{pmatrix}1\\\
u^{*}_{j+\frac{1}{2}}\\\\[3.87495pt]
{\frac{1}{2}}\big{(}u^{*}_{j+\frac{1}{2}}\big{)}^{2}\end{pmatrix}.$

Here, $\rho^{*}_{j+\frac{1}{2}}$ and $(\rho u)^{*}_{j+\frac{1}{2}}$ are the
first and second components of

$\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}}-\left\\{\bm{F}(\bm{U}^{+}_{j+\frac{1}{2}})-\bm{F}(\bm{U}^{-}_{j+\frac{1}{2}})\right\\}}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}},$

$u^{*}_{j+\frac{1}{2}}=(\rho u)^{*}_{j+\frac{1}{2}}/\rho^{*}_{j+\frac{1}{2}}$,
and the minmod function is defined by

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

####  2.1.1 Nonlinear Limiters

As mentioned before, the point values $\bm{U}^{\pm}_{j+\frac{1}{2}}$ are
obtained with the help of a conservative piecewise linear reconstruction,
designed using a proper nonlinear limiter. In this paper, we use a family of
the SBM limiters (introduced in [30]) applied to the local characteristic
variables. To this end, we introduce the matrices
$\widehat{A}_{j+\frac{1}{2}}=A\big{(}(\,\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\big{)}$ and 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}}$ are
diagonal matrices. We then introduce the local characteristic variables
$\bm{\Gamma}$ 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,\,j,\,j+1,\,j+2.$

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

$(\bm{\Gamma}_{x})_{j}=\phi^{\rm
SBM}_{\theta,\tau}\left(\frac{\bm{\Gamma}_{j+1}-\bm{\Gamma}_{j}}{\bm{\Gamma}_{j}-\bm{\Gamma}_{j-1}}\right)\frac{\bm{\Gamma}_{j}-\bm{\Gamma}_{j-1}}{\Delta
x}$ (2.4)

and

$(\bm{\Gamma}_{x})_{j+1}=\phi^{\rm
SBM}_{\theta,\tau}\left(\frac{\bm{\Gamma}_{j+2}-\bm{\Gamma}_{j+1}}{\bm{\Gamma}_{j+1}-\bm{\Gamma}_{j}}\right)\frac{\bm{\Gamma}_{j+1}-\bm{\Gamma}_{j}}{\Delta
x},$ (2.5)

where the two-parameter SBM function

$\phi^{\rm SBM}_{\theta,\tau}(r):=\begin{cases}0&\mbox{if $r<0$},\\\
\min\\{r\theta,1+\tau(r-1)\\}&\mbox{if}~{}0<r\leq 1,\\\ r\phi^{\rm
SBM}_{\theta,\tau}(\frac{1}{r})&\text{otherwise,}\end{cases}$ (2.6)

is applied in the component-wise manner.

The parameters $\theta\in[1,2]$ and $\tau$ in (2.6) can be used to control the
amount of numerical dissipation present in the resulting scheme. First, larger
$\theta$’s correspond to less dissipative but, in general, more oscillatory
reconstructions. In all of the numerical examples reported in §4, we have
taken $\theta=2$. Second, according to [30], if $\tau\geq 0.5$, then the SBM
limiter is dissipative and its use typically causes contact discontinuities to
be severely smeared in time. If $0\leq\tau<0.5$, then the SBM limiter is
compressive and in this case, contact waves are usually resolved sharply
within few points, but smooth extrema might be slightly compressed resulting
in continuous solution profiles having a kink. If $\tau<0$, then the limiter
is overcompressive so that contact discontinuities typically stay very sharp
for long time, while smooth solutions become overcompressed as time evolves
resulting in the appearance of artificial ${\cal O}(1)$ jump discontinuities.

Equipped with (2.4) and (2.5), we 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 then 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}}.$

###### Remark 2.1

For detailed explanations on how the 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, we refer the reader to [3, Appendix A].

###  2.2 One-Dimensional Adaptive Schemes

We now turn to the description of the proposed adaptive schemes. The key
ingredient of the new schemes is the use of the different limiters from the
family (2.6) in different parts of the computational domain. In particular, we
use an overcompressive limiter with $\tau=-0.25$ in the “rough” parts of the
computed solution and a dissipative limiter with $\tau=0.5$ elsewhere. The
latter limiter is, in fact, the Minmod2 limiter, which can be written in a
simpler form since (2.4) and (2.5) with $\theta=2$ and $\tau=0.5$ reduce to

$(\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),$

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

In order to implement this simple scheme adaption approach, we need to
automatically detect “rough” parts of the computed solution. This is done
using either the MM- or WLR-based SIs briefly described in §2.2.1 and §2.2.2
below.

####  2.2.1 Minmod-Based Smoothness Indicator

We first compute the MM-based quantities

$s_{j}={\rm minmod}\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j+1}-\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}-\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern
0.0pt}}}_{j-1}\big{)},$

and then we say that the cell $C_{j}$ is “rough” if
$|s_{j}|>\max\big{\\{}|s_{j-1}|,|s_{j+1}|\big{\\}}+\delta$, where $\delta$ is
a small positive number. The parameter $\delta$ has to be selected for each
problem at hand and it should indicate a size of a jump in $\rho$, which we
neglect when detecting “rough” parts of the solution. In fact, this SI is not
very sensitive to the choice of $\delta$ and in all of the numerical examples
reported in §4, we have taken $\delta=10^{-4}$.

####  2.2.2 Weak Local Residual-Based Smoothness Indicator

In order to detect “rough” areas, one can also use the WLR-based SI, which we
obtain as follows. First, we assume that the cell averages $\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j}$
are available at a certain time level $t=t^{n}$ and the two previous time
levels $t=t^{n-1}$ (with $t^{n}-t^{n-1}=\Delta t^{n-1}$) and $t^{n-2}$ (with
$t^{n-1}-t^{n-2}=\Delta t^{n-2}$). In addition, we assume that the solution
has been reconstructed at $t=t^{n-1}$ and $t^{n-2}$ and the corresponding
point values at the cell interfaces $x=x_{j+\frac{1}{2}}$ are available. We
will denote these reconstructed point values by $\bm{U}^{n-1}_{j+\frac{1}{2}}$
and $\bm{U}^{n-2}_{j+\frac{1}{2}}$. Recall that we, in fact, obtain two point
values at each cell interface ($\bm{U}^{\pm}_{j+\frac{1}{2}}$) and any of them
can be used to evaluate the WLRs for the density equation $\rho_{t}+(\rho
u)_{x}=0$. We denote these WLRs by
$\varepsilon^{n-\frac{3}{2}}_{j+\frac{1}{2}}$ and compute them according to
[23]:

$\displaystyle\varepsilon^{n-\frac{3}{2}}_{j+\frac{1}{2}}$
$\displaystyle=\frac{\Delta
x}{6}\left[\rho^{n-1}_{j+\frac{3}{2}}-\rho^{n-2}_{j+\frac{3}{2}}+4\big{(}\rho^{n-1}_{j+\frac{1}{2}}-\rho^{n-2}_{j+\frac{1}{2}}\big{)}+\rho^{n-1}_{j-\frac{1}{2}}-\rho^{n-2}_{j-\frac{1}{2}}\right]$
(2.7) $\displaystyle+\frac{\Delta t^{n-2}}{4}\left[(\rho
u)^{n-1}_{j+\frac{3}{2}}-(\rho u)^{n-1}_{j-\frac{1}{2}}+(\rho
u)^{n-2}_{j+\frac{3}{2}}-(\rho u)^{n-2}_{j-\frac{1}{2}}\right].$

The desired SIs are then obtained at each cell interface $x=x_{j+\frac{1}{2}}$
by setting

$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\varepsilon$\kern
0.0pt}}}^{\,n-\frac{3}{2}}_{j+\frac{1}{2}}:=\frac{1}{6}\left[\varepsilon^{\,n-\frac{3}{2}}_{j-\frac{1}{2}}+4\varepsilon^{\,n-\frac{3}{2}}_{j+\frac{1}{2}}+\varepsilon^{\,n-\frac{3}{2}}_{j+\frac{3}{2}}\right].$

As mentioned in [23], the size of the WLRs and thus of the SIs for second-
order schemes are expected to be

$||\varepsilon^{n-\frac{3}{2}}||_{\infty}\sim\begin{cases}\Delta&\mbox{near
shock waves,}\\\ \Delta^{\alpha}&\mbox{near contact,}\\\ \Delta^{4}&\mbox{in
smooth regions},\end{cases}$ (2.8)

where $\Delta:=\max\\{\Delta t,\Delta x\\}$ and $1<\alpha\leq 2$; see also
[20].

Finally, we take advantage of (2.8), which suggests that the size of SIs
ranges from ${\cal O}(\Delta)$ near shocks to ${\cal O}(\Delta^{4})$ in the
smooth regions and develop the following simple strategy for the automatic
detection of “rough” parts of the computed solution $\\{\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j}\\}$. We mark the cell $C_{j}$ as “rough” as long as

$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\varepsilon$\kern
0.0pt}}}^{\,n-\frac{3}{2}}_{j}>\texttt{C}(\Delta x)^{2},$ (2.9)

where C is a positive tunable constant to be selected for each problem at
hand. The robustness of this shock detection strategy depends on the
sensitivity of the proposed algorithm to the choice of C. In principle, C can
be tuned on a coarse mesh and then used for fine mesh computations, but as we
demonstrate in Examples 2 in §4, this approach may fail. Therefore, even
though the use of the WLR-based SI may lead to extremely sharp results (like
in, for instance, Example 3 in §4), the adaption strategy that relies on this
SI may not be robust.

###### Remark 2.2

While implementing (2.7), we have used
$\bm{U}_{j+\frac{1}{2}}=\bm{U}^{-}_{j+\frac{1}{2}}$ in all the numerical
examples reported in §4. Other choices like
$\bm{U}_{j+\frac{1}{2}}=\bm{U}^{+}_{j+\frac{1}{2}}$ or
$\bm{U}_{j+\frac{1}{2}}=\big{(}\bm{U}^{-}_{j+\frac{1}{2}}+\bm{U}^{+}_{j+\frac{1}{2}}\big{)}/2$
can also be used with no visible advantages or disadvantages of any of them.

##  3 Two-Dimensional Scheme Adaption Algorithm

In this section, we extend the 1-D adaptive strategy introduced in §2 to the
2-D hyperbolic systems of conservation laws (1.2).

###  3.1 2-D Low-Dissipation Central-Upwind (LDCU) Schemes

Let the computational domain be covered with uniform 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}}]$
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}=(y_{k-\frac{1}{2}}+y_{k+\frac{1}{2}})/2$. We assume that 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\limits_{C_{j,k}}\bm{U}(x,y,t)\,{\rm d}y\,{\rm d}x,$

have been computed at a certain time $t\geq 0$.

According to the semi-discrete LDCU scheme from [27], the computed cell
averages are evolved in time by numerically solving the following system of
ODEs:

$\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}}_{{j+\frac{1}{2}},k}-\bm{\mathcal{F}}_{{j-\frac{1}{2}},k}}{\Delta
x}-\frac{\bm{\mathcal{G}}_{j,{k+\frac{1}{2}}}-\bm{\mathcal{G}}_{j,{k-\frac{1}{2}}}}{\Delta
y},$ (3.1)

where $\bm{{\cal F}}_{{j+\frac{1}{2}},k}=\bm{{\cal
F}}_{{j+\frac{1}{2}},k}\big{(}\bm{U}_{{j+\frac{1}{2}},k}^{-},\bm{U}_{{j+\frac{1}{2}},k}^{+}\big{)}$
and $\bm{{\cal G}}_{j,{k+\frac{1}{2}}}=\bm{{\cal
G}}_{j,{k+\frac{1}{2}}}\big{(}\bm{U}_{j,{k+\frac{1}{2}}}^{-},\bm{U}_{j,{k+\frac{1}{2}}}^{+}\big{)}$
are the LDCU numerical fluxes defined by

$\displaystyle\bm{\mathcal{F}}_{{j+\frac{1}{2}},k}$
$\displaystyle=\frac{a^{+}_{{j+\frac{1}{2}},k}\bm{F}^{-}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}\bm{F}^{+}_{{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}^{x}_{{j+\frac{1}{2}},k},$
$\displaystyle\bm{\mathcal{G}}_{j,{k+\frac{1}{2}}}$
$\displaystyle=\frac{b^{+}_{j,{k+\frac{1}{2}}}\bm{G}^{-}_{j,{k+\frac{1}{2}}}-b^{-}_{j,{k+\frac{1}{2}}}\bm{G}^{+}_{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}^{y}_{j,{k+\frac{1}{2}}}.$

Here,
$\bm{F}^{\pm}_{{j+\frac{1}{2}},k}:=\bm{F}\big{(}\bm{U}^{\pm}_{{j+\frac{1}{2}},k}\big{)}$
and
$\bm{G}^{\pm}_{j,{k+\frac{1}{2}}}:=\bm{G}\big{(}\bm{U}^{\pm}_{j,{k+\frac{1}{2}}}\big{)}$,
and $\bm{U}^{\pm}_{{j+\frac{1}{2}},k}$ and $\bm{U}^{\pm}_{j,{k+\frac{1}{2}}}$
are the one-sided point values of $\bm{U}$ at the cell interfaces
$(x_{j+\frac{1}{2}},y_{k})$ and $(x_{j},y_{k+\frac{1}{2}})$, respectively. We
reconstruct the point values $\bm{U}^{\pm}_{{j+\frac{1}{2}},k}$ and
$\bm{U}^{\pm}_{j,{k+\frac{1}{2}}}$ using the LCD; see Appendix A for details.
The one-sided local speeds of propagation in the $x$\- and $y$-directions,
$a^{\pm}_{{j+\frac{1}{2}},k}$ and $b^{\pm}_{j,{k+\frac{1}{2}}}$, can be
estimated by the largest and smallest eigenvalues of the Jacobians
$A(\bm{U}):=\frac{\partial\bm{F}}{\partial\bm{U}}(\bm{U})$ and
$B(\bm{U}):=\frac{\partial\bm{G}}{\partial\bm{U}}(\bm{U})$, for example, by
setting

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

Finally, $\bm{q}^{x}_{{j+\frac{1}{2}},k}$ and $\bm{q}^{y}_{j,{k+\frac{1}{2}}}$
are built-in “anti-diffusion” terms, which can be derived for a particular
system (1.2) at hand. For instance, we consider the 2-D Euler equations of gas
dynamics, which read as (1.2) with

$\bm{U}=\big{(}\rho,\rho u,\rho v,E\big{)}^{\top},~{}~{}\bm{F}=\big{(}\rho
u,\rho u^{2}+p,\rho uv,u(E+p)\big{)}^{\top},~{}~{}\bm{G}=\big{(}\rho v,\rho
uv,\rho v^{2}+p,v(E+p)\big{)}^{\top},$ (3.2)

where $v$ is the $y$-velocity and the other variables are as the same as in
(2.2). The system (3.2) is completed through the following EOS for ideal
gases:

$p=(\gamma-1)\Big{[}E-\frac{\rho}{2}(u^{2}+v^{2})\Big{]}.$ (3.3)

For the sake of brevity, we omit the details on the built-in “anti-diffusion”
terms $\bm{q}^{x}_{{j+\frac{1}{2}},k}$ and $\bm{q}^{y}_{j,{k+\frac{1}{2}}}$.
For the system (1.2), (3.2)–(3.3), they have been derived in [27].

###  3.2 Two-Dimensional Adaptive Schemes

We now turn to the description of the proposed adaptive schemes for the 2-D
system. As in the 1-D case, we use an overcompressive SBM limiter with
$\tau=-0.25$ in the “rough” parts of the computed solution and a dissipative
Minmod2 limiter elsewhere. To this end, we detect the “rough” parts of the
numerical solution using either the MM- or WLR-based SIs described in §3.2.1
and §3.2.2 below.

####  3.2.1 Two-Dimensional Minmod-Based Smoothness Indicator

The 1-D MM-based SI introduced in §2.2.1, is extended to the 2-D case in the
“dimension-by-dimension” manner. We first compute the MM-based quantities in
the $x$-direction,

$s^{x}_{j,k}={\rm
minmod}\big{(}\rho_{j+1,k}-\rho_{j,k},\rho_{j,k}-\rho_{j-1,k}\big{)},$

and use the overcompressive SBM limiter to compute the slopes in the
$x$-direction only in those cell $C_{j,k}$, where
$|s^{x}_{j,k}|>\max\big{\\{}|s^{x}_{j-1,k}|,|s^{x}_{j+1,k}|\big{\\}}+\delta$.
Similarly, we compute the MM-based quantities in the $y$-direction,

$s^{y}_{j,k}={\rm
minmod}\big{(}\rho_{j,k+1}-\rho_{j,k},\rho_{j,k}-\rho_{j,k-1}\big{)},$

and use the overcompressive SBM limiter to compute the slopes in the
$y$-direction only in those cell $C_{j,k}$, where
$|s^{y}_{j,k}|>\max\big{\\{}|s^{y}_{j,k-1}|,|s^{y}_{j,k+1}|\big{\\}}+\delta$.

####  3.2.2 Two-Dimensional Weak Local Residual-Based Smoothness Indicator

One can also detect the “rough” parts of the numerical solution using the 2-D
WLR-based SI, which we obtain as follows. First, we compute the WLRs
introduced in [23]. For the 2-D density equation $\rho_{t}+(\rho u)_{x}+(\rho
v)_{y}=0$, these WLRs are

$\varepsilon^{n-\frac{3}{2}}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}=\frac{1}{36\Delta}\Delta
x\Delta y\,{\cal
U}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}^{n-\frac{3}{2}}+\frac{1}{12\Delta}\Big{(}\Delta
y\Delta t{\cal F}_{j,k}^{n-\frac{3}{2}}+\Delta x\Delta t{\cal
G}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}^{n-\frac{3}{2}}\Big{)},$

where $\Delta:=\max\\{\Delta t,\Delta x,\Delta y\\}$ and

$\displaystyle{\cal U}_{{j+\frac{1}{2}},\,{k+\frac{1}{2}}}^{n-\frac{3}{2}}$
$\displaystyle=\left[\rho^{n-1}_{j+\frac{3}{2},k+\frac{3}{2}}-\rho^{n-2}_{j+\frac{3}{2},k+\frac{3}{2}}+\rho^{n-1}_{j+\frac{3}{2},{k-\frac{1}{2}}}-\rho^{n-2}_{j+\frac{3}{2},{k-\frac{1}{2}}}+\rho^{n-1}_{{j-\frac{1}{2}},k+\frac{3}{2}}-\rho^{n-2}_{{j-\frac{1}{2}},k+\frac{3}{2}}\right.$
$\displaystyle\left.+\rho^{n-1}_{{j-\frac{1}{2}},{k-\frac{1}{2}}}-\rho^{n-2}_{{j-\frac{1}{2}},{k-\frac{1}{2}}}\right]+4\left[\rho^{n-1}_{j+\frac{3}{2},{k+\frac{1}{2}}}-\rho^{n-2}_{j+\frac{3}{2},{k+\frac{1}{2}}}+\rho^{n-1}_{{j-\frac{1}{2}},{k+\frac{1}{2}}}-\rho^{n-2}_{{j-\frac{1}{2}},{k+\frac{1}{2}}}\right.$
$\displaystyle\left.+\rho^{n-1}_{{j+\frac{1}{2}},k+\frac{3}{2}}-\rho^{n-2}_{{j+\frac{1}{2}},k+\frac{3}{2}}+\rho^{n-1}_{{j+\frac{1}{2}},{k-\frac{1}{2}}}-\rho^{n-2}_{{j+\frac{1}{2}},{k-\frac{1}{2}}}\right]+16\left[\rho^{n-1}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}-\rho^{n-2}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}\right],$
$\displaystyle{\cal F}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}^{n-\frac{3}{2}}$
$\displaystyle=\left[(\rho u)^{n-1}_{j+\frac{3}{2},k+\frac{3}{2}}-(\rho
u)^{n-1}_{{j-\frac{1}{2}},k+\frac{3}{2}}+(\rho
u)^{n-1}_{j+\frac{3}{2},{k-\frac{1}{2}}}-(\rho
u)^{n-1}_{{j-\frac{1}{2}},{k-\frac{1}{2}}}\right.$ $\displaystyle\left.+(\rho
u)^{n-2}_{j+\frac{3}{2},k+\frac{3}{2}}-(\rho
u)^{n-2}_{{j-\frac{1}{2}},k+\frac{3}{2}}+(\rho
u)^{n-2}_{j+\frac{3}{2},{k-\frac{1}{2}}}-(\rho
u)^{n-2}_{{j-\frac{1}{2}},{k-\frac{1}{2}}}\right]$ $\displaystyle+4\left[(\rho
u)^{n-1}_{j+\frac{3}{2},{k+\frac{1}{2}}}-(\rho
u)^{n-1}_{{j-\frac{1}{2}},{k+\frac{1}{2}}}-(\rho
u)^{n-2}_{j+\frac{3}{2},{k+\frac{1}{2}}}-(\rho
u)^{n-2}_{{j-\frac{1}{2}},{k+\frac{1}{2}}}\right],$ $\displaystyle{\cal
G}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}^{n-\frac{3}{2}}$
$\displaystyle=\left[(\rho v)^{n-1}_{j+\frac{3}{2},k+\frac{3}{2}}-(\rho
v)^{n-1}_{j+\frac{3}{2},{k-\frac{1}{2}}}+(\rho
v)^{n-1}_{{j-\frac{1}{2}},k+\frac{3}{2}}-(\rho
v)^{n-1}_{{j-\frac{1}{2}},{k-\frac{1}{2}}}\right.$ $\displaystyle\left.+(\rho
v)^{n-2}_{j+\frac{3}{2},k+\frac{3}{2}}-(\rho
v)^{n-2}_{j+\frac{3}{2},{k-\frac{1}{2}}}+(\rho
v)^{n-2}_{{j-\frac{1}{2}},k+\frac{3}{2}}-(\rho
v)^{n-2}_{{j-\frac{1}{2}},{k-\frac{1}{2}}}\right]$ $\displaystyle+4\left[(\rho
v)^{n-1}_{{j+\frac{1}{2}},k+\frac{3}{2}}-(\rho
v)^{n-1}_{{j+\frac{1}{2}},{k-\frac{1}{2}}}-(\rho
v)^{n-2}_{{j+\frac{1}{2}},k+\frac{3}{2}}-(\rho
v)^{n-2}_{{j+\frac{1}{2}},{k-\frac{1}{2}}}\right].$

Note that no values at the cell corners
$(x_{j+\frac{1}{2}},y_{k+\frac{1}{2}})$ are reconstructed. We therefore set

$\bm{U}^{\ell}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}:=\frac{1}{4}\left[\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,\ell}_{j,k}+\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,\ell}_{j+1,k}+\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,\ell}_{j,k+1}+\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,\ell}_{j+1,k+1}\right],\quad\ell=n-2,\,n-1,$

where $\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,n-2}_{j,k}$ and
$\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,n-1}_{j,k}$ denote the
cell averages computed at times $t=t^{n-2}$ and $t=t^{n-1}$, respectively.

The desired SIs are then obtained at each cell corner
$(x_{j+\frac{1}{2}},y_{k+\frac{1}{2}})$ using a proper averaging, for
instance, by

$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\varepsilon$\kern
0.0pt}}}^{\,n-\frac{3}{2}}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}:$
$\displaystyle=\frac{1}{36}\left[\varepsilon^{n-\frac{3}{2}}_{{j-\frac{1}{2}},{k-\frac{1}{2}}}+\varepsilon^{n-\frac{3}{2}}_{{j-\frac{1}{2}},k+\frac{3}{2}}+\varepsilon^{n-\frac{3}{2}}_{j+\frac{3}{2},{k-\frac{1}{2}}}+\varepsilon^{n-\frac{3}{2}}_{j+\frac{3}{2},k+\frac{3}{2}}\right.$
$\displaystyle\left.+\,4\left(\varepsilon^{n-\frac{3}{2}}_{{j-\frac{1}{2}},{k+\frac{1}{2}}}+\varepsilon^{n-\frac{3}{2}}_{{j+\frac{1}{2}},{k-\frac{1}{2}}}+\varepsilon^{n-\frac{3}{2}}_{{j+\frac{1}{2}},k+\frac{3}{2}}+\varepsilon^{n-\frac{3}{2}}_{j+\frac{3}{2},{k+\frac{1}{2}}}\right)+16\,\varepsilon^{n-\frac{3}{2}}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}\right].$

As mentioned in [23], the size of the WLRs and thus of the SIs for second-
order schemes are expected to be the same as in (2.8). We therefore act
similarly to (2.9) and mark the cell $C_{j,k}$ as “rough” as long as
$\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\varepsilon$\kern
0.0pt}}}^{\,n-\frac{3}{2}}_{{j+\frac{1}{2}},{k+\frac{1}{2}}}>\texttt{C}\max\\{(\Delta
x)^{2},(\Delta y)^{2}\\}$, where C is a positive tunable constant to be
selected for each problem at hand.

##  4 Numerical Examples

In this section, we test the developed adaptive schemes on several numerical
examples. To this end, we compare the performance of the original LDCU and the
adaptive LDCU schemes by applying them to a number of initial-boundary value
problems for the 1-D and 2-D Euler equations of gas dynamics. The adaptive
LDCU schemes with the MM- and WLR-based SIs used to detect the “rough” areas
will be referred to as the _A-MM_ and _A-WLR_ schemes, respectively.

In all of the numerical examples reported, we have solved the ODE systems
(2.1) and (3.1) using the three-stage third-order strong stability preserving
(SSP) Runge-Kutta method; see, e.g., [14, 15]. We take $\gamma=1.4$ in Example
1–5 and $\gamma=5/3$ in Example 6. The CFL number is 0.4 in all of the
examples.

###  4.1 One-Dimensional Examples

##### Example 1—Shock-Entropy Wave Interaction Problem.

In the first example taken from [40], we consider the shock-entropy wave
interaction problem. 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 number 1.1 interacting with
high-frequency density perturbations, that is, as the shock wave moves, the
perturbations spread ahead. We set the free boundary condition at the both
ends of the computational domain $[-10,5]$.

We compute the numerical solution until the final time $t=5$ by the LDCU,
A-MM, and A-WLR schemes on a uniform mesh with $\Delta x=1/80$. We use the
adaption constant $\texttt{C}=0.1$ in the A-WLR scheme. The numerical results
at time $t=5$ are presented in Figure 4.1 along with the reference solution
computed by the LDCU scheme on a much finer mesh with $\Delta x=1/1600$. On
the right panel of Figure 4.1, we zoom the obtained solutions at the interval
$[-1,0]$, at which the exact solution is smooth but has an oscillatory nature.
As one can see, this part of the solution is resolved much more accurately by
the two adaptive schemes, especially by the A-WLR one.

Figure 4.1: Example 1: Density $\rho$ computed by the LDCU, A-MM, and A-WLR
schemes (left) and zoom at $x\in[-1,0]$ (right).

##### Example 2—Shock-Density Wave Interaction Problem.

In the second example taken from [41], we consider the shock-density wave
interaction problem. The initial data,

$(\rho,u,p)(x,0)=\begin{cases}\Big{(}\dfrac{27}{7},\dfrac{4\sqrt{35}}{9},\dfrac{31}{3}\Big{)},&x<-4,\\\\[3.44444pt]
(1+0.2\sin(5x),0,1),&x>-4,\end{cases}$

are prescribed in the computational domain $[-5,15]$ subject to the free
boundary conditions.

We compute the numerical solutions by the LDCU, A-MM, and A-WLR schemes on the
uniform mesh with $\Delta x=1/40$ until the final time $t=5$. The A-WLR scheme
is used with the adaption constant $\texttt{C}=0.35$. We present the obtained
numerical results in Figure 4.2 together with the reference solution computed
by the LDCU scheme on a much finer mesh with $\Delta x=1/400$. It can be
clearly seen in Figure 4.2 (right) that both of the adaptive schemes produce
more accurate results compared to those obtained by the LDCU scheme. One can
also observe that in this example, unlike the previous one, the A-MM scheme
achieves higher resolution of the smooth parts of the solution compared with
its A-WLR counterpart. This is attributed to a relatively large value of C
used in this example.

Figure 4.2: Example 2: Density $\rho$ computed by the LDCU, A-MM, and A-WLR
schemes (left) and zoom at $x\in[11.8,13.6]$ (right).

Recall that one of the key points in the proposed WLR-based scheme adaption
strategy is tuning the adaption constant C. In [23], where a WLR-based
adaptive artificial viscosity was introduced and studied, the viscosity
coefficient, which is directly related to C, was first adjusted on a coarse
mesh and then used for the high-resolution computations on finer meshes.
However, this strategy does not seem to be robust in the A-WLR scheme as the
numerical results computed by the A-WLR scheme may still have staircase-like
structure in the areas where the coarse mesh solution is smooth. In order to
illustrate this, we use a slightly smaller adaption constant $\texttt{C}=0.2$
and compute the results with $\texttt{C}=0.2$ and $\texttt{C}=0.35$ on the
uniform meshes with $\Delta x=1/40$ and $1/200$. We present the obtained
densities (zoomed at $x\in[9,11]$) in Figure 4.3, where one can see that even
though the solution computed on the coarse mesh with $\texttt{C}=0.2$ is
smooth, it develops clear staircase-like structures when the mesh is refined.

Figure 4.3: Example 2: Density $\rho$ computed by the A-WLR scheme with
$\texttt{C}=0.2$ and 0.35 on a coarse (left) and fine (right) meshes. Zoom at
$x\in[9,11]$.

##### Example 3—Blast Wave Problem.

In the last 1-D example, we consider the strong shocks interaction problem
from [48], 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 the LDCU,
A-MM, and A-WLR (with the adaption constant $\texttt{C}=0.1$) schemes on a
uniform mesh with $\Delta x=1/400$ and implement the LDCU scheme on a much
finer grid with $\Delta x=1/4000$ to compute the reference solution. The
obtained results, presented in Figure 4.4, demonstrate that while the A-MM
scheme outperforms the LDCU one, the A-WLR results are even more accurate and
the A-WLR scheme is capable of achieving a superb resolution of the contact
wave located at about $x=0.6$. It is well-known that this contact wave is the
one, which is hardest to get sharply resolved, and the A-WLR scheme is, to
best of our knowledge, the first Riemann-problem solver-free scheme that can
achieve this goal.

Figure 4.4: Example 3: Density $\rho$ computed by the LDCU, A-MM, and A-WLR
schemes (left) and zoom at $x\in[0.55,0.85]$.

###  4.2 Two-Dimensional Examples

##### Example 4—2-D Riemann Problem.

In the first 2-D example, we consider Configuration 3 of the 2-D Riemann
problems from [26] (see also [36, 37, 49]) with the initial conditions

$(\rho,u,v,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}$

prescribed in the computational domain $[0,1.2]\times[0,1.2]$ subject to the
free boundary conditions.

We compute the numerical solution until the final time $t=1$ by the LDCU,
A-MM, and A-WLR (with the adaption constant $\texttt{C}=4$) schemes on the
uniform mesh with $\Delta x=\Delta y=3/2500$ and present the obtained results
in Figure 4.5, where one can see that both the A-MM and A-WLR schemes
outperform the LDCU scheme in capturing a sideband instability of the jet in
the zones of strong along-jet velocity shear and the instability along the
jets neck.

Figure 4.5: Example 4: Density $\rho$ computed by the LDCU (left), A-MM
(middle), and A-WLR (right) schemes.

In Figure 4.6, we show the regions which have been detected as “rough” by the
SIs at the final time. We first indicate (in the left and middle panels) the
regions in which the MM-based SI detected large $x$\- and $y$-directional
derivatives, respectively, and where the overcompressive SBM directional
limiters have been used. In the right panel, we show the “rough” regions
indicated by the WLR-based SI. As one can see, when the A-WLR scheme is used,
a sharper SBM limiter is implemented only in a small part of the computational
domain, mostly around the shocks.

Figure 4.6: Example 4: Areas detected as having large $x$\- (left) and
$y$-derivatives (middle) by the MM-based SI and the “rough” areas detected by
the WLR-based SI (right).

##### Example 5—Implosion Problem.

In this example, we consider the implosion problem taken from [31]. The
initial conditions,

$(\rho,u,v,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 $[0,0.5]\times[0,0.5]$ with the solid wall boundary
conditions imposed at all of the four sides. This example was designed to test
the amount of numerical diffusion present in different schemes as there is a
jet forming near the origin and propagating along the diagonal $y=x$
direction, and schemes containing large numerical diffusion may not resolve
the jet at all or the jet propagation velocity may be affected by the
numerical diffusion.

We compute the numerical solution until the final time $t=2.5$ by the LDCU,
A-MM, and A-WLR (with the adaption constant $\texttt{C}=5$) schemes on the
uniform mesh with $\Delta x=\Delta y=1/2000$ and present the obtained
numerical results in Figure 4.7. As one can observe, while the jet is
generated by all of the studied schemes, it propagates much further in the
diagonal direction when the solution is computed by one of the adaptive
schemes. In fact, the A-MM scheme seems to contain even smaller amount of
numerical diffusion than the A-WLR one. This can be confirmed by the results
presented in Figure 4.8, where we see that the “rough” areas detected at the
final time by the A-WLR scheme are concentrated along the shock waves only.
This means that the A-WLR scheme uses the overcompressive SBM limiter in a
smaller part of the computational domain compared with the A-MM scheme. It is
worth noting that the A-WLR scheme can be made less dissipative by decreasing
the adaption constant C, but it might be difficult to tune C in this example
as for smaller values of C the A-WLR scheme may produce reasonably sharp
results on a coarse mesh, but it may develop severe instabilities when the
mesh is refined to $\Delta x=\Delta y=1/2000$.

Figure 4.7: Example 5: Density $\rho$ computed by the LDCU (left), A-MM
(middle), and A-WLR (right) schemes.

Figure 4.8: Example 5: Areas detected as having large $x$\- (left) and
$y$-derivatives (middle) by the MM-based SI and the “rough” areas detected by
the WLR-based SI (right).

###### Remark 4.1

The solution of the studied initial-boundary value problem is symmetric with
respect to the axis $y=x$, but this symmetry may be destroyed by the roundoff
errors when the solution is computed by the studied low-dissipation schemes.
In order to prevent the loss of symmetry, we have used a very simple strategy
introduced in [46]: upon completion of each time evolution step, we replace
the computed values $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,\,k}$ with
$\widehat{\bm{U}}_{j,k}$, where

$\widehat{\rho}_{j,k}:=\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}}}_{k,j}}{2},~{}~{}(\widehat{\rho 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 v$\kern
0.0pt}}})_{k,j}}{2},~{}~{}(\widehat{\rho 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 u$\kern
0.0pt}}})_{k,j}}{2},~{}~{}\widehat{E}_{j,k}:=\frac{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}_{j,k}+\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}_{k,j}}{2},$

for all $j,k$. For more sophisticated symmetry enforcement techniques, we
refer the reader to, e.g., [5, 8, 6, 45].

##### Example 6—RT Instability.

In the last example taken from [38], 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. To this end, we first modify the 2-D Euler
equations of gas dynamics (1.2), (3.2)–(3.3) by adding the gravitational
source terms acting in the positive direction of the $y$-axis into the RHS of
the system:

$\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}+\left[u(E+p)\right]_{x}+\left[v(E+p)\right]_{y}=\rho v,$

and then use the following initial conditions:

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

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

$(\rho,u,v,p)(x,1,t)=(1,0,0,2.5),\quad(\rho,u,v,p)(x,0,t)=(2,0,0,1).$

We compute the numerical solution until the final time $t=2.95$ by the LDCU,
A-MM, and A-WLR (with the adaption constant $\texttt{C}=3$) schemes on the
uniform mesh with $\Delta x=\Delta y=1/1024$ in the computational domain
$[0,0.25]\times[0,1]$ and then present the numerical results obtained at the
times $t=1.95$ and 2.95 in Figure 4.9. As one can see, the A-MM and A-WLR
schemes resolve more small structures than the LDCU scheme, which again
demonstrates that using a sharper SBM limiter can produce sharper numerical
results. As in the previous examples, one can observe that the A-MM scheme
uses the overcompressive SBM limiter in a larger number of cells compared with
the A-WLR scheme; see Figure 4.10, where the indicated “rough” areas are
plotted.

Figure 4.9: Example 6: Density $\rho$ computed by the LDCU (left), A-MM
(middle), and A-WLR (right) schemes at $t=1.95$ (top row) and 2.95 (bottom
row).

Figure 4.10: Example 6: Areas detected as having large $x$\- (left) and
$y$-derivatives (middle) by the MM-based SI and the “rough” areas detected by
the WLR-based SI (right) at $t=1.95$ (top row) and 2.95 (bottom row).

###### Remark 4.2

In this example, the solution is symmetric with respect to the vertical axis
$x=0.125$. In order to enforce this symmetry, we have applied the strategy
from [46]: upon completion of each time evolution step, we replace the
computed cell averages $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}$ with
$\widehat{\bm{U}}_{j,k}$, where

$\displaystyle\widehat{\rho}_{j,k}$ $\displaystyle=\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}}}_{M-j,k}}{2},$
$\displaystyle(\widehat{\rho u})_{j,k}$
$\displaystyle=\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 u$\kern 0.0pt}}})_{M-j,k}}{2},$
$\displaystyle(\widehat{\rho v})_{j,k}$
$\displaystyle=\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 v$\kern 0.0pt}}})_{M-j,k}}{2},$
$\displaystyle\widehat{E}_{j,k}$ $\displaystyle=\frac{\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$E$\kern
0.0pt}}}_{j,k}+\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}_{M-j,k}}{2},$

for all $j=1,\ldots,M$ and for all $k$. Alternative symmetry enforcement
techniques can be found in, e.g., [5, 6, 8, 45].

##  5 Conclusion

In this paper, we have introduced new second-order adaptive low-dissipation
central-upwind schemes for the one- (1-D) and two-dimensional (2-D) hyperbolic
systems of conservation laws. The new adaptive schemes are based on the
recently proposed low-dissipation central-upwind (LDCU) fluxes and two
smoothness indicators (SIs) (the minmod (MM)- and weak local residual
(WLR)-based ones) used to automatically detect “rough” areas of the computed
solutions. We then use the overcompressive SBM limiters in the “rough” areas
and the dissipative Minmod2 limiters elsewhere to achieve higher resolution of
the computed shocks and contact discontinues and, at the same time, to avoid
the staircase-like overcompressed structures in the computed results. We have
applied the developed adaptive schemes to the 1-D and 2-D Euler equations of
gas dynamics and the obtained numerical results clearly demonstrate that both
of the adaptive schemes outperform the original LDCU scheme.

We have also compared the performance of the two proposed adaptive schemes. It
turns out that even though the use of the WLR-based SI may be advantageous in
some examples, this SI relies on an adaption constant, which may be hard to
tune: this affects the robustness of the resulting adaptive scheme. The use of
the MM-based SI, on the other hand, leads to a robust adaption strategy. Other
SIs may be tested and they may turn out to be even more robust and sharp, but
we leave this study for the future work.

## Acknowledgments

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 2-D LCD-Based Piecewise Linear Reconstruction

In this appendix, we describe how to reconstruct the one-sided point values
$\bm{U}^{\pm}_{{j+\frac{1}{2}},k}$ (the point values
$\bm{U}^{\pm}_{j,{k+\frac{1}{2}}}$ can be computed in a similar manner and we
omit the details for the sake of brevity). To this end, as in the 1-D case, we
first 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,\,j,\,j+1,\,j+2,$

where the matrix $R_{{j+\frac{1}{2}},k}$ is such that
$R^{-1}_{{j+\frac{1}{2}},k}\widehat{A}_{{j+\frac{1}{2}},k}R_{{j+\frac{1}{2}},k}$
is diagonal and a locally linearized Jacobian is
$\widehat{A}_{{j+\frac{1}{2}},k}:=A\big{(}(\,\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\big{)}$.

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}=\phi^{\rm
SBM}_{\theta,\tau}\left(\frac{\bm{\Gamma}_{j+1,k}-\bm{\Gamma}_{j,k}}{\bm{\Gamma}_{j,k}-\bm{\Gamma}_{j-1,k}}\right)\frac{\bm{\Gamma}_{j,k}-\bm{\Gamma}_{j-1,k}}{\Delta
x},$

and

$(\bm{\Gamma}_{x})_{j+1,k}=\phi^{\rm
SBM}_{\theta,\tau}\left(\frac{\bm{\Gamma}_{j+2,k}-\bm{\Gamma}_{j+1,k}}{\bm{\Gamma}_{j+1,k}-\bm{\Gamma}_{j,k}}\right)\frac{\bm{\Gamma}_{j+1,k}-\bm{\Gamma}_{j,k}}{\Delta
x},$

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

$\bm{\Gamma}^{-}_{{j+\frac{1}{2}},k}=\bm{\Gamma}_{j,k}+\frac{\Delta
x}{2}(\bm{\Gamma}_{x})_{j,k}\quad\mbox{and}\quad\bm{\Gamma}^{+}_{{j+\frac{1}{2}},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}^{\pm}_{j+\frac{1}{2}}=R_{{j+\frac{1}{2}},k}\bm{\Gamma}^{\pm}_{{j+\frac{1}{2}},k}.$

###### Remark A.1

The matrices $R_{{j+\frac{1}{2}},k}$ and $R^{-1}_{{j+\frac{1}{2}},k}$ for the
2-D Euler equation of gas dynamics (3.2)–(3.3) can be found in [3, Appendix
C].

## References

  * [1] F. Aràndiga, A. Baeza, and R. Donat, Vector cell-average multiresolution based on Hermite interpolation, Adv. Comput. Math., 28 (2008), pp. 1–22.
  * [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] A. Chertock, S. Chu, M. Herty, A. Kurganov, and M. Lukáčová-Medviďová, Local characteristic decomposition based central-upwind scheme, J. Comput. Phys., 473 (2023).  Paper No. 111718, 24 pp.
  * [4] J. Dewar, A. Kurganov, and M. Leopold, Pressure-based adaption indicator for compressible Euler equations, Numer. Methods Partial Differential Equations, 31 (2015), pp. 1844–1874.
  * [5] 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.
  * [6] W. S. Don, P. Li, K. Y. Wong, and Z. Gao, Improved symmetry property of high order weighted essentially non-oscillatory finite difference schemes for hyperbolic conservation laws, Adv. Appl. Math. Mech., 10 (2018), pp. 1418–1439.
  * [7] M. Dumbser, O. Zanotti, R. Loubère, and S. Diot, A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolic conservation laws, J. Comput. Phys., 278 (2014), pp. 47–75.
  * [8] N. Fleischmann, S. Adami, and N. A. Adams, Numerical symmetry-preserving techniques for low-dissipation shock-capturing schemes, Comput. & Fluids, 189 (2019), pp. 94–107.
  * [9] K. O. Friedrichs, Symmetric hyperbolic linear differential equations, Comm. Pure Appl. Math., 7 (1954), pp. 345–392.
  * [10] G. Fu and C.-W. Shu, A new troubled-cell indicator for discontinuous Galerkin methods for hyperbolic conservation laws, J. Comput. Phys., 347 (2017), pp. 305–327.
  * [11] A. Gelb and E. Tadmor, Spectral reconstruction of piecewise smooth functions from their discrete data, M2AN Math. Model. Numer. Anal., 36 (2002), pp. 155–175.
  * [12] , Adaptive edge detectors for piecewise smooth data based on the minmod limiter, J. Sci. Comput., 28 (2006), pp. 279–306.
  * [13] 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.
  * [14] 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.
  * [15] S. Gottlieb, C.-W. Shu, and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev., 43 (2001), pp. 89–112.
  * [16] J.-L. Guermond, R. Pasquetti, and B. Popov, Entropy viscosity method for nonlinear conservation laws, J. Comput. Phys., 230 (2011), pp. 4248–4267.
  * [17] A. Harten, ENO schemes with subcell resolution, J. Comput. Phys., 83 (1989), pp. 148–184.
  * [18] J. S. Hesthaven, Numerical methods for conservation laws: From analysis to algorithms, Comput. Sci. Eng. 18, SIAM, Philadelphia, 2018.
  * [19] E. Johnsen, On the treatment of contact discontinuities using WENO schemes, J. Comput. Phys., 230 (2011), pp. 8665–8668.
  * [20] S. Karni and A. Kurganov, Local error analysis for approximate solutions of hyperbolic conservation laws, Adv. Comput. Math., 22 (2005), pp. 79–99.
  * [21] S. Karni, A. Kurganov, and G. Petrova, A smoothness indicator for adaptive algorithms for hyperbolic systems, J. Comput. Phys., 178 (2002), pp. 323–341.
  * [22] 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.
  * [23] A. Kuganov and Y. Liu, New adaptive artificial viscosity method for hyperbolic systems of conservation laws, J. Comput. Phys., 231 (2012), pp. 8114–8132.
  * [24] 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.
  * [25] A. Kurganov and C.-T. Lin, On the reduction of numerical dissipation in central-upwind schemes, Commun. Comput. Phys., 2 (2007), pp. 141–163.
  * [26] A. Kurganov and E. Tadmor, Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers, Numer. Methods Partial Differential Equations, 18 (2002), pp. 584–608.
  * [27] A. Kurganov and R. Xin, New low-dissipation central-upwind schemes, J. Sci. Comput., 96 (2023).  Paper No. 56, 33 pp.
  * [28] P. D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical computation, Comm. Pure. Appl. Math., 7 (1954), pp. 159–193.
  * [29] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Appl. Math., Cambridge University Press, Cambridge, UK, 2002\.
  * [30] K.-A. Lie and S. Noelle, 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.
  * [31] 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.
  * [32] H. Nessyahu and E. Tadmor, Nonoscillatory central differencing for hyperbolic conservation laws, J. Comput. Phys., 87 (1990), pp. 408–463.
  * [33] G. Puppo and M. Semplice, Numerical entropy and adaptivity for finite volume schemes, Commun. Comput. Phys., 10 (2011), pp. 1132–1160.
  * [34] 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.
  * [35] , A comparison of troubled-cell indicators for Runge-Kutta discontinuous Galerkin methods using weighted essentially nonoscillatory limiters, SIAM J. Sci. Comput., 27 (2005), pp. 995–1013.
  * [36] C. W. Schulz-Rinne, Classifcation of the riemann problem for two-dimensional gas dynamics, SIAM J. Math. Anal., 24 (1993), pp. 76–88.
  * [37] 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–1414.
  * [38] J. Shi, Y.-T. Zhang, and C.-W. Shu, Resolution of high order WENO schemes for complicated fow structures, J. Comput. Phys., 186 (2003), pp. 690–696.
  * [39] C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes, Acta Numer., 5 (2020), pp. 701–762.
  * [40] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys., 77 (1988), pp. 439–471.
  * [41] C.-W. Shu and S. Osher, Efficient implementation of essentially nonoscillatory shock-capturing schemes. II, J. Comput. Phys., 83 (1989), pp. 32–78.
  * [42] P. K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal., 21 (1984), pp. 995–1011.
  * [43] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics: A practical introduction, Springer-Verlag, Berlin, Heidelberg, third ed., 2009.
  * [44] M. J. Vuik and J. K. Ryan, Automated parameters for troubled-cell indicators using outlier detection, SIAM J. Sci. Comput., 38 (2016), pp. A84–A104.
  * [45] H. Wakimura, S. Takagi, and F. Xiao, Symmetry-preserving enforcement of low-dissipation method based on boundary variation diminishing principle, Comput. & Fluids, 233 (2022).  Paper No. 105227, 17 pp.
  * [46] 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.
  * [47] W. Wang, C.-W. Shu, H. C. Yee, D. V. Kotov, and B. Sjögreen, High order finite difference methods with subcell resolution for stiff multispecies discontinuity capturing, Commun. Comput. Phys., 17 (2015), pp. 317–336.
  * [48] P. Woodward and P. Colella, The numerical solution of two-dimensional fluid flow with strong shocks, J. Comput. Phys., 54 (1988), pp. 115–173.
  * [49] 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.