A Quantitative Study of Social Organisation in Open Source Software
Communities

Marcelo Serrano Zanetti, Emre Sarigöl, Ingo Scholtes, Claudio Juan Tessone,
Frank Schweitzer

## References

Marcelo Serrano Zanetti  Emre Sarigöl  Ingo Scholtes  Claudio Juan Tessone
Frank Schweitzer  Chair of Systems Design, ETH Zurich,
Zurich, Switzerland<EMAIL_ADDRESS>
  * [1]
  * [2] http://www.sg.ethz.ch
  * [3]
  * [4]

## A Quantitative Study of Social Organisation in Open Source Software
Communities

Marcelo Serrano Zanetti  Emre Sarigöl  Ingo Scholtes  Claudio Juan Tessone
Frank Schweitzer  Chair of Systems Design, ETH Zurich,
Zurich, Switzerland<EMAIL_ADDRESS>
###### Abstract

The success of open source projects crucially depends on the voluntary
contributions of a sufficiently large community of users. Apart from the mere
size of the community, interesting questions arise when looking at the
_evolution of structural features_ of collaborations between community
members. In this article, we discuss several network analytic proxies that can
be used to quantify different aspects of the social organisation in social
collaboration networks. We particularly focus on measures that can be related
to the cohesiveness of the communities, the distribution of responsibilities
and the resilience against turnover of community members. We present a
comparative analysis on a large-scale dataset that covers the full history of
collaborations between users of $14$ major open source software communities.
Our analysis covers both aggregate and time-evolving measures and highlights
differences in the social organisation across communities. We argue that our
results are a promising step towards the definition of suitable, potentially
multi-dimensional, resilience and risk indicators for open source software
communities.

* [5]
* [6]
* [7]

##  1 Introduction

* [8] What are the most important social factors that lead to successful and sustainable open source software projects? According to _Linus’ Law_ \- which states that “given enough eyeballs, all bugs are shallow” [bazaar1999] \- the quality and success of open source software (OSS) critically depends on the existence of a sufficiently large community of developers who review, modify and improve the publicly available source code. Apart from development efforts, another important success factor is the existence of a stable community of users who report software defects, request and inspire new features, reproduce bugs or comment on issues reported by other users. By employing the collective knowledge and diverse experiences of many contributors, most OSS communities manage to provide technical assistance to less experienced users, often on a time scale that is competitive to commercial software support.
* [9] Depending on the distribution of competencies and responsibilities of contributors, largely different patterns of collaborations may arise. While it is generally difficult to assess these social factors of OSS projects, the availability of large scale data on community dynamics increasingly allows to study the _social dimension of OSS projects_ from a quantitative perspective [roblesturnover2006, mockuslongtermcontributions2012]. Previous studies have mainly focused on rather simple proxies of social dynamics like the evolution of the number of contributors and contributions or the time span of a user’s activity and were mostly based on a rather limited set of snapshots of a single project. Using a large scale dataset of time-stamped social interactions that has been collected from the Bugzilla bug-tracker installations of $14$ major OSS projects, in this paper we study the _fine-grained evolution of structural features of networks of user collaborations_. We thus take a _network perspective on OSS communities_ and highlight differences in the social organisation of software projects that can be related to their activity, their cohesion as well as their resilience against fluctuations in the community members. By applying standard measures from social network analysis we particularly quantify how tightly community members collaborate, how equal responsibilities are distributed and how resilient collaboration topologies are against the loss of (central) community members. While similar tools have been applied to OSS projects before [howisonsocial2006, niavaliditynet2010], to the best of our knowledge, the present paper is the first to study these network analytic measures on a dataset that covers the full, fine-grained history of $14$ well-established and successful OSS communities.
* [10]
* [11]

##  2 Social Organisation in OSS Communities: A Network Perspective

* [12] In order to make substantiated statements about the structure and dynamics of the social organisation of OSS communities, we recently completed collecting data on the history of user collaborations recorded by the Bugzilla installation of $14$ well-established OSS projects. Bugzilla[serranobugzilla2005] is an open source bug tracking system which is utilised by users and developers alike to report bugs, keep track of open issues and feature requests and comment on issues reported by others. Since the Bugzilla installations of OSS projects are used to foster collaboration between community members, it constitutes a valuable source of data that allows us to track social interactions between developers and users.
* [13]
* [14]

###  2.1 Building Social Networks from Bug-Reports

* [15] Data in the Bugzilla database are arranged around the notion of _bug reports_. Each bug report has a set of fields describing aspects like the user who initially filed the bug report, its current status (e.g. _pending_ , _reproduced_ , _solved_ , etc), to whom the responsibility to provide a fix has been assigned, attachments which may be used to reproduce or resolve the issue, comments and hints by other community members, or a list of community members which shall be informed about future updates. Apart from an initial bug report, Bugzilla additionally stores the full history of all updates to any of the fields of a bug report. Each of these change records includes a time stamp, the ID of the user performing the change as well as the new values of the changed fields. While our dataset comprises change records for all possible fields, in this article we focus on those that indicate changes in the users that are assigned responsibility to fix an issue (henceforth called the _ASSIGNEE_ field) and changes to the list of users to whom future updates of the bug shall be sent via E-Mail (henceforth called the _CC_ field). We consider any updates in the CC and ASSIGNEE field of a bug report as a time-stamped edge from the user who performed the update to the user(s) who were added to the CC field or the ASSIGNEE list of responsible developers respectively.
* [16] Based on the data extraction procedure described above, we obtain a large time-aggregated network of nodes representing community members and time-stamped edges representing a particular interaction between two users. For most of the projects considered, the Bugzilla history from which we extract the network is longer than ten years. The fact that - in social networks aggregated over such long periods of time - most of the users represented by nodes have never been active within the same time period limits the expressiveness of the network structure in terms of a project’s “social organisation”. In order to overcome this issue, we perform a _dynamic network analysis_ by defining a sequence of _monthly collaboration networks_ based on the time stamps of edges. In particular, we define a 30 day sliding time window and filter out those edges whose time stamps are outside the window and those nodes who did not have any interactions in the corresponding time period. By progressively advancing the start date of the sliding 30 day time window by one day increments we obtain a sequence of collaboration networks that allows us to study the structure of the community’s social organisation as well as its evolution over time. Naturally, most of the monthly networks obtained in the way described above will not be fully connected. Since the network analytic measures we intend to apply assume connected topologies, we perform a component analysis on all snapshots and restrict our quantitative analysis to the largest connected component (LCC). In order to test the significance of our findings we further compute the fraction of those nodes who are part of the largest connected component. Table 1 shows the $14$ OSS projects that are included in our dataset along with the time period and the total number of bug reports and updates that we included in our analysis. The column _LCC/TOTAL_ furthermore indicates the fraction of users in the LCC, averaged over all monthly snapshots of the corresponding project. Here one observes that our data shows a rather large degree of variation with respect to this fraction, which may be seen as an argument that this measure is an interesting indicator for the _cohesiveness_ of OSS communities by itself. Nevertheless, we argue that for all projects the fraction of users in the LCC is sufficiently large to make substantiated statements about the project’s social organisation.
* [17]
Table 1: Aggregated measures for the studied projects. From column _LCC/Total_
to the last on the right, the numbers indicate the mean value $\pm$ standard
deviation.

Project | Bugs | Updates | Period | LCC/Total | Nodes in LCC | Edges | Mean | Assortativity | Closeness | Clustering
---|---|---|---|---|---|---|---|---|---|---
Name | Degree | Central. | Coefficient
xamarin | 4552 | 20721 | 2011-2012 | 0.93$\pm$0.05 | 46.76$\pm$8.12 | 98.15$\pm$22.70 | 2.07$\pm$0.29 | -0.14$\pm$0.11 | 0.40$\pm$0.07 | 0.22$\pm$0.05
thunderbird | 35388 | 313957 | 2000-2012 | 0.53$\pm$0.26 | 64.82$\pm$53.49 | 86.44$\pm$80.05 | 1.05$\pm$0.42 | -0.23$\pm$0.17 | 0.40$\pm$0.27 | 0.04$\pm$0.05
libreoffice | 8916 | 78341 | 2010-2012 | 0.78$\pm$0.11 | 73.83$\pm$32.06 | 114.41$\pm$49.10 | 1.56$\pm$0.26 | -0.20$\pm$0.10 | 0.40$\pm$0.09 | 0.13$\pm$0.06
mageia | 6600 | 46921 | 2006-2012 | 0.93$\pm$0.07 | 77.54$\pm$21.80 | 156.00$\pm$59.24 | 1.95$\pm$0.30 | -0.37$\pm$0.12 | 0.54$\pm$0.09 | 0.14$\pm$0.04
mandriva | 60546 | 368463 | 2002-2012 | 0.70$\pm$0.18 | 88.15$\pm$60.70 | 142.16$\pm$118.44 | 1.41$\pm$0.38 | -0.29$\pm$0.15 | 0.40$\pm$0.14 | 0.07$\pm$0.05
firefox | 112953 | 1067914 | 1999-2012 | 0.58$\pm$0.23 | 171.77$\pm$117.79 | 240.79$\pm$180.44 | 1.16$\pm$0.44 | -0.15$\pm$0.11 | 0.32$\pm$0.23 | 0.04$\pm$0.04
seamonkey | 90040 | 993392 | 1998-2012 | 0.67$\pm$0.15 | 210.39$\pm$251.43 | 364.42$\pm$482.54 | 1.48$\pm$0.48 | -0.19$\pm$0.13 | 0.34$\pm$0.11 | 0.08$\pm$0.06
netbeans | 210921 | 1875878 | 2000-2012 | 0.96$\pm$0.05 | 269.71$\pm$292.07 | 1069.72$\pm$1509.12 | 3.39$\pm$1.13 | -0.12$\pm$0.08 | 0.37$\pm$0.05 | 0.23$\pm$0.08
openoffice | 118135 | 915749 | 2000-2012 | 0.88$\pm$0.19 | 319.01$\pm$169.88 | 931.35$\pm$591.80 | 2.52$\pm$0.84 | -0.12$\pm$0.10 | 0.34$\pm$0.15 | 0.12$\pm$0.06
gentoo | 140216 | 661783 | 2002-2012 | 0.80$\pm$0.07 | 338.97$\pm$110.86 | 617.73$\pm$211.92 | 1.82$\pm$0.27 | -0.29$\pm$0.10 | 0.49$\pm$0.13 | 0.04$\pm$0.03
kde | 179470 | 648331 | 2002-2012 | 0.75$\pm$0.12 | 361.16$\pm$246.16 | 424.61$\pm$301.20 | 1.15$\pm$0.07 | -0.16$\pm$0.07 | 0.32$\pm$0.07 | 0.01$\pm$0.01
eclipse | 356415 | 2594385 | 2001-2012 | 0.78$\pm$0.08 | 472.58$\pm$180.71 | 964.47$\pm$411.94 | 2.06$\pm$0.38 | 0.05$\pm$0.08 | 0.25$\pm$0.05 | 0.13$\pm$0.03
gnome | 550722 | 2751441 | 2000-2012 | 0.67$\pm$0.12 | 523.76$\pm$585.26 | 610.16$\pm$616.81 | 1.25$\pm$0.22 | -0.17$\pm$0.09 | 0.25$\pm$0.08 | 0.03$\pm$0.04
redhat | 414163 | 3777634 | 2006-2012 | 0.45$\pm$0.26 | 658.06$\pm$865.97 | 983.58$\pm$1297.18 | 1.19$\pm$0.35 | -0.12$\pm$0.20 | 0.30$\pm$0.23 | 0.00$\pm$0.01

* [18]
* [19]

###  2.2 Network Measures

While the literature is rich in terms of measures able to quantify structural
features of networks [stanleysocial1994, newmannetsbook2010], due to space
limitations here we focus on three measures which are able to capture basic
network qualities that relate to the _cohesiveness_ of a community, the
distribution of responsibilities among its members and its resilience against
fluctuations in the user base. The first network measure is based on the
_closeness centrality_ of a node, which is defined as the inverse of the sum
of the shortest path length to all other nodes in the network.

$Cc(n_{i})=\sum_{j=1,j\neq i}^{N}\frac{N-1}{d(n_{i},n_{j})}\in[0,1]$ (1)

where $Cc(n_{i})$ corresponds to the _closeness centrality_ score of node
$n_{i}$, $d(n_{i},n_{j})$ is the length of the shortest path between nodes
$n_{i}$ and $n_{j}$, while $N$ corresponds to the total number of nodes in a
given network. Finally, the factor $N-1$ is a normalisation constant
[freemancentrality1979]. Based on this, the _closeness centralisation_ of a
network $(Cc_{global})$ can be calculated by taking the sum of the differences
between the node with the highest value of closeness centrality $(n^{*})$ and
the closeness centrality scores of all other nodes. This quantity is then
normalised to the range of $0$ to $1$ using the theoretical value that results
from a (maximally centralised) star network. Equation (2) presents the formal
definition, while more details can be found in [freemancentrality1979,
stanleysocial1994]. In the context of OSS collaboration networks, closeness
centralisation captures to what degree responsibilities, collaboration and
communication are distributed equally across community members.

* [20] $Cc_{global}=\sum_{i=1}^{N}\frac{Cc(n^{*})-Cc(n_{i})}{\frac{(N-2)(N-1)}{2N-3}}\in[0,1]$ (2)
* [21] The second measure, the _clustering coefficient_ of a network $(C)$, measures how closely community members interact with each other in the sense that an interaction between a user $X$ and $Y$, as well as an interaction between user $Y$ and $Z$ will also entail a direct interaction between the users $Y$ and $Z$. The formal definition is presented in equations (3) and (4).
* [22] $C(n_{i})=\frac{2L_{D_{n_{i}}}}{D_{n_{i}}(D_{n_{i}}-1)}\in[0,1]$ (3) $C=\frac{1}{N}\sum_{i=1}^{N}C(n_{i})\in[0,1]$ (4) where $D_{n_{i}}$ is the number of nodes directly connected to the node $n_{i}$, while $L_{D_{n_{i}}}$ is the number of edges between them. Therefore, the clustering coefficient $C(n_{i})$ of node $n_{i}$ expresses the fraction of edges that were realised from the possible $\frac{D_{n_{i}}(D_{n_{i}}-1)}{2}$ edges which are expected in a fully connected network with $D_{n_{i}}$ nodes. We obtain the clustering coefficient of a network by averaging the clustering coefficient scores of all existing nodes (see equation (4)). This procedure can be seen as measuring how _cohesive_ the community is in terms of nodes being embedded in collaborating clusters [stanleysocial1994].
* [23] Finally, the _assortativity_ $(r)$ measures an individual’s preference to connect to other individuals that have a similar or different degree of connectivity (the degree being a node’s number of connections to different nodes). Networks in which nodes are preferentially connected to nodes with similar degree are called assortative. In this case a positive degree assortativity $(0\ll r\leq 1)$ indicates a positive correlation between the degrees of neighbouring nodes. Networks in which nodes are preferentially connected to nodes with different degree are called disassortative and in this case degree assortativity is negative $(0\gg r\geq-1)$. In networks with zero degree assortativity, there is no correlation between the degrees of connected nodes, i.e. nodes do not exhibit a preference for one or the other. Formally,  $r=\frac{\sum_{ij}ij(e_{i,j}-q(i)q(j))}{\sigma(q)^{2}}\in[-1,1]$ (5) where $e_{ij}$ is the fraction of all links in the network that join together nodes with degrees $i$ and $j$, $q(i)=\sum_{j}e_{i,j}$, $q(j)=\sum_{i}e_{i,j}$ and $\sigma(q)$ is the standard deviation of the distribution of $q$. The term $q(i)q(j)$ is the equivalent to the expected value of $e_{i,j}$ inferred from a random network. Therefore, if $r=0$ the pattern of interconnection between nodes is also random [newmanR2003].
* [24]
* [25]

##  3 Comparative Analysis of OSS Communities

* [26] As described above, the preliminary results presented here have been obtained for the LCC of the network of monthly collaborations in terms of CC and ASSIGNEE interactions. While Table 1 shows the aggregate measures averaged over all time windows for every project in our database, due to space constraints we limit the presentation of the dynamics of the social organisation to the projects Gentoo and KDE (both Gnu/Linux related projects) as well as Eclipse and NetBeans (both Java IDEs). These have been chosen because a) their communities are of comparable size and age, b) the respective pairs of projects address similar problem domains and c) they represent contrasting examples with respect to the measures studied in this paper.
* [27] Figure 1 shows the evolution of the number of nodes in the LCC, its assortativity, clustering coefficient and closeness centralisation for these four projects. For all projects, the fraction of nodes in the LCC is rather stable with values between $0.7$ and $1$ consistent with the aggregate values given in Table 1. The same is true for the evolution of the mean degree. We thus omit these plots. The four projects show significant differences in the evolution of the clustering coefficient that cannot be explained by mere size effects. In the particular time frame between $2006$ and $2008$, the clustering coefficient of the Eclipse community ($\approx 0.15$) was roughly ten times higher than that of the Gentoo community ($\approx 0.01$), although the LCCs of both communities were of comparable size ($\approx 500$ nodes). In addition, the clustering coefficient of the Gentoo community shows an interesting dynamics, dropping to a very small value between $2006$ and $2008$ and increasing thereafter.
* [28]

(a) Nodes in LCC (b) Assortativity (c) Clustering coeff. (d) Closeness
central.

(e) Nodes in LCC (f) Assortativity (g) Clustering coeff. (h) Closeness
central.

Figure 1: Evolution of structural measures of the LCC in the monthly Bugzilla
collaboration networks. (a-d): Gnu/Linux related projects Gentoo (black) and
KDE (green), (e-h): IDEs Eclipse (black) and NetBeans (green).

* [29]

(a) Gentoo (Jan/2006)
$nodes=535,edges=785$ (b) KDE (Feb/2011)
$nodes=543,edges=630$ (c) Eclipse (Jan/2010)
$nodes=502,edges=868$ (d) Netbeans (Sep/2008)
$nodes=566,edges=2753$

Figure 2: Four monthly collaboration networks with comparable size showing
largely different social organisation (the visualisation was generated by
Gephi [gephi]).

* [30] A different perspective of the structural change the Gentoo community was undergoing is given in Figure 1(d) which displays a visible plateau in the closeness centralisation of the network within the same period. In fact, as can be seen in the network depicted in Figure 2(a), in the period between $2006$ and $2008$ most of the collaborations were mediated by a single central community member, while the social organisation of the Eclipse community depicted in 2(c) was structured in a much more homogeneous way. The evolution of degree assortativity is captured in Figures 1(b) and 1(f). Both the level of degree assortativity as well as its dynamics differ across the projects. The collaboration network of Eclipse exhibits a tendency towards assortative structures (meaning that high degree nodes are preferentially connected to high degree nodes). The opposite is true for the KDE and the Gentoo communities which show a tendency towards disassortativity. We thus argue that assortativity is suitable to further differentiate the social organisation of OSS communities.
* [31]
* [32]

##  4 Conclusions and Future Work

We have studied measures that capture different structural dimensions in the
social organisation of OSS projects. Our analysis is based on a comprehensive
dataset collected from the bug tracking communities of $14$ major OSS
projects. We view the social organisation from the perspective of time-
evolving networks and highlight how projects, although similar in terms of
size, problem domain and age, a) largely differ in terms of clustering
coefficient, assortativity and closeness centralisation and b) that some
projects show interesting dynamics with respect to these measures that cannot
be explained by mere size effects. We argue that the phase of high closeness
centralisation and low clustering coefficient observed in the Gentoo community
between $2006$ and $2008$ may be interpreted as a lack of _social cohesion_
which can possibly pose a risk for the project.

* [33] While our results are necessarily preliminary, we currently extend our work by adding spectral measures like algebraic connectivity and inequality measures like the Gini index that can highlight further differences in the social organisation [zanettiACS2013]. A detailed case study is under preparation [zanettiICSE2013] and further includes community performance indicators (e.g. response times, bug fixing times and fraction of open issues) that can be mined from our dataset. The eventual goal of our project is the provision of multi-dimensional indicators for the social and technical organisation of OSS projects that are correlated with performance and that can be considered in the management and evaluation of OSS projects [zanettiDOCSYMPICSE2012, zanettiQmodularity12, sustainablegrowthtessone2011]. Such indicators can be useful when taking informed decisions about which OSS project to invest in or rely on. Furthermore, due to the distributed nature of collaborations, individuals often lack a global perspective on evolving communication and coordination structures, even though these can influence long-term success. An inclusion of suitable indicators in community platforms like e.g. Bugzilla can assist in determining risks and allow project managers to timely react by shifting responsibilities, fostering information flow or changing organisational procedures.
* [34]
* [35]

##  5 Acknowledgement

* [36] This work was supported by the SNF through grant CR12I1_125298. C.J.T. acknowledges financial support from the SBF through grant C09.0055.
* [37]
* [38]


# New Low-Dissipation Central-Upwind Schemes. Part II

Shaoshuai Chu Department of Mathematics and Shenzhen International Center for
Mathematics, Southern University of Science and Technology, Shenzhen, 518055,
China<EMAIL_ADDRESS>Alexander Kurganov Department 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>and Ruixiao Xin Department of Mathematics, Southern University of Science and
Technology, Shenzhen, 518055, China<EMAIL_ADDRESS>
###### Abstract

The low-dissipation central-upwind (LDCU) schemes have been recently
introduced in [A. Kurganov and R. Xin, J. Sci. Comput., 96 (2023), Paper No.
56] as a modification of the central-upwind (CU) schemes from [A. Kurganov and
C. T. Lin, Commun. Comput. Phys., 2 (2007), pp. 141-163]. The LDCU schemes
achieve much higher resolution of contact waves and many (two-dimensional)
structures resulting from complicated wave interaction. However, the LDCU
schemes sometimes produce more oscillatory results compared with the CU
schemes, especially near the computational domain boundaries.

In this paper, we propose a very simple—yet systematic—modification of the
LDCU schemes, which completely eliminates the aforementioned oscillations
almost without affecting the quality of the computed solution.

Key words: Hyperbolic systems of conservation laws, low-dissipation central-
upwind schemes, Euler equations of gas dynamics.

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

##  1 Introduction

We consider hyperbolic systems of conservation laws, which in the one-
dimensional (1-D) case read as

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

and 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.2)

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

We focus on finite-volume central-upwind (CU) schemes, which were introduced
in [4, 6, 7] as a “black-box” solver for general hyperbolic systems (1.1) and
(1.2). Even though the CU schemes are quite accurate, efficient, and robust
for a wide variety of hyperbolic systems, higher resolution of the numerical
solutions can be achieved by further reducing numerical dissipation present in
the CU schemes; see, e.g., [5, 1, 8]. The low-dissipation CU (LDCU) schemes
have been recently proposed in [8] as a modification of the CU scheme from
[5]. The LDCU schemes achieve much higher resolution of contact waves and many
(two-dimensional) structures resulting from complicated wave interaction.
However, the LDCU schemes sometimes produce more oscillatory results compared
with the CU schemes, especially near the computational domain boundaries.

In order to suppress these spurious oscillations, we modify both 1-D and 2-D
LDCU schemes. The proposed modifications are very simple, yet systematic as
they are based on a more accurate projection of the evolved solution onto the
original (uniform) mesh. The new LDCU schemes are developed for the 1-D and
2-D Euler equations of gas dynamics and tested on several numerical examples.
The obtained numerical results demonstrate that the new schemes contain almost
the same small amount of numerical dissipation as the LDCU schemes from [8]
but produce substantially “cleaner”, non-oscillatory computed solutions.

##  2 New LDCU Schemes

We follow the derivation of the LDCU scheme in [8] and use precisely the same
notation as in [8].

We cover the computational domain with the finite volume cells
$C_{j}=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$, $j=1,\ldots,N$, which are
assumed to be uniform, that is,
$x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}\equiv\Delta x$. We then assume that the
solution, realized in terms of its cell averages $\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,n}_{j}$, is available at a certain time level $t=t^{n}$ and
reconstruct a second-order piecewise linear interpolant
$\sum_{j}[\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\,n}_{j}+(\bm{U}_{x})^{n}_{j}(x-x_{j})]{\cal X}_{C_{j}}$, where
$\cal X$ denotes the characteristic function of the corresponding intervals
and $(\bm{U}_{x})^{n}_{j}$ are the slopes obtained using a nonlinear limiter.
In the numerical experiments reported in §3, we have used a generalized minmod
limiter with the minmod parameter $\theta=1.3$; see, e.g., [9, 11, 15, 16].

We then evaluate the local speeds of propagation $a^{\pm}_{j+\frac{1}{2}}$,
introduce the corresponding points
$x^{n}_{{j+\frac{1}{2}},\ell}:=x_{j+\frac{1}{2}}+a^{-}_{j+\frac{1}{2}}\Delta
t^{n}$ and
$x^{n}_{{j+\frac{1}{2}},r}:=x_{j+\frac{1}{2}}+a^{+}_{j+\frac{1}{2}}\Delta
t^{n}$, and integrate the system (1.1) over the space-time control volumes
consisting of the “smooth”
$[x_{{j-\frac{1}{2}},r},x_{{j+\frac{1}{2}},\ell}]\times[t^{n},t^{n+1}]$ and
“nonsmooth”
$[x_{{j+\frac{1}{2}},\ell},x_{{j+\frac{1}{2}},r}]\times[t^{n},t^{n+1}]$ areas,
where $t^{n+1}:=t^{n}+\Delta t^{n}$. This way the solution is evolved in time
and upon the completion of the evolution step, we obtain the intermediate cell
averages $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\rm\,int}_{j+\frac{1}{2}}$ (see [8, Eq. (2.7)]) and
$\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\rm\,int}_{j}$ (see [8,
Eq. (2.10)]).

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

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

In order to develop the original LDCU scheme in [8], we have used

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

where $\widetilde{\bm{U}}^{\,\rm int}_{j+\frac{1}{2}}(x)$ is discontinuous at
$x=x_{j+\frac{1}{2}}$. The values $\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$ and $\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$ are determined based on the conservation requirement

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

and $d$ additional degrees of freedom, which can be used to design an accurate
projection. The way these degrees freedom are utilized depends on the problem
at hand. One, however, can introduce an additional degree of freedom, which
may be used to further improve the accuracy of the projection step.

###  2.1 Modification of the Projection Step

We propose to replace (2.2) with

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

where
$\widetilde{x}_{j+\frac{1}{2}}=x_{j+\frac{1}{2}}+\widetilde{u}_{j+\frac{1}{2}}\Delta
t^{n}$ and
$\widetilde{u}_{j+\frac{1}{2}}\in(a^{-}_{j+\frac{1}{2}},a^{+}_{j+\frac{1}{2}})$
represents an additional degree of freedom. The conservation requirements now
give

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

compare with (2.3), which is obtained from (2.5) if
$\widetilde{u}_{j+\frac{1}{2}}=0$. The modified projection step is outlined in
Figure 2.1; compare it with [8, Fig. 3].

Figure 2.1: The new projection step.

In order to complete the derivation of the new LDCU scheme, we need to
consider a particular hyperbolic system. As in [8], we study the Euler
equation of gas dynamics.

###  2.2 New LDCU Scheme for the 1-D Euler Equations of Gas Dynamics

The 1-D Euler equations of gas dynamics reads as (1.1) with $\bm{U}=(\rho,\rho
u,E)^{\top}$ and $\bm{F}(\bm{U})=(\rho u,\rho u^{2}+p,u(E+p))^{\top}$, where
$\rho$ is the density, $u$ is the velocity, $p$ is the pressure, and $E$ is
the total energy. The system is closed using the equation of states (EOS),
which is in the case of ideal gas is $p=(\gamma-1)\big{[}E-\frac{1}{2}\rho
u^{2}\big{]}$, $\gamma={\rm Const}$.

In order to derive the new LDCU scheme for the 1-D Euler equations, we first
follow [8, §3.1] to obtain the point values $\rho^{\pm}_{j+\frac{1}{2}}$,
$(\rho u)^{\pm}_{j+\frac{1}{2}}$, and $E^{\pm}_{j+\frac{1}{2}}$, evaluate the
corresponding point values $u^{\pm}_{j+\frac{1}{2}}$ and
$p^{\pm}_{j+\frac{1}{2}}$, and estimate the one-sided local speeds of
propagation $a^{\pm}_{j+\frac{1}{2}}$.

We then proceed with the evolution of the subcell averages
$\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}=\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}},(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\,\rm
int,L}_{j+\frac{1}{2}},\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}\big{)}^{\top}$ and $\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}=\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}},(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\,\rm
int,R}_{j+\frac{1}{2}},\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}\big{)}^{\top}$ required in (2.4). Similarly to [8, §3],
we enforce the continuity of $u$ and $p$ across the cell interfaces by setting

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

Next, (2.6) together with the conservation requirement (2.5) applied to
$\rho$, $\rho u$, and $E$,

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

yield a system of five algebraic equations (2.6)–(2.7), which we solve for
$(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho
u$\kern 0.0pt}}})^{\,\rm int,L}_{j+\frac{1}{2}}$, $(\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\,\rm
int,R}_{j+\frac{1}{2}}$, $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$, and $\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$E$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$, and express these quantities in terms of
$\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$ and $\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$:

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

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

Note that the first two equations in (2.8) give $u^{\rm
int,L}_{j+\frac{1}{2}}=u^{\rm int}_{j+\frac{1}{2}}=u^{\rm
int,R}_{j+\frac{1}{2}}$, which suggests that in this piecewise constant
solution approximation, the velocity is constant in
$[x_{{j+\frac{1}{2}},\ell},x_{{j+\frac{1}{2}},r}]$. We therefore use the same
velocity value and set $\widetilde{u}_{j+\frac{1}{2}}=u^{\rm
int}_{j+\frac{1}{2}}$.

Finally, we follow the steps in [8, §3], where we make the difference
$\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$ as large as possible without creating any new local
extrema. To this end, we denote by

$S^{-}_{j+\frac{1}{2}}:=(u^{\rm
int}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}})\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\rho^{\rm
int}_{{j+\frac{1}{2}},\ell}\big{)}\quad\mbox{and}\quad
S^{+}_{j+\frac{1}{2}}:=(a^{+}_{j+\frac{1}{2}}-u^{\rm
int}_{j+\frac{1}{2}})\big{(}\rho^{\rm
int}_{{j+\frac{1}{2}},r}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}\big{)},$

where the point values $\rho^{\rm int}_{{j+\frac{1}{2}},\ell}$ and $\rho^{\rm
int}_{{j+\frac{1}{2}},\ell}$ were introduced in [8, Eq. (2.14)], and then
determine $\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}$ and $\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}$ by

$\begin{aligned} &\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+\frac{\delta_{j+\frac{1}{2}}}{a^{-}_{j+\frac{1}{2}}-u^{\rm
int}_{j+\frac{1}{2}}},\\\ &\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}}=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}+\frac{\delta_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}-u^{\rm
int}_{j+\frac{1}{2}}},\end{aligned}\qquad\delta_{j+\frac{1}{2}}:={\rm
minmod}\big{(}S^{-}_{j+\frac{1}{2}},S^{+}_{j+\frac{1}{2}}\big{)};$ (2.9)

compare these formulae with [8, Eq. (3.14)]. Here, ${\rm
minmod}(a,b):=\frac{1}{2}({\rm sgn}(a)+{\rm sgn}(b))\min(|a|,|b|)$.

We then substitute (2.9) into (2.8) and obtain

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

####  2.2.1 Fully Discrete Scheme

We now derive a new fully discrete LDCU scheme based on the new projection
step. To this end, we integrate the piecewise constant interpolant (2.1),
(2.4) over the cell $C_{j}$ and obtain

$\displaystyle\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,n+1}_{j}$
$\displaystyle=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j}+\frac{\Delta t^{n}}{\Delta
x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{j-\frac{1}{2}}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm int}_{j}\big{)}$
(2.11) $\displaystyle+\max(u^{\rm int}_{j-\frac{1}{2}},0)(\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{j-\frac{1}{2}}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{j-\frac{1}{2}})-\max(u^{\rm
int}_{j+\frac{1}{2}},0)(\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,L}_{j+\frac{1}{2}}-\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int,R}_{j+\frac{1}{2}})\Big{]}$
$\displaystyle\stackrel{{\scriptstyle(\ref{3.14}),(\ref{3.10f})}}{{=}}\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}^{\rm\,int}_{j}+\frac{\Delta t^{n}}{\Delta
x}\Big{[}a^{+}_{j-\frac{1}{2}}\big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j-\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}-a^{-}_{j+\frac{1}{2}}\big{(}\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j+\frac{1}{2}}-\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}^{\,\rm
int}_{j}\big{)}+\alpha^{\rm
int}_{j-\frac{1}{2}}{\bm{\delta}}_{j-\frac{1}{2}}-\alpha^{\rm
int}_{j+\frac{1}{2}}{\bm{\delta}}_{j+\frac{1}{2}}\Big{]},$

where

$\alpha^{\rm int}_{j+\frac{1}{2}}=\left\\{\begin{aligned}
&\frac{a^{+}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}-u^{\rm
int}_{j+\frac{1}{2}}}&&\mbox{if}~{}u^{\rm int}_{j+\frac{1}{2}}<0,\\\
&\frac{a^{-}_{j+\frac{1}{2}}}{a^{-}_{j+\frac{1}{2}}-u^{\rm
int}_{j+\frac{1}{2}}}&&\mbox{otherwise},\end{aligned}\right.\qquad\bm{\delta}_{j+\frac{1}{2}}=\delta_{j+\frac{1}{2}}\begin{pmatrix}1\\\
u^{\rm int}_{j+\frac{1}{2}}\\\\[3.44444pt] \dfrac{1}{2}\big{(}u^{\rm
int}_{j+\frac{1}{2}}\big{)}^{2}\end{pmatrix}.$ (2.12)

####  2.2.2 Semi-Discrete Scheme

We now pass to the semi-discrete limit $\Delta t^{n}\to 0$ in (2.11)–(2.12)
and proceed as in [8, §3.1.2] to end up with the following semi-discretization

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

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

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

and

$\bm{q}_{j+\frac{1}{2}}=\alpha^{*}_{j+\frac{1}{2}}q^{\rho}_{j+\frac{1}{2}}\begin{pmatrix}1\\\
u^{*}_{j+\frac{1}{2}}\\\\[3.44444pt]
\dfrac{1}{2}\big{(}u^{*}_{j+\frac{1}{2}}\big{)}^{2}\end{pmatrix}$ (2.15)

is a modified “built-in” anti-diffusion term with

$\displaystyle\rho^{*}_{j+\frac{1}{2}}=\frac{a^{+}_{j+\frac{1}{2}}\rho^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\rho^{-}_{j+\frac{1}{2}}-\left[(\rho
u)^{+}_{j+\frac{1}{2}}-(\rho
u)^{-}_{j+\frac{1}{2}}\right]}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}},$
(2.16) $\displaystyle(\rho
u)^{*}_{j+\frac{1}{2}}=\frac{a^{+}_{j+\frac{1}{2}}(\rho
u)^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}(\rho
u)^{-}_{j+\frac{1}{2}}-\left[\rho^{+}_{j+\frac{1}{2}}\big{(}u^{+}_{j+\frac{1}{2}}\big{)}^{2}+p^{+}_{j+\frac{1}{2}}-\rho^{-}_{j+\frac{1}{2}}\big{(}u^{-}_{j+\frac{1}{2}}\big{)}^{2}-p^{-}_{j+\frac{1}{2}}\right]}{a^{+}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}},$
$\displaystyle u^{*}_{j+\frac{1}{2}}=\frac{(\rho
u)^{*}_{j+\frac{1}{2}}}{\rho^{*}_{j+\frac{1}{2}}},\quad\alpha^{*}_{j+\frac{1}{2}}=\left\\{\begin{aligned}
&\frac{a^{+}_{j+\frac{1}{2}}}{a^{+}_{j+\frac{1}{2}}-u^{*}_{j+\frac{1}{2}}}&&\mbox{if}~{}u^{*}_{j+\frac{1}{2}}<0,\\\
&\frac{a^{-}_{j+\frac{1}{2}}}{a^{-}_{j+\frac{1}{2}}-u^{*}_{j+\frac{1}{2}}}&&\mbox{otherwise},\end{aligned}\right.$
$\displaystyle q^{\rho}_{j+\frac{1}{2}}={\rm
minmod}\Big{(}\big{(}u^{*}_{j+\frac{1}{2}}-a^{-}_{j+\frac{1}{2}}\big{)}\big{(}\rho^{*}_{j+\frac{1}{2}}-\rho^{-}_{j+\frac{1}{2}}\big{)},\big{(}a^{+}_{j+\frac{1}{2}}-u^{*}_{j+\frac{1}{2}}\big{)}\big{(}\rho^{+}_{j+\frac{1}{2}}-\rho^{*}_{j+\frac{1}{2}}\big{)}\Big{)}.$

Note that all of the indexed quantities in (2.14)–(2.16) are time dependent,
but from now on we will omit this dependence for the sake of brevity.

###### Remark 2.1

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

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

###  2.3 New LDCU Scheme for the 2-D Euler Equations of Gas Dynamics

In this section, we extend the modified semi-discrete LDCU scheme from §2.2.2
to the 2-D Euler equations of gas dynamics, which read as (1.2) with
$\bm{U}=(\rho,\rho u,\rho v,E)^{\top}$, $\bm{F}(\bm{U})=(\rho u,\rho
u^{2}+p,\rho uv,u(E+p))^{\top}$, and $\bm{G}(\bm{U})=(\rho v,\rho uv,\rho
v^{2}+p,v(E+p))^{\top}$, where the notations are the same as in the 1-D case
except for that now $u$ and $v$ are the $x$\- and $y$-velocities,
respectively. The system is closed using the EOS for the ideal gas
$\,p=(\gamma-1)\big{[}E-\frac{\rho}{2}(u^{2}+v^{2})\big{]}$.

We first introduce a uniform mesh consisting of the finite-volume cells
$C_{j,k}:=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]\times[y_{k-\frac{1}{2}},y_{k+\frac{1}{2}}]$
of the uniform size 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$, $j=1,\ldots,N_{x}$,
$k=1,\ldots,N_{y}$. We assume that at certain time level $t$, an approximate
solution, realized in terms of the cell averages $\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}$,
is available. These cell averages are then evolved in time by solving the
following system of ODEs:

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

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

$\displaystyle\begin{aligned} \bm{{\cal
F}}_{{j+\frac{1}{2}},k}&=\frac{a^{+}_{{j+\frac{1}{2}},k}\bm{F}(\bm{U}^{-}_{{j+\frac{1}{2}},k})-a^{-}_{{j+\frac{1}{2}},k}\bm{F}(\bm{U}^{+}_{{j+\frac{1}{2}},k})}{a^{+}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}}+\frac{a^{+}_{{j+\frac{1}{2}},k}a^{-}_{{j+\frac{1}{2}},k}}{a^{+}_{{j+\frac{1}{2}},k}-a^{-}_{{j+\frac{1}{2}},k}}\left(\bm{U}^{+}_{{j+\frac{1}{2}},k}-\bm{U}^{-}_{{j+\frac{1}{2}},k}\right)\\\
&+\bm{q}_{{j+\frac{1}{2}},k},\end{aligned}$ (2.18)
$\displaystyle\begin{aligned} \bm{{\cal
G}}_{j,{k+\frac{1}{2}}}&=\frac{b^{+}_{j,{k+\frac{1}{2}}}\bm{G}(\bm{U}^{-}_{j,{k+\frac{1}{2}}})-b^{-}_{j,{k+\frac{1}{2}}}\bm{G}(\bm{U}^{+}_{j,{k+\frac{1}{2}}})}{b^{+}_{j,{k+\frac{1}{2}}}-b^{-}_{j,{k+\frac{1}{2}}}}+\frac{b^{+}_{j,{k+\frac{1}{2}}}b^{-}_{j,{k+\frac{1}{2}}}}{b^{+}_{j,{k+\frac{1}{2}}}-b^{-}_{j,{k+\frac{1}{2}}}}\left(\bm{U}^{+}_{j,{k+\frac{1}{2}}}-\bm{U}^{-}_{j,{k+\frac{1}{2}}}\right)\\\
&+\bm{q}_{j,{k+\frac{1}{2}}}.\end{aligned}$ (2.19)

To obtain $\bm{U}^{\pm}_{{j+\frac{1}{2}},k}$ and
$\bm{U}^{\pm}_{j,{k+\frac{1}{2}}}$ in (2.18)–(2.19), we reconstruct the
second-order piecewise linear interpolant
$\sum_{j,k}\left[\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{j,k}+(\bm{U}_{x})_{j,k}(x-x_{j})+(\bm{U}_{y})_{j,k}(y-y_{k})\right]{\cal
X}_{C_{j,k}}(x,y)$ where $(\bm{U}_{x})_{j,k}$ and $(\bm{U}_{y})_{j,k}$ are the
slopes which are supposed to be computed using a nonlinear limiter to ensure a
non-oscillatory nature of the reconstruction. In the numerical experiments
reported in §3, we have used the generalized minmod limiter with the minmod
parameter $\theta=1.3$. We then follow [8, §3.2] to evaluate the corresponding
point values $u^{\pm}_{{j+\frac{1}{2}},k}$, $u^{\pm}_{j,{k+\frac{1}{2}}}$,
$v^{\pm}_{{j+\frac{1}{2}},k}$, $v^{\pm}_{j,{k+\frac{1}{2}}}$,
$p^{\pm}_{{j+\frac{1}{2}},k}$, and $p^{\pm}_{j,{k+\frac{1}{2}}}$, and estimate
the one-sided local speeds of propagation $a^{\pm}_{{j+\frac{1}{2}},k}$ and
$b^{\pm}_{j,{k+\frac{1}{2}}}$.

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

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

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

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

We then go through all of the steps in the derivation of the 1-D fully
discrete scheme for the systems in (2.20) following [8, §3.2] and §2 up to
(2.4), which now reads as

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

and the corresponding local conservation requirements (2.5) become

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

where we take $\widetilde{u}_{{j+\frac{1}{2}},k}=u^{\rm
int}_{{j+\frac{1}{2}},k}=(\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})^{\rm
int}_{{j+\frac{1}{2}},k}/\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}^{\,\rm
int}_{{j+\frac{1}{2}},k}$. We then proceed as in [8, §3.2], where we enforce
the continuity of $u$ and $p$ across the cell interfaces $x=x_{j+\frac{1}{2}}$
by setting

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

and enforce sharp (yet, non-oscillatory) jumps of the $\rho$\- and $\rho
v$-components. This leads to the following formulae analogous to (2.9):

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

where $a^{\rm int,\pm}_{{j+\frac{1}{2}},k}:=a^{\pm}_{{j+\frac{1}{2}},k}-u^{\rm
int}_{{j+\frac{1}{2}},k}$, and

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

and the point values $\big{(}\rho^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{\ell}$, $\big{(}\rho^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{r}$, $\big{(}(\rho v)^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{\ell}$, and $\big{(}(\rho v)^{\rm
int}_{{j+\frac{1}{2}},k}\big{)}_{r}$ were introduced in [8, Eq. (2.14)].

Next, we proceed as in §2.2 and complete the derivation of the fully discrete
scheme (not shown here for the sake of brevity), and after this, we pass to
the semi-discrete limit and end up with the new LDCU flux (2.18) with the
following “built-in” anti-diffusion term:

$\bm{q}_{{j+\frac{1}{2}},k}=\alpha^{*}_{{j+\frac{1}{2}},k}\big{(}q^{\rho}_{{j+\frac{1}{2}},k},u^{*}_{{j+\frac{1}{2}},k}q^{\rho}_{{j+\frac{1}{2}},k},q^{\rho
v}_{{j+\frac{1}{2}},k},q^{E}_{{j+\frac{1}{2}},k}\big{)}^{\top}.$

Here,

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

and

$\alpha^{*}_{{j+\frac{1}{2}},k}=\left\\{\begin{aligned}
&\frac{a^{+}_{{j+\frac{1}{2}},k}}{a^{*,+}_{{j+\frac{1}{2}},k}}&&\mbox{if}~{}u^{*}_{{j+\frac{1}{2}},k}<0,\\\
&\frac{a^{-}_{{j+\frac{1}{2}},k}}{a^{*,-}_{{j+\frac{1}{2}},k}}&&\mbox{otherwise},\end{aligned}\right.\qquad
a^{*,\pm}_{{j+\frac{1}{2}},k}=a^{\pm}_{{j+\frac{1}{2}},k}-u^{*}_{{j+\frac{1}{2}},k}.$

Similarly, the “built-in” anti-diffusion term in the new $y$-directional LDCU
flux (2.19) is

$\bm{q}_{j,{k+\frac{1}{2}}}=\alpha^{*}_{j,{k+\frac{1}{2}}}\big{(}q^{\rho}_{j,{k+\frac{1}{2}}},q^{\rho
u}_{j,{k+\frac{1}{2}}},v^{*}_{j,{k+\frac{1}{2}}}q^{\rho}_{j,{k+\frac{1}{2}}},q^{E}_{j,{k+\frac{1}{2}}}\big{)}^{\top}$

with

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

and

$\alpha^{*}_{j,{k+\frac{1}{2}}}=\left\\{\begin{aligned}
&\frac{b^{+}_{j,{k+\frac{1}{2}}}}{b^{*,+}_{j,{k+\frac{1}{2}}}}&&\mbox{if}~{}v^{*}_{j,{k+\frac{1}{2}}}<0,\\\
&\frac{b^{-}_{j,{k+\frac{1}{2}}}}{b^{*,-}_{j,{k+\frac{1}{2}}}}&&\mbox{otherwise},\end{aligned}\right.\qquad
b^{*,\pm}_{j,{k+\frac{1}{2}}}=b^{\pm}_{j,{k+\frac{1}{2}}}-v^{*}_{j,{k+\frac{1}{2}}}.$

###### Remark 2.2

As in the 1-D case, the computation of numerical fluxes in (2.18) and (2.19)
should be desingularized to avoid division by zero or very small numbers:

$\bullet$ If $a^{+}_{{j+\frac{1}{2}},k}<\varepsilon$ and
$a^{-}_{{j+\frac{1}{2}},k}>-\varepsilon$ for a small positive $\varepsilon$,
we replace the flux $\bm{{\cal F}}_{{j+\frac{1}{2}},k}$ with

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

$\bullet$ If $b^{+}_{j,{k+\frac{1}{2}}}<\varepsilon$ and
$b^{-}_{j,{k+\frac{1}{2}}}>-\varepsilon$, we replace the flux $\bm{{\cal
G}}_{j,{k+\frac{1}{2}}}$ with

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

##  3 Numerical Examples

In this section, we apply the developed schemes to several initial-boundary
value problems for the 1-D and 2-D Euler equations of gas dynamics (with
$\gamma=1.4$) and compare the performance of the new and original second-order
LDCU schemes, which will be referred to as the NEW and OLD schemes.

For time integration of the ODE systems (2.13) and (2.17), we have used the
three-stage third-order strong stability preserving (SSP) Runge-Kutta method
(see, e.g.,[2, 3]) with the CFL number 0.475. We have taken the small
desingularization parameter $\varepsilon=10^{-12}$.

#### Example 1—Shock-Entropy Problem

In the first example taken from [14], we consider the shock-entropy wave
interaction problem with the following initial condition:

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

which corresponds 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 compute the numerical solution using both the NEW and OLD schemes in the
computational domain $[-5,5]$ on a uniform mesh with $\Delta x=1/80$. We
impose the free boundary conditions by simply setting $\,\hbox{\vbox{\hrule
height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{0}:=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{1}$ and
$\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern
0.0pt}}}_{N+1}:=\,\hbox{\vbox{\hrule height=0.5pt\kern
1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{N}$ in the ghost cells
$C_{0}$ and $C_{N+1}$ on the left and on the right, respectively. The
numerical results at time $t=5$ are presented in Figure 3.1 along with the
corresponding reference solution computed by the NEW scheme on a much finer
mesh with $\Delta x=1/800$. As one can see, the numerical solution computed by
the OLD scheme is very inaccurate near the left boundary of the computational
domain, while the NEW scheme accurately captures the solution throughout the
entire computational domain.

Figure 3.1: Example 1: Density $\rho$ computed by the OLD and NEW schemes
(left) and zoom at $x\in[-5,-2]$ (right).

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

In the second example taken from [5], the initial conditions,

$\big{(}\rho(x,0),u(x,0),p(x,0)\big{)}=\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 $[0,1]$ subject to the free
boundary conditions.

We compute the numerical solutions until the final time $t=0.012$ using both
the NEW and OLD schemes on two uniform meshes, the coarse and fine ones with
$\Delta x=1/200$ and $1/8000$, respectively. The numerical results, plotted in
Figure 3.2, show that the NEW and OLD schemes achieve similar resolutions of
both shock and contact waves. At the same time, the numerical results computed
by the NEW scheme is non-oscillatory, while the OLD scheme solutions contain
small oscillations, whose magnitude does not decay when the mesh is refined.

Figure 3.2: Example 2: Density $\rho$ computed by the OLD and NEW schemes on
two uniform meshes (left) and zoom at $x\in[0.7,0.9]$ and $x\in[0.8,0.801]$
(right).

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

In the first 2-D example, we consider Configuration 3 of the 2-D Riemann
problems taken from [7]; see also [18, 12, 13]. The initial conditions,

$(\rho(x,y,0),u(x,y,0),v(x,y,0),p(x,y,0))=\left\\{\begin{aligned}
&(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{aligned}\right.$

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

We compute the solution using both the NEW and OLD schemes on a uniform mesh
with $\Delta x=\Delta y=0.001$ until the final time $t=1$ and present the
obtained densities in Figure 3.3. As one can see, the NEW solution is
substantially less oscillatory than the OLD one, and the resolution achieved
by the NEW scheme seems to be comparable and even a little better in some
parts of the computed solutions.

Figure 3.3: Example 3: Density ($\rho$) computed by the OLD (left) and NEW
(right) schemes.

#### Example 4—Explosion Problem

In this example, we consider the explosion problem taken from [10] with the
following initial conditions:

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

which are prescribed in the computational domain $[0,1.5]\times[0,1.5]$
subject to the solid wall boundary conditions at $x=0$ and $y=0$ and free
boundary conditions at $x=1.5$ and $y=1.5$.

We compute the solution using the NEW and OLD schemes on a uniform mesh with
$\Delta x=\Delta y=3/800$ until the final time $t=3.2$. The obtained densities
are presented in Figure 3.4, where one can clearly see that there are obvious
oscillations along the boundaries $x=1.5$ and $y=1.5$ in the numerical results
computed by the OLD scheme, while the oscillations are substantially smaller
in the numerical results computed by the NEW scheme. At the same time, in this
example, the OLD scheme achieves a slightly better resolution.

Figure 3.4: Example 4: Density ($\rho$) computed by the OLD (left) and NEW
(right) schemes.

##### Example 5—Implosion Problem

In the last example, we test the implosion problem also taken from [10]. The
initial conditions,

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

are prescribed in the computational domain $[0,0.3]\times[0,0.3]$ subject to
the solid wall boundary conditions.

We compute the solution using the NEW and OLD schemes on a uniform mesh with
$\Delta x=\Delta y=1/2000$ until the final time $t=2.5$. The obtained
densities are presented in Figure 3.5. As one can see, the NEW solution is
slightly less oscillatory than the OLD one. At the same time, the jet
generated by the NEW scheme propagates to a larger extent than the jet
produced by the OLD scheme, which demonstrates that in this example, the NEW
scheme achieves slightly higher resolution than the OLD scheme.

Figure 3.5: Example 5: Density ($\rho$) computed by the OLD (left) and NEW
(right) schemes.

Declarations:

##### Funding.

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

##### Conflicts of interest.

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

##### Data and software availability.

The data that support the findings of this study and FORTRAN codes developed
by the authors and used to obtain all of the presented numerical results are
available from the corresponding author upon reasonable request.

## References

  * [1] 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.
  * [2] 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.
  * [3] S. Gottlieb, C.-W. Shu, and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev., 43 (2001), pp. 89–112.
  * [4] A. Kuganov, S. Noelle, and G. Petrova, Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations, SIAM J. Sci. Comput., 23 (2001), pp. 707–740.
  * [5] A. Kurganov and C.-T. Lin, On the reduction of numerical dissipation in central-upwind schemes, Commun. Comput. Phys., 2 (2007), pp. 141–163.
  * [6] A. Kurganov and E. Tadmor, New high-resolution semi-discrete central schemes for Hamilton-Jacobi equations, J. Comput. Phys., 160 (2000), pp. 720–742.
  * [7] 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.
  * [8] A. Kurganov and R. Xin, New low-dissipation central-upwind scheme, J. Sci. Comput., 96 (2023).  Paper No. 56, 33 pp.
  * [9] 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.
  * [10] R. Liska and B. Wendroff, Comparison of several difference schemes on 1D and 2D test problems for the Euler equations, SIAM J. Sci. Comput., 25 (2003), pp. 995–1017.
  * [11] H. Nessyahu and E. Tadmor, Nonoscillatory central differencing for hyperbolic conservation laws, J. Comput. Phys., 87 (1990), pp. 408–463.
  * [12] C. W. Schulz-Rinne, Classification of the Riemann problem for two-dimensional gas dynamics, SIAM J. Math. Anal., 24 (1993), pp. 76–88.
  * [13] 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.
  * [14] C.-W. Shu and S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys., 77 (1988), pp. 439–471.
  * [15] P. K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal., 21 (1984), pp. 995–1011.
  * [16] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method, J. Comput. Phys., 32 (1979), pp. 101–136.
  * [17] P. Woodward and P. Colella, The numerical solution of two-dimensional fluid flow with strong shocks, J. Comput. Phys., 54 (1988), pp. 115–173.
  * [18] Y. Zheng, Systems of conservation laws: Two-dimensional Riemann problems, Progress in Nonlinear Differential Equations and their Applications, 38, Birkhäuser Boston, Inc., Boston, MA, 2001.