# On general and complete multidimensional Riemann solvers for nonlinear systems of hyperbolic conservation laws

Elena Gaburro, Mario Ricchiuto, Michael Dumbser

► To cite this version:

Elena Gaburro, Mario Ricchiuto, Michael Dumbser. On general and complete multidimensional Riemann solvers for nonlinear systems of hyperbolic conservation laws. 2025. hal-05124553

HAL Id: hal-05124553

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

Preprint submitted on 22 Jun 2025

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
# On general and complete multidimensional Riemann solvers for nonlinear systems of hyperbolic conservation laws

Elena Gaburro<sup>a,*</sup>, Mario Ricchiuto<sup>b</sup>, Michael Dumbser<sup>c</sup>

<sup>a</sup>Department of Computer Science, University of Verona, Strada le Grazie 15, Verona, 37134, Italy

<sup>b</sup>Inria, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, 200 Avenue de la Vieille Tour, 33405 Talence, France

<sup>c</sup>Laboratory of Applied Mathematics, DICAM, University of Trento, via Mesiano 77, 38123 Trento, Italy

## Abstract

In this work, we introduce a framework to design multidimensional Riemann solvers for nonlinear systems of hyperbolic conservation laws on general unstructured polygonal Voronoi-like tessellations. In this framework we propose two simple but complete solvers. The first method is a direct extension of the Osher-Solomon Riemann solver to multiple space dimensions. Here, the multidimensional numerical dissipation is obtained by integrating the absolute value of the flux Jacobians over a dual triangular mesh around each node of the primal polygonal grid. The required nodal gradient is then evaluated on a local computational simplex involving the $d+1$ neighbors meeting at each corner. The second method is a genuinely multidimensional upwind flux. By introducing a fluctuation form of finite volume methods with corner fluxes, we show an equivalence with residual distribution schemes (RD). This naturally allows to construct genuinely multidimensional upwind corner fluxes starting from existing and well-known RD fluctuations. In particular, we explore the use of corner fluxes obtained from the so-called N scheme, i.e. the Roe's original optimal multidimensional upwind advection scheme.

Both methods use the full eigenstructure of the underlying hyperbolic system and are therefore complete by construction. A simple higher order extension up to fourth order in space and time is then introduced at the aid of a CWENO reconstruction in space and an ADER approach in time, leading to a family of high order accurate fully-discrete one-step schemes based on genuinely multidimensional Riemann solvers.

We present applications of our new numerical schemes to several test problems for the compressible Euler equations of gas-dyanamics. In particular, we show that the proposed schemes are at the same time carbuncle-free and preserve certain stationary shear waves exactly.

**Keywords:** finite volume (FV) and residual distribution (RD) schemes for hyperbolic conservation laws, genuinely multidimensional Riemann solvers, multidimensional Osher scheme, N scheme applied in the context of FV schemes, unstructured polygonal meshes, compressible Euler equations,

## 1. Introduction

The groundbreaking work of Godunov [51] on finite volume schemes in 1959 allowed to solve non-linear hyperbolic systems of PDE with discontinuous solutions numerically making for the first time use of the exact solution of the Riemann problem (RP) between two adjacent states at a cell interface in order to obtain a stable numerical flux for systems of nonlinear hyperbolic conservation laws via nonlinear upwinding. Later, Godunov-type finite volume schemes were developed that allowed also the use of only

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: elena.gaburro@univr.it (Elena Gaburro), mario.ricchiuto@inria.fr (Mario Ricchiuto), michael.dumbser@unitn.it (Michael Dumbser)</footnote>
approximate solutions of the RP to get the numerical fluxes, so called approximate Riemann solvers, see e.g. [83, 53, 36, 73, 99, 98]. Since then, several attempts have been made to extend the idea of nonlinear upwinding via the exact or approximate solution of Riemann problems also to the multidimensional case. First, developments based on genuinely multidimensional Riemann solvers have been made by Roe *et al.* in [85, 88], Colella [26] and Saltzman [89] and later also in the context of finite volume evolution Galerkin methods [66]. Some thorough studies of the solution of two-dimensional Riemann problems have been provided in the work of Schulz-Rinne *et al.* in [91, 92].

A significant step forward in the development of genuinely multidimensional Riemann solvers for high order Godunov-type finite volume schemes was made by Balsara *et al.* in a series of fundamental papers [14, 15, 35, 16, 18, 17], where a new family of genuinely multidimensional Riemann solvers of the HLL-type has been developed, including also multidimensional HLLC-type solvers. In all the aforementioned work the multidimensional Riemann solvers were not only interesting from a mathematical and conceptual point of view, but also a necessary key building block for the construction of exactly divergence-free numerical schemes for magneto-hydrodynamics (MHD). Along these lines, later, multidimensional Riemann solvers have also been successfully employed in the context of semi-implicit schemes for divergence-free MHD, see e.g. [37, 45, 103]. We also would like to mention a very recent contribution made in [49, 33], porting, for the first time ever, commonly used nodal solvers from the context of Lagrangian hydrodynamics [34, 24, 68, 69, 67, 71, 70], where the use of multidimensional node solvers is mandatory to guarantee geometric compatibility on moving unstructured meshes, to the Eulerian description of compressible gas dynamics on fixed grids.

Parallely to the efforts to formulate genuinely multidimensional Riemann solvers, there have been many works over the years dedicated to devise formulations which do not directly rely on the Riemann Problem, and may allow more naturally to go beyond the use of one dimensional wave propagation ideas. These works have led first to the design of the schemes now known as residual distribution (RD) [8]. These methods started as generalizations to multiple dimensions P.L. Roe's upwind splitting [83, 84]. Over a decade, joint efforts of groups revolving around P.L. Roe and H. Deconinck have managed to propose numerical methods for the Euler equations combining a multidimensional generalization of Roe's linearization, a decomposition of the steady state equations in scalar components, which allowed a unique and efficient approximation of steady two-dimensional flows [80, 81, 82, 86, 29, 30, 72]. The methods proposed relied on splitting flux integrals to cell nodes. In this new setting upwinding is not anymore related to choosing information between two states, but allows to select as target for the flux update among three states at once on a triangle, four on a tetrahedron or quadrangle, etc. Among these methods, the so called N scheme emerged as the optimal positive coefficient first order scheme: featuring the minimal cross-wind numerical dissipation among positivity preserving first order upwind schemes, and providing exact preservation of steady data for mesh aligned advection [80, 81, 82].

The methods have seen since then several evolutions. Recasting system equations as coupled scalar wave equations has been often replaced with a full matrix generalization only requiring the knowledge of the eigenstructure of the Jacobian matrices of the underlying hyperbolic equations [102]. The requirement of an exact linearization of the latter matrices to evaluate multidimensional upwind distribution, required e.g. for the N scheme, was removed in [28, 78] allowing the use of a more practical direct quadrature of the fluxes, instead of the quadrature of the linearized quasi-linear form. Ad-hoc fixes for unwanted features as Carbuncle's and expansion shocks arising in high speed applications have been also proposed [50, 93]. The appearance of such anomalies is somewhat balanced for the low cross-stream dissipation of the N scheme, which is on the other hand appealing in the approximation of complex flows with contact waves and boundary layers on coarse meshes [60]. Higher order extensions for steady state hyperbolic equations have been proposed in [10, 5], and the high order extension to time dependent problems clarified in [6, 78, 77, 76, 2]. Recently, the RD setting has been shown to allow naturally the satisfaction of several additional conservation constraints [11, 7, 3]. Many of these contributions are summarized and compared in the reviews [32, 74, 9, 8]. These developments have required somewhat steering away from some of the
initial ideas which motivated the design of the methods. This has led to more recent work on a different setting called Active Flux methods (see e.g. the recent lecture [87] and [19, 12, 13]).

In this paper we proceed differently and combine ideas from the multidimensional upwind RD and multidimensional Riemann solver frameworks. This provides a new setting to design finite volume methods with multiple state corner numerical fluxes, and genuinely multidimensional dissipation. Within this setting we propose and test two specific definitions of the corner fluxes. The first is a direct extension of the Osher-Solomon Riemann solver to multiple space dimensions. Here, the multidimensional numerical dissipation is obtained by integrating the absolute value of the flux Jacobians over a dual triangular mesh around each node of the primal polygonal grid. The required nodal gradient is then evaluated on a local computational simplex involving the $d+1$ neighbors meeting at each corner. The second method is a genuinely multidimensional upwind flux, obtained by an equivalence with RD schemes on the mentioned local computational simplex. In particular, the corner flux is defined as the internal normal flux plus a fluctuation obtained from the N scheme.

The rest of the paper is organized as follows. In Section 2 we briefly introduce the overall setting of the paper, including an introduction to the notation and the class of governing PDEs under consideration. Subsequently, in Section 3 we present a new genuinely multidimensional Osher-type Riemann solver, which is a straightforward extension of [44, 43] to the multidimensional case. A peculiar feature of the new multidimensional Osher-type flux proposed in this paper is the fact that it provides an expression for the entire multidimensional numerical flux tensor, and not just a numerical flux projected into a particular normal direction. In Section 5 we present a fluctuation formulation of multidimensional Riemann solvers with corner fluxes, allowing to establish an analogy with the so-called residual distribution (RD) framework This provides a whole new setting to design corner fluxes for finite volume (FV) schemes, embedding a genuinely multidimensional upwind flavor. As an example we will use in the tests numerical fluxes based on a generalization of the N scheme. All multidimensional fluxes presented before can then be naturally extended to higher order in space and time via a CWENO reconstruction in space [63, 64, 65, 27, 39, 47] and a fully-discrete one-step ADER discretization in time [100, 96, 38, 23], see Section 6. Comparative computational results are presented for the different schemes in Section 7. The paper closes with some concluding remarks and an outlook to future research given in Section 8.

## 2. Governing PDEs and domain discretization

We aim at solving general nonlinear hyperbolic systems of conservation laws in multiple space dimensions that read

$$ \partial_t \mathbf{Q} + \nabla \cdot \mathbf{F}(\mathbf{Q}) = 0 \quad \text{on} \quad \Omega \times [0, t_f] \subset \mathbb{R}^d \times \mathbb{R}^+ \qquad (1) $$

with $\mathbf{Q} \in \Omega_Q \subset \mathbb{R}^m$ the vector of $m$ conservative variables, $\mathbf{F} = \mathbf{F}(\mathbf{Q}) = (f_1(\mathbf{Q}), f_2(\mathbf{Q}), \dots, f_d(\mathbf{Q}))$ the flux tensor and $d$ the number of space dimensions. The system can be also written in the following quasi-linear form as

$$ \partial_t \mathbf{Q} + \mathbf{A}(\mathbf{Q}) \cdot \nabla \mathbf{Q} = 0, \qquad (2) $$

with $\mathbf{A}(\mathbf{Q}) = (A_1(\mathbf{Q}), A_2(\mathbf{Q}), \dots, A_d(\mathbf{Q}))$ the Jacobian matrices of the fluxes $A_i(\mathbf{Q}) = \partial f_i(\mathbf{Q}) / \partial \mathbf{Q}$. The above system is assumed to verify all the classical properties of a hyperbolic system, i.e. it is endowed with a (convex) entropy pair, the components of $\mathbf{A}(\mathbf{Q}) \cdot \mathbf{n}$ have a full set of real eigenvalues and linearly independent eigenvectors for all normal vectors $\mathbf{n} \neq 0$. For later use we define, given a vector $\mathbf{n} = (n_1, n_2, \dots, n_d) \in \mathbb{R}^d$ the projected matrix

$$ \mathbf{K}_n := \mathbf{A}(\mathbf{Q}) \cdot \mathbf{n} = \sum_{i=1}^{d} A_i(\mathbf{Q}) n_i. \qquad (3) $$
Since the system is hyperbolic $K_n$ also has a full set of real eigenvalues and linearly independent eigenvectors.

The computational domain $\Omega$ is discretized via a set of Voronoi-like polygons/polyhedra $\Omega_c$. To construct the polygonal/polyhedral tessellation we first cover our domain with a set of generator points which are connected via a Delaunay triangulation in 2D or tetrahedralization in 3D; then, each polygon/polyhedron is constructed around a generator point by connecting the centroids (i.e. barycenters) of the dual triangles/tetrahedra sharing that generator. Since we employ the centroids of the triangles/tetrahedra, and not the circumcenters as in the classical definition of Voronoi elements, we should always refer to our elements as Voronoi-like polygons/polyhedra, but sometimes, to lighten the notation, we will use just the word Voronoi. The use of centroid-based Voronoi-like elements avoids small edges/faces and in particular zero length/area ones, hence it increases the homogeneity and the quality of the obtained tessellation. Moreover, it guarantees that each vertex of the polygons (respectively polyhedra), a part for those on a boundary, is attached to exactly 3 Delaunay triangles (respectively 4 Delaunay tetrahedra).

In what follows, we restrict our notation for the geometry to the two dimensional case $d = 2$ and we will make explicit the three-dimensional $d = 3$ extension only in the multidimensional flux definitions.

A considered generator, which is also a vertex of a Delaunay triangle, is indicated with the letter $c$, and eventually, if needed, with the letters $a, b$. It is depicted with a point symbol in the Figures and its coordinates are indicated with $\mathbf{x}_a$ ($\mathbf{x}_b$). Correspondingly, the Voronoi element around this point is indicated with $\Omega_c$ (respectively $\Omega_a, \Omega_b$).

Two elements $\Omega_a$ and $\Omega_c$ which share the same edge/face $\partial\Omega_{ac} = \Omega_c \cap \Omega_a$ are called neighbors. The set of neighbors of the polygonal element $\Omega_c$ is denoted by $\mathcal{N}_c$ and contains all polygons that share an edge with $\Omega_c$. The unit normal vector pointing from $\Omega_c$ to $\Omega_a$ is denoted by $\hat{\mathbf{n}}_{ac}$, while the vector including the length

$$|\partial\Omega_{ac}| = \int_{\partial\Omega_{ac}} dS$$

is denoted by $\mathbf{n}_{ac} = \hat{\mathbf{n}}_{ac}|\partial\Omega_{ac}|$.

A generic vertex of the Voronoi tessellation is denoted by $p$ (and eventually, if needed, with the letters $r, q$). It is depicted with a star symbol in the Figures and its coordinates are denoted with $\mathbf{x}_p$. These points are indeed the barycenters of the Delaunay triangles: we denote with $T_p$ the Delaunay triangle around each $p$ whose vertexes consequently are the generators of the 3 polygons which are attached to the point $p$. Furthermore, the set of polygonal cells around a point $p$ is denoted by $C_p$ and the set of vertexes $p$ that compose a polygon $\Omega_c$ is denoted by $\mathcal{P}_c$.

A sketch of the notation used in this paper is provided in Figures 1 and 2 where the set of elements attached to node $p$ corresponds to $C_p = \{a, b, c\}$ (which also defines the triangle $T_p$) and the set of points which are vertexes of $\Omega_c$ is $\mathcal{P}_c = \{p, q, r, \dots\}$.

Finally, we define the corner normal pointing from cell $\Omega_c$ to point $p$ as

$$\mathbf{n}_{pc} = \frac{1}{2}(\mathbf{n}_{ac} + \mathbf{n}_{bc}) = -\mathbf{n}_{cp}, \qquad (4)$$

and we refer to Figure 2 for a visual representation.

### 3. First order finite volume scheme using multidimensional corner fluxes

As usual in cell-centered finite volume schemes, the discrete data $Q_c^n$ are assumed to represent the cell averages of the conserved variables $Q$

$$Q_c^n = \frac{1}{|\Omega_c|} \int_{\Omega_c} Q(\mathbf{x}, t^n) d\mathbf{x}, \quad \text{with} \quad |\Omega_c| = \int_{\Omega_c} d\mathbf{x} \qquad (5)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: In this Figure we depict the polygonal Voronoi-like cell $\Omega_c$ and, on the right, we highlight the edge $\partial\Omega_{ac}$ between the neighboring Voronoi cells $\Omega_c$ (colored in light blue) and $\Omega_a$ (colored in blue). Note that the vertexes of the Delaunay mesh, i.e. the Voronoi generators $a, b, c$, are depicted with point symbols and the vertexes of the Voronoi $p, q, r$, i.e. the barycenters of the dual triangles, are depicted with star symbols.</figcaption></figure>

which we discretize on the primal Voronoi tessellation made by the polygons $\Omega_c$. Moreover, a classical first order finite volume scheme for (1) reads

$$Q_c^{n+1} = Q_c^n - \frac{\Delta t}{|\Omega_c|} \sum_{\Omega_a \in \mathcal{N}_c} |\partial \Omega_{ac}| \hat{\mathbf{f}}_{ac}, \quad (6)$$

where, in addition to the previously defined notation, $\hat{\mathbf{f}}_{ac}$ is the numerical flux in the normal direction $\hat{\mathbf{n}}_{ac}$ that approximates the integral

$$\hat{\mathbf{f}}_{ac} \approx \frac{1}{\Delta t |\partial \Omega_{ac}|} \int_{\Omega_a}^{p} \int_{\partial \Omega_{ac}} \mathbf{F}(\mathbf{Q}(\mathbf{x}, t)) \cdot \hat{\mathbf{n}}_{ac} \, dS \, dt. \quad (7)$$

In the context of classical edge-based two-point fluxes and first order explicit Euler time integration the numerical flux in normal direction would take the simple form

$$\hat{\mathbf{f}}_{ac} = \hat{\mathbf{f}}_{ac}(Q_c^n, Q_a^n, \hat{\mathbf{n}}_{ac}), \quad (8)$$

where as flux function $\hat{\mathbf{f}}_{ac}$ any classical 1D Riemann solver can be used (see e.g. [98]) and it is evaluated at the midpoint of the edge $\partial\Omega_{ac}$.

In this paper, instead, we employ point fluxes, i.e. genuinely multidimensional Riemann solvers, as in [14, 15, 35, 90, 16, 18] and thus the use of the midpoint rule is no more convenient. Thus, in order to approximate the integral in (7), one can employ the trapezoidal rule in such a way to include data at the vertexes $p$ of the control volume $\Omega_c$. Thus, the flux across the edge becomes

$$\hat{\mathbf{f}}_{ac} = \frac{1}{2} (\hat{\mathbf{F}}_{ac}^{-} + \hat{\mathbf{F}}_{ac}^{+}) \cdot \hat{\mathbf{n}}_{ac} = \frac{1}{2} (\hat{\mathbf{F}}_{q} + \hat{\mathbf{F}}_{p}) \cdot \hat{\mathbf{n}}_{ac}, \quad (9)$$

with $\hat{\mathbf{F}}_{ac}^{-} = \hat{\mathbf{F}}_{q}$ and $\hat{\mathbf{F}}_{ac}^{+} = \hat{\mathbf{F}}_{p}$ the numerical flux tensors evaluated in the two end points $q$ and $p$ of the edge $\partial\Omega_{ac}$. Using the corner normal defined by (4), we can rewrite the finite volume scheme (6) also using the
<figure><img src="image_3.png" /><figcaption>Figure 2: In this Figure we sketch the edge normal vectors $\mathbf{n}_{ac}$ and $\mathbf{n}_{bc}$ and the corner normal pointing from cell $\Omega_c$ to point $p$ which is denoted by $\mathbf{n}_{pc}$. Note that in general the intersection of edge $\overline{ac}$ with $\overline{pq}$ is not in the mid point of the edge $\overline{ac}$. Likewise for the intersection of the edge $\overline{bc}$ with $\overline{pr}$.</figcaption></figure>

corner normals and the corner fluxes instead of the edge normals and the edge fluxes as

$$ \mathbf{Q}_c^{n+1} = \mathbf{Q}_c^n - \frac{\Delta t}{|\Omega_c|} \sum_{\Omega_a \in \mathcal{N}_c} |\partial \Omega_{ac}| \frac{1}{2} (\hat{\mathbf{F}}_{ac}^{-} + \hat{\mathbf{F}}_{ac}^{+}) \cdot \mathbf{n}_{ac} = \mathbf{Q}_c^n - \frac{\Delta t}{|\Omega_c|} \sum_{p \in \mathcal{P}_c} \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc}. \quad (10) $$

The last equation allows to introduce a multidimensional framework generalizing the classical two states finite volume setting (8). We will directly define the corner flux $\hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc}$, which is now a function of $d+1$ states

$$ \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} = \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc}(\mathbf{Q}_a, \mathbf{Q}_b, \mathbf{Q}_c) \quad (11) $$

The above multidimensional numerical flux can be characterized with the following generalization of standard consistency and stability conditions.

**Consistency** A numerical corner flux is consistent if

$$ \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc}(\mathbf{Q}, \mathbf{Q}, \mathbf{Q}) = \mathbf{F}(\mathbf{Q}) \cdot \mathbf{n}_{pc}, \quad (12) $$

with $\mathbf{F}$ the physical flux.

**Conservation** Local conservation at each corner is equivalent to

$$ \sum_{c \in \mathcal{C}_p} \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} = 0. \quad (13) $$

**Entropy stability** Corner fluxes verify a generalization of Tadmor's entropy conservation condition [95] which reads

$$ \sum_{c \in \mathcal{C}_p} \mathbf{W}_c^t \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} = \sum_{c \in \mathcal{C}_p} \mathbf{\Psi}_c \cdot \mathbf{n}_{pc}, \quad (14) $$

having denoted by $\mathbf{W}$ the vector of entropy variables, and with $\mathbf{\Psi}$ the entropy potential vector. The derivation of this condition is provided in appendix Appendix A. Any corner flux with more dissipa-
tion than the one of an entropy conservative flux, defines and entropy stable scheme.

In the following we will provide two examples:the first can be seen as a multidimensional generalizations of the Osher-Solomon flux and the second of the 1d upwind flux difference splitting.

4. A new genuinely multidimensional Osher-type flux

Here, the nodal numerical flux tensor $\hat{\mathbf{F}}_p$ in (10) will be obtained via a novel extension of the Osher-Solomon flux [73] to the genuinely multidimensional case. Inspired by the ideas outlined in [43, 44], the genuinely multidimensional Osher-type flux, which provides a full numerical flux tensor instead of a usual numerical flux in a special normal direction, is assumed to have the following form

$$ \hat{\mathbf{F}}_p = \frac{1}{d+1} \sum_{c \in C_p} \mathbf{F}(\mathbf{Q}_c^n) - \frac{h}{d+1} \left( \int_{T_0} \mathrm{diag}(|A_i(\psi(\xi))| d\xi) \cdot \nabla_p^c \mathbf{Q}_c^n \right), \quad (15) $$

where the first term with the sum contains the centered part of the flux, while the term with the integral contains the numerical dissipation. In the above formula (15), $C_p$ is the set of cells $\Omega_c$ attached to a point $p$ (for example in Figure 2 $C_p = \{c, a, b\}$) corresponding to the vertexes of the simplex $T_p$. Instead, the unit simplex is denoted by $T_0$ and has the canonical vertexes 0 and 1 in 1D, (0, 0), (1, 0) and (0, 1) in 2D, and (0, 0, 0), (1, 0, 0), (0, 1, 0) and (0, 0, 1) in 3D. Moreover, $h$ is a characteristic length scale and $\nabla_p^c \mathbf{Q}_c^n$ is a discrete gradient operator whose precise definitions will be given later.

We now recall that the absolute value of a generic matrix $\mathbf{M}$ is defined as

$$ |\mathbf{M}| = \mathbf{R} |\mathbf{A}| \mathbf{R}^{-1}, \quad (16) $$

with $\mathbf{A}$ the diagonal matrix of the eigenvalues of the matrix $\mathbf{M}$, and $\mathbf{R}$ the associated matrix of right eigenvectors. Since the system (1) is assumed to be hyperbolic, all matrices $A_i(\psi(\xi))$ in (15) can be diagonalized, hence $|A_i(\psi(\xi))|$ can be computed.

While in the 1D Riemann solvers [43, 44] it is employed a linear 1D segment path $\psi = \psi(\xi_1)$, here we use a piecewise linear surface (manifold) inside the simplex $T_p$ given by

$$ \psi(\xi) = \sum_{k=1}^{d+1} \varphi_k(\xi) \mathbf{Q}_{c(k,p)}^n, \quad (17) $$

with $k \in \{1, 2, \dots, d+1\}$ a local vertex index in the simplex element $T_p$, $c = c(k, p)$ the global element number associated to node $p$ and vertex $k$ of the simplex $T_p$, and the $\varphi_k(\xi)$ being the classical P1 Lagrange basis functions given inside the unit simplex and expressed in terms of reference coordinates $\xi$. In the $d=2$ case, they read

$$ \varphi_1(\xi) = 1 - \xi_1 - \xi_2, \quad \varphi_2(\xi) = \xi_1, \quad \varphi_3(\xi) = \xi_2, \quad (18) $$

and for $d=3$ they read

$$ \varphi_1(\xi) = 1 - \xi_1 - \xi_2 - \xi_3, \quad \varphi_2(\xi) = \xi_1, \quad \varphi_3(\xi) = \xi_2, \quad \varphi_4(\xi) = \xi_3. \quad (19) $$

In the following we will denote the physical coordinates of the simplex $T_p$ by $\mathbf{X}_k = \mathbf{x}_{c(k,p)}$, $k \in \{1, 2, \dots, d+1\}$.

The discrete gradient $\nabla_p^c \mathbf{Q}_c^n$ (which is a constant inside the simplex $T_p$) in physical coordinates $\mathbf{x}$ can be immediately computed from the linear surface $\psi(\xi)$ as

$$ \nabla_p^c \mathbf{Q}_c^n = \frac{\partial \psi(\xi)}{\partial \xi} \mathbf{J}^{-1} = \sum_{k=1}^{d+1} \frac{\partial \varphi_k(\xi)}{\partial \xi} \mathbf{Q}_{c(k,p)}^n \mathbf{J}^{-1} = \frac{1}{|\Omega_p|} \sum_{c \in C_p} \mathbf{n}_{cp} \mathbf{Q}_c^n. \quad (20) $$
Here, **J** is the Jacobian matrix associated with the mapping from the unit simplex $T_0$ to the simplex $T_p$

$$ \mathbf{x} = \mathbf{x}(\xi) = \sum_{k=1}^{d+1} \varphi_k(\xi) \mathbf{X}_k, \qquad (21) $$

and reads

$$ \mathbf{J} = \frac{\partial \mathbf{x}}{\partial \xi} = \sum_{k=1}^{d+1} \frac{\partial \varphi_k(\xi)}{\partial \xi} \mathbf{X}_k = (\mathbf{X}_{d+1} - \mathbf{X}_1, \mathbf{X}_d - \mathbf{X}_1, \dots, \mathbf{X}_2 - \mathbf{X}_1). \qquad (22) $$

For the sake of simplicity we define the characteristic length scale as $h = (\frac{1}{d}|\mathbf{J}|)^{1/d}$. The integral in (15) is approximated at the aid of numerical quadrature formulae of suitable accuracy, see e.g. [94].

Inserting (20) into (15) yields the following final expression

$$ \hat{\mathbf{F}}_p = \frac{1}{d+1} \sum_{c \in C_p} \mathbf{F}(\mathbf{Q}_c^n) - \frac{h}{d+1} \left( \int_{T_0} \mathrm{diag}(|\mathbf{A}_i(\psi(\xi))|) d\xi \right) \cdot \left( \frac{1}{|\Omega_p|} \sum_{c \in C_p} \mathbf{n}_{cp} \mathbf{Q}_c^n \right), \qquad (23) $$

which together with the definition of the manifold (17) provides a closed expression for the numerical flux tensor $\hat{\mathbf{F}}_p$ in terms of the states $\mathbf{Q}_c^n$ adjacent to the vertex $\mathbf{x}_p$. This completes the description of the new genuinely multidimensional Osher-type Riemann solver.

The one-dimensional case as a special case of the general multidimensional setting. The 1D case already treated in [44, 43] is immediately included as a special case of the more general framework outlined in this Section. Indeed, it is enough to take as basis function $\varphi_1(\xi) = 1 - \xi_1$ and $\varphi_2(\xi) = \xi_1$, with the simplex $T_p$ for node $p = i + \frac{1}{2}$ being the interval $[x_i, x_i + 1]$, i.e. $C_p = \{i, i + 1\}$, and the 1D unit simplex being $T_0 = [0, 1]$. Thus the integration manifold reduces to the segment $\psi(\xi_1) = (1 - \xi_1)\mathbf{Q}_i^n + \xi_1\mathbf{Q}_{i+1}^n = \mathbf{Q}_i^n + \xi_1(\mathbf{Q}_{i+1}^n - \mathbf{Q}_i^n)$. Moreover, we simply have $h = \Delta x$ hence

$$ h \nabla_p^c \mathbf{Q}_c = \mathbf{Q}_{i+1}^n - \mathbf{Q}_i^n, \qquad (24) $$

which leads to the 1D numerical flux

$$ \hat{\mathbf{f}}_{i+\frac{1}{2}} = \frac{1}{2}(\mathbf{f}_1(\mathbf{Q}_i^n) + \mathbf{f}_1(\mathbf{Q}_{i+1}^n)) - \frac{1}{2}\left(\int_{T_0} |\mathbf{A}_1(\psi(\xi_1))| d\xi_1\right)(\mathbf{Q}_{i+1}^n - \mathbf{Q}_i^n), \qquad (25) $$

already proposed in [44, 43], where the integral was also approximated at the aid of Gaussian quadrature formulae of suitable accuracy.

## 5. Corner fluxes coming from the residual distribution framework

The construction of the multidimensional Osher-type solver exploits two geometrical properties of our Voronoi-like cells. In particular, it exploits the fact that every cell corner identifies uniquely $d + 1$ cells which can be used to define a unique local gradient using classical linear interpolation on the triangle $T_p$. Here we take this idea one step further. Let us start again from the multidimensional prototype with corner fluxes already given in (10)

$$ \mathbf{Q}_c^{n+1} = \mathbf{Q}_c^n - \frac{\Delta t}{|\Omega_c|} \sum_{p \in \mathcal{P}_c} \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc}, \qquad (26) $$

and focus on the definition of $\hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc}$ for a fixed $p$. We will now exploit the point-wise conservation constraint on the numerical flux (13). In particular, we require that for each $c \in C_p$ the numerical flux
should be given by the internal one plus a fluctuation:

$$ \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} = \mathbf{F}_c \cdot \mathbf{n}_{pc} + \boldsymbol{\phi}_{pc}, \qquad (27) $$

where $\mathbf{F}_c = \mathbf{F}(\mathbf{Q}_c^n)$ and the $\boldsymbol{\phi}_{pc}$ are nodal fluctuations which must be defined. The conservation condition (13) requires that these fluctuations verify the constraint

$$ 0 = \sum_{c \in C_p} \mathbf{F}_c \cdot \mathbf{n}_{pc} + \sum_{c \in C_p} \boldsymbol{\phi}_{pc}. \qquad (28) $$

Setting

$$ \boldsymbol{\phi}_p := \sum_{c \in C_p} \mathbf{F}_c \cdot \mathbf{n}_{cp}, \qquad (29) $$

Conservation at each node is equivalent to

$$ \sum_{c \in C_p} \boldsymbol{\phi}_{pc} = \boldsymbol{\phi}_p. \qquad (30) $$

Conservation at each corner is thus equivalent to the requirement that the fluctuations $\boldsymbol{\phi}_{pc}$ should be a splitting of the residual $\boldsymbol{\phi}_p$ defined in (29). The latter can be seen as the exact integral of the divergence of a linear approximation of the flux, over the triangle $\hat{T}_p$ of Figure 3, obtained by mirroring the triangle $rq$ w.r.t. the segment $[r, q]$. For this equivalence to be correct, the values of the unknown at the nodes of $\hat{T}_p$ should be the traces of the $d+1$ discrete finite volume solution at node $p$. In the first order case this trivially reduced to requiring that the nodal values of $\hat{T}_p$ should be $\mathbf{Q}_a, \mathbf{Q}_b$, and $\mathbf{Q}_c$. Note, that in general $\hat{T}_p$ does not coincide exactly with $T_p$, unless the segments $\overline{pq}$ and $\overline{ac}$ in Figure 2 intersect in correspondence of their mid-points, and similarly for $\overline{pr}$ and $\overline{bc}$.

The analogy presented above is a multidimensional generalization of the classical duality between finite volume methods and residual distributions (see e.g. [8] and references therein). Besides conservation, all the properties discussed previously can be generalized to this setting. In particular, the consistency requirement (12), reduces in this case to the condition

$$ \boldsymbol{\phi}_p(\mathbf{Q}, \mathbf{Q}, \mathbf{Q}) = 0. \qquad (31) $$

In Appendix A we also show that the generalized Tadmor shuffle condition (32), is equivalent to the entropy conservation condition

$$ \sum_{c \in C_p} \mathbf{W}_c^t \boldsymbol{\phi}_{pc} = \sum_{c \in C_p} \mathbf{G}_c \cdot \mathbf{n}_{cp}, \qquad (32) $$

where the right hand side in the last expression is the integral over $\hat{T}_p$ of the divergence of a linear approximation of the entropy flux $\mathbf{G}(\mathbf{Q})$.

The problem of finding a splitting (29) is classical in the context of second order Residual Distribution (RD) schemes [32, 9, 75]. We will recall hereafter some examples. Before that let us note that in
absence of a high order reconstruction we have

$$
\begin{align}
\mathbf{Q}_c^{n+1} &= \mathbf{Q}_c^n - \frac{\Delta t}{|\Omega_c|} \sum_{p \in \mathcal{P}_c} \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} \nonumber \\
&= \mathbf{Q}_c^n - \frac{\Delta t}{|\Omega_c|} \mathbf{F}_c \cdot \sum_{p \in \mathcal{P}_c} \mathbf{n}_{pc} - \frac{\Delta t}{|\Omega_c|} \sum_{p \in \mathcal{P}_c} \boldsymbol{\phi}_{pc} \tag{33} \\
&= \mathbf{Q}_c^n - \frac{\Delta t}{|\Omega_c|} \sum_{p \in \mathcal{P}_c} \boldsymbol{\phi}_{pc}, \nonumber
\end{align}
$$

having used the fact that $\sum_{p \in \mathcal{P}_c} \mathbf{n}_{pc} = 0$. This shows that *at first order* the above is not only an analogy, but we have a full equivalence between finite volumes with corner fluxes and RD schemes.

<figure><img src="image_4.png" /><figcaption>Figure 3: In this Figure we show the triangle $\hat{T}_p$ used to compute the residual distribution fluctuations $\boldsymbol{\phi}_p$ around each vertex $p$ of a polygonal element $c$.</figcaption></figure>

### 5.1. Example: multidimensional Osher-Solomon flux in fluctuation form

Using (27), we can readily deduce a fluctuation corresponding to the multidimensional Osher-type flux (15), that we have introduced in the previous section, which in 2D, and with the notation of Figures 2 and 3, reads

$$ \boldsymbol{\phi}_{pc}^{\text{OS}} = \frac{1}{3}(\mathbf{F}_a - \mathbf{F}_c) \cdot \mathbf{n}_{pc} + \frac{1}{3}(\mathbf{F}_b - \mathbf{F}_c) \cdot \mathbf{n}_{pc} - \frac{h}{3} \left[ \left( \int_{T_0} \operatorname{diag}(|A_i(\psi(\xi)|) d\xi \right) \cdot \nabla_h \mathbf{Q}_p \right] \cdot \mathbf{n}_{pc}. \quad (34) $$

Note that for piecewise constant data we can remove from the update (33)

### 5.2. Example: Rusanov method and an entropy stable splitting

A Rusanov flux, or multidimensional variant of local the Lax-Friedrich's scheme, is obtained with the following definition

$$ \boldsymbol{\phi}_{cp}^{\text{Rv}} = \frac{1}{3} \boldsymbol{\phi}_p + \alpha_p (\mathbf{Q}_c - \bar{\mathbf{Q}}_p), \quad (35) $$

where $\bar{\mathbf{Q}}_p$ is the average

$$ \bar{\mathbf{Q}}_p = \frac{1}{d+1} \sum_{c \in C_p} \mathbf{Q}_c. \quad (36) $$
Usually the constant $\alpha_p$ is chosen as the mesh size times the largest spectral radius of the system matrix projected in the normal directions. In other words

$$ \alpha_p = h \times \max(\rho(K_{n_p}(Q_c)), \rho(K_{n_{ap}}(Q_a)), \rho(K_{n_{bp}}(Q_b))). \quad (37) $$

The Rusanov splitting has many interesting properties, including a discrete maximum principle when using (convex combinations of) the explicit Euler time integration, and energy stability for linear symmetric systems (see e.g. [75, 32]). It is a very dissipative method, but its form is useful to construct schemes with desirable properties, or design corrections for existing schemes.

As an example, we can use the above splitting to obtain entropy conservative or stable schemes. Inspired by the previous definition, and by recent work by R. Abgrall [11], we can easily provide an extension of the previous definition which is strictly entropy conservative/stable. In particular, instead of using (37) to evaluated $\alpha_p$, we now compute it by explicitly imposing the entropy conservation condition (32):

$$ \alpha_p = \frac{\sum_{c \in C_p} \mathbf{G} \cdot \mathbf{n}_{cp} - \overline{\mathbf{W}}_p^t \phi_p}{(\sum_{c \in C_p} \mathbf{W}_c^t (\mathbf{Q}_c - \overline{\mathbf{Q}}_p))} \quad (38) $$

where $\overline{\mathbf{W}}_p = \sum_{c \in C_p} \mathbf{W}_c/(d + 1)$, and where the denominator can be easily shown to be always positive definite for all convex entropies.

An entropy stable flux can be obtained by adding to the above any positive contribution, or by adding to $\phi_{pc}$ an extra dissipation operator as e.g. the one of the multidimensional Osher-Solomon method.

### 5.3. Example: a multidimensional upwind flux using the N scheme

In the scalar case, the N scheme is Roe's multidimensional upwind optimum first order positive coefficient method on triangles [80, 81, 82]. Its extension to systems of conservation laws can be performed using an exact linearization of the flux Jacobians, as in Roe's one dimensional method, and either seeking a decomposition of the system in scalar waves, or via a compact matrix formulation [31, 102]. For cases in which an exact mean value linearization of the flux Jacobians is not available or impractical, the authors of [28, 78] have introduced a more general formulation compatible with the use of arbitrary linearized states for the system Jacobians. Using the notation introduced in (3), the N scheme splitting is defined by

$$ \phi_{cp}^{N} = \mathbf{K}_{n_{cp}}^{+} (\mathbf{Q}_{c} - \tilde{\mathbf{Q}}_{p}). \quad (39) $$

As in [28, 78] the tilde state $\tilde{\mathbf{Q}}_p$ is obtained by imposing the consistency condition (29):

$$ \tilde{\mathbf{Q}}_p = \mathbf{Q}_p^+ - \mathbf{N}^{-1} \phi_p, \quad (40) $$

where N denotes the matrix

$$ \mathbf{N} = \sum_{c \in C_p} \mathbf{K}_{n_{cp}}^{+}, \quad (41) $$

and where

$$ \mathbf{Q}_p^+ := \mathbf{N}^{-1} \sum_{c \in C_p} \mathbf{K}_{n_{cp}}^+ \mathbf{Q}_c. \quad (42) $$

The multidimensional upwind N scheme is the optimum first order scheme, where “optimum” means it has the least amount of numerical dissipation among first order positive coefficient discretizations for scalar advection [80, 81, 82]. The N scheme verifies a discrete maximum principle for scalar conservation laws, when using (convex combinations of) the explicit Euler time integration. For linear symmetric systems, it can also be shown to be energy stable (see e.g. [75, 4, 32] for all the previous statements). For nonlinear problems, only asymptotic entropy inequalities have been shown (see e.g. [4]), which do not preclude
violations of such condition. Indeed, when considering high Mach number flows, while the low dissipation of the N scheme has shown high potential to resolve complex flows with contact waves and boundary layers on coarse meshes [60], ad-hoc fixes for unwanted features as Carbuncle's and expansion shocks may be necessary [50, 93]. In this work, we also introduce a correction inspired by the Rusanov splitting:

$$ \Phi_{cp}^{N-C} = \Phi_{cp}^{N} + \mu_{C}(Q)(Q_{c} - \bar{Q}_{p}) \quad (43) $$

where $\mu_C(Q)$ is a small extra artificial viscosity only active in strong shocks. This modification has in particular no effect on contact waves. The precise definition of $\mu_C$, as well as the shock detection follows almost exactly [79].

Roe's one-dimensional upwind flux difference splitter as a special case of the multidimensional setting. In one space dimension, looking at a cell $i$, the computational triangle $\hat{T}_p$ reduces to a segment at the interface $i \pm 1/2$. The corner normals reduce to $\pm 1$. If we focus on the interface $i + 1/2$, we can easily see that the $K_{n,cp}$ matrices reduce to $\pm A_1(Q_{i+1/2})$, and in particular for cell $i$ the matrix $K_{n,cp}^+$ is replaced by $[-A_1(Q_{i+1/2})]^+ = -A_1^-(Q_{i+1/2})$. One easily shows that the nodal fluctuation becomes (see [8] for example)

$$ \Phi_{i+1/2} = \mathbf{f}_1(Q_{i+1}^n) - \mathbf{f}_1(Q_i^n) \quad (44) $$

while the matrix $N$ in (41) reduces to

$$ N_{i+1/2} = -A_1^-(Q_{i+1/2}) + A_1^+(Q_{i+1/2}) = |A_1(Q_{i+1/2})|. \quad (45) $$

All things assembled, and after few manipulations we obtain the numerical fluxes

$$ \begin{aligned} \hat{\mathbf{f}}_{i+\frac{1}{2}} &= \mathbf{f}_1(Q_i^n) - \frac{A_1^-(Q_{i+1/2})|A_1(Q_{i+1/2})|^{-1}}{2}(\mathbf{f}_1(Q_{i+1}^n) - \mathbf{f}_1(Q_i^n)) \\ &= \frac{1}{2}(\mathbf{f}_1(Q_i^n) + \mathbf{f}_1(Q_{i+1}^n)) - \frac{A_1(Q_{i+1/2})|A_1(Q_{i+1/2})|^{-1}}{2}(\mathbf{f}_1(Q_{i+1}^n) - \mathbf{f}_1(Q_i^n)). \end{aligned} \quad (46) $$

When $A_1(Q_{i+1/2})$ is evaluated averaging the Roe parameter, the above is precisely Roe's flux [83].

# 6. Higher order extension using the ADER approach with CWENO reconstruction

## 6.1. Spatial CWENO reconstruction

In the context of better than second order accurate finite volume schemes a reconstruction procedure is necessary to obtain better than piecewise constant approximations of the data within each cell. For this purpose one needs to compute a high order accurate but at the same time non-oscillatory polynomial representation $w_h^n(x, t^n)$ of the solution $Q(x, t^n)$ for each polygon $\Omega_c$, starting from the values of the cell averages $Q_a^n$ in $\Omega_c$ and its neighbors contained in a so-called reconstruction stencil $S_c^s$. For a polynomial of degree $M \ge 0$ in two space dimensions the number of degrees of freedom per conserved variable is $M = (M+1)(M+2)/2$. Moreover, as already mentioned in [20], on general unstructured meshes it is necessary to use more than $M$ cells to get a stable numerical scheme. Furthermore, one reconstruction stencil alone is not enough in order to obtain an essentially non-oscillatory scheme, but ENO/WENO reconstruction techniques [57, 52] relying on multiple stencils, as those introduced in [1, 55, 42, 104, 97, 101, 39] for unstructured meshes, are needed in order to circumvent the Godunov theorem which states that there are no better than first order accurate linear monotone schemes [51]. In this paper we follow in particular the successful approach of CWENO schemes introduced by Lewy, Puppo, Semplice and Russo et al. in a series of papers in [63, 64, 65, 27, 39], and which was for the first time applied on general unstructured polygonal meshes in [47, 46]. Here, we closely follow the work outlined in [39, 47], but for the sake of completeness and to make the paper self-contained for the reader, we report here the entire algorithm.
The reconstruction starts from the computation of a so-called central polynomial $\mathbf{P}_{\text{opt}}$ of degree $M$. In order to define $\mathbf{P}_{\text{opt}}$ in a robust manner, following [20, 42, 39], we consider a stencil $S_c^0$ which contains a total number of $n_e = f \cdot M$ control volume and which includes cell $\Omega_c$ and its neighbors

$$S_c^0 = \bigcup_{k=1}^{n_e} \Omega_{a(k)}, \qquad (47)$$

with the safety factor $f \ge 1.5$. Stencil $S_c^0$ is constructed by starting from the central polygon $\Omega_c$ and then adding recursively of neighbors and neighbors of elements that have been already included in the stencil, until the desired total number of $n_e$ elements is reached. The polynomial $\mathbf{P}_{\text{opt}}(\mathbf{x}, t^n)$ is then defined by imposing that its average on each cell $\Omega_a \in S_c^0$ matches the known cell average $\mathbf{Q}_a^n$. Since $n_e > M$, this of obviously leads to an overdetermined linear system, which is solved using a constrained least-squares technique introduced in [41]. Here, the polynomial $\mathbf{P}_{\text{opt}}$ has exactly the cell average $\mathbf{Q}_c^n$ on the polygon $\Omega_c$ and matches all other cell averages $\mathbf{Q}_a^n$ for $a \ne c$ only in a least-square sense. The polynomial $\mathbf{P}_{\text{opt}}$ is expressed in terms of spatial reconstruction basis functions $\psi_{\ell}(\mathbf{x}, t^n)$ of degree $M$,

$$\mathbf{P}_{\text{opt}}(\mathbf{x}, t^n) = \sum_{\ell=0}^{M-1} \psi_{\ell}(\mathbf{x}, t^n) \hat{\mathbf{p}}_{\ell,i}^n. \qquad (48)$$

Since there is no simple mapping to universal reference elements for general polygonal meshes, as in [47] we use the following reconstruction basis functions that employ rescaled and shifted Taylor monomials

$$\psi_{\ell}(\mathbf{x}, t^n)|_{P_i^n} = \frac{(x - x_c)^{p_{\ell}}}{p_{\ell}! h_c^{p_{\ell}}} \frac{(y - y_c)^{q_{\ell}}}{q_{\ell}! h_c^{q_{\ell}}}, \quad \ell = 0, \dots, M-1, \quad 0 \le p_{\ell} + q_{\ell} \le M, \qquad (49)$$

with $h_c$ the radius of the circumcircle and $\mathbf{x}_c$ the barycenter of the polygon $\Omega_c$. The expansion coefficients $\hat{\mathbf{p}}_{\ell,i}^n$ in (48) therefore represent the rescaled derivatives of the Taylor expansion about $\mathbf{x}_c$.

According to what was written above we now require integral conservation in the central stencil $S_c^0$ as follows:

$$\frac{1}{|\Omega_a|} \int_{\Omega_a} \mathbf{P}_{\text{opt}}(\mathbf{x}, t^n) d\mathbf{x} = \mathbf{Q}_a^n, \quad \forall \Omega_a \in S_c^0. \qquad (50)$$

and the integrals appearing in (50) are conveniently computed in each polygon by simply summing up the contributions over all sub-triangles of a polygon and on the sub-triangles we employ suitable conical products of the one-dimensional Gauss-Jacobi formula, see [94] for details.

In order to obtain a nonlinear reconstruction procedure, which is needed to avoid spurious oscillations in the vicinity of shock waves, the CWENO method uses also other polynomials of lower degree. In particular, in this work we employ piecewise linear reconstructions for the lower order reconstruction polynomials. Given a polygon $\Omega_c$ with $N_c$ Voronoi neighbors contained in the set $V(\Omega_c)$, we construct $N_c$ interpolating polynomials of degree $M^s = 1$ referred to as sectorial polynomials. In order to do this we consider $N_c$ stencils $S_c^s$ with $s \in [1, N_c]$, each of them containing three cells, namely the central polygon $\Omega_c$ and two consecutive neighbors in the set $V(\Omega_c)$. An example of such reconstruction stencils is reported for the case $M = 2$ in Figure 4.

For each stencil $S_c^s$ we compute a linear polynomial $\mathbf{P}_s(\mathbf{x}, t^n) \in \mathbb{P}_1$ by solving the reconstruction systems

$$\frac{1}{|\Omega_a|} \int_{\Omega_a} \mathbf{P}_s(\mathbf{x}, t^n) d\mathbf{x} = \mathbf{Q}_a^n, \quad \forall \Omega_a \in S_c^s \qquad (51)$$

on all sectorial stencils $S_c^s$. These reconstruction equations are not overdetermined and thus do not require the constrained least squares procedure. Following the general framework introduced in [27], we select a
set of positive coefficients $\lambda_0, \dots, \lambda_{N_p}$ such that

$$ \sum_{s=0}^{N_c} \lambda_s = 1 \qquad (52) $$

and we define a new auxiliary polynomial

$$ \mathbf{P}_0(\mathbf{x}, t^n) = \frac{1}{\lambda_0} \left( \left( \sum_{s=0}^{N_c} \lambda_s \right) \mathbf{P}_{\text{opt}}(\mathbf{x}, t^n) - \sum_{s=1}^{N_c} \lambda_s \mathbf{P}_s(\mathbf{x}, t^n) \right) \in \mathbb{P}_M, \qquad (53) $$

so that the linear combination of the polynomials $\mathbf{P}_0, \dots, \mathbf{P}_{N_c}$ with the coefficients $\lambda_0, \dots, \lambda_{N_c}$ is equal to $\mathbf{P}_{\text{opt}}$ and conservation is ensured. Specifically, we consider the linear weights used in [39], namely $\lambda_0 = 10^5$ for $S_i^0$ and $\lambda_s = 1$ for the sectorial stencils. These weights are then automatically renormalized in order to sum to unity, according to the requirement (52), see also eqn. (53) which contains the automatic rescaling. All polynomials $\mathbf{P}_s$ with $s \in [0, N_c]$ are then nonlinearly combined with each other as in classical WENO schemes, see e.g. [57, 55, 41]. We therefore obtain the final nonlinear CWENO reconstruction polynomial $\mathbf{w}_h(\mathbf{x}, t^n)$ in $\Omega_c$ as

$$ \mathbf{w}_h(\mathbf{x}, t^n) = \sum_{s=0}^{N_c} \omega_s \mathbf{P}_s(\mathbf{x}, t^n), \quad \mathbf{x} \in \Omega_c, \qquad (54) $$

<figure><img src="image_5.png" /><figcaption>Figure 4: Illustration of the different reconstruction stencils used for the CWENO reconstruction of order three ($M=2$) with a safety factor of $f=1.5$ for a pentagonal central element $\Omega_c$ in blue. Top-left: central stencil needed for the reconstruction of $\mathbf{P}_{\text{opt}}$ (in light blue). In the other panels we report the 5 sectorial stencils containing the element itself and two consecutive neighbors needed to reconstruct piecewise linear sectorial polynomials.</figcaption></figure>
where the normalized nonlinear weights $\omega_s$ read

$$ \omega_s = \frac{\tilde{\omega}_s}{\sum_{m=0}^{N_c} \tilde{\omega}_m}, \quad \text{with} \quad \tilde{\omega}_s = \frac{\lambda_s}{(\sigma_s + \epsilon)^r}. \qquad (55) $$

In the equation above the non-normalized weights $\tilde{\omega}_s$ depend on the linear weights $\lambda_s$ and the oscillation indicators $\sigma_s$. The parameters $\epsilon = 10^{-14}$ and $r = 4$ are chosen according to [41, 39, 47]. As in [47], the oscillation indicators $\sigma_s$ appearing in (55) are simply given by the quadratic form

$$ \sigma_s = \sum_l (\hat{\mathbf{p}}_{l,i}^{n,s})^2. \qquad (56) $$

Note that in smooth regions of the domain the central stencil with the optimal order polynomial is preferred due to the large linear weight $\lambda_0$, while in the vicinity of discontinuities the reconstruction switches to the lower order sectorial polynomials to reduce oscillations.

## 6.2. High order element-local space-time ADER predictor

Unlike the original ADER schemes of Toro and Titarev [96, 100] here we employ the element-local space-time DG predictor as introduced in [40, 38]. This local space-time predictor solution is valid locally inside each $\Omega_c$ for the current time interval $[t^n, t^{n+1}]$ and is given by high order piecewise polynomials $\mathbf{q}_h(\mathbf{x}, t)$ of degree $M$ in space-time which read

$$ \mathbf{q}_h = \mathbf{q}_h(\mathbf{x}, t) = \sum_{\ell=0}^{Q-1} \theta_{\ell}(\mathbf{x}, t) \hat{\mathbf{q}}_{\ell}^{n}, \quad (\mathbf{x}, t) \in \Omega_c \times [t^n, t^{n+1}], \quad Q = M(M+1), \qquad (57) $$

with $\theta_{\ell}(\mathbf{x}, t)$ the modal space-time basis already introduced in [47, 48]:

$$ \theta_{\ell}(\mathbf{x}, y, t)|_{C_i^n} = \frac{(x - x_c)^{p_{\ell}}}{p_{\ell}! h_c^{p_{\ell}}} \frac{(y - y_c)^{q_{\ell}}}{q_{\ell}! h_c^{q_{\ell}}} \frac{(t - t^n)^{q_{\ell}}}{q_{\ell}! h_c^{q_{\ell}}}, \quad \ell = 0, \dots, Q-1, \quad 0 \le p_{\ell} + q_{\ell} + r_{\ell} \le M. \qquad (58) $$

To obtain a weak formulation of the governing PDE system in space-time we multiply the system of conservation laws (1) by a test function $\theta_k$, integrate over the space-time control volume $C_c = \Omega_c \times [t^n, t^{n+1}]$ and substitute the discrete solution $\mathbf{q}_h(\mathbf{x}, t)$ for $\mathbf{Q}$. This yields

$$ \int_{C_c} \theta_k(\mathbf{x}, t) \frac{\partial \mathbf{q}_h}{\partial t} d\mathbf{x}dt + \int_{C_c} \theta_k(\mathbf{x}, t) \nabla \cdot \mathbf{F}(\mathbf{q}_h) d\mathbf{x}dt = 0. \qquad (59) $$

In order to introduce the initial condition in this weak formulation of the local Cauchy problem, given by the spatial reconstruction polynomial obtained via the CWENO procedure outlined in the previous section, we integrate the first term by parts in time and then use upwinding in time (causality principle), thus obtaining

$$ \int_{\Omega_c} \theta_k(\mathbf{x}, t^{n+1}) \mathbf{q}_h(\mathbf{x}, t^{n+1}) d\mathbf{x} - \int_{C_c} \frac{\partial}{\partial t} \theta_k(\mathbf{x}, t) \mathbf{q}_h d\mathbf{x}dt + \int_{C_c} \theta_k(\mathbf{x}, t) \nabla \cdot \mathbf{F}(\mathbf{q}_h) d\mathbf{x}dt = \int_{\Omega_c} \theta_k(\mathbf{x}, t^n) \mathbf{w}_h(\mathbf{x}, t^n) d\mathbf{x}. \qquad (60) $$

The element-local space-time predictor $\mathbf{q}_h^n$ can then be easily computed via a simple fixed-point iteration, the convergence of which has been shown in [23]. For details, the reader is referred to [38, 56, 23].

## 6.3. Final high order ADER-CWENO scheme based on multidimensional Riemann solvers

For better than second order accurate schemes in space and time we need higher order accurate numerical quadrature formulas to integrate the fluxes along the element edges. As done in [35, 17, 22] we use the
Keplerian barrel rule<fnref n="1" />, which is a Gauss-Lobatto quadrature formula with 3 quadrature points, two being located at the vertexes of an edge and one in the edge barycenter. As such, the high order version of the ADER-CWENO scheme based on vertex fluxes reads

$$Q_c^{n+1} = Q_c^n - \frac{1}{3} \frac{\Delta t}{|\Omega_c|} \sum_{p \in \mathcal{P}_c} \mathbf{n}_{pc} \cdot \hat{\mathbf{F}}_p - \frac{2}{3} \frac{\Delta t}{|\Omega_c|} \sum_{\Omega_a \in \mathcal{N}_c} |\partial \Omega_{ac}| \hat{\mathbf{f}}_{ac}, \quad (61)$$

with $\hat{\mathbf{F}}_p$ the time-averaged vertex fluxes computed in the vertex-extrapolated states of the high order space-time predictor solution $\mathbf{q}_h$ and $\hat{\mathbf{f}}_{ac}$ the time-averaged edge fluxes computed at the aid of a conventional 1D Riemann solver. Throughout this paper we employ the 1D Riemann solvers to which the multidimensional fluxes reduce to: the 1D Osher solver [44, 43] for the flux of Section 4, and the Roe method for the multidimensional upwind N scheme solver. Both choices are consistent with the aim of this paper of using complete Riemann solvers anywhere. More precisely, to compute the time average flux $\hat{\mathbf{f}}_{ac}$ with the appropriate order of accuracy, we evaluate it as

$$\hat{\mathbf{f}}_{ac} = \sum_{j} \omega_j \hat{\mathbf{f}}_{ac} (\mathbf{q}_h(\mathbf{x}_{ac}^{-}, t^n + \tau_j \Delta t), \mathbf{q}_h(\mathbf{x}_{ac}^{+}, t^n + \tau_j \Delta t), \hat{\mathbf{n}}_{ac}), \quad (62)$$

where $\tau_j$ and $\omega_j$ are the points and weights of a classical Gauss-Legendre quadrature formula of suitable order in 1D. Two quadrature points would be enough both for third and fourth order of accuracy, however we employ $M + 1$ points to respect the underlying original structure of our ADER solver. Furthermore, $\mathbf{q}_{\bar{h}} = \mathbf{q}_h(\mathbf{x}_{ac}^{-}, t^n + \tau_j \Delta t)$ and $\mathbf{q}_h(\mathbf{x}_{ac}^{+}, t^n + \tau_j \Delta t)$ are the boundary-extrapolated states that enter the 1D Riemann solver and $\mathbf{x}_{ac}^{-}$ and $\mathbf{x}_{ac}^{+}$ is the edge barycenter seen from the left element $\Omega_c$ and from the right element $\Omega_a$, respectively. Likewise, the vertex flux must be computed in the vertex-extrapolated states. In case of the new multidimensional Osher flux we obtain

$$\hat{\mathbf{F}}_p = \sum_j \omega_j \frac{1}{d+1} \sum_{c \in C_p} \mathbf{F}(\mathbf{q}_h(\mathbf{x}_{cp}, t^n + \tau_j \Delta t)) - \sum_j \omega_j \frac{h}{d+1} \left( \int_{T_0} \mathrm{diag}(|A_i(\psi(\xi))|) d\xi \right) \cdot \left( \frac{1}{|\Omega_p|} \sum_{c \in C_p} \mathbf{n}_{cp} \mathbf{q}_h(\mathbf{x}_{cp}, t^n + \tau_j \Delta t) \right). \quad (63)$$

Here, $\mathbf{x}_{cp}$ is the vertex coordinate $\mathbf{x}_p$ taken from inside control volume $\Omega_c$, i.e. the predictor $\mathbf{q}_h(\mathbf{x}_{cp}, t^n + \tau_j \Delta t)$ is also taken from within $\Omega_c$ and is computed in the location $\mathbf{x}_p$ of the vertex.

The other node fluxes presented in this paper can be computed in a similar manner, simply replacing $Q_c^n$ with the vertex-extrapolated value $\mathbf{q}_h(\mathbf{x}_{cp}, t^n + \tau_j \Delta t)$ and using Gaussian quadrature in time to compute the time-averaged flux.

## 7. Numerical results

In order to numerically verify the capabilities of our multidimensional Riemann solvers and the applicability of our strategy for their high order extension in space and time, we test them on a series of benchmarks typical of the compressible Euler equations of gas-dynamics.

<footnote n="1">¹This quadrature formula is widely known under the name Simpson rule, but it had already been found more than 100 years earlier by Johannes Kepler in 1615 in his work *Nova stereometria doliorium vinariorum, in primis Austriaci*, see [59], hence we deliberately prefer to refer to the original reference.</footnote>
<table><caption>Table 1: Isentropic Shu-type vortex. Numerical convergence results for our ADER-CWENO scheme based on the multidimensional Osher and the N scheme as Riemann solvers and employed from the first order ($M = 0$) to the fourth order ($M = 3$) of accuracy. The error norms refer to the variable $\rho$ at time $t = 1.0$ in the $L_2$ norm.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">h</th><th colspan="2">MultiD Osher</th><th colspan="2">N scheme</th></tr><tr><th>ε<sub>L2</sub>(ρ)</th><th>O(ρ)</th><th>ε<sub>L2</sub>(ρ)</th><th>O(ρ)</th></tr></thead><tbody><tr><td rowspan="4">M = 0 → O1</td><td>3.07E-2</td><td>5.07E-2</td><td>-</td><td>2.18E-2</td><td>-</td></tr><tr><td>2.47E-2</td><td>4.11E-2</td><td>0.95</td><td>1.74E-2</td><td>1.02</td></tr><tr><td>2.06E-2</td><td>3.46E-2</td><td>0.96</td><td>1.45E-2</td><td>1.01</td></tr><tr><td>1.77E-2</td><td>2.99E-2</td><td>0.96</td><td>1.24E-2</td><td>1.01</td></tr><tr><td rowspan="4">M = 1 → O2</td><td>4.89E-2</td><td>1.65E-3</td><td>-</td><td>1.38E-3</td><td>-</td></tr><tr><td>4.10E-2</td><td>1.08E-3</td><td>2.40</td><td>9.58E-4</td><td>2.03</td></tr><tr><td>3.07E-2</td><td>5.63E-4</td><td>2.26</td><td>5.37E-4</td><td>2.03</td></tr><tr><td>2.47E-2</td><td>3.55E-4</td><td>2.10</td><td>3.48E-4</td><td>1.98</td></tr><tr><td rowspan="4">M = 2 → O3</td><td>6.11E-2</td><td>2.62E-3</td><td>-</td><td>2.25E-3</td><td>-</td></tr><tr><td>4.89E-2</td><td>1.29E-3</td><td>3.16</td><td>1.09E-3</td><td>3.26</td></tr><tr><td>4.10E-2</td><td>7.60E-4</td><td>2.99</td><td>6.41E-4</td><td>2.96</td></tr><tr><td>3.07E-2</td><td>3.17E-4</td><td>3.05</td><td>2.67E-4</td><td>3.07</td></tr><tr><td rowspan="4">M = 3 → O4</td><td>1.19E-1</td><td>1.54E-2</td><td>-</td><td>1.56E-2</td><td>-</td></tr><tr><td>6.11E-2</td><td>7.30E-4</td><td>4.56</td><td>6.83E-4</td><td>4.68</td></tr><tr><td>4.89E-2</td><td>3.07E-4</td><td>3.87</td><td>2.87E-4</td><td>3.88</td></tr><tr><td>4.10E-2</td><td>1.34E-4</td><td>3.64</td><td>1.44E-4</td><td>3.87</td></tr></tbody></table>

The Euler equations can be cast in the form (1) in the two-dimensional framework by choosing

$$ Q = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho E \end{pmatrix}, \quad F = \begin{pmatrix} \rho u & \rho v \\ \rho u^2 + p & \rho uv \\ \rho uv & \rho v^2 + p \\ u(\rho E + p) & v(\rho E + p) \end{pmatrix}. \qquad (64) $$

Here, the vector of conserved variables $Q$ involves the fluid density $\rho$, the momentum density vector $\rho v = (\rho u, \rho v)$ and the total energy density $\rho E$. The fluid pressure $p$ is related to the conserved quantities $Q$ using the equation of state for an ideal gas

$$ p = (\gamma - 1) \left( \rho E - \frac{1}{2} \rho v^2 \right), \qquad (65) $$

where $\gamma$ ($\gamma = 1.4$ in our tests) is the ratio of specific heats so that the speed of sound takes the form $c = \sqrt{\frac{\gamma p}{\rho}}$.

Let us notice that for all the benchmarks we have deliberately chosen to work with coarse meshes in order to better appreciate the influence of the Riemann solvers on the final result and to make the results quickly replicable.

## 7.1. Order of convergence of our high order ADER-CWENO schemes based on multidimensional Riemann solvers

To check the order of accuracy of our multidimensional Riemann solvers embedded in the strategy presented in Section 6, we employ a smooth isentropic vortex test case inspired to the one proposed in the example 3.3. of [54], which is a stationary equilibrium of the Euler system of equations.

The computational domain is the square $\Omega = [-10; 20] \times [-10; 20]$ with wall boundary conditions. The initial condition is given by some perturbations $\delta$ that are superimposed onto a homogeneous background field $V_0 = (\rho, u, v, p) = (1, 0, 0, 1)$ (note that w.r.t. [54] here we do not translate the vortex). The
perturbations for density and pressure are

$$ \delta\rho = (1 + \delta T)^{\frac{1}{\gamma-1}} - 1, \quad \delta p = (1 + \delta T)^{\frac{\gamma}{\gamma-1}} - 1, \qquad (66) $$

with the temperature fluctuation $\delta T = -\frac{(\gamma-1)\epsilon^2}{8\gamma\pi^2}e^{1-r^2}$ and the vortex strength is $\epsilon = 5$. The velocity field is affected by the following perturbations

$$ \left( \begin{array}{c} \delta u \\ \delta v \end{array} \right) = \frac{\epsilon}{2\pi} e^{\frac{1-r^2}{2}} \left( \begin{array}{c} -(y-5) \\ (x-5) \end{array} \right). \qquad (67) $$

When the simulation reaches the final time $t = 1$, we compute the discrete $L^2$ error norms by evaluating the difference between the reconstructed numerical solution $w_h$ and the exact solution $Q_{\text{ref}}$ as

$$ \epsilon_{L_2} = \|w_h - Q_{\text{ref}}\|_{L^2} = \sum_{c=1}^{N_{\omega}} \left( \int_{\Omega_c} |w_h(\mathbf{x}, t_f) - Q_{\text{ref}}(\mathbf{x}, t_f)| dxdy \right)^{\frac{1}{2}}, $$

using Gaussian quadrature rules of appropriate order. The obtained orders of convergence, which are achieved as expected, are reported in Table 1 for both the multidimensional Osher and the N scheme Riemann solvers employed in our CWENO-ADER scheme of order from one to four.

## 7.2. Contact and shear waves

In this section we concentrate on the study of isolated contact waves and shear waves.

**Steady contact wave.** We first notice that both our novel multidimensional Riemann solvers as well as the classical complete Riemann solvers that we have used for comparison reasons (Osher and Roe solvers) preserve steady contact waves up to machine even when the position of the initial discontinuity is not aligned with the mesh. As numerical evidence of this property, we show in Figure 5 the discontinuous initial density profile of our test case

$$ (\rho, u, v, p)(\mathbf{x}) = \begin{cases} (1.0, 0, 0, 1.0) & \text{if } y > 0 \ \& \ y < -5x + 0.5, \\ (0.1, 0, 0, 1.0) & \text{if } y > 0 \ \& \ y \geq -5x + 0.5, \\ (1.0, 0, 0, 1.0) & \text{if } y < 0 \ \& \ y \geq 5x - 0.5, \\ (0.1, 0, 0, 1.0) & \text{if } y < 0 \ \& \ y < 5x - 0.5, \end{cases} \qquad (68) $$

<figure><img src="image_6.png" /><figcaption>Figure 5: Steady contact wave. On the left we show the discontinuous initial density profile, which is not aligned with the mesh. Then we show the 1d cut along $y = 0$ of the final density (middle) and the final total velocity (right) obtained with our novel multidimensional solvers and some classical 1d solvers. Note that the not aligned steady contact wave is maintained with machine precision.</figcaption></figure>
simulated on the square [-0.5; 0.5] x [-0.5; 0.5] with periodic boundary conditions on top and bottom boundaries, and the 1d cut along y = 0 of the density and the total velocity (V = $\sqrt{u^2 + v^2}$) obtained at time t=1.0 with our first order numerical schemes, which coincide with the initial condition up to machine precision.

*Isolated shear waves.* We then consider an isolated shear wave characterized by the following initial conditions

$$ (\rho, u, v, p)(\mathbf{x}) = \begin{cases} (1.0, 0, 0.1, 1.0) & \text{if } x < 0, \\ (1.0, 0, -0.1, 1.0) & \text{if } x \ge 0, \end{cases} \qquad (69) $$

simulated on a square of dimension [-0.5; 0.5] x [-0.5; 0.5] with periodic boundary on top and bottom boundaries. We remark that, despite the mathematical discontinuity is aligned with $x = 0$, our unstructured polygonal mesh, if not constructed by imposing some constraints on the position of the initial generator points c, thus it does not maintain this alignment property. Hence, the results in Figure 6 refer to a shear wave not aligned with the mesh and demonstrate that our novel multidimensional solvers significantly reduce the numerical dissipation.

In addition, we would like to underline that, if the generator points c are chosen to be aligned with the discontinuity, then our multidimensional Osher Riemann solvers, with the dissipation part of 23 computed using primitive variables (circles), is able to simulate the isolated shear waves with machine precision, as can be noticed in Figure 7.

### 7.3. Circular Sod explosion problem

This circular explosion problem is a multidimensional extension of the classical Sod test case simulated on a square of dimension [-1; 1] x [-1; 1] with the following initial conditions composed of two different states, separated by a discontinuity at radius $r_d = 0.5$

$$ (\rho, u, v, p)(\mathbf{x}) = \begin{cases} (1.0, 0, 0, 1) & \text{if } r \le r_d, \\ (0.125, 0, 0, 0.1) & \text{if } r > r_d. \end{cases} \qquad (70) $$

<figure><img src="image_7.png" /><figcaption>Figure 6: Shear wave not aligned with the tassellation. Here, we show the numerical results obtained with our first order schemes and our novel multidimensional Rieman solvers (squares) compared with those obtained with standard complete solvers (dashed lines). Our novel solvers significantly reduce the numerical dissipation.</figcaption></figure>
Here, we numerically solve this problem with our novel multidimensional solvers and our first and third order finite volume schemes on a mesh with resolution $h \approx 1/50$. The final time is chosen to be $t_f = 0.25$, so that the shock wave does not cross the external boundary of the domain, where wall boundary conditions are imposed. We have obtained a reference solution thanks to the rotational symmetry of the problem which reduces to a one-dimensional problem with geometric source terms, which we have solved by using

<figure><img src="image_8.png" /><figcaption>Figure 7: Shear wave simulated on a tassellation where the generator points $c_i$ are aligned with the discontinuity. One can notice that our multidimensional Osher solver in primitive variables (red circles) exhibits zero numerical dissipation for the simulation of shear waves which are aligned with the mesh.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Circular Sod explosion. We show the numerical results obtained with the first order (M=0) and the third order (M=2) CWENO-ADER schemes on the same mesh and our multidimensional Oscher scheme. (The results obtained with the N scheme are really similar so we omit them).</figcaption></figure>
a classical second order TVD scheme on a very fine mesh.

The obtained results are reported in Figure 8 and show the increased resolution power of our high order approach.

### 7.3.1. Lax shock tube

We continue with our set of Riemann problems by solving the Lax shock tube, originally presented in [62], which in addition to one shock wave also shows a contact wave and a rarefaction fun. The initial conditions for this problem are

$$ (\rho, u, v, p)(\mathbf{x}) = \begin{cases} (0.5, 0, 0, 0.571) & \text{if } x > 0.5, \\ (0.445, 0.698, 0, 3.528) & \text{if } x \le 0.5, \end{cases} \qquad (71) $$

and it is solved on the domain $\Omega = [0, 1] \times [0, 0.1]$ discretized with a coarse mesh with resolution $h \approx 1/50$. The final simulation time is $t_f = 0.14$. We report the numerical results obtained with our multidimensional Riemann solvers and our first and third order finite volume schemes in Figure 9.

## 7.4. Two dimensional Riemann Problem

In this Section, we consider some of the truly multidimensional Riemann problems presented and studied in [61].

**Configuration 12.** We start with the configuration number 12 which involves both steady contact discontinuities, and strong moving shocks whose interaction produces two triple points from which parallel contact lines emerge. From the initial singularity a strong fluid jet emerges whose resolution is enhanced by high order and less dissipative methods as our multidimensional solvers.

<figure><img src="image_10.png" /><figcaption>Figure 9: Lax shock tube. We show the numerical results obtained with the first order (M=0) and the third order (M=2) CWENO-ADER schemes on the same mesh and our N scheme Riemann solver. (The results obtained with the multidimensional Osher solver are really similar so we omit them).</figcaption></figure>
Here, the computational domain is the square $\Omega = [-2, 2] \times [-2, 2]$ discretized with a mesh with resolution $h \approx 1/110$, and the discontinuous initial conditions read as follows

$$ (\rho, u, v, p)(\mathbf{x}) = \begin{cases} (1.0, 0.7276, 0.0, 1.0) & \text{if } x < 0 \ \& \ y > 0, \\ (0.8, 0.0, 0.0, 1.0) & \text{if } x < 0 \ \& \ y < 0, \\ (0.5313, 0.0, 0.0, 0.4) & \text{if } x > 0 \ \& \ y > 0, \\ (1.0, 0.0, 0.7276, 1.0) & \text{if } x > 0 \ \& \ y < 0. \end{cases} \qquad (72) $$

<figure><img src="image_11.png" /><figcaption>Figure 10: Two dimensional Riemann Problem (configuration 12 of [61]). Here we report the numerical results obtained with our first order P0, second order P1, and third order P2 schemes. In particular, we show the density profile with 40 contours line in the interval [0.54, 1.7]. We refer to the test case description for detailed comments on the obtained results.</figcaption></figure>
In Figure 10, we report the numerical results obtained at time $t_f = 0.25$ with a standard 1D Osher solver compared with those achieved with our multidimensional solvers of first, second, and third order of accuracy. Moreover, in Figure 11 we show the cut of the density profile along $x = -1.2$ and $y = -1.2$. We can clearly notice that both our multidimensional solvers exhibit an increased resolution w.r.t. the 1D solver showing many more vortexes along the traveling shock discontinuities. We would also like to remark that the third order multidimensional schemes (in the last two panel of the Figure), as presented in Section 6, which also involve an edge contribution computed with a 1D solver, have less resolution power w.r.t. the second order schemes, which only use corner fluxes. Moreover, a variant of the third order scheme (left-bottom panel), employing a third order space-time polynomial reconstruction but a second order order flux computation, based only on corner fluxes, effectively provides numerical results more accurate w.r.t the second order schemes. However, due to the lack of sufficient accuracy of the numerical quadrature, this variant is only second order accurate. This feature will be object of future investigations.

**Configuration 3.** We now consider the configuration number 3 of the same reference. The initial solution involves again four states corresponding to two oblique shocks and two normal shocks converging in one point. The interaction leads to the appearance of a strong normal shock connecting the two oblique shocks, and of two lambda shocks stemming from the normal ones. The lambda shocks lead to the appearance of two strong moving contact discontinuities, and of a strong jet of fluid whose resolution is enhanced by high order and less dissipative methods as our multidimensional solvers.

Here, the computational domain is the square $\Omega = [0, 1.2] \times [0, 1.2]$ discretized with meshes of resolution $h \approx 1/200$ or $h \approx 1/400$, and the discontinuous initial conditions read as follows

$$ (\rho, u, v, p)(\mathbf{x}) = \begin{cases} (0.5323, 1.206, 0.0, 0.3) & \text{if } x < 1 \ \& \ y > 1, \\ (0.138, 1.206, 1.206, 0.029) & \text{if } x < 1 \ \& \ y < 1, \\ (1.5, 0.0, 0.0, 1.5) & \text{if } x > 1 \ \& \ y > 1, \\ (0.5323, 0.0, 1.206, 0.3) & \text{if } x > 1 \ \& \ y < 1. \end{cases} \qquad (73) $$

<figure><img src="image_12.png" /><figcaption>Figure 11: Cut along $x = -1.2$ (left) and $y = -1.2$ (right) for the density values of the solution of the two dimensional Riemann Problem (configuration 12 of [61]). One can easily notice that the discontinuity between the two steady states is captured in a sharper manner by our multidimensional solvers.</figcaption></figure>
The final computational time is $t_f = 1.0$.

We report the results obtained with our first and second order finite volume schemes in Figure 12 so to compare the different resolutions achieved with the classical Roe edge solver and with our multidimensional corner fluxes.

<figure><img src="image_13.png" /><figcaption>Figure 12: Two dimensional Riemann problem (configuration 3 of [61]). Here we report the numerical results obtained with our first order P0 and second order P1 schemes over a coarse mesh with $h \approx 1/100$ and a finer one with $h \approx 1/400$. In particular, we show the density profile with 40 contours line in the interval $[0.2, 1.6]$. We can notice that the N scheme, even if the carbuncle fix is active, captures clearly the shear instabilities, since the carbuncle fix is active only across shock waves. Finally, we recall that the asymmetry of the results is due to the use of fully unstructured meshes.</figcaption></figure>
7.5. Hypersonic flow past a blunt body

We consider now a hypersonic (Mach = 5 and Mach = 10) flow past the forebody of a circular cylin-
der [58, 49, 25]. As computational domain we take the half ring between two circles of internal radius
r₀ = 1 and external r₁ = 3 (centered at the origin) and with x < 0. We discretize it with a coarse mesh
of characteristic size h ≈ 1/70. The initial condition is given by an incoming hypersonic flow such that
(ρ, u, v, p)(x) = (√γ, Mach, 0, 1) and the boundary conditions are of wall type on r₀, of transmissive type

<figure><img src="image_14.png" /><figcaption>Figure 13: Here, we show the density profile of the blunt body test case at Mach=5 (top line) and Mach=10 (bottom line) with 21 contours line in [1.2, 8.3]. In particular, we compare the behavior of our multidimensional solvers with the classical 1d edge solvers inserted in a P0 first order finite volume scheme. We first notice that the Roe scheme, which is usually affected by the carbuncle phenomenon on triangular meshes, here on polygonal meshes does not suffer of this problem. The incomplete simple Rusanov scheme is, as usual, carbuncle-free but quite dissipative. Instead, we do not report any results for the 1D Osher scheme because it crashes at the very beginning of the simulation for the early apparition of the carbuncle phenomenon. On the contrary, our multidimensional extension of the Osher scheme, without the needing of any trick, provides high resolved results both at Mach=5 and Mach=10 with just a small carbuncle effect at high Mach which however do not affect the simulation. Finally, the N scheme, which we recall that it is implemented with a well-known carbuncle fix, is completely carbuncle free and shows a reduced dissipation despite the employed fix.</figcaption></figure>
on the right vertical boundaries and of Dirichlet type on $r_1$.

We report the obtained results in Figure 13 together with detailed comments regarding in particular the behavior of our schemes w.r.t the well-known appearance of the carbuncle phenomenon.

### 7.6. Circular two-dimensional explosion with boundary reflections

We conclude our sets of benchmarks with an explosion problem which has an initial configuration similar to the Sod test case but it is performed on a larger domain and for a longer time, so to see many more reflections and interactions around the origin. The final computational time is $t_f = 3.2$. The initial

<figure><img src="image_15.png" /><figcaption>Figure 14: Circular two-dimensional explosion with boundary reflections, comparison of Osher type fluxes at first and second order of accuracy. In particular, we show the density contours at time $t_f = 3.2$ with 16 contours level in [0.1, 0, 2].</figcaption></figure>
conditions for this problem are

$$ (\rho, u, v, p)(\mathbf{x}) = \begin{cases} (1.0, 0.0, 0.0, 1.0) & \text{if } x > 0.4, \\ (0.125, 0.0, 0.0, 0.1) & \text{if } x \le 0.4, \end{cases} \qquad (74) $$

and the domain of interest is $\Omega = [0, 1.5] \times [0, 1.5]$. (We have solved the problem on a larger domain $\Omega = [-3.5, 3.5] \times [-3.5, 3.5]$ to avoid the effect of external boundary conditions).

We have simulated this test case with our first and second order finite volume schemes, and we report in Figure 14 the results obtained with the Osher solvers and in Figure 15 those obtained with the Roe scheme and the N scheme.

<figure><img src="image_16.png" /><figcaption>Figure 15: Circular two-dimensional explosion with boundary reflections, comparison of Roe and N scheme solvers at first and second order of accuracy. In particular, we show the density contours at time $t_f = 3.2$ with 16 contours level in [0.1, 0, 2].</figcaption></figure>
## 8. Conclusions and outlook to future work

In this paper we have designed new high order accurate fully-discrete one-step finite volume schemes that are cell centered and employ nodal fluxes. High order of accuracy in space is achieved via the CWENO reconstruction procedure and high order of accuracy in time is obtained at the aid of the ADER approach, making use of a weak formulation of the PDE system in space-time. With respect to existing high order cell-centered finite volume schemes that employ vertex-based numerical fluxes this paper introduces two novelties:

1. the formulation of a new genuinely multidimensional Osher-type Riemann solver for general non-linear systems of hyperbolic conservation laws, see eqn. (23) for the algebraic expression of the numerical flux and (17) for the generalization of the simple 1D integration path used in [44, 43] to a manifold in multiple space dimensions. The very peculiar feature of this new vertex-based Riemann solver is the fact that it produces an entire numerical flux tensor rather than just a numerical flux projected into a particular direction. The mathematical structure of the flux is a simple and direct extension of the one-dimensional Osher-type solvers introduced in [44, 43]. It consists in a central part of the flux in which simply the arithmetic averages of the flux tensors surrounding a vertex is computed. The multidimensional numerical dissipation is then achieved by computing a surface integral of the matrix absolute value operator applied to the flux Jacobians in each spatial direction over a virtual simplex element constituted by the cell centers around each vertex. The integral is computed numerically via suitable quadrature rules. These integrals are then multiplied by a discrete multidimensional gradient that is simply obtained via the Gauss theorem, as already proposed in [21] to obtain discrete gradients on general unstructured meshes.

2. the careful rewriting of well-known fluctuations from the residual distribution (RD) framework allows to construct suitable approximate multidimensional Riemann solvers that can then be incorporated inside classical high order cell-centered WENO finite volume schemes. It is therefore very interesting to note that RD fluctuations are not limited to the RD context, but can be used in a much more general setting. The most prominent example of genuinely multidimensional upwinding that accounts for all characteristic fields present in the nonlinear hyperbolic system is perhaps the N scheme. Its rewriting as vertex flux in the context of multidimensional Riemann solvers has been achieved in Equations (39)-(42) in combination with the general definition (27) that relates existing RD fluctuations to numerical point fluxes.

Future work will concern both a more in-depth analysis of our approach for the extension to higher than second order of accuracy and the application of these new multidimensional Riemann solvers in the context of exactly divergence-free schemes for MHD and more complex hyperbolic PDE systems, as well as the extension to non-conservative hyperbolic systems.

## 9. Acknowledgments

E. Gaburro gratefully acknowledges the support received from the European Union with the ERC Starting Grant ALcHyMiA (grant agreement No. 101114995).

M. Dumbser was funded by the Italian Ministry of Education, University and Research (MIUR) in the framework of the PRIN 2022 project *High order structure-preserving semi-implicit schemes for hyperbolic equations*, via the Departments of Excellence Initiative 2018–2027 attributed to DICAM of the University of Trento (grant L. 232/2016) and by the European Union Next Generation EU projects PNRR Spoke 7 CN HPC and SMART. M. Dumbser has also received funding by the Fondazione Caritro via the project SOPHOS.

Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.
E. Gaburro and M. Dumbser are members of the INdAM GNCS group in Italy. M. Ricchiuto is member of the CARDAMOM team at the Inria center of the University of Bordeaux.

## In memoriam

This paper is dedicated to the memory of Prof. Arturo Hidalgo López (*July 3rd 1966 - †August 26th 2024) of the Universidad Politecnica de Madrid, organizer of HONOM 2019, active participant in many other editions of HONOM and, most of all, a very dear friend. Our thoughts and wishes go to his wife Lourdes and his sister María Jesús, whom he left behind.

## Appendix A. Entropy conservation: Tadmor's shuffle relation and point entropy fluctuation

We now consider the entropy evolution for the multidimensional scheme with corner fluxes

$$|\Omega_c| \frac{d\mathbf{Q}_c}{dt} + \sum_p \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} = 0.$$

Let $(\eta(\mathbf{Q}), \mathbf{G}(\mathbf{Q}))$ denote the entropy/entropy flux pair, with as usual

$$\mathbf{G}(\mathbf{Q}) = \mathbf{W}^t \hat{\mathbf{F}}(\mathbf{Q}) - \mathbf{\Psi}(\mathbf{Q})$$

with $\mathbf{W}^t = \nabla_Q \eta$ the entropy variables, and $\mathbf{\Psi}$ the entropy potential vector.

We first generalize Tadmor's shuffle condition. To this end we compute

$$\sum_c |\Omega_c| \mathbf{W}_c^t \frac{d\mathbf{Q}_c}{dt} + \sum_c \sum_{p \in \mathcal{P}_c} \mathbf{W}_c^t \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} = 0$$

with $\mathbf{W}_c = \mathbf{W}(\mathbf{Q}_c)$ the entropy variables vector. By definition we have

$$\sum_c |\Omega_c| \frac{d\eta(\mathbf{Q}_c)}{dt} + \sum_c \sum_{p \in \mathcal{P}_c} \mathbf{W}_c^t \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} = 0$$

Following [95] we now write at a given $p \in \mathcal{P}_c$

$$\mathbf{W}_c = \overline{\mathbf{W}}_p + (\mathbf{W}_c - \overline{\mathbf{W}}_p) = \overline{\mathbf{W}}_p + \frac{1}{d+1} \sum_{\ell \in C_p} (\mathbf{W}_c - \mathbf{W}_\ell)$$

where $\overline{\mathbf{W}}_p = \sum_{\ell \in C_p} \mathbf{W}_\ell / (d+1)$. As a consequence we have

$$\sum_c |\Omega_c| \frac{d\eta(\mathbf{Q}_c)}{dt} + \sum_c \sum_{p \in \mathcal{P}_c} \overline{\mathbf{W}}_p^t \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} + \sum_c \sum_{p \in \mathcal{P}_c} (\mathbf{W}_c - \overline{\mathbf{W}}_p)^t \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} = 0$$

Assume now that $\forall c \in C_p$ we have

$$\sum_{c \in C_p} \mathbf{W}_c \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} = \sum_{c \in C_p} \mathbf{\Psi}_c \cdot \mathbf{n}_{pc} \qquad (A.1)$$

which is equivalent to

$$\sum_{c \in C_p} (\mathbf{W}_c - \overline{\mathbf{W}}_p)^t \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} = \sum_{c \in C_p} (\mathbf{\Psi}_c - \overline{\mathbf{\Psi}}_p) \cdot \mathbf{n}_{pc} \qquad (A.2)$$
Under this condition we can write

$$ \sum_c |\Omega_c| \frac{d\eta(Q_c)}{dt} + \sum_c \sum_{p \in \mathcal{P}_c} (\overline{\mathbf{W}}_p^t \hat{\mathbf{F}}_p - \overline{\mathbf{\Psi}}_p) \cdot \mathbf{n}_{pc} = - \sum_c \mathbf{\Psi}_c \cdot \sum_{p \in \mathcal{P}_c} \mathbf{n}_{pc} = 0 $$

which shows that under the generalized shuffle condition (A.1) entropy conservation is verified with the multidimensional numerical entropy flux

$$ \hat{\mathbf{G}}_p := \overline{\mathbf{W}}_p^t \hat{\mathbf{F}}_p - \overline{\mathbf{\Psi}}_p. $$

As argued in [95] any scheme with more numerical viscosity than an entropy conservative scheme will be entropy stable.

As an alternative approach we now consider the schemes in fluctuation form. In absence of reconstruc- tion, and proceeding as before this gives

$$ \sum_c |\Omega_c| \frac{d\eta(Q_c)}{dt} + \sum_c \sum_{p \in \mathcal{P}_c} \mathbf{W}_c^t \boldsymbol{\phi}_{pc} = 0 $$

We can add the sum of the internal entropy value flux along the nodal normals, which is null due to the properties of the normals:

$$ \sum_c |\Omega_c| \frac{d\eta(Q_c)}{dt} + \sum_c \sum_{p \in \mathcal{P}_c} \mathbf{G}_c \cdot \mathbf{n}_{pc} + \sum_c \sum_{p \in \mathcal{P}_c} \mathbf{W}_c^t \boldsymbol{\phi}_{pc} = 0 $$

with $\mathbf{G}_c = \mathbf{G}(Q_c)$. We can easily see that a scheme is entropy conservative if for any corner $p$ we have

$$ \sum_{c \in C_p} \mathbf{W}_c^t \boldsymbol{\phi}_{pc} = \mathbf{\Phi}_\eta := \sum_{c \in C_p} \mathbf{G}_c \cdot \mathbf{n}_{cp} \qquad (A.3) $$

with associated numerical entropy flux verifying

$$ \hat{\mathbf{G}}_p \cdot \mathbf{n}_{pc} = \mathbf{G}_c \cdot \mathbf{n}_{pc} + \mathbf{W}_c^t \boldsymbol{\phi}_{pc} $$

**Proposition 1** (Equivalence of the entropy condition in potential/fluctuation form). *The potential form of the shuffle condition (A.1) and the fluctuation condition (A.3) are equivalent.*

*Proof.* We prove that one condition implies the other and vice-versa. This first implication be easily shown by a direct computation, summing up the shuffle conditions around node $p$:

$$ \sum_{c \in C_p} (\mathbf{W}_c - \overline{\mathbf{W}}_p)^t \mathbf{F}_c \cdot \mathbf{n}_{pc} + \sum_{c \in C_p} (\mathbf{W}_c - \overline{\mathbf{W}}_p)^t \boldsymbol{\phi}_{pc} = \sum_{c \in C_p} (\mathbf{\Psi}_c - \overline{\mathbf{\Psi}}_p) \cdot \mathbf{n}_{pc} $$

$$ \sum_{c \in C_p} \mathbf{W}_c^t \mathbf{F}_c \cdot \mathbf{n}_{pc} + \overline{\mathbf{W}}_p^t \boldsymbol{\phi}_p + \sum_{c \in C_p} \mathbf{W}_c^t \boldsymbol{\phi}_{pc} - \overline{\mathbf{W}}_p^t \boldsymbol{\phi}_p = \sum_{c \in C_p} \mathbf{\Psi}_c \cdot \mathbf{n}_{pc} $$

which readily shows that

$$ \sum_{c \in C_p} \mathbf{W}_c^t \boldsymbol{\phi}_{pc} = - \sum_{c \in C_p} \mathbf{W}_c^t \mathbf{F}_c \cdot \mathbf{n}_{cp} - \sum_{c \in C_p} \mathbf{\Psi}_c \cdot \mathbf{n}_{cp} = \sum_{c \in C_p} \mathbf{G}_c \cdot \mathbf{n}_{cp} $$
The reverse implication is also easy to show using the definition of the $\phi_{pc}$:

$$ 
\begin{align*}
\sum_{c \in C_p} \mathbf{W}_c^t (\hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} - \mathbf{F}_c \cdot \mathbf{n}_{pc}) &= \sum_{c \in C_p} \mathbf{G}_c \cdot \mathbf{n}_{cp} \\
\sum_{c \in C_p} \mathbf{W}_c^t (\hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} - \mathbf{F}_c \cdot \mathbf{n}_{pc}) &= - \sum_{c \in C_p} (\mathbf{W}_c^t \mathbf{F}_c \cdot \mathbf{n}_{pc} - \mathbf{\Psi}_c \cdot \mathbf{n}_{pc}) \\
\sum_{c \in C_p} \mathbf{W}_c^t \hat{\mathbf{F}}_p \cdot \mathbf{n}_{pc} &= \sum_{c \in C_p} \mathbf{\Psi}_c \cdot \mathbf{n}_{pc}
\end{align*}
$$

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

## References

[1] R. Abgrall. On essentially non-oscillatory schemes on unstructured meshes: analysis and implementation. *J. Comput. Phys.*, 144:45–58, 1994.

[2] R. Abgrall. High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices. *J. Sci. Comput.*, 73:461–494, 2017.

[3] R. Abgrall. A personal discussion on conservation, and how to formulate it. In Emmanuel Franck, Jürgen Fuhrmann, Victor Michel-Dansac, and Laurent Navoret, editors, *Finite Volumes for Complex Applications X—Volume 1, Elliptic and Parabolic Problems*, pages 3–19, Cham, 2023. Springer Nature Switzerland.

[4] R. Abgrall and T.J. Barth. Residual distribution schemes for conservation laws via adaptive quadrature. *SIAM J. Sci. Comput.*, 24(3):732–769, 2002.

[5] R. Abgrall, A. Larat, and M. Ricchiuto. Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes. *Journal of Computational Physics*, 230(11):4103–4136, 2011.

[6] R. Abgrall and M. Mezine. Construction of second-order accurate monotone and stable residual distribution schemes for unsteady flow problems. *J. Comput. Phys.*, 188:16–55, 2003.

[7] R. Abgrall, P. Öffner, and H. Ranocha. Reinterpretation and extension of entropy correction terms for residual distribution and discontinuous Galerkin schemes: Application to structure preserving discretization. *Journal of Computational Physics*, 453:110955, 2022.

[8] R. Abgrall and M. Ricchiuto. Hyperbolic balance laws: residual distribution, local and global fluxes. *Numerical Fluid Dynamics*, pages 177–222, 2022.

[9] R. Abgrall and Mario Ricchiuto. High-order methods for cfd. In *Encyclopedia of Computational Mechanics Second Edition*, pages 1–54. John Wiley & Sons, Ltd, 2017.

[10] R. Abgrall and P.L. Roe. High-order fluctuation schemes on triangular meshes. *J. Sci. Comput.*, 19(3):3–36, 2003.

[11] Remi Abgrall. A general framework to construct schemes satisfying additional conservation relations. Application to entropy conservative and entropy dissipative schemes. *J. Comput. Phys.*, 372:640–666, 2018.

[12] Rémi Abgrall, Wasilij Barsukow, and Christian Klingenberg. The Active Flux method for the Euler equations on Cartesian grids. *arXiv preprint arXiv:2310.00683*, 2023.

[13] Rémi Abgrall, Wasilij Barsukow, and Christian Klingenberg. A Semi-discrete Active Flux Method for the Euler Equations on Cartesian Grids. *Journal of Scientific Computing*, 102(2):36, 2025.

[14] D.S. Balsara. Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows. *J. Comput. Phys.*, 229:1970–1993, 2010.

[15] D.S. Balsara. A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows. *J. Comput. Phys.*, 231:7476–7503, 2012.

[16] D.S. Balsara. Multidimensional Riemann problem with self-similar internal structure – Part I – Application to hyperbolic conservation laws on structured meshes. *J. Comput. Phys.*, 277:163–200, 2014.

[17] D.S. Balsara and M. Dumbser. Divergence-free MHD on unstructured meshes using high order finite volume schemes based on multidimensional Riemann solvers. *Journal of Computational Physics*, 299:687–715, 2015.

[18] D.S. Balsara and M. Dumbser. Multidimensional Riemann problem with self-similar internal structure – Part II – Application to hyperbolic conservation laws on unstructured meshes. *J. Comput. Phys.*, 287:269–292, 2015.

[19] Wasilij Barsukow, Jonathan Hohm, Christian Klingenberg, and Philip L Roe. The active flux scheme on Cartesian grids and its low Mach number limit. *Journal of Scientific Computing*, 81:594–622, 2019.

[20] T.J. Barth and P.O. Frederickson. Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction. *AIAA paper no. 90-0013*, 28th Aerospace Sciences Meeting January 1990.

[21] T.J. Barth and D.C. Jespersen. The design and application of upwind schemes on unstructured meshes. *AIAA Paper 89-0366*, pages 1–12, 1989.

[22] W. Boscheri, M. Dumbser, and D.S. Balsara. High Order Lagrangian ADER-WENO Schemes on Unstructured Meshes – Application of Several Node Solvers to Hydrodynamics and Magnetohydrodynamics. *International Journal for Numerical Methods in Fluids*, 76:737–778, 2014.
[23] S. Busto, S. Chiocchetti, M. Dumbser, E. Gaburro, and I. Peshkov. High order ADER schemes for continuum mechanics. *Frontiers in Physics*, 8:32, 2020.

[24] G. Carré, S. Del Pino, B. Després, and E. Labourasse. A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension. *Journal of Computational Physics*, 228:5160–5183, 2009.

[25] M. Ciallella, E. Gaburro, M. Lorini, and M. Ricchiuto. Shifted boundary polynomial corrections for compressible flows: high order on curved domains using linear meshes. *Applied Mathematics and Computation*, 441:127698, 2023.

[26] P. Colella. Multidimensional upwind methods for hyperbolic conservation laws. *Journal of Computational Physics*, 87:171–200, 1990.

[27] I. Cravero, G. Puppo, M. Semplice, and G. Visconti. Cweno: uniformly accurate reconstructions for balance laws. *Math. Comput.*, 87(312):1689–1719, 2018.

[28] Á. Csík, M. Ricchiuto, and H. Deconinck. A Conservative Formulation of the Multidimensional Upwind Residual Distribution Schemes for General Nonlinear Conservation Laws. *J. Comput. Phys.*, 179(2):286–312, 2002.

[29] H. Deconinck. Beyond the riemann problem, part ii. In M. Y. Hussaini, A. Kumar, and M. D. Salas, editors, *Algorithmic Trends in Computational Fluid Dynamics*, pages 341–367, New York, NY, 1993. Springer New York.

[30] H. Deconinck, P.L. Roe, and R. Struijs. A multidimensional generalization of Roe's flux difference splitter for the Euler equations. *Computers and Fluids*, 22:215–222, 1993.

[31] H. Deconinck, P.L. Roe, and R. Struijs. A multidimensional generalization of Roe's difference splitter for the Euler equations. *Computers and Fluids*, 22(2/3):215–222, 1993.

[32] Herman Deconinck and Mario Ricchiuto. Residual distribution schemes: Foundations and analysis. In *Encyclopedia of Computational Mechanics Second Edition*, pages 1–53. John Wiley & Sons, Ltd, 2017.

[33] A. Del Grosso, M.J. Castro, A. Chan, G. Gallice, R. Loubère, and P.H. Maire. A well-balanced, positive, entropy-stable, and multi-dimensional-aware finite volume scheme for 2d shallow-water equations with unstructured grids. *Journal of Computational Physics*, 503:112829, 2024.

[34] B. Després and C. Mazeran. Lagrangian gas dynamics in two-dimensions and Lagrangian systems. *Archive for Rational Mechanics and Analysis*, 178:327–372, 2005.

[35] M. Dumbser D.S. Balsara and R. Abgrall. Multidimensional HLLC Riemann Solver for Unstructured Meshes - With Application to Euler and MHD Flows. *J. Comput. Phys.*, 261:172–208, 2014.

[36] M. Dumbser and D. S. Balsara. A new efficient formulation of the HLLEM Riemann solver for general conservative and non-conservative hyperbolic systems. *Journal of Computational Physics*, 304(C):275–319, 2016.

[37] M. Dumbser, D.S. Balsara, M. Tavelli, and F. Fambri. A divergence-free semi-implicit finite volume scheme for ideal, viscous and resistive magnetohydrodynamics. *International Journal for Numerical Methods in Fluids*, 89:16–42, 2019.

[38] M. Dumbser, D.S. Balsara, E.F. Toro, and C.-D. Munz. A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes. *Journal of Computational Physics*, 227(18):8209–8253, 2008.

[39] M. Dumbser, W. Boscheri, M. Semplice, and G. Russo. Central weighted eno schemes for hyperbolic conservation laws on fixed and moving unstructured meshes. *SIAM J. Sci. Comput.*, 39(6):A2564–A2591, 2017.

[40] M. Dumbser, C. Eaux, and E.F. Toro. Finite Volume Schemes of Very High Order of Accuracy for Stiff Hyperbolic Balance Laws. *J. Comput. Phys.*, 227:3971–4001, 2008.

[41] M. Dumbser and M. Käser. Arbitrary high order non-oscillatory Finite Volume schemes on unstructured meshes for linear hyperbolic systems. *J. Comput. Phys.*, 221:693–723, 2007.

[42] M. Dumbser, M. Käser, V.A Titarev, and E.F. Toro. Quadrature-Free Non-Oscillatory Finite Volume Schemes on Unstructured Meshes for Nonlinear Hyperbolic Systems. *J. Comput. Phys.*, 226:204–243, 2007.

[43] M. Dumbser and E. F. Toro. On Universal Osher-Type Schemes for General Nonlinear Hyperbolic Conservation Laws. *Communications in Computational Physics*, 10:635–671, 2011.

[44] M. Dumbser and E. F. Toro. A simple extension of the Osher Riemann solver to non-conservative hyperbolic systems. *J. Sci. Comput.*, 48:70–88, 2011.

[45] F. Fambri, E. Zampa, S. Busto, L. Río-Martín, F. Hindenlang, E. Sonnendrücker, and M. Dumbser. A well-balanced and exactly divergence-free staggered semi-implicit hybrid finite volume / finite element scheme for the incompressible MHD equations. *Journal of Computational Physics*, 493:112493, 2023.

[46] E. Gaburro. A unified framework for the solution of hyperbolic PDE systems using high order direct Arbitrary-Lagrangian-Eulerian schemes on moving unstructured meshes with topology change. *Archives of Computational Methods in Engineering*, 28(3):1249–1321, 2021.

[47] E. Gaburro, W. Boscheri, S. Chiocchetti, C. Klingenberg, V. Springel, and M. Dumbser. High order direct arbitrary-lagrangian-eulerian schemes on moving voronoi meshes with topology changes. *Journal of Computational Physics*, 407:109167, 2020.

[48] Elena Gaburro, Walter Boscheri, Simone Chiocchetti, and Mario Ricchiuto. Discontinuous Galerkin schemes for hyperbolic systems in non-conservative variables: quasi-conservative formulation with subcell finite volume corrections. *Computer Methods in Applied Mechanics and Engineering*, 431:117311, 2024.

[49] G. Gallice, A. Chan, R. Loubère, and P.H. Maire. Entropy stable and positivity preserving godunov-type schemes for multidimensional hyperbolic systems on unstructured grid. *Journal of Computational Physics*, 468:111493, 2022.

[50] J. Garicano-Mena, A. Lani, and H. Deconinck. An energy-dissipative remedy against carbuncle: Application to hypersonic flows around blunt bodies. *Computers and Fluids*, 133:43–54, 2016.

[51] S.K. Godunov. Finite Difference Methods for the Computation of Discontinuous Solutions of the Equations of Fluid Dynamics. *Mathematics of the USSR: Sbornik*, 47:271–306, 1959.
[52] A. Harten, B. Engquist, S. Osher, and S.R. Chakravarthy. Uniformly High Order Accurate Essentially Non-oscillatory Schemes III. *J. Comput. Phys.*, 71:231–303, 1987.

[53] A. Harten, P. D. Lax, and B. van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. *SIAM Review*, 25(1):35–61, 1983.

[54] C. Hu and C.-W. Shu. Weighted essentially non-oscillatory schemes on triangular meshes. *Journal of Computational Physics*, 150(1):97–127, 1999.

[55] C. Hu and C.W. Shu. A high-order WENO finite difference scheme for the equations of ideal magnetohydrodynamics. *Journal of Computational Physics*, 150:561–594, 1999.

[56] H. Jackson. On the eigenvalues of the ADER-WENO Galerkin predictor. *J. Comput. Phys.*, 333:409–413, 2017.

[57] G.S. Jiang and C.W. Shu. Efficient implementation of weighted ENO schemes. *J. Comput. Phys.*, 126:202–228, 1996.

[58] F. Kemm, E. Gaburro, F. Thein, and M. Dumbser. A simple diffuse interface approach for compressible flows around moving solids of arbitrary shape based on a reduced Baer–Nunziato model. *Computers & fluids*, 204:104536, 2020.

[59] J. Kepler. *Nova stereometria doliorium vinariorum, in primis Austriaci*. Plancus, Lincii, Austria, 1615.

[60] A. Lani, M. Panesi, and H. Deconinck. Conservative Residual Distribution Method for Viscous Double Cone Flows in Thermochemical Nonequilibrium. *Communications in Computational Physics*, 13(2):479–501, 2013.

[61] P. Lax and X.D. Liu. Solution of two-dimensional Riemann problems of gas dynamics by positive schemes. *SIAM Journal of Scientific Computing*, 19:319–340, 1998.

[62] P.D. Lax. Weak solutions of nonlinear hyperbolic equations and their numerical approximation. *Comm. Pure Appl. Math.*, 7:159–193, 1954.

[63] D. Levy, G. Puppo, and G. Russo. Central WENO schemes for hyperbolic systems of conservation laws. *M2AN Math. Model. Numer. Anal.*, 33(3):547–571, 1999.

[64] D. Levy, G. Puppo, and G. Russo. A third order central WENO scheme for 2D conservation laws. *Applied Numerical Mathematics*, 33:415–421, 2000.

[65] D. Levy, G. Puppo, and G. Russo. A fourth-order central WENO scheme for multidimensional hyperbolic systems of conservation laws. *SIAM J. Sci. Comput.*, 24:480–506, 2002.

[66] M. Lukacova-á-Medvidova, G. Warnecke, and Y. Zahaykah. Finite volume evolution Galerkin (FVEG) methods for three-dimensional wave equation system. *Applied Numerical Mathematics*, 57:1050–1064, 2007.

[67] P.-H. Maire. A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes. *J. Comput. Phys.*, 228:2391–2425, 2009.

[68] P.-H. Maire, R. Abgrall, J. Breil, and J. Ovadia. A cell-centered lagrangian scheme for two-dimensional compressible flow problems. *SIAM J. Sci. Comput.*, 29:1781–1824, 2007.

[69] P.-H. Maire and J. Breil. A second-order cell-centered lagrangian scheme for two-dimensional compressible flow problems. *Int. J. Numer. Methods Fluids*, 56:1417–1423, 2007.

[70] P.H. Maire. A high-order one-step sub-cell force-based discretization for cell-centered Lagrangian hydrodynamics on polygonal grids. *Computers and Fluids*, 46(1):341–347, 2011.

[71] P.H. Maire. A unified sub-cell force-based discretization for cell-centered Lagrangian hydrodynamics on polygonal grids. *International Journal for Numerical Methods in Fluids*, 65:1281–1294, 2011.

[72] H. Nishikawa, M. Rad, and P.L. Roe. A third-order fluctuation splitting scheme that preserves potential flow. 15th AIAA Computational Fluid Dynamics Conference, Anaheim, CA, USA, June 2001.

[73] S. Osher and F. Solomon. Upwind Difference Schemes for Hyperbolic Conservation Laws. *Math. Comput.*, 38:339–374, 1982.

[74] M. Ricchiuto. *Contributions to the development of residual discretizations for hyperbolic conservation laws with application to shallow water flows*. Habilitation à diriger des recherches, Université Sciences et Technologies - Bordeaux I, December 2011. PDF available at https://theses.hal.science/tel-00651688/.

[75] M. Ricchiuto. *Contributions to the development of residual discretizations for hyperbolic conservation laws with application to shallow water flows*. Habilitation à diriger des recherches, Université Sciences et Technologies - Bordeaux I, December 2011. PDF available at https://theses.hal.science/tel-00651688/.

[76] M. Ricchiuto and R. Abgrall. Explicit runge-kutta residual distribution schemes for time dependent problems: Second order case. *Journal of Computational Physics*, 229(16):5653–5691, 2010.

[77] M. Ricchiuto, R. Abgrall, and H. Deconinck. Application of conservative residual distribution schemes to the solution of the shallow water equations on unstructured meshes. *J. Comput. Phys.*, 222:287–331, 2007.

[78] M. Ricchiuto, Á. Csík, and H. Deconinck. Residual distribution for general time dependent conservation laws. *J. Comput. Phys*, 209(1):249–289, 2005.

[79] A.V. Rodionov. Artificial viscosity in godunov-type schemes to cure the carbuncle phenomenon. *Journal of Computational Physics*, 345:308–329, 2017.

[80] P. L. Roe. Linear advection schemes on triangular meshes. Technical Report CoA 8720, Cranfield Institute of Technology, 1987.

[81] P. L. Roe. "Optimum" upwind advection on a triangular mesh. Technical Report ICASE 90-75, ICASE, 1990.

[82] P. L. Roe and D. Sidilkover. Optimum positive linear schemes for advection in two and three dimensions. *SIAM Journal on Numerical Analysis*, 29(6):1542–1568, 1992.

[83] P.L. Roe. Approximate Riemann Solvers, Parameter vectors, and Difference schemes. *J. Comput. Phys.*, 43:357–372, 1981.

[84] P.L. Roe. Fluctuations and signals - a framework for numerical evolution problems. In K.W. Morton and M.J. Baines, editors, *Numerical Methods for Fluids Dynamics*, pages 219–257. Academic Press, 1982.
[85] P.L. Roe. Discrete models for the numerical analysis of time-dependent multidimensional gas dynamics. *Journal of Computational Physics*, 63:458-476, 1986.

[86] P.L. Roe. Beyond the riemann problem, part i. In M. Y. Hussaini, A. Kumar, and M. D. Salas, editors, *Algorithmic Trends in Computational Fluid Dynamics*, pages 341–367, New York, NY, 1993. Springer New York.

[87] P.L. Roe. Musings of a computational philosopher. In James C. Tyacke and Nagabhushana Rao Vadlamani, editors, *Proceedings of the Cambridge Unsteady Flow Symposium 2024*, pages 1–35, Cham, 2025. Springer Nature Switzerland.

[88] C.L. Rumsey, B. van Leer, and P.L. Roe. A multidimensional flux function with applications to the Euler and Navier-Stokes equations. *Journal of Computational Physics*, 105:306–323, 1993.

[89] J. Saltzman. An unsplit 3D upwind method for hyperbolic conservation laws. *Journal of Computational Physics*, 115:153–168, 1994.

[90] K. A Schneider, J.M. Gallardo, D.S. Balsara, B. Nkonga, and C. Parés. Multidimensional approximate Riemann solvers for hyperbolic nonconservative systems. Applications to shallow water systems. *Journal of Computational Physics*, 444:110547, 2021.

[91] C. W. Schulz-Rinne. Classification of the Riemann problem for two-dimensional gas dynamics. *SIAM J. Math. Anal.*, 24:76–88, 1993.

[92] C.W. Schulz-Rinne, J.P. Collins, and H.M. Glaz. Numerical solution of the Riemann problem for two-dimensional gasdynamics. *SIAM Journal on Scientific Computing*, 14:1394–1414, 1993.

[93] K. Sermeus and H. Deconinck. An entropy fix for multidimensional upwind residual distribution schemes. *Computers and Fluids*, 34(4):617–640, 2005.

[94] A.H. Stroud. *Approximate Calculation of Multiple Integrals*. Prentice-Hall Inc., Englewood Cliffs, New Jersey, 1971.

[95] E. Tadmor. Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems. *Acta Numerica*, 12:451–512, 2003.

[96] V.A. Titarev and E.F. Toro. ADER: Arbitrary high order Godunov approach. *Journal of Scientific Computing*, 17(1-4):609–618, December 2002.

[97] V.A. Titarev, P. Tsoutsanis, and D. Drikakis. WENO schemes for mixed-element unstructured meshes. *Communications in Computational Physics*, 8:585–609, 2010.

[98] E. F. Toro. *Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction, Third edition*. Springer-Verlag, Berlin, 2009.

[99] E.F. Toro, M. Spruce, and W. Speares. Restoration of the contact surface in the Harten-Lax-van Leer Riemann solver. *Journal of Shock Waves*, 4:25–34, 1994.

[100] E.F. Toro and V.A. Titarev. ADER schemes for scalar non-linear hyperbolic conservation laws with source terms in three-space dimensions. *Journal of Computational Physics*, 202(1):196–215, 2005.

[101] P. Tsoutsanis, V.A. Titarev, and D. Drikakis. WENO schemes on arbitrary mixed-element unstructured meshes in three space dimensions. *J. Comput. Phys.*, 230:1585–1601, 2011.

[102] E. van der Weide and H. Deconinck. Positive matrix distribution schemes for hyperbolic systems. In *Computational Fluid Dynamics*, pages 747–753, New York, 1996. Wiley.

[103] E. Zampa, S. Busto, and M. Dumbser. A divergence-free hybrid finite volume / finite element scheme for the incompressible mhd equations based on compatible finite element spaces. *Applied Numerical Mathematics*, 198:346–374, 2024.

[104] Y.T. Zhang and C.W. Shu. Third Order WENO Scheme on Three Dimensional Tetrahedral Meshes. *Communications in Computational Physics*, 5:836–848, 2009.