Jeaniffer Vides

Boniface Nkonga

Edouard Audit

A Simple Two-Dimensional Extension of the HLL Riemann Solver for Gas Dynamics

Keywords: Multidimensional Riemann solvers, Godunov-type scheme, conservation laws, gas dynamics, aerodynamics Solveurs de Riemann multidimensionnel, schémas de type Godunov, lois de conservation, dynamique des gaz, aérodynamique

We report on our study aimed at deriving a simple method to numerically approximate the solution of the two-dimensional Riemann problem for gas dynamics, using the literal extension of the well-known HLL formalism as its basis. Essentially, any strategy attempting to extend the three-state HLL Riemann solver to multiple space dimensions will by some means involve a piecewise constant approximation of the complex two-dimensional interaction of waves, and our numerical scheme is not the exception. In order to determine closed form expressions for the involved fluxes, we rely on the equivalence between the consistency condition and the use of Rankine-Hugoniot conditions that hold across the outermost planar waves emerging from the Riemann problem's initial discontinuities. The proposed scheme is then carefully designed to simplify its eventual numerical implementation and its advantages are analytically attested. We also present first numerical results that put into evidence its robustness and stability.

Introduction

In the context of standard finite volume approximations, the conservative form of the governing equations is especially important when dealing with problems admitting discontinuities (such as shocks and contact discontinuities) in the solution. Non-conservative formulations lead not only to incorrect estimations of wave speeds, but also to inconsistency of the numerical approximation with the weak form [START_REF] Hou | Why nonconservative schemes converge to wrong solutions: Error analysis[END_REF], whereas convergent conservative approximations are known to be always consistent with the weak formulation [START_REF] Lax | Systems of conservation laws[END_REF] as Rankine-Hugoniot relations are satisfied.

It has been pointed out, since the early fifties [START_REF] Richtmyer | Difference methods for initial-value problems[END_REF][START_REF] Vonneumann | A method for the numerical calculation of hydrodynamic shocks[END_REF], that artificial viscosity is unavoidable when designing convergent approximations for hyperbolic systems, and in this context, we mention the pioneering work of Godunov [START_REF] Godunov | A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF], who formulated a conservative scheme that employs the solutions of Riemann problems at cell interfaces to approximate the local numerical flux. Indeed, it is well known that numerical dissipation is the effect of fine unresolved scales on the coarse resolved ones (refer to the variational multiscale, abbreviated VMS, framework), and interface Riemann problems are ways to describe the evolution of these fine scales for given coarse scale data.

For the one-dimensional case, the previously mentioned interfaces are simply vertices connecting two coarse cells and the corresponding Riemann problem can be solved exactly for the Euler equations with the ideal gas equation of state, even though approximated solutions are more commonly used in practice. Thanks to strict hyperbolicity and entropy dissipation in shocks, a key concept used in the one-dimensional analysis is the fact that weak solutions evolve in time toward a non-interacting scattering state (decay of Glimm's interaction potential, see [START_REF] Glimm | Decay of solutions of systems of hyperbolic conservation laws[END_REF][START_REF] Liu | Well-posedness theory for hyperbolic conservation laws[END_REF]). Then, consider that in two dimensions the interfaces are defined as collections of vertices connected by edges, and in three dimensions, by edges and surfaces. The Riemann problems at the vertices are genuinely multidimensional (MultiD), involving interactions of more than two coarse data, whereas the Riemann problems associated with the cell edges are locally one-dimensional (in the 2D situations). In principle, a proper multidimensional approximation should take into account interactions of both 1D and MultiD Riemann problems. Unfortunately, even for the Euler equations with the ideal gas equation of state, an adequate approximation of multidimensional Riemann problems is a challenging problem [START_REF] Abgrall | Approximation du problème de Riemann vraiment multidimensionnel des équations d'Euler par une méthode de type Roe (I): La linéarisation[END_REF][START_REF] Abgrall | Approximation du problème de Riemann vraiment multidimensionnel des équations d'Euler par une méthode de type Roe (II): Solution du problème de Riemann approché[END_REF], especially given the complexity of the nonlinear interaction of wave patterns [START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF][START_REF] Li | Characteristic decompositions and interactions of rarefaction waves of 2-D Euler equations[END_REF][START_REF] Zhang | Conjecture on the structure of solutions of the Riemann problem for two-dimensional gas dynamics systems[END_REF]. As a matter of fact, the corresponding Mul-tiD solutions do not systematically provide a non-interacting scattering state because of these complicated effects.

Put in the simplified context of the Euler equations with the ideal equation of state, for the two-dimensional Riemann problem, even under the premise that each jump between neighboring initial states projects one planar wave consisting of a single shock, rarefaction or contact discontinuity, the number of allowable distinct self-similar configurations can be up to 77 [START_REF] Schulz-Rinne | Numerical solution of the Riemann problem for two-dimensional gas dynamics[END_REF]. Moreover, non-self-similar bounded admissible solutions exist and recent results rigorously show the ill-posedness of entropy solutions for rather elementary multidimensional systems. It is then unrealistic to expect that an exact MultiD Riemann solver can be used as a building block for numerical schemes. Yet, we can still view the numerical flux as the one-dimensional flux across cell boundaries plus multidimensional corrections emanating from the corners (in the 2D case). In the Lagrangian context, the approximation of corner interactions has already been combined with one-dimensional fluxes to obtain robust approximations satisfying a discrete entropy inequality [START_REF] Després | Lagrangian gas dynamics in two dimensions and Lagrangian systems[END_REF][START_REF] Maire | A cell-centered Lagrangian scheme for two-dimensional compressible flow problems[END_REF][START_REF] Maire | Multi-scale Godunov-type method for cell-centered discrete Lagrangian hydrodynamics[END_REF]. However, it is known that for most of the current flux based numerical strategies, approximations often neglect the corner corrections and only use one-dimensional wave characteristics, even when an operator splitting technique is not employed. We mention that several alternative strategies consisting of Riemann solver-free formulations are available, such as residual distribution schemes [START_REF] Deconinck | Status of multidimensional upwind residual distribution schemes and applications in aeronautics[END_REF][START_REF] Deconinck | Introduction to residual distribution methods[END_REF], VMS methods applied to finite element solutions [START_REF] Hughes | Large eddy simulation and the variational multiscale method[END_REF], and Riemann-solvers-free central schemes [START_REF] Kurganov | Semi-discrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations[END_REF][START_REF] Kurganov | Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers[END_REF], among others, but are not within the scope of this paper.

Our focus here is on Riemann based unsplit formulations taking into account the interactions associated to the "corner" boundaries. A nine state Riemann solver was formulated in [START_REF] Wendroff | A two-dimensional HLLE Riemann solver and associated Godunov-type difference scheme for gas dynamics[END_REF] to obtain numerical approximations that include these interactions as constant states, extending the one-dimensional HLL theory [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] with Einfeldt's wave speed estimates [START_REF] Einfeldt | On Godunov-type methods for gas dynamics[END_REF] to two dimensions. Although his approach includes a valuable interpretation of the approximate structure of 2D solutions at a given time, it regrettably lacks explicit expressions that would enable a direct implementation. One year later, in another line of development, Brio et al. [START_REF] Brio | Two-dimensional Riemann solver for Euler equations of gas dynamics[END_REF] proposed a multistate Riemann solver (defined at the corner) as a linear hyperbolic propagation of acoustic waves, which can be regarded as a partial correction to the 1D solver applied at the interfaces such that the final numerical flux results from a convex combination of purely one-dimensional and corner fluxes. However, the solver was solely developed for the Euler equations of gas dynamics and, given its linear nature, requires considerable reformulation for application to other systems of conservation laws. Recently, Balsara re-examined Wendroff's contribution and formulated a multidimensional solver in [START_REF] Balsara | Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows[END_REF] and a more robust version in [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF], which include, among other things, the closed-form approximate expressions that Wendroff failed to provide. The calculation of the involved states and fluxes at a corner is done by means of the integral form of the conservation laws over a space-time volume that is essentially different from Wendroff's to facilitate the computation of the resulting equations in the subsonic case; to handle supersonic cases, slight modifications must be performed to the fluxes and/or signal speeds, nonetheless.

We propose a formulation that utilizes the advantageous space-time structure proposed by Wendroff [START_REF] Wendroff | A two-dimensional HLLE Riemann solver and associated Godunov-type difference scheme for gas dynamics[END_REF] to devise an approximation that naturally includes all possible flow configurations (subsonic and supersonic). Our solver, suitably built as an extension of the HLL formalism to multidimensions, inevitably leads to the construction of an approximated profile of propagation consisting of constant states and representative of the complexity of the waves associated with the multidimensional Riemann problem. We make use of the consistency with the integral formulation through the Rankine-Hugoniot relations, which hold across planar waves separating these constant states, to derive general closed-form expressions (in the sense of [START_REF] Balsara | Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows[END_REF][START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF]) for the fluxes. All expressions are, in fact, the solution of an overdetermined linear system resolved by the method of ordinary least squares and provide a straightforward implementation of our robust and stable scheme. Although we will restrict our attention to the case of the Euler equations for inviscid compressible gas flows in two space variables, all formalisms developed in this paper can be extended to higher dimensions and applied to any system of conservation laws.

In the subsequent section, we present the general form of a first-order hyperbolic system, particularly the Euler equations for gas dynamics, as well as background information about Riemann solvers (both one-and two-dimensional). In Section 3, we combine ideas from existing solvers with the enforcing of jump conditions, to design a strategy for the construction of simple MultiD Riemann solvers. Next, the developed strategy is validated through applications to test problems in Section 4 and finally, concluding remarks are given in the last section.

Governing Equations and HLL Riemann Solvers

We begin by considering a system of ϑ conservation laws in d spatial dimensions

∂ t w + ∇ • F (w) = 0, in R d × (0, ∞), (2.1a)
Inria with w = (w [1] , . . . , w [ϑ] ) T the state vector of conservative variables and F = (f 1 , . . . , f d ) the flux tensor. The unknown w : R d × [0, ∞) → V is a function from space x ∈ R d and time t to the system's state space V, and each flux in the mth spatial dimension is defined as f m : V → R ϑ , for m = 1, . . . , d. The numerical solution of such system, complemented with initial conditions of the form

w = w 0 , on R d × {t = 0}, (2.1b)
is of considerable interest for modeling diverse physical phenomena, such as in gas dynamics. Due to their known importance in applications, we focus on the Euler equations for inviscid compressible gas flows, given by the following system of nonlinear hyperbolic partial differential equations:

   ∂ t (ρ) + ∇ • (ρu) = 0, ∂ t (ρu) + ∇ • (ρu ⊗ u) + ∇p = 0, ∂ t (ρe) + ∇ • ((ρe + p)u) = 0, (2.2)
where ρ > 0 is the density, u ∈ R d the velocity, and e = ǫ + |u| 2 /2 the specific total energy. In order to close system (2.2), we consider an equation of state of the form

p = p(ρ, ǫ), ( 2.3)
to relate the thermodynamic pressure p with both density and specific internal energy ǫ. Unless stated otherwise, the ideal equation of state p = (γ -1)ρǫ, with the specific heat capacity ratio γ = 1.4, is assumed. Let us note that system (2.2) can be easily be put in compact form (2.1a), by casting the pressure gradient as a divergence, i.e., ∇p = ∇ • (pI) with I the identity matrix. The need of obtaining accurate numerical solutions of the previous equations has led to the development of several schemes. In this section, we aim to describe and review several existing approximate HLL solvers, developed both for the one-and two-dimensional cases.

One-Dimensional Systems

We first consider the nonlinear system of conservation laws (2.2) in one dimension. Denoting the physical flux function by f , we write

∂ t w + ∂ x f (w) = 0, (2.4)
with w = (ρ, ρu, ρe) T and f (w) = (ρu, ρu 2 + p, (ρe + p)u) T . It is well-known that, regardless of the initial data's nature in w 0 , a smooth solution to system (2.4) can break down at a finite time t > 0, such that it no longer satisfies the differential equations in the classical sense. Thus, it is necessary to introduce the associated integral form to account for possible discontinuities. For any control volume of dimension (x a , x b ) × (t 1 , t 2 ), the integral form of (2.4) over it becomes

x b xa w(x, t 2 ) dx = x b xa w(x, t 1 ) dx + t2 t1 f (w(x a , t)) dt - t2 t1 f (w(x b , t)) dt.
(2.5)

Let us now set a uniform numerical mesh with

N x cells C i = (x i-1 2 , x i+ 1 2
) of a determined width ∆x, where x i± 1 2 = x i ± ∆x/2. Henceforth, subscripts refer to spatial location, with cell centers denoted by integer subscripts i = 1, . . . , N x and interfaces denoted by half integers. The time increment is given by ∆t such that t n+1 = t n + ∆t, for n ∈ N.

A direct evaluation of the integral form (2.5) over a rectangle in this computational space-time domain yields

Ci w(x, t n+1 ) dx = Ci w(x, t n ) dx + t n+1 t n f (w(x i-1 2 , t)) dt - t n+1 t n f (w(x i+ 1 2 , t)) dt, (2.6)
which can then be expressed as

w n+1 i = w n i - ∆t ∆x φ i+ 1 2 -φ i-1 2 , ( 2.7)
where w n i is a cell-averaged value of w at time level t n and φ i± 1 2 are time-averaged numerical fluxes at x = x i± 1 2 . Specifically,

w n i = 1 ∆x x i+ 1 2 x i-1 2 w(x, t n ) dx and φ i± 1 2 = 1 ∆t t n+1 t n f (w(x i± 1 2 , t)) dt. (2.8)
We note that the mathematical formulation (2.7) cannot yet be considered as a numerical scheme since we have not specified how to compute φ i± 1 2 . In his seminal paper [START_REF] Godunov | A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF], Godunov introduced a novel first-order numerical approach that forms the basis of numerous interesting schemes, several of them presented in what follows.

Riemann Problems and Godunov-Type Schemes

A Riemann problem centered at x = x 0 is simply an initial-value problem

∂ t w + ∂ x f (w) = 0, w(x, t 0 ) = w l if x < x 0 , w r if x > x 0 , (2.9)
which has a solution that depends only on the initial left and right states, given respectively by w l and w r , and on the value ξ = (x -x 0 )/(t -t 0 ), 0 ≤ t 0 < t. Thus, we denote an exact solution of (2.9) by w(ξ; w l , w r ). For convenience, from this point on, we assume x 0 = 0 and t 0 = 0. Within the finite-volume framework, Godunov's first-order method [START_REF] Godunov | A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF] assumes a piecewise constant distribution of the conserved quantities over each cell, see the first equation in (2.8), and evolves it in time by solving a one-dimensional Riemann problem in the normal direction at each cell interface. Consequently, the value w n+1 i is calculated in terms of the exact solutions of local Riemann problems in the following way:

w n+1 i = 1 ∆x   xi x i-1 2 w x -x i-1 2 ∆t ; w n i-1 , w n i dx + x i+ 1 2 xi w x -x i+ 1 2 ∆t ; w n i , w n i+1 dx   ,
(2.10) which can be written in conservative form (2.7) by defining the intercell numerical fluxes as

φ i-1 2 = f (W(0; w n i-1 , w n i )) and φ i+ 1 2 = f (W(0; w n i , w n i+1 )). (2.11)
For the Euler equations, the main drawback of Godunov's scheme results from computing the exact solution of each nonlinear Riemann problem, which has a direct impact on calculation cost. It is therefore necessary to consider an approximation W(ξ; w l , w r ) to the Riemann problem (2.9), which satisfies the consistency condition

xr x l W(ξ; w l , w r ) dx = x r w r -x l w l + ∆t(f l -f r ), (2.12)
Inria as long as the complicated structure of the exact solution w(ξ; w l , w r ) is contained in the control volume (x l , x r ) × (0, ∆t), x l ≤ 0 ≤ x r . Note that f l = f (w l ) and f r = f (w r ). Using then these approximations W(ξ; w l , w r ), Harten et. al (refer to [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]) define a Godunov-type scheme as

w n+1 i = 1 ∆x xi x i-1 2 W x -x i-1 2 ∆t ; w n i-1 , w n i dx + 1 ∆x x i+ 1 2 xi W x -x i+ 1 2 ∆t ; w n i , w n i+1 dx.
(2.13)

HLL Riemann Solver

One of the simplest Godunov-type schemes is the so-called HLL Riemann solver proposed by Harten et. al [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], where the exact Riemann fan is approximated by two waves containing a single constant state in between (see Figure 1). These waves propagate with speeds s l and s r denoting the smallest and largest signal speeds, the estimation of which will be detailed later in this section. As pointed out in [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], any scheme (2. Thus, considering the previously mentioned configuration, the HLL approximate Riemann solver is defined as

W hll (ξ; w l , w r ) =      w l if ξ < s l , w * if s l < ξ < s r , w r if ξ > s r , (2.15)
with ξ = x/t. The term w * represents the average intermediate state between the two waves and can easily be derived from the conservation laws (2.12), i.e.,

w * = s r w r -s l w l + f l -f r s r -s l . (2.16)
Moreover, we are interested in the determination of the associated numerical flux. For this purpose, by applying the integral conservation laws (2.5) over two distinct rectangles (x l , 0) × (0, ∆t), (0, x r ) × (0, ∆t), two fluxes along the t-axis are obtained

f |l = f l - x l ∆t w l - 1 ∆t 0 x l W hll x ∆t ; w l , w r dx, ( 2.17a
)

f |r = f r - x r ∆t w r + 1 ∆t xr 0 W hll x ∆t ; w l , w r dx, ( 2.17b)
respectively. Clearly, with the equality f |l = f |r , the consistency condition (2.12) is recovered.

Let us now introduce some useful notation. For any constant h ∈ R, we define

h + = max(0, h) and h -= min(0, h), ( 2.18)
recalling that h = h + + h -. Employing this notation, we are able to rewrite (2.17) as

f |l = f l + s l -(w * -w l ), f |r = f r + s r + (w * -w r ), (2.19)
which are both useful to obtain an approximation to the numerical flux along the t-axis. Hence, we substitute the state (2.16), with both signal speeds replaced by s l -and s r + , into any of the previous equations (2.19) to get

φ hll (w l , w r ) = s r + f l -s l -f r + s l -s r + (w r -w l ) s r + -s l - , ( 2.20)
so the scheme can be written in the conservative form (2.7) simply by defining

φ i-1 2 = φ hll (w n i-1 , w n i ) and φ i+ 1 2 = φ hll (w n i , w n i+1 ). (2.21)
Yet another equivalent and simpler way to construct the HLL intermediate state vector and flux relies on applying the Rankine-Hugoniot jump conditions across each of the waves. Specifically, we may think of these intermediate quantities as solutions of the linear system

f * = f l + s l (w * -w l ), (2.22a
)

f * = f r + s r (w * -w r ), (2.22b)
i.e., the Rankine-Hugoniot conditions across the left and right waves, respectively. These conditions hold across curves of discontinuities and are mentioned here briefly for later reference in Section 3. Solving system (2.22) then yields the state w * (2.16) and the HLL flux

f * = s r f l -s l f r + s l s r (w r -w l ) s r -s l . (2.23)
Under the assumption of a subsonic solution where s l < 0 < s r , it is evident that φ hll = f * . With a slight modification to the speeds, we then obtain the intercell flux (2.20). Now, in order to completely determine the numerical fluxes previously described, an adequate choice of the wave speeds s l and s r is needed. In [START_REF] Einfeldt | On Godunov-type methods for gas dynamics[END_REF], Einfeldt derived approximations for the minimum and maximum physical signal velocities of the exact Riemann problem, generalized to

s l = min 1≤θ≤ϑ min λ θ (w l ), λθ (w l , w r )
and s r = max 1≤θ≤ϑ max λ θ (w r ), λθ (w l , w r ) ,

(2.24) where λ θ is the θ-th eigenvalue of the Jacobian matrix ∇ w f associated with system (2.4) and λθ is the θ-th eigenvalue of the Roe matrix [START_REF] Einfeldt | On Godunov-type methods for gas dynamics[END_REF][START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF].

The HLL approach [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] together with Einfeldt's wave speed estimates (2.24) is not only effective and robust but also rather easy to implement. Several details regarding the scheme's ability to preserve the positivity of the internal energy and density throughout the computational process are given in [START_REF] Einfeldt | On Godunov-type methods near low densities[END_REF].

Inria

Two-Dimensional Systems

Due to our specific interest in two-dimensional gas dynamics, we henceforth only consider system (2.2) in d = 2 dimensions with x = (x, y) and denote f and g the fluxes f 1 and f 2 , respectively. We then write

∂ t w + ∂ x f (w) + ∂ y g(w) = 0, (2.25)
with

w =     ρ ρu ρv ρe     , f (w) =     ρu ρu 2 + p ρuv (ρe + p)u     , g(w) =     ρv ρvu ρv 2 + p (ρe + p)v     .
(2.26)

After Harten et al.'s contribution in [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF], several extensions of their HLL scheme have been proposed to find approximate solutions to the above system and in this section we will carefully review two of them. However, following the developments in the one-dimensional case, we first establish the integral form of the conservation laws (2.25). Thus, for all control volumes (x a , x b )× (y c , y d ) × (t 1 , t 2 ), the following integral form should hold: (2.27) As noted before, any finite-volume method based on a Godunov-type approach strongly depends on the exact or approximate solution of the Riemann problem. Conceptually speaking, there exists an infinite number of configurations for a two-dimensional (2D) Riemann problem, as shown later in Section 3.3.1, but for a rectangular mesh, we are only interested in the one centered at x 0 = (x 0 , y 0 ) where the initial data for system (2.25) are piecewise constant such that

w(x, y, t 0 ) =          w sw if x < x 0 , y < y 0 , w se if x > x 0 , y < y 0 , w nw if x < x 0 , y > y 0 , w ne if x > x 0 , y > y 0 .
(2.28)

It is known that system (2.25) is invariant under scaling of the form (x, y, t) → (κx, κy, κt), for any constant κ > 0, and with the initial data given by (2.28), we expect the solution to be a similarity solution of the form w(ξ, η; w ne , w nw , w sw , w se ), hereafter w(ξ, η), i.e., a function that is constant along ξ = (x -x 0 )/(t -t 0 ) and η = (y -y 0 )/(t -t 0 ) and which is self-similar. Despite the reduction in the number of dimensions (from three to two in this case), the solution is expected to be fairly complex nonetheless, as demonstrated by Schulz-Rinne et al. [START_REF] Schulz-Rinne | Numerical solution of the Riemann problem for two-dimensional gas dynamics[END_REF].

To find a numerical approximation, we break the spatial domain into rectangular grid cells with centers indexed as i, j, where i refers to the x-coordinate direction and j to the y-coordinate direction. Here, ∆x = 1/N x and ∆y = 1/N y are the grid spacing such that x i = (i -1)∆x and y j = (j -1)∆y, with i = 1, . . . , N x and j = 1, . . . , N y , and as before, the corresponding cell interfaces are denoted by half integers. Once such a grid has been constructed, we can then define the average value of the gas dynamic state at time level t n over a particular cell

C i,j = (x i-1 2 , x i+ 1 2 ) × (y j-1 2 , y j+ 1 2 ) as w n i,j = 1 |C i,j | Ci,j w(x, y, t n ) dx, ( 2.29)
where |C i,j | = ∆x∆y, and describe the solution at each time step as a piecewise constant function

w h (x, y, t n ) = i=1,...,Nx j=1,...,Ny w n i,j χ i,j (x, y) with χ i,j (x, y) = 1 if (x, y) ∈ C i,j , 0 if (x, y) ∈ C i,j . (2.30)
It is worth noticing that in view of the above-mentioned finite volume grid configuration, the evolution of solutions w h (x, y, t n ) over a time step results from solving a series of one-and twodimensional Riemann problems at the cells' edges and vertices, as depicted in Figure 2(a). Since the one-dimensional theory was already covered in Section 2.1, we now focus on approximately solving the local 2D Riemann problem formed at the vertex

x i+ 1 2 ,j+ 1 2 = (x i+ 1 2 , y j+ 1 2 ), i.e., ∂ t w + ∂ x f (w) + ∂ y g(w) = 0, w 0 (x, y) = w h (x, y, t n ), (x, y) ∈ D i,j , (2.31)
interpreting D i,j = (x i , x i+1 ) × (y j , y j+1 ) as a staggered cell centered in x i+ 1 2 ,j+ 1 2 , (see Figure 2(b)). Before proceeding, some remarks are in order. As previously mentioned, the solution of the initial value problem (2.31) is assumed to have a rather complex structure and one of the reasons for this stems from wave interactions taking place. A proper study requires considering the numerous combinations of initial data that are possible for this two-dimensional problem, and even under the premise that each jump between neighboring initial states projects one planar wave consisting of a single shock, rarefaction or contact discontinuity, the number of possible combinations for a polytropic gas reduces to nineteen [START_REF] Chang | On the 2-D Riemann problem for the compressible Euler equations. I. Interaction of shocks and rarefaction waves[END_REF][START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF] (ignoring the sign of the slip lines, to fifteen [START_REF] Schulz-Rinne | Numerical solution of the Riemann problem for two-dimensional gas dynamics[END_REF]), and for each combination, the Inria solution's complexity is evident in the conjectures of Zhang and Zheng [START_REF] Zhang | Conjecture on the structure of solutions of the Riemann problem for two-dimensional gas dynamics systems[END_REF] and in the numerical experiments [START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF][START_REF] Schulz-Rinne | Numerical solution of the Riemann problem for two-dimensional gas dynamics[END_REF]. In both predicted and numerical results, we are able to perceive that the two-dimensional Riemann problem gives rise to a region of strong interaction consisting of a complex similarity solution. This interaction region can then be approximated in a way conceptually similar to that of the one-dimensional intermediate state in the context of the HLL method.

y x C i,j C i-1,j C i+1,j C i,j-1 C i,j+1 D i,j x i-1 2 x i+ 1 2 y j-
Thus, considering the definition of the local Riemann problem (2.31), it is apparent that the jump discontinuities at the cell's edges lead to two one-dimensional Riemann problems in the x-direction and two one-dimensional Riemann problems in the y-direction, and their effects on one another at the vertex give rise to the region of strong interaction. Since waves propagate with finite velocities, we can then approximately delimit this interaction region by means of the wave model detailed in the subsequent paragraphs of this section.

Computing the one-dimensional HLL smallest and largest wave speed estimates with (2.24) for each of the previously discussed one-dimensional problems yields eight approximate signal velocities that determine the following vectors:

s sw = (
g * ν = s ν n g sν -s ν s g nν + s ν s s ν n (w nν -w sν ) s ν n -s ν s , ν = e or w.
(2.35b) Denote now by O = (x 0 , t 0 ) a local origin in the space-time domain (x, y, t). For a 2D Riemann problem centered at this origin, the extent of its strong interaction region on the xyplane at a time t > t 0 lies by construction within a quadrilateral with time-dependent vertices q sw (t) = x 0 + s sw t, q se (t) = x 0 + s se t, q nw (t) = x 0 + s nw t, q ne (t) = x 0 + s ne t, (2.36) having assumed that t 0 = 0, as in Section 2.1.1. For later convenience, at this stage we define the four points at some fixed small time ∆t > 0

Q sw = (q sw (∆t), ∆t), Q se = (q se (∆t), ∆t), Q nw = (q nw (∆t), ∆t), Q ne = (q ne (∆t), ∆t), (2.37)
and specify a rectangular space-time control volume Q = R × (0, ∆t) that contains these points (2.37), as well as the local origin O. We define the rectangle R = (x w , x e ) × (y s , y n ), with x w ≤ x 0 ≤ x e and y s ≤ y 0 ≤ y n , and identify its four corners as

x sw = (x w , y s ), x se = (x e , y s ), x nw = (x w , y n ), x ne = (x e , y n ).
(2.38)

Wendroff's Nine-State Solver

Wendroff formulated in [START_REF] Wendroff | A two-dimensional HLLE Riemann solver and associated Godunov-type difference scheme for gas dynamics[END_REF] a nine-state two-dimensional version of the three-state HLL Riemann solver, using the literal extension of Godunov's formulation (2.10) to two dimensions as its basis. However, given the absence of an exact solution to the initial value problem (2.31), he employed an approximation W(ξ, η; w ne , w nw , w sw , w se ), or W(ξ, η) for short, to obtain the two-dimensional analogue of (2.13), given by

w n+1 i,j = 1 |C i,j | xi x i-1 2 yj y j-1 2 W x -x i-1 2 ∆t , y -y j-1 2 ∆t ; w n i,j , w n i-1,j , w n i-1,j-1 , w n i,j-1 dy dx + 1 |C i,j | xi x i-1 2 y j+ 1 2 yj W x -x i-1 2 ∆t , y -y j+ 1 2 ∆t ; w n i,j+1 , w n i-1,j+1 , w n i-1,j , w n i,j dy dx + 1 |C i,j | x i+ 1 2 xi yj y j-1 2 W x -x i+ 1 2 ∆t , y -y j-1 2 ∆t ; w n i+1,j , w n i,j , w n i,j-1 , w n i+1,j-1 dy dx + 1 |C i,j | x i+ 1 2 xi y j+ 1 2 yj W x -x i+ 1 2 ∆t , y -y j+ 1 2 ∆t ; w n i+1,j+1 , w n i,j+1 , w n i,j , w n i+1,j
dy dx,

(2.39) under the assumption that the approximated solutions at the vertices do not interact with each other during the time interval ∆t provided the condition ∆t max R=Di,j max µ=s,n ν=e,w

|s µ ν |, |s ν µ | ≤ min ∆x 2 , ∆y 2
.
(2.40)

With our attention directed towards Wendroff's derivation of W(ξ, η), we begin by analyzing the structure of a solution to the 2D Riemann problem (2.25, 2.28) at the planar faces of the space-time control volume Q, illustrated in Figure 3(a). The top surface of Q constitutes nine regions resulting from the finite propagation of waves, in accordance with the wave model introduced in the previous section, up to a small time ∆t. For reference, we show this flat surface in Figure 3(b) and acknowledge the following: the central extent corresponds to the strong interaction region; the four corner zones are simply rectangles containing the undisturbed initial data w ne , w nw , w sw , and w se , written in a counterclockwise order starting from the top right quadrant in the xy-plane; and each of the remaining regions represents the total area covered at time ∆t by the intermediate state (2.34) obtained from the application of a three-state HLL solver at the underlying edge. Let us here summarize the technique used by Wendroff to obtain an approximation for the former region. The central idea is to lump together all of the region's complicated structure into a constant state w * * , in agreement with the one-dimensional approach presented in Section 2.1.2, and make use of the integral conservation laws (2.27) over Q to obtain its specific value.

Performing a simple geometric analysis of the top surface described in the previous paragraph, hereafter defined as S = R × {t = ∆t}, we note that the two-dimensional interactions are contained in the quadrilateral with vertices q sw (∆t), q se (∆t), q nw (∆t), and q ne (∆t), which can be located anywhere on S. For this reason, to simplify the developments and reduce coding difficulties, Wendroff adopts the notation (2.18) and proposes to reformulate (2.32) as

Inria y x s s w ∆t s s e ∆t s e s ∆ t s e n ∆ t t (xsw, 0) (xse, 0) (xne, 0) (xsw, ∆t) (xnw, ∆t) (xne, ∆t) O Q nw Q ne Q se Q sw (a) Space-time control volume Q = R × (0,
s sw = (s s- w , s w- s ), s se = (s s+ e , s e- s ), s nw = (s n- w , s w+ n ), s ne = (s n+ e , s e+ n ), (2.41)
and as a result, relaxes every solution to the subsonic case in two-dimensions, which accounts for the introduction of additional numerical dissipation.

We are now in the position to formulate his explicit approximation of the solution for (2.31). At the small fixed time ∆t, let R µν (µ = n, * , s, ν = e, * , w) be the nine regions with corresponding areas a µν , and then define

W x -x 0 ∆t , y -y 0 ∆t = w µν for all (x, y) ∈ R µν , ( 2
a µν w µν = δx e-0 δy n-0 w ne + δx 0-w δy n-0 w nw + δx 0-w δy 0-s w sw + δx e-0 δy 0-s w se -(f e -f w ) -(g n -g s ) , (2.43)
as long as R * * = 0. Each flux on the right-hand side of the previous equation is obtained from a time-surface integral at the control volume's outer face α, see Figure 3(a), e.g.,

f e = ∆t 2 ((2y u -y 0 -s e n ∆t) f ne + (s e n -s e s ) ∆tf * e + (s e s ∆t + y 0 -2y d ) f se ) . (2.44)
We note that the author in [START_REF] Wendroff | A two-dimensional HLLE Riemann solver and associated Godunov-type difference scheme for gas dynamics[END_REF] does not explicitly mention the procedure to define the transverse fluxes f * e , f * w , g n * and g s * , but instead writes that f * e = f (w * e ) and g n * = g(w n * ), assuming analogous expressions for the other fluxes. Moreover, in (2.43), we used the abbreviations δx α-β and δy α-β to indicate the differences

δx α-β = x α -x β and δy α-β = y α -y β , α, β ∈ {n, s, e, w, 0}. (2.45)
Once the intermediate constant state is found from solving equation (2.43), the approximation (2.42) is properly defined in the subdomain R = D i,j . If we repeat the process for each vertex of C i,j , we obtain the four approximations needed to calculate the value w n+1 i,j (2.39). However, despite the valuable wave model and the consistency with the integral form introduced by Wendroff, the resulting scheme is mainly first-order in both space and time and a higher-order version is not straightforward considering its general formulation. For some details regarding the difficulty in proving the positivity and stability of the method, as well as the behavior of entropy, we refer the reader to [START_REF] Wendroff | A two-dimensional HLLE Riemann solver and associated Godunov-type difference scheme for gas dynamics[END_REF].

Balsara's Multidimensional HLL Solver

In 2010, eleven years after Wendroff's contribution to the recently growing collection of multidimensional solvers, Balsara [START_REF] Balsara | Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows[END_REF] formulated a two-dimensional HLL solver that included closedform, approximate expressions for the fluxes, thus providing a relatively simple implementation of the resulting scheme. Two years later, the same author succeeded in constructing a more robust version of his own solver and presented it in the first pages of [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF]. Here, we concisely detail the fundamental ideas behind the latter, in a way that will be useful for future comparisons with the method proposed in this paper.

′ sw , 0) (q ′ se , 0) (q ′ sw , ∆t) (q ′ nw , ∆t) (q ′ ne , ∆t) σ y sw σ y s * σ y se σ x ne σ x * e σ x se (b) Rectangular prism Q ′ = R ′ × (0, ∆t)
Figure 4: Specific choice for the rectangle proposed by Balsara [START_REF] Balsara | Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows[END_REF][START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF], which bounds the strong interaction region arising from the two-dimensional interaction of waves at the origin O.

Once again, we are interested in finding an approximate solution to the 2D Riemann problem (2.25, 2.28), now conveniently centered at the origin O = (x 0 , t 0 ) with x 0 = (0, 0) and t 0 = 0. It is evident that any strategy that attempts to extend the one-dimensional HLL Riemann solver to multidimensions will somehow involve a constant approximation of the two-dimensional interaction region's composite structure, and Balsara's proposal [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF] is not an exception. Under a subsonic condition, he assumes this region to be bounded at time ∆t > 0 by the rectangle R ′ = (s w ∆t, s e ∆t) × (s s ∆t, s n ∆t), s w , s s < 0 and s n , s e > 0, (2.46) with its vertices q ′ ne , q ′ nw , q ′ sw , and q ′ se respectively located in the four known quadrants of the xy-plane. He then chooses the control volume Q ′ to be the rectangular prism formed with R ′ (2.46) as its base to make the forthcoming integral evaluations easier. Figure 4(b) aims to show Inria this three-dimensional element graphically and Figure 4(a) might assist in visualizing how the interaction region is chosen.

Considering the rectangular prism 2 (s e -s w ) (s n -s s ), an algebraic expression for the constant state w * * can be found based on the integral form (2.27), namely

Q ′ = R ′ × (0, ∆t) with |R ′ | = ∆t
|R ′ |w * * = s s s w ∆t 2 w sw -s s s e ∆t 2 w se -s n s w ∆t 2 w nw + s n s e ∆t 2 w ne + σ x sw f sw + σ x nw f nw -σ x ne f ne -σ x se f se + σ y sw g sw + σ y se g se -σ y nw g nw -σ y ne g ne + σ x * w f * w -σ x * e f * e + σ y s * g s * -σ y n * g n * ,
(2.47) denoting by σ x µe and σ x µw , µ = n, * , s, the areas of the trapezoidal or triangular sections generated by the slowest and fastest waves arising from each of the one-dimensional Riemann problems at the corresponding planar faces x = s e ∆t and x = s w ∆t, as can be appreciated in Figure 4(b). In an analogous manner, the areas of the zones formed at the outer surfaces y = s n ∆t and y = s s ∆t are respectively represented by σ y nν and σ y sν , with ν = e, * , w. To solve for the state w * * , we note that all variables in equation (2.47) are known, with the exception of the fluxes appearing in the last line, i.e., the transverse fluxes introduced briefly in the prior section. Momentarily focusing on the one-dimensional Riemann problem above the x-axis with initial states w l = w nw and w r = w ne , we realize that its solution provided by the HLL approximate Riemann solver yields the intermediate constant state w n * (2.35) and normal flux f n * (2.34), but not the transverse flux g n * . A similar scenario holds for each of the other one-dimensional Riemann problems. However, Balsara offers a solution in [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF], which will be carefully summarized in the following paragraph.

Roughly, each transverse flux can be constructed using values extracted from the associated intermediate state and normal flux. Using the notation introduced in the first paragraph of Section 2, where vector elements are designated by superscripts placed in brackets to avoid confusion with exponents, one constructs the transverse fluxes as

f * ν =     w [2]     , g µ * =    
w [3] µ * w [2] µ * w [3] µ * / w [1] µ * f [2] µ * + ((w [3] µ * ) 2 -(w [2] µ * ) 2 )/ w [1] µ * w [3] µ * f [4] µ * / w [2] µ *

    , (2.48)
where ν = e or w and µ = n or s. Now, the only unknown in equation (2.47) is w * * , which can be expanded by substituting particular values for each of the zone areas at the four faces normal to the main directions x and y such that

|R ′ | ∆t 2 w * * = s s s w w sw -s s s e w se -s n s w w nw + s n s e w ne + s s (f se -f sw ) -s n (f ne -f nw ) + s e (g se -g ne ) -s w (g sw -g nw ) + 1 2 [s w s (f sw -f * w ) -s e s (f se -f * e ) -s w n (f nw -f * w ) + s e n (f ne -f * e )] + 1 2 [s s w (g sw -g s * ) -s s e (g se -g s * ) -s n w (g nw -g n * ) + s n e (g ne -g n * )] = d1.
(2.49)

The focus now shifts from obtaining the interaction state to determining the associated fluxes f * * and g * * . Balsara's approach to derive them is based on the arguments used to define the one-dimensional HLL flux along the t-axis by means of any of the two equations (2.17). In two dimensions, equivalent formulations can be found by employing the integral form (2.27) over sub-rectangular prisms obtained by partitioning the principal space-time control volume so that the time axis is positioned at a face. Among other choices, it can for instance be divided along the x-or y-axis to get four possible volumes Q ′ w = (s w ∆t, 0) × (s s ∆t, s n ∆t) × (0, ∆t), Q ′ s = (s w ∆t, s e ∆t) × (s s ∆t, 0) × (0, ∆t), Q ′ e = (0, s e ∆t) × (s s ∆t, s n ∆t) × (0, ∆t) and lastly, Q ′ n = (s w ∆t, s e ∆t) × (0, s n ∆t) × (0, ∆t). Before proceeding, we remark that the assumed subsonic case (s w , s s < 0 and s n , s e > 0) guarantees that the inverted pyramidal structure, a consequence for the evolution of the rectangular interaction region from time 0 to ∆t, contains the vertical time axis. One can therefore determine the two unknowns f * * and g * * by performing space-time integrations over any two of the above-mentioned volumes. Choosing first Q ′ e , Balsara [START_REF] Balsara | Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows[END_REF][START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF] applies the integral conservation law (2.27) over this spacetime volume to obtain the equality (2.23). In order to obtain the numerical y-flux g * * , the integration of the conservation law (2.27) is performed over the volume Q ′ n shown in Figure 5(a), yielding the expression

y x t w * * O (q ′ sw , 0) (q ′ se , 0) (q ′ se , ∆t) (q ′ sw , ∆t) (a) Sub-rectangular prism Q ′ n y x t w * * O (q ′ nw , 0) (q ′ nw , ∆t) (q ′ ne , 0) (q ′ ne , ∆t) (q ′ sw , 0) (b) Sub-rectangular prism Q ′ e
s e (s n -s s ) w * * -1 2 (s n -s s ) f * * = -
s n (s e -s w ) w * * -1 2 (s e -s w ) g * * = -s n s w w nw + s n s e w ne + 1 2 [s w+ s f sw -s e+ s f se + s e φ hll e -s w φ hll w ] + 1 2 [(2s n -s w+ n ) f nw -(2s n -s e+ n ) f ne -(s e+ n -s e+ s )f * e + (s w+ n -s w+ s )f * w ] -1 2 [(-2s w + s n w ) g nw + (2s e -s n e ) g ne + (s n e -s n w ) g n * ] = d 3 , (2.51)
written in terminology analogous to that described above.

Inria

Using Figure 5 as reference, we are able to understand why the wave speeds involving a plus sign, following the notation established in (2.18), are needed in the previous equations. Basically, their introduction allows us to handle supersonic cases provided they appear. For example, consider the situation where the x-directional Riemann problem with initial states w nw and w ne admits speed estimates s n w < 0 and s n e > 0, so that part of the required subsonic condition (2.46) is satisfied. Now, let us suppose that only the flow below the x-axis is supersonic with strictly positive wave speeds s s w , s s e . Note then that the terms s s+ w g sw and (s s+ e -s s+ w )g s * in equation (2.50) are non-zero, as they should be, providing the needed contributions to the appropriate estimation of f * * .

We merely wish to comment that the system of linear equations (2.49, 2.50, 2.51) can be easily expressed in matrix form Cy = d, specifically

  1 ∆t 2 |R ′ | 0 0 s e (s n -s s ) -1 2 (s n -s s ) 0 s n (s e -s w ) 0 -1 2 (s e -s w )     w * * f * * g * *   =   d 1 d 2 d 3   , (2.52)
such that we are readily able to retrieve det C = 1 4∆t 4 |R ′ | 2 , which without any doubt is strictly positive as long as |R ′ | = 0 (recall that ∆t = 0). From standard linear algebra, we thus find the inverse matrix

C -1 = ∆t 2 |R ′ |   1 0 0 2 s e -2 (s e -s w ) 0 2 s n 0 -2 (s n -s s )   , (2.53)
and compute the unique solution y = C -1 d for the subsonic case.

In the event that the underlying flow is supersonic in both x and y directions, Balsara solves for the intermediate state directly from equation (2.49) and explicitly defines f * * and g * * at point (x 0 , ∆t) as the properly upwinded fluxes

F * * = (f * * , g * * ) =          ( f sw , g sw )
if s s ≥ 0 and s w ≥ 0, ( f se , g se ) if s s ≥ 0 and s e ≤ 0,

( f nw , g nw ) if s n ≤ 0 and s w ≥ 0, ( f ne , g ne ) if s n ≤ 0 and s e ≤ 0, (2.54)
as given in [4, p. 7483]. On the same page, we also find the expressions for the fluxes that are meant to be used in the remaining situations where the flow is fully supersonic in one of the two spatial directions, but subsonic in the other. For specific details concerning the appropriate use of F * * at the cells' interfaces, see Section 3.3.2.

Simple Two-Dimensional HLL Riemann Solver

In this section, a simple method is developed to numerically approximate the solution of the two-dimensional Riemann problem (2.25, 2.28). Suitably built as an extension of the well-known HLL formalism to two dimensions, the scheme relies heavily on the proper utilization of Rankine-Hugoniot relations, which hold across the surfaces of discontinuity that emerge from the origin O, to estimate the constant flux F * * . Hence, before embarking on the details, we must understand what these conditions are like in two dimensions.

Rankine-Hugoniot Relation in Two Dimensions

We are interested in the derivation of the Rankine-Hugoniot relation in two dimensions from a general point of view. The system of conservation laws in (2.25) can be rewritten as

∇ • U [θ] = 0, θ = 1, . . . , ϑ, (3.1)
with ∇ the nabla operator in the physical space (x, y, t) and U [θ] = (f [θ] , g [θ] , w [θ] ). We assume that U [θ] is a bounded measurable function and divergence-free in a weak sense over an open region Ω ⊂ R 2 × (0, ∞) such that

Ω ∇ϕ • U [θ] dx dt = 0, (3.2)
for all continuously differentiable test functions ϕ : Ω → R with compact support and for every θ ∈ {1, . . . , ϑ}.

Let us consider a two-dimensional smooth surface S that splits the region Ω into two open parts Ω l and Ω r . Suppose that U [θ] is smooth in each of the parts' interiors and uniformly continuous up to S, and also that it has limits along S from the left and from the right, denoted by U [θ] l and U [θ] r . Thus, based on the Rankine-Hugoniot relation theorem presented in [START_REF] Zheng | Systems of conservation laws: Two-dimensional Riemann problems[END_REF], equation (3.2) along S is equivalent to

n • U [θ] r -U [θ] l = 0, θ = 1, . . . , ϑ, (3.3)
where n = (n x , n y , n t ) is the unit normal of the surface, pointing from Ω l to Ω r . The previous equation in expanded form is

n x (f [θ] l -f [θ] r ) + n y (g [θ] l -g [θ] r ) + n t (w [θ] l -w [θ]
r ) = 0, but is most commonly expressed as [n x f [θ] + n y g [θ] ] = σ [w [θ] ] , θ = 1, . . . , ϑ, (3.4) having defined σ = -n t and the jump

[(•)] = (•) l -(•) r . (3.5)
For self-similar solutions, the discontinuity surface S described by an equation of the form (a, b, c) • (x, y, t) = 0, may be identified by a similarity curve in the (ξ, η) plane having the form Γ(ξ, η) = aξ + bη + c = 0 and its normal can easily be obtained as the gradient ∇Γ(ξ, η) in the physical space (x, y, t), i.e.,

n = ∂Γ ∂ξ ∂ξ ∂x , ∂Γ ∂η ∂η ∂y , ∂Γ ∂ξ ∂ξ ∂t + ∂Γ ∂η ∂η ∂t = 1 t (a, b, c) . (3.6)
Now that all the theoretical notions have been formally introduced, we can proceed with the complete description of the proposed scheme.

Derivation of Intermediate States and Fluxes

Any approach that deliberately aims to extend the approximate Riemann solver of HLL type to two dimensions involves in some way the constant approximation of intermediate states, as is our case. Furthermore, we expect the solution of the two-dimensional Riemann problem (2.25, 2.28) to be self-similar, as was indicated in Section 2.2, and we need to ensure that an implementation of our method invariably satisfies this property. We therefore set the wave model and space-time control volume Q that will be used in our forthcoming developments to the ones suggested by Wendroff in [START_REF] Wendroff | A two-dimensional HLLE Riemann solver and associated Godunov-type difference scheme for gas dynamics[END_REF], accurately portrayed in Figure 3 and carefully detailed in Section 2.2.1 of this paper. Our strong preference for his setup over other possible ones comes from the fact that it encloses not only subsonic but also supersonic flow structures, e.g., those shown in Figure 6, and most importantly, it provides the correct profile for the eventual use of jump conditions, which is readily seen from the fact that the associated discontinuity planes and their normals are self-similar.

Strongly Interacting State

Regardless of the type, a flow characterized by the considered wave model forms specific patterns at the control volume's faces that can be generalized as follows: the bottom flat surface t = t 0 comprises the four rectangles determined by the initial states w ne , w nw , w sw , and w se , each in its corresponding quadrant; the top face t = ∆t consists of four trapezoids that result from applying the one-dimensional HLL [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF] solver at the underlying edges as described in Section 2.2.1, four rectangles containing unaltered initial data, and a quadrilateral defined by the points Q sw , Q se , Q nw and Q ne (with speeds given in equation (2.32) for the general case) which bounds all two-dimensional interactions at time ∆t that originated from O; and each of the lateral faces retains a structure similar to that of Figure 1 created by a two-wave Riemann fan.

Therefore, all of the individual regions found at the rectangular prism's surfaces have simple geometric shapes and a computation of their areas does not pose a challenge. This facilitates the application of the integral conservation law (2.27) over it, i.e., over the discussed control volume Q = R × (0, ∆t), providing a straightforward expression to obtain w * * given by µ=s, * ,n ν=e, * ,w

a µν w µν = x r y u w ne -x l y u w nw + x l y d w sw -x r y d w se -∆t 2 [(2y u -s e n ∆t) f ne -(2y u -s w n ∆t) f nw + (2y d -s w s ∆t) f sw -(2y d -s e s ∆t) f se ] -∆t 2 [(2x r -s n e ∆t) g ne -(2y l -s n w ∆t) g nw + (2y l -s s w ∆t) g sw -(2x r -s s e ∆t) g se ] -∆t 2 2 [s n e (s n e -s n w ) g n * -(s s e -s s w ) g s * + (s e n -s e s ) f * e -(s w n -s w s ) f * w ] , (3.7)
where the origin O was set to (0, 0, 0), taking x 0 = (0, 0) and t 0 = 0, to make later computations simpler. We wish to note that the nine areas a µν (µ = n, * , s, ν = e, * , w) correspond to the previously identified regions located at the top face, which in turn has a total surface area |R| = δx e-w δy n-s since R = (x w , x e ) × (y s , y n ), and it is rather easy to define the main quadrilateral's extent a * * as a remaining value, notably Before turning to the determination of the flux F * * , we would like to point out that the strong interaction region in the xy-plane at time ∆t is most likely delimited by a sonic line consisting of a circle or ellipse (see [START_REF] Zheng | Systems of conservation laws: Two-dimensional Riemann problems[END_REF]) and well contained in the quadrilateral with vertices q sw (∆t), q se (∆t), q nw (∆t), and q ne (∆t) and area a * * . If a larger interaction region is considered, more dissipation is certainly introduced. For instance, in the case that all wave speeds relax to the minimal and maximal ones defined in equation (2.33), as in [START_REF] Balsara | Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows[END_REF], the volume's top and bottom surfaces will be formed entirely of rectangular regions and the lateral faces of triangular parts, and equation (3.9) would reduce to

a * * = |R| -a ne -a n * -a nw -a * w -a sw -a s * -a se -a * e = ∆t
w * * = ∆t 2 2|R ′ | * [ 2s n s e w ne -2s n s w w nw + 2s s s w w sw -2s s s e w se -s n (f ne -f nw ) -s s (f sw -f se ) -s e (g ne -g se ) -s w (g sw -g nw ) -(s n -s s ) (f * e -f * w ) -(s e -s w ) (g n * -g s * ) ] ,
(3.10) with R ′ = (s w ∆t, s e ∆t) × (s s ∆t, s n ∆t). By comparing this expression with the one derived by Balsara in [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF], it is apparent that the dissimilarities arise from considering different wave configurations at the control volume's lateral faces. To be precise, using Figure 4(b) as reference, the approach described in Section 2.2.2 does not require each corner of the two-dimensional interaction region to exactly coincide with the sides of the two neighboring one-dimensional intermediate scopes, whereas the proposed technique essentially does to allow for the use of jump conditions, as will be explained in the immediate section. We remark that in the limit s β α → s α for α, β ∈ {n, s, e, w}, both equations (3.10) and (2.49) are nevertheless the same.

Inria

Fluxes from Jump Conditions

We will now restrict ourselves to the derivation of closed form expressions for both fluxes f * * and g * * . Recall that for the one-dimensional case, we exposed in Section 2.1.2 the equivalence between the integral relation (2.12), also known as the consistency condition, and the use of Rankine-Hugoniot conditions across the outermost waves, as a means to recover the flux in the intermediate region. We wish to extend this latter concept to two dimensions. For the remaining directions, a similar procedure is employed to get the conditions across the corresponding planes, and a summary of all, including (3.12), is presented below: Relations (3.13) form a system of linear equations and, since the specific value of the strongly interacting state w * * is completely determined by equation (3.9), we opt to rewrite it as

y x t (xsw, 0) (xse, 0) (xne, 0) (xsw, ∆t) (xnw, ∆t) (xne, ∆t) ne O Q nw Q ne Q se Q sw
N δ w 1 (s n w s e n -s w n s n e )(w n * -w * * ) + δ f 1 (s w n -s e n )(f n * -f * * ) + δ g 1 (s n e -s n w )(g n * -g * * )= 0, (3.13a) W δ w 2 (s s w s w n -s w s s n w )(w * w -w * * ) + δ f 2 (s w s -s w n )(f * w -f * * ) + δ g 2 (s n w -s s w )(g * w -g * * )= 0, ( 3
δ f 1 f * * + δ g 1 g * * = δ w 1 (w n * -w * * ) + δ f 1 f n * + δ g 1 g n * = b 1 , (3.14a) δ f 2 f * * + δ g 2 g * * = δ w 2 (w * w -w * * ) + δ f 2 f * w + δ g 2 g * w = b 2 , (3.14b) δ f 3 f * * + δ g 3 g * * = δ w 3 (w s * -w * * ) + δ f 3 f s * + δ g 3 g s * = b 3 , (3.14c) δ f 4 f * * + δ g 4 g * * = δ w 4 (w * e -w * * ) + δ f 4 f * e + δ g 4 g * e = b 4 , (3.14d)
with the unknown fluxes on the left hand-side. It is evident that (3.14) is overdetermined, seeing that there are four equations to solve for two unknowns, and the method of ordinary least squares can be utilized to find a solution. Hence, we express the linear system (3.14) in the form Ax = b, by defining having substituted the quantities defined in (3.13). Equations (3.19) aid to confirm that our proposed approach is able to pick out the right ingredients for the determination of the numerical flux F * * . It is worth observing that for this particular case, f * * (respectively, g * * ) is simply the average of the jump conditions across the eastern and western (respectively, northern and southern) planes of the inverted rectangular pyramid, as expected.

A =     δ f 1 δ g 1 δ f 2 δ g 2 δ f 3 δ g 3 δ f 4 δ g 4     , x = f * * g * * , and b =     δ w 1 (w n * -w * * ) + δ f 1 f n * + δ g 1 g n * δ w 2 (w * w -w * * ) + δ f 2 f * w + δ g 2 g * w δ w 3 (w s * -w * * ) + δ f 3 f s * + δ g 3 g s * δ w 4 (
The above analysis inspired us to develop two alternative formulations to (3.17), which will be duly justified in the subsequent section. For all details pertaining the appropriate use of the resolved fluxes at the primary cells' interfaces, refer to Section 3.3.

Alternative Formulations

As the linear system (3.14) is mathematically overdetermined, we could theoretically propose infinitely many formulations to estimate F * * , not all of which would be sensible. However, in this spirit, we detail two of which will give reasonable solutions and yield shorter expressions than in (3.17), for later interpretation and implementation. The first method gives fluxes that are dependent on the intermediate state w * * , as opposed to the ones provided by the second.

Form I

We first calculate the difference between equations (3.14a) and (3.14c), and separately, the one between (3.14d) and (3.14b), to recover a system of two, not four, linear equations that can be written in condensed form as

δ f 1 -δ f 3 δ g 1 -δ g 3 δ f 4 -δ f 2 δ g 4 -δ g 2 f * * g * * = b 1 -b 3 b 4 -b 2 , ( 3.20)
where the 2×2 real matrix on the left-hand side is denoted by A I . A straightforward substitution of the terms introduced in (3.13) into this matrix allows us to compute its determinant as det A I = -4 ∆t 2 a * * , which is certainly less than zero on the assumption that both ∆t and a * * are positive quantities. The unique solution of (3.20)

Consistency

In the next few paragraphs, we give various statements concerning the consistency of the proposed numerical scheme, where w * * is defined by equation (3.9) and F * * by (3.17). For this, let us define a state w constant in x ∈ R 2 , as well as we , wn , ww , and ws constant in x > 0, y > 0, x < 0 and y < 0, respectively, with the corresponding fluxes denoted in a similar way.

Given that the scheme is in conservative form, we need to verify if the numerical fluxes are consistent with the physical ones, i.e., if f * * ( w, w, w, w) = f ( w) and g * * ( w, w, w, w) = g( w). Making use of equations (3.19) with the speeds defined as s e = s n = s = -s s = -s w , and recalling the fact that the one-dimensional HLL fluxes are consistent, e.g., f * e ( w, w) = f ( w), we surely recover that these basic consistency equalities are satisfied.

In addition, if all variations occur in one spatial direction, as depicted in Figure 8, equation (3.9) reduces to an analogue of (2.16). For instance, if the variation is restricted to the y-direction, we are sure that wn = w ne = w nw , ws = w se = w sw , F n = F ne = F nw and F s = F se = F sw , and the equality (3.9) changes into

w * * = s n wn -s s ws -(ḡ n -ḡs ) s n -s s . ( 3

Extensions and Computational Remarks

The goal of this section is to exhibit a general representation of the proposed solver for twodimensional meshes satisfying certain properties and simultaneously provide insight into its implementation. We then elaborate on how to assemble the total flux across the mesh edges or interfaces and on how to obtain a second-order version of the scheme.

Non-Rectangular Meshes

In the framework of finite volume methods, the technique presented in Section 3 naturally relies on integral relations that can be easily applied over any right prism in the physical space (x, y, t). So, let T be an admissible mesh defined over an open bounded region Ω ∈ R 2 in the sense of [START_REF] Eymard | Finite volume methods[END_REF] (Definition 9.1) consisting of polygonal cells C c ∈ T (with c a unique index), edges, and vertices; the latter belonging to a family denoted by P. Moreover, for each vertex p ∈ P, we construct a dual convex cell D p by connecting the centers of the polygons that share this point and establish the space-time control volume Q p = D p × (0, ∆t). An integration of (2.25) over this prism Q p then yields the expression

Dp w(x, ∆t) dx - Dp w(x, 0) dx + Qp ∇ • F (w(x, t)) dx dt = 0, (3.28)
with F = (f , g) and x = (x, y), which is known to be equivalent to

Dp w(x, ∆t) dx = Dp w(x, 0) dx - Sp F (w(s, t)) • dS, (3.29)
by application of the divergence theorem, given S p = ∂D p × (0, ∆t). For presentation purposes, we shall hereafter limit ourselves to the case where D p is a basic triangle. However, we urge the reader to keep in mind that all subsequent developments can be readily generalized to any convex polygonal dual cell. Thus, we consider the volume Q p to be a triangular prism, provided D p is adequately defined by joining together the three distinct centers

x α = (x α , y α ), α = 1, 2, 3, (3.30) Dp C 2 C 3 C 1 y x (a) Global view of possible structures y x w 3 w 1 w * 1 w * 2 w 2 w * * w * 3 p x 3 x 2 x 1 L l 1 L 1 L r 1 s 1 r ∆ t s 1 l ∆ t s 2 l ∆ t s 2 r ∆ t s 3 l ∆t s 3 r ∆t ν 1
(b) Zoom on a triangular dual cell of the corresponding primary cells C α , as exemplified in Figure 9. Owing to its construction, the dual cell contains the vertex p that is evidently the intersection point of three edges L α , for which we know their respective formulations and unit normals ν α = (ν αx , ν αy ). We employ the notation L α to designate the line segment having p as one of its endpoints and lying between x α and x α+1 , with the index α following a cyclic order such that for α = 3, α + 1 = 1.

As explained in Section 2.2, a constant state w α is assumed within each cell C α at the initial time t 0 ≥ 0, for α = 1, 2, 3. Using now Figure 9(b) as reference, it is evident that in order to obtain the desired values at the face L α × (t 0 , t 0 + ∆t) in the interior of Q p , we must solve not only a one-dimensional Riemann problem in its normal direction, but also a local two-dimensional Riemann problem, involving initial data w 1 , w 2 and w 3 , at the vertex p.

For the former, we identify as s α l and s α r the left and right HLL wave speeds in the direction of ν α , after making use of equation (2.24) with w l = w α+1 and w r = w α . Moreover, solving the one-dimensional Riemann problems at the edges allows us to derive expressions for the constant state w * α (2.16) and HLL fluxes (2.23) associated with each problem, i.e.,

w * α = s α r w α -s α l w α+1 + ν α • (F α+1 -F α ) s α r -s α l , (3.31a) f n * α = ν α • (s α r F α+1 -s α l F α ) + s α l s α r (w α -w α+1 ) s α r -s α l , (3.31b)
where f n * α = ν α • F * α is the flux perpendicular to the edge α ∈ {1, 2, 3}. We recall that the transverse fluxes, denoted here as

f n⊥ * α = (-ν αy , ν αx )•(f * α , g * α )
, can be constructed using values extracted from the associated normal flux (3.31b) and intermediate state (3.31a), as was done in Section 2.2.2. In addition, we easily recover the two parallel lines that bound the extent of w * α on the xy-plane at some time t > t 0 , which are given by

L r α (t) = {(x, y) | ν αx x + ν αy y = s α r (t -t 0 )} , (3.32a) L l α (t) = {(x, y) | ν αx x + ν αy y = s α l (t -t 0 )} , ( 3

.32b)

Inria having conveniently set p = (0, 0) and recalling that ν α = 1. From now on, we also consider the time t = ∆t and t 0 = 0. With regard to the two-dimensional Riemann problem found at the origin O = (p, 0), we specify a triangular region by connecting the vertices q α = (x ′ α ∆t, y ′ α ∆t), for α = 1, 2, 3, which are essentially the intersection points of the bounding lines L l α (∆t) and L r α+1 (∆t) such that

x ′ α = s α+1 r ν αy -s α l ν (α+1)y ν (α+1)x ν αy -ν αx ν (α+1)y and y ′ α = s α l ν (α+1)x -s α+1 r ν αx ν (α+1)x ν αy -ν αx ν (α+1)y . (3.33)
Let us note that all two-dimensional complex interactions emanating from O and projected onto the planar surface t = ∆t, are well-contained in the previously mentioned triangle, which in turn generates an inverted triangular pyramid during its time evolution from 0 to ∆t. The previous statements facilitate the derivation of an algebraic equation to compute the constant state w * * , when applying the integral conservation law (3.29) over the control volume Q p . However, for the determination of the flux f * * , we still need to provide the Rankine-Hugoniot relations that hold across the surfaces of discontinuities emerging from O, following the approach carefully detailed in Section 3.2.2. We therefore start by identifying the points Q α = (q α , ∆t) and the normals

n α = --→ OQ α+1 × --→ OQ α = ∆t 2 (y ′ α+1 -y ′ α ) i + (x ′ α -x ′ α+1 ) j + (x ′ α+1 y ′ α -x ′ α y ′ α+1 ) t , (3.34)
that will allow us to deduce the jump conditions for α = 1, 2, 3 as

δ w α (x ′ α+1 y ′ α -x ′ α y ′ α+1 )(w * α -w * * ) + δ f α (y ′ α+1 -y ′ α )(f * α -f * * ) + δ g α (x ′ α -x ′ α+1 )(g * α -g * * ) = 0, (3.35)
which can be written in the following form:

δ f α f * * + δ g α g * * = δ w α (w * α -w * * ) + δ f α f * α + δ g α g * α = b α , α = 1, 2, 3. (3.36)
Clearly, relations (3.36) form once again an overdetermined system of linear equations, since there are more equations (three) than unknowns (two). In view of the method of ordinary least squares, we then express this system as A T x = b T by defining

A T =   δ f 1 δ g 1 δ f 2 δ g 2 δ f 3 δ g 3   , x = f * * g * * , and b T =   δ w 1 (w * 1 -w * * ) + δ f 1 f * 1 + δ g 1 g * 1 δ w 2 (w * 2 -w * * ) + δ f 2 f * 2 + δ g 2 g * 2 δ w 3 (w * 3 -w * * ) + δ f 3 f * 3 + δ g 3 g * 3   , (3.37)
and look for the least squares solution

x = M T -1 A T T b T , provided M T is invertible.
It is worth mentioning that in the general case when three or more edges (ε ≥ 3) meet at the vertex p, the matrix

M T =        α=1,...,ε δ f α 2 α=1,...,ε δ f α δ g α α=1,...,ε δ f α δ g α α=1,...,ε δ g α 2        , (3.38)
has a determinant that can be expressed as the sum of square binomials of the form (δ f α δ g β -δ f β δ g α ) 2 , for all α, β ∈ {1, 2, . . . , ε}, which is geometrically proven to be strictly positive as long as the area of the interaction region a * * is not zero.

Hence, we confirm that the proposed approach is simple enough to be applied on any admissible mesh, yielding useful closed expressions for the intermediate state w * * and flux F * * .

Flux Assembling at Faces

The assembling of the total flux at the cells' faces is an important aspect that has to be carefully considered. Although we perform the subsequent study focusing on the rectangular mesh used in Section 3, we note that all formalisms developed here can be easily adapted to other configurations such as the ones presented in the previous section.

After analyzing the structures shown in Figures 2(a) and 9(a), we can infer that the final expression will be a convex combination of one-and two-dimensional fluxes, as in [START_REF] Balsara | Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows[END_REF][START_REF] Brio | Two-dimensional Riemann solver for Euler equations of gas dynamics[END_REF][START_REF] Kurganov | Semi-discrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations[END_REF]. Each flux at an edge of the primary cell C i,j can be obtained from the surface average of the flux normal to the underlying space-time face, e.g., φ x,i+ 1 2 ,j = 1 ∆y ∆t

y j+ 1 2 y j-1 2 ∆t 0 F (w(x i+ 1 2 , t)) • ni+ 1 2 dt dy, (3.39)
being ni+ 1 2 = (1, 0) the unit normal vector of S i+ 1 2 = (y j-1 2 , y j+ 1 2 ) × (0, ∆t). In the subsonic case, the above definite integral can be evaluated by considering the contributions at the face coming from the one-and two-dimensional Riemann solvers initially applied at the cell's edge and vertices, respectively. The averaged flux becomes

φ x,i+ 1 2 ,j = θ s,i+ 1 2 ,j+ 1 2 φ hll2D x,i+ 1 2 ,j+ 1 2 + θ y,i+ 1 2 ,j φ hll x,i+ 1 2 ,j + θ n,i+ 1 2 ,j-1 2 φ hll2D x,i+ 1 2 ,j-1 2 , ( 3.40)
with θ s,i+ 1 2 ,j+ 1 2 = ∆t 2 ∆y |ŝ s,i+ 1 2 ,j+ 1 2 |, θ n,i+ 1 2 ,j-1 2 = ∆t 2 ∆y |ŝ n,i+ 1 2 ,j-1 2 |, and θ y,i+ 1 2 ,j = 1 -θ s,i+ 1 2 ,j+ 1 2 -θ n,i+ 1 2 ,j-1 2 .
It is obvious that the weights θ s,i+ 1 2 ,j+ 1 2 and θ n,i+ 1 2 ,j-1 2 determine the amount of two-dimensional contributions that the total flux at the cell's face will possess, and they are directly linked through the time step with the CFL number used in practice (for details regarding the associated CFL condition, we refer the reader to [3, p. 1977] where it is fully explained). We remark that a simpler alternative to (3.40) is fixing the weights to the coefficients in Simpson's rule, used to numerically integrate the average of F (w(x i+ 1 2 , ∆t)) • ni+ 1 2 over the edge (y j-1 2 , y j+ 1 2 ), yielding

φ x,i+ 1 2 ,j = 1 6 φ hll2D x,i+ 1 2 ,j+ 1 2 + 4 6 φ hll x,i+ 1 2 ,j + 1 6 φ hll2D x,i+ 1 2 ,j-1 2 . ( 3.41)
Let us note that in equations (3.40) and (3.41), we have employed new variables that need to be appropriately defined. First, the element having the superscript "hll" is essentially the one-dimensional numerical flux described in equation (2.21), but with a two-dimensional index instead, such that φ hll x,i+ 1 2 ,j = φ hll (w n i,j , w n i+1,j ). In a similar manner, we specify

φ hll2D x,i+ 1 2 ,j+ 1 2 = φ hll2D x (w n i+1,j+1 , w n i,j+1 , w n i,j , w n i+1,j ), (3.42a
)

φ hll2D x,i+ 1 2 ,j-1 2 = φ hll2D x (w n i+1,j , w n i,j , w n i,j-1 , w n i+1,j-1 ), (3.42b)
where φ hll2D x,m,n is the local flux f * * of the two-dimensional Riemann problem defined at the vertex x m,n , with analogous considerations for the y direction. Next, the wave speeds denoted by ŝα , for α ∈ {n, s, e, w}, are determined from the intersection of the interaction region with the xand y-axes at time ∆t (notice the circles pictured in Figure 3

Inria and in the limit s β α → s α for α, β ∈ {n, s, e, w}, ŝα = s α . Equations (3.40) and (3.41) are strictly valid only for a flow that is subsonic everywhere. However, if we reformulate our signal speeds defined in (2.32) as s n+ e , s e+ n , s n- w , s w+ n , s s- w , s w- s s s+ e , and s e- s , and utilize them to recover all one-and two-dimensional states and fluxes associated with our multidimensional Riemann solver (2.34, 2.35, 2.48, 3.9, 3.17), we are able to account for supersonic situations and employ both (3.40) and (3.41) without any additional adjustments. Another way to deal with supersonic situations is to approximately define the fluxes φ hll2D

x,m,n and φ hll2D y,m,n as the upwinded fluxes at the space-time point (x m,n , ∆t), as given in [4, p. 7483] and partially described in equation (2.54). Henceforth, we use the term Simpson assembling to denominate equation (3.41) together with these upwinded 2D fluxes.

The most accurate approach to handle supersonic flows is to manually track the position of the interaction region and deduce the actual elements that contribute to the total flux at each of the cell's faces. For instance, on examining the configuration shown in Figure 6(a), it is clear that for this example, f * * does not have an effect in the total flux assembled at the vertical edge, and choosing f * w when assembling would certainly be more accurate. For this, we choose to introduce a new expression

φ x,i+ 1 2 ,j = θ s,i+ 1 2 ,j+ 1 2 φ hll2D x,i+ 1 2 ,j+ 1 2 + θ y,i+ 1 2 ,j φ hll x,i+ 1 2 ,j + θ n,i+ 1 2 ,j-1 2 φ hll2D x,i+ 1 2 ,j-1 2 , ( 3.44)
where

θ s,i+ 1 2 ,j+ 1 2 = ∆t 2 ∆y |s s,i+ 1 2 ,j+ 1 2 |, θ n,i+ 1 2 ,j-1 2 = ∆t 2 ∆y |s n,i+ 1 2 ,j-1 2 |, and θ y,i+ 1 2 ,j = 1 -θ s,i+ 1 2 ,j+ 1 2 -θ n,i+ 1 2 ,j-1 2 .
Here, the x-directional flux φ hll2D x,m,n needs to be regarded as a convex combination of local fluxes associated to the nine states, characterized by the set L m,n = { * * , n * , s * , * e, * w, ne, nw, sw, se}, of the two-dimensional Riemann solver centered at x m,n , i.e., φ hll2D x,i+ 1 2 ,j+

1 2 = φ hll2D x (w n i+1,j+1 , w n i,j+1 , w n i,j , w n i+1,j ) = (µν) ∈ L i+ 1 2 ,j+ 1 2 β y,µν f µν , (3.45a) φ hll2D x,i+ 1 2 ,j-1 2 = φ hll2D x (w n i+1,j , w n i,j , w n i,j-1 , w n i+1,j-1 ) = (µν) ∈ L i+ 1 2 ,j-1 2 β y,µν f µν , (3.45b)
and β y,µν ≥ 0, β y,µν = 1, with similar equations and notation for the y-direction. Each coefficient β y,µν corresponds to the fraction of the interaction surface associated to the local state µν. In a subsonic configuration, notice that only β y, * * = 1 and all others are zero. Equation (3.44) will now be referred to as manual assembling and in Annex B, we provide the specific details for the implementation of (3.44, 3.45) with the help of pseudocode. There, we observe that the third algorithm corresponds to the approximate (due to some simplifications) two-dimensional analogue of having written equation (2.20) as

φ hll (w l , w r ) =      f l if 0 < s l , f * if s l < 0 < s r , f r if 0 > s r .
(3.46)

Predictor-Corrector Scheme of Second-Order Accuracy

The first order scheme is now complete and the next step is to define an appropriate second-order extension. This is achieved by means of a simple predictor-corrector approach similar to the one proposed by Balsara in [START_REF] Balsara | Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows[END_REF][START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF], which roughly consists of using the two-dimensional Riemann solver for both steps involved in the algorithm. A consequence of employing this solver in the corrector step is the introduction of more isotropy into the simulation. Instead, "our use of the multidimensional Riemann solver in the predictor step has the happy consequence of raising the maximal CFL number", in the words of the author of [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF], and we add the phrase "when needed".

In the predictor step, the vector w is spatially reconstructed from the center of the primary cell to its corners or vertices following a MUSCL-type approach, i.e., using piecewise linear interpolations with slope limiters such as the MC limiter [START_REF] Van Leer | Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection[END_REF], minmod [START_REF] Roe | Characteristic-based schemes for the Euler equations[END_REF] or the positive preserving limiter [START_REF] Suresh | Positivity-preserving schemes in multidimensions[END_REF]. The reconstruction provides the four states that are necessary to solve the multidimensional Riemann solver (2.34, 2.35, 2.48, 3.9, 3.17) at each corner, yielding the xand y-directional fluxes that will contribute to the assembling at the cell's faces. Now, to avoid the appearance of spurious solutions for certain second-order simulations, we propose employing

φ x,i+ 1 2 ,j = θ y,i+ 1 2 ,j+ 1 2 φ hll2D xs,i+ 1 2 ,j+ 1 2 + θ y,i+ 1 2 ,j φ hll x,i+ 1 2 ,j + θ y,i+ 1 2 ,j-1 2 φ hll2D xn,i+ 1 2 ,j-1 2 , ( 3.47)
where

θ y,i+ 1 2 ,j± 1 2 = ζ∆t 2 ∆y max(|s s,i+ 1 2 ,j± 1 2 |, |s n,i+ 1 2 ,j± 1 2 |) and θ y,i+ 1 2 ,j = 1 -θ y,i+ 1 2 ,j+ 1 2 -θ y,i+ 1 2 ,j-1
2 , with analogous expressions for the other main direction. We note that we have introduced new variables "hll2D" which are defined in a way similar to (3.45) but with a different subscript that denotes the location relative to the local vertices (we refer the reader to Algorithm 4 in Annex B for more details regarding their precise definition). Moreover, the flux φ hll x,i+ 1 2 ,j at the center of the face is simply the average of the nearest two x-directional HLL fluxes, one initially computed below the vertex x i+ 1 2 ,j+ 1 2 and the other above x i+ 1 2 ,j-1 2 , as partial inputs for our two-dimensional Riemann solver. We then utilize the total assembled fluxes to estimate the vector of conservative variables at the half time step t n + ζ ∆t, with ζ = 1 2 . For the corrector step, we repeat the operations carried out in the predictor step but having set ζ = 1. The slopes computed at time t n are now applied to spatially reconstruct variables that are centered both in time and space, i.e., the output of the predictor step.

In the coming section, we will present several numerical tests that were obtained with this simple yet second-order accurate predictor-corrector approach. First-order approximations can be obtained using the same computer code, by performing only one step with unreconstructed states and ζ = 1.

Numerical Results

The purpose of this section is to validate our scheme with several multidimensional test problems on a uniform rectangular mesh. The numerical implementation of our solver and Balsara's [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF] (hereafter referred to as BAL2012) has been done in the HERACLES code [START_REF] González | HERACLES: A three-dimensional radiation hydrodynamics code[END_REF] for astrophysical fluid flows. By having a common computational framework, we can fairly compare the accuracy and robustness of both methods.

We note that employing Simpson assembling at the cells' faces for our solver yields almost the exact same results as BAL2012, for which this type of assembling is the default for second-order approximations, and thus, unless stated otherwise, we opt to display our method's solutions with the manual assembling (3.47). All tests were run with a CFL number of 0.9 and making use of the predictor-corrector scheme mentioned in Section 3.3.3. As for the choice of slope limiters, we applied the MC limiter [START_REF] Van Leer | Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection[END_REF] for all except the last (Section 4.4), where minmod [START_REF] Roe | Characteristic-based schemes for the Euler equations[END_REF] was utilized instead.

Inria

Accuracy Analysis

We wish to estimate the rate at which the L 1 error for the proposed scheme decreases as the numerical grid is refined. For this, we consider the initial density profile [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF] ρ 0 (x, y) = 1 + 0.2 sin(π(x + y)), (4.1) together with the velocities and pressure defined in Section 4.1 of [START_REF] Liska | Comparison of several difference schemes on 1D and 2D test problems for the Euler equations[END_REF], i.e., u 0 = 1, v 0 = -0.5 and p 0 = 1. The simulation is run to time t = 4, which corresponds to the time it takes for the wave to be advected once around the periodic domain spanning [0, 2] × [0, 2]. The final state is then compared with the analytical one. Table 1: L 1 density errors and orders of accuracy for the wave advection test, using our proposed scheme with different slope limiters.

In Table 1, the accuracy results for our scheme are summarized. For both assembling methods being compared, the L 1 density errors decrease as the numerical resolution increases and we are able to see that the lowest values are obtained when the MC limiter and manual assembling are present. In addition, we observe that second-order accuracy is reached when any of the three selected limiters are used, as was previously mentioned in Section 3.3.3. We mention that our method with Simpson assembling not only gives roughly the same results as BAL2012 (e.g., differences after the seventh or eight decimal place for the MC limiter results) but also takes the same amount of CPU time to complete with HERACLES. The advantage, in terms of accuracy, of taking into account the multidimensional contributions for this problem is evident when we compare the rows corresponding to the MC limiter in Table 1 with those of Table 2. The latter were obtained utilizing only 1D HLL Riemann solvers in both steps of the predictor-corrector algorithm. Firstly, we consider the two-dimensional problem described in [5, p. 183], with initial data given in Table 3 (left). This Riemann problem (RP) initially consists of two contact discontinuities J and two forward shock waves -→ S , specifically -→ S 21 J 32 J 34 -→ S 41 borrowing the notation used in [START_REF] Schulz-Rinne | Numerical solution of the Riemann problem for two-dimensional gas dynamics[END_REF]. In general, we expect that both slip lines encounter the sonic circle of the constant state in the third quadrant of the xy-plane and bend to end in spirals inside the subsonic area of the circle's portion lying in this quadrant. In addition, from the interaction of the shocks -→ S 21 and -→ S 41 , we await the appearance of a pair of three-shock configurations, such that part of the subsonic area is bounded by two joining Mach shocks and two reflected shocks. Table 3: Initial data for the first MultiD RP described in [START_REF] Brio | Two-dimensional Riemann solver for Euler equations of gas dynamics[END_REF] and the sixth of [START_REF] Kurganov | Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers[END_REF].

Multidimensional Riemann Problems

Multidimensional Riemann Problem 3

Quadrant ρ 0 (x, y) u 0 (x, y) v 0 (x, y) p 0 (x, y)

x > 0, y > 0 1.5 0.0 0.0 1.5

x < 0, y > 0 0.5323 [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF].

By means of the contour plots shown in Figures 10 and11, we can analyze the solutions of this MultiD Riemann problem, computed on uniform grids of 400 2 and 1000 2 cells, for both methods being compared. All results follow the expected behavior described in the above paragraph, with the additional property of being symmetric with respect to the x = y line, as was also anticipated. There is no significant difference between the plots obtained with the manual assembling and those of BAL2012 for this particular test. Let us note that the low resolution of the contact discontinuities is not surprising given that both methods are based on the two-wave model of the HLL Riemann solver. Next, we study the multidimensional Riemann problem ---------→ J 21 J 32 J 34 J 41 that involves only slip line initial data, summarized on the right of Table 3. Particularly, we expect the solutions to have a vortex-type structure that turns in clockwise direction, with contact discontinuities spiraling around the center, and this is the case for the numerical results presented in Figure 12. We can observe that the ripples that are created in the first and third quadrants of the left plot have comparable resolution to those found in [START_REF] Kurganov | Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers[END_REF][START_REF] Lax | Solution of two-dimensional Riemann problems of gas dynamics by positive schemes[END_REF][START_REF] Schulz-Rinne | Numerical solution of the Riemann problem for two-dimensional gas dynamics[END_REF]. The detail of the ripples and slip lines is greatly improved when the number of zones is increased, e.g., to one million cells, as shown in the right image. We wish to mention that once again there is no visible difference between the contours obtained with our scheme and the ones acquired with BAL2012, and this is the reason why we only displayed the former. 4. The expected behavior of this problem is properly specified in Configuration 3 of [START_REF] Schulz-Rinne | Numerical solution of the Riemann problem for two-dimensional gas dynamics[END_REF]. Here, we briefly mention that during its time evolution, the solution develops a double Mach reflection and a shock propagates in the southwest direction Inria at a 45-degree angle to the grid lines. This can be appreciated in Figures 13 and14, where we display the density variable ρ at time t = 1.1 by means of contour plots obtained with our scheme and BAL2012, on two different meshes. Clearly, the region of strong and complex interactions associated with the problem is located in the third quadrant, where we are able to observe a well-resolved mushroom cap (especially in the finer grid where it is sharp), which is consistent with expectations. Once more, all solutions are visually symmetric with respect to x = y. Figure 14: Density estimated with BAL2012 (left) and our scheme (right) for MultiD RP3, using 1000 × 1000 cells on 64 proc.; contours chosen as in Fig. 13 (ρ:0.14-1.75).

We wish to note that if we do not properly assemble the flux at each cell's faces for this problem, second-order computations will suffer from the spurious solution known as the carbuncle phenomenon. In fact, prevention of this instability is what inspired us to derive and recommend equation (3.47) instead of (3.40) for the assembling (Algorithms 3 and 2 in Annex B, respectively). For first-order approximations, both mentioned expressions are equivalent; the results shown in Figure 15 were created using extremely fine meshes of 16 million and 100 million cells, proving the robustness of the first-order scheme.

Double Mach Reflection

The double Mach reflection problem proposed by Woodward and Colella [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF] starts off as a Mach 10 oblique shock in air encountering a reflecting wall. Using the set-up originally given in [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF], we run the simulation until the final time t = 0.2 for the range of resolutions considered by Balsara [START_REF] Balsara | Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows[END_REF][START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF] Figure 16 shows twenty-five density contours obtained with the two methods being compared, on the above-mentioned meshes. We are able to see that the jet formed by the double Mach reflection is well captured, especially on the finest grid; in all plots, we can observe the slipping contact line that leads around to the forward moving Mach stem, which rolls-up creating a vortex head. For the single-step Eulerian MUSCL results found in Figure 9e of [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF], Woodward and Colella provided an explanation for the oscillations and noise present near the slowly moving shock, which we now quote as it directly applies to our case: "the shocks are extremely thin, but this thinness has permitted a numerical instability to generate noise where the shocks move slowly and are nearly aligned with the mesh". However, despite this noise, all results are satisfactory.

Sedov Explosion

The multidimensional blast test presented by Sedov in 1946 [START_REF] Sedov | Propagation of strong blast waves[END_REF] comprises an intense explosion resulting from a punctual quantity of energy placed in the center of the numerical domain.

Inria

We then expect the solution to be a strong spherical shock propagating outwards towards the boundaries, which are set to be periodic. For this problem, both the ambient gas density and the explosion energy are initially set to unity. The latter is deposited in the central cell of a 65 × 65 or 129 × 129 grid covering the computational domain [-0.5, 0.5] × [-0.5, 0.5] and the simulations are run until a time t = 0.2, i.e., before the shock reaches the boundaries. The grid is purposely chosen coarse to be able to easily detect the anisotropic behavior commonly observed when performing this test with traditional Godunov codes. The results shown in Figures 17 and18 for BAL2012 and our scheme were obtained using Simpson assembling, in order to reiterate the fact that both methods yield almost identical numerical solutions when this type of assembling is employed, and they are surely more isotropic than those obtained with the conventional second order HLL scheme. In addition, Figure 18 aids in understanding that the more we refine the grid, the more the corresponding solutions will resemble the analytical ones.

Conclusions

We have described a simple multidimensional Riemann solver for compressible homogeneous flows governed by the Euler system of equations. The associated numerical strategy defines an approximated profile of 2D Riemann problems composed of plane waves and makes use of Rankine-Hugoniot conditions as a guideline to adequately derive constant state approximations on both sides of the discontinuities. The MultiD solver is a two-dimensional extension of the well-known HLL scheme for the four-quadrant Riemann problem that generalizes the 2D solver proposed by Balsara [START_REF] Balsara | Multidimensional HLLE Riemann solver: Application to Euler and magnetohydrodynamic flows[END_REF][START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF]. For the considered approximated profile consisting of nine constant states, jump conditions lead to an overdetermined system that we solved using a least squares approximation. Notwithstanding, the derived numerical 2D fluxes look remarkably similar to the typical HLL flux and all formulations reduce to those of the 1D solver when the initial Riemann data model a one-dimensional flow.

Sample numerical results presented in this paper show the effectiveness and robustness of the proposed methodology when applied to subsonic and supersonic flows. For the latter, particular attention must be paid when assembling the total flux at the cells' faces with varying weights; therefore, we have provided a straightforward and robust assembling approach, comparable to that which uses weights fixed to the coefficients in Simpson's rule for all time steps. In addition to the simplicity, we also propose a generalization to unstructured grids with a formulation that is mostly algebraic rather than geometrical and, following this line, we argue that there is a way to derive an HLL solver for Riemann problems with an arbitrary number of initial constant states connected at a single point.

Jump conditions can be improved by designing complex profiles so that the Rankine-Hugoniot relations define an invertible system with intermediate states and associated fluxes as unknowns. Moreover, given that the consistency with the integral formulation through these relations holds in three dimensions as well, we believe that a genuine three-dimensional solver can readily be obtained in future work.

Inria

After obtaining the eight approximate signal velocities that determine the vectors (2.32), we propose to restrict in some cases the diagonal crossing of the interaction region corners (2.36), in order to reduce further coding difficulties. Basically, we want to avoid the type of situations where a corner lies in the quadrant diagonal to that where its subsonic counterpart would be and, at the same time, at least two of the remaining corners are in their respective "subsonic" quadrants. For this, we use Algorithm 1 Restricting certain crossings We then recover all one-and two-dimensional states and fluxes with equations (2.34), (2.35), (2.48), (3.9), and (3.17). To estimate the values of the speeds and fluxes mentioned in the first paragraph of this annex, we suggest considering the following pieces of pseudocode: Algorithm 2 Defining the speeds sα , for α ∈ {n, s, e, w}

Inria

Algorithm 3 serves to understand how the "hll2D" fluxes found in equation (3.44) are to be defined. However, we introduce a more general algorithm associated with the manual assembling (3.47) that serves to obtain robust first-and second-order approximations. Observing closely, we note that Algorithm 3 is somehow contained in what follows:

Figure 1 :

 1 Figure 1: One-dimensional HLL(E) Riemann problems.

  y, t 2 ) dy dx = (x a , y, t)) dy dtx, y c , t)) dx dtx, y d , t)) dx dt.

Figure 2 :

 2 Figure 2: Example of a solution's structure at time t = ∆t, resulting from a series of one-and two-dimensional HLL Riemann problems on a rectangular mesh.

Figure 3 :

 3 Figure 3: Structures formed by the outward propagation of waves from the staggered cell's origin O and edges, as suggested by Wendroff in [37].

  Surface at time t = ∆t of Q

Figure 5 :

 5 Figure 5: Sub-rectangular prisms used by Balsara [3, 4] to determine the x-flux f * * and y-flux g * * associated with the strong interaction region.

Figure 6 :

 6 Figure 6: Two examples of structures formed at time t = ∆t by outward propagating waves related to flows that are supersonic in at least one of the spatial directions.

Figure 7 :

 7 Figure 7: Planes in the interior of the space-time volume Q. Let us denote by ne = n e / n e the unit normal of the plane S e that contains the three points Q se , Q ne , and O, as illustrated in Figure 7. Without difficulty, we identify the distinct vectors --→ OQ ne = (s n e ∆t, s e n ∆t, ∆t) and --→ OQ se = (s s e ∆t, s e s ∆t, ∆t) that lie on the flat surface and determine their cross product to compute the normal n e , namely n e = --→ OQ ne × --→ OQ se = ∆t 2 [(s e n -s e s ) i -(s n e -s s e ) j + (s n e s e s -s e n s s e ) t] , (3.11) with i, j and t the standard basis in the physical space (x, y, t). The Rankine-Hugoniot condition across the discontinuity surface S e is then n e • (f * *f * e , g * *g * e , w * *w * e ) = 0, obtained by means of equation (3.4), which can be rewritten as (s n e s e s -s e n s s e )(w * ew * * ) + (s e n -s e s )(f * ef * * ) + (s s e -s n e )(g * eg * * ) = 0. (3.12)

(δ g and δ f 1 = δ f 3 = δ g 2 = δ g 4 = 0 .δ g 1 2 = -k 14 , k 21 = 2δ g 1 δ f 2 2=

 401142 w * ew * * ) + δ f 4 f * e + δ g 4 g * e     , (3.15) and write the normal equations in matrix notation as A T Ax = A T b. The least squares solution x = M -1 A T b = Kb, provided M = A T A can be inverted, is the exact one if it exists or an approximate one if it does not. Considering that M has in fact a strictly positive determinant (A.3) and is consequently nonsingular (see Annex A), we are thus able to get explicit forms for the fluxes in the interaction region as f * * = [k 11 b 1 + k 12 b 2 + k 13 b 3 + k 14 b 4 ] / det M , g * * = [k 21 b 1 + k 22 b 2 + k 23 b 3 + k 24 b 4 ] / det M , of the matrix elements k 11 = δ f 1 This further implies that the four mentioned weights become zero and k 12 = 2δ f 2 -k 23 , allowing us to find f * * = 1 2 [(s e + s w ) w * * -s e w * e -s w w * w + f * e + f * w ] , (3.19a) g * * = 1 2 [(s n + s s ) w * * -s n w n * -s s w s * + g n * + g s * ] , (3.19b)

  Variation in the y-direction

Figure 8 :

 8 Figure 8: Two-dimensional variations restricted to a single spatial direction.

Figure 9 :

 9 Figure 9: A solution example at time t = ∆t, resulting from the application of HLL Riemann solvers at the interfaces and at the vertices of the underlying non-rectangular mesh.

Figure 10 :

 10 Figure 10: Density variable ρ obtained using BAL2012 (left, ρ:0.53-1.72) and our scheme (right, ρ:0.53-1.71) for the MultiD RP1; computations performed on a 400 × 400 grid and 30 contour lines displayed from 0.54 to 1.70 with a step of 0.04.

Inria Multidimensional Riemann Problems 1 & 2

 2 Quadrantρ 0 (x, y) u 0 (x, y) v 0 (x, y) p 0 (x, y) ρ 0 (x, y) u 0 (x, y) v 0 (x, y) p 0 (x, Computational domain: [-1, 1] × [-1, 1]; Free-flow boundary conditions

Figure 11 :

 11 Figure 11: The density computed with BAL2012 (left) and our scheme (right) for MultiD RP1, using 1000 × 1000 cells on 64 processors; contour lines chosen as in Fig. 10 (ρ:0.53-1.73).

Figure 15 :

 15 Figure15: First-order density results for MultiD RP3 obtained with our scheme on 4000 × 4000 (left, 512 processors) and 10000 × 10000 (right, 1024 processors) cells; contours chosen as in Fig.13(ρ:0.14-1.75).

Figure 16 :

 16 Figure 16: Results for the double Mach reflection problem obtained with BAL2012 (left) and our scheme (right), using 25 density contours ranging from 1.77 to 22.44 with a constant step; computations performed on 960×240 (top), 1920×480 (middle) and 2400×600 (bottom) meshes, 64 processors. All results have been plotted up to x = 3.

Figure 17 :

 17 Figure 17: Scatter plots for three different methods compared with the analytical time dependent solution of the density (top) and pressure (bottom); computations performed on a 65 × 65 grid.

Figure 18 :

 18 Figure 18: Density (top) and pressure (bottom) scatter plots compared with analytical solutions represented as solid lines, using three methods; computations performed on a 129 × 129 grid.

Algorithm 3 the fluxes φ hll2D x and φ hll2D y 1 :←←←←←

 31 Defining if (sw ≥ 0 and ss ≥ 0) then 2: φ hll2D x ← [(sn -ss) f * w + ss f sw ]/sn 3: φ hll2D y ← [(se -sw) g s * + sw g sw ]/se 4: else if (sw ≥ 0 and sn ≤ 0) then 5: φ hll2D x ← [(ss -sn) f * w + sn f nw ][(se -sw) g n * + sw g nw ]/se 7: else if (se ≤ 0 and ss ≥ 0) then 8: φ hll2D x ← [(sn -ss) f * e + ss f se ][(sw -se) g s * + se g se ]/sw 10: else if (se ≤ 0 and sn ≤ 0) then 11: φ hll2D x ← [(ss -sn) f * e + sn f ne ]/ss 12: φ hll2D y ← [(sw -se) g n * + se g ne ]/sw 13: else if sw ≥ 0 then [(se -sw) g * * + sw g * w ]/se 16: else if se ≤ 0 then [(sw -se) g * * + se g * e ]/sw 19: else if ss ≥ 0 then 20: φ hll2D x ← [(sn -ss) f * * + ss f s * ][(ss -sn) f * * + sn f n * ]

Algorithm 4 ←: φ hll2D yw ← g nw 14 :←

 414 Defining the fluxes φ = max(|sn|, |ss|) 2: sx = max(|se|, |sw|) 3: 4: if (sw ≥ 0 and ss ≥ 0) then 5: φ hll2D xn← [(sn -ss) f * w + ss f sw ]se -sw) g s * + sw g sw ]/se 8: φ hll2D yw ← g sw 9: else if (sw ≥ 0 and sn ≤ 0) [(ss -sn) f * w + sn f nw ]/ss 12: φ hll2D ye ← [(se -sw) g n * + sw g nw ]/se 13else if (se ≤ 0 and ss ≥ 0) then 15: φ hll2D xn ← [(sn -ss) f * e + ss f se ]sw -se) g s * + se g se ]/sw 19: else if (se ≤ 0 and sn ≤ 0) ss -sn) f * e + sn f ne ]sw -se) g n * + se g ne ]/sw 24: else if sw ≥ 0 then 25: φ hll2D xn ← [(sy -sn) f nw + sn f * w ]/sy 26: φ hll2D xs ← [(sy + ss) f sw -ss f * w ]/sy 27: φ hll2D ye ← [(se -sw) g * * + sw g * w ]sy -sn) f ne + sn f * e ]/sy 31: φ hll2D xs ← [(sy + ss) f se -ss f * e ][(sw -se) g * * + se g * e ]/sw 34: else if ss ≥ 0 then 35: φ hll2D xn ← [(sn -ss) f * * + ss f s * ]sx -se) g se + se g s * ]/sx 38: φ hll2D yw ← [(sx + sw) g sw -sw g s * ]/sx 39: else if sn ≤ 0 then ss -sn) f * * + sn f n * )/ss 42: φ hll2D ye ← [(sx -se) g ne + se g n * )/sx 43: φ hll2D yw ← [(sx + sw) g nw -sw g n * )/sx 44: else 45: φ hll2D xn ← [(sy -sn) f n * + sn f * * )/sy 46: φ hll2D xs ← [(sy + ss) f s * -ss f * * )/sy 47: φ hll2D ye ← [(sx -se) g * e + se g * * )/sx 48: φ hll2D yw ← [(sx + sw) g * w -sw g * * )/sx RR n°8540

  Solving the one-dimensional Riemann problems at the edges not only allows us to deduce the speeds (2.32) but also lets us derive expressions for the constant state w * (2.16) and HLL fluxes (2.23) associated with each problem. We denote by w µ * , µ = n or s, the state in the intermediate constant region of the one-dimensional Riemann problem with initial data set to w l = w µw and w r = w µe . Analogous notations are used for w * ν , ν = e or w, with initial states given by w l = w sν and w r = w nν . One then has

	s s w , s w s ), s se = (s s e , s e s ), s nw = (s n w , s w n ), s ne = (s n e , s e n ),	(2.32)
	where s n w and s n e (respectively, s s w and s s e ) are the left and right speeds for the one-dimensional
	Riemann problem above (respectively, below) the x-axis, with equivalent definitions for the y-
	direction. Additionally, one can specify the minimal and maximal wave speeds
	s n = max(s e n , s w n ), s s = min(s e s , s w s ), s e = max(s n e , s s e ), s w = min(s n w , s s w ).	(2.33)
	w µ * =	s µ e w µe -s µ w w µw + f µw -f µe s µ w e -s µ	, µ = n or s,	(2.34a)
	w * ν =	s ν n w nν -s ν s w sν + g sν -g nν s ν s n -s ν	,		ν = e or w,	(2.34b)
	with corresponding HLL fluxes
	f µ * =	s µ e f µw -s µ w f µe + s µ w s µ e (w µe -w µw ) s µ w e -s µ	, µ = n or s,	(2.35a)

  Substituting then the constant one-dimensional intermediate states w n * , w s * , w * e and w * w defined in (2.34), in conjunction with the particular values of all areas having the form a µν such as a * * (3.8), into the left-hand side of equation (3.7) yields

		2 2 [(s n e -s s w ) (s w n -s e s ) + (s e n -s w s ) (s s e -s n w )] .	(3.8)
	w * * =	1 s ) + (s e n -s e w ) (s w e -s s (s n n -s w s ) (s s e -s n w )	*
		[ (s w n s n e + s s e s e n )w ne -(s e n s n w + s s w s w n )w nw + (s e s s s w + s n w s w s )w sw -(s w s s s e + s n e s e s )w se
		-s w

n f ne + s e n f nw -s e s f sw + s w s f se -(s e n -s e s ) f * e + (s w n -s w s ) f * w -s s e g ne + s s w g nw -s n w g sw + s n e g se -(s n e -s n w ) g n * + (s s e -s s w ) g s * ] , (3.9) after considerable algebraic manipulation. This formulation is only complete after appropriate definitions for the fluxes f * e , f * w , g n * and g s * have been furnished. For this end, we adopt the solution (2.48) proposed by Balsara, thoroughly detailed in Section 2.2.2, and construct each of these unknown transverse fluxes with values obtained from the corresponding normal flux and intermediate state.

Scheme with Manual Asmb. Scheme with Simpson Asmb.

	Limiter	Resolution	L 1 error	L 1 order	L 1 error	L 1 order
	None	25 × 25	4.8975e-01		4.9536e-01
		50 × 50	4.1098e-01	0.25	4.2679e-01	0.21
		100 × 100	2.8650e-01	0.52	3.0515e-01	0.48
		200 × 200	1.7279e-01	0.73	1.8704e-01	0.71
		400 × 400	9.5487e-02	0.86	1.0429e-01	0.84
		800 × 800	5.0278e-02	0.93	5.5167e-02	0.92
		1600 × 1600	2.5808e-02	0.96	2.8383e-02	0.96
	MM [31]	25 × 25	2.0198e-01		2.1682e-01
		50 × 50	6.5074e-02	1.63	6.7657e-02	1.68
		100 × 100	2.8358e-02	1.20	3.0418e-02	1.15
		200 × 200	7.8803e-03	1.85	8.4898e-03	1.84
		400 × 400	2.1739e-03	1.86	2.3589e-03	1.85
		800 × 800	5.9648e-04	1.87	6.4860e-04	1.86
		1600 × 1600	1.5876e-04	1.91	1.7275e-04	1.91
	PP [34]	25 × 25	7.7416e-02		7.9507e-02
		50 × 50	2.5379e-02	1.61	2.3234e-02	1.77
		100 × 100	5.5017e-03	2.21	5.1137e-03	2.18
		200 × 200	1.1486e-03	2.26	1.0969e-03	2.22
		400 × 400	2.3008e-04	2.32	2.2597e-04	2.28
		800 × 800	4.5417e-05	2.34	4.5906e-05	2.30
		1600 × 1600	8.9790e-06	2.34	9.4001e-06	2.29
	MC [35]	25 × 25	3.2846e-02		4.2989e-02
		50 × 50	4.4552e-03	2.88	5.8172e-03	2.89
		100 × 100	8.8114e-04	2.34	1.0172e-03	2.52
		200 × 200	2.0742e-04	2.09	2.1494e-04	2.24
		400 × 400	4.8755e-05	2.09	4.8707e-05	2.14
		800 × 800	1.1452e-05	2.09	1.1523e-05	2.08
		1600 × 1600	2.7050e-06	2.08	2.7568e-06	2.06

Table 2 :

 2 Density errors measured in the L 1 norm for the wave advection test using the MC limiter and one-dimensional HLL Riemann solvers in both the predictor and corrector steps.

	1D HLL Solvers Only
		25 × 25	50 × 50	100 × 100	200 × 200	400 × 400	800 × 800	1600 × 1600
	L 1 error	4.2162e-02 6.5289e-03 1.3566e-03 3.5122e-04 8.7975e-05 2.1608e-05	5.2697e-06
	L 1 order		2.69	2.27	1.95	2.00	2.03	2.04

Table 4 :

 4 Initial data for the second MultiD RP described in

  , i.e., grids consisting of 960 × 240, 1920 × 480 and 2400 × 600 zones, spanning the domain [0, 4] × [0, 1].

					Density at t =0.20							Density at t =0.20
								22.44						22.44
								19.86						19.86
	0.8							0.8
								17.27						17.27
	0.6							14.69 0.6						14.69
	y						y	12.11						12.11
	0.4							9.52 0.4						9.52
								6.94						6.94
	0.2							0.2
								4.35						4.35
								1.77						1.77
	0.0	0.5	1.0	1.5	2.0	2.5	3.0	0.0	0.5	1.0	1.5	2.0	2.5	3.0
				x							x
					Density at t =0.20							Density at t =0.20
								22.44						22.44
								19.86						19.86
	0.8							0.8
								17.27						17.27
	0.6							14.69 0.6						14.69
	y						y	12.11						12.11
	0.4							9.52 0.4						9.52
								6.94						6.94
	0.2							0.2
								4.35						4.35
								1.77						1.77
	0.0	0.5	1.0	1.5	2.0	2.5	3.0	0.0	0.5	1.0	1.5	2.0	2.5	3.0
				x							x
					Density at t =0.20							Density at t =0.20
								22.44						22.44
								19.86						19.86
	0.8							0.8
								17.27						17.27
	0.6							14.69 0.6						14.69
	y						y	12.11						12.11
	0.4							9.52 0.4						9.52
								6.94						6.94
	0.2							0.2
								4.35						4.35
								1.77						1.77
	0.0	0.5	1.0	1.5	2.0	2.5	3.0	0.0	0.5	1.0	1.5	2.0	2.5	3.0
				x							x

RR n°8540

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

The advantage of the suggested formulation over existing ones is that it efficiently encloses all feasible subsonic or supersonic configurations for the two-dimensional interaction of waves associated with the Riemann problem (2.25, 2.28), while providing a single and perspicuous implementation of the approximated variables w * * (3.9) and F * * (3.17).

If we regard the elements of the matrix K as weights, we notice that k 11 , k 13 , k 22 and k 24 become smaller as the strongly interaction region in the t = ∆t plane turns rectangular. In fact, studying once more the situation discussed at the end of the previous section where this region is a rectangle R ′ , we perceive that in the limit s β α → s α for α, β ∈ {n, s, e, w}, δ

Inria

Form II

We shall now describe a method that is built with the specific intention of eliminating the contribution of the resolved state w * * to the flux tensor F * * . We start by summing equation (3.13a) multiplied by δ w 3 and equation (3.13c) multiplied by -δ w 1 , to get

and in an analogous manner, we multiply equation (3.13d) by δ w 2 and equation (3.13b) by -δ w 4 so that their sum gives

(3.24)

Using the same methodology as in Form I, we employ matrix notation to write both linear equations as

Here, we are interested in finding the determinant of M , introduced in Section 3.2.2, and analyzing the result. We begin by evaluating the matrix product A T A in order to get

which is a two-dimensional square matrix with a straightforward determinant

2) after expansion and some simplification. It is well-known from basic linear algebra that M is nonsingular if and only if its determinant is nonzero. Thus, we propose to rewrite equation (A.2) as a sum of squared binomials

that is obviously greater or equal to zero. It suffices then to find a term that is greater than zero to prove the determinant is strictly positive. For this, we consider the squared binomial r ne = (δ f 4 δ g 1 -δ f 1 δ g 4 )

2 and acknowledge the following: δ f 4 = s e n -s e s > 0, recalling that s e n and s e s are the distinct maximum and minimum signal velocities (2.24) associated with the one-dimensional Riemann problem on the right of the y-axis, and δ g 1 = s n e -s n w > 0, using an analogous reasoning. Therefore, we are certain that the product δ f 4 δ g 1 > 0 and the only way that r ne would become zero is if

Next, we denote by r sw the term (δ f 3 δ g 2 -δ f 2 δ g 3 ) 2 and using an approach similar to the previous one, we find that r sw

given that δ f 2 δ g 3 > 0 since both δ g 3 , δ f 2 < 0. Two simple but useful properties that hold for the absolute value are |ab| = |a||b|, for any a, b ∈ R, and ab ≤ |a||b|, so that we can obtain

By performing some algebra, we find that s = 2 ∆t 2 a * * , and, provided the assumption ∆t > 0 is satisfied, we recover r en + r sw > 0 as long as a * * = 0. In other words, if a * * = 0, the determinant of M is strictly positive and the matrix is invertible

B Pseudocode for Manual Assembling

With regard to equation (3.44) and its analogue in the y direction, we present three snippets of pseudocode that will help the reader compute the values for the speeds sα , with α ∈ {n, s, e, w}, and for the "hll2D" fluxes.