Efficient conservative ADER schemes based on WENO reconstruction and space-time predictor in primitive variables

We present a new version of conservative ADER-WENO finite volume schemes, in which both the high order spatial reconstruction as well as the time evolution of the reconstruction polynomials in the local space-time predictor stage are performed in primitive variables, rather than in conserved ones. To obtain a conservative method, the underlying finite volume scheme is still written in terms of the cell averages of the conserved quantities. Therefore, our new approach performs the spatial WENO reconstruction twice: the first WENO reconstruction is carried out on the known cell averages of the conservative variables. The WENO polynomials are then used at the cell centers to compute point values of the conserved variables, which are subsequently converted into point values of the primitive variables. This is the only place where the conversion from conservative to primitive variables is needed in the new scheme. Then, a second WENO reconstruction is performed on the point values of the primitive variables to obtain piecewise high order reconstruction polynomials of the primitive variables. The reconstruction polynomials are subsequently evolved in time with a novel space-time finite element predictor that is directly applied to the governing PDE written in primitive form. The resulting space-time polynomials of the primitive variables can then be directly used as input for the numerical fluxes at the cell boundaries in the underlying conservative finite volume scheme. Hence, the number of necessary conversions from the conserved to the primitive variables is reduced to just one single conversion at each cell center. We have verified the validity of the new approach over a wide range of hyperbolic systems, including the classical Euler equations of gas dynamics, the special relativistic hydrodynamics (RHD) and ideal magnetohydrodynamics (RMHD) equations, as well as the Baer-Nunziato model for compressible two-phase flows. In all cases we have noticed that the new ADER schemes provide less oscillatory solutions when compared to ADER finite volume schemes based on the reconstruction in conserved variables, especially for the RMHD and the Baer-Nunziato equations. For the RHD and RMHD equations, the overall accuracy is improved and the CPU time is reduced by about 25 %. Because of its increased accuracy and due to the reduced computational cost, we recommend to use this version of ADER as the standard one in the relativistic framework. At the end of the paper, the new approach has also been extended to ADER-DG schemes on space-time adaptive grids (AMR).


Introduction
Since their introduction by Toro and Titarev (Toro et al. ; Titarev and Toro ; Toro and Titarev ; Titarev and Toro ; Toro and Titarev ), ADER (arbitrary high order derivatives) schemes for hyperbolic partial differential equations (PDE) have been improved and developed along different directions. A key feature of these methods is their ability to achieve uniformly high order of accuracy in space and time in a single step, without the need of intermediate Runge-Kutta stages (Pareschi et al. ; Pidatella et al. ), by exploiting the approximate solution of a Generalized Riemann Problem (GRP) at cell boundaries. ADER schemes have been first conceived within the finite volume (FV) framework, but they were soon extended also to the discontinuous Galerkin (DG) finite element framework (Dumbser and Munz ; Taube et al. ) and to a unified formulation of FV and DG schemes, namely the so-called P N P M approach (Dumbser et al. a). In the original ADER approach by Toro and Titarev, the approximate solution of the GRP is obtained through the solution of a conventional Riemann problem between the boundary-extrapolated values, and a sequence of linearized Riemann problems for the spatial derivatives. The required time derivatives in the GRP are obtained via the so-called Cauchy-Kowalevski procedure, which consists in replacing the time derivatives of the Taylor expansion at each interface with spatial derivatives of appropriate order, by resorting to the strong differential form of the PDE. Such an approach, though formally elegant, becomes prohibitive or even impossible as the complexity of the equations increases, especially for multidimensional problems and for relativistic hydrodynamics and magneto-hydrodynamics. On the contrary, in the modern reformulation of ADER (Dumbser et al. b; Dumbser et al. a; Balsara et al. ), the approximate solution of the GRP is achieved by first evolving the data locally inside each cell through a local space-time discontinuous Galerkin predictor (LSDG) step that is based on a weak form of the PDE, and, second, by solving a sequence of classical Riemann problems along the time axis at each element interface. This approach has the additional benefit that it can successfully cope with stiff source terms in the equations, a fact which is often encountered in physical applications. For these reasons, ADER schemes have been applied to real physical problems mostly in their modern version. Notable examples of applications include the study of Navier-Stokes equations, with or without chemical reactions (Hidalgo and Dumbser ; Dumbser ), geophysical flows (Dumbser et al. ), complex three-dimensional free surface flows (Dumbser ), relativistic magnetic reconnection (Dumbser and Zanotti ; Zanotti and Dumbser ), and the study of the Richtmyer-Meshkov instability in the relativistic regime (Zanotti and Dumbser ). In the last few years, ADER schemes have been enriched with several additional properties, reaching a high level of flexibility. First of all, ADER schemes have been soon extended to deal with non-conservative systems of hyperbolic PDE (Toro and Hidalgo ; Dumbser et al. ; Dumbser et al. ), by resorting to path-conservative methods (Parés and Castro ; Pares ). ADER schemes have also been extended to the Lagrangian framework, in which they are currently applied to the solution of multidimensional problems on unstructured meshes for various systems of equations, (Boscheri and Dumbser ; Dumbser and Boscheri ; Boscheri et al. a; Boscheri et al. b; Boscheri and Dumbser ). On another side, ADER schemes have been combined with Adaptive Mesh Refinement (AMR) techniques (Dumbser et al. ; Zanotti and Dumbser ), exploiting the local properties of the discontinuous Galerkin predictor step, which is applied cell-by-cell irrespective of the level of refinement of the neighbour cells. Moreover, ADER schemes have also been used in combination with discontinuous Galerkin methods, even in the presence of shock waves and other discontinuities within the flow, thanks to a novel a posteriori sub-cell finite volume limiter technique based on the MOOD approach (Clain et al. ; Diot et al. ), that is designed to stabilize the discrete solution wherever the DG approach fails and produces spurious oscillations or negative densities and pressures (Dumbser et al. ; Zanotti et al. a; Zanotti et al. b).
The various implementations of ADER schemes mentioned so far differ under several aspects, but they all share the following common features: they apply the local spacetime discontinuous Galerkin predictor to the conserved variables, which in turn implies that, if a WENO finite volume scheme is used, the spatial WENO reconstruction is also performed in terms of the conserved variables. Although this may be regarded as a reasonable choice, it has two fundamental drawbacks. The first one has to do with the fact that, as shown by Munz (), the reconstruction in conserved variables provides the worst shock capturing fidelity when compared to the reconstruction performed either in primitive or in characteristic variables. The second drawback is instead related to computational performance. Since the computation of the numerical fluxes requires the calculation of integrals via Gaussian quadrature, the physical fluxes must necessarily be computed at each space-time Gauss-Legendre quadrature point. However, there are systems of equations (e.g. the relativistic hydrodynamics or magnetohydrodynamics equations) for which the physical fluxes can only be written in terms of the primitive variables. As a result, a conversion from the conserved to the primitive variables is necessary for the calculations of the fluxes, and this operation, which is never analytic for such systems of equations, is rather expensive. For these reasons it would be very desirable to have an ADER scheme in which both the reconstruction and the subsequent local space-time discontinuous Galerkin predictor are performed in primitive variables. It is the aim of the present paper to explore this possibility. It is also worth stressing that in the context of high order finite difference Godunov methods, based on traditional Runge-Kutta discretization in time, the reconstruction in primitive variables has been proved to be very successful by Del Zanna et al. () in their ECHO general relativistic code (see also Bucciantini and Del Zanna ; Zanotti et al. ). In spite of the obvious differences among the numerical schemes adopted, the approach that we propose here and the ECHO-approach share the common feature of requiring a single (per cell) conversion from the conserved to the primitive variables.
The plan of the paper is the following: in Section  we describe the numerical method, with particular emphasis on Section . and on Section ., where the spatial reconstruction strategy and the local space-time discontinuous Galerkin predictor in primitive variable are described. The results of our new approach are presented in Section  for a set of four different systems of equations. In Section  we show that the new strategy can also be extended to pure discontinuous Galerkin schemes, even in the presence of space-time adaptive meshes (AMR). Finally, Section  is devoted to the conclusions of the work.

Numerical method
We present our new approach for purely regular Cartesian meshes, although there is no conceptual reason preventing the extension to general curvilinear or unstructured meshes, which may be considered in future studies.

Formulation of the equations
We consider hyperbolic systems of balance laws that contain both conservative and non-conservative terms, i.e.

∂Q ∂t
where Q ∈ Q ⊂ R ν is the state vector of the ν conserved variables, which, for the typical gas dynamics equations, are related to the conservation of mass, momentum and energy.
is the flux tensor a for the conservative part of the PDE system, while represents the nonconservative part of it. Finally, S(Q) is the vector of the source terms, which may or may not be present. In the follow up of our discussion it is convenient to recast the system () in quasilinear form as where A(Q) = [A x , A y , A z ] = ∂F(Q)/∂Q + B(Q) accounts for both the conservative and the non-conservative contributions. As we shall see below, a proper discretization of Eq. () can provide the time evolution of the conserved variables Q, but when the primitive variables V are adopted instead, Eq. () translates into In the following we suppose that the conserved variables Q can always be written analytically in terms of the primitive variables V, i.e. the functions are supposed to be analytic for all PDE systems under consideration. On the contrary, the conversion from the conserved to the primitive variables, henceforth the cons-toprim conversion, is not always available in closed form, i.e. the functions may not be analytic (e.g. for relativistic hydrodynamics and magnetohydrodynamics to be discussed in Section .), thus requiring an approximate numerical solution. As a result, the matrix ( ∂Q ∂V ) - , which in principle could be simply computed as in practice it cannot be obtained in this manner, but it must be computed as where we have introduced the notation which will be used repeatedly below. Since Q(V) is supposed to be analytic, the matrix M can be easily computed. Equation () will serve us as the master equation to evolve the cell averages of the conserved variables Q via a standard finite volume scheme. However, both the spatial WENO reconstruction and the subsequent LSDG predictor will act on the primitive variables V, hence relying on the alternative formulation given by Eq. (). The necessary steps to obtain such a scheme are described in the Sections .-. below.

The finite volume scheme
In Cartesian coordinates, we discretize the computational domain through space-time control volumes and t = t n+t n . Integration of Eq. () over I ijk yields the usual finite volume discretization where the cell averagē is the spatial average of the vector of conserved quantities at time t n . In Eq. () we recognize two different sets of terms, namely those due to the conservative part of the system (), and those coming from the non-conservative part of it. In the former set we include the three time-averaged fluxes and the space-time averaged source term We emphasize that the terms v h in Eqs. ()-(), as well as in the few equations below, are piecewise space-time polynomials of degree M in primitive variables, computed according to a suitable LSDG predictor based on the formulation (), as we will discuss in Section .. This marks a striking difference with respect to traditional ADER schemes, in which such polynomials are instead computed in conserved variables and are denoted as q h (see, e.g. Hidalgo and Dumbser ). The integrals over the smooth part of the non-conservative terms in Eq. () yield the following contribution, According to this approach, the following path integrals must be prescribed where (s) is a path joining the left and right boundary extrapolated states vh and v + h in state space of the primitive variables. The simplest option is to use a straight-line segment path Pragmatic as it is, b the choice of the path () allows to evaluate the terms D i in () as that we compute through a three-point Gauss-Legendre formula (Dumbser et al. ; Dumbser and Toro a; Dumbser and Toro b). The computation of the numerical fluxesf i in Eq. () requires the use of an approximate Riemann solver, see Toro (). In this work we have limited our attention to a local Lax-Friedrichs flux (Rusanov flux) and to the Osher-type flux proposed in Dumbser and Toro (b), Dumbser and Toro (a), Castro et al. (). Both of them can be written formally as where D i ≥  is a positive-definite dissipation matrix that depends on the chosen Riemann solver. For the Rusanov flux it simply reads where |s max | is the maximum absolute value of the eigenvalues admitted by the PDE and I is the identity matrix. The matrix M is a Roe matrix that allows to write the jumps in the conserved variables in terms of the jump in the primitive variables, i.e.
Since M = ∂Q/∂V, the Roe matrix M can be easily defined by a path integral as which in the case of the simple straight-line segment path () leads to the expression In the case of the Osher-type flux, on the other hand, the dissipation matrix reads with the usual definition of the matrix absolute value operator The path in Eqs. () and () is the same segment path adopted in () for the computation of the jumps D i .

A novel WENO reconstruction in primitive variables
Since we want to compute the time averaged fluxes [cf.
Eqs. ()-()] and the space-time averaged sources [cf. Eq. ()] directly from the primitive variables V, it is necessary to reconstruct a WENO polynomial in primitive variables. However, the underlying finite volume scheme () will still advance in time the cell averages of the conserved variablesQ n ijk , which are the only known input quantities at the reference time level t n . Hence, the whole procedure is performed through the following three simple steps: . We perform a first standard spatial WENO reconstruction of the conserved variables starting from the cell averagesQ n ijk . This allows to obtain a reconstructed polynomial w h (x, y, z, t n ) in conserved variables valid within each cell. . Since w h (x, y, z, t n ) is defined at any point inside the cell, we simply evaluate it at the cell center in order to obtain the point value Q n ijk = w h (x i , y j , z k , t n ). This conversion from cell averagesQ n ijk to point values Q n ijk is the main key idea of our new method, since the simple identity Q n ijk =Q n ijk is valid only up to second order of accuracy! After that, we perform a conversion from the point-values of the conserved variables to the point-values in primitive variables, i.e. we apply Eq. (), thus obtaining the corresponding primitive variables V n ijk = V(Q n ijk ) at each cell center. This is the only step in the entire algorithm that needs a conversion from the conservative to the primitive variables. . Finally, from the point-values of the primitive variables at the cell centers, we perform a second WENO reconstruction to obtain a reconstruction polynomial in primitive variables, denoted as p h (x, y, z, t n ). This polynomial is then used as the initial condition for the new local space-time DG predictor in primitive variables described in Section .. As for the choice of the spatial WENO reconstruction, we have adopted a dimension-by-dimension reconstruction strategy, discussed in full details in our previous works (see Dumbser et al. ; Dumbser et al. ; Zanotti and Dumbser ). Briefly, we first introduce space-time reference coordinates ξ , η, ζ , τ ∈ [, ], defined by and, along each spatial direction, we define a basis of polynomials {ψ l (λ)} M+ l= , each of degree M, formed by the M +  Lagrange interpolating polynomials, that pass through the M +  Gauss-Legendre quadrature nodes {μ k } M+ k= . According to the WENO philosophy, a number of stencils is introduced such that the final polynomial is a data-dependent nonlinear combination of the polynomials computed from each stencil. Here, we use a fixed number N s of onedimensional stencils, namely N s =  for odd order schemes (even polynomials of degree M), and N s =  for even order schemes (odd polynomials of degree M). For example, focusing on the x direction for convenience, every stencil along x is formed by the union of M +  adjacent cells, i.e.
where L = L(M, s) and R = R(M, s) are the spatial extension of the stencil to the left and to the right. c Now, an important difference emerges depending on whether we are reconstructing the conserved or the primitive variables. In the former case, corresponding to the computation of w h (x, y, z, t n ) at step  above, we require that the reconstructed polynomial must preserve the cellaverages of the conserved variables over each element I ijk .
Since the polynomials reconstructed along the x direction can be written as Equations () provide a system of M +  linear equations for the unknown coefficientsŵ n,s ijk,r , which is conveniently solved through linear algebra packages. Once this operation has been performed for each stencil, we construct a data-dependent nonlinear combination of the resulting polynomials, i.e. . The whole procedure must be repeated along the two directions y and z. Hence, although each direction is treated separately, the net effect provides a genuine multidimensional reconstruction. We now proceed with the key step of the new algorithm presented in this paper and compute the point values of the conserved quantities at the cell centers, simply by evaluating the reconstruction polynomials in the barycenter of each control volume: These point values of the conserved quantities Q n ijk are now converted into point values of the primitive variables V n ijk , which requires only a single cons-to-prim conversion per cell. In RHD and RMHD, this is one of the most expensive and most delicate parts of the entire algorithm: The reconstruction polynomials in primitive variables are spanned by the same basis functions ψ r (ξ ) used for w h , hence According to step  listed above, we now require that the reconstructed polynomial must interpolate the pointvalues of the primitive variables at the centers of the cells forming each stencil, i.e.
The reconstruction equations () will also generate a system of M +  linear equations for the unknown coefficientsp n,s ijk,r . The rest of the WENO logic applies in the same way, leading to We emphasize that thanks to our polynomial WENO reconstruction (instead of the original point-wise WENO reconstruction of Jiang and Shu ), the point-value of w h (x, y, z, t n ) at each cell center, which is required at step  above, is promptly available after evaluating the basis functions at the cell center. In other words, there is no need to perform any special transformation from cell averages to point-values via Taylor series expansions, like in Buchmüller and Helzel (), Buchmüller et al. (). On the other hand, since the WENO reconstruction is performed twice, once for the conserved variables and once for the primitive variables, we expect that our new approach will become convenient in terms of computational efficiency only for those systems of equations characterized by relations V(Q) that cannot be written in closed form. In such circumstances, in fact, reducing the number of cons-toprim conversions from M(M + ) d+ + d(M + ) d in d space dimensions (due to the space-time predictor and the numerical flux computation in the finite volume scheme) to just one single conversion per cell will compensate for the double WENO reconstruction in space that we must perform. On the contrary, for systems of equations, such as the compressible Euler, for which the cons-to-prim conversion is analytic, no benefit will be reported in terms of computational efficiency, but still a significant benefit will be reported in terms of numerical accuracy. All these comments will be made quantitative in Section .

.. Description of the predictor
As already remarked, the computation of the fluxes through the integrals ()-() is more conveniently performed if the primitive variables are available at each space-time quadrature point. In such a case, in fact, no conversion from the conserved to the primitive variables is required. According to the discussion of the previous Section, it is possible to obtain a polynomial p h (x, y, z, t n ) in primitive variables at the reference time t n . This is however not enough for a high accurate computation of the numerical fluxes, and p h (x, y, z, t n ) must be evolved in time, locally for each cell, in order to obtain a polynomial v h (x, y, z, t) approximating the solution at any time in the range [t n ; t n+ ].
To this extent, we need an operation, to be performed locally for each cell, which uses as input the high order polynomial v h obtained from the WENO reconstruction, and gives as output its evolution in time, namely This can be obtained through an element-local space-time discontinuous Galerkin predictor that is based on the weak integral form of Eq. (). From a mathematical point of view, Eq. () is a hyperbolic system in non-conservative form. Therefore, the implementation of the space-time discontinuous Galerkin predictor follows strictly the strategy already outlined in Dumbser et al. () for nonconservative systems. Here we recall briefly the main ideas, focusing on the novel aspects implied by the formulation of Eq. (). The sought polynomial v h (x, y, z, t) is supposed to be expanded in space and time as where the degrees of freedomv n l are the unknowns. The space-time basis functions θ l are given by a dyadic product of the Lagrange interpolation polynomials that pass through the Gauss-Legendre quadrature points, i.e. the tensor-product quadrature points on the hypercube [, ] d+ , see Stroud (). The system () is first rephrased in terms of the reference coordinates τ and ξ = (ξ , η, ζ ), yielding Expression () is then multiplied by the piecewise spacetime polynomials θ k (ξ , η, ζ , τ ) and integrated over the space-time reference control volume, thus providing where we have replaced V with its discrete representation v h . Integrating the first term by parts in time yields Equation () is an element-local nonlinear algebraic equation that must be solved locally for each grid-cell in the unknownsv n l . In practice, we solve the system of equations () through a discrete Picard iteration, see Dumbser and Zanotti (), Hidalgo and Dumbser (), where additional comments about its solution can be found.

.. An efficient initial guess for the predictor
A proper choice of the initial guess for each of the spacetime degrees of freedomv l can improve the convergence of the Picard process. The easiest strategy is to set v h (x, t) = p h (x, t n ) i.e. the reconstruction polynomial is simply extended as a constant in time. This is, however, not the best approach. A better strategy for obtaining a good initial guess for the LSDG predictor was presented in Hidalgo and Dumbser (), and it is based on the implementation of a MUSCL scheme for the explicit terms, plus a second-order Crank-Nicholson scheme in case stiff source terms are present. In the following, we refer to this version of the initial guess for the LSDG predictor as the MUSCL-CN initial guess. If the source terms are not stiff, however, an even more efficient approach is possible which is based on a space-time extension of multi-level Adams-Bashforth-type ODE integrators. For that purpose, the space-time polynomial denoted by v n- h (x, t) obtained during the previous time step [t n- , t n ] is simply extrapolated in time to the new time step [t n , t n+ ] by simple L projection: The comparison has been performed for the isentropic vortex solution and the numbers have been normalized to the value obtained with the traditional MUSCL-CN initial guess (see Section 2.4.2 for more details).
In terms of the degrees of freedomv n l andv n- l this relation becomes with τ =  + τ t n t n- and t n- = t nt n- . In the following, we refer to this second version of the initial guess for the LSDG predictor as the Adams-Bashforth (AB) initial guess. In Table  we show a comparison among the performances of the LSDG predictor with these two different implementations of the initial guess.

Numerical tests with the new ADER-WENO finite volume scheme in primitive variables
In the following we explore the properties of the new ADER-WENO finite volume scheme by solving a wide set of test problems belonging to four different systems of equations: the classical Euler equations, the relativistic hydrodynamics (RHD) and magnetohydrodynamics (RMHD) equations and the Baer-Nunziato equations for compressible two-phase flows. For the sake of clarity, we introduce the notation ' ADER-Prim' to refer to the novel approach of this work for which both the spatial WENO reconstruction and the subsequent LSDG predictor are performed on the primitive variables. On the contrary, we denote the traditional ADER implementation, for which both the spatial WENO reconstruction and the LSDG predictor are performed on the conserved variables, as ' ADER-Cons' . In a few circumstances, we have also compared with the ' ADER-Char' scheme, namely a traditional ADER scheme in which, however, the spatial reconstruction is performed on the characteristic variables. In this Section we focus our attention on finite volume schemes, which, according to the notation introduced in Dumbser et al. (a), are denoted as P  P M methods, where M is the degree of the approximating polynomial. In Section  a brief account is given to discontinuous Galerkin methods, referred to as P N P N methods, for which an ADER-Prim version is also possible.

Euler equations
First of all we consider the solution of the classical Euler equations of compressible gas dynamics, for which the vectors of the conserved variables Q and of the fluxes f x , f y and f z are given respectively by Here v x , v y and v z are the velocity components, p is the pressure, ρ is the mass density, / is the total energy density, while γ is the adiabatic index of the supposed ideal gas equation of state, which is of the kind p = ρ (γ -), being the specific internal energy.

.. D isentropic vortex
It is important to assess the convergence properties of the new scheme, in particular comparing with the traditional ADER scheme in conserved and in characteristic variables.
To this extent, we have studied the two-dimensional isentropic vortex, see e.g. Hu and Shu (). The initial conditions are given by a uniform mean flow, to which a per-turbation is added, such that Whatever the perturbation δT in the temperature is, it is easy to verify that there is not any variation in the specific entropy s = p/ρ γ , and the flow is advected smoothly and isentropically with velocity v = (, ,  with r  = (x -)  + (y -)  , vortex strength =  and adiabatic index γ = .. Table  contains the results of our calculation, in which we have compared the convergence properties of three different finite volume ADER schemes: ADER-Prim, ADER-Cons and ADER-Char, obtained with the Osher-type Riemann solver, see Dumbser and Toro (b). While all the schemes converge to the nominal order, it is interesting to note that the smallest L  error is obtained for the new ADER finite volume scheme in primitive variables, and that the difference with respect to the other two reconstructions increases with the order of the method.
In addition to the convergence properties, we have compared the performances of the Adams-Bashforth version of the initial guess for the LSDG predictor with the traditional version based on the MUSCL-CN algorithm. The comparison has been performed over a  ×  uniform grid. The results are shown in Table , from which we conclude that the Adams-Bashforth initial guess is indeed computationally more efficient in terms of CPU time. However, we have also experienced that it is typically less robust, and in some of the most challenging numerical tests discussed in the rest of the paper we had to use the more traditional MUSCL-CN initial guess.

.. Sod's Riemann problem
We have then solved the classical Riemann problem named after Sod (Sod ), assuming an adiabatic index γ = ., and evolved until t final = .. In spite of the fact that this is a one-dimensional test, we have evolved this problem in two spatial dimensions over the domain [, ] × [-., .], using periodic boundary conditions along the passive y direction. In Figure  we show the comparison among the solutions obtained with ADER-Prim, ADER-Cons and ADER-Char, together with the exact solution provided in Toro (). We have adopted the finite volume scheme at the fourth order of accuracy, namely the P  P  scheme, in combination with the Rusanov numerical flux and using  cells along the x-direction. Although  all of the ADER implementations show a very good agreement with the exact solution, a closer look at the tail of the rarefaction, highlighted in the bottom right panel, reveals that the ADER-Cons scheme is actually the worst one, while the solution obtained with ADER-Prim is more similar to the reconstruction in characteristic variables.
On the contrary, in terms of CPU-time, ADER-Prim is not convenient for this system of equations because the price paid for performing the double WENO reconstruction in space is not significantly compensated by the reduced number of conversions that are needed from the conserved to the primitive variables. Table  reports the CPU times, normalized with respect to the ADER-Prim implementation, for different orders of accuracy, showing that the ADER-Prim scheme is ∼ % slower than the traditional ADER-Cons scheme. As we will see in Table  of Section ., the comparison will change in favor of ADER-Prim schemes, when the relativistic equations are solved instead.

.. Interacting blast waves
The interaction between two blast waves was first proposed by Woodward and Colella () and it is now a standard test for computational fluid dynamics. The initial conditions are given by where the adiabatic index is γ = .. We have evolved this problem in two spatial dimensions over the domain [-., .] × [-., .], using reflecting boundary conditions in x direction and periodic boundary conditions along the y direction. The results of our calculations, obtained with the P  P  scheme, are reported in Figure , where only the one-dimensional cuts are shown. The number of cells chosen along the x-direction, namely N x = , is not particularly large, at least for this kind of challenging problem. This has been intentionally done to better highlight potential differences among the two alternative ADER-Prim and ADER-Cons schemes. As it turns out from the figure, the two methods are very similar in terms of accuracy: the sharp peak in the density at time t = . (left panel) is somewhat better resolved through the ADER-Prim, while the opposite is true for the highest peak at time t = . (right panel). On the overall, however, the two schemes perform equally well for this test.

.. Double Mach reflection problem
As a representative test for the Euler equations in two space dimensions, we have considered the double Mach reflection problem, which implies the interaction of several waves. The dynamics of this problem is triggered by a shock wave propagating towards the right with a Mach number M = , and intersecting the x-axis at x = / with  As a tentative conclusion about the performances of ADER-Prim for the Euler equations, we may say that, although it is the most accurate on smooth solutions (see Table ), and comparable to a traditional ADER with reconstruction in characteristic variables, it is computationally more expensive than ADER-Cons and ADER-Char. Hence, ADER-Prim will rarely become the preferred choice in standard applications for the Euler equations.

Relativistic hydrodynamics and magnetohydrodynamics
From a formal point of view, the equations of special relativistic hydrodynamics and magnetohydrodynamics can be written in conservative form like the classical Euler equations (see, however, the comments below), namely as in Eq. (), with the vectors of the conserved variables and of the corresponding fluxes given by where the conserved variables (D, S j , U, B j ) can be expressed as d while the spatial projection of the energy-momentum tensor of the fluid is (Del Zanna et al. ) Here ijk is the Levi-Civita tensor and δ ij is the Kronecker symbol. We have used the symbol h =  + + p/ρ to denote the specific enthalpy of the plasma and in all our calculations the usual ideal gas equation of state has been assumed.
The components of the electric and of the magnetic field in the laboratory frame are denoted by E i and B i , while the Lorentz factor of the fluid with respect to this reference frame is W = (v  ) -/ . We emphasize that the electric field does not need to be evolved in time under the assumption of infinite electrical conductivity, since it can always be computed in terms of the velocity and of the magnetic field as E =v × B.
Although formally very similar to the classical gas dynamics equations, their relativistic counterpart present two fundamental differences. The first one is that, while the physical fluxes f i of the classical gas dynamics equations can be written analytically in terms of the conserved variables, i.e. f i = f i (Q), those of the relativistic hydrodynamics (or magnetohydrodynamics) equations need the knowledge of the primitive variables, i.e. f i = f i (V) for RMHD. The second difference is that, in the relativistic case, the conversion from the conserved to the primitive variables, i.e. the operation (D, S j , U, B j ) → (ρ, v i , p, B i ), is not analytic, and it must be performed numerically through some appropriate iterative procedure. Since in an ADER scheme such a conversion must be performed in each space-time degree of freedom of the space-time DG predictor and at each Gaussian quadrature point for the computation of the fluxes in the finite volume scheme, we may expect a significant computational advantage by performing the WENO reconstruction and the LSDG predictor directly on the primitive variables. In this way, in fact, the conversion (D, S j , U, B j ) → (ρ, v i , p, B i ) is required only once at the cell center (see Section .), and not in each space-time degree of freedom of the predictor and at each Gaussian point for the quadrature of the numerical fluxes. We emphasize that the choice of the variables to reconstruct for the relativistic velocity is still a matter of debate. The velocity v i may seem the most natural one, but, as first noticed by Komissarov (), reconstructing Wv i can increase the robustness of the scheme. However, this is not always the case (see Section .. below) and in our tests we have favored either the first or the second choice according to convenience. Concerning the specific strategy adopted to recover the primitive variables, in our numerical code we have used the third method reported in Section . of Del Zanna et al. (). Alternative methods can be found in Noble et al. (), Rezzolla and Zanotti ().
Finally, there is an important formal change in the transition from purely hydrodynamics systems to genuinely magnetohydrodynamics systems. As already noticed by Londrillo and Del Zanna (), the RMHD equations should not be regarded as a mere extension of the RHD ones, with just a larger number of variables to evolve. Rather, their formal structure is better described in terms of a coupled system of conservation laws (the five equations for the dynamics of the plasma) and a set of Hamilton-Jacobi equations, those for the evolution of the vector potential of the magnetic field (Jin and Xin ). The different mathematical structure of the RMHD equations reflects the existence of the divergence-free property of the magnetic field, which must be ensured at all times during the evolution. Numerically, we have adopted a simplified and well known approach, which consists of augmenting the system () with an additional equation for a scalar field , aimed at propagating away the deviations from ∇ · B = . We therefore need to solve while the fluxes for the evolution of the magnetic field are also changed, namely In the following, we first limit our attention to a few physical systems for which B i = E i = , hence to relativistic hydrodynamics, and then we consider truly magnetohydrodynamics tests with B i = . Table  reports the initial conditions of the two onedimensional Riemann problems that we have considered, and whose wave-patterns at the final time t f = . are shown in Figure  and Figure , respectively. In order to appreciate the differences among the available ADER implementations, we have again solved each problem with
In the first Riemann problem, which was also analyzed by Mignone and Bodo (), two rarefaction waves are produced, separated by a contact discontinuity. It has been solved through a fourth order P  P  scheme, using the Ru-sanov Riemann solver over a uniform grid with  cells. As it is clear from Figure , the ADER-Prim scheme performs significantly better than the ADER-Cons. In particular, the overshoot and undershoot at the tail of the right rarefaction is absent. In general, the results obtained with ADER-Prim are essentially equivalent to those of ADER-Char, namely when the reconstruction in characteristic variables is adopted. This is manifest after looking at the bottom right panel of Figure , where a magnification of the rest mass density at the contact discontinuity is shown. Additional interesting comparisons can be made about the second Riemann problem, which can be found in Radice and Rezzolla (), and which is displayed in Figure . In this case a third order P  P  scheme has been used, again with the Rusanov Riemann solver over a uniform grid with  cells. The right propagating shock has a strong jump  Table 4) with the fourth order ADER-WENO scheme at time t = 0.4. The bottom right panel shows a magnification around the contact discontinuity. Table 4) with the third order ADER-WENO scheme at time t = 0.4. The bottom right panel shows a magnification around the right propagating shock.

Figure 5 Solution of RHD-RP2 (see
in the rest mass density, as it is visible from the bottom right panel of the figure, and the position of the shock front is better captured by the two schemes ADER-Prim and ADER-Char.
It is particularly interesting to address the issue of CPU time comparison among different implementations of ADER, as already done for the Euler equations. The result of such a comparison, performed for the RHD-RP problem, are reported in Table , which should be read in synopsis with Table . Clearly, ADER-Prim is not only more accurate than ADER-Cons, but it is also more efficient. As anticipated, this is in agreement with our expectations, since in the ADER-Prim implementation a single cons-to-prim operation is needed within the cell, rather than at each Gaussian quadrature point and at each spacetime degree of freedom. For other tests, see for instance Section .., the CPU time reduction implied by ADER-Prim is even more evident, but the numbers shown in Table  describe with good fidelity the relative performances of the different ADER in a large number of relativistic tests.

.. RHD Kelvin-Helmholtz instability
In the relativistic regime, the Kelvin-Helmholtz (KH) instability is likely to be responsible for a variety of physical effects, which are encountered in the dynamics of extragalactic relativistic jets (Bodo et al. ; Perucho et al. ; Perucho et al. ). As an academic test, we simulate the linear growth phase of the KH instability in two spatial dimensions, taking the initial conditions from Mignone et al. () (see also Beckwith and Stone  and Radice and Rezzolla ). In particular, the rest-mass density is chosen as with ρ  = . and ρ  = .. Assuming that the shear layer has a velocity v s = . and a characteristic size a = ., the velocity along the x-direction is modulated as It is convenient to add a perturbation in the transverse velocity, i.e.
where η  = . is the amplitude of the perturbation, while σ = . is its length scale. The adiabatic index is γ = / and the pressure is uniform, p = . The problem has been solved over the computational domain [-., .] × [-, ], covered by a uniform mesh with  ×  cells, using the P  P  scheme and the Osher-type numerical flux. Periodic boundary conditions are fixed both in x and in y directions. Figure  shows the results of the calculations: in the left, in the central and in the right panels we have reported the solution obtained with the ADER-Prim, with the ADER-Cons and with the ADER-Char scheme, respectively, while the top and the bottom panels correspond to two different times during the evolution, namely t = . and t = .. Interestingly, two secondary vortices are visible when the reconstruction is performed in primitive and characteristic variables (see left the right panels), but only one is present in the simulation using the reconstruction in conserved variables. In Zanotti and Dumbser () we have already commented about the elusive character of these details in the solution, which depend both on the resolution and on the Riemann solver adopted. Based on our results, we infer that the ADER-Cons scheme is the most diffusive, while ADER-Prim and ADER-Char seem to produce the same level of accuracy in the solution. However, if we look at the CPU times in the two cases, we find that ADER-Prim is a factor . faster than ADER-Cons and a factor  faster than ADER-Char, and therefore should be preferred in all relevant applications of RHD.

.. RMHD Alfvén wave
In Table  of Section .. we have reported the comparison of the convergence rates among three different implementations of ADER for the Euler equations. We believe it is important to verify the convergence of the new ADER-Prim scheme also for the RMHD equations, which indeed admits an exact, smooth unsteady solution, namely the propagation of a circularly polarized Alfvén wave (see Komissarov ; Del Zanna et al.  for a full account). The wave is assumed to propagate along the x direction in a constant density and constant pressure background, say ρ = p = . The magnetic field, on the other hand, is given by where η =  is the amplitude of the wave, B  =  is the uniform magnetic field, k is the wave number, while v A is speed of propagation of the wave. We have solved this problem over the computational domain = [; π] × [; π], using periodic boundary conditions, the Rusanov Riemann solver and the Adams-Bashforth version for the initial guess of the LSDG predictor. We have compared the numerical solution with the analytic one after one period T = L/v A = π/v A . Table  contains the results of our analysis, showing the L  and the L  norms of the error of B y . As apparent from the table, the nominal order of convergence of the new ADER-Prim scheme is recovered with very good accuracy.

.. RMHD Riemann problems
Riemann problems are very relevant also in RMHD, admitting a larger number of waves than in hydrodynamics. The exact solution was provided by Giacomazzo and Rezzolla () already ten years ago, making them very popular as a precise tool to validate numerical codes. We have selected Test  and Test  in Table  of Balsara (), with initial left and right states that are reported in Table . Both the tests have been solved using a fourth order ADER-WENO scheme, the Rusanov Riemann solver and over a uniform grid composed of  cells. The damping factor for the divergence-cleaning procedure is set to κ = . Figure  and Figure  allow to compare the exact solution with the results obtained through the ADER-Prim and the ADER-Cons schemes. Especially for RMHD-RP, the solution obtained with the traditional ADER-Cons scheme is significantly more oscillatory than that produced by ADER-Prim. This is particularly evident in the rest-mass density and in the velocity v x . We have here a good indication that the ADER-Prim scheme behaves better than the ADER-Cons scheme when applied to the equations of special relativistic magnetohydrodynamics.

.. RMHD rotor problem
The relativistic version of the magnetic rotor problem, originally proposed by Balsara and Spicer (), has by now become a standard numerical test in RMHD. It describes the evolution of a high density plasma which, at time t = , rotates rapidly with angular velocity ω and is surrounded by a low density plasma at rest: p = , γ = /. Due to rotation, a sequence of torsional Alfvén waves are launched outside the cylinder, with the net effect of reducing the angular velocity of the rotor. We have solved this problem over a computational domain = [-., .] × [-., .], discretized by  ×  numerical cells and using a fourth order finite volume scheme with the Rusanov Riemann solver. No taper has been applied to the initial conditions, thus producing true discontinuities right at the beginning. Figure  shows the rest-mass density, the thermal pressure, the relativistic Mach number and the magnetic pressure at time t = .. We obtain results which are in good qualitative agreement with those available in the literature (see, for instance, Del Zanna et al. ; Dumbser and Zanotti ; Loubère et al.  and Kim and Balsara ). We emphasize that for this test the reconstruction of the primitive variables v i turns out to be more robust than that achieved through the reconstruction of the products Wv i .

The Baer-Nunziato equations
As a genuinely non-conservative system of hyperbolic equations we consider the Baer-Nunziato model for compressible two-phase flow (see also Baer and Nunziato ; Saurel and Abgrall ; Andrianov and Warnecke ; Schwendeman et al. ; Deledicque and Papalexandris ; Murrone and Guillard ). In the rest of the paper we define the first phase as the solid phase and the second phase as the gas phase. As a result, we will use the subscripts  and s as well as  and g as synonyms. Sticking to Baer and Nunziato (), we prescribe the interface velocity v I and the pressure p I as v I = v  and p I = p  , respectively, although other choices are also possible (Saurel and Abgrall ). With these definitions, the system of Baer-Nunziato equations can be cast in the form prescribed by () after defining the state vector Q as where φ k is the volume fraction of phase k, with the condition that φ  + φ  = . On the other hand, the fluxes f i , the sources S and the non-conservative matrices B i are expressed by   Table 7) with the fourth order ADER-WENO scheme at time t = 0.4. The Rusanov Riemann solver has been used over a 400 cells uniform grid.
where e i is the unit vector pointing in direction i (i ∈ {x, y, z}) and ν and μ are two parameters related to the friction between the phases and to the pressure relaxation. e The equation of state is the so-called stiffened gas equation of state, which is a simple modification of the ideal gas EOS and where π k expresses a reference pressure. For brevity, we have solved this system of equations only for a set of onedimensional Riemann problems, with initial conditions reported in Table . The name of the models, BNRP, BNRP, etc., respects the numeration adopted in Dumbser et al. (). A reference solution is available for these tests, In all the tests, with the exception of BNRP, the ADER-Prim scheme behaves significantly better than the ADER-Cons scheme. On several occasions, such as for v s and v g in BNRP, or for most of the quantities in BNRP, the solution provided through ADER-Cons manifest evident oscillations, which are instead strongly reduced, or even absent, when the ADER-Prim scheme is used. The CPU time overhead implied by ADER-Prim is comparatively limited, and never larger than ∼ %.  Table 7) with the fourth order ADER-WENO scheme at time t = 0.55. The Rusanov Riemann solver has been used over a 400 cells uniform grid.

Extension to discontinuous Galerkin and adaptive mesh refinement
Although we have so far concentrated on the implementation of the new ADER-Prim scheme in the context of finite volume methods, the same idea can be extended to discontinuous Galerkin (DG) schemes as well. Incidentally, we note that the interest of computational astrophysics towards DG methods is increasing (Radice and Rezzolla ; Teukolsky ), and, especially in the relativistic context, they are expected to play a crucial role in the years to come. In a sequence of papers, we have recently developed a class of robust DG schemes which are able to cope even with discontinuous solutions, by incorporating an a posteriori subcell limiter (Dumbser et al. ; Zanotti et al. b; Zanotti et al. a). The whole logic can be briefly summarized as follows. First we assume a discrete representation of the solution, in conserved variables, at any given time t n as in which the polynomials are built using the spatial Lagrange interpolation polynomials already adopted for the WENO reconstruction. The time evolution of the degrees of freedomû n l is then obtained after considering the weak form of the governing Figure 9 Solution of the RMHD rotor problem at time t = 0.4, obtained with the P 0 P 3 scheme on a uniform grid with 300 × 300 cells. Top panels: rest-mass density (left) and thermal pressure (right). Bottom panels: Mach number (left) and magnetic pressure (right).  are then used as initial conditions for the LSDG predictor, i.e.

Table 8 Initial states left (L) and right (R) for the Riemann problems for the Baer-Nunziato equations
In those cells in which the main scheme of Eq. () fails, either because unphysical values of any quantity are encountered, or because strong oscillations appear in the solution which violate the discrete maximum principle, the computation within the troubled cell goes back to the time level t n and it proceeds to a complete re-calculation. In practice, a suitable subgrid is generated just within the troubled cell, and a traditional finite volume scheme is used on the subgrid using an alternative data representation in terms of cell averages defined for each cell of the subgrid.

RMHD blast wave problem
At time t = , the rest-mass density and the pressure are ρ = . and p = , respectively, within a cylinder of radius R = ., while outside the cylinder ρ =  - and p =  ×  - . Moreover, there is a constant magnetic field B  along the x-direction and the plasma is at rest, while a smooth ramp function between r = . and r =  modulates the initial jump between inner and outer values, similarly to Komissarov () and Del Zanna et al. ().
The computational domain is = [-, ] × [-, ], and the problem has been solved over an initial coarse mesh with  ×  elements. During the evolution the mesh is adaptively refined using a refinement factor along each direction r =  and two levels of refinement. A simple Rusanov Riemann solver has been adopted, in combination with the P  P  version of the ADER-DG scheme. On the subgrid we are free to choose any finite volume scheme that we wish, and for this specific test we have found convenient to adopt a second-order TVD scheme. The results for B  = . are shown in Figure , which reports the rest-mass density, the thermal pressure, the Lorentz factor and the magnetic pressure at time t = .. At this time, the solution is composed by an external circular fast shock wave, which is hardly visible in the rest mass density, and a reverse shock wave, which is compressed along the ydirection. The magnetic field is mostly confined between these two waves, as it can be appreciated from the contour plot of the magnetic pressure. The two bottom panels of the figure show the AMR grid (bottom left) and the map of the limiter (bottom right). In the latter we have used the red color to highlight those cells which required the activa- tion of the limiter over the subgrid, while the blue color is for the regular cells. In practice, the limiter is only needed at the inner shock front, while the external shock front is so weak that the limiter is only occasionally activated. These results confirm the ability of the new ADER-Prim scheme to work also in combination with discontinuous Galerkin methods, and with complex systems of equations like RMHD.

Leblanc, Sedov and Noh problem
Here we solve again the classical Euler equations of compressible gas dynamics on a rectangular domain for the Leblanc problem and on a circular domain in the case of  Boscheri and Dumbser (). For the low pressure region that is present in the above test problems, we use p =  - for the Leblanc and the Noh problem. The computational results obtained with very high order ADER-DG P  P  schemes are depicted in Figures ,  and , showing an excellent agreement with the exact solution in all cases, apart from the overshoot in the case of the Leblanc shock tube. We stress that all test problems are extremely severe and therefore clearly demonstrate the robustness of the new approach.

Conclusions
The new version of ADER schemes introduced in Dumbser et al. (b) relies on a local space-time discontinuous Galerkin predictor, which is then used for the computation of high order accurate fluxes and sources. This approach has the advantage over classical Cauchy-Kovalewski based ADER schemes (Toro et al. ; Titarev and Toro ; Toro and Titarev ; Titarev and Toro ; Toro and Titarev ; Dumbser and Munz ; Taube et al. ) that it is in principle applicable to general nonlinear systems of conservation laws. However, for hyperbolic systems in which the conversion from conservative to primi- Figure 13 Results for the Baer-Nunziato Riemann problem BNRP5. The Rusanov Riemann solver has been used over a 300 cells uniform grid. tive variables is not analytic but only available numerically, a large number of such expensive conversions must be performed, namely one for each space-time quadrature point for the integration of the numerical fluxes over the element interfaces and one for each space-time degree of freedom in the local space-time DG predictor.
Motivated by this limitation, we have designed a new version of ADER schemes, valid primarily for finite volume schemes but extendible also to the discontinuous Galerkin finite element framework, in which both the spatial WENO reconstruction and the subsequent local space-time DG predictor act on the primitive variables. In the finite volume context this can be done by performing a double WENO reconstruction for each cell. In the first WENO step, piece-wise polynomials of the conserved variables are computed from the cell averages in the usual way. Then, these reconstruction polynomials are simply evaluated in the cell centers, in order to obtain point val-ues of the conserved variables. After that, a single conversion from the conserved to the primitive variable is needed in each cell. Finally, a second WENO reconstruction acts on these point values and provides piece-wise polynomials of the primitive variables. The local space-time discontinuous Galerkin predictor must then be reformulated in a Figure 15 Solution of the RMHD blast wave at time t = 4.0, obtained with the ADER-DG P 3 P 3 scheme supplemented with the a posteriori second order TVD subcell finite volume limiter. Top panels: rest-mass density (left) and thermal pressure (right). Central panels: Lorentz factor (left) and magnetic pressure (right), with magnetic field lines reported. Bottom panels: AMR grid (left) and limiter map (right) with troubled cells marked in red and regular unlimited cells marked in blue. non-conservative fashion, supplying the time evolution of the reconstructed polynomials for the primitive variables.
For all systems of equations that we have explored, classical Euler, relativistic hydrodynamics (RHD) and magnetohydrodynamics (RMHD) and the Baer-Nunziato equations, we have noticed a significant reduction of spurious oscillations provided by the new reconstruction in primitive variables with respect to traditional reconstruction in conserved variables. This effect is particularly evident for the Baer-Nunziato equations. In the relativistic regime, there is also an improvement in the ability of capturing the position of shock waves (see Figure ). To a large extent, the new primitive formulation provides results that are comparable to reconstruction in characteristic variables.
Moreover, for systems of equations in which the conversion from the conserved to the primitive variables can-not be obtained in closed form, such as for the RHD and RMHD equations, there is an advantage in terms of computational efficiency, with reductions of the CPU time around ∼ %, or more. We have also introduced an additional improvement, namely the implementation of a new initial guess for the LSDG predictor, which is based on an extrapolation in time, similar to Adams-Bashforth-type ODE integrators. This new initial guess is typically faster than those traditionally available, but it is also less robust in the presence of strong shocks.
We predict that the new version of ADER based on primitive variables will become the standard ADER scheme in the relativistic framework. This may become particularly advantageous for high energy astrophysics, in which both high accuracy and high computational efficiency are required.  c See Appendix A of  for a graphical representation. d We note that, since the spacetime is flat and we are using Cartesian coordinates, the covariant and the contravariant components of spatial vectors can be used interchangeably, namely A i = A i , for the generic vector A.