# An Efficient Implementation of Flux Formulae in Multidimensional Relativistic Hydrodynamical Codes

## Abstract

We derive and analyze a simplified formulation of the numerical viscosity terms appearing in the expression of the numerical fluxes associated to several High-Resolution Shock-Capturing schemes. After some algebraic pre-processing, we give explicit expressions for the numerical viscosity terms of two of the most widely used flux formulae, which implementation saves computational time in multidimensional simulations of relativistic flows. Additionally, such treatment explicitely cancells and factorizes a number of terms helping to amortiguate the growing of round-off errors. We have checked the performance of our formulation running a 3D relativistic hydrodynamical code to solve a standard test-bed problem and found that the improvement in efficiency is of high practical interest in numerical simulations of relativistic flows in Astrophysics.

thanks: This work has been supported by the Spanish DGES (grant PB97-1432). The calculations were carried out on a SGI Origin 2000, at the Centre de Informàtica de la Universitat de València.  

PACS 47.11.+j, 47.75.+f, 95.30.Lz

Key words: Non-linear Systems of Conservation Laws. High Resolution Shock Capturing methods. Special Relativistic Hydrodynamics. General Relativistic Hydrodynamics.

The numerical study of the evolution of multidimensional relativistic flows turns out to be a topic of crucial interest in, at least, two different scientific fields: Nuclear Physics (studies of the properties of the equation of state for nuclear matter via comparison of simulations and experiments of heavy ion collisions) and Relativistic Astrophysics. The field of Numerical Relativistic Astrophysics is recently undergoing an extraordinary developement after the important efforts of people working in building up robust codes able to describe many different astrophysical scenarios, such that relativistic jets in quasars and microquasars, accretion onto compact objects, collision of compact objects, stellar core collapse and recent models of Gamma-Ray bursts (see, e.g., the recent review in  and references therein). Thus, the improvement in the efficiency of multidimensional hydro-codes becomes a necessity.

It is well known the performance of modern high-resolution shock-capturing techniques (HRSC) in simulations of complex classical flows. Most of the HRSC methods are based on the solution of local Riemann problems (i.e., initial value problems with discontiuous initial data) and since 1991  several Riemann Solvers or Flux Formulae have been specifically designed in relativistic fluid dynamics (see, e.g.,  for a review on Riemann solvers in Relativistic Astrophysics). In addition, in a recent paper  we showed the way for applying special relativistic Riemann solvers in General Relativistic Hydrodynamics, hence any future new Riemann solver, exhaustively analyzed in Special Relativistic Hydrodynamics (SRH), could be applied to get the numerical solution of local Riemann problems in General Relativistic Hydrodynamics. Consequently, the interest of the results we obtain in this note goes beyond the domain of SRH and can be easily extended to General Relativistic Hydrodynamics.

For consistency, we start by summarizing the basics of the HRSC techniques. A system of conservation laws  is

$$\frac{𝐮}{t}+\underset{i=1}{\overset{3}{}}\frac{𝐟^i(𝐮)}{x_i}=0$$
(1)   
where u$`\mathrm{}^d`$ is the vector of unknowns and $`𝐟^i`$(u) is the flux in the $`i`$-direction. In the above system (1) we can define a $`5\times 5`$-Jacobian matrix $`^i`$(u) associated to the flux in the $`i`$-direction as:

$$^i=\frac{𝐟^i(𝐮)}{𝐮}.$$
(2)   
The system is said to be hyperbolic if the Jacobian matrices have real eigenvalues.

The main ingredients of a HRSC algorithm are:

i) A finite discretization of the equations in conservation form (1). Using a method of lines, this discretization reads:

$$\frac{d𝐮_{i,j,k}(t)}{dt}=\frac{\widehat{𝐟}_{i+\frac{1}{2},j,k}\widehat{𝐟}_{i\frac{1}{2},j,k}}{\mathrm{\Delta }x}\frac{\widehat{𝐠}_{i,j+\frac{1}{2},k}\widehat{𝐠}_{i,j\frac{1}{2},k}}{\mathrm{\Delta }y}\frac{\widehat{𝐡}_{i,j,k+\frac{1}{2}}\widehat{𝐡}_{i,j,k\frac{1}{2}}}{\mathrm{\Delta }z}$$
(3)   
where subscripts $`i,j,k`$ are related, respectively, with $`x`$, $`y`$ and $`z`$-discretizations, and refer to cell-centered quantities. The cell width, in the three coordinate directions are, respectively, $`\mathrm{\Delta }x`$, $`\mathrm{\Delta }y`$ and $`\mathrm{\Delta }z`$.

ii) Quantities $`\widehat{𝐟}_{i+\frac{1}{2},j,k}`$, $`\widehat{𝐠}_{i,j+\frac{1}{2},k}`$ and $`\widehat{𝐡}_{i,j,k+\frac{1}{2}}`$ are called the numerical fluxes at the cell interfaces. These numerical fluxes are, in general, functions of the states of the system at each side of the cell interface. Some HRSC methods derive expressions for the numerical fluxes by giving a consistent flux formulae or solving local Riemann problems, with an exact  or approximate Riemann solver, after a cell reconstruction procedure that gives the state at both sides of the interface, denoted by $`L`$ (left state) and $`R`$ (right state). Several monotonic cell reconstruction prescriptions have been given in the scientific literature to achieve different orders of spatial accuracy .

For clarity, from now on we will omit the indexes relative to the grid and restrict our study to the $`x_1`$-splitting of the above system (1), assuming that the vector of unknowns satisfies $`𝐮=𝐮(x_1,t)`$.

We have focussed our analysis to some of the most popular HRSC algorithms, and analyzed their expressions for the numerical fluxes. Hence, the sample considered is: HLLE , Roe , Marquina (M) , and a modified Marquina’s flux formula (MM) . The above selection gathers the most fundamental differences among the large sample of HRSC flux formulae. HLLE is the simplest one, it does not need the full spectral decomposition of the Jacobian matrices. Roe’s solver linearizes the information contained in the spectral decomposition into an average state. Marquina’s (and its modified version) flux formula considers the information coming from each side of a given interface (it is not a Riemann solver) and, in some astrophysical applications, has produced better results in modelling complex flows.

After some algebraic work, all these flux formulae can be cast into the following general form:

$$\widehat{𝐟}(𝐮^L,𝐮^R)=\frac{1}{2}\left((+\stackrel{~}{}^L)𝐟^L+(\stackrel{~}{}^R)𝐟^R+(𝒬^L𝐮^L𝒬^R𝐮^R)\right)$$
(4)   
where $`𝐟^{L,R}`$ stands for the fluxes evaluated at the states $`𝐮^{L,R}`$ and $``$ is the unit matrix. Following Harten , the $`𝒬^,`$ terms in the above equation will be referred as the numerical viscosity matrix.

Matrices $`\stackrel{~}{}^{L,R}`$ and $`𝒬^{L,R}`$ can be expressed as linear combinations of the projectors onto each eigenspace, i.e., the direct product of the corresponding left and right eigenvectors $`𝐥_p,𝐫_p`$ associated to the p-th characteristic field (p=1,…,d),

$$\stackrel{~}{}^{L,R}=\underset{p=1}{\overset{d}{}}b_p𝐥_p^{L,R}𝐫_p^{L,R}$$
(5)   
$$𝒬^{L,R}=\underset{p=1}{\overset{d}{}}c_p𝐥_p^{L,R}𝐫_p^{L,R}$$
(6)   
where superscripts $`L,R`$ indicate that the eigenvectors are evaluated at the state $`𝐮^{L,R}`$. The values of the coefficients $`b_p`$ and $`c_p`$ appearing in the above definitions of matrices $`\stackrel{~}{}^{L,R}`$ and $`𝒬^{L,R}`$ depend on the eigenvalues $`\lambda _p`$ as shown in Table I, for the four flux formulae analyzed.

Several comments concerning Table I are in order:

i) If we take into account the orthonormality relations between the right and left eigenvectors

$$\underset{p=1}{\overset{d}{}}𝐥_p𝐫_p=$$
(7)   
and the fact that the coefficients $`b_p`$ and $`c_p`$ are, in the case of HLLE, independents of $`p`$, then matrices $`\stackrel{~}{}^{L,R}`$ and $`𝒬^{L,R}`$ are, trivially, the unit matrix multiplied by the corresponding factors.

ii) For HLLE’s and Roe’s flux formulae their corresponding matrices $`\stackrel{~}{}^{L,R}`$ and $`𝒬^{L,R}`$ satisfy the relations: $`\stackrel{~}{}^L`$ = $`\stackrel{~}{}^R`$, $`𝒬^L`$ = $`𝒬^R`$

iii) As it is well known, the knowledge of the spectral decomposition of the Jacobian matrices is a basic ingredient to build up Riemann solvers or many flux formulae. Nevertheless, while HLLE’s flux formula only needs the values of the maximum and minimum speeds of propagation of the signals, Roe’s and Marquina’s flux formulae need explicitly the full knowledge of the spectral decomposition, including right and left eigenvectors.

The system governing the evolution of multidimensional relativistic perfect fluids can be written in Cartesian coordinates in the form (1), with $`d=5`$, where, in units such that the speed of light $`c=1`$, the vector of unknowns $`𝐮`$ is given by

$$𝐮=(D,S^1,S^2,S^3,\tau )^T,$$
(8)   
the fluxes are defined by

$$𝐟^i=(Dv^i,S^1v^i+p\delta ^{1i},S^2v^i+p\delta ^{2i},S^3v^i+p\delta ^{3i},S^iDv^i)^T$$
(9)   
where $`D(=\rho W)`$ is the rest mass density, $`S^j(=\rho hW^2v^j)`$ is the j-component of the momentum density, and $`\tau (=\rho hW^2p\rho W)`$ is the energy density, $`W=(1\mathrm{v}^2)^{1/2}`$ is the Lorentz factor, $`\rho `$ is the rest–mass density, $`p`$ the pressure and $`h`$ the specific enthalpy given by $`h=1+\epsilon +p/\rho `$ with $`\epsilon `$ being the specific internal energy. The system of SRH is closed with an equation of state $`p=p(\rho ,\epsilon )`$ from which the local sound speed, $`c_s`$, can be obtained

$$hc_s^2=\frac{p}{\rho }+(p/\rho ^2)\frac{p}{ϵ},$$
(10)   
In a previous paper  we derived the explicit analytical expressions for the full (right and left) spectral decomposition. We denote the five characteristic fields by $`p=1,\mathrm{},5,0,0,0,+`$, in the standard ordenation.

A very worthy simplification on the calculation of matrices $`𝒬`$ arises when some eigenvalue is degenerate, i.e., when the system is not strictly hyperbolic. In SRH, like in multidimensional Newtonian hydrodynamics, there is a linearly degenerate field, $`p=0`$, such that the corresponding eigenvalue $`\lambda _0`$ is triple (the system in the $`j`$-direction splitting is not strictly hyperbolic, although the set of eigenvectors is complete). According to equation (6), and using the orthonormality relations between the right and left eigenvectors

$$\underset{k=1}{\overset{3}{}}r_{0,k}^ml_{0,k}^n=\delta ^{mn}r_+^ml_+^nr_{}^ml_{}^n$$
(11)   
where $`m,n=1,\mathrm{}5`$ denote the components of the 5-vector, it is possible to eliminate the three eigenvectors associated to the degenerate field and to write down the following simplified expression (omitting $`L`$,$`R`$ superscripts)

$$𝒬^{mn}=c_0\delta ^{mn}+(c_+c_0)r_+^ml_+^n+(c_{}c_0)r_{}^ml_{}^n.$$
(12)   
Notice that only $`r_\pm `$ and $`l_\pm `$ are needed to evaluate the numerical viscosity. The same procedure can be applied to any system of conservation laws where one of the eigenvectors has multiple degeneracy, because orthogonality relations always allow us to skip the explicit dependence on one of the vector subspaces of the spectral decomposition. In particular, it could be of great interest in the case of the equations of relativistic magnetohydrodynamics where, in the ansatz of a directional splitting, similar degeneracy arises in the structure of the characteristic fields associated to each one of the fluxes.

The explicit formulae for the numerical viscosity term corresponding to the system of equations of special relativistic hydrodynamics are:

HLLE’s flux formulae. Since the numerical viscosity matrix is proportional to the identity, the application of the above recipes is obvious.

Roe’s flux formulae. The numerical flux across some given interface can be written

$$\widehat{𝐟}(𝐮^L,𝐮^R)=\frac{1}{2}[𝐟^L+𝐟^R+𝐪]$$
(13)   
$`𝐪`$ being the five–vector calculated from the corresponding numerical viscosity matrices of Table I:

$$𝐪=𝒬(𝐮^L𝐮^R)𝒬\mathrm{\Delta }𝐮$$
(14)   
In Roe’s Riemann solver the quantities relative to the spectral decomposition are evaluated using the corresponding Roe-averages of the left and right states, denoted by $`\stackrel{~}{𝐮}`$ (see , for the Newtonian case and  for the relativistic case). In practice, other particular averaging (e.g., arithmetic means) have also been used. Note that in the following expressions (15) all quantities are evaluated using Roe’s average, except for $`\mathrm{\Delta }u_m`$. After some algebra, the viscosity vector associated to the numerical flux in the $`j`$-direction is

$`q_1`$ $`=`$ $`\lambda _0\mathrm{\Delta }u_1+\chi _a`$   
$`q_2`$ $`=`$ $`\lambda _0\mathrm{\Delta }u_2+hW(v_x\chi _a+\chi _b\delta _{jx})`$   
$`q_3`$ $`=`$ $`\lambda _0\mathrm{\Delta }u_3+hW(v_y\chi _a+\chi _b\delta _{jy})`$   
$`q_4`$ $`=`$ $`\lambda _0\mathrm{\Delta }u_4+hW(v_z\chi _a+\chi _b\delta _{jz})`$   
$`q_5`$ $`=`$ $`\lambda _0\mathrm{\Delta }u_5+hW(\chi _a+v_j\chi _b)\chi _a`$ (15)   
where

$`\chi _a`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{5}{}}}\left[(\lambda _+\lambda _0)l_+^m+(\lambda _{}\lambda _0)l_{}^m\right]\mathrm{\Delta }u_m`$ (16)   
$`\chi _b`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{5}{}}}\left[(\lambda _+\lambda _0)V_+^jl_+^m+(\lambda _{}\lambda _0)V_{}^jl_{}^m\right]\mathrm{\Delta }u_m`$ (17)   
$`V_\pm ^j`$ $`=`$ $`{\displaystyle \frac{\lambda _\pm v^j}{1v^j\lambda _\pm }}`$ (18)   
M and MM- flux formulae. The numerical flux across a given interface can be written like equation (13) with

$$𝐪=𝐪^L𝐪^R$$
(19)   
$$𝐪^{L,R}=𝒬^{L,R}𝐮^{L,R}$$
(20)   
Omitting the superscripts $`L,R`$ and taken into account the expressions in Table I for MM and the results in , the viscosity vector in the $`x`$-splitting is:

$`q_1^{L,R}`$ $`=`$ $`{\displaystyle \frac{h^2}{\mathrm{\Delta }}}\left\{M\left[𝒜_{}\mathrm{\Omega }_+𝒜_+\mathrm{\Omega }_{}\right]+p(c_+\mathrm{}_+c_{}\mathrm{}_{})\right\}+`$   
$`c_0p{\displaystyle \frac{W}{h}}\left\{{\displaystyle \frac{𝒦}{𝒦1}}+{\displaystyle \frac{v_y^2+v_z^2}{1v_x^2}}\right\}`$   
$`q_2^{L,R}`$ $`=`$ $`{\displaystyle \frac{h^2W}{\mathrm{\Delta }}}\left\{M𝒜_+𝒜_{}\left[\mathrm{\Omega }_+\lambda _+\mathrm{\Omega }_{}\lambda _{}\right]+p(c_+\lambda _+𝒜_+\mathrm{}_+c_{}\lambda _{}𝒜_{}\mathrm{}_{})\right\}+`$   
$`c_0pW^2v_x\left\{{\displaystyle \frac{1}{𝒦1}}+2{\displaystyle \frac{v_y^2+v_z^2}{1v_x^2}}\right\}`$   
$`q_3^{L,R}`$ $`=`$ $`{\displaystyle \frac{h^2W}{\mathrm{\Delta }}}v_y\left\{M\left[\mathrm{\Omega }_+𝒜_{}\mathrm{\Omega }_{}𝒜_+\right]+p(c_+\mathrm{}_+c_{}\mathrm{}_{})\right\}+`$   
$`c_0p\left\{{\displaystyle \frac{W^2}{𝒦1}}+{\displaystyle \frac{1+2W^2(v_y^2+v_z^2)}{1v_x^2}}\right\}`$   
$`q_4^{L,R}`$ $`=`$ $`{\displaystyle \frac{h^2W}{\mathrm{\Delta }}}v_z\left\{M\left[\mathrm{\Omega }_+𝒜_{}\mathrm{\Omega }_{}𝒜_+\right]+p(c_+\mathrm{}_+c_{}\mathrm{}_{})\right\}+`$   
$`c_0p\left\{{\displaystyle \frac{W^2}{𝒦1}}+{\displaystyle \frac{1+2W^2(v_y^2+v_z^2)}{1v_x^2}}\right\}`$   
$`q_5^{L,R}`$ $`=`$ $`{\displaystyle \frac{h^2}{\mathrm{\Delta }}}\left\{M\left[𝒜_{}\mathrm{\Omega }_+𝒟_+𝒜_+\mathrm{\Omega }_{}𝒟_{}\right]+p[c_+\mathrm{}_+𝒟_+c_{}\mathrm{}_{}𝒟_{}]\right\}+`$ (21)   
$`c_0p{\displaystyle \frac{W}{h}}\left\{{\displaystyle \frac{hW𝒦}{𝒦1}}+{\displaystyle \frac{(2hW1)(v_y^2+v_z^2)}{1v_x^2}}\right\},`$   
with the following auxiliary quantities

$`M=\rho hW^2(𝒦1),\mathrm{\Omega }_\pm =c_\pm (v_x\lambda _{}),𝒟_\pm =hW𝒜_\pm 1,`$   
$`𝒦{\displaystyle \frac{\stackrel{~}{\kappa }}{\stackrel{~}{\kappa }c_s^2}},\stackrel{~}{\kappa }={\displaystyle \frac{1}{\rho }}{\displaystyle \frac{p}{\epsilon }}|_\rho ,𝒜_\pm {\displaystyle \frac{1v_xv_x}{1v_x\lambda _\pm }}`$   
$`\mathrm{\Delta }=h^3W(𝒦1)(1v_xv_x)(𝒜_+\lambda _+𝒜_{}\lambda _{})`$   
$`\mathrm{}_\pm =\pm \left\{v_xW^2(\mathrm{v}^2v_xv_x)(2𝒦1)(v_x𝒜_\pm \lambda _\pm )+𝒦𝒜_\pm \lambda _\pm \right\}`$ (22)   
where quantities $`c_{\pm ,0}`$ are given in Table I and $`\mathrm{\Delta }`$ is the determinant of the matrix of right-eigenvectors.

The corresponding viscosity vectors in the other directions are trivially obtained by a cyclic permutation of subindices $`x,y,z`$.

We have tested the efficiency of our numerical proposal, for Roe’s and MM’s flux formulae, running GENESIS (a 3D special relativistic hydro-code ), without any optimization at compilation level, in a SGI Origin 2000. A standard initial value problem has been chosen: $`\rho _L=10`$, $`ϵ_L=2`$, $`v_L=0`$, $`\gamma _L=5/3`$, $`\rho _R=1`$, $`ϵ_R=10^6`$, $`v_R=0`$ and $`\gamma _R=5/3`$, where the subscript $`L`$ ($`R`$) denotes the state to the left (right) of the initial discontinuity. This test problem has been considered by several authors in the past (see  for details in 1D, 2D and 3D).

We have compared the performance of two different implementations of the numerical flux subroutine:

i) Case A, stands for the results obtained using our analytical prescription. This means to write down, in the numerical flux routine, just the expressions derived here for the viscosity vector $`𝐪`$.

ii) Case B, stands for the results obtained running the code with a standard high-efficiency subroutine for inverting matrices (we use a LU decomposition plus an implicit pivoting which is, for general matrices, $`O(N^3)`$). This subroutine is called to get the left eigenvectors from the matrix of right eigenvectors and is adapted to the particular dimensions of the matrices ($`3\times 3`$, in 1D, $`4\times 4`$, in 2D and $`5\times 5`$, in 3D). Hence, unlike case A, now we have to calculate numerically the following quantities: the matrix of left eigenvectors, the characteristic variables and, finally, the components of the viscosity vector $`𝐪`$.

Table II summarizes the results: the direct implementation of our numerical viscosity formulae leads to an improvement of the efficiency (in terms of CPU time) of the numerical fluxes subroutine in a factor which, in 3D calculations, ranges between about eleven and twelve depending on the particular flux formula used. When comparing Roe’s and MM’s cases a factor two –in favour of Roe– arises due to the fact that MM’s flux formulae needs to compute two viscosity vectors (one per each side of a given interface), unlike Roe’s flux formula which needs only one viscosity vector evaluated at the average state. As it must be, the efficiency increases with the number of spatial dimensions involved in the problem due to the computationally expensive matrix inverting operations performed at each interface to get the numerical fluxes. Since the numerical flux routine is, typically, one of the most time-consuming, it translates into a speed up factor between two and four in the total execution time, depending on the specific weight of the flux formulae routine in each particular application.

Our formulation also gives a unified description of the numerical fluxes (4), permitting a unique implementation with the possibility of switching in cases when the utilisation of a specific flux formula is more appropriate. In addition, due to the fact that we have eliminated, in the generalized MM’s flux formula, all the conditional clauses, the efficiency is ensured either for scalar or vectorial processors.

Another worthy by-products of our algebraic pre-processing concerns with the significant reduction of round-off errors, as a consequence of the number of operations suppressed and factorization. One of the important issues in designing a multidimensional hydro-code is the accurate preservation of any symmetries present in a physical problem. A numerical violation of these symmetries could arise as a consequence of accumulation of round-off errors in the calculation of the numerical fluxes, as we have explained in a previous paper . The algebraic simplifications, shown in the present paper, reduce the number of operations and cure such problem.

Two last additional consequences arise from our work. First is that similar expressions can be worked out for any non-linear hyperbolic system of conservation laws for which the full spectral decomposition is known. In particular, when some of the vectorial subspaces has multiple degeneracy, a similar algebraic preprocessing is very convenient. The other important consequence is that an appropriate combination of a simplified formulation of the numerical viscosity together with the use of special relativistic Riemann solvers in General Relativistic Hydrodynamics , should allow a very easy and efficient extension to General Relativistic Hydrodynamics.