A direct Primitive Variable Recovery Scheme for hyperbolic conservative equations: The case of relativistic hydrodynamics

In this article we develop a Primitive Variable Recovery Scheme (PVRS) to solve any system of coupled differential conservative equations. This method obtains directly the primitive variables applying the chain rule to the time term of the conservative equations. With this, a traditional finite volume method for the flux is applied in order avoid violation of both, the entropy and “Rankine-Hugoniot” jump conditions. The time evolution is then computed using a forward finite difference scheme. This numerical technique evades the recovery of the primitive vector by solving an algebraic system of equations as it is often used and so, it generalises standard techniques to solve these kind of coupled systems. The article is presented bearing in mind special relativistic hydrodynamic numerical schemes with an added pedagogical view in the appendix section in order to easily comprehend the PVRS. We present the convergence of the method for standard shock-tube problems of special relativistic hydrodynamics and a graphical visualisation of the errors using the fluctuations of the numerical values with respect to exact analytic solutions. The PVRS circumvents the sometimes arduous computation that arises from standard numerical methods techniques, which obtain the desired primitive vector solution through an algebraic polynomial of the charges.


Introduction
The use of numerical methods to solve differential equations has constituted a substantial amount of work since the conception of approximate solutions to a given set of equations. In the last few decades, digital computers have been a great help to heavily iterate complicated partial differential equations using extensive numerical, parallel and adaptive mesh techniques in personal computers and large clusters.
Physical laws are often written in a set of conservative differential equations, for which there are many well established convergent numerical techniques to obtain accurate solutions. In spite of this, there is an intermediate step that is often, depending on the nature of the problem, extremely cumbersome to deal with. This appears since the general solution to the problem is obtained as a set of vector charges q at every point or cell on a given domain of space at a particular time in the iteration. However, physical phenomena are described and measured by means of a set of vector primitive variables u. Depending on the nature of the physical problem, the function u(q) may not have an analytic form and so, at every point or cell of the integration space a cumbersome technique requires to be performed for each time step. No matter how fast this routine may be, it introduces an extra computational time that can heavily grow when the space-time resolution increases. In problems of special relativistic hydrodynamics this fact appears and, at each time step, a 10th degree algebraic polynomial has to be solved for a unique given value of each component of the vector u (see e.g., for an excellent account on this, [1]).
To make things even more complicated, for each particular physical problem it is necessary to have either an analytic solution u(q) or a specific numerical technique to obtain it.
In this article we show how it is possible to construct a general numerical iteration method, using a combination of finite differences and finite volume integration techniques for the time and spatial evolutions respectively, to directly find the solutions u avoiding any middle cumbersome step such as the ones mentioned above. This technique is so general that requires no analytical knowledge whatsoever of u(q). The method developed is general and valid to any set of coupled conservative equations. We also show how this method can be applied in the particular case of 1D special relativistic hydrodynamics (1DRHD). For this particular case, we construct convergence tests.
The article is organised as follows. In the Appendix, we briefly mention some (mostly used in relativistic hydrodynamics for shock capturing) of the traditional methods to solve a set of conservative equations. In Section 1 we construct our "Primitive Variable Recovery Scheme (PVRS)" which can directly obtain the primitive variables from quite a standard numerical procedure. Section 4 deals with different convergence relativistic Sod [2] shock-tube tests and error estimates are given using a standard L 1 -norm. Also, the errors are graphically interpreted using the fluctuations of the solution with respect to analytical known values is presented. Finally, in Section 5 we discuss and conclude our results.

Primitive Variable Recovery Scheme (PVRS)
In the appendix we discuss some of the standard techniques for discretising any set of scalar and coupled conservative equations. This is done in order to easy understand the further developments of the article for the less expert reader, and not to interrupt the experienced one with such well known methods. However, we note that in the appendix and in what follows Einstein's summation convention will be used throughout the equations displayed in this article, something that does not usually appear in the literature.
The usual way to solve a system of hyperbolic equations (cf. Eq (15)): is by implementing Finite Difference and Finite Volume Methods (FDM & FVM) in order to obtain solutions for the conservative charges q. In the particular case of relativistic and nonrelativistic hydrodynamics, these charges are the linear momentum along the three dimensions S i , the energy τ and the particle density D. In order to compare the numerical solution with experiments and/or observations, a set of primitive physical measurable variables u needs to be constructed. For this particular case, this primitive variable set is given by by the pressure p, the velocity along three spacial dimensions v i and, the particle number density n. Some authors prefer to find the particle mass density ρ rather than the particle number density n. For most practical proposes, both variables are related by ρ = mn where m is the average mass per particle. In here and in what follows all thermodynamical quantities (pressure p, particle number density n and energy density e and so on, are measured on its proper reference frame following the convention in [3,4]). The explicit dependences q = q(u(x, t)) and f = f(u(x, t)) for 1D flow in the special relativistic case are given by (see e.g. [3,4]): where e is the total (rest plus internal) proper energy per unit volume which can be related with the density and pressure via a state equation e = e(n, p) like the one derived by Tooper [5] for a polytropic relativistic gas: where κ is the polytropic index. In the previous equations and in what follows we choose a system of units in which the velocity of light is set to unity. As we can see from Eqs (2)(3)(4), obtaining the inverse function u = u(q(x, t)) results in quite a completed algebraic problem. In fact, the solution to this problem leads to a system of transcendental algebraic equations that have been deeply studied by Riccardi [1]. One way of solving this system is by using a Newton-Raphson method (cf. [6]) but this or any other numerical solution to obtain u(q(x, t)) will carry an extra error besides the proper numerical error of the FDM or FVM. This procedure also adds a bit of computational processing time since an iteration loop to find the solution needs to be carried out at each cell every time step. In order to avoid this cumbersome task, we show now how it is possible to obtain a direct numerical solution of the primitive variables, which is valid for all conservative equation systems (cf. Eq (15)).

PVRS attempts with finite difference methods
Let us begin by writing the system of m hyperbolic equations showing the explicit dependence on m primitive variables, i.e.: A necessary and sufficient condition for the existence of the solution u 1 , . . ., u m is that a = 1, . . ., m. Now, using the chain rule, the above equation can be written in the following quasilinear form: where @q/@u and @f/@u are the Jacobian matrixes of the vectors q and f respectively. Multiplying the previous equation by the inverse matrix (@q/@u) −1 we get If we perform a discretisation of Eq (8) using a FDM (see e.g. Section Finite differences approach of the appendix), we obtain the following numerical expression: No matter how complicated the functional representations of q(u) and f(u), it is possible (if not by hand, using a Computer Algebra System) to compute the matrix M ab only once before implementing a discretisation scheme. In what follows we show how to implement a numerical scheme to find directly the primitive variables u solving Eq (8). By doing this, the cumbersome step of recovering u from q at every cell for each time step is not needed anymore.
The discretisation Eq (10) is accurate to the first-order and yields quite good results on smooth solutions. When the solution contains a shock wave, the method is stable but not consistent and so no convergent. This could be understood because Eq (10) is mathematically similar to Eq (16) of the appendix [7] with the substitution of the vector u instead of q. Furthermore, Eq (10) is written in a non-conservative form and so, the entropy and Rankine-Hugoniot jump conditions are not satisfied across the shock waves. Due to this fact, the obtained solution converges to a different weak solution as compared to the one obtained by a conservative method (see e.g. [8]). In other words, this FDM scheme does not work and the approach to follow is to consider flux contributions as in standard FVM.

Primitive Variable Recovery Scheme using combined FDM and FVM
We now show how to implement a Primitive Variable Recovery Scheme (PVRS) using both a FDM and a FVM schemes for the time and spatial evolution of the equations. As mentioned at the end of the previous section, the fluxes contribution in the method must not be altered because the entropy and Rankine-Hugoniot jump conditions must be accomplished. To do so, the spatial derivative term must be evolved using a Godunov-type method (e.g. an HLL-type Riemann solver).
In the appendix it is shown that the conservative set of Eq (1) can be discretised in the form of relation Eq (48), which can be written in a semi-discrete form as: where F HLL stands for the HLL-type Riemann solver approximation for the spatial fluxes (see appendix). Using the chain rule on the left hand side of the previous equation, it follows that: ð½F HLL a n iþ1=2 À ½F HLL a n iÀ 1=2 Þ: ð12Þ where A ¼ ð@q=@uÞ À 1 . By applying a forward-difference formula scheme on the left hand side of Eq (12), we get In Eq (13), we take a numerical flux approach as in standard FVM and a finite difference of the time derivative over the primitive variables u. The approximate solution to the Riemann problem, where Rankine-Hugoniot's condition take place, is the same as the one presented in the appendix HLL Riemann solver section. Furthermore, the characteristic velocities used in the HLL solver which correspond to the the eigenvalues of the Jacobian @f/@q, can be computed either from matrix M ab Eq (9) or from @f a /@q b since both matrixes are similar [7]. All matrixes and vectors ðM ab ; A ab ; f a ; q a Þ are computed using a piecewise reconstructionũ of the primitive variables, except for matrix A ab which is evaluated on the midpoint x i of the cell C i .
It is important to note that in Eqs (8) and (13) the second term on the right hand side has an implicit sum over the repeated index a.
Note that, although it seems that the PVRS discretisation Eq (13) arises directly from discretising the hybrid quasilinear equation @u=@t þ A @f =@x ¼ 0 -which can be directly obtained by using the chain rule on Eq (1), it is impossible to obtain the PVRS discretisation shown in Eq (13) using a standard conservative FVM as presented in the appendix, and which satisfies the entropy and Rankine-Hugoniot jump conditions. By using Eq (13) on a numerical code, it would no longer be a concern to recover the primitive variables from the computed conservative charges; they would instead be solved directly! Therefore, it would not be necessary to create a module in the code to obtain the final required solution u(x, t). In general terms, this procedure works out for any kind of conservative system in which q(u(x, t)) and f(u(x, t)) are at least given at some initial time.
The time step evolution of the Eq (13) that we use for our numerical simulations is given by the Method of Lines (MoL): where L(u(x i )) is the right hand side of Eq (13) (see e.g. [9]), which can be further implemented with a Runge-Kutta integration.

Convergence test for PVRS in relativistic hydrodynamics
In this section we are going to show how this new method handles the evolution of a relativistic gas in a particular Riemann problem namely the shock tube (see e.g. [9]). This relativistic Sod [2] shock tube problem is a standard test that any code must fulfil for its validation. It has an exact analytical solution for both special relativistic and non-relativistic hydrodynamics and it is used for comparisons with numerical methods. We calculated the numerical solution using PVRS discretisation Eq (13) with an approximate HLL Riemann solver, a minmod limiter for the reconstructionũ and a 4th order Runge-Kutta Method of Lines (MoL-RK4) for the integration. The problem was solved in the domain [0, 1] with N = 800 identical grid cells. We made three relativistic Sod tests with the initial discontinuity located at x = 0.5 and with initial states shown on Table 1. Furthermore, we compared the numerical results with the exact solution [9]. Also, we have estimated the usual L 1norm error for the following different resolutions: Δx 1 = 1/200, Δx 2 = 1/400, Δx 3 = 1/800, Δx 4 = 1/1600, Δx 5 = 3200 and Δx 6 = 1/6400.
The time-step condition used in this method is different from the commonly used by many authors (cf. [10]). A general CFL-condition applied to this numerical scheme was constructed by us and used in the set of examples presented. The exact condition and its derivation is a subject beyond the scope of this article and will be published elsewhere. For practical purposes, the time step interval can be chosen as a sufficiently smaller number than the corresponding CFL condition (cf. Eq (40)). For the examples presented below, we have chosen a fixed time step for each simulation.

Test 1: Weak relativistic blast wave
The first test corresponds to a lowly relativistic blast wave explosion. The results can be seen in Fig 1, where we compare the numerical solution (points) with the exact solution (lines). It is clear that for both, smooth parts and discontinuities, the numerical solution converges quite well to the exact one.

Test 2: Mildly relativistic blast wave
The second test corresponds to a mildly relativistic blast wave explosion. The results can be seen in Fig 2, where we compare the numerical solution (points) with the exact one (lines). The importance of this test is to see if, with a relative high difference in pressure between both states, the numerical method is capable of solving the density function at the contact discontinuity.

Test 3: Strong relativistic blast wave
Finally, the last test corresponds to a strongly relativistic blast wave explosion. In this case, the density discontinuity is produced by a a 5 orders of magnitude difference between right and left initial detonation pressure, creating a thin shell which numerically is harder to resolve at low resolutions. However, with a relatively small number of cells and a weak variable reconstruction, the results shown on Fig 3 are as good as the ones obtained by other codes (cf. [10,11]).

Error estimates
We have calculated the error of each test using the traditional L 1 -norm value. The convergence order of this test is given by log(error i /error i−1 )/log(1/2), where error j is the L 1 -norm of the Δx j resolution. As we can see from Table 2, the error decreases when the resolution increases, as expected. Also, we obtain first order convergence for all test in at least one  resolution. Additionally, we made an experiment following [6] of a static Gaussian curve in order to estimate the order of convergence of a smooth static profile which, for this case, reaches a convergence value of about 2 in all the tested resolutions for a fixed time step of 0.01. As expected, this means that the important error of the relativistic Sod shock tube test relays on the discontinuities. This is the reason as to why we consider that taking the L 1 -norm is not a clear indicator of the "real" error at the shock waves, so we propose a more relevant useful visual interpretation of this estimation as follows.
In  Table 2. The L 1 -norm for the error in the numerical density for the minmod limiter with different numerical resolutions. The L 1 -norm is computed for all shock-tube and Gaussian tests at time t = 0.35. We also show the order of convergence between different resolutions. Since the error decreases when the resolution increases, the PVRS constructed in the article is stable and converges to the exact solution. The data of this is presented in S1 File.

Error
Order of Convergence   tends to zero as the resolution increases. Working with the fluctuation of the numerical solution about the exact solution is a much better way to easily see the convergence of a numerical method, rather than the traditional L 1 -norm for which smoothing of the errors can be wrongly interpreted as a positive convergence test.

Discussion
In this article we have developed a new numerical algorithm to solve any set of coupled differential conservative equations for which the primitive variable vector u is directly obtained. This is a forward step in numerical methods, since it avoids any intermediate step reconstruction of the primitive variable vector from a previously obtained charge vector q at all points or cells in space at each time. In principle, this means that numerical codes can be written in a more direct form. Also, depending on the nature of the physical problem to solve, the computational time may be reduced with this technique. For practical purposes, we always had in mind special relativistic hydrodynamical problems and for this reason the specific techniques used throughout the article deal with hydrodynamical shock capturing schemes. We demonstrated in the article that the Primitive Variable Recovery Scheme (PVRS) showed good convergence for three shock-tube and one Gaussian tests. Further explorations in other directions, such as a non-static Gaussian test [12] need to be investigated. We will explore more details in future works.
The PVRS presented in this article can be implemented straightforward to any standard hydrodynamical code that already uses HLL Riemann solvers given by Eq (13).
In summary, the PVRS is a numerical maneuver to circumvent the embroiling construction of the primitive vector once the charge vector is obtained from any standard procedure used to solve a set of coupled conservative equations in physical systems.

Traditional approach for numerically solving conservative equations
In this appendix, we deal with traditional well known methods for solving conservative equations. Our intention is to briefly introduce the less versed reader to this topics using Einstein's summation convention.
A system of m conservative equations in one dimension is usually written as: where the subindex a takes values from 1 to m, q ≔ q(u(x, t)) is the vector of conservative charges and f ≔ f(q(u(x, t))) is the corresponding flux vector along the x axis at a given time t. The vector u corresponds to the primitive variables for which its number of entries and functional form of q(u) depends on the particular problem to solve. From this point onwards, we are going to use f(x, t) instead of the cumbersome notation f (q(u(x, t))), bearing in mind that both, charges and flux vectors, depend on the primitive variables u(x, t). As it is shown in section 1, the fluxes also have an explicit dependence on the primitive variables but are usually written in terms of the conservative charges.
We can rewrite Eq (15) in the quasilinear following form where J ab is the Jacobian matrix of f(q). From now on, we use Einstein implicit sum convention over two repeated subindexes contained in the set {a, b, c, d}. If the Jacobian matrix satisfies the conditions of having real eigenvalues and a set of independent eigenvectors, then we say that the system Eq (15) is hyperbolic (see e.g. [8]). In the linear cases (when f is a linear function of q), there exists an analytical solution for (15), but many physical cases give rise to nonlinear conservative systems which are required to be solved using numerical methods.
In the following subsections we briefly mention two of the main numerical methods used to solve 1D conservative systems such as the one written in Eq (15).

Finite differences approach
The finite differences method (FDM) is one of the most useful and simple numerical methods for solving ordinary and partial differential equations. It consists of an approximation of the derivatives of fluxes and charges based on approximations of their values on sufficiently small intervals of space and time. The space is divided in a grid of N centred points spaced by equal length Δx intervals in which the equation is evaluated.
Using Taylor expansions of the involved quantities, it is possible to work out the finite difference form of Eq (15) to find the value of q in all the grid at time t + Δt ≕ t n+1 based on its value at t ≕ t n : where x i is the i-th point on the grid. This is the Forward Time Central Space (FTCS) Euler method [8]. In Eq (17), the derivative @f a /@x at a given time t n was written using a central approximation value given by (f a (x i+1 ) − f a (x i−1 ))/(2Δx). For the left and right boundary points this derivative can be written using a right or left derivative approximation given by: (f a (x 1 ) − f a (x 0 ))/Δx and (f a (x N−1 ) − f a (x N ))/Δx respectively. Unfortunately, Eq (17) leads to numerical unstable solutions [13]. To overcome this problem, many higher order methods have been developed and successfully implemented over time [14]. When a second-order finite differences approximation method is used, additional source artificial viscosity terms appear in Eq (17). Those additional terms are either due to the second derivative approximation in Taylor series or to second differences approximation of the first derivatives (see e.g. [14]). The artificial viscosity name was given by von Neumann [13] since it resembles the viscosity term of the Navier-Stokes equation, but has nothing to do with any physical viscosity.
The general form of the artificial viscosity can be written as [14]: where AE a are the coefficients of second-order explicit artificial viscosity and Dq AE a ðx i ; t n Þ ¼ AEqðx iAE1 ; t n Þ Ç qðx i ; t n Þ. The choice AE a ¼ 2Dx=Dt simplifies the above equation to: which is known as the Lax-Friedrich method. Other second-order-two-step methods, such as the Lax-Wendroff method, have been developed and successfully implemented in many numerical codes. One such favourite two-step method was proposed by MacCormack [15]. It makes a forward-prediction of q and with it, a backward-correction: wheref ≔ f ðqÞ. This method has been proved to be consistent, convergent and stable which is the requirement for any numerical method used in a computational code. Nevertheless, in discontinuities and regions with high pressure gradients, such as regions with shock-waves, this algorithm introduces a dispersive error called the Gibbs phenomenon, which consists on the presence of large spurious oscillations near the finite-jump, such as the example shown in Fig 5. To solve this problem, it is common to apply a corrective diffusion in the regions where the non-physical oscillations appear. The correction presented by Book [16] is where η is the antidiffusion coefficient at space-time points x i and t n :

Finite volume approach
A more natural way of obtaining the discretisation form of (15) is the Finite Volume Method (FVM) which is based on a subdivision of the spatial domain into intervals (also called control volumes or grid cells) The integration of Eq (15) over C i between times t n and t n+1 yields: The integral of @ t q over time and the integral of @ x f over space can be solved exactly and so, the next integral form of the previous equation is found: At this point, both integrals in the previous equation cannot be integrated unless we have the exact form of q, which is precisely the solution to the problem. In order to overcome this, we define each integration as a new numerical vector in the following form: where ½Q a n i is the average charge vector of q over C i at time t n and ½F a n iAE1=2 is the average flux vector across the boundaries of C i . From now on, the square brackets notation [ ] around any numerical function is used to denote the corresponding (space or time) average related to that specific numerical function.
If q(u(x, t)) is a smooth function, then the integral Eq (26) agrees with the value of q at the midpoint of the interval to OðDx 2 Þ [8].
The indexes outside the square bracket do not denote the spatial and time evaluation of the average vector, they are just labels that refer to the time and grid positions of the corresponding numerical values. Substituting the definitions Eqs (26) and (27) in Eq (25) we obtain the main discretisation for the finite volume scheme usually presented in the literature (cf. [8]): ½F a n iþ1=2 À ½F a n iÀ 1=2 : ð28Þ Eq (28) is a numerical recipe of how to compute the mean value ½Q a nþ1 i using the average flux and charge values one time-step backwards for each grid cell C i . This discretisation has the same exact form as Eq (15) except for the choice of the values Eqs (26) and (27).
The advantage of this method over any finite difference scheme is that the conservative nature of the system is preserved, even across strong discontinuities such as shock waves. This is the reason as to why a finite volume scheme is often used when dealing with the physics of high energy flows where discontinuities may appear.

Numerical flux
The flux f at Eq (27) depends on the value of q at every time. This is why it is impossible to integrate the average flux. Somehow, we have to find a good approximation for this integral. Moreover, the flux f inside the integral is evaluated on the boundaries x i±1/2 of the grid cell which, numerically speaking, has no sense because we can only approximate the values of the average charges on the midpoint of the finite volume. This set of midpoints can be "safely" considered the ones used in the finite difference mesh mentioned in the Finite differences approach section.
One way to approximate ½F a n iAE1=2 is to assume that it can be obtained as a function of the cell average values of q on either side of the interface x i±1/2 , i.e., ½Q a n iAE1 and ½Q a n i : ½F a n iAE1=2 ¼ F a ð½Q a n iAE1 ; ½Q a n i Þ: The previous result is expected since in a hyperbolic problem the information of how q change on every cell propagates at a finite characteristic speed (see e.g. [3,14]). The function F a can be thought as a numerical flux function for which its functional form will depend on the problem or the particular numerical scheme used to solve it. Substitution of Eq (29) into Eq (28) yields: Dx ½F a ð½Q a n iþ1 ; ½Q a n i Þ À F a ð½Q a n iÀ 1 ; ½Q a n i Þ: ð30Þ The numerical flux function is then determined by the evolution of the solution in each interface. A good first guess for the function F a is to relate it to the corresponding average flux function of a local (for each cell) Riemann problem [9] with two constant states on each side of the boundary. In order to obtain an accurate numerical flux function, is important to study the behaviour of the solution based on the form and properties of the governing equation at these particular initial conditions.

Riemann problem
Let us now consider a single conservative equation (i.e. Eq (15) with a = 1 only) in which the flux is written as f ðqÞ ¼ũq whereũ is a constant value: This is the advection equation in whichũ corresponds to the propagation velocity of q. Note that, since f 0 ðqÞ ¼ũ, Eq (31) is also its own quasilinear version. The function qðx; tÞ ¼qðx ÀũtÞ satisfies Eq (31) for any functionq. However, it is more useful for us to describe the problem observing the behaviour of the solution q along characteristic curves in the t − x plane. To do so, we perform the time derivative of q(X(t), t) and equate the result to zero, i.e.: d dt qðXðtÞ; tÞ ¼ @q @t þ X 0 ðtÞ @q @x ¼ 0: Direct comparison of the above equation with Eq (31), means that that the solution q(X(t), t) is constant all along the ray XðtÞ ¼ x 0 þũt, where x 0 is some initial value. In the most general case, the set of all rays X(t) are called the characteristics of the equation.
If we consider the particular case in which the initial conditions of the problem consists on two constant states where q l and q r are the left and right states respectively, the characteristics X(t) of (31) are then rays with slopeũ in the t − x plane. With this, the solution can be written as Let us consider now a system of m conservative equations (i.e. a = 1, 2, . . ., m in Eq (15)), where f a = A ab q b , i.e.: where A ab is a constant m × m matrix and so, the system of conservative equations is linear. If A ab is diagonalisable such that: where R ac is the matrix of eigenvectors, with r p a the p-th eigenvector, R À 1 db its inverse and L cd ¼ diagðl 1 ; . . . ; l m Þ, for λ p the p-th eigenvalue. If we define the characteristic variables w a as it is then possible to rewrite Eq (35) as the following system of m advective equations: In the case of the Riemann problem, the solution for the p-th advective equation is w p ðx; tÞ ¼w p ðx À l p t; 0Þ, and the solution q a (x, t) is obtained using the definition of w a : q a ðx; tÞ ¼ R abwb ðx; tÞ: In this way one can think that q a is a superposition of m waves moving with characteristic velocities λ 1 , λ 2 , . . . and λ m , respectively [14].
Another way to see this is by comparing Eq (35) with the time derivative of q(X(t), t) in (32). From this, it follows that the characteristics are curves for which their corresponding slopes are exactly the eigenvalues of the matrix A ab .
In order to obtain a real contribution of one of these waves to the evolution of a contiguous grid cell, the size of the control volume must be larger than the distance travelled by the wave, moving at its characteristic velocity, at a certain fixed time Δt, i.e., The quantity λΔt/Δx is know as the Courant number and the fulfilment of Eq (40) is called Courant-Friedrich-Levy (CFL) condition. This is a convergence requirement for several numerical methods that solve conservative equations. The Riemann problem discussed in this subsection, is used to accurately estimate the value of the numerical fluxes at the boundaries of two contiguous grid cells as will be seen in the following section. The problem with applying Godunov's scheme on non-linear systems and considering wave propagation of characteristic waves on all interfaces, is that the characteristic velocities are not constant at all times and also they change values at different grid cells. For the case of a quasilinear system such as the one of Eq (16), an approximation has to be made. Many methods for obtaining an approximate Riemann solution have been developed and successfully implemented in classical and relativistic magnetohydrodynamic codes (see e.g. [11,19]).

HLL Riemann solver
One of the most popular approximate Riemann solvers is the called HLL solver [20]. This Godunov's base method considers a Riemann problem with constant states q L and q R on each side of the interface in a space-time grid cell [x L , x R ] × [0, T] as shown on Fig 6. Instead of following the solution of all the characteristic variables along their own characteristic velocities, the idea of the HLL approximation consists on considering the larger eigenvalues λ R and λ L moving across the interface to the right and left respectively. The region delimited by these characteristic rays is denoted by the state q HLL .
Note that, since we are working with a system of m conservative equations, 2m characteristic rays will emerge from each interface. The values λ L and λ R are to be chosen taking into account all 2m characteristic velocities.
The approximate solution to the Riemann problem derived by this scheme has the following form (see e.g. [8] or [21]): where where f R,L ≔ f(q R,L ). One can work out the approximate solution to the flux through the interface by integrating the hyperbolic equation over the space-time domain outlined in Fig 6 and using the Rankine-Hugoniot jump condition at each characteristic ray (λ R,L ). The final result is that [21]: Notice that f HLL 6 ¼ f(q HLL ). The flux Eq (45) can be used along with the Godunov scheme to solve the local Riemann problem of to contiguous grid cells.
Let us now consider the boundary x i−1/2 between two control volumes C i and C i−1 and suppose that a constant reconstructionq from the average values of q has been made. With this, letq a L ðx iÀ 1=2 ; t n Þ ≔ ½Q a n iÀ 1 andq a R ðx iÀ 1=2 ; t n Þ ≔ ½Q a n i to be the reconstruction points that lay at the interface x i−1/2 . Note that these values are going to be different if a polynomial reconstruction is made. With this, we can write the numerical flux at x i−1/2 used in the Godunov scheme in the following form: The flux through x i+1/2 is obtained in an analogous way. So, by substituting these numerical fluxes in the discretisation Eq (30), we finally get the numerical solution for the hyperbolic Eq (15) in the finite volume scheme using Godunov's algorithm with a high resolution [22] approximate Riemann HLL solver: ½F HLL a n iþ1=2 À ½F HLL a n iÀ 1=2 : ð47Þ A simple way of computing ½Q a n i is by considering that this average value match the magnitude of q evaluated at the midpoint of the grid cell x i . If q(x, t) is smooth, the error introduced by this approximation is of order OðDx 2 Þ [8]. In other words: q a ðx i ; t nþ1 Þ ¼ q a ðx i ; t n Þ À Dt Dx ½F HLL a n iþ1=2 À ½F HLL a n iÀ 1=2 : ð48Þ Many other HLL-type Riemann solvers have been developed (cf. [21]) and successfully implemented (cf. [19]) but they are beyond the scope of the present article.

Limiters
At first approximation, the reconstruction of q over the grid cell was made considering a constant value ½Q n i which is taken as the midpoint value of q of the corresponding control volume C i . A better way of improving the precision of the above procedure is by considering a piecewise polynomial approximation for this variable.
In the linear case, the reconstruction of q over C i is given bỹ where s n i is the slope of the linear reconstruction. To use the limiters together with a HLL-type Riemann solver, all we need to consider are those points ofq in each contiguous grid cells, evaluated at the interfaces x i±1/2 . In this respect, it is not important to do a complete reconstruction of q. The knowledge of q at the boundaries is sufficient for this approximation, and so the values required to effectively evolve the solution of the hyperbolic equation over the grid cell C i are:q L ðx iÀ 1=2 ; t n Þ ¼ qðx iÀ 1 ; t n Þ þ 1 2 s n iÀ 1 Dx; ð50Þ q R ðx iÀ 1=2 ; t n Þ ¼ qðx i ; t n Þ À 1 2 s n i Dx; ð51Þ q L ðx iþ1=2 ; t n Þ ¼ qðx i ; t n Þ þ 1 2 s n i Dx; ð52Þ Each pair Eqs (50 and 51) and Eqs (52 and 53), constitute a Riemann problem to be solved at the interface x i−1/2 and x i+1/2 , respectively. The polynomial reconstruction are useful to accurate capture discontinuities such as shock-waves. Eqs (50)-(53) are also valid for each component of the vector q when a coupled system of conservative equations is required. The usual way of computing σ is by considering some useful function based on finite derivatives of q over C i . The most used but dissipative reconstruction (also called limiter [8]) is the minmod limiter (MM) introduced in [23]: where the function m i±1/2 is the average slope (or the finite derivative) of q centred at x i±1/2 : The minmod function of two values a and b stands for: minmodða; bÞ ≔ 0; if ab 0; a; if jaj < jbj; b; if jbj < jaj: This limiter has been successfully implemented in the case of relativistic hydrodynamics (cf. [11,18]).
The monotonic centred limiter MC, proposed by van Leer [24], has less dissipation than minmod near discontinuities, but has been proved to create spurious oscillations in the strong shock cases [11]. Nevertheless, it produces relatively well damped solutions that capture not too strong shock waves. The slope σ is written as in Eq (54)  where c ≔ (a + b)/2. Another piecewise linear reconstruction is the superbee limiter, also proposed by Roe in 1986 [23]. This one has a better shock-wave capture than the previous scheme as shown in Fig  7, where comparisons of the superbee limiter with the previous ones and with the piecewise constant reconstruction (godunov) is made. For this slope, the function is slightly more complicated than the previous ones and is given by: b; if jaj < jbj: Colella in 1984 [25] developed a piecewise parabolic reconstruction (PPM), that have been successfully used by many authors in both relativistic [11] and non-relativistic hydrodynamics (cf. [26]) but for the purposes of this paper, it will not be considered.
Supporting information S1 File. Numerical vs. exact data. In the file "S1_File.tar.gz", we attached the relevant data of the comparison between numerical simulations and exact solutions used to obtain the L 1norm. (GZ)