# Relativistic primitive variable reconstruction in Gkeyll¶

In conservative formulations of non-relativistic hydrodynamics (as described, for instance, by the non-relativistic Euler equations), the recovery of primitive variables from conservative variables, as required for the computation of fluxes, is a purely algebraic operation: it consists simply of dividing various quantities through by the mass density \(\rho\), and then applying the equation of state to recover the pressure \(P\) from the total energy \(E\). In special and general relativistic hydrodynamics, however, this recovery operation becomes considerably more complicated, due to the appearance of the Lorentz factor:

which guarantees that the (spatial) three-momenta \(p_k\), as measured by a normal observer:

are no longer algebraically independent of one another, in contrast to the non-relativistic case. Thus, for a generic equation of state, it is argued by [Marti2003] that the conservative to primitive variable reconstruction in relativistic hydrodynamics cannot be represented as a closed-form algebraic operation, and instead requires one to perform some kind of (potentially higher-dimensional) root-finding operation. To this end, Gkeyll employs a certain “robustified” variant of the algorithm proposed by [Eulderink1995], which has been specifically modified to accommodate low densities, low pressures, and high values of the Lorentz factor, without going unstable, whilst still retaining the favorable convergence properties of the original Eulderink and Mellema algorithm in less extreme cases. In this short technical note, we will follow the notational and terminological conventions introduced within the description of the GRHD equations in Gkeyll. In particular, we assume a \({3 + 1}\) “ADM” decomposition of our 4-dimensional spacetime into 3-dimensional spacelike hypersurfaces with induced metric tensor \(\gamma_{i j}\), lapse function \(\alpha\), and shift vector \(\beta^i\). We also assume a perfect relativistic fluid with (rest) mass density \(\rho\), spatial velocity components (as perceived by normal observers) \(v^i\), pressure \(P\), specific enthalpy \(h\), and local sound speed \(c_s\). Einstein summation convention (in which repeated tensor indices are implicitly summed over) is assumed throughout. The Greek indices \(\mu, \nu\) range over the full spacetime coordinate directions \(\left\lbrace 0, \dots, 3 \right\rbrace\), while the Latin indices \(i, j, k, l\) range over the spatial coordinate directions \(\left\lbrace 1, \dots, 3 \right\rbrace\) only.

See also this note on the complete eigensystem for the GRHD equations (as implemented as part of Gkeyll’s stable time-step calculation), and this note on Gkeyll’s handling of general black hole spacetimes and the Kerr-Schild coordinate system.

## Reconstruction algorithm for an ideal gas equation of state¶

Here we assume an ideal gas equation of state, with adiabatic index \(\Gamma\), for which the specific relativistic enthalpy:

with specific internal energy \(\varepsilon \left( \rho, P \right)\), is now given by:

such that the local fluid sound speed:

now becomes:

The approach advocated by [Eulderink1995] is then to use a non-linear root-finding algorithm (namely the one-dimensional Newton-Raphson method) to find the roots of the following quartic polynomial in \(\xi\):

where we have defined the variable \(\xi\) and the constant term \(\eta\) to be:

and

respectively, and where the coefficients \(\alpha_4\), \(\alpha_2\), \(\alpha_1\) and \(\alpha_0\) in the quartic are given by:

and:

respectively. Once an appropriately accurate value of \(\xi\) has been calculated in this way, it is then possible to compute an updated value for the Lorentz factor \(W_{new}\) using:

and, from this, an updated value for the fluid (rest) mass density \(\rho_{new}\) using:

It is also possible to compute, from the numerical value of \(\xi\), an updated value for the specific relativistic enthalpy \(h_{new}\) using:

from which, in turn, an updated value for the hydrostatic pressure \(P_{new}\) can also be recovered by means of the ideal gas equation of state:

Finally, using the updated values of all three quantities (i.e. \(\rho_{new}\), \(h_{new}\), and \(W_{new}\)) in conjunction, we are able to reconstruct the updated components of the spatial velocity (co)vector \(v_{k}^{new}\), as measured by observers moving normal to the spacelike hypersurfaces in the foliation, as:

One of the key modifications that we have made to the original algorithm of Eulderink and Mellema within Gkeyll is the extraction of the prefactor:

outside of the brackets in the expression for the updated value of the Lorentz factor \(W_{new}\), thus ensuring that the term inside the brackets remains of order unity, even for very large values of the Lorentz factor, which guarantees much greater numerical stability in the strongly relativistic fluid regime. In contrast to the prescription of [Marti2003] for ensuring numerical stability in the presence of low densities, low pressures, and high fluid velocities, instead of imposing a minimum pressure \(P_{min}\) of the general form:

for some small \(\varepsilon > 0\), we opt rather to restrict the fluid velocity so as to ensure that the physicality condition:

remains consistently satisfied throughout the computational domain. This avoids the kind of finite-precision arithmetic problems described by [Goldberg1991] that one would otherwise experience when computing the updated Lorentz factor \(W_{new}\) using the minimum-pressure prescription. In Gkeyll we have opted to use a rapidly-converging Newton-Raphson iteration scheme to construct an approximate solution for \(\xi\), rather than solving the quartic analytically, since we have found that the slight increase in numerical accuracy obtained from an analytical solution is not worth the considerable increase in computational complexity associated with calculating it.

## Simplified reconstruction for a stiff equation of state¶

We start from the general form of the ultra-relativistic equation of state, with adiabatic index \(\Gamma\):

and we set \(\Gamma = 2\), yielding \(P = \rho\), known as a *stiff* relativistic
equation of state. The equivalence of fluid density and pressure simplifies the equations
of fluid evolution considerably, and renders the conservative to primitive reconstruction
step a purely algebraic operation, just as it is for non-relativistic fluids. Following
the approach of [Neilsen1999], we are able to recover the updated fluid density
\(\rho_{new}\) as:

and the updated components of the spatial velocity (co)vector \(v_{k}^{new}\), as measured by normal observers, as:

As a consequence, the Gkeyll fluid solver for stiff relativistic fluids is considerably more computationally efficient than it is for general ideal or ultra-relativistic fluids, due to the lack of need for a Newton-Raphson iteration scheme as part of the reconstruction step.

## References¶

J. M. Martí and E. Müller, “Numerical Hydrodynamics in Special
Relativity”, *Living Reviews in Relativity* **6** (7). 2003.

F. Eulderink and G. Mellema, “General Relativistic Hydrodynamics
with a Roe solver”, *Astronomy and Astrophysics Supplement Series* **110**: 587-623.
1995.

D. Golderberg, “What every computer scientist should know about
floating-point arithmetic”, *ACM Computing Surveys* **23** (1): 5-48. 1991.

D. W. Neilsen, *Extremely Relativistic Fluids in Strong-Field Gravity*,
PhD Thesis, University of Texas at Austin. 1999.