{'fraction_non_alphanumeric': 0.0437509527646508, 'fraction_numerical': 0.015854003789281125, 'mean_word_length': 4.295202952029521, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On a class of two-dimensional incomplete Riemann solvers

We propose a general class of genuinely two-dimensional incomplete Riemann solvers for systems of conservation laws. In particular, extensions of Balsara's multidimensional HLL scheme [J. Comput. Phys. 231 (2012) 7476-7503] to two-dimensional PVM/RVM (Polynomial/Rational Viscosity Matrix) finite volume methods are considered. The numerical flux is constructed by assembling, at each edge of the computational mesh, a one-dimensional PVM/RVM flux with two purely two-dimensional PVM/RVM fluxes at vertices, which take into account transversal features of the flow. The proposed methods are applicable to general hyperbolic systems, although in this paper we focus on applications to magnetohydrodynamics. In particular, we propose an efficient technique for divergence cleaning of the magnetic field that provides good results when combined with our two-dimensional solvers. Several numerical tests including genuinely two-dimensional effects are presented to test the performances of the proposed schemes.


Introduction
Since the pioneering work of Godunov [14], Riemann solvers have been an important ingredient in the design of robust and accurate numerical methods for hyperbolic conservation laws. Usually, the exact solution of a Riemann problem contains many complex features which makes it unsuitable in practice. For this reason, a number of incomplete Riemann solvers have been devised in the literature (Lax-Friedrichs, Rusanov, HLL, FORCE, etc.; see [7] and the references therein), which takes into account only some of the waves appearing in the Riemann fan.
In the paper [7] the authors introduced the so-called PVM methods, which constitute a general class of incomplete Riemann solvers defined in terms of viscosity matrices based on the polynomial evaluation of a given Roe matrix or the Jacobian of the flux at some average state. Indeed, they showed that many classical schemes in the literature can be viewed as particular cases of PVM solvers. The idea behind PVM methods was later extended in [8] to the case of rational functions, giving rise to the class of RVM methods, which provide a much more precise representation of internal waves. In what follows, we will use the term AVM (Approximate Viscosity Matrix) to refer to both PVM and RVM solvers. In [8] it was also shown that the choice of a precise first-order solver is important when designing higher order methods in terms of computational efficiency, at least for solutions involving complex structures. In particular, the choice of appropriate underlying functions for an AVM solver allows to control the amount of numerical difussion of the resulting scheme. An important feature of AVM methods is that they do not need the whole spectral decomposition of the system, but only a bound on its spectral radius. This makes this kind of solvers particularly efficient for treating complex problems: see [7,8] for applications to magnetohydrodynamics (MHD) and multilayer shallow water systems. Further applications to relativistic MHD where considered in [10], where a Jacobian-free version of some PVM solvers was also introduced. Furthermore, in [9] the idea of AVM solvers was applied to build approximate versions of the DOT method [12], which in turn can be viewed as approximations of the classical Osher-Solomon scheme.
When solving multidimensional problems, it is a common practice to consider one-dimensional projected Riemann solvers on the edges of the cells of the computational domain. However, many researchers consider that within this approach, one-dimensional solvers lose much of their efficiency, mainly due to the fact that they do not take into account features of the solution propagating transversally to the cells boundaries. Moreover, the biasing in-troduced by one-dimensional solvers produces a reduction of the permissible Courant number when simulating multidimensional flows. For these reasons, there has been many attempts in the literature to build genuinely multidimensional Riemann solvers (see [2,3,34] and the references therein; see also [29] for a more detailed account on multidimensional upwinding), although most of them are specifically designed for the Euler or MHD equations.
In [35] Wendroff proposed a two-dimensional Riemann solver for the Euler equations which extended the one-dimensional HLL method [17], in which the two-dimensional interactions at cell corners were taken into account through the approximate solutions of two-dimensional Riemann problems. However, a drawback of Wendroff's solver is that it does not possess explicit expressions allowing its direct implementation; also, high-order extensions of the solver do not seem straighforward to build. Some years later, Balsara [2] modified Wendroff's formulation and proposed a two-dimensional HLL solver for the Euler and MHD equations on structured meshes, which included closed forms of the fluxes and allowed for an easier high-order extension. A more robust version of Balsara's solver was later proposed in [3], and an extension to unstructured meshes was presented in [4]. Another interesting reformulation of Wendroff's approach has been recently proposed in [34].
In this paper we introduce new types of genuinely two-dimensional incomplete Riemann solvers, which can be viewed as general AVM extensions of Balsara's HLL solver [3]. To achieve this, first we propose a reinterpretation of Balsara's solver in the particular case in which a simple four-wave model is considered for the two-dimensional corrections at cell vertices. These corrections can be viewed as suitable combinations of onedimensional numerical fluxes of HLL type. Once the role of the underlying one-dimensional HLL solver is clearly identified, it can be changed by an appropriate one-dimensional AVM numerical flux. This leads to twodimensional contributions of AVM type at the cell corners, which take into account the transversal features of the flow. It is noteworthy that the proposed solvers are applicable to general hyperbolic systems of conservation laws; in particular, they are not restricted to the Euler or MHD equations. An additional feature of these multidimensional AVM solvers is that they are theoretically stable up to a CFL number of unity.
Several families of basis functions for AVM methods were considered in [7,8]; in particular, good results were obtained with Chebyshev polynomials and Newman rational functions. In this paper we will consider PVM solvers based on the internal polynomials previously studied in [10], together with new families of RVM solvers based on suitable Padé approximations, which enjoy good stability properties. An additional advantage of these solvers is that no entropy fix is needed in the presence of sonic points.
We will focus in this paper on the comparison of the proposed methods in the first-order case, as done previously for one-dimensional RVM methods in [8]. As commented before, the choice of accurate first-order solvers is important when analizing complex scenarios. On the other hand, the use of our solvers as building blocks for higher order methods will be considered in a forthcoming paper.
In order to compare with other methods in the literature, applications to MHD equations will be considered. An additional difficulty that appears in the numerical resolution of MHD is the divergence-free constraint on the magnetic field, that has to be imposed in order to ensure the accuracy and stability of the numerical schemes. Several divergence cleaning techniques has been proposed in the literature: see [33] for a complete account on this topic. Regarding this, we present a novel technique for divergence cleaning based on Powell's [26] formulation of MHD equations. This technique has been included in our two-dimensional AVM solvers, modifying the numerical fluxes accordingly to [22]. Numerical experiments show that the results of our methods are comparable to those obtained with the well-known projection method [5], both in precision and computational time.
As remarked in [15], the inclusion of source terms in Balsara's solver seems to be a nontrivial matter. A way to handle source or coupling terms is to reformulate the system in nonconservative form, and then to apply the theory of path-conservative numerical schemes [25]. Taking into account the form of our multidimensional AVM methods, it is possible to extend them to the nonconservative case following the guidelines in [7,8]. This extension will be analized in a future paper, where we will construct highorder multidimensional AVM schemes that will be applied, in particular, for solving multilayer shallow water systems including depth variations.
The paper is organized as follows. After some preliminaries in Section 2, Balsara's multidimensional HLL solver is recalled in Section 3. Then, a brief review of AVM solvers is given in Section 4, together with the definitions of internal polynomials and Padé approximants. The core of the paper is Section 5, where multidimensional AVM solvers are introduced. Applications to MHD, including the proposed divergence cleaning technique, is considered in Section 6. Numerical experiments are then presented in Section 7. Finally, some conclusions are drawn in Section 8.

Preliminaries
Consider a two-dimensional hyperbolic system of conservation laws of the form We are interested in the numerical solution of (1) by means of finite volume methods on structured meshes. The spatial domain is divided into rectangular cells C ij = [x i−1/2 , x i+1/2 ] × [y j−1/2 , y j+1/2 ] of size ∆x × ∆y and center (x i , y j ), and the time step is denoted as t n = n∆t. where and the numerical fluxes verify In what follows, time dependence will be dropped unless necessary. Usual definitions of the numerical fluxes are based on one-dimensional Riemann solvers in the coordinate directions, which lead to expressions of the form F i+1/2,j = F(U ij , U i+1,j ) and G i,j+1/2 = G(U ij , U i,j+1 ) for certain flux functions F(U, V ) and G(U, V ). However, this approach does not take into account the genuinely two-dimensional effects that can be produced at each of the four vertices (x i±1/2 , y j±1/2 ) of a cell. A way to introduce these contributions in the definition of the numerical fluxes was proposed by Balsara in [2] and [3] for a HLL-type scheme. Specifically, the two-dimensional interactions at a given vertex are taken into account by considering the approximate solution of the two-dimensional Riemann problem that has as initial condition the four constant neighboring states at the vertex. Assuming that the strong interaction area in the solution of the two-dimensional Riemann problem is rectangular, Balsara was able to compute closed-form expressions for the numerical fluxes, thus improving previous work done in this direction by Wendroff [35]. For the sake of completeness, and to introduce some notation, Balsara's solver will be reviewed in the next section.

Balsara's multidimensional HLL solver
The purpose of this section is to give a brief review of Balsara's HLL multidimensional solver. Although all the details can be found in [2,3] (see also [34]), we give here the information that will be needed later to construct our schemes in Section 5.
To define the numerical flux at a given edge, Balsara's solver considers a linear convex combination of a one-dimensional HLL flux and two multidimensional contributions at the vertices. To be more precise, let us consider for a given cell C ij the six-cell stencil depicted in Figure 1. The flux in the x-direction through the right edge can be expressed as where F * i+1/2,j is the usual one-dimensional HLL flux between the cells C ij and C i+1,j , and F * i+1/2,j±1/2 represent the two-dimensional contributions at the vertices (x i+1/2 , y j±1/2 ); proper choices of the coefficients α, β and γ will be given later. Similar expressions can be stated for the y-flux G i,j+1/2 .
To define the two-dimensional flux F * i+1/2,j+1/2 , Balsara considered a two-dimensional Riemann problem with initial condition given by the constant states at each of the four cells which have (x i+1/2 , y j+1/2 ) as a common vertex. To simplify the notation, let us consider the local stencil depicted in Figure 2 (left); for example, for O ≡ (x i+1/2 , y j+1/2 ) the cell values would be given by U LD = U ij , U RD = U i+1,j , U LU = U i,j+1 and U RU = U i+1,j+1 . The initial condition for the local Riemann problem is then It is known that problem (1)-(4) has a self-similar solution roughly consisting of four directional one-dimensional Riemann problems arising at common edges, along with a region of strong interaction where complex structures appear; see, for example, [19,31,37] for more detailed descriptions.  In the design of numerical methods, however, it is impractical to resolve all the fine structures arising in the strong interaction region. Thus, this region is averaged into a single constant state, in a similar way to the intermediate state arising in the one-dimensional HLL solver. In Balsara's formulation, the region of strong interaction is assumed to be rectangular, although other options are possible (see, e. g., [34]).
To build a two-dimensional Riemann solver of HLL type, first we need to identify the maximal speeds in the coordinate directions. Following the notation in [3], let S U L and S U R denote the fastest left and right speeds at the edge between U LU and U RU (see Figure 2, right). An usual choice for S U R and S U L are given by where λ 1 x (U LU ) and λ N x (U RU ) denote the maximal left-and right-going wave speeds associated to the states U LU and U RU , respectively, whilē λ 1 x (U LU , U RU ) andλ N x (U LU , U RU ) are the maximal left-and right-going speeds arising from a linearized Riemann solver between the states U LU and U RU . Similar definitions apply for S D L and S D R . Speeds in the y-direction, S L D , S R D , S L U and S R U , are defined in an analogous way. Finally, the strong interaction region is assumed to be bounded by the maximal left and right wave speeds given by and the maximal downward and upward speeds A simplifying assumption can be made at this point, which consists in considering only the four speeds S L , S R , S D and S U in the wave structure of the solution. That is, we assume that S α = S D α = S U α for α = L, R, and S β = S L β = S R β for β = D, U . This leads to compact expressions for the strongly interacting state and its associated numerical fluxes. As later in Section 5 we will be interested on a four-wave model, this assumption will be made in what follows.
Assume now that we are in the subsonic case in both directions, that is, S L < 0 < S R and S D < 0 < S U . Figure 2 (right) shows the proposed wave structure at a time t = ∆t. The one-dimensional Riemann problems arising at the edges give rise to four HLL states; thus, for α = D, U , the states U Lα and U Rα give while for β = L, R, the states U βD and U βU are where F βα and G βα denote, respectively, the fluxes F (U βα ) and G(U βα ). The corresponding HLL fluxes are then given by Now, integration of (1) on the control volume (S L ∆t, S R ∆t)×(S D ∆t, S U ∆t)× (0, ∆t) leads to an explicit expression for the constant resolved state U * , namely Figure 3 shows the considered structure of the solution of the two-dimensional Riemann problem at a vertex. Notice that all the terms appearing in the definition of U * are known, except the transverse fluxes appearing in the last two terms, which have to be defined taking into account the corresponding intermediate state and the normal flux. For example, F * R is constructed using the state U * R and the flux G * R previously defined; in particular, notice that F * R = F (U * R ) in general. Balsara offers an explicit solution in [3] for the Euler equations, although similar constructions may be done for other systems: see Section 6 for the case of MHD equations.
To compute the resolved flux F * , equation (1) is integrated on the sub- Notice again that, in general, F * = F (U * ). Similarly, integration of (1) on (S L ∆t, S R ∆t) × (0, S U ∆t) × (0, ∆t) leads to As stated before, the fluxes (8) and (9) are used when the flow is subsonic in both spatial directions. When the flow is supersonic in some of the spatial directions, the corresponding expressions for F * and G * may be defined as the corresponding upwinded fluxes; the details can be found in [3, Sec. 2.1]. Finally, the assembling (3) is obtained from a time-space average of the normal flux through the edge (see Figure 4). This leads to the coefficients α = −S D,i+1/2,j+1/2 ∆t 2∆y , γ = S U,i+1/2,j−1/2 ∆t 2∆y and with the obvious notations. Strictly speaking, these values are valid only for the subsonic case in both directions. However, supersonic cases can be accommodated by simply substituting S U,R by S + U,R = max(S U,R , 0) and S D,L by S − D,L = min(S D,L , 0). Another option is to use Simpson's rule to approximate the integral, which leads to the simpler choice We have chosen this approach in the numerical experiments, although no significative differences have been found when the first option is used. More accurate approaches are possible: see, e. g., [34,Sect. 3.3.2]. With respect to the CFL condition needed to ensure stability, the multidimensional corrections allow for a maximum CFL number of unity: see [2] for a complete account on this point.

Some notes on AVM solvers
In Section 5 we will see how to construct more precise versions of Balsara's multidimensional scheme by changing the underlying HLL solver in a suitable way. To this end, and for the sake of completeness, we give in this section a brief survey on AVM solvers.
For simplicity, in this part will be enough to focus on one-dimensional systems. Therefore, we are interested in the numerical solution of the hyperbolic system by means of a finite volume method of the form with numerical flux given by Here Q i+1/2 denotes the numerical viscosity matrix, which determines the numerical diffusion of the scheme. It is worth noticing that Roe's method [28] can be written in the form (10) by considering the viscosity matrix Q i+1/2 = |A i+1/2 |, being A i+1/2 a Roe matrix for the system. The idea of PVM (Polynomial Viscosity Matrix) solvers introduced in [7] consists in approximating |A i+1/2 | by means of a suitable polynomial evaluation of A i+1/2 . More precisely, let p(x) be a polynomial approximation of |x| in the interval [−1, 1], and let λ i+1/2,max be the eigenvalue of A i+1/2 with maximum modulus (or an upper bound of it). Then the viscosity matrix of the PVM method associated to p(x) is given by Notice that the best p(x) approaches |x|, the closer the behaviour of the associated PVM scheme will be to that of Roe's method. A fundamental issue is that the spectral decomposition of A i+1/2 is not needed in the construction of a PVM method, but only a bound on its spectral radius. This feature makes PVM methods greatly efficient and applicable to systems in which the eigenstructure is not known or is difficult to obtain. In the cases in which a Roe matrix is not available or is costly to compute, A i+1/2 can be taken as the Jacobian matrix of the system evaluated at some average state. Several well-known schemes in the literature belong to the family of PVM methods: Lax-Friedrichs, Rusanov, HLL, FORCE, Roe, etc. (see [7]).
It is known that rational functions may provide much more precise approximations to |x| than polynomials. The RVM methods introduced in [8] take advantage of this fact, following the PVM idea but using rational functions as basis instead of polynomials. In what follows we will encompass both kind of solvers, PVM and RVM, under the common term AVM (Approximate Viscosity Matrix) solvers. AVM solvers have also been used in [9] to construct efficient versions of the Osher-Solomon method, following the ideas in [12]. Recently, in [10] some types of AVM solvers have been implemented in Jacobian-free form; this is particularly interesting when solving complex systems as, for example, relativistic MHD.
Finally, we must remark that the stability of an AVM scheme strongly depends on the properties of the underlying function f (x). In particular, it must verify the stability condition With respect to the choice of solvers, in [8] the authors proposed a family of PVM methods based on Chebyshev polynomials, which provide optimal approximations to the absolute value function, and could also be implemented in Jacobian-free form. In the same paper, RVM solvers based on Newman rational functions [23] where found to be the more efficient choice: they provided similar results as Roe's method for complex problems in MHD and multilayer shallow water systems, but with a much smaller computational cost. Despite their efficiency, a drawback of Chebyshev and Newman AVM solvers is that they do not satisfy the stability condition (12) strictly, so a slight modification has to be performed which may perturb the computation of the maximal speeds of propagation. On the other hand, in [10] a new family of internal polynomial approximations to |x| verifying (12) was proposed. Here we will consider this family of internal polynomials together with a class of rational Padé approximations based on [18], which also satisfy (12). For the sake of completeness, we briefly describe them in what follows.
Regarding efficiency issues, we have found that a direct implementation of p n (x) is a better option than using its recursive form. In Table 1 are given the coefficients of the polynomials Figure 5 (top left) shows the polynomials p n (x), for n = 1, 2, 3, 4. On the other hand, for given integers k, m ≥ 0, we consider the polynomials  and where we have used the notation (a) n = a(a + 1) . . . (a + n − 1) and (a) 0 = 1. Then, we follow [18] to build a sequence of Padé approximations of the form (the superscript [m/k] will be dropped unless necessary). It is not difficult to prove that each r n (x) verifies the properties 1-5 listed before for p n (x).  Table 2 the coefficients of the rational function r n+1 (x) written in the form for 0 ≤ k, m ≤ 3. For example, the family of Padé approximants of order [1/1] would be given by while those of order [2/2] would be r n+1 (x) = r n (x) r n (x) 4 + 10r n (x) 2 x 2 + 5x 4 5r n (x) 4 + 10r n (x) 2 x 2 + x 4 , n = 0, 1, 2, . . .   Finally, we would like to remark that an additional advantage of the considered approximations p n (x) and r n (x) is that they provide an automatic entropy fix to handle sonic flow. The reason for this is that neither p n (x) nor r n (x) crosses the origin: see [10,Sect. 4].
Remark 1. After extensive numerical investigation, our conclusion is that for relatively simple problems it is sufficient to consider a simple AVM solver, as for example HLL or the one based in the internal polynomial p 1 (x). However, when complex structures appear, the choice of more precise AVM solvers is a determinant factor in terms of efficiency (see [8]). This is even true when AVM solvers are used as building blocks in the design of highorder schemes, as it was demonstrated in [9]. 2

Multidimensional AVM solvers
The purpose of this section is to extend Balsara's multidimensional solver by changing, in an appropriate way, the underlying HLL flux by a more precise AVM solver. The key point to achieve this is the reinterpretation of Balsara's solver in such a way that the two-dimensional corrections are expressed as suitable combinations of one-dimensional numerical fluxes of HLL type, which in turn can be changed by another AVM flux. We will use the same notations as in Section 3 for the speeds, states and fluxes. In particular, we are interested in a simple four-wave model, so only the speeds S L , S R , S D and S U will be taken into account.
To begin with, consider the local stencil shown in Figure 2, and let us see how the two-dimensional resolved flux F * is defined. Let U * L and U * R be the HLL states given by (7), and let F * L and F * R be the transverse fluxes considered in Section 3. Notice that the explicit form of these transverse fluxes is problem-dependent: see Section 6 for the case of the MHD equations. Next, consider the AVM-type flux (10) in the x-direction given by is a function verifying the stability condition (12), A * is a Roe-like matrix such that and λ max is a bound on its spectral radius.
Remark 2. As stated in Section 3, in general F * α = F (U * α ), α = L, R, so A * is not exactly a Roe matrix. Its form will depend on how the transverse fluxes are built. In the cases in which such a matrix cannot be constructed, A * can be taken as the Jacobian of the flux F evaluated at some average state, so (14) holds only approximately. It was shown in [7] that the HLL flux can be viewed as an AVM solver based on the first-order polynomial with Thus, substituting in (13) and using (14), we get the following expression: where the explicit form (7) for the HLL states has been applied. Notice that the first two lines of (13) can be interpreted as two one-dimensional HLL fluxes associated, respectively, to the upper and lower parts of the stencil, while the last line accounts for the cross derivative of the flux G.
Assuming now that we are in the subsonic case in both directions, so substitution in (16) leads to exactly the same formula (8) for Balsara's resolved flux. This means that the two-dimensional resolved flux (8) can be reinterpreted as an AVM-type flux of the form (13). Analogous conclusions are obtained when the flux is supersonic in one or both directions. Similar computations can be performed to express the resolved flux (9) in the y-direction as where Q * y = g(B * ); here, g(x) = β 0 + β 1 x with Once the two-dimensional fluxes (8) and (9) have been expressed as the AVM-type fluxes (13) and (17), it is possible to extend these to build more precise two-dimensional solvers. This is achieved by changing the underlying HLL polynomial (15) by another suitable function satisfying (12), thus obtaining an AVM-type two-dimensional flux. It is important to remark that to build an AVM solver, only the maximum wave speeds S α are needed; despite simplicity, this is another reason why we consider a four-wave model. In particular, we always assume that the wave structure of the solution is that shown in Figure 2 (right). On the other hand, notice that in (13)-(17) the resolved states U * α are precisely the HLL states (6)-(7). We will conserve these resolved states even though the basis function f (x) used to build the fluxes (13)-(17) is changed. A possible way to introduce more structure in the resolved states, in the case of a PVM solver, could be through the equivalence of PVM solvers and simple Riemann solvers stated in [22]; this idea is currently under investigation.
An important consequence of the above considerations is that our twodimensional AVM solvers could be adapted to define path-conservative schemes for nonconservative systems, in the spirit of [25]. In particular, this allows to extend the two-dimensional AVM solvers to the case of hyperbolic systems with coupling and source terms. This will be the subject of a forthcoming paper, so we will not go into details here.
For the sake of completeness, we end this section by summarizing the construction of the horizontal numerical flux F i+1/2,j for a 2D AVM method. The vertical flux G i,j+1/2 is defined in a totally analogous way.
First of all, for generic cell values U 0 and U 1 , and maximal speeds of propagation S 0 and S 1 , define the HLL states Given a function f (x) satisfying (12), define the one-dimensional AVM fluxes where the viscosity matrices Q x and Q y are defined as in (11), in terms of suitable Roe or Jacobian matrices A x and A y , respectively. Consider now the six-cell stencil around the vertical edge depicted in Figure 1. Then: • The 1d flux through the edge is defined by • To build the 2d flux F * i+1/2,j+1/2 at the vertex (x i+1/2 , y j+1/2 ), first compute the four local speeds of propagation S L , S R , S U and S D (we omit the indexes for clarity) as explained in Section 3.
Subsonic case in both directions: S L < 0 < S R and S D < 0 < S U . Define the states and the 1d AVM fluxes Then build the transverse flux F * L in terms of U * L and G * L , and the transverse flux F * R in terms of U * R and G * R . Finally, the 2d flux is given by where Q * x is defined as explained in Section 5 in terms of a Roe-like matrix A * x satisfying (14); if such a matrix is not available, simply take the Jacobian of the physical flux F evaluated at some average state. Supersonic case in some or both directions: The flux F * i+1/2,j+1/2 is defined as the corresponding upwinded flux, exactly as in [3, Sec. 2.1].
• The 2d flux F * i+1/2,j−1/2 at the lower vertex is defined similarly. • Finally, the 1d and 2d fluxes are assembled as where, for simplicity, Simpson's rule has been considered (see the end of Section 3).
In short, the two-dimensional flux F * can be interpreted as follows. First, from the right states U RD and U RU in Figure 2, we build the associated state U * R and the AVM vertical flux G * R ; from these we define the transverse flux F * R , which may be viewed as a flux in the x-direction involving the states U RD and U RU and the contributions G RD and G RU in the y-direction. Next, the flux F * L is similarly defined from the left states U LD and U LU . Once U * L , U * R , F * L and F * R have been constructed, we use them as building blocks to define F * as an AVM-type flux: F * ≡ F AVM (U * L , U * R ), where the fluxes F (U * L ) and F (U * R ) are substituted by F * L and F * R respectively. The AVM solver used is determined by the choice of the viscosity matrix Q * x . Notice that the two-dimensional fluxes at vertices have been constructed as combination of several one-dimensional fluxes. Thus, for practical implementation only the definition of a one-dimensional AVM flux is really needed, what makes the computer coding simple and clear.
Finally, the choice of the time step is given by the CFL condition where each local time step ∆t ij is defined as where λ α ij denotes the maximal speed of propagation in the cell C ij in the α-direction, taking into account both the speeds at the edges and those arising in the two-dimensional Riemann problems at vertices. As remarked in Section 3 (see also [2]), the two-dimensional contributions allow the use of a maximal CFL number δ of unity. This represents an additional advantage with respect to standard methods in which a one-dimensional solver is applied dimension by dimension, for which the maximum CFL number is 0.5.

Application to ideal magnetohydrodynamics
In this section we focus on applications to ideal Newtonian MHD, although we must remark that the two-dimensional AVM schemes introduced in the previous section are applicable to general hyperbolic systems. We will use the following notation: ρ is density, v = (v x , v y , v z ) t is the velocity field, B = (B x , B y , B z ) t is the magnetic field, E is the total energy and P is the hydrostatic pressure. Then, the equations of MHD read as together with the divergence-free condition To close the system, an ideal equation of state P = (γ − 1)ρε is considered, where γ is the adiabatic constant and ε is the specific internal energy, which is related to the total energy by E = 1 2 ρv 2 + 1 2 B 2 + ρε. The total pressure is defined as P * = P + 1 2 B 2 , where 1 2 B 2 represents the magnetic pressure. On the other hand, define b = B/ √ ρ and the acoustic sound speed as a = γP/ρ. Then the fastest speeds of propagation in the α-direction, for α = x, y, are given by where For a detailed description of the spectral structure of system (22) see [6,30]. We remark that in the derivation of (22) from physical principles, the condition (23) is thoroughly used. However, it is also possible to derive a MHD model without imposing the divergence constraint in the process, as proposed in [26] (see also [13,27]). This leads to the nonconservative form of the MHD equations: Notice that if condition (23) is imposed, then (25) reduces to (22). The spectral structure of problem (25) was fully analized in [26]; in particular, the maximal speeds of propagation are the same as in the conservative case.
In the design of numerical methods for solving (22) it is fundamental to handle the divergence constraint (23) in a proper way. In particular, in the presence of shocks, standard numerical methods can produce large divergence errors with may lead to negative densities or pressures. Several methods have been proposed in the literature to impose the divergencefree condition numerically: see [33] for an in-depth review on this topic. A technique of particular interest to handle the divergence constraint was proposed by Powell in [26], which is based in the eight-waves model (25). Notice that taking divergence of the magnetic field equation in (25) gives so any possible non-zero divergence produced by the numerical scheme is advected by the velocity field and eventually leaves the domain through the boundaries. Thus, the nonconservative form (25) of the MHD system include a mechanism to handle the divergence constraint; we will back to this point later.
From now on we will concentrate on the two-dimensional case. Considering the vector of conserved variables U = (ρ, ρv, B, E) t , system (22) can be written in conservative form as where the physical fluxes are given by The transverse fluxes (see Remark 3) can be defined using values from the associated intermediate state and the normal fluxes. Thus, denoting the components of the states (18) as U * α = (u * α,i ) 8 i=1 and the fluxes (19) as taking into account the form of the fluxes F and G. The transverse fluxes G * β , for β = D, U , can be constructed in an analogous way from U * β and F * β . As commented before, some kind of divergence cleaning mechanism has to be added to the numerical scheme. For doing this we consider the nonconservative system (25), that can be expressed as where the source term is given by Alternatively, (26) can be written in the form where B 1 (U ) and B 2 (U ) are the 8 × 8 matrices with 0 the zero column vector in R 8 . Now, following [22], both the HLL states and the numerical fluxes must be properly modified in order to take into account the extra terms B 1 (U ) and B 2 (U ) appearing in (27). In particular, the HLL states should be modified as being U some average state. Accordingly, the numerical fluxes now should read as

Remark 4. Notice that for MHD the vectors
Finally, the definition of a two-dimensional AVM flux for the nonconservative system (27) follows exactly the same guidelines explained at the end of Section 5, substituting U HLL Notice that the proposed schemes are designed for solving the eight-waves model (25), so they automatically include a divergence cleaning procedure. Thus, if both the initial and boundary conditions satisfy the divergence constraint, it is expected that the scheme will keep it up to the accuracy of the truncation error. Notice also that the proposed nonconservative divergence cleaning technique could be applied, in a similar way, to any arbitrary numerical method.

Numerical results
The performances of the proposed multidimensional AVM schemes are analyzed in this section. To this end, we have chosen a number of tests that constitute standard references in MHD, in order that our schemes could be confronted with other methods in the literature.
An important issue in the construction of the schemes is the proper choice of the speeds of propagation. For this we have followed the guidelines stated in Section 3, with the expressions (24) for the maximal wave speeds. In particular, for the intermediate speedsλ α (U 0 , U 1 ) appearing in (5)  Another possible but more expensive choice is to consider U as the Roe state associated to the states U 0 and U 1 . In our computational experiments, both choices have led to similar results. Once the maximal speeds of propagation have been computed, the time step is determined from the CFL condition (21). In particular, we remark that for the 2D methods the theoretical admissible CFL number is 1, while for the 1D methods is 0.5.
We will denote as Int-n and Padé-[m/k] the AVM methods based on the internal polynomial p n (x) and the Padé approximation r [m/k] 1 (x), respectively. For each AVM method, 2D denotes its multidimensional version, where the numerical flux at an edge is given by (20); on the other hand, 1D refers to the corresponding one-dimensional solver applied dimension by dimension, which is equivalent to taking F i+1/2,j+1/2 = F * i+1/2,j in (20). Finally, unless otherwise stated, the divergence cleaning technique introduced in Section 6 will be used in all the numerical experiments.

Accuracy test
To study the accuracy of the proposed schemes, we consider a test from [21] for the ideal Euler equations with exact smooth solution given by ρ(x, y, t) = 1 + 0.2 sin(π(x + y − t/2)), v x = 1, v y = −1/2, P = 1,   Figure 6 shows a cut along the main diagonal of the final densities obtained on a 100 × 100 mesh with the 1D and 2D versions of HLL, Int-3 and Padé- [4/4] schemes. Notice that the Int-3 method is based on an polynomial of eighth degree, while the Padé- [4/4] method is based on a rational function of degree [8/8]. All the computations have been done using a common CFL number of 0.5, although the 2D methods perform well with CFL numbers up to 0.8. As it can be observed, for each scheme the two-dimensional contributions provide a more precise computation of the solution than their one-dimensional counterparts.
To check the behavior of the solutions in terms of the basis function chosen for an AVM 2D method, we have represented in Figure 7 the results obtained with the Int-n 2D schemes, for n = 1, 2, 3, 4, 5. It is clear that the precision of the solutions increases as the degree of the basis polynomial p n (x) increases. A similar remark is also valid for the Padé-[m/k] 2D schemes.  On the other hand, to check the convergence to the exact solution as the mesh is refined, in Figure 8 are represented the solutions obtained with different meshes of size N × N , for N = 25, 50, 100, 200. Table 3 shows the L 1 errors in the density variable, while the corresponding error curves are depicted in Figure 9. Notice that the Int-3 2D scheme provides very similar result as the Padé- [4/4]

Orszag-Tang vortex
The Orszag-Tang vortex [24] is a well-known model of complex MHD flow containing many significant features. Departing from an initial smooth state, the system develops a series of complex interactions between the different shock waves generated as the system evolves in the transition to turbulence.
Specifically, we consider the initial conditions proposed in [13]: We use this problem to test the robustness of several AVM 2D schemes equipped with the divergence cleaning technique introduced in Section 6. As it is well-known (see, e.g., [13]), negative pressures may appear if the divergence is not properly controlled. In order to compare, we have considered the HLL 2D, Int-3 2D and Padé- [4/4] 2D methods, with a common CFL number of 0.8. Figure 10 shows the densities and pressures obtained at time t = π. Although there is no accepted reference solution for this test, our results are in good agreement with those found in the literature: see, e.g., [2,13,20]. As it can be observed, the Int-3 2D and Padé-[4/4] 2D schemes provide a sharper resolution of the solution than Balsara's HLL 2D method. This can also be seen in Figure 11, where cuts along the main diagonal at different times are shown. As time evolves and the solution becomes more complex, the Padé-[4/4] scheme provides significantly the best approximations, followed by Int-3 2D.
On the other hand, we have compared the technique for imposing the divergence constraint defined in Section 6 with the well-known projection method introduced in [5]. Figure 12 shows cuts along the main diagonal of the solutions obtained with the Int-3 2D scheme at time t = 0.5 with both divergence cleaning techniques. As it can be observed, the results are comparable. Also, Table 4 shows the relative CPU times obtained when using both divergence cleaning methods. The results are again comparable, although our method seems to be slighly faster.

The rotor problem
The rotor problem, initially proposed in [1], is considered in this section. The initial condition consists of a dense rotating disk at the center of the domain, which is connected by means of a taper function with the ambient fluid at rest. The rotor is not in equilibrium, as the centrifugal forces are not balanced. Eventually, the rotating dense fluid will be confined into an oblate shape due to the action of the magnetic field. Given r 0 = 0.1, r 1 = 0.115, f = (r 1 − r)/(r 1 − r 0 ) and r = [(x − 0.5) 2 + (y − 0.5) 2 ] 1/2 , the initial conditions are given by  with B x = 2.5/ √ 4π, B y = 0 and P = 0.5. We take γ = 5/3. The computational domain is [0, 1] × [0, 1] and transmissive boundary conditions are considered. Figure 13 shows the solution obtained with the Padé-[4/4] 2D scheme at time t = 0.295 on a 200 × 200 mesh with CFL= 0.8. The results are in good agreement with those presented in [1,9,20,33]. In particular, it is worth noticing that the results produced with the first-order Padé-[4/4] 2D scheme have a similar quality as those obtained in [9] with a one-dimensional third-order approximate Chebyshev-DOT scheme on the same mesh.

Two-dimensional Riemann problem
In this section we consider a two-dimensional Riemann problem originally proposed in [11]. The initial condition, given in   are simple waves. Denoting the quadrants by Roman numbers as in Table  5, there is a rarefaction wave for problem I↔II and shocks for II↔III and III↔IV. The problem has been solved in the computational domain [−1, 1] ×    Figure 14 shows the contours of B x and B y computed with the Padé-[4/4] 2D scheme using a CFL of 0.8; for this test, the Int-3 2D scheme gives similar results. It was reported in [11] that some schemes present problems to keep B y constant across the shock in II↔III, and they can also produce strong distortions in B x and B y behind the rarefaction wave in I↔II. As it can be observed in Figure 14, none of those pathologies appear in the computed solutions, that can be directly compared with those presented in [11,20].
Finally, Figure 15 shows a comparison between the solution of the twodimensional Riemann problem obtained with the Padé-[4/4] 2D scheme and a one-dimensional reference solution. This serves to compare the quality of the solution of the I↔IV Riemann problem at x = 0.93, in a similar way as it was done in [11,20].

Spherical explosion
This test concerns the evolution of an explosion produced by an overpressured spherical region located at the center of the domain. It was originally proposed in [36], although we have considered here the data proposed in [32]. Specifically, the initial pressure is P = 100 inside a circle of radius r = 10 and P = 1 outside. The density is ρ = 1 on the whole domain and the velocity field is set to zero. As in [32], three different values of the initial magnetic field are considered, with a common value of the adiabatic con- propagates symmetrically in the radial direction, as expected for the Euler equations. This behavior can be seen in Figure 16, which shows the results of the simulation at time t = 3. If the initial magnetic field is taken as (B x , B y , B z ) = (0, 5/ √ π, 0), the increasing of the strength of the magnetic field produces a slightly elongated shock in the y-direction: see Figure 17.
For an even stronger initial magnetic field (B x , B y , B z ) = (0, 50/ √ π, 0), the solution becomes highly anisotropic, with essentially no displacement of fluid in the orthogonal direction to the magnetic field. Figure 18 shows the results obtained at time t = 1.05, which is the time considered in [32] for which the perturbation has not yet reached the boundary.
In all the three cases so far considered, our results are in good agreement with those presented in [32,36].

Conclusions
We have introduced a general class of genuinely two-dimensional incomplete Riemann solvers for systems of conservation laws. Departing from a reinterpretation of Balsara's two-dimensional HLL scheme, we have extended it by substituting the underlying one-dimensional HLL solver with an arbitrary AVM-type numerical flux. This allows, in particular, to control the numerical diffusion of the scheme by using suitable aproximations to the absolute value function used as basis of the AVM flux. The numerical flux at an edge is built as a linear convex combination of a one-dimensional AVM flux and two multidimensional AVM corrections at the corners which take into account transversal features of the flow. In principle, the twodimensional AVM solvers could be extended to three-dimensional problems, although this point is not addressed in this paper. The proposed multidimensional AVM methods can be applied to arbitrary hyperbolic systems. In this paper we have focused on applications to the MHD equations, in order to make comparisons with other existing methods in the literature. Additionally, we have proposed a new technique for imposing the divergence-free condition on the magnetic field, which gives comparable results as the standard projection method. The performances of the proposed schemes have been assessed in a number of challenging problems in MHD.
Future work in this line of research include the extension of the multidimensional AVM solvers to the case of nonconservative hyperbolic systems, the construction of higher order versions of the schemes, as well as applica- tions to different hyperbolic systems.