A direct primitive variable recovery scheme for hyperbolic conservative equations: the case of relativistic hydrodynamics

TL;DR

This paper introduces a Primitive Variable Recovery Scheme (PVRS) for hyperbolic conservative equations, especially in relativistic hydrodynamics, which directly computes primitive variables avoiding complex algebraic solutions.

Contribution

The PVRS method generalizes existing techniques by directly recovering primitive variables through the chain rule, simplifying computations in relativistic hydrodynamics.

Findings

Demonstrates convergence on shock-tube problems

Provides graphical error visualization

Avoids algebraic polynomial solutions for primitive variables

Abstract

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.

Tables

Table 1. Table 1: Initial parameters used for the relativistic Sod ( sod, ) shock tube tests described in the article. κ 𝜅 \kappa stands for the polytropic index.

Test	$p_L$	$v_L$	$n_L$	$p_R$	$v_R$	$n_L$	$\kappa$
1	1.0	0.0	1.0	0.1	0.0	0.125	4/3
2	13.33	0.0	10.00	0.1	0.0	1.0	4/3
3	1000	0.0	1.0	0.01	0.0	1.0	5/3

Table 2. 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 𝑡 0.35 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.

	Error					Order of Convergence
Resolution		Test 1	Test 2	Test 3	Smooth		Test 1	Test 2	Test 3	Smooth
$\Delta x_1$		3.83e-3	9.00e-2	1.93e-1	4.74e-4		-	-	-	-
$\Delta x_2$		2.12e-3	5.04e-2	1.60e-1	1.28e-4		0.85	0.84	0.27	1.89
$\Delta x_3$		1.21e-3	2.61e-2	1.21e-1	0.34e-4		0.81	0.95	0.40	1.91
$\Delta x_4$		6.68e-4	1.51e-2	8.03e-2	0.09e-4		0.85	0.79	0.59	1.92
$\Delta x_5$		4.01e-4	1.02e-2	4.56e-2	0.02e-4		0.74	0.57	0.81	2.17
$\Delta x_6$		2.22e-4	5.07e-3	2.62e-2	-		0.83	1.00	1.05	-

Equations

$$\frac{\partial \mathbf{q}}{\partial t}+\frac{\partial \mathbf{f}(\mathbf{q})}{\partial x}=0,$$

$$q_1=D=\frac{n}{\sqrt{1-v^2}}\text{and}{f}_1=v\frac{n}{\sqrt{1-v^2}},$$

$$q_2=S^x=v\frac{e+p}{1-v^2}\text{and}{f}_2=v^2\frac{e+p}{1-v^2}+p,$$

$$q_3=\tau =\frac{e+v^{2}p}{1-v^2}\text{and}{f}_3=v\frac{e+p}{1-v^2}.$$

$$e=nm+\frac{p}{\kappa -1},$$

$$\frac{\partial q_a(u_1,...,u_m)}{\partial t}+\frac{\partial f_a(u_1,...,u_m)}{\partial x}=0.$$

$$\frac{\partial q_a}{\partial u_b}\frac{\partial u_b}{\partial t}+\frac{\partial f_a}{\partial u_c}\frac{\partial u_c}{\partial x}=0,$$

$$\frac{\partial u_a}{\partial t}+M_{ab}\frac{\partial u_b}{\partial x}=0,$$

$$M_{ab}:={\left(\frac{\partial q_c}{\partial u_a}\right)}^{-1}\left(\frac{\partial f_c}{\partial u_b}\right).$$

$$\begin{array}{rl}{u}_a(x_i,t_{n+1})=u_a(x_i,t_n)-\frac{\Delta t}{2\Delta x}{M}_{ab}[& u_b(x_{i+1},t_n)\\ -& u_b(x_{i-1},t_n)].\end{array}$$

$$\begin{split}\frac{\partial q_{a}(x_{i})}{\partial t}=-\frac{1}{\Delta x}\Big{(}&[F_{a}^{\scriptstyle HLL}]_{i+1/2}^{n}\\ -&[F_{a}^{\scriptstyle HLL}]_{i-1/2}^{n}\Big{)},\end{split}$$

$$\begin{split}\frac{\partial u_{a}(x_{i})}{\partial t}=-\frac{\mathcal{A}_{ab}(x_{i},t_{n})}{\Delta x}\Big{(}&[F_{a}^{\scriptstyle HLL}]_{i+1/2}^{n}\\ -&[F_{a}^{\scriptstyle HLL}]_{i-1/2}^{n}\Big{)}.\end{split}$$

$$\begin{split}{u}_{a}(x_{i},t_{n+1})={u}_{a}(x_{i},t_{n})-\frac{\Delta t}{\Delta x}\mathcal{A}_{ab}(x_{i},t_{n})\Big{(}[&F_{b}^{\scriptstyle HLL}]_{i+1/2}^{n}\\ -&[F_{b}^{\scriptstyle HLL}]_{i-1/2}^{n}\Big{)}.\end{split}$$

$$\frac{\partial \mathbf{u}(x_i)}{\partial t}=-\mathbf{L}(\mathbf{u}(x_i)),$$

$$\frac{\partial q_a}{\partial t}+\frac{\partial f_a(q_1,...,q_m)}{\partial x}=0,$$

$$\frac{\partial q_a}{\partial t}+J_{ab}\frac{\partial q_b}{\partial x}=0,$$

$$\begin{split}q_{a}(x_{i},t_{n+1})=q_{a}(x_{i},t_{n})-\frac{\Delta t}{2\Delta x}\big{[}&f_{a}(x_{i+1},t_{n})\\ -&f_{a}(x_{i-1},t_{n})\big{]},\end{split}$$

$$\begin{split}q_{a}(x_{i},t_{n+1})=&q_{a}(x_{i},t_{n})\\ -\frac{\Delta t}{2\Delta x}\big{[}&f_{a}(x_{i+1},t_{n})-f_{a}(x_{i-1},t_{n})\big{]}\\ +\frac{\Delta t}{2\Delta x}\big{[}&\epsilon_{a}^{+}\Delta q_{a}^{+}(x_{i},t_{n})-\epsilon_{a}^{-}\Delta q_{a}^{-}(x_{i},t_{n})\big{]},\end{split}$$

$$\begin{split}q_{a}(x_{i},t_{n+1})=\frac{1}{2}\big{(}&q_{a}(x_{i+1},t_{n})+q_{a}(x_{i-1},t_{n})\big{)}\\ -\frac{\Delta t}{2\Delta x}\big{[}&f_{a}(x_{i+1},t_{n})-f_{a}(x_{i-1},t_{n})\big{]},\end{split}$$

$$\begin{split}\tilde{q}_{a}(x_{i},t_{n}):=q_{a}(x_{i},t_{n})-\frac{\Delta t}{\Delta x}\big{[}&f_{a}(x_{i+1},t_{n})\\ -&f_{a}(x_{i},t_{n})\big{]},\end{split}$$

$$\begin{split}q_{a}(x_{i},t_{n+1})=\frac{1}{2}\bigg{\{}&q_{a}(x_{i},t_{n})+\tilde{q}_{a}(x_{i},t_{n})\\ -\frac{\Delta t}{\Delta x}\big{[}&\tilde{f}_{a}(x_{i},t_{n})-\tilde{f}_{a}(x_{i-1},t_{n})\big{]}\bigg{\}}.\end{split}$$

$$\begin{array}{rl}{q}_a^{\ast }(x_i,t_n)=q_a(x_i,t_n)& +\eta [q_a(x_{i+1},t_n)\\ & -2q_a(x_i,t_n)+q_a(x_{i-1},t_n)],\end{array}$$

$$\eta=\left\{\begin{array}[]{ll}\eta_{0}\leq 1/4,&\textup{if}\quad(\Delta q_{a}^{+})(\Delta q_{a}^{-})<0,\\ 0,&\textup{if}\quad(\Delta q_{a}^{+})(\Delta q_{a}^{-})>0.\end{array}\right.$$

$$\int_{t_{n}}^{t_{n+1}}\int_{C_{i}}\bigg{[}\frac{\partial q_{a}(x,t)}{\partial t}+\frac{\partial f_{a}(x,t)}{\partial x}\bigg{]}\textup{d}x\,\textup{d}t=0.$$

$$\begin{array}{rl}& {\int }_{C_i}(q_a(x,t_{n+1})-q_a(x,t_n))\text{d}x\\ +& {\int }_{t_n}^{t_{n+1}}(f_a(x_{i+1/2},t)-f_a(x_{i-1/2},t))\text{d}t=0.\end{array}$$

$$[Q_a{]}_i^n=\frac{1}{\Delta x}{\int }_{C_i}{q}_a(x,t_n)\text{d}x,$$

$$[F_a{]}_{i\pm 1/2}^n=\frac{1}{\Delta t}{\int }_{t_n}^{t_{n+1}}{f}_a(x_{i\pm 1/2},t)\text{d}t,$$

$$[Q_{a}]_{i}^{n+1}=[Q_{a}]_{i}^{n}-\frac{\Delta t}{\Delta x}\big{(}[F_{a}]_{i+1/2}^{n}-[F_{a}]_{i-1/2}^{n}\big{)}.$$

$$[F_{a}]_{i\pm 1/2}^{n}=\mathcal{F}_{a}\big{(}[Q_{a}]_{i\pm 1}^{n},[Q_{a}]_{i}^{n}\big{)}.$$

$$\begin{split}[Q_{a}]_{i}^{n+1}=[Q_{a}]_{i}^{n}-\frac{\Delta t}{\Delta x}\big{[}&\mathcal{F}_{a}\big{(}[Q_{a}]_{i+1}^{n},[Q_{a}]_{i}^{n}\big{)}\\ -&\mathcal{F}_{a}\big{(}[Q_{a}]_{i-1}^{n},[Q_{a}]_{i}^{n}\big{)}\big{]}.\end{split}$$

$$\frac{\partial q}{\partial t}+\tilde{u}\frac{\partial q}{\partial x}=0.$$

Full text

Special relativity; General theory in fluid dynamics; Conservation laws and constitutive relations; Computational methods in fluid dynamics

I 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 $\boldsymbol{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 $\boldsymbol{u}$. Depending on the nature of the physical problem, the function $\boldsymbol{u}(\boldsymbol{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 $\boldsymbol{u}$ (see e.g. riccardi2008, , for an excellent account on this).

To make things even more complicated, for each particular physical problem it is necessary to have either an analytic solution $\boldsymbol{u}(\boldsymbol{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 $\boldsymbol{u}$ avoiding any middle cumbersome step such as the ones mentioned above. This technique is so general that requires no analytical knowledge whatsoever of $\boldsymbol{u}(\boldsymbol{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 section A, 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 II we construct our “Primitive Variable Recovery Scheme (PVRS)” which can directly obtain the primitive variables from quite a standard numerical procedure. Section V deals with different convergence relativistic Sod (sod, ) shock-tube tests and error estimates are given using a standard L1-norm. Also, the errors are graphically interpreted using the fluctuations of the solution with respect to analytical known values is presented. Finally, in section VI we discuss and conclude our results.

II 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. equation (15)):

[TABLE]

is by implementing Finite Difference and Finite Volume Methods (FDM & FVM) in order to obtain solutions for the conservative charges $\boldsymbol{q}$. In the particular case of relativistic and non-relativistic hydrodynamics, these charges are the linear momentum along the three dimensions $S^{i}$, the energy $\tau$ 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$111Some authors prefer to find the particle mass density $\rho$ rather than the particle number density $n$. For most practical proposes, both variables are related by $\rho=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 by (landau1987, ; weinberg, )). The explicit dependences $\boldsymbol{q}=\boldsymbol{q}(\boldsymbol{u}(x,t))$ and $\boldsymbol{f}=\boldsymbol{f}(\boldsymbol{u}(x,t))$ for 1D flow in the special relativistic case are given by (see e.g. (weinberg, ; landau1987, )):

[TABLE]

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 (tooper1965, ) for a polytropic relativistic gas:

[TABLE]

where $\kappa$ 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 relations (2)-(4), obtaining the inverse function $\boldsymbol{u}=\boldsymbol{u}(\boldsymbol{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 (riccardi2008, ). One way of solving this system is by using a Newton-Raphson method (cf. yousaf2015, ) but this or any other numerical solution to obtain $\boldsymbol{u}(\boldsymbol{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 (15).

III 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.:

[TABLE]

A necessary and sufficient condition for the existence of the solution $u_{1},\ldots,u_{m}$ is that $a=1,\ldots,m$. Now, using the chain rule, the above equation can be written in the following quasilinear form:

[TABLE]

where $\partial\boldsymbol{q}/\partial\boldsymbol{u}$ and $\partial\boldsymbol{f}/\partial\boldsymbol{u}$ are the Jacobian matrixes of the vectors $\boldsymbol{q}$ and $\boldsymbol{f}$ respectively. Multiplying the previous equation by the inverse matrix $(\partial\boldsymbol{q}/\partial\boldsymbol{u})^{-1}$ we get

[TABLE]

where

[TABLE]

If we perform a discretisation of equation (8) using a FDM (see e.g. section A.1 of the appendix), we obtain the following numerical expression:

[TABLE]

No matter how complicated the functional representations of $\boldsymbol{q}(\boldsymbol{u})$ and $\boldsymbol{f}(\boldsymbol{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 $\boldsymbol{u}$ solving equation (8). By doing this, the cumbersome step of recovering $\boldsymbol{u}$ from $\boldsymbol{q}$ at every cell for each time step is not needed anymore.

The discretisation (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 equation (10) is mathematically similar to relation (16) of the appendix (font1994, ) with the substitution of the vector $\boldsymbol{u}$ instead of $\boldsymbol{q}$. Furthermore, equation (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. (leveque2002, )). In other words, this FDM scheme does not work and the approach to follow is to consider flux contributions as in standard FVM.

IV 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 equations (1) can be discretised in the form of relation (48), which can be written in a semi-discrete form as:

[TABLE]

where $\boldsymbol{F}^{\scriptstyle 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:

[TABLE]

where $\mathcal{A}=(\partial q/\partial u)^{-1}$. By applying a forward-difference formula scheme on the left hand side of equation (12), we get

[TABLE]

In equation (13), we take a numerical flux approach as in standard FVM and a finite difference of the time derivative over the primitive variables $\boldsymbol{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 section A.2.4. Furthermore, the characteristic velocities used in the HLL solver which correspond to the the eigenvalues of the Jacobian $\partial\boldsymbol{f}/\partial\boldsymbol{q}$, can be computed either from matrix $M_{ab}$ (9) or from $\partial f_{a}/\partial q_{b}$ since both matrixes are similar (font1994, ). All matrixes and vectors $(M_{ab},\mathcal{A}_{ab},f_{a},q_{a})$ are computed using a piecewise reconstruction $\tilde{u}$ 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 equations (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 (13) arises directly from discretising the hybrid quasilinear equation $\partial\boldsymbol{u}/\partial t+\mathcal{A}\,\partial\boldsymbol{f}/\partial x=0$ –which can be directly obtained by using the chain rule on equation (1), it is impossible to obtain the PVRS discretisation shown in equation (13) using a standard conservative FVM as presented in the appendix, and which satisfies the entropy and Rankine-Hugoniot jump conditions.

By using discretisation (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 $\boldsymbol{u}(x,t)$. In general terms, this procedure works out for any kind of conservative system in which $\boldsymbol{q}(\boldsymbol{u}(x,t))$ and $\boldsymbol{f}(\boldsymbol{u}(x,t))$ are at least given at some initial time.

The time step evolution of the discretisation (13) that we use for our numerical simulations is given by the Method of Lines (MoL):

[TABLE]

where $\boldsymbol{L}(\boldsymbol{u}(x_{i}))$ is the right hand side of equation (13) (see e.g. (lora2013, )), which can be further implemented with a Runge-Kutta integration.

V 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. lora2013, ). This relativistic Sod (sod, ) 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 (13) with an approximate HLL Riemann solver, a minmod limiter for the reconstruction $\tilde{u}$ 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 shown by (lora2013, ). Also, we have estimated the usual $\textup{L}_{1}$-norm error for the following different resolutions: $\Delta x_{1}=1/200$, $\Delta x_{2}=1/400$, $\Delta x_{3}=1/800$, $\Delta x_{4}=1/1600$, $\Delta x_{5}=3200$ and $\Delta x_{6}=1/6400$.

The time-step condition used in this method is different from the commonly used by many authors (cf. delzanna2002, ). 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 elsewhere222For practical purposes, the time step interval can be chosen as a sufficiently smaller number than the corresponding CFL condition (cf. equation (40)).. For the examples presented below, we have chosen a fixed time step for each simulation.

V.0.1 Test 1: Weak relativistic blast wave

The first test corresponds to a lowly relativistic blast wave explosion. The results can be seen in Figure 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.

V.0.2 Test 2: Mildly relativistic blast wave

The second test corresponds to a mildly relativistic blast wave explosion. The results can be seen in Figure 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.

V.0.3 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 Figure 3 are as good as the ones obtained by other codes (cf. delzanna2002, ; lora2015, ).

V.1 Error estimates

We have calculated the error of each test using the traditional L${}_{1}-\text{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}-\text{norm}$ of the $\Delta 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 yousaf2015 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}-\text{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 Figure 4 we show both exact (red dashed-line) and numerical (blue dashed-line) solution vs. the fluctuation $|u_{\text{num}}-u_{\text{exact}}|/u_{\text{exact}}$ at each point (black line), for the density in Test 3 at every resolution. We can see how the Full Width at Half Maximum (FWHM) of the fluctuation 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 L1-norm for which smoothing of the errors can be wrongly interpreted as a positive convergence test.

VI 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 $\boldsymbol{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 $\boldsymbol{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 (e.g. rezzolla, ) 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 equation (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.

We are constructing a GNU Public Licensed (GPL) free software (http://www.gnu.org) called “aztekas” (http://www.aztekas.org) that deals with relativistic hydrodynamics using this PVRS technique.

VII Acknowledgements

This work was supported by DGAPA-UNAM and CONACyT grants (PAPIIT IN112616 and CB-2014-1 #240512). AAO, SM and DO acknowledge economic support from CONACyT (788898, 26344, 255602).

Appendix A 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:

[TABLE]

where the subindex $a$ takes values from 1 to $m$, $\boldsymbol{q}:=\boldsymbol{q}(\boldsymbol{u}(x,t))$ is the vector of conservative charges and $\boldsymbol{f}:=\boldsymbol{f}(\boldsymbol{q}(\boldsymbol{u}(x,t)))$ is the corresponding flux vector along the $x$ axis at a given time $t$. The vector $\boldsymbol{u}$ corresponds to the primitive variables for which its number of entries and functional form of $\boldsymbol{q}(\boldsymbol{u})$ depends on the particular problem to solve333From this point onwards, we are going to use $\boldsymbol{f}(x,t)$ instead of the cumbersome notation $\boldsymbol{f}(\boldsymbol{q}(\boldsymbol{u}(x,t)))$, bearing in mind that both, charges and flux vectors, depend on the primitive variables $\boldsymbol{u}(x,t)$.. As it is shown in section II, the fluxes also have an explicit dependence on the primitive variables but are usually written in terms of the conservative charges.

We can rewrite equation (15) in the quasilinear following form

[TABLE]

where $J_{ab}$ is the Jacobian matrix of $\boldsymbol{f}(\boldsymbol{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 (15) is hyperbolic (see e.g. leveque2002, ).

In the linear cases (when $\boldsymbol{f}$ is a linear function of $\boldsymbol{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 equation (15).

A.1 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 $\Delta 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 equation (15) to find the value of $\boldsymbol{q}$ in all the grid at time $t+\Delta t=:t_{n+1}$ based on its value at $t=:t_{n}$:

[TABLE]

where $x_{i}$ is the $i$-th point on the grid. This is the Forward Time Central Space (FTCS) Euler method (leveque2002, )444In equation (17), the derivative $\partial f_{a}/\partial 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\Delta 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}))/\Delta x$ and $(f_{a}(x_{N-1})-f_{a}(x_{N}))/\Delta x$ respectively.. Unfortunately, discretisation (17) leads to numerical unstable solutions vonneumann1950 . To overcome this problem, many higher order methods have been developed and successfully implemented over time (laney1998, ).

When a second-order finite differences approximation method is used, additional source artificial viscosity terms appear in (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. laney1998, ). The artificial viscosity name was given by von Neumann (vonneumann1950, ) 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 (laney1998, ):

[TABLE]

where $\epsilon_{a}^{\pm}$ are the coefficients of second-order explicit artificial viscosity and $\Delta q_{a}^{\pm}(x_{i},t_{n})=\pm q(x_{i\pm 1},t_{n})\mp q(x_{i},t_{n})$. The choice $\epsilon_{a}^{\pm}=2\Delta x/\Delta t$ simplifies the above equation to:

[TABLE]

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 (maccormack1972, ). It makes a forward-prediction of $\boldsymbol{q}$ and with it, a backward-correction:

[TABLE]

[TABLE]

where $\tilde{\boldsymbol{f}}:=\boldsymbol{f}(\tilde{\boldsymbol{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 Figure 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 (book1975, ) is

[TABLE]

where $\eta$ is the antidiffusion coefficient at space-time points $x_{i}$ and $t_{n}$:

[TABLE]

A.2 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) $C_{i}:=[x_{i-1/2},x_{i+1/2}]$. The integration of (15) over $C_{i}$ between times $t_{n}$ and $t_{n+1}$ yields:

[TABLE]

The integral of $\partial_{t}q$ over time and the integral of $\partial_{x}f$ over space can be solved exactly and so, the next integral form of the previous equation is found:

[TABLE]

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:

[TABLE]

[TABLE]

where $[Q_{a}]_{i}^{n}$555From 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. is the average charge vector of $q$ over $C_{i}$ at time $t_{n}$ and $[F_{a}]_{i\pm 1/2}^{n}$ is the average flux vector across the boundaries of $C_{i}$.

If $\boldsymbol{q}(\boldsymbol{u}(x,t))$ is a smooth function, then the integral (26) agrees with the value of $\boldsymbol{q}$ at the midpoint of the interval to $\mathcal{O}(\Delta x^{2})$ (leveque2002, ).

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 (26) and (27) in equation (25) we obtain the main discretisation for the finite volume scheme usually presented in the literature (cf. leveque2002, ):

[TABLE]

Equation (28) is a numerical recipe of how to compute the mean value $[Q_{a}]_{i}^{n+1}$ 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 (15) except for the choice of the values (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.

A.2.1 Numerical flux

The flux $\boldsymbol{f}$ at (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 $\boldsymbol{f}$ inside the integral is evaluated on the boundaries $x_{i\pm 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 volume666This set of midpoints can be ”safely” considered the ones used in the finite difference mesh mentioned in section A.1..

One way to approximate $[F_{a}]_{i\pm 1/2}^{n}$ is to assume that it can be obtained as a function of the cell average values of $\boldsymbol{q}$ on either side of the interface $x_{i\pm 1/2}$, i.e., $[Q_{a}]_{i\pm 1}^{n}$ and $[Q_{a}]_{i}^{n}$:

[TABLE]

The previous result is expected since in a hyperbolic problem the information of how $\boldsymbol{q}$ change on every cell propagates at a finite characteristic speed (see e.g. laney1998, ; landau1987, ). The function $\mathcal{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 equation (29) into (28) yields:

[TABLE]

The numerical flux function is then determined by the evolution of the solution in each interface. A good first guess for the function $\mathcal{F}_{a}$ is to relate it to the corresponding average flux function of a local (for each cell) Riemann problem (lora2013, ) 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.

A.2.2 Riemann problem

Let us now consider a single conservative equation (i.e. relation (15) with $a=1$ only) in which the flux is written as $f(q)=\tilde{u}q$ where $\tilde{u}$ is a constant value:

[TABLE]

This is the advection equation in which $\tilde{u}$ corresponds to the propagation velocity of $q$. Note that, since $f^{\prime}(q)=\tilde{u}$, equation (31) is also its own quasilinear version.

The function $q(x,t)=\tilde{q}(x-\tilde{u}t)$ satisfies equation (31) for any function $\tilde{q}$. 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 $\boldsymbol{q}(X(t),t)$ and equate the result to zero, i.e.:

[TABLE]

Direct comparison of the above equation with (31), means that that the solution $q(X(t),t)$ is constant all along the ray $X(t)=x_{0}+\tilde{u}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

[TABLE]

where $q_{l}$ and $q_{r}$ are the left and right states respectively, the characteristics $X(t)$ of (31) are then rays with slope $\tilde{u}$ in the $t-x$ plane. With this, the solution can be written as

[TABLE]

Let us consider now a system of $m$ conservative equations (i.e. $a=1,\ 2,\ \ldots,m$ in relation (15)), where $f_{a}=A_{ab}q_{b}$, i.e.:

[TABLE]

where $A_{ab}$ is a constant $m\times m$ matrix and so, the system of conservative equations is linear. If $A_{ab}$ is diagonalisable such that:

[TABLE]

where $R_{ac}$ is the matrix of eigenvectors, with $r_{a}^{p}$ the $p$-th eigenvector, $R_{db}^{-1}$ its inverse and $\Lambda_{cd}=\textup{diag}(\lambda^{1},...,\lambda^{m})$, for $\lambda^{p}$ the $p$-th eigenvalue. If we define the characteristic variables $w_{a}$ as

[TABLE]

it is then possible to rewrite equation (35) as the following system of $m$ advective equations:

[TABLE]

In the case of the Riemann problem, the solution for the $p$-th advective equation is $w_{p}(x,t)=\tilde{w}_{p}(x-\lambda^{p}t,0)$, and the solution $q_{a}(x,t)$ is obtained using the definition of $w_{a}$:

[TABLE]

In this way one can think that $q_{a}$ is a superposition of $m$ waves moving with characteristic velocities $\lambda^{1}$, $\lambda^{2}$, …and $\lambda^{m}$, respectively (laney1998, ).

Another way to see this is by comparing equation (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 $\Delta t$, i.e.,

[TABLE]

The quantity $\lambda\Delta t/\Delta x$ is know as the Courant number and the fulfilment of relation (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.

A.2.3 Godunov scheme

(godunov1959, ) proposed a numerical scheme for solving conservative equations and this method can be used in terms of the Riemann problem as follows. Consider the single equation (31). The algorithm proposed by Godunov has the following recipe:

Compute the average values of the charges $q$ at the time $t=t_{n}$ using equation (26) for $a=1$ only:

[TABLE] 2. 2.

Reconstruct from $[Q]_{i}^{n}$ a polynomial function $\tilde{q}(x,t_{n})$ for every value of $x$. The simplest case for this is to take a constant function:

[TABLE]

In practice (eg. rezzolla2013, ), the value $[Q]_{i}^{n}$ is consider to be $q$ evaluated at the midpoint of the grid cell. 3. 3.

Evolve the hyperbolic equation in an exact or approximate way by a time $\Delta t$ to obtain $\tilde{q}(x,t_{n+1})$. 4. 4.

Take the average of $\tilde{q}(x,t_{n+1})$ over $C_{i}$ to obtain $[Q]_{i}^{n+1}$. 5. 5.

Go back to the first item on the list and iterate until a final time is reached.

As we discuss above, it is impossible to compute exactly the average flux $[F]_{i\pm 1/2}^{n}$ because we do not know the value of $q$ at all times. However, if we consider a Riemann problem in the interface $x_{i\pm 1/2}$ between the grid cells $C_{i}$ and $C_{i\pm 1}$ and apply step 3 of Godunov’s algorithm, we get that $\tilde{q}(x_{i\pm 1/2},t)$ is constant along the curves that satisfies $(x-x_{i\pm 1/2})/t=\text{const.}$

In summary, if we denote by $q^{\downarrow}([Q]_{i}^{n},[Q]_{i\pm 1}^{n})$ the solution to the Riemann problem at $x_{i\pm 1/2}$, the computation of the average fluxes reduces on computing an integral over a constant function (leveque2002, ). In this way, the Godunov’s algorithm can be expressed in terms of average fluxes using the following recipe:

Solve the Riemann problem in the interfaces $x_{i\pm 1/2}$ of the $C_{i}$ grid cell in order to obtain $q^{\downarrow}([Q]_{i}^{n},[Q]_{i\pm 1}^{n})$. 2. 2.

Define $\mathcal{F}([Q]_{i}^{n},[Q]_{i\pm 1}^{n})=f(q^{\downarrow}([Q]_{i}^{n},[Q]_{i\pm 1}^{n}))$. 3. 3.

Apply discretisation (30).

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 equation (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. (miyoshi2005, ),  (lora2015, )).

A.2.4 HLL Riemann solver

One of the most popular approximate Riemann solvers is the one proposed by (HLL1983, ). This Godunov’s base method considers a Riemann problem with constant states $\boldsymbol{q}^{L}$ and $\boldsymbol{q}^{R}$ on each side of the interface in a space-time grid cell $[x_{L},x_{R}]\times[0,T]$ as shown on Figure 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 $\lambda_{R}$ and $\lambda_{L}$ moving across the interface to the right and left respectively. The region delimited by these characteristic rays is denoted by the state $\boldsymbol{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 $\lambda_{L}$ and $\lambda_{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.  (leveque2002, ) or (toro2009, )):

[TABLE]

where

[TABLE]

where $\boldsymbol{f}^{R,L}:=\boldsymbol{f}(\boldsymbol{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 Figure 6 and using the Rankine-Hugoniot jump condition at each characteristic ray ($\lambda_{R,L}$). The final result is that (toro2009, ):

[TABLE]

Notice that $\boldsymbol{f}^{HLL}\neq\boldsymbol{f}(\boldsymbol{q}^{HLL})$. The flux (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 reconstruction $\tilde{\boldsymbol{q}}$ from the average values of $\boldsymbol{q}$ has been made. With this, let $\tilde{{q}_{a}}^{L}(x_{i-1/2},t_{n}):=[Q_{a}]_{i-1}^{n}$ and $\tilde{{q}_{a}}^{R}(x_{i-1/2},t_{n}):=[Q_{a}]_{i}^{n}$ 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:

[TABLE]

The flux through $x_{i+1/2}$ is obtained in an analogous way. So, by substituting these numerical fluxes in the discretisation (30), we finally get the numerical solution for the hyperbolic equation (15) in the finite volume scheme using Godunov’s algorithm with a high resolution(marti2003, ) approximate Riemann HLL solver:

[TABLE]

A simple way of computing $[Q_{a}]_{i}^{n}$ is by considering that this average value match the magnitude of $\boldsymbol{q}$ evaluated at the midpoint of the grid cell $x_{i}$. If $\boldsymbol{q}(x,t)$ is smooth, the error introduced by this approximation is of order $\mathcal{O}(\Delta x^{2})$ (leveque2002, ). In other words:

[TABLE]

Many other HLL-type Riemann solvers have been developed (cf. toro2009, ) and successfully implemented (cf. miyoshi2005, ) but they are beyond the scope of the present article.

A.2.5 Limiters

At first approximation, the reconstruction of $q$ over the grid cell was made considering a constant value $[Q]_{i}^{n}$ 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 by

[TABLE]

where $\sigma_{i}^{n}$ 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 of $\tilde{q}$ in each contiguous grid cells, evaluated at the interfaces $x_{i\pm 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:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Each pair (50-51) and (52-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. Equations (50)-(53) are also valid for each component of the vector $\boldsymbol{q}$ when a coupled system of conservative equations is required.

The usual way of computing $\sigma$ is by considering some useful function based on finite derivatives of $q$ over $C_{i}$. The most used but dissipative reconstruction (also called limiter (leveque2002, )) is the minmod limiter (MM) introduced by (roe1986, ):

[TABLE]

where the function $m_{i\pm 1/2}$ is the average slope (or the finite derivative) of $q$ centred at $x_{i\pm 1/2}$:

[TABLE]

[TABLE]

The minmod function of two values $a$ and $b$ stands for:

[TABLE]

This limiter has been successfully implemented in the case of relativistic hydrodynamics (cf. lora2015, ; rezzolla2013, ).

The monotonic centred limiter MC, proposed by (vanleer1977, ), has less dissipation than minmod near discontinuities, but has been proved to create spurious oscillations in the strong shock cases (lora2015, ). Nevertheless, it produces relatively well damped solutions that capture not too strong shock waves. The slope $\sigma$ is written as in (54) but the MC function has the following form:

[TABLE]

where $c:=(a+b)/2$.

Another piecewise linear reconstruction is the superbee limiter, also proposed by (roe1986, ). This one has a better shock-wave capture than the previous scheme as shown in Figure 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:

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

(colella1984, ) developed a piecewise parabolic reconstruction (PPM), that have been successfully used by many authors in both relativistic (cf. lora2015, ) and non-relativistic hydrodynamics (cf. amik2007, ) but for the purposes of this paper, it will not be considered.