The GRHydro code: three-dimensional relativistic hydrodynamics

Original authors: Luca Baiotti, Ian Hawke, Pedro Montero
Contributors: Sebastiano Bernuzzi, Giovanni Corvino, Toni Font, Joachim Frieben,
Roberto De Pietri, Thorsten Kellermann, Frank Löffler, Christian D. Ott,
Luciano Rezzolla, Nikolaos Stergioulas, Aaryn Tonita,
and many others,
especially those who were in the European Network “Sources of Gravitational Waves”

Date: 2009-12-07 19:20:47 -0600 (Mon, 07 Dec 2009)


GRHydro is a fully general-relativistic three-dimensional hydrodynamics code. It evolves the hydrodynamics using High Resolution Shock Capturing methods and can work with realistic equations of state. The evolution of the spacetime can be done by any other “appropriate” thorn, such as the CCATIE code, maintained and developed at the Albert-Einstein-Institut (Potsdam).

1 Introduction

The GRHydro1 code is based upon the public version of the Whisky code which is a product of the EU Network on Sources of Gravitational Radiation2 , and was initially written by Luca Baiotti, Ian Hawke and Pedro Montero, based on the publicly available GR3D code and with many other important contributors. With the help of large parts of the community, the GRHydro code got improved, extended and included into the Einstein Toolkit.

2 Using This Thorn

What follows is a brief introduction to using GRHydro. It assumes that you know the required physics and numerical methods, and also the basics of Cactus3 . If you don’t, then skip this section and come back to it after reading the rest of this ThornGuide of Cactus. For more details such as thornlists and parameter files, take a look at the Einstein Toolkit web page which is currently stored at

GRHydro uses the hydro variables defined in HydroBase and provides own “conserved” hydro variables and methods to evolve them. It does not provide any information about initial data or equations of state. For these, other thorns are required. A minimal list of thorns for performing a shock-tube test is given in the shock-tube test parameter file, found at


and will include the essential thorns

GRHydro EOS_Omni ADMBase ADMCoupling MoL

Initial data for shocks can be set using


Initial data for spherically symmetric static stars (with perturbations or multiple “glued” stars) can be set by


The actual evolution in time is controlled by the Method of Lines thorn MoL. For full details see the relevant ThornGuide. For the purposes of GRHydro at least two parameters are relevant; ode_method and mol_timesteps. If second-order accuracy is all that is required then ode_method can be set to either "rk2" (second-order TVD Runge-Kutta evolution) or "icn" (Iterative Crank Nicholson, number of iterations controlled by mol_timesteps, defaults to three), and mol_timesteps can be ignored. A more generic (and hence less efficient) method can be chosen by setting ode_method to "genrk" which is a Shu-Osher type TVD Runge-Kutta evolution. Then the parameter mol_timesteps controls the number of intermediate steps and hence the order of accuracy. First to seventh order are currently supported.

GRHydro currently uses a Reconstruction-Evolution type method. The type of reconstruction is controlled by the parameter recon_method. The currently supported values are "tvd" for slope limited TVD reconstruction, "ppm" for the Colella-Woodward PPM method, and "eno" for the Essentially Non-Oscillatory method of Harten et al. Each of these has further controlling parameters. For example there are a number of slope limiters controlled by the keyword tvd_limiter, the PPM method supports shock detection by the Boolean ppm_detect, and the ENO method can have various orders of accuracy controlled by eno_order. Note that the higher-order methods such as PPM and ENO require the stencil size to be increased using GRHydro_stencil and driver::ghost_size.

For the evolution various approximate Riemann solvers are available, controlled by riemann_solver. Currently implemented are "HLLE", "Roe" and "Marquina". For the Roe and Marquina methods there are added options to choose which method is used for calculating the left eigenvectors. This now defaults to the efficient methods of the Valencia group, but the explicit matrix inversion is still there for reference.

For the equations of state, two “types” are recognized, controlled by the parameter GRHydro_eos_type. These are "Polytype" where the pressure is a function of the rest-mass density, P = P(ρ), and the more generic "General" type where the pressure is a function of the rest-mass density and the specific internal energy, P = P(ρ,𝜖). For the Polytype equations of state one fewer equation is evolved and the specific internal energy is set directly from the rest-mass density. The actual equation of state used is controlled by the parameter GRHydro_eos_table. Polytype equations of state include "2D_Polytrope" and general equations of state include "Ideal_Fluid".

2.1 Obtaining This Thorn

The public version of GRHydro can be found on the website

2.2 Basic Usage

The simplest way to start using GRHydro would be to download some example parameter files from the web page to try it. There are a number of essential parameters which might reward some experimentation. These include:

2.3 Special Behaviour

Although in theory GRHydro can deal with conformal metrics as well as physical metrics, this part of the code is completely untested as we don’t have conformal initial data (although this would not be hard - we just haven’t had the incentive).

2.4 Interaction With Other Thorns

GRHydro provides the appropriate contribution to the stress energy through the TmunuBase interface. Those spacetime evolvers that use this interface can use GRHydro without change.

GRHydro uses the MoL thorn to perform the actual time evolution. This means that if all other evolution thorns are also using MoL then the complete evolution will have the accuracy of the MoL evolution method without change. This (currently) allows for up to fourth-order accuracy in time without any changes to GRHydro.

For the general equations of state GRHydro uses the EOS_Omni interface. This returns the necessary hydrodynamical quantities, such as the pressure and derivatives with general function calls. The parameter GRHydro_eos_table controls which equation of state is used during evolution.

For the metric quantities GRHydro uses the standard CactusEinstein arrangement, especially ADMBase. This allows the standard thorns to be used for the calculation of constraint violations, emission of gravitational waves, location of the apparent horizon, and more.

2.5 Support and Feedback

GRHydro is part of the Einstein Toolkit, with its web page located at

It contains information on obtaining the code, together with thornlists and sample parameter files. For questions, send an email to the Einstein Toolkit mailing list

3 Physical System

The equations of general relativistic hydrodynamics can be written in the flux conservative form

tq + xif(i)(q) = s(q), (1)

where q is a set of conserved variables, f(i)(q) the fluxes and s(q) the source terms. The five conserved variables are labeled D, Si, and τ, where D is the generalized particle number density, Si are the generalized momenta in each direction, and τ is an internal energy term. These conserved variables are composed from a set of primitive variables, which are ρ, the rest-mass density, p, the pressure, vi, the fluid 3-velocities, 𝜖, the specific internal energy, and W, the Lorentz factor, via the following relations  

D = γWρ Si = γρhW2 vi τ = γ ρhW2 p D, (2)

where γ is the determinant of the spatial 3-metric γij and h 1 + 𝜖 + pρ. Only five of the primitive variables are independent. Usually the Lorentz factor is defined in terms of the velocities and the metric as W = (1 γijvivj)12. Also one of the pressure, rest-mass density or specific internal energy terms is given in terms of the other two by an equation of state.

The fluxes are usually defined in terms of both the conserved variables and the primitive variables:

Fi(U) = [D(αvi βi),S j(αvi βi) + pδ ji,τ(αvi βi) + pvi]T. (3)

The source terms are

s(U) = 0,Tμν μgνj + Γμνδg δj,αTμ0 μ lnα TμνΓ νμ0. (4)

Note that the source terms do not contain differential operators acting on the stress-energy tensor and that this is important for the numerical preservation of the hyperbolicity character of the system. Also note that in a curved spacetime the equations are not in a strictly-homogeneous conservative form, which is recovered only in flat spacetime. This conservative form of the relativistic Euler equations was first derived by the group at the Universidad de Valencia [3] and therefore was named the Valencia formulation.

For more detail see the review of Font [9].

4 Numerical Implementation

5 High Resolution Shock Capturing methods

The numerical evolution of a hydrodynamical problem is complicated by the occurrence of shocks, i.e. genuine nonlinear discontinuities that will generically form. It is also complicated by the conservation constraint. In a High Resolution Shock Capturing (HRSC) method the requirement of conservation is used to ensure the correct evolution of a shock. A HRSC method also avoids spurious oscillations at shocks which are known as Gibbs’ phenomena, while retaining a high order of accuracy over the majority of the domain.

For a full introduction to HRSC methods see [14], [20], [15], [16] and [9].

In the GRHydro code it was decided to use the method of lines as a base for the HRSC scheme. The method of lines is a way of turning a partial differential equation such as (1) into an ordinary differential equation. For the GRHydro code the following steps are required.

This ordinary differential equation can be solved by the Cactus thorn MoL. All that GRHydro has to do is to return the values of the discrete spatial differential operator

L(q) = sd3x +xj122xj+122 xk123xk+123 f(1)(q(x i121,y,z))dydz xj122xj+122 xk123xk+123 f(1)(q(x i+121,y,z))dydz + (6)

given the data q, and to supply the boundary conditions that will be required to calculate this right hand side at the next time level. We note that in the current implementation of MoL the solution to the ordinary differential equation (5) will be Nth-order accurate provided that the time integrator used by MoL is Nth-order accurate in time, and that the discrete operator L is Nth-order accurate in space and first-order (or better) accurate in time. For more details on the method of lines, and the options given with the time integration for MoL, see the relevant chapter in the ThornGuide.

In this implementation of GRHydro the right hand side operator L will be simplified considerably by approximating the integrals by the midpoint rule to get

L(q) = si,j,k + fi12,j,k(1) f i+12,j,k(1) + (7)

where the dependence of the flux on q and spatial position is implicit in the notation. Given this simplification, the calculation of the right hand side operator splits simply into the following two parts:

  1. Calculate the source terms s(q(xi1,xj2,xk3)):

    Given that q̄ is a second-order accurate approximation to q at the midpoint of the cell, which is precisely (xi1,xj2,xk3), for second-order accuracy it is sufficient to use s(q̄i,j,k).

  2. In each coordinate direction xa, calculate the intercell flux f(a)(qi+12,j,k):

    From the initial data q̄ given at time tn we can reconstruct the data at the cell boundary, (qi+12,j,k) to any required accuracy in space (except in the vicinity of a shock, where only first-order accuracy is guaranteed). However this will only be zeroth-order accurate in time. To ensure first-order accuracy in time, we have to find (qi+12,j,k)(t) while retaining the high spatial order of accuracy. This requires two steps:

    1. Reconstruct the data q over the cells adjacent to the required cell boundary. This reconstruction should use the high spatial order of accuracy. This gives two values of (qi+12,j,k), which we call qL and qR, where qL is obtained from cell i (left cell) and qR from cell i + 1 (right cell).
    2. The values qL,R are used as initial data for the Riemann problem. This is the initial value problem given by the partial differential equation (1) with semi-infinite piecewise constant initial data qL,R. As the true function q is probably not piecewise constant we will not get the exact solution of the general problem even if we solve the local Riemann problem exactly. However, it will be first-order accurate in time and retain the high order of accuracy in space which, as explained in the documentation to the MoL thorn, is sufficient to ensure that the method as a whole has a high order of accuracy.

So, the difficult part of GRHydro is expressed in two routines. One that reconstructs the function q at the boundaries of a computational cell given the cell average data q̄, and another that calculates the intercell flux f at this cell boundary.

6 Reconstruction

In the reduction of all of GRHydro to two routines in the last section one point was glossed over. That is, in order for the numerical method to be consistent and convergent it must retain conservation and not introduce spurious oscillations. Up to this point all the steps have either been exact or have neither changed the conservation properties or the profile of the function. Also, the calculation of the intercell flux from the Riemann problem can be made to be “exactly correct”. That is, even though as explained above it may not be the true flux for the real function q, it will be the exact physical solution for the values qL,R given by the reconstruction routine, so the intercell flux cannot violate conservation or introduce oscillations. Unphysical effects such as these can only be introduced by an incorrect reconstruction.

For a full explanation of reconstruction methods see Laney [14], Toro [20], Leveque [15]. For the moment we will concentrate on the simplest methods that are better than first-order accurate in space, the TVD slope-limited methods. More complex methods such as ENO will be introduced later.

In the late 1950’s Godunov proved a theorem that, in this context, says

Any linear reconstruction method of higher-than-first-order accuracy may introduce spurious oscillations.

For this theorem linear meant that the reconstruction method was independent of the data it was reconstructing. Simple centred differencing is a linear second-order method and is oscillatory near a shock. Instead we must find a reconstruction method that depends on the data q being reconstructed.

Switching our attention to conservation, we note that there is precisely one conservative first-order reconstruction method,

qFirst(x) = q̄ i,x [xi12,xi+12], (8)

and that any second-order conservative method can be written in terms of a slope or rather difference Δi as

qSecond(x) = q̄ i + x xi xi+12 xi12Δi,x [xi12,xi+12]. (9)

6.1 TVD Reconstruction

As we want a method that is accurate (i.e., at least to second order) while being stable (i.e., only first order or nonlinear at shocks) then the obvious thing to do is to use some second-order method in the form of equation (9) in smooth regions but which changes to the form of equation (8) smoothly near shocks.

In the articles describing the GRAstro_Hydro code4 , this was described as an average of the two reconstructions,

qTV D(x) = ϕ(q)qSecond + (1 ϕ(q))qFirst, (10)

where ϕ [0,1] varies smoothly in some sense, and is zero near a shock and 1 in smooth regions. In Toro’s notation [20] (which we usually adopt here) the slope limiter function ϕ (having the same attributes as above) directly multiplies the slope, giving

qTV D(x) = q̄ i + x xi xi+12 xi12ϕ(q)Δi,x [xi12,xi+12]. (11)

Equations (10) and (11) are equivalent.

For details on how to construct a limiter, on their stability regions and on the explicit expressions for the limiters used here, see [20]. The GRHydro code implements the minmod limiter (the most diffusive and the default), the Van Leer Monotonized Centred (MC) (VanLeerMC) limiter in a number of forms (which should give equivalent results), and the Superbee limiter. The limiter specified by the parameter value VanLeerMC2 is the recommended one.

As an example, we show how TVD with the minmod limiter is implemented in the code. First, we define the minmod function:

minmod(a,b) = min(|a|,|b|)if (ab > 0) 0 otherwise. (12)

For reconstructing q we choose two differences

Δupw =qi qi1 Δloc =qi+1 qi, (13)

and write

qTV D,minmod(x) = q̄ i + x xi xi+12 xi12minmod(Δupw,Δloc),x [xi12,xi+12]. (14)

6.2 PPM reconstruction

The piecewise parabolic method (PPM) of Colella and Woodward [4] is a rather more complex method that requires a number of steps. The implementation in the GRHydro code is specialized to use evenly spaced grids. Also, some of the more complex features are not implemented; in particular, the dissipation algorithm is only the simplest given in the original article. Here we just give the implementation details. For more details on the method we refer to the original article.

Again we assume we are reconstructing a scalar function q as a function of x in one dimension on an evenly spaced grid, with spacing Δx. The first step is to interpolate a quadratic polynomial to the cell boundary,

qi+12 = 1 2 qi+1 + qi + 1 6 δmqi δmqi+1 , (15)


δmqi = min(|δqi|,2|qi+1 qi|,2|qi qi1|) sign(δqi)if (qi+1 qi)(qi qi1) > 0, 0 otherwise. , (16)


δqi = 1 2(qi+1 qi1). (17)

At this point we set both left and right states at the interface to be equal to this,

qiR = q i+1L = q 1+12. (18)

This reconstruction will be oscillatory near shocks. To preserve monotonicity, the following replacements are made:

qiL = q iR = q i if (qiR q i)(qi qiL) 0 (19) qiL = 3q i 2qiR if (q iR q iL) q i 1 2(qiL + q iR) > 1 6(qiR q iL)2 (20) qiR = 3q i 2qiL if (q iR q iL) q i 1 2(qiL + q iR) < 1 6(qiR q iL)2. (21)

However, before applying the monotonicity preservation two other steps may be applied. Firstly we may steepen discontinuities. This is to ensure sharp profiles and is only applied to contact discontinuities. This may be switched on or off using the parameter ppm_detect. This part of the method replaces the cell boundary reconstructions of the density with

ρiL = ρ iL(1 η) + ρ i1 + 1 2δmρi1 η (22) ρiR = ρ iR(1 η) + ρ i+1 1 2δmρi+1 η (23)

where η is only applied if the discontinuity is mostly a contact (see [4] for the details) and is defined as

η = max(0,min(1,η1(η̃ η2))), (24)

where η1,η2 are positive constants and

η̃ = ρi2ρi+2+4δρi 12δρi  if  δ2ρi+1δ2ρi1 < 0 (ρi+1 ρi1) 𝜖min(|ρi+1|,|ρi1|) > 0 0 otherwise ,

with 𝜖 another positive constant and

δ2ρ i = ρi+1 2ρi + ρi1 6Δx2 . (25)

Another step that is performed before monotonicity enforcement is to flatten the zone structure near shocks. This adds simple dissipation and is always in the code. In short, the reconstructions are again altered to

qiL,R = ν iqiL,R + (1 ν i)qi, (26)


νi = max(0,1 max(0,ω2(pi+1pi1 pi+2pi2 ω1))) if ω0min(pi1,pi+1) < |pi+1 pi1| and vi1x vi+1x > 0 1 otherwise (27)

and ω0,ω1,ω2 are constants.

The above flattening procedure is not the one in the original article of Colella and Woodward, but it has been adapted from it in order to have a stencil of three points. The original flattening procedure is also implemented in GRHydro. Instead of 26, it consists in the formula

qiL,R = ν̃ iqiL,R + (1 ν̃ i)qi, (28)


ν̃i = maxνi,νi+sign(pi1pi+1) (29)

and νi is given by 27. This can be activated by setting the parameter ppm_flatten to stencil_4. Formula 28, despite requiring more computational resources (especially when mesh refinement is used), usually gives very similar results to 26, so we routinely use 26.

6.3 ENO Reconstruction

An alternative way of getting higher-than-second-order accuracy is the implementation of the Essentially Non-Oscillatory methods of Harten [11]. The essential idea is to alter the stencil to use those points giving the smoothest reconstruction. The only restriction is that the stencil must include the cell to be reconstructed (for stability). Here we describe the simplest ENO type reconstruction: piecewise polynomial reconstruction using the (un)divided differences to measure the smoothness.

Let k be the order of the reconstruction. Suppose we are reconstructing the scalar function q in cell i. We start with the cell Ii. We then add to the stencil cell Ij, where j = i ± 1, where we choose j to minimize the Newton divided differences

q xi1,xi = qi qi1 xi+12 xi32 (30) q xi,xi+1 = qi+1 qi xi+32 xi12. (31)

We then recursively add more cells, minimizing the higher-order Newton divided differences q xik,,xi+j defined by

q xik,,xi+j = q xik+1,,xi+j q xik,,xi+j1 xi+j xik . (32)

The reconstruction at the cell boundary is given by a standard kth-order polynomial interpolation on the chosen stencil.

[19] has outlined an elegant way of calculating the cell boundary values solely in terms of the stencil and the known data. If the stencil is given by

S(i) = Iir,,Ii+kr1 , (33)

for some integer r, then there exist constants crj depending only on the grid xi such that the boundary values for cell Ii are given by

qi+12 = j=0k1c rjqir+j,qi12 = j=0k1c r1,jqir+j. (34)

The constants crj are given by the rather complicated formula

crj = m=j+1k l=0,lmk q=0,qm,lk xi+12 xir+q12 l=0,lmk xir+m12 xir+L12 Δxir+j. (35)

This simplifies considerably if the grid is even. The coefficients for an even grid are given (up to seventh order) by [19].

7 Riemann Problems

Given the reconstructed data, we then solve a local Riemann problem in order to get the intercell flux. The Riemann problem is specified by an equation in flux conservative homogeneous form,

tq + xif(i)(q) = 0 (36)

with piecewise constant initial data qL,R separated by a discontinuity at x(1) = 0. Flux terms for the other directions are given similarly. There is no intrinsic scale to this problem and so the solution must be a self similar solution with similarity variable ξ = x(1)t. The solution is given in terms of waves which separate different states, with each state being constant. The waves are either shocks, across which all hydrodynamical variables change discontinuously, rarefactions, across which all the variables change continuously (the wave is not a single value of ξ for a rarefaction, but spreads across a finite range), or contact or tangential discontinuities, across which some but not all of the hydrodynamical variables change discontinuously and the rest are constant. The characteristics of the matter evolution converge and break at a shock, diverge at a rarefaction and are parallel at the other linear discontinuities.

The best references for solving the Riemann problem either exactly or approximately are [15], [20], [14]. Here, we start by giving a simple outline. We start by considering the N dimensional linear problem in one dimension given by

tq + Axq = 0, (37)

where A is a N × N matrix with constant coefficients. We define the eigenvalues λj with associated right and left eigenvectors rj,lj, where the eigenvectors are normalized to each other (i.e., their dot product is δji). We shall assume that the eigenvectors span the space. The characteristic variables wi are defined by

wi = li q. (38)

Then equation (37) when written in terms of the characteristic variables becomes

tw + Λxw = 0, (39)

where Λ is the matrix containing the eigenvalues λi on the diagonals and zeros elsewhere. Hence each characteristic variable wi obeys the linear advection equation with velocity a = λi. This solves the Riemann problem in terms of characteristic variables.

In order to write the solution in terms of the original variables q we order the variables in such a way that λ1 λN. We also define the differences in the characteristic variables Δwi = (wi)L (wi)R across the ith characteristic wave. These differences are single numbers (‘scalars’). We note that these differences can easily be found from the initial data using

Δwi = li qL qR . (40)

As the change in the solution across each wave is precisely the difference in the associated characteristic variable, the solution of the Riemann problem in terms of characteristic variables can be written as either

wi = (wi)L + j=1MΔw jejifλ M < ξ < λM+1, (41)


w = (wi)R j=M+1NΔw jejifλ M < ξ < λM+1, (42)

or as the average

wi = 1 2 (wi)L + (wi)R + j=1MΔw jej j=M+1NΔw jej ifλ M < ξ < λM+1, (43)

where ei is the column vector (ei)j = δji.

Converting back to the original variables q we have the solution

q = 1 2 qL + qR + i=1MΔw iri i=M+1NΔw iri ifλ M < ξ < λM+1. (44)

In the case where we are only interested in the flux along the characteristic ξ = 0 we can write the solution in the simple form

f(q) = 1 2 f(qL) + f(qR) i=1N|λ i|Δwiri . (45)

All exact Riemann solvers have to solve at least an implicit equation and so are computationally very expensive. As the solution of Riemann problems takes a large portion of the time to run in a HRSC code, approximations that speed the calculation of the intercell flux are often essential. This is especially true in higher dimensions (¿1), where the solution of the ordinary differential equation to give the relation across a rarefaction wave makes the use of an exact Riemann solver impractical.

Approximate Riemann solvers can have problems, as shown in depth by Quirk [17]. Hence it is best to compare the results of as many different solvers as possible. Here we shall describe the three approximate solvers used in this code, starting with the simplest.

7.1 HLLE solver

The Harten-Lax-van Leer-Einfeldt (HLLE) solver of Einfeldt [7] is a simple two-wave approximation. We assume that the maximum and minimum wave speeds ξ± are known. The solution of the Riemann problem is then given by requiring conservation to hold across the waves. The solution takes the form

q = qLifξ < ξ qifξ < ξ < ξ+ qRifξ > ξ+, (46)

and the intermediate state q is given by

q = ξ+qR ξqL f(qR) + f(qL) ξ+ ξ . (47)

If we just want the numerical flux along the boundary then this takes the form

f(q) = ξ̂+f(qL) ξ̂f(qR) + ξ̂+ξ̂(qR qL) ξ̂+ ξ̂ , (48)


ξ̂ = min(0,ξ),ξ̂+ = max(0,ξ+). (49)

Knowledge of the precise minimum and maximum characteristic velocities ξ± requires knowing the solution of the Riemann problem. Instead, the characteristic velocities are usually found from the eigenvalues of the Jacobian matrix fq evaluated at some intermediate state. To ensure that the maximum and minimum eigenvalues over the entire range between the left and right states are found, we evaluate the Jacobian in both the left and right states and take the maximum and minimum over all eigenvalues. This ensures, for the systems of equations considered here, that the real maximum and minimum characteristic velocities are contained within [ξ,ξ+].

If we set α = max(|ξ|,|ξ+|) and replace the characteristic velocities ξ± with ± α, we find the Lax–Friedrichs flux (cf. also Tadmor’s semi-discrete scheme [13])

f(q) = 1 2 f(qL) + f(qR) + α(qL qR). (50)

This is very diffusive, but also very stable.

7.2 Roe solver

The linearized solver of Roe [18] is probably the most popular approximate Riemann solver. The simplest interpretation is that the Jacobian fq is linearized about some intermediate state. Then the conservation form reduces to the linear equation

tq + Axq = 0, (51)

where A is a constant coefficient matrix. This is identical to equation (37) and so all the results of section 7 on linear systems hold. We reiterate that the standard form for the flux along the characteristic ray ξ = 0 is

f(q) = 1 2 f(qL) + f(qR) i=1N|λ i|Δwiri . (52)

There is a choice of which intermediate state the Jacobian should be evaluated at. Roe gives three criteria that ensure the consistency and stability of the numerical flux:

  1. A(qRoe) qR qL = f(qR) f(qL),
  2. A(qRoe) is diagonalizable with real eigenvalues,
  3. A(qRoe) fq smoothly as qRoe q.

A true Roe average for relativistic hydrodynamics, i.e., an intermediate state that satisfies all these conditions, has been constructed by Eulderink [8]. However, frequently it is sufficient to use

qRoe = 1 2 qR + qL , (53)

which satisfies only the last two conditions. For simplicity we have implemented this arithmetic average.

7.3 Marquina solver

Unlike all the other Riemann solvers introduced so far, the Marquina solver as outlined in [6] does not solve the Riemann problem completely. Instead, only the flux along the characteristic ray ξ = 0 is given. It can be seen as a generalized Roe solver, as the results are the same except at sonic points. These points are where the fluid velocity is equal to the speed of sound of the fluid. In the context of Riemann problems, sonic points are found when the ray ξ = 0 is within a rarefaction wave.

Firstly define the left l(qL,R) and right r(qL,R) eigenvectors and the eigenvalues λ(qL,R) of the Jacobian matrix fq evaluated at the left and right states. Next define left and right characteristic variables wL,R and fluxes ϕL,R by

(wi)L,R = li(qL,R) qL,R,(ϕi)L,R = li(qL,R) f(qL,R). (54)

Then the algorithm chooses the correct-sided characteristic flux if the eigenvalues λi(qL), λi(qR) have the same sign, and uses a Lax–Friedrichs type flux if they change sign. In full, the algorithm is given in figure 1.

Fori = 1,,Ndo Ifλi(qL)λi(qR) > 0then Ifλi(qL) > 0then ϕ+i =ϕLi ϕi =0 else ϕ+i =0 ϕi =ϕRi endif else αi =max(|λi(qL),λi(qR)|) ϕ+i =1 2 ϕLi + αiw Li ϕi =1 2 ϕRi αiw Ri endif enddo (55)

Figure 1: The algorithm to calculate the Marquina flux.

Then the numerical flux is given by

f(q) = i=1N ϕ +iri(q L) + ϕiri(q R). (56)

The above implementation is based on [1].

8 Other points in GRHydro

There are a number of other things done by GRHydro which, whilst not as important as reconstruction and evolution, are still essential.

8.1 Source terms

In a curved spacetime the equations are not in homogeneous conservation-law form but also contain source terms. These are actually calculated first, before the flux terms (it simplifies the loop very slightly). There are a few points to note about the calculation of the sources.

In what follows Greek letters range from 0 to 3 and roman letters from 1 to 3.

For the following computations, we need the expression of some of the 4-Christoffel symbols (4)Γμνρ applied to the 3+1 decomposed variables. In order to remove time derivatives we will frequently make use of the ADM evolution equation for the 3-metric in the form

tγij = 2 αKij + (iβj) (3)Γ ijkβ k . (57)

As it is used in what follows, we also recall that is the covariant derivative associated with the spatial 3-surface and we note that it is compatible with the 3-metric:

iγjk = iγjk + 2(3)Γ iljγlk = 0. (58)

We start from the (4)Γ000 symbol:

(4)Γ 000 = 1 2α2 tβkβk + 2α tα + 2βi tβi βi iβkβk + 2αβi iα (59)

and we expand the derivatives as

tβkβk = tγjkβjβk = 2γ jkβj tβk + βjβk tγjk = = 2βktβk 2αK jkβjβk + 2βjβk jβk 2(3)Γ kjiβ iβjβk (60)


iβkβk = iγjkβ jβk = 2γjkβ jiβk + βjβkiγjk = 2β kiβk 2(3)Γ ikjβ jβk, (61)

where we have used (57) and (58), respectively. Inserting (60) and (61), equation (59) becomes

(4)Γ 000 = 1 αtα + βi iα + Kjkβjβk. (62)

With the same strategy we then compute

(4)Γ i00 = 1 2α2i(βkβ k α2) βj( iβj jβi + tγij) = 1 αiα βjK ij (63)


(4)Γ ij0 = 1 2α2iβj + jβi tγij βk( iγkj + jγki kγij) = 1 αKij. (64)

Other more straightforward calculations give

(4)Γ 00j = (4)Γ 0jνg ν0 = 1 2j βkβk α2 , (65) (4)Γ l0j = (4)Γ ljνg ν0 = αKlj + lβj + jβl βk(3)Γ ljk, (66) (4)Γ 0lj = (4)Γ 0jνg νl = αKjl + lβj βk(3)Γ ljk, (67) (4)Γ lmj = (4)Γ ljνg νm = (3)Γ lmj, (68)

where (57) was used to derive (66) and (67).

8.1.1 Source term for Sk

Now we have all the expressions for calculating the source terms. The ones for the variables Sk are

𝒮Skj = TνμΓ μjν = TμνΓ μνj. (69)

After expanding the derivative in (65), the coefficient of the T00 term in (69) becomes

(4)Γ 00j = 1 2βlβm jγlm αjα + βmjβm. (70)

The coefficient of the T0i term is

(4)Γ 0ij + (4)Γ i0j = jβi = βl iγjl + γiljβl. (71)

The coefficient of the Tlm term is simply

(3)Γ lmj = 1 2jγml + mγjl lγmj. (72)

Finally, summing (70)–(72) we find

𝒮Skj = T00 1 2βlβm jγlm αjα + T0iβl jγil + Ti0 jβi + 1 2Tlm jγlm, (73)

which is the expression implemented in the code.

8.1.2 Source term for τ

The source term for τ is [cf. (4)]

𝒮τ = α Tμ0 μα αTμν(4)Γ μν0 . (74)

For clarity, again we consider separately the terms containing as a factor the different components of Tμν. From (62) we find the coefficient of T00 to be

αtα α(4)Γ 000 = αβi iα + βkβlK kl. (75)

The coefficient of the term T0i is given by

αiα 2α(4)Γ i00 = 2αβjK ij αiα (76)

and, finally, the coefficient for Tij is

α2(4)Γ ij0 = αK ij. (77)

The final expression implemented in the code is thus

𝒮τ = αT00 βiβjK ij βi iα + T0i iα + 2βjK ij + TijK ij. (78)

8.2 Conversion from conservative to primitive variables

As noted in section 3 the variables that are evolved are the conserved variables D,Si,τ. But in order to calculate the fluxes and sources we require the primitive variables ρ,vi,P. Conversion from primitive to conservative is given analytically by equation (2). Converting in the other direction is not possible in a closed form except in certain special circumstances.

There are a number of methods for converting from conservative to primitive variables; see [16]. Here we use a Newton-Raphson type iteration. If we are using a general equation of state such as an ideal gas, then we find a root of the pressure equation. Given an initial guess for the pressure P̄ we find the root of the function

f = P̄ P(ρ̄,𝜖̄), (79)

where the approximate density and specific internal energy are given by

ρ̄ = D̃ τ̃ + P̄ + D̃ (τ̃ + P̄ + D̃)2 S2, (80) W̄ = τ̃ + P̄ + D̃ (τ ̃ + P ̄ + D ̃ )2 S2, (81) 𝜖̄ = D̃1 (τ ̃ + P ̄ + D ̃ )2 S2 P̄W̄ D̃. (82)

Here the conserved variables are all “undensitized”, e.g.,

D̃ = γ12D, (83)

where γ is the determinant of the 3-metric, and S2 is given by

S2 = γijS̃ iS̃j. (84)

In order to perform a Newton-Raphson iteration we need the derivative of the function with respect to the dependent variable, here the pressure. This is given by

f = 1 P ρ ρ P P 𝜖 𝜖 P, (85)

where P ρ and P 𝜖 given by calls to EOS_Base, and

ρ P = D̃S2 (τ ̃ + P ̄ + D ̃ )2 S2(τ̃ + P̄ + D̃)2, (86) 𝜖 P = P̄S2 ρ (τ̃ + P̄ + D̃)2 S2 (τ̃ + P̄ + D̃). (87) (88)

For a polytropic type equation of state, the function is given by

f = ρ̄W̄ D̃, (89)

where ρ̄ is the variable solved for, the pressure, specific internal energy and enthalpy h̄ are set from the EOS and the Lorentz factor is found from

W̄ = 1 + S2 (D̃h̄)2. (90)

The derivative is given by

f = W̄ ρ̄S2h̄ W̄D̃2h̄3, (91)


h̄ = ρ̄1P ρ . (92)

8.3 A note on the Roe and Marquina Riemann Solvers

Finding the Roe or Marquina fluxes as given is section 7 requires the left eigenvectors to either be supplied analytically or calculated numerically.

When this is done by inverting the matrix of right eigenvectors, in the actual code this is combined with the calculation of, e.g., the characteristic jumps Δw. Normally the eigenvalues and vectors are ordered lexicographically. However for the polytropic equation of state one of the equations is redundant, so the matrix formed by these eigenvectors is linearly dependent and hence singular. It turns out that this is only a minor problem; by rearranging the order of the eigenvalues and vectors it is possible to numerically invert the matrix. This means that no specific ordering of the eigenvalues should be assumed. It also explains the slightly strange ordering in the routines GRHydro_EigenProblem*.F90.

The current default is that the left eigenvectors are calculated analytically - for the expressions see Font [9]. For both the Roe and the Marquina solvers an optimized version of the flux calculation has been implemented. For more details on the analytical form and the optimized flux calculation see Ibáñez et al. [12], Aloy et al. [2] and Frieben et al. [10].

8.4 The atmosphere

In simulations of compact objects, often the matter is located only on a (small) portion of the numerical grid. In fact, over much of the evolved domain the physical situation is likely to be sufficiently well approximated by vacuum. However, in the vacuum limit the continuity equations describing the fluid break down. The speed of sound tends to the speed of light and everything fails (especially the conversion from conserved to primitive variables).

To avoid this problem it is customary to introduce an atmosphere. In our implementation, this is a low-density region surrounding the compact objects and initially it has no velocity and is in equilibrium. The introduction of an atmosphere is managed by the initial data thorns.

However GRHydro itself also knows about the atmosphere, of course. If the conserved variables D or τ are beneath some minimum value, or an evolution step might push them below such a value, then the relevant cell is not evolved. Also, if the density should fall below a minimum value in the routine that converts from conservative to primitive variables, all the variables are reset to the values adopted for the atmosphere.

Probably the hardest part of using GRHydro is to correctly set these atmosphere values. In the current implementation the atmosphere is used in three separate places. These are

  1. Set up of the initial data. Initial-data routines should set an atmosphere consistent with the one that will be evolved.
  2. In the routine that converts from conserved variables to primitive variables. This is where the majority of the atmosphere resets will occur.

    If the equation of state is polytropic then an attempt is made to convert to primitive variables. If the iterative algorithm returns a negative (and hence unphysical) value of ρ, then ρ is reset to the atmosphere value, the velocities are set to zero, and P, 𝜖, Si and τ are reset to be consistent with ρ (and D). Note that even though the polytropic equation of state gives us sufficient information to calculate a consistent value of D, this is not done.

    If the equation of state is the more general type (such as that of an ideal fluid) and if ρ is less than the specified minimum, then, as above, we assume we are in the atmosphere and call the routine that changes from the conserved to the primitive variables for the polytrope.

  3. When applying the update. If the calculated update terms for a specific cell would lead to either D or τ becoming negative, then two steps are taken. First, we do not update this specific cell. Second, the data in this cell is reset to be the atmosphere.

The reason why the routine that converts to the primitive variables does not ensure that D is consistent with the other variables is practical rather than accurate. If the value of the variables is set such that they all lie precisely on the atmosphere, then small errors (typically initially of the order of 1025 for a 643-point TOV star in octant symmetry) would move certain cells above the atmosphere values. Combined with the necessary atmosphere treatment this leads to high-frequency noise. This will lead to waves of matter falling onto the star. Despite their extremely low density (typically only an order of magnitude higher than the floor) they will result in visible secondary overtones in the oscillations of, e.g., the central density.

The parameters controlling the atmosphere are the following.

The motivation for these parameters referring only to the initial data is that it is sometimes best to set the initial atmosphere to slightly below the atmosphere cutoff used during evolution, as this avoids truncation error and metric evolution leading to low density waves travelling across the atmosphere.

The routines essential to the atmosphere are contained in GRHydro_Minima.F90, GRHydro_Con2Prim.F90, GRHydro_UpdateMask.F90.

8.5 Advection of passive scalars (’tracers’)

For some astrophysical problems it is necessary to advect passive compositional scalars such as the electron fraction Y e (number of electrons per baryon). For a generic tracer Xk, the evolution equation looks like

t(DXk) + xjf(j)(DX k) = 0, (93)

where D is the generalized particle number density as defined in Eq. (2). GRHydro currently supports any number of independent tracer variables. The following parameters have to be set to use the tracers:

Note, that your initial data thorn must set initial data for GRHydro::tracer[k] and GRHydro::cons_tracer[k] for all tracers you want to advect. GRHydro::cons_tracer[k] stores DXk.

8.5.1 Implementation and Limitations

9 History

The approximate time line is something like this:

This is necessarily only a sketch; many people have contributed to the history of this code, and the present authors were not around for most of it...

9.1 Thorn Source Code

This was initially written by Luca Baiotti, Ian Hawke and Pedro Montero with considerable assistance from Luciano Rezzolla, Toni Font, Nick Stergioulas and Ed Seidel. This led to the basic GRHydro thorns, GRHydro itself, GRHydro_Init_Data and GRHydro_RNSID.

Since then most of the maintenance has been done by Ian Hawke, Luca Baiotti and Frank Löffler. Various people have contributed to the development. In particular

9.2 Thorn Documentation

This documentation was first written largely by Ian Hawke and Scott Hawley in 2002. Long due, rather necessary and considerably large updates were made in 2008 by Luca Baiotti.

9.3 Acknowledgements

As already mentioned, the history behind this code leads to a long list of people to be acknowledged.

Firstly, without the work of the Valencia group this sort of code would be impossible.

Secondly, the incomparable work of Mark Miller and the Washington University - AEI Collaboration in producing the GR3D and GRAstro_Hydro codes gave an essential benchmark to aim for, and encouragement that it was possible!

Thirdly, the support of the Cactus team, especially Tom Goodale, Gabrielle Allen and Thomas Radke made life a lot easier.

Finally, for their work in coding, ideas and suggestions, or just plain encouragement, we would like to thank all at the AEI and in the EU Network, especially Toni Font, Luciano Rezzolla, Nick Stergioulas, Ed Seidel, Carsten Gundlach and José-Maria Ibáñez.

Originally Ed Seidel and then Luciano Rezzolla and Gabrielle Allen and many others have been granting (in addition to valuable scientific advice) financial support and human resources to the development of the code.


[1]   Aloy M.A., Ibánez J.M., MartíJ.M., Müller E. Astroph. J. Supp., 122: 151 (1999).

[2]   M. A. Aloy, J. A. Pons, and J. M. Ibáñez. An efficient implementation of flux formulae in multidimensional relativistic hydrodynamical codes. Comput. Phys. Commun., 120:115–121, 1999.

[3]   Banyuls F., Font J.A., Ibánez J.M., Martí J.M., Miralles J.A. Astrophys. J., 476: 221 (1997).

[4]   P. Colella and P. R. Woodward. The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations. J. Comput. Phys., 54, 174–201, 1984.

[5]   G. Cook Initial Data for Numerical Relativity Living Rev. Relativity, 3, 2000. [Article in on-line journal], cited on 31/08/02, Articles/Volume3/2000-5cook/index.html.

[6]   R. Donat and A. Marquina. Capturing shock reflections: An improved flux formula. J. Comput. Phys., 125:42–58, 1996.

[7]   Einfeldt B. Journal of Computational Physics, 25: 294 (1988).

[8]   Eulderink F., Mellema G. Astron. Astrophys., 284: 652 (1994).

[9]   J. A. Font. Numerical hydrodynamics in General Relativity. Living Rev. Relativity, 3, 2000. [Article in on-line journal], cited on 31/07/01, Articles/Volume3/2000-2font/index.html.

[10]   J. Frieben, J. M. Ibáñez, and J. Pons. in preparation

[11]   A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy. Uniformly high order accurate essentially non-oscillatory schemes, III. J. Comput. Phys., 71:231–303, 1987.

[12]   J. M. Ibáñez et al. in Godunov Methods: Theory and Applications. New York, 485–503, (2001)

[13]   A. Kurganov and E. Tadmor. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations. J. Comput. Phys., 160:241, 2000.

[14]   C. B. Laney. Computational Gasdynamics. Cambridge University Press, 1998.

[15]   R. J. LeVeque. Nonlinear conservation laws and finite volume methods for astrophysical fluid flow. In O. Steiner and A. Gautschy, editors, Computational Methods for Astrophysical Fluid Flow. Springer-Verlag, 1998.

[16]   J. M. Martí and E. Müller. Numerical hydrodynamics in Special Relativity. Living Rev. Relativity, 2, 1999. [Article in on-line journal], cited on 31/7/01,

[17]   J. J. Quirk. A contribution to the great Riemann solver debate. Int. J. Numer. Methods Fluids, 18:555–574, 1994.

[18]   Roe P.L. J. Comput. Phy., 43: 357 (1981).