# An asymptotic-preserving multidimensionality-aware finite volume numerical scheme for Euler equations

Alessia Del Grosso, Wasilij Barsukow, Raphaël Loubère, Pierre-Henri Maire

► To cite this version:

Alessia Del Grosso, Wasilij Barsukow, Raphaël Loubère, Pierre-Henri Maire. An asymptotic-preserving multidimensionality-aware finite volume numerical scheme for Euler equations. ICCFD12 - International Conference on Computational Fluid Dynamics, Jul 2024, Kobe, Japan. hal-04736912

HAL Id: hal-04736912

https://inria.hal.science/hal-04736912v1

Submitted on 15 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An asymptotic-preserving multidimensionality-aware finite volume numerical scheme for Euler equations

A. Del Grosso*, W. Barsukow**, R. Loubère** and P.-H. Maire***

Corresponding author: alessia.del-grosso@inria.fr

* Centre Inria de l'université de Bordeaux, Team MEMPHIS, 33400, Talence, France.

** Institut de Mathématiques de Bordeaux (IMB), Université de Bordeaux, UMR 5152, F33400, Talence, France.

*** CEA Cesta, 15 avenue des sablières, Le Barp, France.

**Abstract:** In the context of the numerical approximation of Euler equations, great efforts have been devoted to developing schemes that can accurately reproduce solutions in low Mach number flows. Solutions of classic Finite Volume (FV) schemes are usually plagued by an excessive diffusion as the numerical scheme is not consistent with the limit equations for the Mach number that tends to zero. Instead, a numerical scheme that satisfies such a property is called Asymptotic-Preserving (AP). In this paper, we propose an AP FV scheme for the multidimensional Euler equations. Contrarily to classic multidimensional FV schemes, our method is not obtained by a dimensional splitting. Instead, we include genuine multidimensional effects by exploiting a particular Riemann Solver (RS). It involves a nodal parameter which depends on all the cells around given node. Such an idea has been exploited by Barsukow et al. (2023) for the linear acoustic equations. Their method performs excellently in the low Mach number regime but its extension to the Euler equations proved to be far from trivial. For such a reason, a change of perspective is needed in the definition of the RS.

**Keywords:** Conservation laws, Finite volume schemes, multidimensionality-aware approximate Riemann solver, Low Mach number flows, Asymptotic-preserving scheme.

## 1 Introduction

In the context of the numerical approximation of the compressible Euler equations, great efforts have been devoted to developing numerical schemes that can also accurately reproduce solutions in low Mach number flows, i.e. when the fluid velocity is much smaller than the sound speed. Indeed, there is a wide range of applications for such a problem, from astrophysical stellar evolution to processes in nuclear reactors [1, 2]. However, solutions of classic Finite Volume (FV) schemes are usually affected by an excessive diffusion [3, 4]. As a consequence, to obtain a useful solution, extreme mesh refinement would be required, making the computational cost unreasonably high. In few words, the excessive diffusion is usually caused by a lack of consistency of the numerical scheme with the limit (incompressible) equations for the Mach number *M* that tends to zero [5, 6, 3, 4, 1]. More specifically, when designing numerical schemes for the compressible Euler equations, we impose the discretization to be consistent with the continuous equations, namely the latter are the formal limit of the scheme for $\Delta \to 0$, where $\Delta$ includes both time and space steps. A similar statement can be made for numerical methods for incompressible equations. However, it is not generally true that the discrete compressible equations tend to a consistent discretization of the incompressible ones for $M \to 0$. A numerical scheme that, on the contrary, satisfies such a property is defined as Asymptotic-Preserving (AP) [7, 1, 8, 2], and leads to satisfying results in the low Mach number regime. This concept is also graphically explained in figure 1, where the numerical scheme usually fails in the red limit. Hence, in the literature, we can find papers that "correct" previously existing numerical schemes by using preconditioning strategies [3, 4, 9] or adding *ad-hoc* fixes [10, 11, 12, 13]. The AP property is also specifically addressed in [1, 7, 14, 15, 16], to name but a few. It is also important to highlight that the issue of accurate numerical simulations in low Mach number flows is strictly related to the dimension of the space. In 1D simulations, numerical schemes have a good behavior without the need of any special treatment as the limit to the incompressible equations is trivial. Problems appear when considering multidimensional equations, where the limit involves truly multi-dimensional operators (i.e. the divergence of the velocity) [6, 3]. Hence, it is not surprising that some AP schemes rely on multidimensional information, for instance for the acoustic equations [17, 18]
<figure><img src="image_2.png" /><figcaption>Figure 1: Graphical explanation of the issue of low Mach number flows. The numerical schemes for which the red limit is valid are called asymptotic-preserving.</figcaption></figure>

or the Euler equations [19]. We also emphasize that classic multidimensional FV schemes rely on their 1D counterpart and 1D (approximate or not) Riemann Solver (RS). Specifically, once the 1D solver is designed, a dimensional splitting is performed [20]. Hence, since the problem of low Mach number flows is multi-dimensional, a natural thought may be that a truly multidimensional RS could solve the issue. However, this has been shown to be not true by Barsukow and Klingenberg in [21]. Indeed, they provided the exact solution of the Cauchy problem combined with the truly multidimensional Godunov method for the acoustic equations, only to observe that also such a scheme was affected by the same diffusion flaw as classic methods. Of course, one may argue that, rather than developing AP numerical schemes for the compressible Euler system, one should directly consider the incompressible equations. However, to separate the two regimes is not always possible, as one might be specifically interested in compressible phenomena at low Mach number or in applications which involves both compressible and incompressible phenomena (e.g. one regime in only one part of the domain) [6]. Furthermore, the issue of low Mach number flows for the Euler equations is also automatically transferred to the complete Navier-Stokes equations, making compelling the need to solve it [6].

In this manuscript, we present a numerical scheme for the multidimensional Euler equations that performs particularly well in low Mach number flows, without adding any *ad-hoc* corrections but relying on the structure of the numerical method. Contrarily to classic multidimensional numerical schemes, our method is not obtained by a dimensional splitting. Instead, we include genuine multidimensional effects inside the numerical flux. Such a feature is obtained by exploiting a very particular RS which couples all cells in the vicinity of the current one [22]. This solver is no longer 1D across one edge. as the numerical flux does not depend on only two states (i.e. left and right to a given face) but on a multiple number of states (i.e. the ones around a given node), whose exact number depends on the specific mesh, see also figure 2 for a graphical explanation. This is also why we distinguish between (classic) two-point schemes and (multidimensionality-aware) multi-point methods. Such idea has first been proposed by [22], where the authors exploited a nodal velocity in order to include multidimensional knowledge. However, while such a numerical scheme behaves very well in supersonic and hypersonic flows, it provides very diffusive results in low Mach number flows, just as classic methods. It is therefore natural to wonder whether a different way of including multidimensional knowledge would help to solve the multidimensional problem of low Mach number flows. In the context of linear acoustic equations, the authors of [18] proposed to use a nodal pressure instead of a nodal velocity inside the RS. Such a method was found to be vorticity-preserving and thus, behaved excellently in the low Mach number regime. However, this scheme was developed for linear acoustic equations, and its extension to the Euler equations proved to be far from trivial. A direct extension leads to a numerical scheme that is no longer consistent nor conservative. Therefore, a change of perspective is needed in the definition of the RS, specifically in the definition of the fluxes. Hence, in this paper, by exploiting a nodal pressure inside the solver, we propose a multidimensionality-aware method for Euler equations which is AP. Furthermore, in addition to being AP, the method can also be viewed as all-speed [14, 19], in the sense that it can be used not only for low Mach number flows but also for $M \approx 1$ or higher. We highlight that, even if the method is developed in the context of the Euler equations, it can be extended to other conservation laws by following the same procedure described in [22], where the authors used a nodal velocity.

In the following section 2, we focus on the formulation of both compressible and incompressible Euler equations. Next, in section 3, we describe the approximate RSs on which the numerical methods are based. Specifically, we first recall the original multi-dimensional RS with nodal velocity [22]. Then, we propose another solver with a nodal pressure. Particular care is given to the flux definition. Then, the complete numerical method based on the new RS is presented in section 4, together with the theoretical
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Stencils for the two-point (left) and multi-point (right) flux.</figcaption></figure>

proof of the good behavior of the scheme in low Mach number flows. In section 5, numerical simulations are provided to numerically validate the new multi-point scheme. Some conclusions and perspectives are given in the last section 6.

## 2 The compressible and incompressible Euler equations

In order to have all the ingredients to present and analyze the numerical method we propose, we first go further into details of the Euler system and the discretized domain. We highlight that, even if we only consider the 2D version of the equations, the numerical scheme we propose can easily be generalized to 3D as has been done in [23].

The 2D compressible Euler equations read

$$ \partial_t \mathbf{U} + \nabla \cdot \mathbb{F}(\mathbf{U}) = \mathbf{0}, \qquad (1) $$

with

$$ \mathbf{U} = \begin{pmatrix} \rho \\ \rho \mathbf{u}^T \\ \rho E \end{pmatrix}, \quad \mathbb{F}(\mathbf{U}) = \begin{pmatrix} \rho \mathbf{u}^T \\ \rho \mathbf{u} \otimes \mathbf{u} + p \mathbb{I}_2 \\ \rho E \mathbf{u}^T + p \mathbf{u}^T \end{pmatrix}, $$

where $t$ and $x$ respectively represent time and space. Then, $\rho$ indicates the gas density, $\mathbf{u}^T = (u, v)$ is the velocity vector, $p$ is the pressure term, $E = \varepsilon + \frac{1}{2}\|\mathbf{u}\|^2$ is the total energy where $\varepsilon = \varepsilon(\tau, \eta)$ is the specific internal energy. The latter is assumed to be a strictly convex function with respect to the specific volume $\tau = \frac{1}{\rho}$ and to the specific entropy $\eta$. Given the temperature $\theta > 0$, for the thermodynamic closure of the system, we use the complete equation of state,

$$ p(\tau, \eta) = - \frac{\partial \varepsilon}{\partial \tau} \Big|_{\eta} \quad \text{and} \quad \theta(\tau, \eta) = \frac{\partial \varepsilon}{\partial \eta} \Big|_{\tau} > 0, $$

which also implies the fundamental Gibbs relation, $\theta d\eta = pd\tau + d\varepsilon$. The associated entropy inequality reads

$$ \partial_t (\rho \eta) + \nabla \cdot (\rho \eta \mathbf{u}) \geq 0 \qquad (2) $$

where the equality to 0, i.e. entropy conservation, is obtained for smooth solutions. Finally, by the convexity assumption of $\varepsilon(\tau, \eta)$, we can also define the isentropic sound speed $a$:

$$ \frac{a^2}{\tau^2} = - \frac{\partial p}{\partial \tau} \Big|_{\eta} = \frac{\partial^2 \varepsilon}{\partial \tau^2} \Big|_{\eta} > 0. $$

With this notation, the Mach number $M$ is defined as $M = \frac{\|\mathbf{u}\|}{a}$.
## 2.1 Quasi-1D Euler Equations

Let us now also discuss the quasi-1D version of the equations in a given normal direction **n**. Our numerical flux will rely on an approximate RS based on the quasi 1D system, as to develop 2D approximate RS is known to be far too complicate. However, unlike classic methods, our aim is to introduce multidimensional knowledge into the solver by including a nodal parameter.

Hence, let **t** be the unit vector such that (n, **t**) is an orthonormal basis attached to a generic interface. The normal and tangential components of the velocity vector **u** respectively read $u_n = \mathbf{u} \cdot \mathbf{n}$ and $u_t = \mathbf{u} \cdot \mathbf{t}$, with $\mathbf{u} = u_n\mathbf{n} + u_t\mathbf{t}$. Then, the Euler system in the normal direction **n** reads

$$ \frac{\partial \mathbf{U}_n}{\partial t} + \frac{\partial \mathbf{F}_n(\mathbf{U})}{\partial x_n} = \mathbf{0}, $$

where the variables vector and physical flux projected onto the normal direction **n** are respectively given by

$$ \mathbf{U}_n = \begin{pmatrix} \rho \\ \rho u_n \\ \rho u_t \\ \rho E \end{pmatrix}, \quad \mathbf{F}_n = \begin{pmatrix} \rho u_n \\ \rho u_n^2 + p \\ \rho u_n u_t \\ \rho u_n E + p u_n \end{pmatrix}, $$

with $x_n = \mathbf{x} \cdot \mathbf{n}$. This model has three eigenvalues given by $\Lambda_{n,l} = u_n - a$, $\Lambda_{n,0} = u_n$ (multiplicity 2), $\Lambda_{n,r} = u_n + a$. As we have previously seen, the Euler system is associated with an entropy inequality (2). In this case, it reads

$$ \partial_t(\rho\eta) + \partial_{x_n}(\rho\eta u_n) \geq 0. $$

**Remark 1 (Flux decomposition).** The physical flux can be decomposed in the following way,

$$ \mathbf{F}_n = u_n \mathbf{U}_n + \mathbf{L}_n, \quad \text{with} \quad \mathbf{L}_n = (0, p, 0, pu_n)^T, \qquad (3) $$

where we separated the transport and acoustic contributions. The term $\mathbf{L}_n$ can also be interpreted as the Lagrangian counterpart of the physical flux. We will see that we mimic this relation also at the discrete level.

## 2.2 Incompressible Euler equations

In order to describe the incompressible Euler equations, we first briefly introduce the non-dimensional version of the 2D Euler equations (1). Hence, following [3], we define the non-dimensional variables

$$ \tilde{\rho} = \frac{\rho}{\rho^*}, \quad \tilde{\mathbf{u}} = \frac{\mathbf{u}}{u^*}, \quad \tilde{p} = \frac{p}{p^*}, \quad \tilde{E} = \frac{E}{E^*}, \quad \tilde{t} = \frac{t}{t^*}, \quad \tilde{\mathbf{x}} = \frac{\mathbf{x}}{\delta^*} $$

where $\delta^*$ is an arbitrary length scale,

$$ \rho^* = \max_{\mathbf{x}}(\rho(\mathbf{x})), \quad u^* = \max_{\mathbf{x}}(\|\mathbf{u}(\mathbf{x})\|), \quad p^* = \rho^*(a^*)^2, \quad E^* = \rho^*p^*, \quad t^* = \frac{\delta^*}{u^*}, $$

$\tilde{p} = (\gamma - 1) \left( \tilde{\rho}\tilde{E} - \frac{M^2}{2}\tilde{\rho}(\tilde{u}^2 + \tilde{v}^2) \right)$, $M$ is the reference Mach number $M = \frac{u^*}{a^*}$, and where the new parameters with symbol * represent the characteristic length, time, velocity, etc. Then, the non-dimensional 2D Euler system is given by

$$ \left\{ \begin{array}{l} \partial_t \tilde{\rho} + \tilde{\nabla} \cdot (\tilde{\rho} \tilde{u}) = 0 \\ \partial_t (\tilde{\rho} \tilde{u}) + \tilde{\nabla} \cdot (\tilde{\rho} \tilde{u} \otimes \tilde{u} + \frac{\tilde{p}}{M^2} \mathbb{I}_2) = 0 \\ \partial_t (\tilde{\rho} \tilde{E}) + \tilde{\nabla} \cdot (\tilde{\rho} \tilde{E} \tilde{u} + \tilde{p} \tilde{u}) = 0. \end{array} \right. \qquad (4) $$

At this stage, one should perform an asymptotic expansion in power of $M$ and collect the terms with equal power of $M$. Then, we can conclude that, when the Mach number tends to zero ($M \to 0$), the
<figure><img src="image_5.png" /><figcaption>Figure 3: Geometrical entities attached to the polygonal cell $\omega_c$. Image from [22].</figcaption></figure>

compressible Euler equations (4) tend to the incompressible ones:

$$ \left\{ \begin{aligned} \partial_t \tilde{\rho} + \tilde{u} \tilde{\nabla} \cdot \tilde{\rho} &= 0 \\ \tilde{\nabla} \cdot \tilde{u} &= 0 \\ \partial_t (\tilde{\rho} \tilde{u}) + \tilde{\nabla} \cdot (\tilde{\rho} \tilde{u} \otimes \tilde{u}) + \tilde{\nabla} p_2 &= 0. \end{aligned} \right. \qquad (5) $$

Here, $p_2$ is the dynamic pressure given by the asymptotic expansion of the pressure in power of the Mach number:

$$ \tilde{p}(\mathbf{x}, t) = P_0(t) + M^2 p_2(\mathbf{x}, t) + \mathcal{O}(M^3). $$

Observe that the pressure is constant in space up to fluctuations of order $M^2$. The details are left to the reader, see [3] for the formal derivation of system (5). Further details about the incompressible equations will be given in section 4.2, for the discussion of the AP property of the method.

## 2.3 Notations

Finally, we also provide some notations about the discretization of the domain, see also [22, 24] for similar notations and figure 3 for a graphical explanation.

First of all, the computational domain is a polygon in $\mathbb{R}^2$, divided into a set of non overlapping polygonal cells $\omega_c$. Observe that $c$ is the general label for a cell, whereas $f$ and $p$ respectively indicate a given face (or subface) and point (node) of the domain. Hence, $\mathbf{x}_p$ denotes the vector position of $p$ and the set of vertices (points) of $\omega_c$ is indicated by $\mathcal{P}(c)$. The points of cell $\omega_c$ are counter-clockwise ordered, where $p^-$ and $p^+$ are respectively the previous and the next points with respect to $p$. By joining the cell centroid, $\mathbf{x}_c$, to the midpoints of $[\mathbf{x}_{p^-}, \mathbf{x}_p]$, $[\mathbf{x}_p, \mathbf{x}_{p^+}]$ and to $\mathbf{x}_p$, we get the subcell $\omega_{pc}$. Their union for $p \in \mathcal{P}(c)$ is a partition of the cell $\omega_c$:

$$ \omega_c = \bigcup_{p \in \mathcal{P}(c)} \omega_{pc}. $$

A similar reasoning can be applied to the faces. A given face can be decomposed into subfaces by means of the partition of $c$ induced by the subcells $pc$ for $p \in \mathcal{P}(c)$. Hence, the set of subfaces attached to the corner $pc$ is denoted as $\mathcal{SF}(pc)$. As such, the set of faces of $\omega_c$ is given by

$$ \mathcal{F}(c) = \bigcup_{p \in \mathcal{P}(c)} \mathcal{SF}(pc). $$

We also denote the measure and unit outward normal of a subface $f$ as $l_{pf}$ and $\mathbf{n}_{pf} = (n_x, n_y)_{pf}$. Given a cell $c$ and one of its face $f$, a unique neighbor cell can be associated and we refer to it as $d(c, f)$ or $d$ to shorten the notation. Obviously $c$ is the neighbor cell of $d$ through edge $f$ as well.

Observe that, throughout the paper, we either use the notation $pf$, $pf$ or $f$ depending on whether or not it is necessary to specify the link towards the cell $c$ and node $p$.
## 3 Approximate Riemann solvers for Euler equations

In this section, we lay the foundations of the numerical schemes. We describe three different approximate RS which are then used to compute the associated numerical flux. The aim is to compare these solvers which, although they can be inscribed in the same framework, can lead to completely different numerical results depending on whether (and which) nodal parameters are present. The first solver is called two-point RS: it has been presented in [25] and is similar to the HLLC one with different wave speeds approximation. The second RS is the multi-point one with nodal velocity [22]. This solver incorporates multi-dimensional knowledge thanks to the presence of a nodal parameter: the nodal velocity. More specifically, it allows to couple all the Riemann problems related to the surfaces attached to the considered node. When the nodal velocity coincides with the classic Godunov velocity, the multi-point RS reduces itself to the more classic two-point solver. However, both solvers produce very diffusion solutions in low Mach number flows. Hence, the novelty of the present paper is the third solver, which is a multi-point one but with a nodal pressure instead of a nodal velocity. The idea is to enhance the numerical results for low Mach number flows as it has been done in the context of the linear acoustic equations in [18]. We will see that this solver also reduces to the two-point one when using the Godunov face-based pressure instead of the nodal one.

### 3.1 Common ground for the approximate Riemann solvers

We start by laying the common ground to the three different approximate RSs and then further specify the differences in the next sections. Observe that this presentation is different from the one used in [25, 22], which was based on a Lagrange-Euler interpretation, see remark 2.

Given a surface $f$ impinging at node $p$, the Eulerian Riemann problem in the outward normal direction $n_{pf} \equiv n$ reads

$$ (RP): \left\{ \begin{array}{l} \frac{\partial \mathbf{U}_n}{\partial t} + \frac{\partial[\mathbf{F}_n(\mathbf{U})]}{\partial x_n} = \mathbf{0}, \\ \mathbf{U}_n(x_n, 0) = \begin{cases} \mathbf{U}_{n,lf} & \text{if } x_n < 0, \\ \mathbf{U}_{n,rf} & \text{if } x_n \ge 0. \end{cases} \end{array} \right. $$

where $\mathbf{U}_{n,lf}$, $\mathbf{U}_{n,rf}$ are the left and right states to the given face $f$ in Eulerian coordinates. Assuming the wave speeds to be a priori ordered $\Lambda_{n,lf} < \Lambda_{n,0f} < \Lambda_{n,rf}$, we assume the approximate solution $\mathbf{W}_{pf}$ of $(RP)$ to be composed of four states separated by the three discontinuities:

$$ \mathbf{W}_{pf}(\mathbf{U}_{n,lf}, \mathbf{U}_{n,rf}, \frac{x_n}{t}, \mathbf{n}_{pf}, \Psi_p) = \begin{cases} \mathbf{U}_{n,lf} & \text{if } \frac{x_n}{t} \le \Lambda_{n,lf}, \\ \mathbf{U}_{n,lf}^* & \text{if } \Lambda_{n,lf} < \frac{x_n}{t} \le \Lambda_{n,0f}, \\ \mathbf{U}_{n,rf}^* & \text{if } \Lambda_{n,0f} < \frac{x_n}{t} \le \Lambda_{n,rf}, \\ \mathbf{U}_{n,rf} & \text{if } \Lambda_{n,rf} \le \frac{x_n}{t}. \end{cases} \quad (6) $$

See also figure 4 for a graphical explanation. $\Psi_p$ indicates the dependence to a nodal parameter which will be specified further. The assumption of ordered wave speeds is actually common in the literature [25]. We will later see that this assumption will be naturally satisfied by our definitions for the wave speeds, which will be based on positivity-preserving and entropy-stable properties. As usual, the intermediate star values $\mathbf{U}_{n,sf}^*$ are given by

$$ \mathbf{U}_{n,sf}^* = \begin{pmatrix} \rho_{sf}^* \\ \rho_{sf}^* u_{n,sf}^* \\ \rho_{sf}^* u_{t,sf}^* \\ \rho_{sf}^* E_{sf}^* \end{pmatrix}, \quad (7) $$

with $s = l, r$. We also observe that $E_{sf}^* = \varepsilon_{sf}^* + \frac{1}{2}((u_{n,sf}^*)^2 + (u_{t,sf}^*)^2)$. Given the left and right fluxes $\mathbf{F}_{n,sf} = \mathbf{F}_n(\mathbf{U}_{n,sf})$, $s = l, r$, we also define the intermediate fluxes $\hat{\mathbf{F}}_{n,sf}$ by analogy with the flux decomposition (3) satisfied by the physical flux:

$$ \hat{\mathbf{F}}_{n,sf} = \hat{u}_{n,sf} \mathbf{U}_{n,sf}^* + \hat{\mathbf{L}}_{n,sf}, \quad (8) $$
<figure><img src="image_6.png" /><figcaption>Figure 4: Graphical explanation of the solution of the Riemann problem.</figcaption></figure>

where

$$ \hat{\mathbf{L}}_{\mathbf{n},sf} = (0, \hat{p}_{sf}, 0, \hat{p}_{sf}\hat{u}_{\mathbf{n},sf})^T, \quad (9) $$

and with $\hat{p}_{sf}, \hat{u}_{\mathbf{n},sf}, s = l, r$, additional unknowns. Hence, definitions (7) and (8) lead to a total of 8 + 4 unknowns that need to be found. We do not consider the wave speeds $\Lambda_{\mathbf{n},lf}, \Lambda_{\mathbf{n},0f}, \Lambda_{\mathbf{n},rf}$ as unknowns, but given parameters whose definition will be given as we move forward. Observe that definitions (8)-(9) are not common. In fact, the following easier formula is often used,

$$ \mathbf{F}_{\mathbf{n},sf}^* = u_{\mathbf{n},sf}^* \mathbf{U}_{\mathbf{n},sf}^* + \mathbf{L}_{\mathbf{n},sf}^*, \quad \text{with} \quad \mathbf{L}_{\mathbf{n},sf}^* = (0, p_{sf}^*, 0, p_{sf}^* u_{\mathbf{n},sf}^*)^T, $$

that is to say $\hat{X} \equiv X^*$ with $X$ any of the given variables (e.g. $u_n, p$), see for instance [25, 22]. Indeed, this choice would lead to only 8 + 2 unknowns instead of the previous 12. However, in this paper and depending on the approximate RS, we may need additional freedom and consider the general case in which $\hat{X}$ and $X^*$ might not coincide. As the reader might suspect, this choice is actually related to the need to include a nodal pressure in the RS. Indeed, the use of the more classic definition $\hat{X} \equiv X^*$ together with a nodal pressure could lead to a RS which is no longer consistent. Hence, a change of paradigm is needed. More specifically, it is necessary to "lighten" the classic assumptions in order to obtain a satisfying RS. Further details will be given later when discussing the RS with nodal pressure, see section 3.4.

At this stage, in order to determine the star values, it is classic to assume the two external waves to be discontinuities and impose the following jump conditions across them:

$$ -\Lambda_{\mathbf{n},lf}(\mathbf{U}_{\mathbf{n},lf}^* - \mathbf{U}_{\mathbf{n},lf}) + \hat{\mathbf{F}}_{\mathbf{n},lf} - \mathbf{F}_{\mathbf{n},lf} = \mathbf{0} \quad \text{and} \quad -\Lambda_{\mathbf{n},rf}(\mathbf{U}_{\mathbf{n},rf} - \mathbf{U}_{\mathbf{n},rf}^*) + \mathbf{F}_{\mathbf{n},rf} - \hat{\mathbf{F}}_{\mathbf{n},rf} = \mathbf{0}, \quad (10) $$

providing 8 relations, which can also be reformulated as

$$ \begin{aligned} & - \mathbf{U}_{\mathbf{n},lf}^*(\Lambda_{\mathbf{n},lf} - \hat{u}_{\mathbf{n},lf}) + \mathbf{U}_{\mathbf{n},lf}(\Lambda_{\mathbf{n},lf} - u_{\mathbf{n},lf}) + \hat{\mathbf{L}}_{\mathbf{n},lf} - \mathbf{L}_{\mathbf{n},lf} = \mathbf{0}, \\ & - \mathbf{U}_{\mathbf{n},rf}(\Lambda_{\mathbf{n},rf} - u_{\mathbf{n},rf}) + \mathbf{U}_{\mathbf{n},rf}^*(\Lambda_{\mathbf{n},rf} - \hat{u}_{\mathbf{n},rf}) + \mathbf{L}_{\mathbf{n},rf} - \hat{\mathbf{L}}_{\mathbf{n},rf} = \mathbf{0}, \end{aligned} \quad (11) $$

by using definition (8). Focusing on the mass equation, (11) simply gives

$$ - \rho_{lf}^*(\Lambda_{\mathbf{n},lf} - \hat{u}_{\mathbf{n},lf}) + \rho_{lf}(\Lambda_{\mathbf{n},lf} - u_{\mathbf{n},lf}) = 0, \quad \text{and} \quad - \rho_{rf}(\Lambda_{\mathbf{n},rf} - u_{\mathbf{n},rf}) + \rho_{rf}^*(\Lambda_{\mathbf{n},rf} - \hat{u}_{\mathbf{n},rf}) = 0. $$

By assuming that there exist two real parameters $\lambda_{lf}, \lambda_{rf} > 0$ such that

$$ - \lambda_{lf} = \rho_{lf}^*(\Lambda_{\mathbf{n},lf} - \hat{u}_{\mathbf{n},lf}) = \rho_{lf}(\Lambda_{\mathbf{n},lf} - u_{\mathbf{n},lf}), \quad \text{and} \quad \lambda_{rf} = \rho_{rf}^*(\Lambda_{\mathbf{n},rf} - \hat{u}_{\mathbf{n},rf}) = \rho_{rf}(\Lambda_{\mathbf{n},rf} - u_{\mathbf{n},rf}), $$

we get a natural definition for the wave speeds $\Lambda_{\mathbf{n},lf}, \Lambda_{\mathbf{n},rf}$:

$$ \Lambda_{\mathbf{n},lf} = \hat{u}_{\mathbf{n},lf} - \frac{\lambda_{lf}}{\rho_{lf}^*} = u_{\mathbf{n},lf} - \frac{\lambda_{lf}}{\rho_{lf}}, \quad \text{and} \quad \Lambda_{\mathbf{n},rf} = \hat{u}_{\mathbf{n},rf} + \frac{\lambda_{rf}}{\rho_{rf}^*} = u_{\mathbf{n},rf} + \frac{\lambda_{rf}}{\rho_{rf}}. \quad (12) $$

Hence, we have transferred the problem of defining the wave speeds $\Lambda_{\mathbf{n},lf}, \Lambda_{\mathbf{n},rf}$ to the two parameters $\lambda_{lf}, \lambda_{rf}$, which will be given according to positivity-preserving and entropy-stability properties. Observe that $\lambda_{lf}, \lambda_{rf}$ are interpreted as approximations of $\rho a$ so that $\Lambda_{\mathbf{n},lf}, \Lambda_{\mathbf{n},rf}$ are indeed approximations of
$u_n \pm a$. Furthermore, the previous formulas also provide the compatibility relations for the mass equation,

$$ \lambda_{lf} \left( \frac{1}{\rho_{lf}^*} - \frac{1}{\rho_{lf}} \right) - (\hat{u}_{n,lf} - u_{n,lf}) = 0, \quad \text{and} \quad -\lambda_{rf} \left( \frac{1}{\rho_{rf}} - \frac{1}{\rho_{rf}^*} \right) - (u_{n,rf} - \hat{u}_{n,rf}) = 0, \quad (13) $$

which are also known as the (H1) condition for the left and right waves, namely the weak form of the Lagrangian volume conservation equation across the left and right discontinuities [26]. At this stage, we are able to reformulate the left/right jump conditions (10) as follows,

$$ \lambda_{lf} \left( \frac{\mathbf{U}_{n,lf}^*}{\rho_{lf}^*} - \frac{\mathbf{U}_{n,lf}}{\rho_{lf}} \right) + \hat{\mathbf{L}}_{n,lf} - \mathbf{L}_{n,lf} = \mathbf{0}, \quad \text{and} \quad -\lambda_{rf} \left( \frac{\mathbf{U}_{n,rf}}{\rho_{rf}} - \frac{\mathbf{U}_{n,rf}^*}{\rho_{rf}^*} \right) + \mathbf{L}_{n,rf} - \hat{\mathbf{L}}_{n,rf} = \mathbf{0}, \quad (14) $$

which read

$$ \left\{ \begin{array}{l} 0 = 0, \\ \lambda_{lf}(u_{n,lf}^* - u_{n,lf}) + \hat{p}_{lf} - p_{lf} = 0, \\ \lambda_{lf}(u_{t,lf}^* - u_{t,lf}) = 0, \\ \lambda_{lf}(E_{lf}^* - E_{lf}) + (\hat{p}_{lf}\hat{u}_{n,f} - p_{lf}u_{n,lf}) = 0, \end{array} \right. \quad \left\{ \begin{array}{l} 0 = 0, \\ \lambda_{rf}(u_{n,rf}^* - u_{n,rf}) - (\hat{p}_{rf} - p_{rf}) = 0, \\ \lambda_{rf}(u_{t,rf}^* - u_{t,rf}) = 0, \\ \lambda_{rf}(E_{rf}^* - E_{rf}) - (\hat{p}_{rf}\hat{u}_{n,f} - p_{rf}u_{n,rf}) = 0. \end{array} \right. \quad (15) $$

Observe that formulas (15) provide only 6 relations for 12 unknowns, the other two being redundant with the definition of $\Lambda_{n,sf}$ and $\lambda_{sf}$, $s = l, r$.

Finally, it remains to impose the jump condition across the 0-wave as we still need to find 4 relations together with a definition for the middle Eulerian wave speed $\Lambda_{0f}$. Depending on what we apply across the middle wave $\Lambda_{0f}$, we obtain different RSs. A special role is played by the balance across the middle wave $\Lambda_{0f}$:

$$ \mathbf{B}_{n,f}^0 = -\Lambda_{n,0f}(\mathbf{U}_{n,rf}^* - \mathbf{U}_{n,lf}^*) + \hat{\mathbf{F}}_{n,rf} - \hat{\mathbf{F}}_{n,lf} \\ = -\mathbf{U}_{n,rf}^*(\Lambda_{n,0f} - \hat{u}_{n,rf}) + \mathbf{U}_{n,lf}^*(\Lambda_{n,0f} - \hat{u}_{n,lf}) + \hat{\mathbf{L}}_{n,rf} - \hat{\mathbf{L}}_{n,lf}. \quad (16) $$

As for now, we only impose the jump condition across the middle wave for the mass equation, which reads

$$ \rho_{rf}^*(\Lambda_{n,0f} - \hat{u}_{n,rf}) = \rho_{lf}^*(\Lambda_{n,0f} - \hat{u}_{n,lf}). \quad (17) $$

Condition (17) is simply satisfied by imposing

$$ \Lambda_{n,0f} = \hat{u}_{n,rf} = \hat{u}_{n,lf} = \hat{u}_{n,f}, \quad (18) $$

which also provides the definition of the middle Eulerian wave speed $\Lambda_{0f}$. This also implies a simplification of the balance term (16) into

$$ \mathbf{B}_{n,f}^0 = \hat{\mathbf{L}}_{n,rf} - \hat{\mathbf{L}}_{n,lf}. $$

However, remark that condition (18) does not necessarily imply $u_{n,lf}^* = u_{n,rf}^*$. Furthermore, thanks to definition (18), formulas (12) for $\Lambda_{n,lf}$, $\Lambda_{n,rf}$ and by assuming the positivity of the star values of the specific volumes, we also obtain that the wave speeds are automatically ordered: $\Lambda_{n,lf} < \Lambda_{n,0f} < \Lambda_{n,rf}$. This justifies the previous definition of the Riemann solution (6). For the condition for the positivity of the intermediate states of the specific volumes, see remark 3. We also refer to [25, 22] for more details. Remark that, at this stage, we still have 3 unknowns and we may write all the star values as functions of $\hat{u}_{n,f}$, $\hat{p}_{lf}$, $\hat{p}_{rf}$ (and, of course, parameters $\lambda_{lf}$, $\lambda_{rf}$).

**Remark 2 (The Lagrange-Euler interpretation).** The RSs in this article could be seen as the Eulerian version of Lagrangian RSs. Such Eulerian RSs would be obtained by exploiting a discrete version of the relations at the continuous level when defining the passage from Lagrangian to Eulerian coordinates. For instance, the terms $\mathbf{L}_{n,sf}$ would correspond to the Lagrangian counterpart of the fluxes, whereas $-\lambda_{lf}$, $\lambda_{rf}$ would be interpreted as the left/right Lagrangian wave speeds. For more details about the Lagrange-Euler transformation, see [25, 22]. In this paper, for simplicity, we directly present the schemes in Eulerian coordinates.

**Remark 3 (Positivity of the star values of the specific volumes).** In order to have a physically-consistent RS, we also ask for positivity of the star values of the specific volumes, namely $\frac{1}{\rho_{sf}^*} > 0$, with
s = l, r. Following [22], to enforce the positivity, we use Dukowicz definition [23, 27]:

$$ \lambda_{lf} \geq \rho_{lf} a_{lf} \left( 1 + G \frac{(\hat{u}_{n,lf} - u_{n,lf})^{(-)}}{a_{lf}} \right), \quad \text{and} \quad \lambda_{r} \geq \rho_{rf} a_{rf} \left( 1 + G \frac{(\hat{u}_{n,rf} - u_{n,rf})^{(+)}}{a_{rf}} \right), \quad (19) $$

with $G = \frac{\gamma+1}{2}$ and, given $X \in \mathbb{R}$, $X^{(-)} = \frac{1}{2}(|X| - X)$, $X^{(+)} = \frac{1}{2}(|X| + X)$. However, this definition is implicit, and it requires an iterative procedure. In this paper, we only use it with two iterations, as it is in practice enough to enforce the positivity of the star values of the specific volumes. Further information about the wave speeds will be given later.

We have completed the presentation of the common parts among the different solvers. Let us now specify how each solver is designed from this point on and highlight their strengths and weaknesses.

## 3.2 Two-point Riemann solver

The classic two-point RS which has been described in [25]. We highlight that this solver is very close to the HLLC one [28], but the main differences rely on the wave speeds definition as well as its derivation.

First of all, for this solver, we consider the usual assumption $\hat{X} \equiv X^*$ with $X = p, u_n$, so that $\hat{L}_{n,sf} = L_{n,sf}^* = (0, p_{sf}^*, 0, p_{sf}^* u_{n,f}^*)^T$. Observe that $p_{sf}^*$ is simply an equivalent notation for $\hat{p}_{sf}$ as $p_{sf}^*$ are not present in the star vectors $U_{n,sf}^*$. Hence, we are left with only 1 missing relation which is simply retrieved by imposing the classic face-based jump conditions for the middle wave as well. Namely, at each (sub)face $f$, we impose

$$ B_{n,f}^0 = 0, \quad (20) $$

which implies

$$ L_{n,lf}^* = L_{n,rf}^* \quad (21) $$

and thus

$$ p_{lf}^* = p_{rf}^* = p_{f}^*. $$

Observe that condition (20) also implies relation $u_{n,lf}^* = u_{n,rf}^* = u_{n,f}^*$, which we already have thanks to the equivalency $\hat{X} \equiv X^*$. We recall that all these jump conditions imply the local conservation through waves, and thus the classic (Eulerian) consistency conditions

$$ -\Lambda_{n,lf}(\mathbf{U}_{n,lf}^* - \mathbf{U}_{n,lf}) - \Lambda_{n,0f}(\mathbf{U}_{n,rf}^* - \mathbf{U}_{n,lf}^*) - \Lambda_{n,rf}(\mathbf{U}_{n,rf} - \mathbf{U}_{n,rf}^*) + \mathbf{F}_{n,rf} - \mathbf{F}_{n,lf} = 0. \quad (22) $$

Hence, after few computations, we obtain the classic Godunov star values for velocity and pressure:

$$ u_{n,f}^* = u_{n,f}^{*,2P} = \frac{\lambda_{lf}u_{n,lf} + \lambda_{rf}u_{n,rf}}{\lambda_{lf} + \lambda_{rf}} - \frac{(p_{rf} - p_{lf})}{\lambda_{rf} + \lambda_{lf}}, \quad (23) $$

$$ p_f^* = p_f^{*,2P} = \frac{\lambda_{lf}p_{rf} + \lambda_{rf}p_{lf}}{\lambda_{lf} + \lambda_{rf}} - \frac{\lambda_{rf}\lambda_{lf}}{\lambda_{rf} + \lambda_{lf}}(u_{n,rf} - u_{n,lf}). $$

Observe that the superscript "2P" is only to highlight the star values used in the two-point version of the method. The other star values $\rho_{sf}^*, E_{sf}^*$ with $s = l, r$ can be found from $u_{n,f}^*$ and $p_f^*$ using jump conditions (13) and (15).

We highlight that the parameters $\lambda_{lf}, \lambda_{rf}$ for the two-point RS are chosen in such a way to ensure positivity of the star values of the specific volumes and internal energy as well as entropy stability. Hence, we need to impose both

$$ \lambda_{sf} \geq \frac{a^2(\tilde{\tau}_s, \eta_s)}{\tilde{\tau}_s^2} \quad \text{with} \quad \tilde{\tau}_s \in [\tau_s, \tau_s^*], \quad (24) $$

with $s = l, r$ and conditions (19), see [25] for more details. Observe that, in the case of the two-point solver, conditions (19) also imply the positivity of the star values of the specific internal energy. However, exploiting them requires an iterative procedure which may be computationally demanding. In [25], the authors providentially observed that it is actually possible to enforce the positivity of $1/\rho_{lf}^*, 1/\rho_{rf}^*$ and
$\epsilon_{lf}^*, \epsilon_{rf}^*$ by directly using the following wave speeds definition:

$$ 
\begin{aligned}
\lambda_{lf} &= \max \left( \rho_{lf} a_{lf}, \sqrt{\max(\rho_{lf}[\Delta p], 0)}, -[\Delta u]\rho_{lf} \right), \\
\lambda_{rf} &= \max \left( \rho_{rf} a_{rf}, \sqrt{\max(-\rho_{rf}[\Delta p], 0)}, -[\Delta u]\rho_{rf} \right),
\end{aligned}
\quad (25) $$

with $\Delta X = X_r - X_l$, with $X$ any of the variables.

## 3.3 Multi-point Riemann solver with nodal velocity

Next, we recall the multi-point RS with nodal velocity as it has been described in [22] for conservation laws and in [24] to handle the presence of source terms.

As for the two-point RS, we assume $\hat{\mathbf{L}}_{n,sf} = \mathbf{L}_{n,sf}^* = (0, p_{sf}^*, 0, p_{sf}^* u_{n,f}^*)^T$. Hence, once again, we are left with only 1 unknown. However, since we aim to include multidimensional knowledge inside the solver (by using a nodal velocity $\mathbf{u}_p$), we do not impose the typical face-based condition (21) for the middle wave. Instead, given a node $p$, we ask for a node-based version of such condition, namely

$$ \sum_{f \in SF(p)} l_{pf} \mathbf{B}_{n,f}^0 = 0, \quad (26) $$

which also reads

$$ \sum_{f \in SF(p)} l_{pf}(\mathbf{L}_{n,rf}^* - \mathbf{L}_{n,lf}^*) = 0. \quad (27) $$

In this way, we are coupling all the solutions of the Riemann problems around the node $p$, i.e. the Riemann problems defined at each subface that impinges on the node $p$. Furthermore, to impose the node-based jump conditions (27) together with the jump conditions (14) is equivalent to impose the node-based consistency conditions,

$$ \sum_{f \in SF(p)} l_{pf} \left( \lambda_{lf} \left( \frac{\mathbf{U}_{n,lf}^*}{\rho_{lf}^*} - \frac{\mathbf{U}_{n,lf}}{\rho_{lf}} \right) - \lambda_{rf} \left( \frac{\mathbf{U}_{n,rf}}{\rho_{rf}} - \frac{\mathbf{U}_{n,rf}^*}{\rho_{rf}^*} \right) + \mathbf{L}_{n,rf} - \mathbf{L}_{n,lf} \right) = 0, \quad (28) $$

which are equivalent to

$$ \sum_{f \in SF(p)} l_{pf} \left( - \Lambda_{n,lf}(\mathbf{U}_{n,lf}^* - \mathbf{U}_{n,lf}) - \Lambda_{n,0f}(\mathbf{U}_{n,rf}^* - \mathbf{U}_{n,lf}^*) - \Lambda_{n,rf}(\mathbf{U}_{n,rf} - \mathbf{U}_{n,rf}^*) \right. \\ \left. + \mathbf{F}_{n,rf} - \mathbf{F}_{n,lf} \right) = 0. \quad (29) $$

Observe that, if we apply the classic jump condition (20) across the middle wave at each subface $f$ around the node $p$, then the nodal jump condition (26) is also verified. However, the vice versa is not true: we are imposing a less restrictive condition. It is also worth pointing out that, since the RS does not satisfy the classic consistency conditions (22), the scheme associated to the RS is not locally conservative at each subface. However, the nodal consistency conditions (29) are sufficient to prove that the associated scheme is locally conservative at each node and, thus, in the entire domain [22]. Few computations show that the node-based conditions (28) lead to the node-based jump condition related to the pressure jump:

$$ \sum_{f \in SF(p)} l_{pf}(p_{rf}^* - p_{lf}^*) \mathbf{n}_{pf} = 0, \quad (30) $$

or equivalently

$$ \sum_{f \in SF(p)} l_{pf}(\lambda_{lf} + \lambda_{rf})(u_{n,f}^* - u_{n,f}^{*,2P}) \mathbf{n}_{pf} = 0, \quad (31) $$

where we have used the jump conditions (15) and (23). Observe that these relations (30) imply that, generally, $p_{rf}^* - p_{lf}^* \neq 0$ at each subface $f$. Hence, contrarily to the two-point solver, there are two distinguished star pressures. At this stage, $u_{n,f}^* \equiv u_{n_{pf}}^*$ is an unknown parameter attached to each subface impinging at node $p$ and the number of unknown parameters at node $p$ is much greater than the number of equations given by the nodal condition (30). To close this system of equations, we follow
[22] and assume that $u_{n,pf}^*$ corresponds to the projection of the unknown nodal vector $\mathbf{u}_p$ onto the unit normal $\mathbf{n}_{pf}$. That is, for all $p$ and $f$, we impose

$$u_{n,pf}^* = \mathbf{u}_p \cdot \mathbf{n}_{pf}.$$

This fundamental assumption drastically reduces the number of unknowns to one vector unknown $\mathbf{u}_p$ per node, which can be interpreted as an approximation of the nodal velocity. Hence, we are able to express conditions (31) as

$$\sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf}(\lambda_{lf} + \lambda_{rf})(\mathbf{u}_p \cdot \mathbf{n}_{pf} - u_{n,f}^{*,2P})\mathbf{n}_{pf} = \mathbf{0}. \qquad (32)$$

Finally, the node-based condition (30) (or equivalently (32)) becomes the linear system $M_u \mathbf{u}_p = \mathbf{w}_u$, where $\mathbf{u}_p$ is the unknown and

$$M_u = \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf}(\lambda_{lf} + \lambda_{rf})(\mathbf{n}_{pf} \otimes \mathbf{n}_{pf}), \quad \mathbf{w}_u = \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf}(\lambda_{lf} + \lambda_{rf})u_{n,f}^{*,2P}\mathbf{n}_{pf}.$$

This system, called a nodal solver [29, 30, 22], always admits a unique solution and, thus, provides an approximation of the nodal velocity $\mathbf{u}_p$, given physical states and wave speeds for all faces impinging at current node. Finally, concerning the parameters $\lambda_{lf}, \lambda_{rf}$, we simply impose (19) and (24). Indeed, as for the two-point scheme, these conditions are sufficient to ensure positivity- and entropy-stability.

This solver does not exhibit a satisfying behavior in low Mach number flows as its solution is affected by an excessive diffusion as for classic FV methods.

## 3.4 Multi-point Riemann solver with nodal pressure

In the previous section, an approximate RS with nodal velocity has been described. Specifically, the RS is equipped with two star pressures and one star velocity which is the projection of the nodal velocity onto the normal direction to the subface. In this section, we develop an approximate RS with a nodal pressure instead of a nodal velocity. We aim at obtaining a numerical scheme suited for low Mach number flows, inspired by what has been done in [18] for the linear acoustic equations.

The idea consists in having a unique (nodal) star pressure but two different star values for the velocity. However, following the previous sections, we immediately remark that, if we assume both $\hat{\mathbf{L}}_{n,sf} = \mathbf{L}_{n,sf}^*$ and $p_{lf}^* = p_{rf}^* = p_f^*$, we retrieve the two-point RS of section 3.2. Thus, it clearly appears that we generally need to distinguish between $\hat{\mathbf{L}}_{n,sf}$ and $\mathbf{L}_{n,sf}$. However, since we impose a unique pressure value

$$\hat{p}_{lf} = \hat{p}_{rf} = \hat{p}_f,$$

together with condition (18), we automatically retrieve the usual face-based jump condition $\hat{\mathbf{L}}_{n,lf} = \hat{\mathbf{L}}_{n,rf}$ across the middle wave.

We are now left with the question of what relation could we impose to find the pressure $\hat{p}_f$. Indeed, we still have two missing relations and the jump condition for the middle wave does not seem to be helpful in this case. This is not surprising as we have introduced two additional unknowns in the RS.

For the pressure, let us mimic what has been done for the previous RS with nodal velocity. There, the star velocity is recovered from the jump of the star pressures around the node (30)-(31). Here, the pressure $\hat{p}_f$ is obtained from the jump of the star velocities around the node:

$$\sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf}(u_{n,rf}^* - u_{n,lf}^*) = 0, \qquad (33)$$

which also reads

$$\sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} \left( \left( \frac{1}{\lambda_{lf}} + \frac{1}{\lambda_{rf}} \right) \hat{p}_f - \left( \frac{p_{rf}}{\lambda_{rf}} + \frac{p_{lf}}{\lambda_{lf}} - (u_{n,rf} - u_{n,lf}) \right) \right) = \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} \left( \frac{1}{\lambda_{lf}} + \frac{1}{\lambda_{rf}} \right) \left( \hat{p}_f - p_f^{*,2P} \right) = 0, \qquad (34)$$
using the left/right jump conditions (14). While, in the RS with nodal velocity, $u_{n,f}^*$ did not generally coincide with $u_{n,f}^{*,2P}$, here $\hat{p}_f$ is not equal to $p_f^{*,2P}$ a priori. Furthermore, likewise (30) could be interpreted as an approximation of $\nabla p$ in the dual cell $\bigcup_{c \in \mathcal{C}(p)} \omega_{pc}$, formula (33) can be seen as an approximation of $\nabla \cdot u$ in the dual cell. However, now we only have one equation for an undefined number of unknowns $p_f^*$, the number of which depends on the number of subfaces around the node $p$. Following [18], we simply propose to define a unique star value for the pressure for all the subfaces around a given node $p$, namely

$$ \hat{p}_f = \pi_p \quad \forall f \in \mathcal{SF}(p). \qquad (35) $$

Therefore, solving (34) with (35), the nodal pressure value $\pi_p$ is simply given by

$$ \pi_p = \frac{w_p}{M_p}, \qquad (36) $$

where

$$ M_p = \sum_{f \in \mathcal{SF}(p)} l_{pf} \left( \frac{1}{\lambda_{lf}} + \frac{1}{\lambda_{rf}} \right), \quad w_p = \sum_{f \in \mathcal{SF}(p)} l_{pf} \left( \frac{p_{rf}}{\lambda_{rf}} + \frac{p_{lf}}{\lambda_{lf}} - (u_{n,rf} - u_{n,lf}) \right). $$

Observe that $M_p$ is always strictly positive provided that $\lambda_{sf} > 0$, $s = l, r$, hence there is no ambiguity in the definition of $\pi_p$.

**Remark 4 (The choice of the nodal solver).** Remark that another possible nodal solver for the pressure may be retrieved by writing down the pressure equation at the continuous level. Then, its discretization would lead to a slightly different nodal solver for the pressure $\pi_p$. This strategy has actually been developed but it produces very similar results. Hence, we only use formula (36) as no advantages were obtained. The idea of exploiting the pressure equation in the RS is not new. See for instance [31], where a relaxation version has been used.

**Remark 5 (Jump of the velocity around a given node).** Following [18], we highlight that the velocity diffusion part of the nodal pressure $\pi_p$ defined by formula (34)-(36) is a weak consistent discretization of the divergence of $u$ at the node $p$. First of all, we observe that the velocity part of the left-hand-side of (34) is

$$ \sum_{f \in \mathcal{SF}(p)} l_{pf} (u_{n,rf} - u_{n,lf}) = - \sum_{c \in \mathcal{C}(p)} \left( \sum_{f \in \mathcal{SF}(pc)} l_{cf} n_{cf} \right) \cdot \mathbf{u}_c = - \sum_{c \in \mathcal{C}(p)} l_{pc} n_{pc} \cdot \mathbf{u}_c, $$

with

$$ l_{pc} n_{pc} = \sum_{f \in \mathcal{SF}(pc)} l_{cf} n_{cf}, $$

and where $l_{cf}$ and $n_{cf}$ are respectively the length and the outward-pointing normal at the face $f$ pointing out of cell $c$. Hence, we get

$$ \left( \sum_{f \in \mathcal{SF}(p)} l_{pf} \left( \frac{1}{\lambda_{lf}} + \frac{1}{\lambda_{rf}} \right) \right) \pi_p = \left( \sum_{f \in \mathcal{SF}(p)} l_{pf} \left( \frac{p_{rf}}{\lambda_{rf}} + \frac{p_{lf}}{\lambda_{lf}} \right) \right) - |\omega_p| D_p(\mathbf{u}), $$

where

$$ \omega_p = \bigcup_{c \in \mathcal{C}(p)} \omega_{pc}, \quad \text{and} \quad D_p(\mathbf{u}) = - \frac{1}{|\omega_p|} \sum_{c \in \mathcal{C}(p)} l_{pc} n_{pc} \cdot \mathbf{u}_c, $$

with $\mathcal{C}(p)$ the set of cells around node $p$. For the proof that $D_p$ is a weak consistent discretization of the divergence of $u$ at the node $p$, refer to [18].

At this stage, the parameters $\lambda_{lf}$ and $\lambda_{rf}$ are still unknown. In order to define them, we rely on conditions (19) for the positivity of the star values of the specific volumes. However, we would also aim to impose the positivity of the star values of the specific internal energies and entropy preservation. Hence, let us compute $\varepsilon_{sf}^*$ with $s = l, r$, and $E = \varepsilon + \frac{1}{2}\mathbf{u}^2$, $E_{sf}^* = \varepsilon_{sf}^* + \frac{1}{2}((u_{n,sf})^2 + u_{t,sf}^2)$. Using the
left/right jump conditions we obtain,

$$
\begin{align}
\lambda_{rf} \varepsilon_{rf}^* &= \lambda_{rf} \varepsilon_{rf} - \frac{\lambda_{rf}}{2} (u_{n,rf}^* + u_{n,rf}) (u_{n,rf}^* - u_{n,rf}) + \pi_p \hat{u}_{n,pf} - p_{rf} u_{n,rf} \nonumber \\
&\stackrel{(15)}{=} \lambda_{rf} \varepsilon_{rf} - \frac{1}{2} (u_{n,rf}^* + u_{n,rf}) (\pi_p - p_{rf}) + \frac{1}{2} (\pi_p + p_{rf}) (\hat{u}_{n,pf} - u_{n,rf}) + \frac{1}{2} (\pi_p - p_{rf}) (\hat{u}_{n,pf} + u_{n,rf}) \nonumber \\
&= \lambda_{rf} \varepsilon_{rf} + \frac{1}{2} (\pi_p + p_{rf}) (\hat{u}_{n,pf} - u_{n,rf}) + \frac{1}{2} (\pi_p - p_{rf}) (\hat{u}_{n,pf} - u_{n,rf}^*) \nonumber \\
&\stackrel{(12),(15)}{=} \lambda_{rf} \left( \varepsilon_{rf} - p_{rf} (\tau_{rf}^* - \tau_{rf}) + \frac{\lambda_{rf}^2}{2} (\tau_{rf}^* - \tau_{rf})^2 - \frac{1}{2} (\hat{u}_{n,pf} - u_{n,rf}^*)^2 \right) \tag{37}
\end{align}
$$

and similarly for the left star value,

$$
\varepsilon_{lf}^{*} = \varepsilon_{lf} - p_{lf}(\tau_{lf}^{*} - \tau_{lf}) + \frac{\lambda_{lf}^{2}}{2}(\tau_{lf}^{*} - \tau_{lf})^{2} - \frac{1}{2}(\hat{u}_{n_{lf}} - u_{n,lf}^{*})^{2}. \quad (38)
$$

The first three terms on the right-hand sides recall the classic conditions, namely they form a quadratic function with respect to $\tau_{sf}^* - \tau_{sf}$. Hence, it is always greater than the minimum $\varepsilon_{sf} - \frac{p_{sf}^2}{2\lambda_{sf}^2}$, that is to say, if

$$
\lambda_{sf} \geq \frac{p_{sf}}{\sqrt{2\varepsilon_{sf}}}, \qquad (39)
$$

then we have

$$
\varepsilon_{lf} - p_{lf}(\tau_{lf}^{*} - \tau_{lf}) + \frac{\lambda_{lf}^{2}}{2}(\tau_{lf}^{*} - \tau_{lf})^{2} \geq 0.
$$

Conditions (39) imply $\varepsilon_{sf}^* \ge 0$ for the two-point scheme and the multi-point method described in sections 3.2-3.3 as $\hat{u}_{n,pf} = u_{n,sf}^*$. However, this is not generally true for the current solver as we have the additional term $-\frac{1}{2}(\hat{u}_{n,pf} - u_{n,sf}^*)^2$. Since the latter is always negative, it is clear that we obtain a more restrictive condition as we need

$$
\varepsilon_{sf} - p_{sf}(\tau_{sf}^* - \tau_{sf}) + \frac{\lambda_{sf}^2}{2}(\tau_{sf}^* - \tau_{sf})^2 \geq \frac{1}{2}(\hat{u}_{n,pf} - u_{n,sf}^*)^2 \geq 0.
$$

Unfortunately, we are not able to provide an explicit condition on $\lambda_{sf}$. It is clear that the closer $\hat{u}_{n,pf}$ is to both $u_{n,lf}^*$ and $u_{n,rf}^*$, the less restrictive the condition on $\lambda_{sf}$ is. In practice, we impose the following weighted average,

$$
\hat{u}_{n,pf} = \frac{\lambda_{lf} u_{n,lf}^* + \lambda_{rf} u_{n,rf}^*}{\lambda_{lf} + \lambda_{rf}} = u_{n,f}^{*,2P}.
$$

Hence, $\hat{u}_{n,pf}$ does not depend on the nodal pressure $\pi_p$.

**Remark 6 (About the positivity of the specific internal energies at the node).** As previously observed, we are not able to prove the positivity of the star values of the specific internal energies at each subface. However, we are able to establish a link between the nodal solver (33) for the pressure and the sign of the specific internal energies when considering all the RSs around a given node. First of all, given formulas (37)-(38), computations (not shown here) provide a different expression for the star values of the specific internal energies:

$$
\lambda_{rf} \varepsilon_{rf}^* = \lambda_{rf} \varepsilon_{rf} + \frac{1}{2\lambda_{rf}} ((\pi_p)^2 - p_{rf}^2) + \pi_p (\hat{u}_{n,pf} - u_{n,rf}^*),
$$

$$
\lambda_{lf} \varepsilon_{lf}^* = \lambda_{lf} \varepsilon_{lf} + \frac{1}{2\lambda_{lf}} ((\pi_p)^2 - p_{lf}^2) - \pi_p (\hat{u}_{n,pf} - u_{n,lf}^*).
$$
Hence, we get

$$ 
\begin{aligned}
\sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} (\lambda_{lf} \varepsilon_{lf}^* + \lambda_{rf} \varepsilon_{rf}^*) &= \\
\sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} \left( \lambda_{lf} \varepsilon_{lf} + \frac{1}{2\lambda_{lf}} ((\pi_p)^2 - p_{lf}^2) + \lambda_{rf} \varepsilon_{rf} + \frac{1}{2\lambda_{rf}} ((\pi_p)^2 - p_{rf}^2) - \pi_p (u_{n,lf}^* - u_{n,rf}^*) \right) &= (40) \\
\sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} \left( \lambda_{lf} \varepsilon_{lf} + \frac{1}{2\lambda_{lf}} ((\pi_p)^2 - p_{lf}^2) + \lambda_{rf} \varepsilon_{rf} + \frac{1}{2\lambda_{rf}} ((\pi_p)^2 - p_{rf}^2) \right) \ge 0,
\end{aligned}
$$

where we have used

$$ \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} \pi_p (u_{n,lf}^* - u_{n,rf}^*) = \pi_p \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} (u_{n,lf}^* - u_{n,rf}^*) = \pi_p \cdot 0 = 0, $$

and conditions (39). It is clear that result (40) is not completely satisfying but it provides some information about the sign of the averaged value of the specific internal energy around the node. Furthermore, to obtain such a result, we did not exploit the definition of $\hat{u}_{n,pf}$ but only the nodal solver (33) for the pressure.

Let us now continue with the entropy analysis following [22]. We aim to impose $\eta_{sf}^* \ge \eta_{sf}$ across the discontinuities of the RS. We decompose the thermodynamic process $(\tau_{sf}, \eta_{sf}) \to (\tau_{sf}^*, \eta_{sf}^*)$ into the isentropic process $(\tau_{sf}, \eta_{sf}) \to (\tau_{sf}^*, \eta_{sf})$ and the isochoric process $(\tau_{sf}^*, \eta_{sf}) \to (\tau_{sf}^*, \eta_{sf}^*)$. Hence, we obtain

$$ \varepsilon(\tau_{sf}^*, \eta_{sf}^*) - \varepsilon(\tau_{sf}, \eta_{sf}) = \int_{\eta_{sf}}^{\eta_{sf}^*} \theta(\tau_{sf}^*, \eta) d\eta - \int_{\tau_{sf}}^{\tau_{sf}^*} p(\tau, \eta_{sf}) d\tau, \quad (41) $$

see [25] for the details. Inserting formulas (37)-(38) into (41), we get

$$ \int_{\eta_{sf}}^{\eta_{sf}^*} \theta(\tau_{sf}^*, \eta) d\eta = \int_{\tau_{sf}}^{\tau_{sf}^*} (p(\tau, \eta_{sf}) - p(\tau_{sf}, \eta_{sf})) d\tau + \frac{\lambda_{sf}^2}{2} (\tau_{sf}^* - \tau_{sf})^2 - \frac{1}{2} (\hat{u}_{n,f} - u_{n,sf}^*)^2 $$

Once again, the first term on the right-hand side is what one can obtain for the two-point scheme or multi-point one with nodal velocity, whereas the last term on the right-hand side is clearly related to nodal solver for the pressure. Once again, the negative terms $-\frac{1}{2}(\hat{u}_{n,f} - u_{n,sf}^*)^2$ appear and make the resulting condition clearly more restrictive. Indeed, we need to impose

$$ \lambda_{sf}^2 \ge -\frac{2}{(\tau_{sf}^* - \tau_{sf})^2} \int_{\tau_{sf}}^{\tau_{sf}^*} (p(\tau, \eta_{sf}) - p(\tau_{sf}, \eta_{sf})) d\tau + \frac{1}{(\tau_{sf}^* - \tau_{sf})^2} (\hat{u}_{n,f} - u_{n,sf}^*)^2. $$

In conclusion, on the one hand we are able to ask for the positivity of the star values of the specific volumes by imposing the same conditions for the three RSs. On the other hand, the RS with nodal pressure presents more restrictive conditions for both the positivity of the intermediate values of the specific internal energies and for entropy-stability.

**Remark 7 (Consistency and conservation properties).** This remark is only to highlight that, even if there is a nodal parameter in the RS, we still have a scheme which is conservative in a more classic way, namely at each subface. Indeed, since the usual face-based jump condition $\hat{\mathbf{L}}_{n,lf} = \hat{\mathbf{L}}_{n,rf}$ across the middle wave is still valid, also the associated consistency conditions (22) are verified at each subface. Hence, contrarily to the previous RS with nodal velocity, the associated numerical scheme with nodal pressure is conservative at each subface and, a fortiori, at each node. Observe that this is possible thanks to the distinction between $\hat{\mathbf{L}}_{n,sf}$ and $\mathbf{L}_{n,sf}^*$.

**Remark 8 (Going in a difficult direction).** As remarked at the beginning of this section, when developing the RS with nodal pressure, we immediately see that having both $\hat{\mathbf{L}}_{n,sf} = \mathbf{L}_{n,sf}^*$ and $p_{lf}^* = p_{rf}^* = p_f^*$ leads to the two-point RS of section 3.2. However, this is true only if we ask for the face-based jump condition across the middle wave for the mass equation, namely condition (17). However, to proceed in this direction would lead to a RS which is not consistent and thus, a scheme which is not conservative. Indeed, if we do not impose condition (17) but $\hat{\mathbf{L}}_{n,sf} = \mathbf{L}_{n,sf}^*$ and $p_{lf}^* = p_{rf}^* = p_f^*$,
the balance term $\mathbf{B}_{n,f}^0$ reads

$$ \mathbf{B}_{n,f}^0 = -\mathbf{U}_{n,rf}^*(\Lambda_{n,0f} - \mathbf{u}_{n,rf}^*) + \mathbf{U}_{n,lf}^*(\Lambda_{n,0f} - \mathbf{u}_{n,lf}^*) + \mathbf{L}_{n,rf}^* - \mathbf{L}_{n,lf}^* $$

When there is a unique star velocity, the term $-\mathbf{U}_{n,rf}^*(\Lambda_{n,0f} - \mathbf{u}_{n,rf}^*) + \mathbf{U}_{n,lf}^*(\Lambda_{n,0f} - \mathbf{u}_{n,lf}^*)$ disappears by simply imposing $\Lambda_{n,0f} = \mathbf{u}_{n,pf}^*$. However, in the particular case of a RS with nodal pressure, this is not true anymore. This also means that neither the face-based nor the node-based consistency conditions are satisfied, the latter being

$$ \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} \Lambda_{n,0f} (\mathbf{U}_{n,rf}^* - \mathbf{U}_{n,lf}^*) = \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} (\mathbf{u}_{n,pf,rf}^* \mathbf{U}_{n,rf}^* - \mathbf{u}_{n,pf,lf}^* \mathbf{U}_{n,lf}^*) $$

Since it did not seem trivial to impose such relations, it appeared that the previous assumptions were far too restrictive, hence the idea of modifying the definition of the intermediate fluxes.

# 4 Numerical schemes and properties

In this section, we describe and analyze the numerical schemes associated to the three different RSs we described in the previous section 3. For the sake of simplicity, we use the following abbreviations:

* 2P: two-point scheme based on the RS of section 3.2 without nodal parameters;

* MPu: multi-point scheme with nodal velocity based on the RS of section 3.3;

* MPp: multi-point scheme with nodal pressure based on the RS of section 3.4.

## 4.1 Finite volume formulation and flux definition

In this section, we describe the FV formulation together with the flux definitions for the three numerical schemes. For details about the 2P and MPu methods, refer respectively to [25] and [22].

We have already specified that the main difference between two-point and multi-point schemes is the number of cells on which the numerical flux depends. On the one hand, the two-point method is classic in the sense that the flux only depends on the left and right states to the given face (of course, at first-order of accuracy). On the other hand, the multi-point RS depends on all the cells around the given subface, whose exact number depends on the specific mesh, refer also to figure 2. To this end, either the nodal velocity $\mathbf{u}_p$ or the nodal pressure $\pi_p$ have been included inside the RS, and will therefore be present in the flux. In 2D, each face of the mesh has two nodes and it can be decomposed into two subfaces. Hence, we define two fluxes for each face, namely subface fluxes. In this sense, the classic face-based FV formulation has to be modified accordingly in order to take into account the subface fluxes. To obtain the relevant formulation, we simply need to integrate the Euler system in time and space [22]:

$$ \mathbf{U}_c^{n+1} = \mathbf{U}_c^n - \frac{\Delta t}{|\omega_c|} \sum_{p \in \mathcal{P}(c)} \sum_{f \in \mathcal{S}\mathcal{F}(pc)} l_{pcf} \bar{\mathbf{F}}_{pcf}, \quad (42) $$

where $\mathbf{U}_c(t) = \frac{1}{|\omega_c|} \int_{\omega_c} \mathbf{U}(\mathbf{x}, t) dv$ is the cell-averaged value of $\mathbf{U}$ over the cell $\omega_c$ and $\mathbf{U}_c^n \equiv \mathbf{U}_c(t^n)$ for any discrete time $t^n$. Then, $\bar{\mathbf{F}}_{pcf}$ is the left-sided (subface) flux with respect to the subface $f$ impinging at node $p$ and cell $c$. If the subface flux does not depend on the nodes of the mesh, clearly formulation (42) is equivalent to the classic face-based one. Given the unit outward normal $\mathbf{n}_{pcf}$, the left- and right-sided fluxes between cells $c$ and $d \equiv d(c, f)$ in the normal direction $\mathbf{n}_{pcf}$ associated to the RSs of section 3 are respectively given by

$$ \bar{\mathbf{F}}_{n_{pcf}} = \mathbf{F}_{n_{pcf}}(\mathbf{U}_c) - (\Lambda_{n,lf}^{(-)}(\mathbf{U}_{n,lf}^* - \mathbf{U}_c) + \Lambda_{n,0f}^{(-)}(\mathbf{U}_{n,rf}^* - \mathbf{U}_{n,lf}^*) + \Lambda_{n,rf}^{(-)}(\mathbf{U}_d - \mathbf{U}_{n,rf}^*)) \quad (43) $$

and

$$ \bar{\mathbf{F}}_{n_{pcf}}^{+} = \mathbf{F}_{n_{pcf}}(\mathbf{U}_d) - (\Lambda_{n,lf}^{(+)}(\mathbf{U}_{n,lf}^* - \mathbf{U}_c) + \Lambda_{n,0f}^{(+)}(\mathbf{U}_{n,rf}^* - \mathbf{U}_{n,lf}^*) + \Lambda_{n,rf}^{(+)}(\mathbf{U}_d - \mathbf{U}_{n,rf}^*)) \quad (44) $$

where we recall $\Lambda^{(-)} = \frac{1}{2}(|\Lambda| - \Lambda)$ and $\Lambda^{(+)} = \frac{1}{2}(|\Lambda| + \Lambda)$. These formulas for $\bar{\mathbf{F}}_{n_{pcf}}$, $\bar{\mathbf{F}}_{n_{pcf}}^+$ are generally valid for a (subface) flux based on a RS with three waves of speed $\Lambda_{n,lf}$, $\Lambda_{n,0f}$, $\Lambda_{n,rf}$. Hence, we may
use these formulas for all the three methods 2P, MPu and MPp. Furthermore, formulas (43)-(44) do not present any difference from the classic ones [32], apart from being related to subfaces and not faces. Of course, these formulas can also be generalized to any number of discontinuities [32, 22]. We can also define the averaged flux on the subface f by taking the average of $\bar{\mathbf{F}}_{\mathbf{n}_{pcf}}$ and $\bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^+$,

$$ 
\begin{aligned}
\bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{AV} &= \frac{1}{2} (\bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{-} + \bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{+}) \\
&= \frac{1}{2} (\mathbf{F}_{\mathbf{n}_{pcf}}(\mathbf{U}_c) + \mathbf{F}_{\mathbf{n}_{pcf}}(\mathbf{U}_d)) \\
&\quad - \frac{1}{2} (|\Lambda_{\mathbf{n},lf}|(\mathbf{U}_{\mathbf{n},lf}^{*} - \mathbf{U}_c) + |\Lambda_{\mathbf{n},0f}|(\mathbf{U}_{\mathbf{n},rf}^{*} - \mathbf{U}_{\mathbf{n},lf}^{*}) + |\Lambda_{\mathbf{n},rf}|(\mathbf{U}_d - \mathbf{U}_{\mathbf{n},rf}^{*})) \\
&= \bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{2P},
\end{aligned} 
\quad (45) $$

which also correspond to the two-point flux. Concerning the MPu method, we can formulate the numerical flux to highlight the difference as follows,

$$ \bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{\mathrm{MPu}} = \bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{\mathrm{AV}} - \frac{1}{2} (p_d^* - p_c^*) \begin{pmatrix} 0 \\ 1 \\ 0 \\ \mathbf{u}_{\mathbf{n}_{pcf}}^* \end{pmatrix} = \bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{\mathrm{AV}} - \frac{\lambda_{lf} + \lambda_{rf}}{2} (\mathbf{u}_p \cdot \mathbf{n}_{pcf} - u_{\mathbf{n}_{pcf}}^{2P}) \begin{pmatrix} 0 \\ 1 \\ 0 \\ \mathbf{u}_{\mathbf{n}_{pcf}}^* \end{pmatrix}, \quad (46) $$

as the difference between the right and left-sided fluxes is given by

$$ \bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{+ , \mathrm{MPu}} - \bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{- , \mathrm{MPu}} = (p_d^* - p_c^*) \begin{pmatrix} 0 \\ 1 \\ 0 \\ \mathbf{u}_{\mathbf{n}_{pcf}}^* \end{pmatrix}. $$

Finally, it only remains to give further details for the MPp flux, which is simpler than the MPu one. Indeed, since $\hat{\mathbf{L}}_{\mathbf{n},f} = (0, \pi_p, 0, \pi_p \hat{u}_{\mathbf{n}_f})^T$, the left- and right-sided fluxes are equal for the MPp method:

$$ \bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{+ , \mathrm{MPp}} - \bar{\mathbf{F}}_{\mathbf{n}_{pcf}}^{- , \mathrm{MPp}} = 0. $$

This also means that the MPp method is actually conservative not only at the nodes but also at the subfaces. Furthermore, it is now clear that the MPp numerical flux is also given by the average formula (45), which implies that the MPp and 2P flux can be formulated in the same way (of course the intermediate values of the RS remain different).

**Remark 9 (Time step condition).** Since the computations are the same as the ones shown in [22], we directly provide the relevant formula for the time step condition:

$$ \Delta t \le \min_c \left[ \frac{|\omega_c|}{\sum_{p \in \mathcal{P}(c)} \sum_{f \in \mathcal{S}\mathcal{F}(pc)} l_{pcf} \max(|\Lambda_{l,pcf}|, |\Lambda_{r,pcf}|)} \right]. \quad (47) $$

**Remark 10 (Positivity and entropy preservation).** As explained in [22, 24], once the RS has been proved to be positivity- and entropy- stable, it is only a matter of computation to show that the associated first-order FV numerical scheme also satisfies such properties under the time step condition (47). Once again, we rely on the previous papers [22, 24] for the details.

## 4.2 Behavior of the schemes at the low Mach limit

One of the last steps before numerical simulations is to prove that the MPp scheme is actually AP, which we do not expect for the 2P and MPu methods. For this purpose, we first briefly reformulate the RSs in non-dimensional form. Then, an asymptotic expansion in power of the Mach number M is performed to show the limit behavior of the methods.
### 4.2.1 Non-dimensional scheme

Let us consider the non-dimensional 2D Euler system (4). The formulas for the non-dimensional RSs and methods are analogous to the ones of the previous sections 3-4. For instance, formulas (15) become

$$ 
\left\{ 
\begin{array}{l} 
\tilde{\lambda}_{lf}(\tilde{u}_{n,lf}^* - \tilde{u}_{n,lf}) + \frac{1}{M^2}(\hat{p}_{lf} - \tilde{p}_{lf}) = 0, \\ 
\tilde{\lambda}_{lf}(\tilde{u}_{t,lf}^* - \tilde{u}_{t,lf}) = 0, \\ 
\tilde{\lambda}_{lf}(\tilde{E}_{lf}^* - \tilde{E}_{lf}) + (\hat{p}_{lf}\tilde{u}_{n,f} - \tilde{p}_{lf}\tilde{u}_{n,lf}) = 0, 
\end{array} 
\right. 
\quad 
\left\{ 
\begin{array}{l} 
\tilde{\lambda}_{rf}(\tilde{u}_{n,rf}^* - \tilde{u}_{n,rf}) - \frac{1}{M^2}(\hat{p}_{rf} - \tilde{p}_{rf}) = 0, \\ 
\tilde{\lambda}_{rf}(\tilde{u}_{t,rf}^* - \tilde{u}_{t,rf}) = 0, \\ 
\tilde{\lambda}_{rf}(\tilde{E}_{rf}^* - \tilde{E}_{rf}) - (\hat{p}_{rf}\tilde{u}_{n,f} - \tilde{p}_{rf}\tilde{u}_{n,rf}) = 0, 
\end{array} 
\right. 
$$

and thus, the star velocity and pressure for the 2P scheme now read

$$ 
\begin{aligned} 
\tilde{u}_{n,f}^{*,2P} &= \frac{\tilde{\lambda}_{lf}\tilde{u}_{n,lf} + \tilde{\lambda}_{rf}\tilde{u}_{n,rf}}{\tilde{\lambda}_{lf} + \tilde{\lambda}_{rf}} - \frac{1}{M^2} \frac{(\tilde{p}_{rf} - \tilde{p}_{lf})}{\tilde{\lambda}_{rf} + \tilde{\lambda}_{lf}}, \\ 
\frac{1}{M^2} \tilde{p}_{f}^{*,2P} &= \frac{1}{M^2} \frac{\tilde{\lambda}_{lf}\tilde{p}_{rf} + \tilde{\lambda}_{rf}\tilde{p}_{lf}}{\tilde{\lambda}_{lf} + \tilde{\lambda}_{rf}} - \frac{\tilde{\lambda}_{rf}\tilde{\lambda}_{lf}}{\tilde{\lambda}_{rf} + \tilde{\lambda}_{lf}} (\tilde{u}_{n,rf} - \tilde{u}_{n,lf}). 
\end{aligned} 
\quad (48) 
$$

Clearly, the nodal velocity and pressure are respectively given by

$$ \sum_{f \in S_f(p)} \tilde{l}_{pf}(\tilde{\lambda}_{lf} + \tilde{\lambda}_{rf})(\tilde{u}_{n,f}^* - \tilde{u}_{n,f}^{*,2P})\mathbf{n}_{pf} = 0, $$

and

$$ \sum_{f \in S_f(p)} \tilde{l}_{pf} \frac{1}{M^2} \left( \frac{1}{\tilde{\lambda}_{lf}} + \frac{1}{\tilde{\lambda}_{rf}} \right) \left( \hat{p}_f - \tilde{p}_f^{*,2P} \right) = 0. \quad (49) $$

Observe that we defined the wave speeds $\tilde{\lambda}_{n,sf}$ and parameters $\tilde{\lambda}_{sf}$, $s=l,r$, as

$$ \tilde{\lambda}_{n,lf} = \tilde{u}_{n,lf} - \frac{\tilde{\lambda}_{lf}}{\tilde{\rho}_{lf}}, \quad \tilde{\lambda}_{n,rf} = \tilde{u}_{n,rf} + \frac{\tilde{\lambda}_{rf}}{\tilde{\rho}_{rf}}, \quad \text{with} \quad \tilde{\lambda} \equiv \frac{\tilde{\rho}\tilde{a}}{M}. $$

Finally, for the non-dimensional numerical schemes, the previous formulas for the fluxes (43)-(44) are still valid. Observe that, in the rest of the section, we drop the symbol $\sim$ for the sake of simplicity.

### 4.2.2 Analysis of the schemes at the low Mach limit

It only remains to analyze the different schemes by applying an asymptotic expansion in powers of the Mach number. We aim to prove that, in the limit $M \to 0$, the MPp scheme is a consistent approximation of the incompressible equations (5) contrarily to the 2P and MPu methods. For this purpose, few definitions are needed.

**Remark 11 (Asymptotic expansion).** Given a variable $X$, its asymptotic expansion in power of $M$ reads:

$$ X = X_0 + M X_1 + M^2 X_2 + \mathcal{O}(M^3). $$

**Definition 4.1 (Well-prepared data).** Consider the vector of variables $(\rho, \mathbf{u}, p)^T(\mathbf{x}, t=0)$ at initial time $t=0$. We define the initial data as *well-prepared* if

$$ \rho(\mathbf{x}, t=0) = \rho_0 + \mathcal{O}(M) \quad \text{with} \quad \rho_0 > 0, $$

$$ \mathbf{u}(\mathbf{x}, t=0) = \mathbf{u}_0 + \mathcal{O}(M) \quad \text{with} \quad \nabla \cdot \mathbf{u}_0 = 0, $$

$$ p(\mathbf{x}, t=0) = P_0 + \mathcal{O}(M^2) \quad \text{with} \quad \nabla P_0 = 0. $$

**Definition 4.2 (Asymptotic-preserving property).** Consider initial well-prepared data $\mathbf{U}(\mathbf{x}, t=0)$. A numerical scheme is asymptotic-preserving when $M$ tends to zero if it is a consistent discretization of the incompressible Euler equations (5). Namely, if $\Delta t$ satisfies the time step condition (47), $\mathbf{U}(\mathbf{x}, \Delta t)$ are still well-prepared data.

It is well-known that classic FV methods based on dimensional splitting and 1D RSs such as the HLLC or Roe ones [20] are not AP in two or more space dimensions. The problem is usually related
to the fact that the RSs are indeed 1D whereas the incompressible equations incorporate truly multi-dimensional operator, i.e. the divergence of **u**. Let us briefly recall the problem by showing that the 2P method is not AP.

**Remark 12.** The two-point scheme (2P) is not asymptotic-preserving.

*Proof.* What prevents the AP property for the 2P scheme is the intermediate value for the pressure (48) in the RS:

$$ \frac{1}{M^2} p_f^{*,2P} = \frac{1}{M^2} \frac{\lambda_{lf} p_{rf} + \lambda_{rf} p_{lf}}{\lambda_{lf} + \lambda_{rf}} - \frac{\lambda_{rf} \lambda_{lf}}{\lambda_{rf} + \lambda_{lf}} (u_{n,rf} - u_{n,lf}) \\ = \frac{P_0}{M^2} + \mathcal{O}(1) + \mathcal{O}\left(\frac{\Delta x_n}{M}\right), $$

as $\lambda \equiv \mathcal{O}\left(\frac{1}{M}\right)$,

$$ u_{n,0rf} - u_{n,0lf} = \Delta x_n \partial_{x_n} u_{n,0f} + \mathcal{O}(\Delta x_n^2) = \mathcal{O}(\Delta x_n), $$

and with a slight abuse of notation for the term $\mathcal{O}\left(\frac{\Delta x_n}{M}\right)$. Hence, the pressure $p_f^{*,2P}$ involves the 1D jump of the velocity in the normal direction. However, the notion of well-prepared data concerns a truly multi-dimensional operator, the divergence of the velocity. Hence, even if we know that $\nabla \cdot \mathbf{u}_0 = 0$, nothing can be said about the value of $\partial_{x_n} u_{n,0f}$. It is then clear that the presence of the term $\mathcal{O}\left(\frac{\Delta x_n}{M}\right)$ makes the convergence when $M \to 0$ impossible unless $\Delta x_n = \mathcal{O}(M)$, which is usually computationally too expensive. We do not include details about the other equations as similar computations can be found in the literature [6]. □

A similar reasoning can be applied to the MPu method.

**Remark 13.** The multi-point scheme with nodal velocity (MPu) is not asymptotic-preserving.

*Proof.* The proof is analogous to the one for the two-point scheme, see remark 12. □

Finally, we are going to explain the reason why the MPp method does not present the same flaw. Clearly, this is related to the presence of the nodal pressure.

**Theorem 1.** The multi-point scheme with nodal pressure (MPp) is asymptotic-preserving.

*Proof.* Here, we propose an argument based on asymptotic expansion in powers of the Mach number to justify the good behavior of the MPp method in low Mach number flows. Once again, we do not include the complete computations as they are quite technical and similar to what can be found in the literature [6]. However, we explain why the presence of the nodal pressure allows to cure the flaw of classic methods. Let us directly consider formula (49):

$$ \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} \frac{1}{M^2} \left( \frac{1}{\lambda_{lf}} + \frac{1}{\lambda_{rf}} \right) \pi_p = \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} \left( \frac{1}{M^2} \left( \frac{p_{rf}}{\lambda_{rf}} + \frac{p_{lf}}{\lambda_{lf}} \right) - (u_{n,rf} - u_{n,lf}) \right) \\ = \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf} \left( \frac{1}{M^2} \left( \frac{p_{rf}}{\lambda_{rf}} + \frac{p_{lf}}{\lambda_{lf}} \right) \right) - |\omega_p| D_p(\mathbf{u}), $$

where we used remark 5. Hence, knowing $D_p(\mathbf{u})$ to be a weak consistent discretization of $\nabla \cdot \mathbf{u}$, at first-order of accuracy in space, we get

$$ \frac{1}{M^2} \pi_p = \frac{P_0}{M^2} + \mathcal{O}(1) + \mathcal{O}\left(\frac{|\omega_p|}{l_p}\right), $$

with $l_p = \sum_{f \in \mathcal{S}\mathcal{F}(p)} l_{pf}$. We immediately remark that pressure $\pi_p$ does not present the typical term of order $\mathcal{O}\left(\frac{\Delta x_n}{M}\right)$ which usually prevents the convergence to the incompressible equations for $M \to 0$.
Then, after few computations, for the other star values we get

$$ 
\begin{aligned}
\frac{1}{\rho_{lf}^*} &= \frac{1}{\rho_{0lf}} + \mathcal{O}(M) \\
u_{n,lf}^* &= u_{n,0lf} + \mathcal{O}(M) + \mathcal{O}\left(M \frac{|\omega_p|}{l_p}\right), \quad \text{and} \\
u_{t,lf}^* &= u_{t,0lf} + \mathcal{O}(M) \\
E_{n,lf}^* &= E_{n,0lf} + \mathcal{O}(M)
\end{aligned}
\quad \frac{1}{\rho_{rf}^*} = \frac{1}{\rho_{0rf}} + \mathcal{O}(M) \\
u_{n,rf}^* &= u_{n,0rf} + \mathcal{O}(M) + \mathcal{O}\left(M \frac{|\omega_p|}{l_p}\right). \\
u_{t,rf}^* &= u_{t,0rf} + \mathcal{O}(M) \\
E_{n,rf}^* &= E_{n,0rf} + \mathcal{O}(M)
\end{aligned}
$$

Hence, we can conclude that the flux formulas are not affected by the usual $\mathcal{O}(\frac{1}{M})$ diffusion flaw.

<figure><img src="image_7.png" /></figure>

**Remark 14.** In the literature, it is common to propose an argument based on asymptotic expansions in powers of Mach number to explain the behavior of a scheme in low Mach number flows, as we have done. However, we stress that it cannot be considered an exhaustive proof of the AP property. Further evidence for the good behavior of the MPp scheme is provided by numerical simulations.

# 5 Numerical results

This section is devoted to the comparison of the numerical results of the different numerical schemes we described so far. More specifically, we consider the following numerical tests.

* Radial Riemann problem: academic test to check the methods are able to correctly reproduce the solution of a Riemann problem;

* Gresho and Kelvin-Helmoltz tests: to check the behavior of the schemes at low Mach number flows;

* Sedov and flow past a half-cylinder problem: to check the behavior of the methods at supersonic and hypersonic flows.

Unless otherwise stated, we use $\gamma = 1.4$, free boundary conditions and CFL = 0.45. Concerning the parameters $\lambda_{sf}$, $s = l, r$, we use an iterative procedure with only two iterations for the multi-point schemes. More specifically, the parameters are initialized with the two-point conditions (25), then at the second iteration we exploit Dukowicz formulas (19). Clearly, for the two-point scheme, we directly use conditions (25) without any iterative procedure.

## 5.1 Radial Sod shock tube

As a first test, we consider the radial Sod shock tube. This academic test case allows us to check the good behaviors of the different methods for Riemann problems. As a domain, we consider a square $[0, 1] \times [0, 1]$ and impose the following Initial Conditions (IC):

$$ (\rho, u, v, p)_{\text{IC}} = \begin{cases} (1, 0, 0, 1) & \text{if } r < 0.25 \\ (0.125, 0, 0, 0.1) & \text{otherwise} \end{cases} $$

For more information about the setup of this test, we refer to [22]. Then, we use $t_{\text{end}} = 0.1$ for the ending time and a Cartesian mesh of either $100 \times 100$ cells or $400 \times 400$ cells. In both left and middle images of figure 5, we show the density solution along the line from point [0.5, 0.5] to point [1, 0.5]. More specifically, in the left image, we compare the three first-order methods 2P, MPu and MPp. From the previous papers [22]-[24], we already knew (and we confirm) that the MPu scheme is more diffusive than the 2P method. This difference may be related to the multidimensional character of the MPu scheme, and more specifically to the additional term in the flux definition (46). In fact, the MPu method appears to be more diffusive than the MPp scheme as well. The difference in diffusion between the MPp and 2P methods is less evident, even if the 2P scheme seems slightly less diffusive than the MPp one. Then, in the middle image, we insert the MPp solution with both meshes of $100 \times 100$ cells and $400 \times 400$ cells in order to verify the correct convergence behavior of the method. Finally, on the right, we insert the scatter plot for the MPp density solution with respect to the radius, in order to show that the MPp does not lose the cylindrical symmetry. Similar results can be obtained with the other two methods, thus we do not insert them here. The reference solution is always obtained by solving the radial equations with a first order method on a grid of 40000 points.
<figure><img src="image_8.png" /><figcaption>Figure 5: Radial Sod shock tube. On the left, comparison of 2P (dashed green line), MPu (dashed magenta line) and MPp (black line) density solutions. In the middle, MPp density solution computed with 100 × 100 (black line) and 400 × 400 (blue line) cells. On the right, scatter plot for MPp density solution. Reference solution in red.</figcaption></figure>

## 5.2 Gresho vortex

Next, we consider two numerical tests to analyze the behavior of the three numerical methods in the low Mach number regime. We start with the Gresho problem, a stationary vortex characterized by the Mach number. Given a square $[-0.5, 0.5] \times [-0.5, 0.5]$, the stationary (and initial) solution is given by

$$ (\rho, u, v, p)_{\text{IC}} = \begin{cases} (1, -5, 5, p_0 + 12.5r^2) & \text{if } r < 0.2 \\ (1, 5 - \frac{2}{r}, \frac{2}{r} - 5, p_0 + 4 + 12.5r^2 - 20r + 4\log(5r)) & \text{if } 0.2 \le r < 0.4 \\ (1, 0, 0, p_0 - 2 + 4\log(2)) & \text{if } 0.4 \le r \end{cases} $$

where $p_0 = \frac{1}{\gamma M^2}$ and $r$ is the radius. We consider the reference Mach number $M = 10^{-2}$, $M = 10^{-3}$ and $M = 10^{-4}$, a Cartesian grid of 128 × 128 cells, and show the solutions at time $t_{\text{end}} = 0.01$ for the three methods 2P, MPu and MPp in figure 6. As expected, only the MPp method is capable of preserving the stationary vortex, while the other two methods produce very diffusive outputs even if the ending time is very small. Of course, if we run the MPp for longer times, e.g. $t_{\text{end}} = 1$, we may observe further diffusion due to the advective term approximation. This is confirmed by figure (7), in which we insert the scatter plots of the Mach number for $M = 10^{-2}$, $M = 10^{-3}$ and $M = 10^{-4}$ computed at initial time and at ending time $t_{\text{end}} = 1$ with the MPp method. We observe that the Mach number decreases of a factor of order 10 from one image to the other, as expected, but the diffusion behavior remains the same. This means that it is not due to the acoustic effects but to the advective ones. Finally, the previous results are further confirmed by the analysis of the kinetic energy, which is theoretically conserved in the incompressible limit. Hence, it is a valuable tool to understand the behavior of the schemes in low Mach number flows. For this reason, we kept track of the ratio of the total kinetic energy at some time steps and the initial one:

$$ \frac{E_{\text{kin}}(t)}{E_{\text{kin}}(t=0)} = \frac{\sum_c |\omega_c| \frac{1}{2} \rho_c(t) \|\mathbf{u}_c(t)\|^2}{\sum_c |\omega_c| \frac{1}{2} \rho_c(0) \|\mathbf{u}_c(0)\|^2} $$

In figure 8, we observe that the MPp scheme better preserves the total kinetic energy than the other two methods as expected. Furthermore, the MPu output appears to be extremely poor in this regard.

## 5.3 Kelvin-Helmholtz instability

Let us consider the second and last test case in the low Mach number regime: the Kelvin-Helmholtz instability [33]. We take a rectangular domain $[0, 2] \times [-\frac{1}{2}, \frac{1}{2}]$ and discretize it with a Cartesian mesh of 128 × 64 cells. The initial condition is as follows, $(\rho, u, v, p)_{\text{IC}} = (\gamma + rH(y), M_aH(y), \delta M_a \sin(2\pi x), 1)$, where

$$ H(y) = \begin{cases} -\sin\left(\pi\frac{y+\frac{1}{4}}{\omega}\right) & \text{if } -\frac{1}{4} - \frac{\omega}{2} \le y < -\frac{1}{4} + \frac{\omega}{2} \\ -1 & \text{if } -\frac{1}{4} + \frac{\omega}{2} \le y < \frac{1}{4} - \frac{\omega}{2} \\ \sin\left(\pi\frac{y-\frac{1}{4}}{\omega}\right) & \text{if } \frac{1}{4} - \frac{\omega}{2} \le y < \frac{1}{4} + \frac{\omega}{2} \\ 1 & \text{otherwise} \end{cases} $$
<figure><img src="image_9.png" /><figcaption>Figure 6: Gresho test. 2D magnitude velocity solution for $M = 10^{-2}$ (up), $M = 10^{-3}$ (middle) and $M = 10^{-4}$ (bottom). Initial conditions (left), 2P (middle-left), MPu (middle-right) and MPp (right) solutions.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Gresho test. Scatter plot of Mach number at initial time (red symbol) and at ending time $t_{end} = 1$ (blue symbol) with MPp method. Simulations with reference Mach number $M = 10^{-2}$ (left), $M = 10^{-3}$ (middle) and $M = 10^{-4}$ (right).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8: Gresho test. Loss of the total kinetic energy during the simulation until ending time $t_{end} = 0.01$. 2P (blue line), MPp (black line) and MPu (magenta line) methods. Simulations with Mach number $M = 10^{-2}$ (left), $M = 10^{-3}$ (middle) and $M = 10^{-4}$ (right).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9: Kelvin-Helmholtz test case. 2D density with Mach number $M = 10^{-2}$ (up) and $M = 10^{-3}$ (bottom). 2P (left), MPu (middle) and MPp (right) solutions.</figcaption></figure>

with $\omega = \frac{1}{16}$, $\delta = 0.1$ and $r = M$. We compute the solution at time $t_{end} = 0.8/M$ with both Mach number $M = 10^{-2}$ and $M = 10^{-3}$, see figures 9. In both cases, we confirm what has been observed in the Gresho test case. The 2P and MPu methods are very diffusive and do not give any information about the solution, contrarily to the MPp method which shows some vortices. However, in the MPp outputs, we can also observe some typical numerical artifacts which are known to affect numerical solutions of Kelvin-Helmholtz instabilities [33].

## 5.4 Flow past a cylinder

At this stage, we have verified that the MPp behaves very well in low Mach number flows contrarily to the 2P and MPu schemes. Next, even if it is not the purpose of this paper, we propose a couple of numerical simulations in supersonic flows starting with the flow past a cylinder. This test is very well-known in the literature as it presents the infamous Carbuncle instability [34, 35]. It is believed that not having sufficient dissipation within the scheme leads to such an instability. On the one hand, in [22], it has been shown that the MPu scheme appears to be Carbuncle-free even though it is based on a complete RS and it is capable of preserving contact discontinuities. On the other hand, the 2P scheme presents the Carbuncle instability as it is very similar to the HLLC method. Hence, we aim to verify if the MPp scheme presents such a flaw or not.

For the description of this problem, we refer to [22]. This test is run at Mach number $M = 20$: the flow regime is hypersonic. In the domain, we insert a cylinder of radius 0.5. The initial conditions are
<figure><img src="image_13.png" /><figcaption>Figure 10: Flow past a cylinder test case. 2P (left), MPu (middle) and MPp (right) density solutions. 10 contour lines for the density.</figcaption></figure>

given by $(\rho, u, v, p)_{IC} = (1, 20\sqrt{\gamma}, 0, 1)$. Regarding the boundary conditions, we impose a wall-type one at the cylinder surface, while free boundary conditions are considered at the bottom and up-side of the domain. Then, at the left boundary, we impose an inflow condition: $h_{inf} = 1$, $u_{inf} = M\sqrt{\gamma}$ with reference Mach number $M = 20$. We consider a structured polar mesh of 12201 cells. Since the numerical solution should converge to a steady state, the simulation is stopped if the residuals are smaller than a tolerance (here, $10^{-6}$), otherwise we use $t = 5$ as ending time. On the one hand, the MPu solution converges towards the steady one and, in less than $t = 2$, the simulation is stopped as the residuals become smaller than the tolerance. On the other hand, both the 2P and the MPp schemes do not converge. Indeed, in figure 10, we can see that some instabilities are present in the 2P and MPp solutions contrarily to the MPu one. Hence, even if the MPp scheme incorporates multidimensional information, it is not Carbuncle-free. We also highlight that negative intermediate values of the specific internal energy are computed in the RS with the MPp method during the simulation.

## 5.5 Sedov blast wave problem

The last test case consists of a point explosion which gives rise to a diverging shock wave: the Sedov blast wave problem [22]. The domain is a square $[-1.2, 1.2] \times [-1.2, 1.2]$ for which we use a Cartesian mesh of $200 \times 200$ cells. We assume the domain to be filled with a perfect gas at rest and the point explosion to be at the center of the domain, namely

$$ (\rho, u, v, p)_{IC} = \begin{cases} (1, 0, 0, 0.397056) & \text{if } r < 0.01 \\ (1, 0, 0, 10^{-6}) & \text{otherwise.} \end{cases} $$

With this test, we further confirm what has previously been observed in the half-cylinder test: the MPp scheme is not Carbuncle-free. Indeed, in figure 11, we can observe that the MPu solution does not present any flaw as opposed to the other two schemes. However, the spurious effects are more evident for the 2P method than for the MPp ones. Furthermore, once again the MPp scheme presents negative intermediate values of the RS for the specific internal energy.

## 6 Concluding remarks and perspectives

In this paper, we have proposed an extension of the multi-point scheme with nodal pressure for the linear acoustic equations [18] to the Euler system. The resulting scheme is AP and thus behaves excellently in the low Mach number regime. It also preserves the positivity of the density by imposing the same conditions on the wave speeds of the 2P and MPu methods. However, conditions for the positivity of the specific internal energies and entropy stability should be more restrictive than the classic ones.
<figure><img src="image_14.png" /><figcaption>Figure 11: Sedov test case. 2D density for 2P (left), MPu (middle) and MPp (right) schemes.</figcaption></figure>

Furthermore, even if this scheme is AP, it is not Carbuncle-free contrarily to the multi-point scheme with nodal velocity [22]. It is a natural question whether it is possible to merge the two methods in order to have an AP numerical scheme which is also Carbuncle-free. Of course, a simple possibility would be to insert a detector for the Mach number in order to switch from a numerical scheme to the other when needed. Such an option is currently being investigated.

The scheme we proposed is first-order accurate but second-order of accuracy can easily be achieved by following papers [22, 24]. However, for higher orders of accuracy, the same difficulties are present and are the subject of current research.

Finally, since in low Mach number flows the time step is limited by the sound speed, we also envisage to consider either implicit or implicit-explicit approximations to fasten numerical simulations.

## References

[1] Andrea Thomann, Markus Zenk, Gabriella Puppo, and Christian Klingenberg. An all speed second order IMEX relaxation scheme for the Euler equations. *Communications in Computational Physics*, 28:591–620, 06 2020.

[2] Sandro Malusà and Alessandro Alaia. A well-balanced all-mach scheme for compressible two-phase flow. *Computer Physics Communications*, 299:109131, 2024.

[3] Hervé Guillard and Cécile Viozat. On the behaviour of upwind schemes in the low Mach number limit. *Computers & Fluids*, 28(1):63–86, 1999.

[4] Hervé Guillard and Angelo Murrone. On the behavior of upwind schemes in the low Mach number limit: II. Godunov type schemes. *Computers & Fluids*, 33(4):655–675, 2004.

[5] G. Volpe. Performance of compressible flow codes at low mach numbers. *AIAA Journal*, 31(1):49–56, 1993.

[6] Stéphane Dellacherie. Analysis of godunov type schemes applied to the compressible euler system at low mach number. *Journal of Computational Physics*, 229(4):978–1016, 2010.

[7] Giacomo Dimarco, Raphaël Loubère, and Marie-Hélène Vignal. Study of a new asymptotic preserving scheme for the euler system in the low mach number limit. *SIAM Journal on Scientific Computing*, 39(5):A2099–A2128, 2017.

[8] Andrea Thomann, Gabriella Puppo, and Christian Klingenberg. An all speed second order well-balanced imex relaxation scheme for the euler equations with gravity. *Journal of Computational Physics*, 420:109723, 2020.

[9] Miczek, F., Röpke, F. K., and Edelmann, P. V. F. New numerical solver for flows at various mach numbers. *A&A*, 576:A50, 2015.

[10] Christophe Chalons, Mathieu Girardin, and Samuel Kokh. An all-regime Lagrange-projection like scheme for the gas dynamics equations on unstructured meshes. *Communications in Computational Physics*, 20, 06 2014.

[11] Christophe Chalons, Mathieu Girardin, and Samuel Kokh. An all-regime Lagrange-Projection like scheme for the gas dynamics equations on unstructured meshes. *Communications in Computational Physics*, 20, 06 2014.
scheme for 2D homogeneous models for two-phase flows on unstructured meshes. *Journal of Computational Physics*, 335:885-904, April 2017.

[12] Ben Thornber and Dimitris Drikakis. Numerical dissipation of upwind schemes in low mach flow. *International Journal for Numerical Methods in Fluids*, 56:1535–1541, 03 2008.

[13] Felix Rieper. A low-mach number fix for roe's approximate riemann solver. *Journal of Computational Physics*, 230:5263–5287, 06 2011.

[14] Emanuela Abbate, Angelo Iollo, and Gabriella Puppo. An all-speed relaxation scheme for gases and compressible materials. *Journal of Computational Physics*, 351:1–24, 2017.

[15] Christophe Berthon, Christian Klingenberg, and Markus Zenk. An all Mach number relaxation upwind scheme. *The SMAI Journal of computational mathematics*, 6:1–31, 2020.

[16] Paola Allegrini. Analyse et développement de schémas volumes finis asymptotiquement préservants dans la limite bas-Mach pour les équations d'Euler et de Navier-Stokes. PhD thesis, 2023. Thèse de doctorat dirigée par Vignal, Marie-Hélène Mathématiques et Applications Toulouse 3 2023.

[17] Wasilij Barsukow. Stationarity preserving schemes for multi-dimensional linear systems. *Math. Comput.*, 88:1621–1645, 2018.

[18] Wasilij Barsukow, Raphaël Loubère, and Pierre-Henri Maire. A node-conservative vorticity-preserving finite volume method for linear acoustics on unstructured grids. 2023.

[19] Wasilij Barsukow. Truly multi-dimensional all-speed schemes for the euler equations on cartesian grids. *Journal of Computational Physics*, 435:110216, 2021.

[20] Eleuterio F. Toro. *Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction*. Springer, third edition, 01 2009.

[21] Barsukow, Wasilij and Klingenberg, Christian. Exact solution and the multidimensional Godunov scheme for the acoustic equations. *ESAIM: M2AN*, 56(1):317–347, 2022.

[22] Gérard Gallice, Agnes Chan, Raphaël Loubère, and Pierre-Henri Maire. Entropy stable and positivity preserving Godunov-type schemes for multidimensional hyperbolic systems on unstructured grid. *Journal of Computational Physics*, 468:111493, 2022.

[23] Agnès Chan. *Innovative numerical schemes for 3D supersonic aerodynamics on unstructured mesh*. PhD thesis, 2022. Thèse de doctorat dirigée par Loubère, Raphaël et Maire, Pierre-Henri Mathématiques appliquées et calcul scientifique Bordeaux 2022.

[24] Alessia Del Grosso, Manuel Castro, Agnes Chan, Gérard Gallice, Raphaël Loubère, and Pierre-Henri Maire. A well-balanced, positive, entropy-stable, and multi-dimensional-aware finite volume scheme for 2d shallow-water equations with unstructured grids, 01 2023.

[25] Agnes Chan, Gérard Gallice, Raphaël Loubère, and Pierre-Henri Maire. Positivity preserving and entropy consistent approximate Riemann solvers dedicated to the high-order MOOD-based finite volume discretization of Lagrangian and Eulerian gas dynamics. *Computers & Fluids*, 229:105056, 2021.

[26] Gérard Gallice. *Approximation numérique de Systèmes Hyperboliques Non-linéaires Conservatifs ou Non-conservatifs*. Habilitation à diriger des recherches, Bordeaux University, 2002.

[27] John K Dukowicz. A general, non-iterative riemann solver for godunov's method. *Journal of Computational Physics*, 61(1):119–137, 1985.

[28] Eleuterio F. Toro, Michelle Spruce, and William Speares. Restoration of the contact surface in the HLL-Riemann solver. *Shock Waves*, 4(1):25–34, 1994.

[29] Pierre-Henri Maire. A high-order cell-centered lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes. *Journal of Computational Physics*, pages 2391–2425, 04 2009.

[30] Zhijun Shen, Wei Yan, and Guangwei Yuan. A robust and contact resolving riemann solver on unstructured mesh, part ii, ale method. *Journal of Computational Physics*, 268:456–484, 2014.

[31] Alessia Del Grosso. *Numerical simulation of geophysical flows using high-order and well-balanced Lagrange-Projection methods*. PhD thesis, 2022. Thèse de doctorat dirigée par Chalons, Christophe Mathématiques appliquées université Paris-Saclay 2022.

[32] François Bouchut. *Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws*. Birkhauser, first edition, 2004.

[33] Giovanni Leidi, Robert Andrassy, Wasilij Barsukow, Johann Higl, Philipp Edelmann, and Friedrich K. Röpke. Performance of high-order godunov-type methods in simulations of astrophysical low mach number flows. *Astronomy & Astrophysics*, 686, 03 2024.

[34] Nico Fleischmann, Stefan Adami, Xiangyu Y. Hu, and Nikolaus A. Adams. A low dissipation method to cure the grid-aligned shock instability. *Journal of Computational Physics*, 401:109004, 2020.

[35] Alexander Rodionov. Artificial viscosity in Godunov-type schemes to cure the carbuncle phenomenon. *Journal of Computational Physics*, 345, 05 2017.