Dinshaw S Balsara
email: dbalsara@nd.edu

Boniface Nkonga
email: boniface.nkonga@unice.fr

Multidimensional Riemann Problem with Self-Similar Internal Structure -Part III-A Multidimensional Analogue of the HLLI Riemann Solver for Conservative Hyperbolic Systems By

Just as the quality of a one-dimensional approximate Riemann solver is improved by the inclusion of internal sub-structure, the quality of a multidimensional Riemann solver is also similarly improved. Such multidimensional Riemann problems arise when multiple states come together at the vertex of a mesh. The interaction of the resulting one-dimensional Riemann problems gives rise to a strongly-interacting state. We wish to endow this strongly-interacting state with physically-motivated sub-structure. The fastest way of endowing such sub-structure consists of making a multidimensional extension of the HLLI Riemann solver for hyperbolic conservation laws. Presenting such a multidimensional analogue of the HLLI Riemann solver with linear sub-structure for use on structured meshes is the goal of this work. The multidimensional MuSIC Riemann solver documented here is universal in the sense that it can be applied to any hyperbolic conservation law.

The multidimensional Riemann solver is made to be consistent with constraints that emerge naturally from the Galerkin projection of the self-similar states within the wave model. When the full eigenstructure in both directions is used in the present Riemann solver, it becomes a complete Riemann solver in a multidimensional sense. I.e., all the intermediate waves are represented in the multidimensional wave model. The work also presents, for the very first time, an important analysis of the dissipation characteristics of multidimensional Riemann solvers. The present Riemann solver results in the most efficient implementation of a multidimensional Riemann solver with sub-structure. Because it preserves stationary linearly degenerate waves, it might also help with well-balancing. Implementation-related details are presented in pointwise fashion for the one-dimensional HLLI Riemann solver as well as the multidimensional MuSIC Riemann solver.

Several stringent test problems drawn from hydrodynamics, MHD and relativistic MHD are presented to show that the method works very well on structured meshes. Our results demonstrate the versatility of our method. The reader is also invited to watch a video

I) Introduction

One-dimensional Riemann solvers are routinely used in the numerical solution of hyperbolic systems of conservation laws. The one-dimensional Riemann problem is a selfsimilar solution that results from a discontinuity between two constant states. In their numerical study of the multidimensional Riemann problem, Schulz-Rinne, Collins & Glaz [START_REF] Schulz-Rinne | Numerical solution of the Riemann problem for two-dimensional gas dynamics[END_REF] initialized four states around the center of a two-dimensional Cartesian mesh. While one-dimensional Riemann problems arise between each pair of states, those authors showed that the onedimensional Riemann problems interact amongst themselves to form a self-similarly evolving strongly-interacting state. This strongly-interacting state arises at the point where the four states come together. The study of the multidimensional Riemann problem is, therefore, the study of the strongly-interacting state. This strongly-interacting state emerges by propagating into the one-dimensional Riemann problems along its boundary. Consequently, the strongly interacting state, as well as the one-dimensional Riemann problems that form its boundary, evolve in a selfsimilar fashion. We refer to this boundary as the boundary of the multidimensional wave model because it contains the strongly-interacting state. The wave models in all the multidimensional Riemann solvers incorporate this concept of self-similarity. Schulz-Rinne, Collins & Glaz [START_REF] Schulz-Rinne | Numerical solution of the Riemann problem for two-dimensional gas dynamics[END_REF] only presented a computational study of the multidimensional Riemann problem. However, Abgrall [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] was the first to formulate multidimensional Riemann solvers that were usable. The self-similarly evolving strongly-interacting state is an inevitable consequence of having a multidimensional wave model that propagates into the one-dimensional Riemann problems. Seizing on this insight, Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF] presented a self-similar formulation of the multidimensional Riemann problem. Incorporating the physics of the strongly-interacting state has shown to be very advantageous in second order calculations (Balsara [4]) and higher order accurate calculations (Balsara [15]). This is the true motivation for our study of the multidimensional Riemann solver reported here.

Following Abgrall [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], further advances were also reported (Fey [40], [START_REF] Fey | Multidimensional upwinding 2. Decomposition of the Euler equation into advection equation[END_REF], Gilquin, Laurens & Rosier [START_REF] Gilquin | Multidimensional Riemann problems for linear hyperbolic systems[END_REF], Brio, Zakharian & Webb [START_REF] Brio | Two-dimensional Riemann solver for Euler equations of gas dynamics[END_REF]). However, these early formulations were cumbersome and did not see much use. Multidimensional Riemann solvers that are very efficient have also been designed and we focus on a certain class of multidimensional Riemann solvers here (Wendroff [70], 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], [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF], [START_REF] Balsara | Three Dimensional HLL Riemann Solver for Structured Meshes; Application to Euler and MHD Flow[END_REF], Balsara, Dumbser & Abgrall [START_REF] Balsara | Multidimensional HLL and HLLC Riemann Solvers for Unstructured Meshes -With Application to Euler and MHD Flows[END_REF], Vides, Nkonga & Audit [START_REF] Wendroff | A two-dimensional HLLE Riemann solver and associated Godunov-type difference scheme for gas dynamics[END_REF], Balsara & Dumbser [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part II -Application to Hyperbolic Conservation Laws on Unstructured Meshes[END_REF], Balsara et al. [START_REF] Balsara | A Two-Dimensional Riemann Solver with Self-Similar Sub-Structure -Alternative Formulation Based on Least Squares Projection[END_REF]). A video introduction to multidimensional Riemann solvers is available on the following website: http://www.nd.edu/~dbalsara/Numerical-PDE-Course . Such Riemann solvers are applied at the vertices of a two-dimensional or three-dimensional mesh. Many states come together at a vertex from different directions, making it possible to communicate the multidimensionality of the flow to the multidimensional Riemann solver. At the vertex, the job of the multidimensional Riemann solver is to approximate the self-similar multidimensional structure that emanates from the vertex.

By this point in time, there has been substantial progress in one-dimensional and multidimensional Riemann solvers. In this paragraph we list the one-dimensional Riemann solvers and juxtapose them with their multidimensional counterparts. Such a juxtaposition can be very useful in building perspective. Several excellent one-dimensional Riemann solvers have been designed. There are exact Riemann solvers from Godunov [START_REF] Godunov | Finite Difference Methods for the Computation of Discontinuous Solutions of the Equations of Fluid Dynamics[END_REF], [START_REF] Godunov | Numerical Solution of Multi-dimensional Problems in Gas Dynamics[END_REF] and van Leer [START_REF] Van Leer | Toward the Ultimate Conservative Difference Scheme. V. A Second-Order Sequel to Godunov's Method[END_REF] and two-shock approximations thereof (Colella [START_REF] Colella | A direct Eulerian MUSCL scheme for gas dynamics[END_REF], Colella & Woodward [START_REF] Colella | The piecewise parabolic method (PPM) for gas-dynamical simulations[END_REF]). See also the work of Chorin [START_REF] Chorin | Random choice solutions of hyperbolic systems[END_REF]. The linearized Riemann solver by Roe [START_REF] Roe | Approximate Riemann solver, parameter vectors and difference schemes[END_REF] has also proved useful. The multidimensional Riemann solver by Abgrall [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] can be viewed as Roe-type Riemann solver that has been extended to multiple dimensions. One-dimensional HLL Riemann solvers (Harten, Lax & van Leer [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]) have now been extended to two-dimensions (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]) and threedimensions (Balsara [18]). The papers by Balsara offer simple closed form expressions for the multidimensional HLL fluxes that are easy to implement. One-dimensional HLLC Riemann solvers (Toro, Spruce and Speares [START_REF] Toro | Restoration of contact surface in the HLL Riemann solver[END_REF] [66], [START_REF] E F Toro | Restoration of the contact surface in the HLL Riemann solver[END_REF], Chakraborty & Toro [START_REF] Chakraborty | Development of an approximate Riemann solver for the steady supersonic Euler equations[END_REF] and Batten et al. [START_REF] Batten | On the choice of wavespeeds for the HLLC Riemann solver[END_REF]) seek to restore the physics of the contact discontinuity. Multidimensional extensions of the HLLC Riemann solver to structured and unstructured meshes have also become available in recent papers (Balsara [4], Balsara, Dumbser & Abgrall [START_REF] Balsara | Multidimensional HLL and HLLC Riemann Solvers for Unstructured Meshes -With Application to Euler and MHD Flows[END_REF]). While HLLC Riemann solvers seek to restore an isolated contact discontinuity in the HLL Riemann solver, it is always interesting to ask if there are other ways to introduce an intermediate wave into the HLL Riemann solver? The one-dimensional HLLE/HLLEM Riemann solver (Einfeldt [38], Einfeldt et al. [START_REF] Einfeldt | On Godunov-type methods near low densities[END_REF]) tried to do that by introducing a linear profile in the Riemann fan. However, because of an error in the formulation, it did not achieve its intended goal. Dumbser & Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF] rectified the prior deficiencies and also introduced another very important advance. Using the self-similar formulation of Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF], they were able to introduce multiple intermediate waves into the HLL Riemann solver, thus giving rise to the HLLI Riemann solver. Here "I" stands for intermediate waves and acknowledges the fact that the HLLI Riemann solver can accommodate any intermediate wave as long as its eigenstructure is known. The result is a one-dimensional HLLI Riemann solver that benefits from all the good properties of the one-dimensional HLL Riemann solver and simultaneously functions as a Riemann solver that retains sub-structure. When all the intermediate waves are included, the one-dimensional HLLI Riemann solver of Dumbser & Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF] becomes a complete Riemann solver. It is also a fully capable replacement for costlier Riemann solvers by Osher and Solomon [START_REF] Osher | Upwind Difference Schemes for Hyperbolic Systems of Conservation Laws[END_REF] and Dumbser and Toro [START_REF] Dumbser | A simple extension of the Osher Riemann solver to nonconservative hyperbolic systems[END_REF]. It is, therefore, very attractive to present a two-dimensional analogue of the HLLI Riemann solver for hyperbolic conservation laws and that is indeed the first goal of this paper. Such a multidimensional Riemann solver can be made complete in a multidimensional sense if all the intermediate waves in all directions are included. This is a very attractive property and we explore it further in this paper.

Self-similarity has not been used much in the design of one-dimensional Riemann solvers; the only real exception being the HLLI Riemann solver of Dumbser & Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF]. However, it is crucially important in the development of multidimensional Riemann solvers (Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF], Balsara & Dumbser [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part II -Application to Hyperbolic Conservation Laws on Unstructured Meshes[END_REF]). This has prompted the name of MuSIC Riemann solvers, where MuSIC stands for "Multidimensional, Self-similar, strongly-Interacting, Consistent". Such Riemann solvers are multidimensional; they draw on the self-similarity of the problem; they focus on the strongly-interacting state that results when multiple one-dimensional Riemann solvers interact; and the design relies on establishing consistency with the conservation law. MuSIC Riemann solvers that rely on a Petrov-Galerkin projection to obtain the self-similar variation in the strongly interacting state have been presented (Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF], Balsara & Dumbser [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part II -Application to Hyperbolic Conservation Laws on Unstructured Meshes[END_REF]). An alternative projection method consists of satisfying the one-dimensional shock jumps at the boundary of the multidimensional wave model. Vides, Nkonga & Audit [START_REF] Wendroff | A two-dimensional HLLE Riemann solver and associated Godunov-type difference scheme for gas dynamics[END_REF] and Balsara et al. [START_REF] Balsara | A Two-Dimensional Riemann Solver with Self-Similar Sub-Structure -Alternative Formulation Based on Least Squares Projection[END_REF] developed a multidimensional Riemann solver without and with sub-structure respectively that uses least squares minimization methods. A study of the dissipation characteristics of MuSIC Riemann solvers has never been presented. The second goal of this paper is to present a thorough study of the dissipation characteristics of the MuSIC Riemann solvers. We first present an analysis of the dissipation characteristics of the one-dimensional HLLI Riemann solver. We then show that when the one-dimensional HLLI Riemann solver is used as a building block for the MuSIC Riemann solver, its dissipation characteristics mirror those of the HLLI Riemann solver for flows that are mesh-aligned.

It is also worth recalling that the one-dimensional HLLI Riemann solver of Dumbser & Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF] is a universal Riemann solver; i.e. it is applicable to any hyperbolic conservation law. It would be very desirable to have a multidimensional Riemann solver that is also applicable to any conservation law. The third goal of this paper is to show that when the one-dimensional HLLI Riemann solver is used as a building block for the MuSIC Riemann solver we indeed get a universal multidimensional Riemann solver that works for any hyperbolic conservation law. This generality implies that multidimensional Riemann solvers with sub-structure can be built and incorporated into any code for any hyperbolic conservation law. Moreover, the same coding strategy can be used for all hyperbolic conservation laws. Magnetohydrodynamics (MHD) is an interesting example of a hyperbolic system with a more complex wave foliation. One-dimensional linearized Riemann solvers for numerical MHD have been designed (Roe & Balsara [60], Cargo and Gallice [START_REF] Cargo | Roe matrices for ideal MHD and systematic construction of Roe matrices for systems of conservation laws[END_REF], Balsara [START_REF] Balsara | Linearized formulation of the Riemann problem for adiabatic and isothermal magnetohydrodynamics[END_REF]). HLLC Riemann solvers, capable of capturing mesh-aligned contact discontinuities, have been presented by Gurski [START_REF] Gurski | An HLLC-type approximate Riemann solver for ideal magnetohydrodynamics[END_REF] and Li [START_REF] Li | An HLLC Riemann solver for magnetohydrodynamics[END_REF]. Miyoshi and Kusano [START_REF] Miyoshi | A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics[END_REF] drew on Gurski's work to design an HLLD Riemann solver for MHD. It is, therefore, interesting to show that MHD can also be accommodated within our formulation. MHD is a system with an involution constraint, where the divergence of the magnetic field is always zero. Balsara & Spicer [START_REF] Balsara | A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations[END_REF] showed that this is assured within the context of a higher order Godunov scheme by using the upwinded fluxes at the edges of the mesh to update the magnetic fields that are collocated at the faces of a mesh. Gardiner & Stone [START_REF] Gardiner | An unsplit Godunov method for ideal MHD via constrained transport[END_REF], [START_REF] Gardiner | An unsplit Godunov method for ideal MHD via constrained transport in three dimensions[END_REF] have claimed that the dissipation in those upwinded fluxes needs to be doubled all the time in order to stabilize the method. A substantial body of work now exists to show that the suggestion of Gardiner & Stone is completely unnecessary when multidimensional Riemann solvers are used to provide a properly upwinded electric field at the edges of the mesh (Balsara [4], Vides, Nkonga & Audit [START_REF] Vides | A simple two-dimensional extension of the HLL Riemann solver for hyperbolic conservation laws[END_REF], Balsara & Dumbser [START_REF] Balsara | Divergence-Free MHD on Unstructured Meshes Using High Order Finite Volume Schemes Based on Multidimensional Riemann Solvers[END_REF]). Indiscriminate doubling of the dissipation, as per Gardiner & Stone's suggestion, can indeed lead to excessive dissipation of the magnetic field in the direction that is transverse to the upwind direction. The present paper reinforces that finding.

As with classical MHD, progress has also been made in relativistic MHD (RMHD). Balsara [START_REF] Balsara | Total Variation Diminishing Scheme for Relativistic Magnetohydrodynamics[END_REF] and Komissarov [START_REF] Komissarov | A Godunov-type scheme for Relativistic MHD[END_REF] have designed Roe-type Riemann solvers for RMHD. HLLC and HLLD type Riemann solvers for RMHD have also been designed by Mignone & Bodo [START_REF] Mignone | An HLLC Riemann Solver for Relativistic Flows II -Magnetohydrodynamics[END_REF], Honkkila & Janhunen [START_REF] Honkkila | HLLC solver for relativistic MHD[END_REF], Mignone, Ugliano and Bodo [START_REF] Mignone | A five-wave HLL Riemann solver for Relativistic MHD[END_REF] and Kim & Balsara [START_REF] Kim | A Stable HLLC Riemann solver for Relativistic Magnetohydrodynamics[END_REF]. Balsara and Kim [START_REF] Balsara | A subluminal relativistic magnetohydrodynamics scheme with ADER-WENO predictor and multidimensional Riemann solver-based corrector[END_REF] have also shown the value of multidimensional Riemann solvers for RMHD calculations. The present paper reinforces the utility of MuSIC Riemann solvers for accurate RMHD simulations.

Section II describes a one-dimensional HLLI Riemann solver for conservation laws that is indeed novel and has some rather nice properties. Section III provides details associated with the construction of the multidimensional Riemann problem on Cartesian meshes. Section IV shows that schemes that use the multidimensional Riemann solver meet their design accuracy. Section V shows the results of several stringent test problems drawn from Euler, MHD and relativistic MHD flow. Section VI presents conclusions.

II) Quick Derivation of the One-Dimensional HLLI Riemann Solver

In any multidimensional Riemann problem, the strongly-interacting state propagates into a sequence of one-dimensional Riemann problems that lie on its boundary. One dimensional Riemann solvers are, therefore, used as building blocks for the multidimensional Riemann problem. Because we wish to show that the dissipation characteristics of the MuSIC Riemann solver strongly mirror those of the one-dimensional HLLI Riemann solver, we first present a quick derivation of the one-dimensional HLLI Riemann solver in Sub-section II.a and study its dissipation characteristics in Sub-section II.b. This study is somewhat different from the one presented in Dumbser & Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF] because the prior work did not use one of the Galerkin constraints that results from the imposition of self-similarity. A compare-and-contrast is presented in Sub-section II.c. Sub-section II.d presents implementation-related details. In Section III we present a multidimensional Riemann solver in two-dimensions that is a close analogue of the one-dimensional HLLI Riemann solver presented here when Cartesian meshes are used.

II.a) Galerkin Formulation in Similarity Variables

In this section we consider an N-component hyperbolic conservation law,

0 t x ∂ ∂ + ∂ ∂ = U F
, which is restricted to one dimension.

c c R L R L S S S S ξ ξ ξ ξ ξ ξ   - ≡ ≡ + ∆ ≡ -   ∆    (2.1)
Since the solution evolves self-similarly within the Riemann fan, the solution within the Riemann fan can be written in terms of similarity variables. Written in these shifted similarity variables, the conservation law becomes ( )

1 0 c ξ ξ ξ ξ ξ   ∂ - + ∆   + = ∆ ∂ F U U    (2.
2)

The tilde on the top of U  is intended to signify a self-similarly evolving solution. The same is true for F  . Because of self-similarity, U  and F  are functions of only one similarity variable ξ . Eqn. (2.2) is then the governing equation written in terms of the similarity variable. We expand our state and flux as

( ) ξ ξ ξ = + U U U  (2.3) and ( ) ( ) with ξ ξ ξ ∂ = + = ∂ F U F F A U A U  (2.4)
Please note that we have evaluated the characteristic matrix A by using the mean state U ; but there is some flexibility in the evaluation of the characteristic matrix. For example, it can be evaluated using Roe-averages or arithmetic averages, as was done in Dumbser and Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF].

Please also note that ξ ξ ≠ ∆ ∆ U U where ( )

R L ∆ ≡ - U U U . Realize that ξ ∆ ∆ U
is indeed an estimate of the full gradient and, therefore, includes contributions from the extremal waves that make up the Riemann fan. In a numerical Riemann problem, we only want to pick out contributions from waves that are internal to the Riemann fan. We will soon show that ξ U will be obtained by a projection of ∆U onto the subset of waves that are interior to the Riemann fan.

Multiplying the conservation law from eqn. (2.2) with the test function ( )

φ ξ gives ( ) ( ) { } ( ) ( ) ( ) 1 1 0 c c φ ξ ξ ξ ξ φ ξ ξ ξ ξ φ ξ ξ ξ ξ ξ   ∂ - + ∆ ∂     - - + ∆ + =   ∆ ∂ ∆ ∂ F U F U U      (2.5)
Now we are ready to make Galerkin projections with different test functions.

Using ( ) ) ( )

1 1 HLL R R R L L L S S ξ ξ = = - - + - ∆ ∆ U U F U F U (2.6)
In practice, one always evaluates HLL U at the start of the calculation because it plays an important role in the rest of the calculation. This could include the construction of the characteristic matrix A . Realize, therefore, that HLL U from the equation above will always be a positivity-preserving state. Using ( ) φ ξ ξ = and making a Galerkin projection gives ( ) ( )

1 1 1 1 6 2 2 6 c HLL R R R L L L HLL S S ξ ξ ξ ξ ξ ξ ξ ξ ∆ - = - - - - - ⇔ = + ∆ ∆ ∆ ∆ U F U F U F U F F U (2.7)
Here in which case all the waves in the hyperbolic system are considered. Consequently, the Riemann solver becomes a complete Riemann solver.) The best characteristic projection we can do gives us

( ) ( ) ( ) int 2 2 i i R L i R L i I l r ξ δ ∈ = -= -     ∑ U U U R δ L U U  (2.8)
Here R is a matrix of right eigenvectors with dimension Also please notice that when the state is endowed with sub-structure F , which is obtained from eqn. (2.7), is not the classical HLL flux. The final numerical flux at the zone boundary, i.e. at 0 ξ =  , is given by ( )

c numerical c ξ ξ ξ ξ ξ ξ = = - ∆ =- ∆ F F F A U  (2.9a) or ( ) ( ) int int = 2 2 6 c c numerical i i R L i i HLL i i R L i i i I i I l r l r ξ ξ ξ ξ δ λ δ λ ξ ξ ∈ ∈       ∆     - - = + + - -               ∆ ∆           ∑ ∑ F F U U F U U U   (2.9b) or ( ) ( ) ( ) 2 1 = 2 2 2 3 c numerical HLL R L ξ ξ ξ   ∆ - - + -   ∆   F F R δ Λ δ L U U (2.
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 2 1 2 2 2 3 numerical R L R L R L R L R L R L R L R L R L

S S S S S S S S S S S S S S

= +

  + - + - - - + -   - - -   F F F R Λ I δ Λ δ L U U (2.10)
The second term in the above equation helps us to identify the viscosity of our Riemann solver. The square bracket in the above equation gives us the eigenvalues of the viscosity matrix and we want these to be bounded by the eigenvalues of the Roe-matrix viscosity (at the lower end) and the eigenvalues of the HLL viscosity (at the upper end). Using the dissipation properties of the underlying HLL Riemann solver we get the condition for i δ as follows Also notice that when

( ) ( ) ( ) ( ) ( )
i R L i R L i R i L i R L i R L i R L R L i R L S S S S S S S S S S S S S S S S φ λ δ λ λ φ λ φ - +    - - + ≤       =  + -           - - +      - - (2
R L S S  or when L R S S 
, we have

0 i δ → so that numerical HLL → F
F . Now recall the very nice design feature of the HLL Riemann solver which says that the subsonic flux retrieves the supersonic fluxes when the Riemann fan is opened up ever so slightly so as to always force it to be minimally subsonic. From the property stated at the beginning of this paragraph we see that our HLLI Riemann solver also retains that very nice design feature.

II.b) Dissipation Properties of the present HLLI Riemann solver

Recall that the Roe-type Riemann solver provides the theoretical minimum dissipation that any Riemann solver should provide to a scheme in order to ensure stability of the numerical method. However, the Roe-type Riemann solver has problems with positivity enforcement, while the HLLI Riemann solver discussed in this Section is free of this problem. The entropy fix is also naturally built into the HLLI Riemann solver. It is, therefore, worth asking the question, "How much excess dissipation is produced by the present HLLI Riemann solver compared to the Roetype Riemann solver?". We answer that question in this paragraph and the next one. To normalize the search space, we can always require 1 R L S S -= . We also require 0

L R S S ≤ ≤ , i.e.
we focus on the subsonic case. We assume that there is only one intermediate wave with wave speed i λ such that L i R S S λ ≤ ≤ . (Since the dissipation is independently determined for each wave family, the number of wave families that we use does not affect our present analysis.) For such a wave, we can use eqn. (2.11) to evaluate i δ . The square bracket in eqn. (2.10) then gives us the dissipation matrix. The diagonal term in the dissipation matrix for the intermediate wave being considered should be greater than or equal to i λ because this is the theoretically minimum amount of dissipation required by the Roe-type Riemann solver. For various subsonic choices of L S and R S , and with the normalizing restriction 1 R L S S -= , we can indeed step through all possible values of i λ . We can then plot the dissipation produced by the present HLLI-style Riemann solver versus i λ . We can also plot i λ , the dissipation from the Roe-type Riemann solver, versus i λ . Such an exercise is undertaken in the next paragraph and it enables us to get an interesting perspective on the dissipation characteristics of the present HLLI Riemann solver vis a vis the Roe-type Riemann solver.

The previous paragraph outlined a strategy for quantifying the dissipation properties of the HLLI Riemann solver and comparing it to the Roe-type Riemann solver. The results of this exercise are shown in Fig. 1. The dashed lines in Fig. 1 show the dissipation from our HLLI Riemann solver whereas the solid lines show the dissipation from the Roe-type Riemann solver. Fig. 1a shows us the dissipation from the HLLI Riemann solver and also the theoretically minimum dissipation, i λ , on the vertical axis as a function of wave speed, i λ , on the horizontal axis when 0.9 ; we see that the dissipation of the two Riemann solvers is almost identical. Based on such an analysis, we conclude that our present HLLI Riemann solver always produces dissipation that is within ten to twenty percent of the Roe-type Riemann solver under all circumstances. In many of the situations shown in Fig. 1, the two Riemann solvers have identical dissipation. This is a very interesting demonstration in light of the versatility, robustness and favorable positivity properties of our HLLI Riemann solver and the lack thereof for the Roe-type Riemann solver. For all the panels in Fig. 1 we see that our HLLI Riemann solver has zero dissipation when 0

i λ = which
shows that it can also capture stationary linearly degenerate waves exactly. Consequently, we see that it offers all the good attributes of the Roe-type Riemann solver while avoiding all its pitfalls. Fig. 1 of this paper can also be compared to Fig. 3.1 of Castro-Díaz and Fernández-Nieto [START_REF] Castro-Díaz | A class of computationally fast first order finite volume solvers: PVM methods[END_REF] if one wants to analyze the dissipation properties of the HLLI Riemann solver through the perspective of polynomial viscosity methods (PVM).

II.c) Comparison with the HLLI-type Riemann solver of Dumbser and Balsara [37]

In this section we have designed an HLLI-type Riemann solver based on endowing substructure to the HLL Riemann solver. The one-dimensional HLLI Riemann solver described here is very useful because it extends more naturally to multidimensions. In Dumbser and Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF] a slightly different HLLI-type Riemann solver had been presented. The difference is primarily in the fact that the Galerkin projection in eqn. (2.7) is not used in the design of the Riemann solver in Dumbser & Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF]. As a result, eqns. (2.9) and (2.11) are also substantially different. It is interesting to compare and contrast the two variants of HLLI Riemann solvers. To that end, it is valuable to write the explicit expressions for ( )

int int int 2 2 2 6 c HLL i i R L i i I c HLL i i R L i i i R L i i i I i I l r l r l r ξ ξ ξ δ ξ ξ ξ ξ ξ δ δ λ ξ ∈ ∈ ∈     -   = + -         ∆            - ∆     = + - + -               ∆           ∑ ∑ ∑ U U U U F F U U U U          ( ) ( ) ( ) ( ) ( )
The two curly brackets in the above two equations only need to be evaluated once. Appendix A provides the corresponding formulation of this Riemann solver for moving meshes, i.e., ALEtype meshes. Notice that the Galerkin formulation from eqn. (2.7) dictates that the inclusion of sub-structure should cause a change in the mean flux in eqn. (2.12). It is also helpful to be able to compare and contrast this Riemann solver with the HLLI Riemann solver from Dumbser and Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF]. That Riemann solver does not use the first moment of the conservation law, i.e. eqn. (2.7), because it is meant to be generally applicable to hyperbolic systems in conservation and non-conservative forms. As a result, the definition of i δ changes to

1 i i i L R S S λ λ δ - + =- - (2.13)
For the Riemann solver from Dumbser and Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF] we then have

( ) ( ) ( ) ( ) int int 2 2 2 2 c HLL i i R L i i I R L HLL i i R L i i I l r S S l r ξ ξ ξ δ ξ ξ ξ δ ξ ∈ ∈     -   = + -         ∆          -   = + -         ∆      ∑ ∑ U U U U F F U U         (2.14)
Only one curly bracket needs to be evaluated in the above equation, therefore, the HLLI Riemann solver from Dumbser and Balsara [START_REF] Dumbser | A New, Efficient Formulation of the HLLEM Riemann Solver for General Conservative and Non-Conservative Hyperbolic Systems[END_REF] has slightly lower computational complexity. However, both flavors of HLLI Riemann solvers require the evaluation of the intermediate eigenvectors. This eigenvector evaluation often constitutes the bulk of the additional computational cost that is added on top of the cost of the HLL Riemann solver. For that reason, both flavors of HLLI Riemann solver have almost the same overall computational complexity. Notice that i δ can assume larger values in eqn. (2.13) compared to eqn. (2.11). However, the mean HLL flux in eqn. (2.12) and (2.14) we see that the flux in the former equation varies linearly with ξ  whereas the flux in the latter equation varies quadratically with ξ  . This is a consequence of the different philosophies that were used in deriving the two variants of the HLLI Riemann solver. In practice, both work equally well. Both can preserve stationary intermediate waves on a mesh without additional dissipation.

II.d) Implementation-Related Details for one-dimensional HLLI Riemann solver:-

The present HLLI Riemann solver can be easily retrofitted to any HLL Riemann solver and usually provides a very palpable improvement in the simulation quality. The steps in the implementation of this one-dimensional HLLI Riemann solver are as follows: 7) The supersonic cases are obvious.

In this section we have provided details for the one-dimensional HLLI Riemann solver on a fixed mesh. But we also realize that some people might want to apply this Riemann solver to an arbitrary Lagrangian-Eulerian (ALE) mesh. Appendix A gives the formulation of the present one-dimensional HLLI Riemann solver on a moving mesh.

III) MuSIC Riemann Solver that is closest to an HLLI Formulation -Focus on Cartesian Meshes

Sub-section III.a presents the formulation of the MuSIC Riemann solver, including a description of the inclusion of sub-structure. Section III.b presents implementation-related details.

III.a) Formulation of the MuSIC Riemann Solver

Consider the N-component hyperbolic conservation law in two-dimensions, given by

0 t x y ∂ ∂ + ∂ ∂ + ∂ ∂ = U F G
. It can give rise to one-dimensional Riemann problems, but it can also give rise to a multidimensional Riemann problem. The multidimensional Riemann problem is most easily understood on a Cartesian mesh, and we focus on that in this paper because it is possible to get exact answers for the multidimensional Riemann solver in Cartesian geometry. We will defer the inclusion of sub-structure in the multidimensional Riemann solver on unstructured meshes for a subsequent paper. As shown schematically in Fig. 2 a multidimensional Riemann problem arises when four states RU U , LU U , LD U and RD U come together at a zone vertex; the vertex is shown as a gray dot in that figure. The four pairs of mutually contiguous states set up four one-dimensional Riemann problems. However, the onedimensional Riemann problems interact in a strongly-interacting state, as shown in Fig. 2a of Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF]. The strongly interacting state is bounded by a multidimensional wave model. In fig.

2a the thick solid line denotes the boundary of the multidimensional wave model; the interior of the wave model is shaded. The four initial states that come together at a vertex "O" of the mesh are also shown. The thin solid lines in Fig. 2a show the extremal speeds of the one-dimensional Riemann problems in the boundary of the multidimensional wave model. The dashed lines in Fig. 2a show the coordinate axes, measured as speeds. The bounding speeds of the multidimensional wave model are also shown. On such a mesh, the extent of the

multidimensional wave model, [ ] [ ] , , L R D U S S S S ×
, is approximated beforehand. See Balsara [START_REF] Balsara | Multidimensional HLLE Riemann solver; Application to Euler and Magnetohydrodynamic Flows[END_REF] and [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF] for advice on how to pick out the extent of the multidimensional wave model on a Cartesian mesh. The strongly-interacting state is bounded by the multidimensional wave model and evolves self-similarly, just like the one-dimensional Riemann problems at its boundary.

We want to predict the self-similar evolution of the multidimensional, stronglyinteracting state, U  . The tilde on the top of U  is intended to signify a self-similarly evolving solution. Let us, therefore, pick similarity variables in two-dimensions and express the stronglyinteracting state in terms of those two variables. The similarity variables are

) ( ) ( ) ( )
2 ; ; 2 ; ; c R L R L c U D U D c c S S S S S S S S ξ ξ ψ ψ ξ ξ ψ ψ ξ ψ ξ ψ ≡ + ∆ ≡ - ≡ + ∆ ≡ -   -  - ≡ ≡     ∆ ∆       (3.2)
Observe that ξ and ψ are still self-similar variables with the main difference that they now

range over [ ] [ ] 1/ 2,1/ 2 1/ 2,1/ 2 - × -
. This makes it easier to achieve concordance with the onedimensional case described in the previous sub-section. With the change of variables in eqn.

(3.2), the N-component conservation law in two-dimensions becomes ( ) ( )

1 1 2 0 c c ψ ψ ψ ξ ξ ξ ξ ξ ψ ψ     ∂ - + ∆ ∂ - + ∆     + + = ∆ ∂ ∆ ∂ G U F U U      (3.3)
Here the strongly-interacting state ( )

, ξ ψ = U U  
is a function of the two similarity variables. The same is true for the fluxes F  and G  .

We can now expand the strongly-interacting state in the similarity variables as ( )

, ξ ψ ξ ψ ξ ψ =+ + U U U U  (3.4)
Because ( ) ( ) , 0, 0 ξ ψ = corresponds to the centroid of our wave model, U is indeed the mean value associated with our wave model. The x-flux is written in similarity variables as

( ) ( ) ( ) , with ξ ψ ξ ψ ξ ψ ∂ = + + = ∂ F U F F A U U A U  (3.5)
It may also prove convenient to integrate eqn. (3.5) in the ψ -direction to write the numerical x- flux as ( )

1/ 2 1/ 2 , c numerical c d ξ ξ ξ ξ ξ ψ ψ ξ -   = = -∆ = -  ∆   ∫ F F F A U  (3.6)
The y-flux is written in similarity variables as

( ) ( ) ( ) , with = ξ ψ ξ ψ ξ ψ ∂ = + + ∂ G U G G B U U B U  (3.7)
It also proves convenient to integrate eqn. (3.7) in the ξ -direction to write the numerical y-flux as ( ) ). We want to make sure that eqns. (

1/2 1/2 , c numerical c d ψ ψ ξ ψ ψ ψ ξ ψ -   = = - ∆ = -  ∆   ∫ G G G B U  ( 3
( ) ( ) { } ( ) ( ) { } ( ) ( ) ( ) ( ) ( ) , , 1 1 ,
, 1 1 2 , 0 c c c c φ ξ ψ ψ ψ ψ φ ξ ψ ξ ξ ξ ξ ξ ψ ψ φ ξ ψ φ ξ ψ ξ ξ ξ ψ ψ ψ φ ξ ψ ξ ξ ψ ψ     ∂ - + ∆ ∂ - + ∆     + ∆ ∂ ∆ ∂ ∂ ∂     - - + ∆ - - + ∆ + =     ∆ ∂ ∆ ∂ G U F U F U G U U          (3.9)
The test functions are chosen from the same set of functions as the trial functions in eqn. (3.4).

From eqn. (3.4) it is easy to see that our trial functions are ( )

, 1 φ ξ ψ = , ( ) , φ ξ ψ ξ = and ( ) , φ ξ ψ ψ = .
Using the test function ( ) , 1 φ ξ ψ = and integrating over the entire wave model gives

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
( )

1/2 1/2 1/2 1/2 1/2 1/2 1/2 1/2 1 1 1/ 2, 1/ 2, 1/ 2, 1/ 2, 1 2 1 1 ,1/ 2 ,1/ 2 , 1/ 2 , 1/ 2 R L U D S d S d S d S d ψ ψ ψ ψ ψ ψ ξ ξ ξ ξ ξ ξ ξ ξ ψ ψ - - - -   - - - - -   ∆ ∆   = -    + - - - - - ∆ ∆     ∫ ∫ ∫ ∫ F U F U U G U G U (3.10)
In practice, one always obtains U (the mean value of the strongly interacting state) as early as possible in the calculation, because its value plays an important role in subsequent equations.

This value of U is used in eqn. (3.4) for the mean value and also in eqns. (3.5) and (3.7) to construct the characteristic matrices. It is also easy to show that when the flow is aligned with the

x-axis we have

D U S S = -, ( ) ( ) ,1/ 2 , 1/ 2 ξ ξ = - G G and
( ) ( ) 1/2 1/2 1/2 1/2 ,1/ 2 , 1/ 2 HLL d d ξ ξ ξ ξ - - = - = ∫ ∫ U U U
. The upshot is that for mesh-aligned flow, HLL = U U . In other words, when the flow is mesh-aligned, the mean value of the strongly- interacting state in the multidimensional Riemann solver matches with the corresponding state from the one-dimensional HLL Riemann solver, see eqn.

( )

1/2 1/2 1/2 1/2 1/2 1/2 1/2 1/2 1 1 1/ 2, 1/ 2, 1/ 2, 1/ 2, 2 2 1 1 ,1/ 2 ,1/ 2 , 1/ 2 , 1/ 2 4 R L c U D S d S d S d S d ξ ψ ψ ψ ψ ψ ψ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ξ ψ ψ ξ - - - -   - + - - -   ∆ ∆   = + ∆     + - - - - - ∆ ∆     ∆ + ∫ ∫ ∫ ∫ F U F U F U G U G U U (3.11)
For the first moment in the y-direction we use the test function ( ) , φ ξ ψ ψ = and integrate over the entire wave model to get

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
( )

1/2 1/2 1/2 1/2 1/2 1/2 1/2 1/2 1 1 1/ 2, 1/ 2, 1/ 2, 1/ 2, + 1 1 ,1/ 2 ,1/ 2 , 1/ 2 , 1/ 2 2 2 4 R L c U D S d S d S d S d ψ ψ ψ ψ ψ ψ ψ ψ ψ ξ ξ ψ ψ ξ ξ ξ ξ ξ ξ ψ ψ ψ - - - -   - - - - -   ∆ ∆   = ∆     + - + - - - ∆ ∆     ∆ + ∫ ∫ ∫ ∫ F U F U G U G U G U U (3.12)
This completes our description of the moments that are taken over the entire wave model,

[ ] [ ] 1/ 2,1/ 2 1/ 2,1/ 2 - × -
. The above three equations were already derived in Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF]. They are, however, used very differently in this paper to derive a MuSIC Riemann solver that is a close analogue of the one-dimensional HLLI Riemann solver. In principle, any one-dimensional Riemann solver can be used as a building block for the multidimensional Riemann solver, as shown in Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF]. However, to make the connection with the HLLI Riemann solver as tight as possible, we want the present multidimensional Riemann solver to be based on the same philosophy that was used for the one-dimensional HLLI Riemann solver in the limit where the flow is mesh-aligned.

Let us first establish a notational similarity between the multidimensional eigenstructure in this section and the one-dimensional eigenstructure from the previous section. We would like to obtain the best possible representation of the linear profile within the strongly interacting region. Let ( ) ξ ∆ U and ( ) ψ ∆ U denote undivided differences. Let us denote the linear profile in multiple dimensions as follows:

( ) ( ) ( ) , unprojected ξ ψ ξ ψ ξ ψ = + ∆ + ∆ U U U U  (3.13)
Typically, we wish to identify these undivided differences from the multidimensional wave model by looking at the solutions from the one-dimensional Riemann problems in the boundary of the multidimensional wave model. Thus we can write ( )

( ) ( ) 1/2 1/2 1/2 1/2 1/ 2, 1/ 2, d d ξ ψ ψ ψ ψ - - ∆ = - - ∫ ∫ U U U (3.14) and ( ) ( ) ( ) 1/2 1/2 1/2 1/2 ,1/ 2 , 1/ 2 d d ψ ξ ξ ξ ξ - - ∆ = - - ∫ ∫ U U U ( 3
( ) ( ) ( ) int 2 2 x x x x x x i i i i I l r ξ ξ ξ δ ∈   = ∆ = ∆   ∑ U U R δ L U  (3.16)
We can also relate ψ U to ( )

ψ ∆ U as ( ) ( ) ( ) int 2 2 y y y y y y i i i i I l r ψ ψ ψ δ ∈   = ∆ = ∆   ∑ U U R δ L U  (3.17)
Notice that we have evaluated the eigenstructure in both the x-and y-directions. As a result, x R and y R are matrices of right eigenvectors with dimension ) and (3.17). These factors can be greater than unity or they can even become less than unity. The amount of additional weight imparted by these factors is designed to ensure that the multidimensional Riemann solver retains favorable properties, as discussed in an ensuing paragraph.

We now ask the important question, which fluxes and states should we use in the integrals in eqns. (3.10), (3.11) and (3.12)? Our first instinct would be to use the linear profiles from eqn. (2.12). In fact, it can be shown that with that linear profile, and the definition for i δ given in eqn. (2.11), the x-flux in eqn. (3.6) will indeed reduce to the x-flux from the onedimensional HLLI Riemann solver when the flow is aligned with the x-axis. While this is proved in Appendix B, the proof steers us false! The fallacy is not in the math in Appendix B; in fact the mathematics is correct. The source of the fallacy is this:-If the logic of that mathematics is followed, it will lead us to a multidimensional Riemann solver that has some very poor entropy generation properties, especially in the vicinity of strong shocks! The source of the fallacy resides in the fact that we wanted the profiles ( )

,1/ 2 ξ U and ( ) , 1/ 2 ξ - U
to match the linear profiles from eqn. (2.12). However, realize that the one-dimensional HLLI Riemann solver produces overly steepened linear profiles. Such an over-steepened linear profile will produce lower than desired entropy in the transverse fluxes. In other words, the Lagrangian fluxes

( ) ( ) ,1/ 2 ,1/ 2 U S ξ ξ - G U and
( ) ( ) , 1/ 2 , 1/ 2 D S ξ ξ - - - G U
will produce less entropy than desired. When strong non-linearities are present in the flow, the resulting multidimensional Riemann solver will be unstable.

Having gained that insight, we draw upon our first goal. The goal is that for problems with strong discontinuities the expressions for numerical F and numerical G generate sufficient entropy to stabilize the problem. In the limit of strong discontinuities, the substructure, represented by ξ U and ψ U is irrelevant and can even be suppressed with the help of a switch that detects the presence of strong shocks. We therefore require numerical We now draw upon our second goal. When the flow is mesh-aligned, we want the expressions to reduce to their one-dimensional forms from Section II. In other words, when the flow is aligned with the x-axis, we want the expression from eqn. (3.6) to have dissipation characteristics that are similar to the one-dimensional HLLI Riemann solver from Section II. As in Section II, this will enable us to put bounds on the slope ξ U . For x-directional flow, we have

RU RD R = = U U U and LU LD L = = U U U . Eqn. (3.11) then give us ( ) ( ) 1 2 4 4 c HLL R R R L L L HLL S S ξ ξ ξ ξ ξ ∆ ∆ = + - + - + = +     F U F U F U U F U ( 3
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 2 1 2 2 2 2 numerical R L R L R L R L x x x x x x R L R L R L R L R L

S S S S S S S S S S S S S S

= +

  + - + - - - + -   - - -   F F F R Λ I δ Λ δ L U U (3.
x x i R L i R L x x x i x R i L i R L i x R L i R L x R L i R L S S S S S S S S S S S S S S S S φ λ δ λ λ φ λ φ - +    - - + ≤       =  + -           - - +      - - (3.20)
An analogous exercise for the y-flux, which is not repeated here for the sake of brevity, gives us (3.12), U , F and G are also fully specified. Eqns. (3.6) and (3.8) can, therefore, be used to obtain the numerical fluxes from the multidimensional Riemann solver. Also notice that we have already evaluated all or part of the eigenstructure so that we make the simplification I.e., the characteristic matrices are built from the sub-portion of the eigenstructure that has been evaluated.

( ) ( ) ( ) ( ) ( ) 2
- +    - - + ≤       =  + -           - - +      - - (3
The one-dimensional HLLI Riemann solver from Section II also preserves stationary contact discontinuities that are mesh-aligned. Consequently, the multidimensional Riemann solver described in this Section will do the same. As a result, the multidimensional Riemann solver described in this Section will inherit all the good well-balancing properties that were described for the one-dimensional HLLI Riemann solver in Section II.

A further observation is warranted. Recall that the two-dimensional HLLC Riemann solver from Balsara [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF] reduces to the one-dimensional HLLC Riemann solver for mesh-aligned flows on a two-dimensional Cartesian mesh. But this reduction becomes exact only if the onedimensional Riemann solvers on the boundary of the multidimensional wave model are indeed HLLC. If another type of one-dimensional Riemann solver is used, the result may still be pretty good. However, the reduction in the one-dimensional limit is not provably exact. A similar situation prevails for the discussion in this Section. The multidimensional MuSIC Riemann solver becomes closely analogous to the one-dimensional HLLI Riemann solver described in Section II for mesh-aligned flows on a two-dimensional Cartesian mesh. However, for this reduction to become philosophically consistent, the one-dimensional Riemann solvers on the boundary of the multidimensional wave model should indeed be based on the one-dimensional HLL Riemann solver. Appendix C of Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF] provides a computer-implementable strategy for reducing the variation in the mean state. The reduction of ξ ∆ U and ψ ∆ U will naturally result in a reduction in ξ U and ψ U via eqns. (3.16) and (3.17). In the limit of 0

ξ ψ = = U U
the present Riemann solver reduces exactly to an HLL Riemann solver, guaranteeing positivity. Also please note that when the variation in the flow is not one-dimensional and mesh-aligned, the fluxes F and G in eqns. (3.6) and (3.8) will indeed have multidimensional contributions. This is true even in the limit where 0 ξ ψ = = U U . I.e., the MuSIC Riemann solver described here is genuinely multidimensional even when sub-structure is absent.

It is also worth pointing out that in this narrative we have described the simplest multidimensional Riemann solver with sub-structure. In other words, one-dimensional HLL Riemann solvers provide the boundary information for the multidimensional wave model. There also exist one-dimensional HLLC (Toro, Spruce and Speares [START_REF] Toro | Restoration of contact surface in the HLL Riemann solver[END_REF], Gurski [START_REF] Gurski | An HLLC-type approximate Riemann solver for ideal magnetohydrodynamics[END_REF], Li [START_REF] Li | An HLLC Riemann solver for magnetohydrodynamics[END_REF]) and HLLD (Miyoshi and Kusano [56]) Riemann solvers which provide multiple states. The transverse fluxes across each of those states can be evaluated. This permits us to provide a more refined linear profile using a least squares minimization procedure (Vides et al. [START_REF] Vides | A simple two-dimensional extension of the HLL Riemann solver for hyperbolic conservation laws[END_REF], Balsara et al. [START_REF] Balsara | A Two-Dimensional Riemann Solver with Self-Similar Sub-Structure -Alternative Formulation Based on Least Squares Projection[END_REF]). Such a least squares procedure would consist of minimizing the linear profile from eqn. (3.13) along the boundary of the wave model. This is done by minimizing

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1/2 1/2 2 2 1/2 1/2 1/2 1/2 2 2 1/2 1/2 1/ 2, 1/ 2, 1/ 2, 1/ 2, 1 2 ,1/ 2 ,1/ 2 , 1/ 2 , 1/ 2 unprojected unprojected unprojected unprojected d d d d ψ ψ ψ ψ ψ ψ ξ ξ ξ ξ ξ ξ - - - -   - + - --         + - + - - -     ∫ ∫ ∫ ∫ U U U U U U U U     (3.23)
When a one-dimensional HLL Riemann solver is used, eqns. (3.14) and (3.15) provide the better strategy because they ensure that for mesh-aligned flow profiles we retrieve the one-dimensional Riemann solver. But it is useful to note that the above equation works well too when a more refined one-dimensional Riemann solver is used. Experience has shown that the major gain comes from having some reasonable sub-structure in the multidimensional Riemann solver. Introducing a more refined sub-structure using eqn. (3.23) does not provide that much incremental benefit.

III.b) Implementation-Related Details for MuSIC Riemann Solver on Cartesian Meshes:-

Before the multidimensional Riemann problem is evaluated, it is assumed that a onedimensional HLL Riemann solver has been applied immediately outside the boundary of the multidimensional wave model. This is inexpensive and easy to do. The steps in the implementation of the MuSIC Riemann solver with multidimensional sub-structure on structured meshes are as follows:

1) Construct the multidimensional wave model. See 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 construction of the multidimensional wave model for Cartesian meshes.

2) Examine the multidimensional wave model for the existence of a supersonic solution. The supersonic state and its corresponding numerical fluxes are obtained exactly as described in Fig. 5 from Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF]. Please also see the text associated with Fig. 5 of Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF]. If a supersonic case is discovered, the subsequent steps, all of which pertain to the subsonic case, are not carried out.

3) Obtain U from eqn. (3.10). The integrals on the right hand side of eqn. (3.10) are performed on the solution vectors and transverse fluxes that are obtained from the one-dimensional HLL Riemann solvers that abut the multidimensional wave model. (3.15). Please notice that these "undivided differences" should not be confused with ξ U and ψ U , which will be evaluated in the next two steps.

7)

In practice, it is also advisable to apply the multidimensional limiter from Appendix C of Balsara [START_REF] Balsara | Multidimensional Riemann Problem with Self-Similar Internal Structure -Part I -Application to Hyperbolic Conservation Laws on Structured Meshes[END_REF] to the linear profile in eqn. (3.13). Consequently, ξ ∆ U and ψ ∆ U undergo multidimensional limiting.

IV) Accuracy Analysis

It is desirable to demonstrate that our MuSIC Riemann solver, when coupled with high order spatial reconstruction and temporal evolution, produces suitably high order accuracy. To achieve high spatial accuracy we use WENO reconstruction (Jiang and Shu [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF], Balsara and Shu [START_REF] Balsara | Monotonicity preserving weighted non-oscillatory schemes with increasingly high order of accuracy[END_REF], Dumbser and Käser [START_REF] Dumbser | Arbitary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems[END_REF], Balsara et al. [START_REF] Balsara | Efficient, high-accuracy ADER-WENO schemes for hydrodynamics and divergence-free magnetohydrodynamics[END_REF], Balsara [START_REF] Balsara | Self-Adjusting, Positivity Preserving High Order Schemes for Hydrodynamics and Magnetohydrodynamics[END_REF], Balsara, Garain and Shu [START_REF] Balsara | An efficient class of WENO schemes with adaptive order[END_REF]) for the spatial reconstruction of conserved variables. The spatial reconstruction of divergencefree magnetic fields uses an extension of the WENO reconstruction strategy (Balsara [START_REF] Balsara | Divergence-free adaptive mesh refinement for magnetohydrodynamics[END_REF], [START_REF] Balsara | Second-order-accurate schemes for magnetohydrodynamics with divergencefree reconstruction[END_REF], [START_REF] Balsara | Divergence-free reconstruction of magnetic fields and WENO schemes for magnetohydrodynamics[END_REF], [START_REF] Balsara | Self-Adjusting, Positivity Preserving High Order Schemes for Hydrodynamics and Magnetohydrodynamics[END_REF], Balsara and Dumbser [START_REF] Balsara | Divergence-Free MHD on Unstructured Meshes Using High Order Finite Volume Schemes Based on Multidimensional Riemann Solvers[END_REF]). We also use a high order one-step ADER (Arbitrary DERivatives in space and time) time-evolution (Titarev & Toro [62], [START_REF] Titarev | ADER schemes for three-dimensional nonlinear hyperbolic systems[END_REF] and Toro & Titarev [START_REF] Toro | Solution of the generalized Riemann problem for advection reaction equations[END_REF], Dumbser et al. [START_REF] Dumbser | A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes[END_REF], Balsara et al. [START_REF] Balsara | Efficient, high-accuracy ADER-WENO schemes for hydrodynamics and divergence-free magnetohydrodynamics[END_REF], Balsara et al. [START_REF] Balsara | Efficient Implementation of ADER schemes for Euler and Magnetohydrodynamic flow on structured meshes -Comparison with Runge-Kutta methods[END_REF], Balsara and Kim [START_REF] Balsara | A subluminal relativistic magnetohydrodynamics scheme with ADER-WENO predictor and multidimensional Riemann solver-based corrector[END_REF]). Specifically, the ADER time-stepping that we use in this work on structured meshes is closest in style to Balsara et al. [START_REF] Balsara | Efficient, high-accuracy ADER-WENO schemes for hydrodynamics and divergence-free magnetohydrodynamics[END_REF] and Balsara et al. [START_REF] Balsara | Efficient Implementation of ADER schemes for Euler and Magnetohydrodynamic flow on structured meshes -Comparison with Runge-Kutta methods[END_REF].

For the present accuracy analysis, the multidimensional Riemann solver used onedimensional HLL Riemann solvers in the side panels. A linear variation of the variables, along with the limiting from eqns (3.20) and (3.21), was retained in the strongly interacting state. The use of eqns (3.20) and (3.21) contributes significantly to the stability of the multidimensional Riemann solver with sub-structure while simultaneously enabling the accurate treatment of slowly moving discontinuities. The multidimensional Riemann solver was integrated into the code using a simplified Simpson rule so that fourth order fluxes were evaluated at the zone faces. The RIEMANN code is three-dimensional so a weight of 1/6 was assigned to each of the four multidimensional Riemann solvers at the edge-centers associated with each face. A fifth onedimensional Riemann solve was done at the zone center using the one-dimensional Riemann solver described in Section II, and it was assigned a weight of 1/3. See eqn. (25) of Balsara [START_REF] Balsara | Multidimensional HLLE Riemann solver; Application to Euler and Magnetohydrodynamic Flows[END_REF] for further detail. Time-averaged states from the ADER predictor step were used as inputs to the Riemann solvers. The implicit assumption in such an approach is that spatial and temporal accuracy can commute.

A CFL number of 0.8 was used for all second and third order hydrodynamical simulations in this section. Up to third order it seems acceptable to assume that the spatial and temporal accuracy can indeed commute so that the advantages of a large timestep are retained. At fourth order, this assumption clearly breaks down, limiting the CFL to 0.4 if the order property is to be preserved. A recently completed von Neumann stability analysis for the numerical solution of the induction equation is now available, Balsara and Käppeli [START_REF] Balsara | von Neumann Stability Analysis of Globally Divergence-free RKDG and PNPM schemes for the Induction Equation Using Multidimensional Riemann Solvers[END_REF]. It shows that finite volume schemes for the induction equation are only stable up to a CFL of 0.666. As a result, all second and third order MHD and RMHD simulations in this paper used a CFL of 0.65. All fourth order MHD and RMHD simulations in this paper used a CFL of 0.4. For the sake of completeness, we also mention that in Balsara, Dumbser & Abgrall [START_REF] Balsara | Multidimensional HLL and HLLC Riemann Solvers for Unstructured Meshes -With Application to Euler and MHD Flows[END_REF] a Gaussian quadrature in time was used along with an ADER-WENO formulation. With that choice, it was shown that larger CFL numbers can be retained at all orders. Because that choice requires extra programming, we made a simpler choice in this paper. It must be appreciated that if the plan in this paragraph is implemented efficiently, each multidimensional Riemann solver contributes to the formation of four facial fluxes. Furthermore, the number of calls to the Riemann solvers does not increase with increasing order of accuracy, which would be the case if a Runge-Kutta timestepping were to be used. For MHD calculations, the Riemann solver acts at the edges of the mesh, therefore, the edge-aligned component of the electric field is directly available at each edge of the mesh.

In this section we demonstrate the accuracy of our methods by using hydrodynamical and MHD vortex flows. Analogous work for relativistic hydrodynamics and relativistic magnetohydrodynamics has also been done. In Balsara and Kim [START_REF] Balsara | A subluminal relativistic magnetohydrodynamics scheme with ADER-WENO predictor and multidimensional Riemann solver-based corrector[END_REF] we have presented relativistic analogues of the hydrodynamical and MHD vortices that are analyzed here. In that paper we also show that those relativistic vortex problems can be used to demonstrate the higher order accuracy of relativistic hydrodynamics and RMHD codes.

IV.a) Accuracy Analysis for Hydrodynamical Vortex Problem

This hydrodynamical vortex problem was presented in Jiang and Shu [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF]. In this problem, an isentropic vortex propagates at 45° to the grid lines in a domain with periodic boundaries given by [-5, 5] x [-5, 5]. Since the vortex stably preserves its form, it returns to its original location. Since the problem is well-known in the literature, we do not describe it in detail here. The velocity and temperature fluctuations for the circular vortex have an exponential falloff with radius, which ensures that the fluctuations are quite close to zero at the domain boundaries. However, for the fourth order scheme the domain is increased to [-10, 10] x [-10, 10] due to the fact that the nonzero values of the exponential function at the boundaries are picked up by the fourth order scheme on the smaller domain. The stopping time was set to 10 time units for the second and third order schemes and to 20 time units for the fourth order scheme because of the bigger domain. The stopping time was chosen so that the vortex has completed one periodic passage through the computational domain. We report on the accuracy of the vortex after it has completed one orbit in the computational domain. Table I shows the accuracies in the L1 and L∞ norms for the density variable for second, third and fourth order schemes at the stated stopping times. All three ADER-WENO methods with the MuSIC Riemann solver permit the use of a large CFL number. They all meet the expected order of accuracy even for a small number of zones, showing that they have uniform convergence properties.

TABLE I shows the accuracy analysis for the hydrodynamical vortex problem as measured

in the density variable.

Method # of zones 1 L Error 1 L Order L ∞ Error L ∞ Order ADER-WENO 2 nd

IV.b) Accuracy Analysis for MHD Vortex Problem

The magnetized isodensity vortex problem described in Balsara [START_REF] Balsara | Second-order-accurate schemes for magnetohydrodynamics with divergencefree reconstruction[END_REF] consists of a magnetized vortex moving across a domain given by [-5, 5] x [-5, 5] at an angle of 45° for a time of 10 units. As before, for the fourth order scheme the domain is increased to [-10, 10] x [-10, 10] and the simulation time is increased to 20 units. Periodic boundaries are used for the domain. Since the problem is well-known in the literature, we do not describe it in detail here. We report on the accuracy of the x-component of the magnetic field of the vortex after it has completed one orbit in the computational domain. Table II shows the accuracies in the L1 and L∞ norms for the x-magnetic field variable for second, third and fourth order schemes at the stated stopping times. All three ADER-WENO methods with the MuSIC Riemann solver permit the use of a large CFL number. They all meet the expected order of accuracy even for a small number of zones, showing that they have uniform convergence properties.

V) Test Problems

A broad range of test problems is presented here. The problems span many different PDE systems showing that the 2D MuSIC Riemann solver, along with the 1D HLLI Riemann solver, works very well. We first show a few results involving Euler flow. Then we show results from classical MHD. Subsequently, we show results involving RMHD flows.

V.a) Euler Flow: Forward Facing Step

This problem was first presented in Woodward and Colella [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF] and Colella and Woodward [START_REF] Colella | The piecewise parabolic method (PPM) for gas-dynamical simulations[END_REF]. It is so well-known in the computational hydrodynamics community that we do not describe its set-up here. The problem was run on a 960×320 zone mesh spanning the domain [ ] [ ] 0,3 0,1 × to a final time of 4.0. A fourth order ADER-WENO scheme with the MuSIC Riemann solver was used. Fig. 3 shows the density variable from the forward facing step problem at a time of 4.0. The problem consists of a Mach 3 wind tunnel with a forward-facing step. As the bow shock reflects off the step and then the top wall of the wind tunnel, it establishes a triple-point structure. A vortex sheet emanates from the triple point. The vortex sheet that emanates from this triple-point structure shows an interesting roll-up due to Kelvin-Helmholtz instability. We see that our simulation has captured the roll-up of the vortex sheet very well.

V.b) Euler Flow: Double Mach Reflection

This problem was first presented in Woodward and Colella [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF] and Colella and Woodward [START_REF] Colella | The piecewise parabolic method (PPM) for gas-dynamical simulations[END_REF]. It is also well-known in the computational hydrodynamics community, with the result that we do not describe its set-up here. The problem consists of a strong Mach 10 shock that is initially incident on a reflecting wedge that makes an angle of 60 o with the plane of the shock. The problem was run on a 1920×480 zone mesh spanning the domain [ ] [ ] 0, 4 0,1 × to a final time of 0.2. A fourth order ADER-WENO scheme with the MuSIC Riemann solver was used. Fig. 4a shows the density variable from the double Mach reflection problem at a time of 0.2. Fig. 4b shows a zoom-in of the roll-up of the Mach stem. We see that we obtain a very wellresolved roll-up of the Mach stem.

V.c) Classical MHD Flow: Rotor Problem

This well-known MHD problem was first documented in Balsara & Spicer [START_REF] Balsara | A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations[END_REF] and also Balsara [START_REF] Balsara | Second-order-accurate schemes for magnetohydrodynamics with divergencefree reconstruction[END_REF]. It consists of a central, uniformly rotating vortex in a non-rotating ambient medium. An initially uniform magnetic field threads through both regions. The details are described in the above-mentioned references. The problem initially has a uniform magnetic field threading the computational domain [

] [ ] 0.5, 0.5 0.5, 0.5 -× -. Within a radius of 0.1 units, we initialize a dense rotating fluid, which forms the rotor. The motion of the rotor initiates outwar-propagating waves. The problem was run on a 1000×1000 zone mesh to a final time of 0.29. A fourth order ADER-WENO scheme with the MuSIC Riemann solver was used. Fig. 5 shows the results from the MHD Rotor test problem. Figs. 5a, 5b, 5c and 5d show the density, pressure, magnitude of the fluid velocity and the magnitude of the magnetic field at the final time. All the requisite MHD flow features are captured nicely in our simulations.

V.d) Classical MHD Flow: 3D Blast Problem

The present test problem is a 3D variant of the well-known 2D MHD blast wave problem that was first documented in Balsara & Spicer [START_REF] Balsara | A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations[END_REF]. The plasma β measures the ratio of the thermal pressure to the magnetic pressure. As the plasma's β becomes smaller, this problem becomes increasingly stringent. The problem consists of a 1.4 γ = gas with unit density and a pressure of 0.1 initialized on a 257 3 zone mesh spanning the unit cube. Initially we have B B B 150 3

x y z = = = .
The pressure is initially reset to a value of 1000 inside a central region with a radius of 0.1. The plasma's β is initially given by 1.117×10 -4 . A CFL number of 0.4 was used. The problem is run up to a time of 0.0075, by which time a strong magnetosonic blast wave propagates through the domain. The problem was run with a third order ADER-WENO scheme with the MuSIC Riemann solver applied at the edges of the mesh. Fig. 6 shows the variables from the 3D blast problem in the z = 0 mid-plane of the computational domain. Fig. 6a shows the plot of the density for the mid-plane in the z-direction. Fig. 6b shows the same for the pressure in the same plane. Figs. 6c and6d show the magnitude of the velocity and the magnitude of the magnetic field, again in the same plane. We see that the densities and pressures are positive as expected.

V.e) Classical MHD Flow : Decay of Finite Amplitude Torsional Alfven Waves

Turbulence studies play an increasingly important role in several fields, like astrophysics or space physics. (The Orzag Tang [START_REF] Orszag | Small-scale structure of two-dimensional magnetohydrodynamic turbulence[END_REF] problem is a first step in the development of turbulence.) The ability to propagate finite amplitude Alfven waves over large distances and long times on a computational mesh is crucial for carrying out simulations of MHD turbulence. If the Alfven waves are damped strongly because of inherent numerical dissipation in a code, the code will fail to capture the resulting turbulence. This is because MHD turbulence is mainly sustained by Alfven waves. The Alfven wave decay test problem, first presented by Balsara [START_REF] Balsara | Second-order-accurate schemes for magnetohydrodynamics with divergencefree reconstruction[END_REF], examines the numerical dissipation of torsional Alfven waves in two dimensions. In this test problem torsional Alfven waves propagate at an angle of 9.462 o to the y-axis through a domain given by [-3, 3] x [-3, 3] . The domain was set up with 120 x 120 zones and has periodic boundary conditions. We do not present further details of the set-up, because the problem is already well-described in the above-mentioned paper. The simulation was stopped at 129 time units by which time the Alfven waves had crossed the domain several times. Depending on the dissipation properties of the scheme, the amplitude of the torsional Alfven wave will, of course, decay. A more dissipative method will cause greater dissipation of the Alfven wave; a less dissipative method will reduce that dissipation.

It is often said that the quality of the Riemann solver is not very important, especially when high order schemes are used. But practitioners have not quantified the precise order of accuracy of the scheme at which the quality of the Riemann solver becomes immaterial. We set out to quantify this order of accuracy for MHD simulations. To that end, we simulated the torsional Alfven wave decay problem with second, third and fourth order schemes with the 1D HLLI Riemann solver along with the 2D MuSIC Riemann solver with sub-structure. Used in this fashion, both the 1D and 2D Riemann solvers are complete; i.e. they fully represent all the waves that arise in the MHD system. We then simulated the same problem again with the same second, third and fourth order schemes. However, this time we used a 1D HLL Riemann solver along with the 2D MuSIC Riemann solver without any sub-structure. In other words, in our second set of simulations both Riemann solvers did not resolve any intermediate waves.

Figs. 7a and7b show the evolution of the maximum z-velocity and maximum zcomponent of the magnetic field in the torsional Alfven wave as a function of time. For the simulations shown in Figs. 7a and 7b we used the 1D HLLI Riemann solver along with the 2D MuSIC Riemann solver with sub-structure. Figs. 7c and7d show the same information as Figs. 7a and 7b, the only difference being that we used the 1D HLL Riemann solver along with the 2D MuSIC Riemann solver without sub-structure. Comparing the two sets of figures, we see that the inferior Riemann solvers produce a six-times larger decay in the amplitude of the Alfven wave at second order. At third order, the inferior Riemann solvers produce a three-times larger decay in the amplitude of the Alfven wave. Notice that the second order scheme with superior Riemann solvers is less dissipative than the third order scheme with inferior Riemann solvers! At fourth order, the difference between the inferior Riemann solvers and the exact Riemann solvers is almost negligible. We, therefore, conclude that second and third order schemes are greatly benefited by the quality of the Riemann solver. It is only at fourth and higher orders of accuracy that the difference between a superior and an inferior Riemann solver begins to become quite small! However, please note that a fourth order scheme has computational complexity that is substantially higher than a second or third order scheme. The Riemann solver with substructure has a computational complexity that is only marginally higher than a Riemann solver without substructure. As a result, it is very advantageous to improve the quality of all schemes at all orders.

V.f) RMHD Flow: Relativistic Orzag Tang Problem

The Orzag Tang test problem (Orzag Tang [START_REF] Orszag | Small-scale structure of two-dimensional magnetohydrodynamic turbulence[END_REF]) is designed to illustrate the transition to turbulence for MHD flows. The RMHD variant of that test problem has been proposed by Beckwith and Stone [START_REF] Beckwith | A Second-order Godunov Method for Multi-dimensional Relativistic Magnetohydrodynamics[END_REF]. We do not repeat the set-up here. The problem was set up on a unit square with 1000×1000 zones and run to a final time of 0.8. The problem was run with a fourth order ADER-WENO scheme with the MuSIC Riemann solver applied at the edges of the mesh. Figs. 8a, 8b, 8c and 8d show the density, pressure, magnitude of the velocity and magnitude of the magnetic field at the final time for the relativistic Orzag Tang problem. All the requisite MHD flow features are captured nicely in our simulations.

V.g) RMHD Flow: Decay of Finite Amplitude Torsional Alfven Waves

In Section V.e we studied the decay of torsional Alfven waves in classical MHD and showed that the importance of a high-quality Riemann solver is only diminished at fourth and higher orders. It is interesting to speculate whether such a result also applies to other hyperbolic systems. In Balsara and Kim [START_REF] Balsara | A subluminal relativistic magnetohydrodynamics scheme with ADER-WENO predictor and multidimensional Riemann solver-based corrector[END_REF] we constructed a relativistic analogue of the torsional Alfven wave decay problem from Balsara [START_REF] Balsara | Second-order-accurate schemes for magnetohydrodynamics with divergencefree reconstruction[END_REF]. We do not repeat the details of the set-up here.

Figs. 9a and9b show the evolution of the maximum z-velocity and maximum zcomponent of the magnetic field in the relativistic torsional Alfven wave as a function of time. For the simulations shown in Figs. 9a and 9b we used the 1D HLLI Riemann solver along with the 2D MuSIC Riemann solver with sub-structure. Figs. 9c and9d show the same information as Figs. 9a and 9b, the only difference being that we used the 1D HLL Riemann solver along with the 2D MuSIC Riemann solver without sub-structure. Comparing the two sets of figures, we see that the inferior Riemann solvers again show substantially larger dissipation at second and third orders. It is only at fourth order that we find a much-reduced difference between a Riemann solver with sub-structure and a Riemann solver that does not resolve any intermediate waves. We, therefore, conclude that a Riemann solver that resolves intermediate waves is very important for reducing dissipation in second and third order schemes. At fourth and higher orders, that importance is diminished. As mentioned in Sub-section V.e, the incremental costs of including sub-structure in a Riemann solver are only slight, making it advantageous to improve the quality of all schemes at all orders.

VI) Conclusions

In this paper we have shown the importance of introducing sub-structure in one-and multi-dimensional Riemann solvers. When all the one-dimensional waves are represented in a one-dimensional Riemann solver, the Riemann solver is said to be complete. In an analogous fashion we say that when all the sub-structure from all the waves that can propagate in multidimensions is correctly represented in a multi-dimensional Riemann solver, then it too is said to be complete. In this paper, we present MuSIC Riemann solvers that can retain sub-structure from all the waves that propagate in all directions of a multi-dimensional hyperbolic system. This paper represents the first time that the dissipation properties of multi-dimensional Riemann solvers has been studied in detail. All implementation-related details are presented in pointwise fashion to facilitate easy implementation.

Several stringent test problems drawn from Euler flow, non-relativistic MHD flow and RMHD flow are presented. The importance of retaining sub-structure in the Riemann solvers is conclusively demonstrated, especially at second and third orders. It is also shown that for fourth and higher orders, the importance of a Riemann solver that retains sub-structure is diminished. Consequently, for fourth and higher order Godunov schemes it might be possible to use less expensive Riemann solvers (ones which do not retain sub-structure) without experiencing a significant degradation in simulation quality. However, the computational complexity of a Riemann solver with substructure is only marginally higher than a Riemann solver without substructure. As a result, it is very advantageous to improve the quality of all schemes at all orders.

Helpful educational videos and codes that can facilitate the reader's understanding of multidimensional Riemann solvers and their use in higher order schemes are available on http://www.nd.edu/~dbalsara/Numerical-PDE-Course.

Appendix A

Here we provide the formulation of the one-dimensional HLLI Riemann solver for an arbitrary ALE mesh. It is very desirable that intermediate waves moving with the same speed as the mesh should be treated on the moving mesh with no dissipation and in this Appendix we provide details of such an HLLI Riemann solver.

( ) ( )( ) ( ) 2 1 = 2 v 2 2 3 m m c numerical HLL m R L ξ ξ ξ   ∆ - - + - -   ∆   F F R δ Λ I δ L U U
where the superscripts "m" denote a moving mesh. Here v

( ) ( ) ( ) ( ) and R R L L R L R L L R R L R L HLL HLL R L R L S S S S S S S S S S - - - - + - = = - - U U F F F F U U U F
Please also compare our present formula for the HLLI numerical flux to the one in eqn. (2.9c) to see that it has undergone only a slight change with ( )

v m → - Λ Λ I .
The only further thing that needs to be specified is the diagonal matrix δ . For making a practical computer implementation, it is useful to define three auxiliary variables

v R R m S S ≡ -  , v L L m S S ≡ -  and v i i m λ λ ≡ - 
. These three terms, with the tilde on top, can then be used in our regular HLLI Riemann solver to evaluate the weights i δ . We, therefore, have

( ) ( ) ( ) ( ) 2
i R L i R L i R i L i R L i R L i R L i S S S S S S S S S S S S φ λ δ λ λ φ λ φ - +    - - + ≤          =    + -       - - +                           ( ) 2 3 = R L R L S S S S - -    
Notice that the above formula does depend quite strongly on the speed of the mesh, v m . This completes our description of the one-dimensional HLLI Riemann solver for ALE meshes.

The previous equation exactly matches eqn. (2.7), showing that the constant parts are identical. Now compare eqn. (3.5) to eqn. (2.4). We see that in the limit where 0 ψ = U , the two expressions for the x-flux are identical. We, therefore, see that there is a one-to-one correspondence between the constant part of the x-flux as well as the ξ -dependent part of the x- flux. The two expressions for the x-flux can be made equivalent when x i δ is defined analogously to eqn.

(2.11) in the one-dimensional limit. We have, therefore, proved that the multidimensional Riemann solver reduces exactly to the one-dimensional Riemann solver on a Cartesian mesh when the flow has become mesh-aligned and one-dimensional. However, as argued in Section III, this proof steers us false.

Figure Captions

Fig. 1 compares the dissipation from the HLLI Riemann solver (dashed lines) to the least possible dissipation, as exemplified by the Roe-type Riemann solver (solid lines). Fig. 1a shows the dissipation from the HLLI Riemann solver and also the theoretically minimum dissipation as a function of the wave speed when SL = -0.9 and SR = 0.1. The dissipation from the HLLI Riemann solver is only 23.2% larger than the Roe Riemann solver. Fig. 1b shows similar information when SL = -0.7 and SR = 0.3. Fig. 1c shows similar information when SL = -0.5 and SR = 0.5. In Figs. 1b and 1c the dissipation properties coincide for both Riemann solvers. Fig. 1d shows analogous information when SL = -0.2 and SR = 0.8; indicating a difference of 17.6% in the dissipation. Fig. 1e shows the situation when SL = -0.01 and SR = 0.99; indicating practically no difference in the dissipation. Fig. 2 shows the multidimensional wave model for a Cartesian mesh. Here the thick solid line denotes the boundary of the multidimensional wave model; the interior of the wave model is shaded. The four initial states that come together at a vertex "O" of the mesh are also shown. The thin solid lines in Fig. 2 show the extremal speeds of the one-dimensional Riemann problems in the boundary of the multidimensional wave model. The dashed lines in Fig. 2 show the coordinate axes, measured as speeds. The bounding speeds of the multidimensional wave model are also shown. The black dot with (ξc,ψc) is the centroid of the wave model. Fig. 6 shows the variables from the 3D blast problem in the z = 0 mid-plane of the computational domain. Fig. 6a shows the plot of the density for the mid-plane in the z-direction. Fig. 6b shows the same for the pressure in the same plane. Figs. 6c and6d show the magnitude of the velocity and the magnitude of the magnetic field, again in the same plane. Fig. 6 shows the variables from the 3D blast problem in the z = 0 mid-plane of the computational domain. Fig. 6a shows the plot of the density for the mid-plane in the z-direction. Fig. 6b shows the same for the pressure in the same plane. Figs. 6c and6d show the magnitude of the velocity and the magnitude of the magnetic field, again in the same plane.

(

×

  the right eigenvectors being considered; L is a corresponding matrix of left eigenvectors with dimension and δ is a diagonal matrix of dimension ( ) ( ) denotes the number of elements in the set " int I ". We will specify the diagonal elements of δ shortly and we will see that each diagonal term i δ in the diagonal matrix δ depend on the structure of the wave model as well as the wave speed i λ . Therefore, in order to be consistent with the Galerkin projection, we should substitute the value of ξ U from eqn. (2.8) in eqn. (2.7) to get the flux F .

1 L 7 L 7 L 2 L

 1772 from Fig.1athat our HLLI-style Riemann solver always produces dissipation that is within 23.2% of the Roe-type Riemann solver. (Please also note that the analogous plot for 0.to Fig.1aafter it is flipped about the vertical axis given by 0 i λ = . This trend extends to all the other panels in Fig.1.) Fig.1bshows similar information when 0.. 1b we see that the dissipation of the HLLI-style Riemann solver coincides with the dissipation of the Roe-type Riemann solver when 0.us that the two Riemann solvers produce identical dissipation. Fig.1dshows similar information when 0.us that our HLLI Riemann solver always produces dissipation that is within 17.6% of the Roe-type Riemann solver. Fig.1eshows analogous information when 0

  the present HLLI Riemann solver for any value of the similarity variable x t ξ =  . The formulae in this paragraph are valid as long as ξ  lies in the range L R S S ξ < <  , i.e. within the Riemann fan. Using i δ from eqn. (2.11), we get

  (2.14) is left unchanged by the inclusion of sub-structure in the Riemann fan. The flux in eqn. (2.14) is based on considering fluctuations. Comparing the fluxes in eqns.



  correspond most naturally to ( ) , x y . We make a scaled and shifted coordinate transformation in the similarity variables with

(

. 8 )

 8 For eqns.(3.6) and (3.8) recall that the time axis corresponds to ( )( )

  when the flow is mesh-aligned, we want the expressions to become analogous to the one-dimensional forms from Section II. In other words, when the flow is aligned with the x-axis, we want the expression from eqn.(3.6) to have dissipation characteristics that are similar to the one-dimensional HLLI Riemann solver from Section II. As in Section II, this will enable us to put bounds on the slope ξ U . A similar consideration for flow that is aligned with the y-axis will enable us to put bounds on the slope ψ U .By multiplying eqn. (3.3) by a test function ( ) , φ ξ ψ , we can make it more ready for the Galerkin projection in similarity variables. Consequently, we get

  (2.6). Having obtained U with the help of zeroth moments, let us now consider the first moments of the governing equation. For the first moment in the x-direction we use the test function ( ) , φ ξ ψ ξ = and integrate over the entire wave model to get

. 15 )

 15 As in Section II, ( ) ξ ∆ U and ( ) ψ ∆ U can be thought of as the unprojected slopes. They are related to ξ U and ψ U respectively by appropriate projections that can be made with the left and right eigenvectors. The weights that are assigned to those projections are designed to bring out certain favorable properties in the multidimensional Riemann solver. To that end, we identify the interior waves in both directions for the state U . Let { } and right-and left-eigenvectors in the x-direction associated with the state U . and right-and lefteigenvectors in the y-direction associated with the state U . We assume that the eigenstates are so ordered that the same set int I labels the intermediate waves in either direction; this is usually possible for most hyperbolic systems. (For example, in MHD we could use the set int I to label a left-going Alfven wave, an entropy wave in the x-direction and a right-going Alfven wave. We can use the same set to label a downward-going Alfven wave, an entropy wave in the y-direction and an upward-going Alfven wave.) It is worth pointing out that since the x-and y-directional eigenvectors are built from the same state U , waves of a given wave family that are moving in any arbitrary direction can be projected in the linear space of the two sets of eigenvectors. We can now relate ξ U to ( ) ξ ∆ U in a fashion that is closely analogous to eqn. (2.8) as follows

  -and y-directions; and please note that the two matrices are not the same. Similar considerations hold for matrices of left eigenvectors, x L and y L , with dimension ( ) eigenvalues in the x-and y-directions are denoted byx Λ and y Λ respectively. The elements of the two diagonal matrices x δ and y we change the eigenvector projection in eqns.(3.16

  extent, the fluxes and states that we put into the integrals in eqns. (3.10),(3.11) and (3.12) are a matter of choice. We choose to use the piecewise-constant fluxes and states that come from the one-dimensional HLL Riemann solver. With that choice, numerical F and numerical G will indeed reduce to the multidimensional HLL values from Balsara [

δ

  fully specified by the above equations, we realize that eqns. (3.14) and (3.16) give us ξ U . Likewise, eqns. (3.15) and (3.17) give us ψ U . The integrals over the side panels of the multidimensional wave model in eqns. (3.10), (3.11) and (3.12) are fully specified by the one-dimensional HLL Riemann solvers in those side panels. From eqns. (3.10), (3.11) and

  Notice that when " int I " is a complete set of intermediate waves, the MuSIC Riemann solver is indeed complete in a multidimensional sense. In other words, it retains all the intermediate waves that can exist in the multidimensional wave model. Positivity can also be ensured by requiring that positivity is enforced at the vertices of our multidimensional wave model. If the initial profile in eqn. (3.13) indicates a loss of positivity at any of the vertices of the multidimensional wave model, we should feel free to reduce ξ ∆ U and ψ ∆ U in that equation.

8) Evaluate x i δ and y i δ

 i from eqns. (3.20) and (3.21). Then obtain ξ U and ψ U from eqns. (3.16) and (3.17). 9) Using U , ξ U and ψ U in eqns. (3.11) and (3.12), obtain F and G . The integrals on the right hand side of eqns. (3.11) and (3.12) are performed by using the solution vectors and transverse fluxes from the one-dimensional HLL Riemann solvers that abut the multidimensional wave model. 10) Now obtain the numerical fluxes numerical F and numerical G from eqns. (3.6) and (3.8).

  Let the zone boundary at which this Riemann solver is applied have a mesh motion given by v m . Here L S and R S still denote the speeds of the extremal left-going and right-going wave as they would have been evaluated on a fixed mesh. We have v Realize that the mesh motion does not change the left and right eigenvectors. However, the corresponding set of eigenvalues that are evaluated relative to the moving mesh becomes { } analogue of the numerical flux from eqn. (2.9c) for a moving mesh is given by

  evaluated on a moving mesh. The HLL state and flux HLL U and HLL F are defined simply relative to a stationary frame of reference in the usual way as

Fig. 3

 3 Fig.3shows the density variable from the forward facing step problem at a time of 4.0.

Fig. 4a shows

  Fig.4ashows the density variable from the double Mach reflection problem at a time of 0.2. Fig.4bshows a zoom-in of the roll-up of the Mach stem.

Fig. 5

 5 Fig. 5 shows the results from the MHD Rotor test problem. Figs. 5a, 5b, 5c and 5d show the density, pressure, magnitude of the fluid velocity and magnitude of the magnetic field at the final time.

Figs. 7a andFig. 1

 1 Figs.7a and 7bshow the evolution of the maximum z-velocity and maximum z-component of the magnetic field in the torsional Alfven wave as a function of time. For the simulations shown in Figs. 7a and 7b we used the 1D HLLI Riemann solver along with the 2D MuSIC Riemann solver with sub-structure. Figs.7c and 7dshow the same information as Figs.7a and 7b, the only difference being that we used the 1D HLL Riemann solver along with the 2D MuSIC Riemann solver without sub-structure.

Fig. 2 Fig. 3

 23 Fig.2shows the multidimensional wave model for a Cartesian mesh. Here the thick solid line denotes the boundary of the multidimensional wave model; the interior of the wave model is shaded. The four initial states that come together at a vertex "O" of the mesh are also shown. The thin solid lines in Fig.2show the extremal speeds of the one-dimensional Riemann problems in the boundary of the multidimensional wave model. The dashed lines in Fig.2show the coordinate axes, measured as speeds. The bounding speeds of the multidimensional wave model are also shown. The black dot with (ξ c ,ψ c ) is the centroid of the wave model.

FigsFigs

  Figs. 7a and 7b show the evolution of the maximum z-velocity and maximum z-component of the magnetic field in the torsionalAlfven wave as a function of time. For the simulations shown in Figs. 7a and 7b we used the 1D HLLI Riemann solver along with the 2D MuSIC Riemann solver with sub-structure. Figs.7c and 7dshow the same information as Figs.7a and 7b, the only difference being that we used the 1D HLL Riemann solver along with the 2D MuSIC Riemann solver without sub-structure.

  3.6) and (3.8) meet two important goals. First, for problems with strong discontinuities in arbitrary directions the expressions for numerical F and numerical G generate sufficient entropy to stabilize the problem. When strong discontinuities are present, the substructure, represented by ξ U and ψ U is irrelevant and can be zeroed out. This can be accomplished with the help of a sensor function that detects the presence of a strong discontinuity. We therefore require numerical

  .18) Compare eqn. (3.18) to eqn. (2.7) to notice that the two equations differ in detail. Consequently, putting eqn. (3.18) into eqn. (3.6) and simplifying gives us

4 )

 4 If the density in U is substantially lower than the minimum density in the states RU U , LU U , we do not provide linear sub-structure. Similarly, if the density in U is substantially greater than the maximum density in the states RU

									U , LU U , LD U and RD U we also
	do not provide linear sub-structure. Similar considerations are made for the pressure. I.e., this is
	just a reasonable and physical criterion for deciding whether it is justified to include sub-
	structure in the Riemann solver. If the decision is made to exclude sub-structure, the eigenvector
	evaluation and projection in the ensuing few steps is not needed.
	5) Using U , obtain the eigenstructure { λ ∈ : x i i I	int	}	, {	: r i I x i ∈	int	}	and { l i I : x i ∈	int	}	in the x-
	direction as well as the eigenstructure { λ ∈ : y i i I	int	}	, {	: r i I y i ∈	int	}	and { l i I : y i ∈	int	}	in the y-
	direction. Note that only the intermediate waves of interest are needed; and these waves are
	usually easier to evaluate than the entire eigenstructure.
	6) Obtain ξ ∆ U and

LD U and RD U ψ ∆ U from eqns. (3.14) and

TABLE II shows the accuracy analysis for the MHD vortex problem as measured in the x- component of the magnetic field.

 II

	Method	# of zones	1 L Error	1 L Order	L ∞ Error	L ∞ Order
	ADER-WENO
	2 nd Order,
	CFL 0.8
		64x64	3.1447E-3		2.5112E-2
		128x128	7.8041E-4	2.01	6.0608E-3	2.05
		256x256	1.9452E-4	2.00	1.5432E-3	1.97
		512x512	4.8703E-5	2.00	3.8818E-4	1.99
	ADER-WENO
	3 rd Order,
	CFL 0.8
		64x64	5.2781E-4		1.0549E-2
		128x128	6.3890E-5	3.05	1.6369E-3	2.69
		256x256	7.8487E-6	3.03	2.7012E-4	2.60
		512x512	1.0024E-6	2.97	5.0893E-5	2.41
	ADER-WENO
	4 th Order,
	CFL 0.4
		64x64	4.5538E-4		3.0456E-2
		128x128	2.3059E-5	4.30	1.4159E-3	4.43
		256x256	1.3014E-6	4.15	6.2495E-5	4.50
		512x512	7.5300E-8	4.11	3.6479E-6	4.10

Acknowledgements

DSB acknowledges support via NSF grants NSF-DMS-1361197, NSF-ACI-1533850, NSF-DMS-1622457. BN acknowledges the financial support from the National French Research Program (ANR): ANEMOS (2011), ANR-11-MONU-002. Several simulations were performed on a cluster at UND that is run by the Center for Research Computing. Computer support on NSF's XSEDE and Blue Waters computing resources is also acknowledged.

Appendix B

In this Appendix we prove that the multidimensional Riemann solver from Section III reduces exactly to the one-dimensional Riemann solver in Section II on a two-dimensional Cartesian mesh where all the variations in the input state are restricted to the x-direction. This happens when RU . We also have

In that limit, eqn. (3.11) for the constant part of the x-flux, F , reduces as follows:-

In the one-dimensional limit, we also get ( )