The GRHydro code: three-dimensional relativistic hydrodynamics

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

Abstract

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 ﬁles, take a look at the Einstein Toolkit web page which is currently stored at

http://www.einsteintoolkit.org

GRHydro uses the hydro variables deﬁned 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 ﬁle, found at

GRHydro/test/GRHydro_test_shock.par

and will include the essential thorns

Initial data for shocks can be set using

GRHydro_Init_Data

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

TOVSolver

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 eﬃcient) 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 eﬃcient 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\left(\rho \right)$, and the more generic "General" type where the pressure is a function of the rest-mass density and the speciﬁc internal energy, $P=P\left(\rho ,𝜖\right)$. For the Polytype equations of state one fewer equation is evolved and the speciﬁc 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 http://www.einsteintoolkit.org.

2.2 Basic Usage

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

• Reconstruction methods:
• recon_method: The type of reconstruction method to use. tvd is the standard. ppm is more accurate but it requires more resources. eno gives in theory arbitrary order of accuracy but it is practically unworthy to go beyond ﬁfth order.
• tvd_limiter: When using tvd reconstruction the choice of limiter can have a large eﬀect. vanleerMC2 is probably the best to use, but the extremes of minmod and Superbee are also interesting.
• ppm_detect: When using ppm reconstruction with shocks, the shock detection algorithm can notably sharpen the proﬁle.
• Riemann solvers: Marquina is the standard solver used. HLLE is signiﬁcantly faster, but sometimes provides cruder approximation.
• Equations of state: These are controlled by the GRHydro_eos_type and GRHydro_eos_table parameters. Changing these parameters will depend on which equation of state thorns you have compiled in.
• Initial data parameters: GRHydro_rho_central is inherited by many initial data thorns to set the central density of compact ﬂuid objects such as single stars. However, this parameter is depreciated.
• Atmosphere parameters: Many of these are listed in section 8.4.

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

http://www.einsteintoolkit.org

It contains information on obtaining the code, together with thornlists and sample parameter ﬁles. For questions, send an email to the Einstein Toolkit mailing list users@einsteintoolkit.org.

3 Physical System

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

 ${\partial }_{t}q+{\partial }_{{x}^{i}}{f}^{\left(i\right)}\left(q\right)=s\left(q\right),$ (1)

where $q$ is a set of conserved variables, ${f}^{\left(i\right)}\left(q\right)$ the ﬂuxes and $s\left(q\right)$ the source terms. The ﬁve conserved variables are labeled $D$, ${S}_{i}$, and $\tau$, where $D$ is the generalized particle number density, ${S}_{i}$ are the generalized momenta in each direction, and $\tau$ is an internal energy term. These conserved variables are composed from a set of primitive variables, which are $\rho$, the rest-mass density, $p$, the pressure, ${v}^{i}$, the ﬂuid 3-velocities, $𝜖$, the speciﬁc internal energy, and $W$, the Lorentz factor, via the following relations

$\begin{array}{rcll}D& =& \sqrt{\gamma }W\rho & \text{}\\ {S}_{i}& =& \sqrt{\gamma }\rho h{W}^{2}{v}_{i}& \text{}\\ \tau & =& \sqrt{\gamma }\left(\rho h{W}^{2}-p\right)-D,& \text{(2)}\text{}\text{}\end{array}$

where $\gamma$ is the determinant of the spatial 3-metric ${\gamma }_{ij}$ and $h\equiv 1+𝜖+p∕\rho$. Only ﬁve of the primitive variables are independent. Usually the Lorentz factor is deﬁned in terms of the velocities and the metric as $W={\left(1-{\gamma }_{ij}{v}^{i}{v}^{j}\right)}^{-1∕2}$. Also one of the pressure, rest-mass density or speciﬁc internal energy terms is given in terms of the other two by an equation of state.

The ﬂuxes are usually deﬁned in terms of both the conserved variables and the primitive variables:

$\begin{array}{rcll}{F}^{i}\left(U\right)& =& {\left[D\left(\alpha {v}^{i}-{\beta }^{i}\right),{S}_{j}\left(\alpha {v}^{i}-{\beta }^{i}\right)+p{\delta }_{j}^{i},\tau \left(\alpha {v}^{i}-{\beta }^{i}\right)+p{v}^{i}\right]}^{T}\phantom{\rule{1em}{0ex}}.& \text{(3)}\text{}\text{}\end{array}$

The source terms are

$\begin{array}{rcll}s\left(U\right)=\left[\right0,{T}^{\mu \nu }\left(\right{\partial }_{\mu }{g}_{\nu j}+{\Gamma }_{\mu \nu }^{\delta }{g}_{\delta j}\left)\right,\alpha \left(\right{T}^{\mu 0}{\partial }_{\mu }ln\alpha -{T}^{\mu \nu }{\Gamma }_{\nu \mu }^{0}\left)\right\left]\right\phantom{\rule{1em}{0ex}}.& & & \text{(4)}\text{}\text{}\end{array}$

Note that the source terms do not contain diﬀerential 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 ﬂat spacetime. This conservative form of the relativistic Euler equations was ﬁrst 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].

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 diﬀerential equation such as (1) into an ordinary diﬀerential equation. For the GRHydro code the following steps are required.

• Partition the domain of interest into cells. For simplicity we shall assume a regular Cartesian partitioning. This is not necessary for the method of lines, but it is for GRHydro.
• Over a given cell with Cartesian coordinates $\left({x}_{i}^{1},{x}_{j}^{2},{x}_{k}^{3}\right)$, integrate equation (1) in space to ﬁnd the ordinary diﬀerential equation $\begin{array}{rcll}\frac{d\phantom{\rule{0.3em}{0ex}}}{dt}q& =& \int \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\int \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\int s\phantom{\rule{0.3em}{0ex}}{d}^{3}x+{\int }_{{x}_{j-1∕2}^{2}}^{{x}_{j+1∕2}^{2}}{\int }_{{x}_{k-1∕2}^{3}}^{{x}_{k+1∕2}^{3}}{f}^{\left(1\right)}\left(q\left({x}_{i-1∕2}^{1},y,z\right)\right)dy\phantom{\rule{0.3em}{0ex}}dz-& \text{}\\ & & {\int }_{{x}_{j-1∕2}^{2}}^{{x}_{j+1∕2}^{2}}{\int }_{{x}_{k-1∕2}^{3}}^{{x}_{k+1∕2}^{3}}{f}^{\left(1\right)}\left(q\left({x}_{i+1∕2}^{1},y,z\right)\right)dy\phantom{\rule{0.3em}{0ex}}dz+\cdots & \text{(5)}\text{}\text{}\end{array}$

where the boundaries of the Cartesian cells are given by ${x}_{i±1∕2}^{1}$ and so on.

• If we deﬁne $\stackrel{̄}{q}$ as the integral average of $q$ over the cell, after dividing (5) by the volume of the cell, we get an ordinary diﬀerential equation for $\stackrel{̄}{q}$, in terms of the ﬂux functions and the source terms as functions of the spatial diﬀerential of $\stackrel{̄}{q}$. We note that, unlike the spatial diﬀerential of $q$, the spatial diﬀerential of $\stackrel{̄}{q}$ is well deﬁned in a cell containing a discontinuity.

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

$\begin{array}{rcll}L\left(q\right)& =& \int \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\int \phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}\int s\phantom{\rule{0.3em}{0ex}}{d}^{3}x+{\int }_{{x}_{j-1∕2}^{2}}^{{x}_{j+1∕2}^{2}}{\int }_{{x}_{k-1∕2}^{3}}^{{x}_{k+1∕2}^{3}}{f}^{\left(1\right)}\left(q\left({x}_{i-1∕2}^{1},y,z\right)\right)dy\phantom{\rule{0.3em}{0ex}}dz-& \text{}\\ & & {\int }_{{x}_{j-1∕2}^{2}}^{{x}_{j+1∕2}^{2}}{\int }_{{x}_{k-1∕2}^{3}}^{{x}_{k+1∕2}^{3}}{f}^{\left(1\right)}\left(q\left({x}_{i+1∕2}^{1},y,z\right)\right)dy\phantom{\rule{0.3em}{0ex}}dz+\cdots & \text{(6)}\text{}\text{}\end{array}$

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 diﬀerential equation (5) will be ${N}^{th}$-order accurate provided that the time integrator used by MoL is ${N}^{th}$-order accurate in time, and that the discrete operator $L$ is ${N}^{th}$-order accurate in space and ﬁrst-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 simpliﬁed considerably by approximating the integrals by the midpoint rule to get

 $L\left(q\right)={s}_{i,j,k}+{f}_{i-1∕2,j,k}^{\left(1\right)}-{f}_{i+1∕2,j,k}^{\left(1\right)}+\cdots$ (7)

where the dependence of the ﬂux on $q$ and spatial position is implicit in the notation. Given this simpliﬁcation, the calculation of the right hand side operator splits simply into the following two parts:

1. Calculate the source terms $s\left(q\left({x}_{i}^{1},{x}_{j}^{2},{x}_{k}^{3}\right)\right)$:

Given that $\stackrel{̄}{q}$ is a second-order accurate approximation to $q$ at the midpoint of the cell, which is precisely $\left({x}_{i}^{1},{x}_{j}^{2},{x}_{k}^{3}\right)$, for second-order accuracy it is suﬃcient to use $s\left({\stackrel{̄}{q}}_{i,j,k}\right)$.

2. In each coordinate direction ${x}^{a}$, calculate the intercell ﬂux ${f}^{\left(a\right)}\left({q}_{i+1∕2,j,k}\right)$:

From the initial data $\stackrel{̄}{q}$ given at time ${t}^{n}$ we can reconstruct the data at the cell boundary, $\left({q}_{i+1∕2,j,k}\right)$ to any required accuracy in space (except in the vicinity of a shock, where only ﬁrst-order accuracy is guaranteed). However this will only be zeroth-order accurate in time. To ensure ﬁrst-order accuracy in time, we have to ﬁnd $\left({q}_{i+1∕2,j,k}\right)\left(t\right)$ 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 $\left({q}_{i+1∕2,j,k}\right)$, which we call ${q}_{L}$ and ${q}_{R}$, where ${q}_{L}$ is obtained from cell $i$ (left cell) and ${q}_{R}$ from cell $i+1$ (right cell).
2. The values ${q}_{L,R}$ are used as initial data for the Riemann problem. This is the initial value problem given by the partial diﬀerential equation (1) with semi-inﬁnite piecewise constant initial data ${q}_{L,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 ﬁrst-order accurate in time and retain the high order of accuracy in space which, as explained in the documentation to the MoL thorn, is suﬃcient to ensure that the method as a whole has a high order of accuracy.

So, the diﬃcult 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 $\stackrel{̄}{q}$, and another that calculates the intercell ﬂux $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 proﬁle of the function. Also, the calculation of the intercell ﬂux from the Riemann problem can be made to be “exactly correct”. That is, even though as explained above it may not be the true ﬂux for the real function $q$, it will be the exact physical solution for the values ${q}_{L,R}$ given by the reconstruction routine, so the intercell ﬂux cannot violate conservation or introduce oscillations. Unphysical eﬀects 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 ﬁrst-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-ﬁrst-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 diﬀerencing is a linear second-order method and is oscillatory near a shock. Instead we must ﬁnd a reconstruction method that depends on the data $q$ being reconstructed.

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

 ${q}^{First}\left(x\right)={\stackrel{̄}{q}}_{i},\phantom{\rule{2em}{0ex}}x\in \left[{x}_{i-1∕2},{x}_{i+1∕2}\right],$ (8)

and that any second-order conservative method can be written in terms of a slope or rather diﬀerence ${\Delta }_{i}$ as

 ${q}^{Second}\left(x\right)={\stackrel{̄}{q}}_{i}+\frac{x-{x}_{i}}{{x}_{i+1∕2}-{x}_{i-1∕2}}{\Delta }_{i},\phantom{\rule{2em}{0ex}}x\in \left[{x}_{i-1∕2},{x}_{i+1∕2}\right].$ (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 ﬁrst 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,

 ${q}^{TVD}\left(x\right)=\varphi \left(q\right){q}^{Second}+\left(1-\varphi \left(q\right)\right){q}^{First},$ (10)

where $\varphi \in \left[0,1\right]$ 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 $\varphi$ (having the same attributes as above) directly multiplies the slope, giving

 ${q}^{TVD}\left(x\right)={\stackrel{̄}{q}}_{i}+\frac{x-{x}_{i}}{{x}_{i+1∕2}-{x}_{i-1∕2}}\varphi \left(q\right){\Delta }_{i},\phantom{\rule{2em}{0ex}}x\in \left[{x}_{i-1∕2},{x}_{i+1∕2}\right].$ (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 diﬀusive 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 speciﬁed 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 deﬁne the minmod function:

 (12)

For reconstructing $q$ we choose two diﬀerences

 $\begin{array}{ccc}{\Delta }_{upw}\hfill & \hfill =\hfill & {q}_{i}-{q}_{i-1}\hfill \\ {\Delta }_{loc}\hfill & \hfill =\hfill & {q}_{i+1}-{q}_{i}\phantom{\rule{0.3em}{0ex}},\hfill \\ \hfill \end{array}$ (13)

and write

 ${q}^{TVD,minmod}\left(x\right)={\stackrel{̄}{q}}_{i}+\frac{x-{x}_{i}}{{x}_{i+1∕2}-{x}_{i-1∕2}}minmod\left({\Delta }_{upw},{\Delta }_{loc}\right),\phantom{\rule{2em}{0ex}}x\in \left[{x}_{i-1∕2},{x}_{i+1∕2}\right].$ (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 $\Delta x$. The ﬁrst step is to interpolate a quadratic polynomial to the cell boundary,

 ${q}_{i+1∕2}=\frac{1}{2}\left({q}_{i+1}+{q}_{i}\right)+\frac{1}{6}\left({\delta }_{m}{q}_{i}-{\delta }_{m}{q}_{i+1}\right),$ (15)

where

 (16)

and

 $\delta {q}_{i}=\frac{1}{2}\left({q}_{i+1}-{q}_{i-1}\right).$ (17)

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

 ${q}_{i}^{R}={q}_{i+1}^{L}={q}_{1+1∕2}.$ (18)

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

$\begin{array}{rcll}{q}_{i}^{L}={q}_{i}^{R}={q}_{i}& \mathrm{\text{if}}& \left({q}_{i}^{R}-{q}_{i}\right)\left({q}_{i}-{q}_{i}^{L}\right)\le 0& \text{(19)}\text{}\text{}\\ {q}_{i}^{L}=3{q}_{i}-2{q}_{i}^{R}& \mathrm{\text{if}}& \left({q}_{i}^{R}-{q}_{i}^{L}\right)\left({q}_{i}-\frac{1}{2}\left({q}_{i}^{L}+{q}_{i}^{R}\right)\right)>\frac{1}{6}{\left({q}_{i}^{R}-{q}_{i}^{L}\right)}^{2}& \text{(20)}\text{}\text{}\\ {q}_{i}^{R}=3{q}_{i}-2{q}_{i}^{L}& \mathrm{\text{if}}& \left({q}_{i}^{R}-{q}_{i}^{L}\right)\left({q}_{i}-\frac{1}{2}\left({q}_{i}^{L}+{q}_{i}^{R}\right)\right)<-\frac{1}{6}{\left({q}_{i}^{R}-{q}_{i}^{L}\right)}^{2}.& \text{(21)}\text{}\text{}\end{array}$

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

$\begin{array}{rcll}{\rho }_{i}^{L}& =& {\rho }_{i}^{L}\left(1-\eta \right)+\left({\rho }_{i-1}+\frac{1}{2}{\delta }_{m}{\rho }_{i-1}\right)\eta & \text{(22)}\text{}\text{}\\ {\rho }_{i}^{R}& =& {\rho }_{i}^{R}\left(1-\eta \right)+\left({\rho }_{i+1}-\frac{1}{2}{\delta }_{m}{\rho }_{i+1}\right)\eta & \text{(23)}\text{}\text{}\end{array}$

where $\eta$ is only applied if the discontinuity is mostly a contact (see [4] for the details) and is deﬁned as

 $\eta =\mathrm{\text{max}}\left(0,\mathrm{\text{min}}\left(1,{\eta }_{1}\left(\stackrel{̃}{\eta }-{\eta }_{2}\right)\right)\right),$ (24)

where ${\eta }_{1},{\eta }_{2}$ are positive constants and

with $𝜖$ another positive constant and

 ${\delta }^{2}{\rho }_{i}=\frac{{\rho }_{i+1}-2{\rho }_{i}+{\rho }_{i-1}}{6\Delta {x}^{2}}.$ (25)

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

 ${q}_{i}^{L,R}={\nu }_{i}{q}_{i}^{L,R}+\left(1-{\nu }_{i}\right){q}_{i},$ (26)

where

 (27)

and ${\omega }_{0},{\omega }_{1},{\omega }_{2}$ are constants.

The above ﬂattening 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 ﬂattening procedure is also implemented in GRHydro. Instead of 26, it consists in the formula

 ${q}_{i}^{L,R}={\stackrel{̃}{\nu }}_{i}{q}_{i}^{L,R}+\left(1-{\stackrel{̃}{\nu }}_{i}\right){q}_{i},$ (28)

where

$\begin{array}{rcll}{\stackrel{̃}{\nu }}_{i}& =& max\left(\right{\nu }_{i},{\nu }_{i+sign\left({p}_{i-1}-{p}_{i+1}\right)}\left)\right& \text{(29)}\text{}\text{}\end{array}$

and ${\nu }_{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 reﬁnement 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 et.al [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 diﬀerences 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 ${I}_{i}$. We then add to the stencil cell ${I}_{j}$, where $j=i±1$, where we choose $j$ to minimize the Newton divided diﬀerences

$\begin{array}{rcll}q\left[{x}_{i-1},{x}_{i}\right]& =& \frac{{q}_{i}-{q}_{i-1}}{{x}_{i+1∕2}-{x}_{i-3∕2}}& \text{(30)}\text{}\text{}\\ q\left[{x}_{i},{x}_{i+1}\right]& =& \frac{{q}_{i+1}-{q}_{i}}{{x}_{i+3∕2}-{x}_{i-1∕2}}.& \text{(31)}\text{}\text{}\end{array}$

We then recursively add more cells, minimizing the higher-order Newton divided diﬀerences $q\left[{x}_{i-k},\dots ,{x}_{i+j}\right]$ deﬁned by

 $q\left[{x}_{i-k},\dots ,{x}_{i+j}\right]=\frac{q\left[{x}_{i-k+1},\dots ,{x}_{i+j}\right]-q\left[{x}_{i-k},\dots ,{x}_{i+j-1}\right]}{{x}_{i+j}-{x}_{i-k}}.$ (32)

The reconstruction at the cell boundary is given by a standard ${k}^{\mathrm{\text{th}}}$-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\left(i\right)=\left\{{I}_{i-r},\dots ,{I}_{i+k-r-1}\right\},$ (33)

for some integer $r$, then there exist constants ${c}_{rj}$ depending only on the grid ${x}_{i}$ such that the boundary values for cell ${I}_{i}$ are given by

 ${q}_{i+1∕2}=\sum _{j=0}^{k-1}{c}_{rj}{q}_{i-r+j},\phantom{\rule{2em}{0ex}}{q}_{i-1∕2}=\sum _{j=0}^{k-1}{c}_{r-1,j}{q}_{i-r+j}.$ (34)

The constants ${c}_{rj}$ are given by the rather complicated formula

 ${c}_{rj}=\left\{\sum _{m=j+1}^{k}\frac{\sum _{l=0,l\ne m}^{k}\prod _{q=0,q\ne m,l}^{k}\left({x}_{i+1∕2}-{x}_{i-r+q-1∕2}\right)}{\prod _{l=0,l\ne m}^{k}\left({x}_{i-r+m-1∕2}-{x}_{i-r+L-1∕2}\right)}\right\}\Delta {x}_{i-r+j}.$ (35)

This simpliﬁes considerably if the grid is even. The coeﬃcients 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 ﬂux. The Riemann problem is speciﬁed by an equation in ﬂux conservative homogeneous form,

 ${\partial }_{t}q+{\partial }_{{x}^{i}}{f}^{\left(i\right)}\left(q\right)=0$ (36)

with piecewise constant initial data ${q}_{{}_{L,R}}$ separated by a discontinuity at ${x}^{\left(1\right)}=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 $\xi ={x}^{\left(1\right)}∕t$. The solution is given in terms of waves which separate diﬀerent 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 $\xi$ for a rarefaction, but spreads across a ﬁnite 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

 ${\partial }_{t}q+A{\partial }_{x}q=0\phantom{\rule{1em}{0ex}},$ (37)

where $A$ is a $N×N$ matrix with constant coeﬃcients. We deﬁne the eigenvalues ${\lambda }^{j}$ with associated right and left eigenvectors ${r}^{j},{l}_{j}$, where the eigenvectors are normalized to each other (i.e., their dot product is ${\delta }_{j}^{i}$). We shall assume that the eigenvectors span the space. The characteristic variables ${w}_{i}$ are deﬁned by

 ${w}_{i}={l}_{i}\cdot q.$ (38)

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

 ${\partial }_{t}w+\Lambda {\partial }_{x}w=0,$ (39)

where $\Lambda$ is the matrix containing the eigenvalues ${\lambda }_{i}$ on the diagonals and zeros elsewhere. Hence each characteristic variable ${w}^{i}$ obeys the linear advection equation with velocity $a={\lambda }_{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 ${\lambda }_{1}\le \cdots \le {\lambda }_{N}$. We also deﬁne the diﬀerences in the characteristic variables $\Delta {w}_{i}={\left({w}_{i}\right)}_{L}-{\left({w}_{i}\right)}_{R}$ across the ${i}^{th}$ characteristic wave. These diﬀerences are single numbers (‘scalars’). We note that these diﬀerences can easily be found from the initial data using

 $\Delta {w}_{i}={l}_{i}\cdot \left({q}_{L}-{q}_{R}\right).$ (40)

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

 ${w}_{i}={\left({w}_{i}\right)}_{L}+\sum _{j=1}^{M}\Delta {w}_{j}{e}^{j}\phantom{\rule{1em}{0ex}}if\phantom{\rule{1em}{0ex}}{\lambda }_{M}<\xi <{\lambda }_{M+1},$ (41)

or

 $w={\left({w}_{i}\right)}_{R}-\sum _{j=M+1}^{N}\Delta {w}_{j}{e}^{j}\phantom{\rule{1em}{0ex}}if\phantom{\rule{1em}{0ex}}{\lambda }_{M}<\xi <{\lambda }_{M+1},$ (42)

or as the average

 ${w}_{i}=\frac{1}{2}\left({\left({w}_{i}\right)}_{L}+{\left({w}_{i}\right)}_{R}+\sum _{j=1}^{M}\Delta {w}_{j}{e}^{j}-\sum _{j=M+1}^{N}\Delta {w}_{j}{e}^{j}\right)\phantom{\rule{1em}{0ex}}if\phantom{\rule{1em}{0ex}}{\lambda }_{M}<\xi <{\lambda }_{M+1},$ (43)

where ${e}^{i}$ is the column vector ${\left({e}^{i}\right)}_{j}={\delta }_{j}^{i}$.

Converting back to the original variables $q$ we have the solution

 $q=\frac{1}{2}\left({q}_{L}+{q}_{R}+\sum _{i=1}^{M}\Delta {w}_{i}{r}^{i}-\sum _{i=M+1}^{N}\Delta {w}_{i}{r}^{i}\right)\phantom{\rule{1em}{0ex}}if\phantom{\rule{1em}{0ex}}{\lambda }_{M}<\xi <{\lambda }_{M+1}.$ (44)

In the case where we are only interested in the ﬂux along the characteristic $\xi =0$ we can write the solution in the simple form

 $f\left(q\right)=\frac{1}{2}\left(f\left({q}_{L}\right)+f\left({q}_{R}\right)-\sum _{i=1}^{N}|{\lambda }_{i}|\Delta {w}_{i}{r}^{i}\right).$ (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 ﬂux are often essential. This is especially true in higher dimensions (¿1), where the solution of the ordinary diﬀerential 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 diﬀerent 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 ${\xi }_{±}$ 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=\left\{\begin{array}{ccc}\hfill {q}_{L}& \hfill if\hfill & \xi <{\xi }_{-}\hfill \\ \hfill {q}_{\ast }& \hfill if\hfill & {\xi }_{-}<\xi <{\xi }_{+}\hfill \\ \hfill {q}_{R}& \hfill if\hfill & \xi >{\xi }_{+},\hfill \end{array}\right\$ (46)

and the intermediate state ${q}_{\ast }$ is given by

 ${q}_{\ast }=\frac{{\xi }_{+}{q}_{R}-{\xi }_{-}{q}_{L}-f\left({q}_{R}\right)+f\left({q}_{L}\right)}{{\xi }_{+}-{\xi }_{-}}.$ (47)

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

 $f\left(q\right)=\frac{{\stackrel{̂}{\xi }}_{+}f\left({q}_{L}\right)-{\stackrel{̂}{\xi }}_{-}f\left({q}_{R}\right)+{\stackrel{̂}{\xi }}_{+}{\stackrel{̂}{\xi }}_{-}\left({q}_{R}-{q}_{L}\right)}{{\stackrel{̂}{\xi }}_{+}-{\stackrel{̂}{\xi }}_{-}},$ (48)

where

 ${\stackrel{̂}{\xi }}_{-}=min\left(0,{\xi }_{-}\right),\phantom{\rule{1em}{0ex}}{\stackrel{̂}{\xi }}_{+}=max\left(0,{\xi }_{+}\right).$ (49)

Knowledge of the precise minimum and maximum characteristic velocities ${\xi }_{±}$ requires knowing the solution of the Riemann problem. Instead, the characteristic velocities are usually found from the eigenvalues of the Jacobian matrix $\partial f∕\partial q$ 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 $\left[{\xi }_{-},{\xi }_{+}\right]$.

If we set $\alpha =max\left(|{\xi }_{-}|,|{\xi }_{+}|\right)$ and replace the characteristic velocities ${\xi }_{±}$ with $±\alpha$, we ﬁnd the Lax–Friedrichs ﬂux (cf. also Tadmor’s semi-discrete scheme [13])

 $f\left(q\right)=\frac{1}{2}\left[f\left({q}_{L}\right)+f\left({q}_{R}\right)+\alpha \left({q}_{L}-{q}_{R}\right)\right].$ (50)

This is very diﬀusive, 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 $\partial f∕\partial q$ is linearized about some intermediate state. Then the conservation form reduces to the linear equation

 ${\partial }_{t}q+A{\partial }_{x}q=0,$ (51)

where $A$ is a constant coeﬃcient 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 ﬂux along the characteristic ray $\xi =0$ is

 $f\left(q\right)=\frac{1}{2}\left(f\left({q}_{L}\right)+f\left({q}_{R}\right)-\sum _{i=1}^{N}|{\lambda }_{i}|\Delta {w}_{i}{r}^{i}\right).$ (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 ﬂux:

1. $A\left({q}_{Roe}\right)\left({q}_{R}-{q}_{L}\right)=f\left({q}_{R}\right)-f\left({q}_{L}\right)$,
2. $A\left({q}_{Roe}\right)$ is diagonalizable with real eigenvalues,
3. $A\left({q}_{Roe}\right)\to \partial f∕\partial q$ smoothly as ${q}_{Roe}\to q$.

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

 ${q}_{Roe}=\frac{1}{2}\left({q}_{R}+{q}_{L}\right),$ (53)

which satisﬁes 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 ﬂux along the characteristic ray $\xi =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 ﬂuid velocity is equal to the speed of sound of the ﬂuid. In the context of Riemann problems, sonic points are found when the ray $\xi =0$ is within a rarefaction wave.

Firstly deﬁne the left $l\left({q}_{L,R}\right)$ and right $r\left({q}_{L,R}\right)$ eigenvectors and the eigenvalues $\lambda \left({q}_{L,R}\right)$ of the Jacobian matrix $\partial f∕\partial q$ evaluated at the left and right states. Next deﬁne left and right characteristic variables ${w}_{L,R}$ and ﬂuxes ${\varphi }_{L,R}$ by

 ${\left({w}_{i}\right)}_{L,R}={l}_{i}\left({q}_{L,R}\right)\cdot {q}_{L,R},\phantom{\rule{1em}{0ex}}{\left({\varphi }_{i}\right)}_{L,R}={l}_{i}\left({q}_{L,R}\right)\cdot f\left({q}_{L,R}\right).$ (54)

Then the algorithm chooses the correct-sided characteristic ﬂux if the eigenvalues ${\lambda }_{i}\left({q}_{L}\right)$, ${\lambda }_{i}\left({q}_{R}\right)$ have the same sign, and uses a Lax–Friedrichs type ﬂux if they change sign. In full, the algorithm is given in ﬁgure 1.

 $\begin{array}{c}For\phantom{\rule{1em}{0ex}}\phantom{\rule{0.3em}{0ex}}i=1,\dots ,N\phantom{\rule{1em}{0ex}}do\hfill \\ \phantom{\rule{2em}{0ex}}\begin{array}{c}If\phantom{\rule{1em}{0ex}}\phantom{\rule{0.3em}{0ex}}{\lambda }_{i}\left({q}_{L}\right){\lambda }_{i}\left({q}_{R}\right)>0\phantom{\rule{1em}{0ex}}\phantom{\rule{0.3em}{0ex}}then\hfill \\ \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\begin{array}{c}If\phantom{\rule{1em}{0ex}}\phantom{\rule{0.3em}{0ex}}{\lambda }_{i}\left({q}_{L}\right)>0\phantom{\rule{1em}{0ex}}\phantom{\rule{0.3em}{0ex}}then\hfill \\ \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\begin{array}{ccc}\hfill {\varphi }_{+}^{i}& \hfill =\hfill & {\varphi }_{L}^{i}\hfill \\ \hfill {\varphi }_{-}^{i}& \hfill =\hfill & 0\hfill \end{array}\hfill \\ else\hfill \\ \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\begin{array}{ccc}\hfill {\varphi }_{+}^{i}& \hfill =\hfill & 0\hfill \\ \hfill {\varphi }_{-}^{i}& \hfill =\hfill & {\varphi }_{R}^{i}\hfill \end{array}\hfill \\ endif\hfill \end{array}\hfill \\ else\hfill \\ \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\begin{array}{ccc}\hfill {\alpha }^{i}& \hfill =\hfill & max\left(|{\lambda }_{i}\left({q}_{L}\right),{\lambda }_{i}\left({q}_{R}\right)|\right)\hfill \\ \hfill {\varphi }_{+}^{i}& \hfill =\hfill & \frac{1}{2}\left({\varphi }_{L}^{i}+{\alpha }^{i}{w}_{L}^{i}\right)\hfill \\ \hfill {\varphi }_{-}^{i}& \hfill =\hfill & \frac{1}{2}\left({\varphi }_{R}^{i}-{\alpha }^{i}{w}_{R}^{i}\right)\hfill \end{array}\hfill \\ endif\hfill \\ \hfill \end{array}\hfill \\ enddo\hfill \end{array}$ (55)

Figure 1: The algorithm to calculate the Marquina ﬂux.

Then the numerical ﬂux is given by

 $f\left(q\right)=\sum _{i=1}^{N}\left[{\varphi }_{+}^{i}{r}^{i}\left({q}_{L}\right)+{\varphi }_{-}^{i}{r}^{i}\left({q}_{R}\right)\right].$ (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 ﬁrst, before the ﬂux terms (it simpliﬁes the loop very slightly). There are a few points to note about the calculation of the sources.

• The calculation used here, taken from the GR3D code, requires both the metric and the extrinsic curvature.
• In order to calculate the Christoﬀel symbols the gauge and metric variables must be diﬀerenced. Currently centred diﬀerencing of second or fourth (we are safe to use this, as GRHydro requires always at least 2 ghost zones) order is hardwired in. The two diﬀerencings can be selected via the parameter GRHydro::sources_spatial_order.
• For numerical reasons, namely in order to avoid the presence of time derivatives in the source-term computation, the implemented form of the source terms is not (4) directly, but it has been modiﬁed as shown in the following paragraphs (following a clever idea by Mark Miller (see the GR3D code)

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-Christoﬀel symbols ${\phantom{\rule{1em}{0ex}}}^{\left(4\right)}{\Gamma }_{\mu \nu }^{\rho }$ 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

 ${\partial }_{t}{\gamma }_{ij}=2\left(-\alpha {K}_{ij}+{\partial }_{\left(i}{\beta }_{j}{-}^{\left(3\right)}{\Gamma }_{ij}^{k}{\beta }_{k}\right)\phantom{\rule{1em}{0ex}}.$ (57)

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

$\begin{array}{rcll}{\nabla }_{i}{\gamma }^{jk}={\partial }_{i}{\gamma }^{jk}+{2}^{\left(3\right)}{\Gamma }_{il}^{j}{\gamma }^{lk}=0\phantom{\rule{1em}{0ex}}.& & & \text{(58)}\text{}\text{}\end{array}$

We start from the ${}^{\left(4\right)}{\Gamma }_{00}^{0}$ symbol:

$\begin{array}{rcl}{}^{\left(4\right)}{\Gamma }_{00}^{0}=\frac{1}{2{\alpha }^{2}}\left[\right-{\partial }_{t}\left(\right{\beta }_{k}{\beta }^{k}\left)\right+2\alpha {\partial }_{t}\alpha +2{\beta }^{i}{\partial }_{t}{\beta }_{i}-{\beta }^{i}{\partial }_{i}\left(\right{\beta }_{k}{\beta }^{k}\left)\right+2\alpha {\beta }^{i}{\partial }_{i}\alpha \left]\right& & \text{(59)}\text{}\text{}\end{array}$

and we expand the derivatives as

$\begin{array}{rcll}{\partial }_{t}\left(\right{\beta }_{k}{\beta }^{k}\left)\right& =& {\partial }_{t}\left(\right{\gamma }_{jk}{\beta }^{j}{\beta }^{k}\left)\right=2{\gamma }_{jk}{\beta }^{j}{\partial }_{t}{\beta }^{k}+{\beta }^{j}{\beta }^{k}{\partial }_{t}{\gamma }_{jk}=& \text{}\\ & =& 2{\beta }_{k}{\partial }_{t}{\beta }^{k}-2\alpha {K}_{jk}{\beta }^{j}{\beta }^{k}+2{\beta }^{j}{\beta }^{k}{\partial }_{j}{\beta }_{k}-{2}^{\left(3\right)}{\Gamma }_{kj}^{i}{\beta }_{i}{\beta }^{j}{\beta }^{k}& \text{(60)}\text{}\text{}\end{array}$

and

$\begin{array}{rcll}{\partial }_{i}\left(\right{\beta }_{k}{\beta }^{k}\left)\right={\partial }_{i}\left(\right{\gamma }^{jk}{\beta }_{j}{\beta }_{k}\left)\right=2{\gamma }^{jk}{\beta }_{j}{\partial }_{i}{\beta }_{k}+{\beta }_{j}{\beta }_{k}{\partial }_{i}{\gamma }^{jk}=2{\beta }_{k}{\partial }_{i}{\beta }_{k}-{2}^{\left(3\right)}{\Gamma }_{ik}^{j}{\beta }_{j}{\beta }^{k}\phantom{\rule{1em}{0ex}},& & & \text{(61)}\text{}\text{}\end{array}$

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

$\begin{array}{rcl}{}^{\left(4\right)}{\Gamma }_{00}^{0}=\frac{1}{\alpha }\left(\right{\partial }_{t}\alpha +{\beta }^{i}{\partial }_{i}\alpha +{K}_{jk}{\beta }^{j}{\beta }^{k}\left)\right\phantom{\rule{1em}{0ex}}.& & \text{(62)}\text{}\text{}\end{array}$

With the same strategy we then compute

$\begin{array}{rcl}{}^{\left(4\right)}{\Gamma }_{i0}^{0}=& -\frac{1}{2{\alpha }^{2}}\left[\right{\partial }_{i}\left({\beta }^{k}{\beta }_{k}-{\alpha }^{2}\right)-{\beta }^{j}\left({\partial }_{i}{\beta }_{j}-{\partial }_{j}{\beta }_{i}+{\partial }_{t}{\gamma }_{ij}\right)\left]\right=-\frac{1}{\alpha }\left(\right{\partial }_{i}\alpha -{\beta }^{j}{K}_{ij}\left)\right& \text{(63)}\text{}\text{}\end{array}$

and

$\begin{array}{rcl}{}^{\left(4\right)}{\Gamma }_{ij}^{0}=& -\frac{1}{2{\alpha }^{2}}\left[\right{\partial }_{i}{\beta }_{j}+{\partial }_{j}{\beta }_{i}-{\partial }_{t}{\gamma }_{ij}-{\beta }^{k}\left({\partial }_{i}{\gamma }_{kj}+{\partial }_{j}{\gamma }_{ki}-{\partial }_{k}{\gamma }_{ij}\right)\left]\right=-\frac{1}{\alpha }{K}_{ij}\phantom{\rule{1em}{0ex}}.& \text{(64)}\text{}\text{}\end{array}$

Other more straightforward calculations give

$\begin{array}{lllllll}\hfill {}^{\left(4\right)}{\Gamma }_{00j}& =\phantom{\rule{2em}{0ex}}& \hfill {}^{\left(4\right)}{\Gamma }_{0j}^{\nu }{g}_{\nu 0}& =\frac{1}{2}{\partial }_{j}\left({\beta }_{k}{\beta }^{k}-{\alpha }^{2}\right),\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}& \hfill \text{(65)}\\ \hfill & \phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}& \hfill \\ \hfill {}^{\left(4\right)}{\Gamma }_{l0j}& =\phantom{\rule{2em}{0ex}}& \hfill {}^{\left(4\right)}{\Gamma }_{lj}^{\nu }{g}_{\nu 0}& =\alpha {K}_{lj}+{\partial }_{l}{\beta }_{j}+{\partial }_{j}{\beta }_{l}-{\beta {}^{k}}_{}\left(3\right){\Gamma }_{lj}^{k}\phantom{\rule{1em}{0ex}},\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}& \hfill \text{(66)}\\ \hfill & \phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}& \hfill \\ \hfill {}^{\left(4\right)}{\Gamma }_{0lj}& =\phantom{\rule{2em}{0ex}}& \hfill {}^{\left(4\right)}{\Gamma }_{0j}^{\nu }{g}_{\nu l}& =-\alpha {K}_{jl}+{\partial }_{l}{\beta }_{j}-{\beta {}^{k}}_{}\left(3\right){\Gamma }_{lj}^{k}\phantom{\rule{1em}{0ex}},\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}& \hfill \text{(67)}\\ \hfill & \phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}& \hfill \\ \hfill {}^{\left(4\right)}{\Gamma }_{lmj}& =\phantom{\rule{2em}{0ex}}& \hfill {}^{\left(4\right)}{\Gamma }_{lj}^{\nu }{g}_{\nu m}& {=}^{\left(3\right)}{\Gamma }_{lmj}\phantom{\rule{1em}{0ex}},\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}& \hfill \text{(68)}\end{array}$

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

8.1.1 Source term for ${S}_{k}$

Now we have all the expressions for calculating the source terms. The ones for the variables ${S}_{\phantom{\rule{0.3em}{0ex}}k}$ are

 $\left(\right{\mathsc{𝒮}}_{{S}_{k}}{\left)\right}_{j}={T}_{\nu }^{\mu }{\Gamma }_{\mu j}^{\nu }={T}^{\mu \nu }{\Gamma }_{\mu \nu j}\phantom{\rule{1em}{0ex}}.$ (69)

After expanding the derivative in (65), the coeﬃcient of the ${T}^{\phantom{\rule{1em}{0ex}}00}$ term in (69) becomes

$\begin{array}{rcl}{}^{\left(4\right)}{\Gamma }_{00j}=& \frac{1}{2}{\beta }^{l}{\beta }^{m}{\partial }_{j}{\gamma }_{lm}-\alpha {\partial }_{j}\alpha +{\beta }_{m}{\partial }_{j}{\beta }^{m}.& \text{(70)}\text{}\text{}\end{array}$

The coeﬃcient of the ${T}^{\phantom{\rule{0.3em}{0ex}}0i}$ term is

$\begin{array}{rcl}{}^{\left(4\right)}{\Gamma }_{0ij}{+}^{\left(4\right)}{\Gamma }_{i0j}={\partial }_{j}{\beta }_{i}={\beta }^{l}{\partial }_{i}{\gamma }_{jl}+{\gamma }_{il}{\partial }_{j}{\beta }^{l}.& & \text{(71)}\text{}\text{}\end{array}$

The coeﬃcient of the ${T}^{\phantom{\rule{0.3em}{0ex}}lm}$ term is simply

$\begin{array}{rcl}{}^{\left(3\right)}{\Gamma }_{lmj}=\frac{1}{2}\left(\right{\partial }_{j}{\gamma }_{ml}+{\partial }_{m}{\gamma }_{jl}-{\partial }_{l}{\gamma }_{mj}\left)\right.& & \text{(72)}\text{}\text{}\end{array}$

Finally, summing (70)–(72) we ﬁnd

$\begin{array}{rcll}\left(\right{\mathsc{𝒮}}_{{S}_{k}}{\left)\right}_{j}& =& {T}^{00}\left(\frac{1}{2}{\beta }^{l}{\beta }^{m}{\partial }_{j}{\gamma }_{lm}-\alpha {\partial }_{j}\alpha \right)+{T}^{0i}{\beta }^{l}{\partial }_{j}{\gamma }_{il}+{T}_{i}^{0}{\partial }_{j}{\beta }^{i}+\frac{1}{2}{T}^{lm}{\partial }_{j}{\gamma }_{lm}\phantom{\rule{1em}{0ex}},& \text{(73)}\text{}\text{}\end{array}$

which is the expression implemented in the code.

8.1.2 Source term for $\tau$

The source term for $\tau$ is [cf. (4)]

 ${\mathsc{𝒮}}_{\tau }=\alpha \left({T}^{\mu 0}{\partial }_{\mu }\alpha -\alpha {T{}^{\mu \nu }}^{\left(4\right)}{\Gamma }_{\mu \nu }^{0}\right).$ (74)

For clarity, again we consider separately the terms containing as a factor the diﬀerent components of ${T}^{\mu \nu }$. From (62) we ﬁnd the coeﬃcient of ${T}^{\phantom{\rule{0.3em}{0ex}}00}$ to be

$\begin{array}{rcll}\alpha \left(\right{\partial }_{t}\alpha -{\alpha }^{\left(4\right)}{\Gamma }_{00}^{0}\left)\right=-\alpha \left(\right{\beta }^{i}{\partial }_{i}\alpha +{\beta }^{k}{\beta }^{l}{K}_{kl}\left)\right\phantom{\rule{1em}{0ex}}.& & & \text{(75)}\text{}\text{}\end{array}$

The coeﬃcient of the term ${T}^{\phantom{\rule{0.3em}{0ex}}0i}$ is given by

$\begin{array}{rcll}\alpha \left(\right{\partial }_{i}\alpha -2{\alpha }^{\left(4\right)}{\Gamma }_{i0}^{0}\left)\right=2\alpha {\beta }^{j}{K}_{ij}-\alpha {\partial }_{i}\alpha & & & \text{(76)}\text{}\text{}\end{array}$

and, ﬁnally, the coeﬃcient for ${T}^{\phantom{\rule{0.3em}{0ex}}ij}$ is

$\begin{array}{rcll}-{\alpha {}^{2}}^{\left(4\right)}{\Gamma }_{ij}^{0}=\alpha {K}_{ij}\phantom{\rule{1em}{0ex}}.& & & \text{(77)}\text{}\text{}\end{array}$

The ﬁnal expression implemented in the code is thus

$\begin{array}{rcll}{\mathsc{𝒮}}_{\tau }=\alpha \left[\right{T}^{00}\left({\beta }^{i}{\beta }^{j}{K}_{ij}-{\beta }^{i}{\partial }_{i}\alpha \right)+{T}^{0i}\left(-{\partial }_{i}\alpha +2{\beta }^{j}{K}_{ij}\right)+{T}^{ij}{K}_{ij}\left]\right\phantom{\rule{1em}{0ex}}.& & & \text{(78)}\text{}\text{}\end{array}$

8.2 Conversion from conservative to primitive variables

As noted in section 3 the variables that are evolved are the conserved variables $D,{S}_{i},\tau$. But in order to calculate the ﬂuxes and sources we require the primitive variables $\rho ,{v}_{i},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 ﬁnd a root of the pressure equation. Given an initial guess for the pressure $\stackrel{̄}{P}$ we ﬁnd the root of the function

 $f=\stackrel{̄}{P}-P\left(\stackrel{̄}{\rho },\stackrel{̄}{𝜖}\right),$ (79)

where the approximate density and speciﬁc internal energy are given by

$\begin{array}{rcll}\stackrel{̄}{\rho }& =& \frac{\stackrel{̃}{D}}{\stackrel{̃}{\tau }+\stackrel{̄}{P}+\stackrel{̃}{D}}\sqrt{{\left(\stackrel{̃}{\tau }+\stackrel{̄}{P}+\stackrel{̃}{D}\right)}^{2}-{S}^{2}},& \text{(80)}\text{}\text{}\\ \stackrel{̄}{W}& =& \frac{\stackrel{̃}{\tau }+\stackrel{̄}{P}+\stackrel{̃}{D}}{\sqrt{{\left(\stackrel{̃}{\tau }+\stackrel{̄}{P}+\stackrel{̃}{D}\right)}^{2}-{S}^{2}}},& \text{(81)}\text{}\text{}\\ \stackrel{̄}{𝜖}& =& {\stackrel{̃}{D}}^{-1}\left(\sqrt{{\left(\stackrel{̃}{\tau }+\stackrel{̄}{P}+\stackrel{̃}{D}\right)}^{2}-{S}^{2}}-\stackrel{̄}{P}\stackrel{̄}{W}-\stackrel{̃}{D}\right).& \text{(82)}\text{}\text{}\end{array}$

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

 $\stackrel{̃}{D}={\gamma }^{-1∕2}D,$ (83)

where $\gamma$ is the determinant of the 3-metric, and ${S}^{2}$ is given by

 ${S}^{2}={\gamma }^{ij}{\stackrel{̃}{S}}_{i}{\stackrel{̃}{S}}_{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}^{\prime }=1-\frac{\partial P}{\partial \rho }\frac{\partial \rho }{\partial P}-\frac{\partial P}{\partial 𝜖}\frac{\partial 𝜖}{\partial P},$ (85)

where $\frac{\partial P}{\partial \rho }$ and $\frac{\partial P}{\partial 𝜖}$ given by calls to EOS_Base, and

$\begin{array}{rcll}\frac{\partial \rho }{\partial P}& =& \frac{\stackrel{̃}{D}{S}^{2}}{\sqrt{{\left(\stackrel{̃}{\tau }+\stackrel{̄}{P}+\stackrel{̃}{D}\right)}^{2}-{S}^{2}}{\left(\stackrel{̃}{\tau }+\stackrel{̄}{P}+\stackrel{̃}{D}\right)}^{2}},& \text{(86)}\text{}\text{}\\ \frac{\partial 𝜖}{\partial P}& =& \frac{\stackrel{̄}{P}{S}^{2}}{\rho \left({\left(\stackrel{̃}{\tau }+\stackrel{̄}{P}+\stackrel{̃}{D}\right)}^{2}-{S}^{2}\right)\left(\stackrel{̃}{\tau }+\stackrel{̄}{P}+\stackrel{̃}{D}\right)}.& \text{(87)}\text{}\text{}\\ & & & \text{(88)}\text{}\text{}\end{array}$

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

 $f=\stackrel{̄}{\rho }\stackrel{̄}{W}-\stackrel{̃}{D},$ (89)

where $\stackrel{̄}{\rho }$ is the variable solved for, the pressure, speciﬁc internal energy and enthalpy $\stackrel{̄}{h}$ are set from the EOS and the Lorentz factor is found from

 $\stackrel{̄}{W}=\sqrt{1+\frac{{S}^{2}}{{\left(\stackrel{̃}{D}\stackrel{̄}{h}\right)}^{2}}}.$ (90)

The derivative is given by

 ${f}^{\prime }=\stackrel{̄}{W}-\frac{\stackrel{̄}{\rho }{S}^{2}{\stackrel{̄}{h}}^{\prime }}{\stackrel{̄}{W}{\stackrel{̃}{D}}^{2}{\stackrel{̄}{h}}^{3}},$ (91)

where

 ${\stackrel{̄}{h}}^{\prime }={\stackrel{̄}{\rho }}^{-1}\frac{\partial P}{\partial \rho }.$ (92)

8.3 A note on the Roe and Marquina Riemann Solvers

Finding the Roe or Marquina ﬂuxes 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 $\Delta 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 speciﬁc 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 ﬂux calculation has been implemented. For more details on the analytical form and the optimized ﬂux 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 suﬃciently well approximated by vacuum. However, in the vacuum limit the continuity equations describing the ﬂuid 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 $\tau$ 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 $\rho$, then $\rho$ is reset to the atmosphere value, the velocities are set to zero, and $P$, $𝜖$, ${S}_{i}$ and $\tau$ are reset to be consistent with $\rho$ (and $D$). Note that even though the polytropic equation of state gives us suﬃcient 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 ﬂuid) and if $\rho$ is less than the speciﬁed 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 speciﬁc cell would lead to either $D$ or $\tau$ becoming negative, then two steps are taken. First, we do not update this speciﬁc 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 $1{0}^{-25}$ for a $6{4}^{3}$-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 ﬂoor) they will result in visible secondary overtones in the oscillations of, e.g., the central density.

The parameters controlling the atmosphere are the following.

• GRHydro::rho_abs_min. An absolute value for $\rho$ in the atmosphere. Defaults to -1. Any negative value will be ignored, and the value of rho_rel_min used instead.
• GRHydro::rho_rel_min. A relative value for $\rho$ in the atmosphere. Defaults to $1{0}^{-7}$. Only used if rho_abs_min is negative, which is the default behaviour. The actual value of the atmosphere will be $\rho =$rho_rel_min$×$GRHydro_rho_max, where GRHydro::GRHydro_rho_max is a variable containing the maximum value of $\rho$ on the numerical grid at time zero.
• initial_rho_abs_min. An absolute value for rho in the initial atmosphere. It is used only by initial data routines. Unused if negative.
• initial_rho_rel_min. A relative (to the initial maximum rest-mass density) value for rho in the atmosphere. It is used only by initial data routines and only if it is positive and initial_rho_abs_min is negative.
• initial_atmosphere_factor. A relative (to the initial atmosphere) value for rho in the atmosphere. It is used only by initial data routines. It multiplies the atmosphere value used by the initial data solver. Unused if negative.
• GRHydro_atmo_tolerance. A parameter useful mostly in mesh-reﬁned simulations. A point is set to the atmosphere values in the conservative to primitive routines if its rest-mass density is such that $\rho <$ GRHydro_rho_min$\ast \left(1+$GRHydro_atmo_tolerance). This avoids occasional spurious oscillations in (Carpet) buﬀer zones lying in the atmosphere (because prolongation happens on conserved variables).

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 cutoﬀ 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 ${X}_{k}$, the evolution equation looks like

 ${\partial }_{t}\left(D{X}_{k}\right)+{\partial }_{{x}^{j}}{f}^{\left(j\right)}\left(D{X}_{k}\right)=0\phantom{\rule{0.3em}{0ex}},$ (93)

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

• GRHydro::evolve_tracer. Boolean. Set to yes if you want the tracers to be active.
• GRHydro::number_of_tracers. Integer. Defaults to 0. To use tracers, set to at least 1.

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 $D{X}_{k}$.

8.5.1 Implementation and Limitations
• Reconstruction: Currently only TVD and PPM reconstruction of the tracers ${X}_{k}$ are implemented. Since for most astrophysical problems one will associate the tracer with some compositional variable it might be better to reconstruct $\rho {X}_{k}$.
• Riemann Solvers: Only HLLE and Marquina are supported. In HLLE, the ﬂuxes as given in Eq. (93) are computed for each tracer. In the Marquina solver, we multiply the $D$-ﬂux by each tracer to get the individual tracer ﬂuxes according to the following prescription:  $\begin{array}{c}If\phantom{\rule{1em}{0ex}}\phantom{\rule{0.3em}{0ex}}{F}_{i+1∕2}\left(D\right)>0\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}then\hfill \\ \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\begin{array}{c}\phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}{F}_{i+1∕2}\left(D{X}_{k}\right)={X}_{{k}_{i}}^{+}{F}_{i+1∕2}\left(D\right)\hfill \end{array}\hfill \\ else\hfill \\ \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}\begin{array}{c}\hfill \phantom{\rule{2em}{0ex}}\phantom{\rule{2em}{0ex}}{F}_{i+1∕2}\left(D{X}_{k}\right)={X}_{{k}_{i+1}}^{-}{F}_{i+1∕2}\left(D\right)\end{array}\hfill \\ endif\hfill \\ \hfill \end{array}$ (94)

The above was suggested by Miguel Aloy and ﬁrst implemented and tested by Harry Dimmelmeier in his code (CoCoNuT), then by Christian D. Ott in GRHydro.

9 History

The approximate time line is something like this:

•  1995: Valencia group hydrodynamics code, ﬁxed spacetime.
•  1997: Ported to Cactus, extensively rewritten for the Binary Neutron Star Grand Challenge. Primarily written by Mark Miller. Released as GR3D as public domain code.
•  1998-: Developed as Cactus thorn MAHC inside the GRAstro_Hydro arrangement at Washington University, primarily by Mark Miller.
• 2002-2008: Whisky written based on GR3D.
• 2008-: GRHydro based on the public version of Whisky

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öﬄer. Various people have contributed to the development. In particular

• The PPM reconstruction methods were written by Ian Hawke heavily based on the code of Toni Font. They were later expanded by Christian D. Ott and Luca Baiotti.
• The Roe and Marquina solvers were made considerably more eﬃcient thanks to Joachim Frieben.
• The current atmosphere algorithm is a mixture of ideas from the GR3D code, Luciano Rezzolla, Toni Font and Nick Stergioulas. The current setup was written by Ian Hawke and Luca Baiotti.
• The 1-dimensional TOV solver GRHydro_TOVSolver was written by Ian Hawke based on a short paper by Thomas Baumgarte.

9.2 Thorn Documentation

This documentation was ﬁrst 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 scientiﬁc advice) ﬁnancial support and human resources to the development of the code.

References

[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 eﬃcient implementation of ﬂux 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, http://www.livingreviews.org/ Articles/Volume3/2000-5cook/index.html.

[6]   R. Donat and A. Marquina. Capturing shock reﬂections: An improved ﬂux 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, http://www.livingreviews.org/ 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-diﬀusion 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 ﬁnite volume methods for astrophysical ﬂuid ﬂow. 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, http://www.livingreviews.org/Articles/Volume2/1999-3marti/index.html.

[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).

[19]   C. Shu. High Order ENO and WENO Schemes for Computational Fluid Dynamics. In T. J. Barth and H. Deconinck, editors High-Order Methods for Computational Physics. Springer, 1999. A related on-line version can be found under Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws at http://www.icase.edu/library/reports/rdp/97/97-65RDP.tex.refer.html.

[20]   E. F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer-Verlag, 2nd edition, 1999.

10 Parameters

 constrain_to_1d Scope: private BOOLEAN Description: Set ﬂuid velocities to zero for non-radial motion Default: no

 use_cxx_code Scope: private BOOLEAN Description: Use experimental C++ code? Default: yes

 verbose Scope: private BOOLEAN Description: Some debug output Default: no

 apply_h_viscosity Scope: restricted BOOLEAN Description: H viscosity is useful to ﬁx the carbuncle instability seen at strong shocks Default: no

 atmo_falloﬀ_power Scope: restricted REAL Description: The power at which the atmosphere level falls oﬀ as (atmo_falloof_radius/r)**atmo_falloﬀ_power Range Default: 0.0 0:* Anything positive

 atmo_falloﬀ_radius Scope: restricted REAL Description: The radius for which the atmosphere starts to falloﬀ as (atmo_falloﬀ_radius/r)**atmo_falloﬀ_power Range Default: 50.0 0:* Anything positive

 atmo_tolerance_power Scope: restricted REAL Description: The power at which the atmosphere tolerance increases as (r/atmo_tolerance_radius)**atmo_tolerance_power Range Default: 0.0 0:* Anything positive

 atmo_tolerance_radius Scope: restricted REAL Description: The radius for which the atmosphere tolerance starts to increase as (r/atmo_tolerance_radius)**atmo_tolerance_power Range Default: 50.0 0:* Anything positive

 avec_gauge Scope: restricted KEYWORD Description: Which gauge condition to use when evolving the vector potential Range Default: lorenz algebraic Algebraic gauge lorenz Lorenz gauge

 bound Scope: restricted KEYWORD Description: Which boundary condition to use - FIXME Range Default: none ﬂat Zero order extrapolation none None static Static, no longer supported scalar Constant

 c2p_reset_pressure Scope: restricted BOOLEAN Description: If the pressure guess is unphysical should we recompute it? Default: no

 c2p_reset_pressure_to_value Scope: restricted REAL Description: The value to which the pressure is reset to when a failure occurrs in C2P Range Default: 1.e-20 0: greater than zero

 c2p_resort_to_bisection Scope: restricted BOOLEAN Description: If the pressure guess is unphysical, should we try with bisection (slower, but more robust) Default: no

 calculate_bcom Scope: restricted BOOLEAN Description: Calculate the comoving contravariant magnetic 4-vector bâ? Default: no

 check_for_trivial_rp Scope: restricted BOOLEAN Description: Do check for trivial Riemann Problem Default: yes

 check_rho_minimum Scope: restricted BOOLEAN Description: Should a check on rho < GRHydro_rho_min be performed and written as WARNING level 2? Default: no

 clean_divergence Scope: restricted BOOLEAN Description: Use hyperbolic divergence cleaning Default: no

 comoving_attenuate Scope: restricted KEYWORD Description: Which attenuation method for the comoving shift Range Default: tanh sqrt Multiply by sqrt(rho/rho_max) tanh ”Multiply by 1/2*(1+tanh(factor*( rho/rho_max - oﬀset)))”

 comoving_factor Scope: restricted REAL Description: Factor multiplying the velocity for the comoving shift Range Default: 0.0 0.0:2.0 [0,2] is allowed, but [0,1] is probably reasonable

 comoving_tanh_factor Scope: restricted REAL Description: The factor in the above tanh attenuation Range Default: 10.0 (0.0:* Any positive number

 comoving_tanh_oﬀset Scope: restricted REAL Description: The oﬀset in the above tanh attenuation Range Default: 0.05 0.0:1.0 Only makes sense in [0,1]

 comoving_v_method Scope: restricted KEYWORD Description: Which method for getting the radial velocity Range Default: projected projected vr = x_i . vî / r components vr = sqrt(v_i vî)

 con2prim_oct_hack Scope: restricted BOOLEAN Description: Disregard c2p failures in oct/rotsym90 boundary cells with xyz < 0 Default: no

 decouple_normal_bﬁeld Scope: restricted BOOLEAN Description: when using divergence cleaning properly decouple Bx,psidc subsystem Default: yes

 enhanced_ppm_c2 Scope: restricted REAL Description: Parameter for enhancecd ppm limiter. Default from McCorquodale & Colella 2011 Range Default: 1.25 *:* must be greater than 1. According to Colella&Sekora 2008, enhanced ppm is insensitive to C in [1.25,5]

 enhanced_ppm_c3 Scope: restricted REAL Description: Parameter for enhancecd ppm limiter. Default from McCorquodale & Colella 2011 Range Default: 0.1 0:* must be greater than 0.

 eno_order Scope: restricted INT Description: The order of accuracy of the ENO reconstruction Range Default: 2 1:* Anything above 1, but above 5 is pointless

 eos_change Scope: restricted BOOLEAN Description: Recalculate ﬂuid quantities if changing the EoS Default: no

 eos_change_type Scope: restricted KEYWORD Description: Change polytropic K or Gamma? Range Default: Gamma K Change the K Gamma Change the Gamma GammaKS Change K and Gamma, Shibata et al. 2004 3-D GR Core Collapse style

 evolve_tracer Scope: restricted BOOLEAN Description: Should we advect tracers? Default: no

 gradient_method Scope: restricted KEYWORD Description: Which method is used to set GRHydro::DiﬀRho? Range Default: First diﬀ First diﬀ Standard ﬁrst diﬀerences Curvature Curvature based method of Paramesh / FLASH Rho weighted Based on the size of rho

 grhydro_atmo_tolerance Scope: restricted REAL Description: A point is set to atmosphere in the Con2Prim’s if its rho < GRHydro_rho_min *(1+GRHydro_atmo_tolerance). This avoids occasional spurious oscillations in carpet buﬀer zones lying in the atmosphere (because prolongation happens on conserved variables) Range Default: 0.0 0.0: Zero or larger. A useful value could be 0.0001

 grhydro_c2p_failed_action Scope: restricted KEYWORD Description: what to do when we detect a c2p failure Range Default: abort abort abort with error terminate request termination

 grhydro_c2p_reset_eps_tau_hot_eos Scope: restricted BOOLEAN Description: As a last resort, reset tau Default: no

 grhydro_c2p_warn_from_reﬂevel Scope: restricted INT Description: Start warning on given reﬁnement level and on higher levels Range Default: (none) 0: Greater or equal to 0

 grhydro_c2p_warnlevel Scope: restricted INT Description: Warnlevel for Con2Prim warnings Range Default: (none) 0:1 Either 0 or 1

 grhydro_countmax Scope: restricted INT Description: Maximum number of iterations for Con2Prim solve Range Default: 100 1:* Must be positive

 grhydro_countmin Scope: restricted INT Description: Minimum number of iterations for Con2Prim solve Range Default: 1 0:* Must be non negative

 grhydro_del_ptol Scope: restricted REAL Description: Tolerance for primitive variable solve (absolute) Range Default: 1.e-18 0: Do we really want both tolerances?

 grhydro_enable_internal_excision Scope: restricted BOOLEAN Description: Set this to ’false’ to disable the thorn-internal excision. Default: true

 grhydro_eos_hot_eps_ﬁx Scope: restricted BOOLEAN Description: Activate quasi-failsafe mode for hot EOSs Default: no

 grhydro_eos_hot_prim2con_warn Scope: restricted BOOLEAN Description: Warn about temperature workaround in prim2con Default: yes

 grhydro_eos_rf_prec Scope: restricted REAL Description: Precision to which root ﬁnding should be carried out Range Default: 1.0e-8 (0.0:* anything larger than 0 goes

 grhydro_eos_table Scope: restricted STRING Description: Name for the Equation of State Range Default: Ideal_Fluid .* Can be anything

 grhydro_eos_type Scope: restricted KEYWORD Description: Type of Equation of State Range Default: General Polytype P = P(rho) or P = P(eps) General P = P(rho, eps)

 grhydro_eps_min Scope: restricted REAL Description: Minimum value of speciﬁc internal energy - this is now only used in GRHydro_InitData’s GRHydro_Only_Atmo routine Range Default: 1.e-10 0: Positive

 grhydro_hot_atmo_temp Scope: restricted REAL Description: Temperature of the hot atmosphere in MeV Range Default: 0.1e0 (0.0:* Larger than 0 MeV

 grhydro_hot_atmo_y_e Scope: restricted REAL Description: Y_e of the hot atmosphere Range Default: 0.5e0 0.0:* Larger than 0

 grhydro_hydro_excision Scope: restricted INT Description: Turns excision automatically on in HydroBase Range Default: 1 1:1 Only ’1’ allowed

 grhydro_lorentz_overshoot_cutoﬀ Scope: restricted REAL Description: Set the Lorentz factor to this value in case it overshoots (1/0) Range Default: 1.e100 0:* Some big value

 grhydro_max_temp Scope: restricted REAL Description: maximum temperature we allow Range Default: 90.0e0 (0.0:* Larger than 0 MeV

 grhydro_maxnumconstrainedvars Scope: restricted INT Description: The maximum number of constrained variables used by GRHydro Range Default: 37 7:48 A small range, depending on testing or not

 grhydro_maxnumevolvedvars Scope: restricted INT Description: The maximum number of evolved variables used by GRHydro Range Default: 5 when using multirate 5:12 dens scon[3] tau (B/A)vec[3] psidc ye entropy Aphi

 grhydro_maxnumevolvedvarsslow Scope: restricted INT Description: The maximum number of evolved variables used by GRHydro Range Default: (none) do not use multirate 5:12 dens scon[3] tau (B/A)vec[3] psidc ye entropy Aphi

 grhydro_maxnumsandrvars Scope: restricted INT Description: The maximum number of save and restore variables used by GRHydro Range Default: 16 0:16 A small range, depending on testing or not

 grhydro_nan_verbose Scope: restricted INT Description: The warning level for NaNs occuring within GRHydro Range Default: 2 0:* The warning level

 grhydro_oppm_reﬂevel Scope: restricted INT Description: Ref level where oPPM is used instead of ePPM (used with use_enhaced_ppm=yes). Range Default: -1 -1:10 0-10 (the reﬂevel number) or -1 (oﬀ)

 grhydro_perc_ptol Scope: restricted REAL Description: Tolerance for primitive variable solve (percent) Range Default: 1.e-8 0: Do we really want both tolerances?

 grhydro_polish Scope: restricted INT Description: Number of extra iterations after root found Range Default: (none) 0:* Must be non negative

 grhydro_rho_central Scope: restricted REAL Description: Central Density for Star Range Default: 1.e-5 :

 grhydro_stencil Scope: restricted INT Description: Width of the stencil Range Default: 2 0: Must be positive

 grhydro_y_e_max Scope: restricted REAL Description: maximum allowed Y_e Range Default: 1.0 0.0:* Larger than or equal to zero; 1 is default

 grhydro_y_e_min Scope: restricted REAL Description: minimum allowed Y_e Range Default: 0.0 0.0:* Larger than or equal to zero

 hlle_type Scope: restricted KEYWORD Description: Which HLLE type to use Range Default: Standard Standard Standard HLLE solver Tadmor Tadmor’s simpliﬁcation of HLLE

 initial_atmosphere_factor Scope: restricted REAL Description: A relative (to the initial atmosphere) value for rho in the atmosphere. This is used at initial time only. Unused if negative. Range Default: -1.0 -1.0:

 initial_gamma Scope: restricted REAL Description: If changing Gamma, what was the value used in the initial data routine? Range Default: 1.3333 (0.0: Positive

 initial_k Scope: restricted REAL Description: If changing K, what was the value used in the initial data routine? Range Default: 100.0 (0.0: Positive

 initial_rho_abs_min Scope: restricted REAL Description: An absolute value for rho in the atmosphere. To be used by initial data routines only. Unused if negative. Range Default: -1.0 -1.0:

 initial_rho_rel_min Scope: restricted REAL Description: A relative (to the central density) value for rho in the atmosphere. To be used by initial data routines only. Unused if negative. Range Default: -1.0 -1.0:

 kap_dc Scope: restricted REAL Description: The kap parameter for divergence cleaning Range Default: 10.0 0:* Any non-negative value, but 1.0 to 10.0 seems preferred

 left_eigenvectors Scope: restricted KEYWORD Description: How to compute the left eigenvectors Range Default: analytical analytical Analytical left eigenvectors numerical Numerical left eigenvectors

 max_magnetic_to_gas_pressure_ratio Scope: restricted REAL Description: consider pressure to be magnetically dominated if magnetic pressure to gas pressure ratio is higher than this Range Default: -1.0 (0:* any positive value, eg. 100. -1.0 disable

 method_type Scope: restricted KEYWORD Description: Which type of method to use Range Default: RSA FV RSA FV ”Reconstruct-Solve-A verage ﬁnite volume method” Flux Split FD Finite diﬀerence using Lax-Friedrichs ﬂux splitting

 min_tracer Scope: restricted REAL Description: The ﬂoor placed on the tracer Range Default: 0.0 *:* Anything

 mp5_adaptive_eps Scope: restricted BOOLEAN Description: Same as WENO adaptive epsilon: adaptively reduce mp5_eps according to norm of stencil. Original algorithm does not use this. Default: no

 mp5_alpha Scope: restricted REAL Description: alpha parameter for MP5 reconstruction Range Default: 4.0 0:* positive

 mp5_eps Scope: restricted REAL Description: epsilon parameter for MP5 reconstruction Range Default: 0.0 0:* 0.0 or very small and positive. 1e-10 is suggested by Suresh&Huynh. TOV star requires 0.0

 myﬂoor Scope: restricted REAL Description: A minimum number for the TVD reconstruction routine Range Default: 1.e-10 0.0: Must be positive

 number_of_arrays Scope: restricted INT Description: Number of arrays to evolve Range Default: (none) 0:3 Either zero or three, depending on the particles

 number_of_particles Scope: restricted INT Description: Number of particles to track Range Default: (none) 0:* 0 switches oﬀ particle tracking

 number_of_tracers Scope: restricted INT Description: Number of tracer variables to be advected Range Default: (none) 0:* positive or zero

 numerical_viscosity Scope: restricted KEYWORD Description: How to compute the numerical viscosity Range Default: fast fast Fast numerical viscosity normal Normal numerical viscosity

 particle_interpolation_order Scope: restricted INT Description: What order should be used for the particle interpolation Range Default: 2 1:* A valid positive interpolation order

 particle_interpolator Scope: restricted STRING Description: What interpolator should be used for the particles Range Default: Lagrange polynomial interpolation .+ A valid interpolator name

 ppm_detect Scope: restricted BOOLEAN Description: Should the PPM solver use shock detection? Default: no

 ppm_epsilon Scope: restricted REAL Description: Epsilon for PPM zone ﬂattening Range Default: 0.33 0.0: Must be positive. Default is from Colella & Woodward

 ppm_epsilon_shock Scope: restricted REAL Description: Epsilon for PPM shock detection Range Default: 0.01 : Anything goes. Default is from Colella & Woodward

 ppm_eta1 Scope: restricted REAL Description: Eta1 for PPM shock detection Range Default: 20.0 : Anything goes. Default is from Colella & Woodward

 ppm_eta2 Scope: restricted REAL Description: Eta2 for PPM shock detection Range Default: 0.05 : Anything goes. Default is from Colella & Woodward

 ppm_ﬂatten Scope: restricted KEYWORD Description: Which ﬂattening procedure should the PPM solver use? Range Default: stencil_3 stencil_3 our ﬂattening procedure, which requires only stencil 3 and which may work stencil_4 original C&W ﬂattening procedure, which requires stencil 4

 ppm_k0 Scope: restricted REAL Description: K0 for PPM shock detection Range Default: 0.2 : Anything goes. Default suggested by Colella & Woodward is: (polytropic constant)/10.0

 ppm_mppm Scope: restricted INT Description: Use modiﬁed (enhanced) ppm scheme Range Default: (none) 0:1 0 (oﬀ, default) or 1 (on)

 ppm_mppm_debug_eigenvalues Scope: restricted INT Description: write eigenvalues into debug grid variables Range Default: (none) 0:1 0 (oﬀ, default) or 1 (on)

 ppm_omega1 Scope: restricted REAL Description: Omega1 for PPM zone ﬂattening Range Default: 0.75 : Anything goes. Default is from Colella & Woodward

 ppm_omega2 Scope: restricted REAL Description: Omega2 for PPM zone ﬂattening Range Default: 10.0 : Anything goes. Default is from Colella & Woodward

 ppm_omega_tracer Scope: restricted REAL Description: Omega for tracer PPM zone ﬂattening Range Default: 0.50 : Anything goes. Default is from Plewa & Mueller

 ppm_small Scope: restricted REAL Description: A ﬂoor used by PPM shock detection Range Default: 1.e-7 0.0:1.0 In [0,1]

 psidcspeed Scope: restricted KEYWORD Description: Which speed to set for psidc Range Default: light speed char speed Based on the characteristic speeds light speed Set the characteristic speeds to speed of light set speed ”Manually set the characteristic speeds: [setcharmin,setcharm ax]”

 recon_method Scope: restricted KEYWORD Description: Which reconstruction method to use Range Default: tvd tvd Slope limited TVD ppm PPM reconstruction eno ENO reconstruction weno WENO reconstruction weno-z WENO-Z reconstruction mp5 MP5 reconstruction

 recon_vars Scope: restricted KEYWORD Description: Which type of variables to reconstruct Range Default: primitive primitive Reconstruct the primitive variables conservative Reconstruct the conserved variables

 reconstruct_temper Scope: restricted BOOLEAN Description: if set to true, temperature will be reconstructed Default: no

 reconstruct_wv Scope: restricted BOOLEAN Description: Reconstruct the primitive velocity W_Lorentz*vel, rather than just vel. Default: no

 rho_abs_min Scope: restricted REAL Description: A minimum rho below which evolution is turned oﬀ (atmosphere). If negative, the relative minimum will be used instead. Range Default: -1.0 -1.0:

 rho_abs_min_after_recovery Scope: restricted REAL Description: A new absolute value for rho in the atmosphere. To be used after recovering. Unused if negative. Range Default: -1.0 -1.0:

 rho_rel_min Scope: restricted REAL Description: A minimum relative rho (rhomin = centden * rho_rel_min) below which evolution is turned oﬀ (atmosphere). Only used if rho_abs_min < 0.0 Range Default: 1.e-9 0:

 riemann_solver Scope: restricted KEYWORD Description: Which Riemann solver to use Range Default: HLLE Roe Standard Roe solver Marquina Marquina ﬂux HLLE HLLE HLLC HLLC LLF Local Lax-Friedrichs (MHD only at the moment)

 set_trivial_rp_grid_function Scope: restricted INT Description: set gf for triv. rp (only for debugging) Range Default: (none) 0:1 0 for no (default), 1 for yes

 setcharmax Scope: restricted REAL Description: Maximum characteristic speed for psidc if psidcspeed is set Range Default: 1.0 0:1 Any value smaller than speed of light

 setcharmin Scope: restricted REAL Description: Minimum characteristic speed for psidc if psidcspeed is set Range Default: -1.0 -1:0 Any value smaller than speed of light - sign should be negative

 sources_spatial_order Scope: restricted INT Description: Order of spatial diﬀerencing of the source terms Range Default: 2 2 2nd order ﬁnite diﬀerencing 4 4th order ﬁnite diﬀerencing

 sqrtdet_thr Scope: restricted REAL Description: Threshold to apply cons rescalings deep inside the horizon Range Default: -1.0 1.0: Larger values guarantees this sort of rescaling only deep inside the horizon -1.0 Do not apply limit

 sync_conserved_only Scope: restricted BOOLEAN Description: Only sync evolved conserved quantities during evolution. Default: no

 tau_rel_min Scope: restricted REAL Description: A minimum relative tau (taumin = maxtau(t=0) * tau_rel_min) below which tau is reschaled Range Default: 1.e-10 0:

 tmunu_damping_radius_max Scope: restricted REAL Description: damping radius at which Tmunu becomes 0 Range Default: -1 -1 damping switched oﬀ 0:* greater than minimum radius above

 tmunu_damping_radius_min Scope: restricted REAL Description: damping radius at which we start to damp with a tanh function Range Default: -1 -1 damping switched oﬀ 0:* damping radius at which we start to damp

 track_divb Scope: restricted BOOLEAN Description: Track the magnetic ﬁeld constraint violations Default: no

 transport_constraints Scope: restricted BOOLEAN Description: Use constraint transport for magnetic ﬁeld Default: no

 tvd_limiter Scope: restricted KEYWORD Description: Which slope limiter to use Range Default: minmod minmod Minmod vanleerMC2 Van Leer MC - Luca Superbee Superbee

 use_enhanced_ppm Scope: restricted BOOLEAN Description: Use the enhanced ppm reconstruction method proposed by Colella & Sekora 2008 and McCorquodale & Colella 2011 Default: no

 use_evolution_mask Scope: restricted KEYWORD Description: Set this to ’always’ to skip validity tests in regions where CarpetEvolutionMask::evolution_mask vanishes. Range Default: never always use the mask auto check if CarpetEvolutionMask is active, then use the mask never do not use the mask

 use_min_tracer Scope: restricted BOOLEAN Description: Should there be a ﬂoor on the tracer? Default: no

 use_mol_slow_multirate_sector Scope: restricted BOOLEAN Description: Whether to make use of MoL’s slow multirate sector Default: no

 use_optimized_ppm Scope: restricted BOOLEAN Description: use C++ templated version of PPM. Experimental Default: no

 use_weighted_ﬂuxes Scope: restricted BOOLEAN Description: Weight the ﬂux terms by the cell surface areas Default: no

 weno_adaptive_epsilon Scope: restricted BOOLEAN Description: use modiﬁed smoothness indicators that take into account scale of solution (adaptive epsilon) Default: yes

 weno_eps Scope: restricted REAL Description: WENO epsilon parameter. For WENO-z, 1e-40 is recommended Range Default: 1e-26 0:* small and positive

 weno_order Scope: restricted INT Description: The order of accuracy of the WENO reconstruction Range Default: 5 5 Fifth-order

 wk_atmosphere Scope: restricted BOOLEAN Description: Use some of Wolfgang Kastauns atmosphere tricks Default: no

11 Interfaces

Implements:

grhydro

Inherits:

boundary

tmunubase

hydrobase

Grid Variables

11.0.1 PRIVATE GROUPS
 Group Names Variable Names Details inlastmolpoststep compact 0 InLastMoLPostStep description Flag to indicate if we are currently in the last MoL_PostStep dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT execute_mol_step compact 0 execute_MoL_Step description Flag indicating whether we use the slow sector of multirate RK time integration dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT execute_mol_poststep compact 0 execute_MoL_PostStep description Flag indicating whether we use the slow sector of multirate RK time integration dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT grhydro_con_bext compact 0 densplus description Conservative variables extended to the cell boundaries sxplus dimensions 3 syplus distribution DEFAULT szplus group type GF tauplus tags Prolongation=”None” checkpoint=”no” densminus timelevels 1 sxminus variable type REAL grhydro_mhd_con_bext compact 0 Bconsxplus description Conservative variables extended to the cell boundaries Bconsyplus dimensions 3 Bconszplus distribution DEFAULT Bconsxminus group type GF Bconsyminus tags Prolongation=”None” checkpoint=”no” Bconszminus timelevels 1 variable type REAL grhydro_mhd_prim_bext compact 0 Bvecxplus description Primitive mhd variables extended to the cell boundaries Bvecyplus dimensions 3 Bveczplus distribution DEFAULT Bvecxminus group type GF Bvecyminus tags Prolongation=”None” checkpoint=”no” Bveczminus timelevels 1 variable type REAL
 Group Names Variable Names Details grhydro_avec_bext compact 0 Avecxplus description Vector potential extended to the cell boundaries Avecyplus dimensions 3 Aveczplus distribution DEFAULT Avecxminus group type GF Avecyminus tags Prolongation=”None” checkpoint=”no” Aveczminus timelevels 1 variable type REAL grhydro_aphi_bext compact 0 Aphiplus description Vector potential phi extended to the cell boundaries Aphiminus dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL grhydro_mhd_psidc_bext compact 0 psidcplus description Divergence cleaning variable extended to the cell boundaries for diverence cleaning psidcminus dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL grhydro_entropy_prim_bext compact 0 entropyplus description Primitive entropy extended to the cell boundaries entropyminus dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL grhydro_entropy_con_bext compact 0 entropyconsplus description Conservative entropy extended to the cell boundaries entropyconsminus dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL whichpsidcspeed compact 0 whichpsidcspeed description Which speed to set for psidc? Set in ParamCheck dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT
 Group Names Variable Names Details grhydro_coords compact 0 GRHydro_x description Coordinates to use with the comoving shift GRHydro_y dimensions 3 GRHydro_z distribution DEFAULT group type GF timelevels 3 variable type REAL grhydro_coords_rhs compact 0 GRHydro_x_rhs description RHS for coordinates to use with the comoving shift GRHydro_y_rhs dimensions 3 GRHydro_z_rhs distribution DEFAULT group type GF tags Prolongation=”None” timelevels 1 variable type REAL grhydro_trivial_rp_gf_group compact 0 GRHydro_trivial_rp_gf_x description set gf for triv. rp (only for debugging) GRHydro_trivial_rp_gf_y dimensions 3 GRHydro_trivial_rp_gf_z distribution DEFAULT group type GF tags Prolongation=”None” timelevels 1 variable type INT ﬂux_splitting compact 0 densfplus description Fluxes for use in the ﬂux splitting densfminus dimensions 3 sxfplus distribution DEFAULT sxfminus group type GF syfplus tags Prolongation=”None” checkpoint=”no” syfminus timelevels 1 szfplus variable type REAL fs_alpha compact 0 fs_alpha1 description Maximum characteristic speeds for the ﬂux splitting fs_alpha2 dimensions 0 fs_alpha3 distribution CONSTANT fs_alpha4 group type SCALAR fs_alpha5 timelevels 1 variable type REAL y_e_plus compact 0 Y_e_plus description Plus state for the electron fraction dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL
 Group Names Variable Names Details y_e_minus compact 0 Y_e_minus description Minus state for the electron fraction dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL tempplus compact 0 tempplus description Plus state for the temperature dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL tempminus compact 0 tempminus description Minus state for the temperature dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL grhydro_tracer_rhs compact 0 cons_tracerrhs description RHS for the tracer dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 vararray_size number_of_tracers variable type REAL grhydro_tracer_ﬂux compact 0 cons_tracerﬂux description Flux for the tracer dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 vararray_size number_of_tracers variable type REAL grhydro_tracer_cons_bext compact 0 cons_tracerplus description Cell boundary values for the tracer cons_tracerminus dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 vararray_size number_of_tracers variable type REAL
 Group Names Variable Names Details grhydro_tracer_prim_bext compact 0 tracerplus description Primitive cell boundary values for the tracer tracerminus dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 vararray_size number_of_tracers variable type REAL grhydro_tracer_ﬂux_splitting compact 0 tracerfplus description Flux splitting for the tracer tracerfminus dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 vararray_size number_of_tracers variable type REAL grhydro_mppm_eigenvalues compact 0 GRHydro_mppm_eigenvalue_x_left description debug variable for ﬂux eigenvalues in mppm GRHydro_mppm_eigenvalue_x_right dimensions 3 GRHydro_mppm_eigenvalue_y_left distribution DEFAULT GRHydro_mppm_eigenvalue_y_right group type GF GRHydro_mppm_eigenvalue_z_left tags Prolongation=”None” checkpoint=”no” GRHydro_mppm_eigenvalue_z_right timelevels 1 GRHydro_mppm_xwind variable type REAL particles compact 0 particle_x description Coordinates of particles to be tracked particle_y dimensions 1 particle_z distribution DEFAULT ghostsize 0 group type ARRAY size NUMBER_OF_PARTICLES timelevels 3 variable type REAL particle_rhs compact 0 particle_x_rhs description RHS functions for particles to be tracked particle_y_rhs dimensions 1 particle_z_rhs distribution DEFAULT ghostsize 0 group type ARRAY size NUMBER_OF_PARTICLES timelevels 1 variable type REAL particle_arrays compact 0 particle_vx description Temporaries to hold interpolated values for particle tracking particle_vy dimensions 1 particle_vz distribution DEFAULT particle_alp ghostsize 0 particle_betax group type ARRAY particle_betay size NUMBER_OF_PARTICLES particle_betaz tags checkpoint=”no” timelevels 1 variable type REAL
 Group Names Variable Names Details grhydro_maxima_location compact 0 maxima_i description The location (point index) of the maximum value of rho maxima_j dimensions 0 maxima_k distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type REAL grhydro_maxima_iteration compact 0 GRHydro_maxima_iteration description Iteration on which maximum was last set dimensions 0 distribution CONSTANT group type SCALAR timelevels 1 variable type INT grhydro_maxima_separation compact 0 GRHydro_separation description The distance between the centres (locations of maximum density) of a binary NS GRHydro_proper_separation dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type REAL diﬀrho compact 0 DiﬀRho description The ﬁrst diﬀerence in rho dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL eos_temps compact 0 eos_cs2_p description Temporaries for the EOS calls eos_cs2_m dimensions 3 eos_dpdeps_p distribution DEFAULT eos_dpdeps_m group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL roeaverage_temps compact 0 rho_ave description Temporaries for the Roe solver velx_ave dimensions 3 vely_ave distribution DEFAULT velz_ave group type GF eps_ave tags Prolongation=”None” checkpoint=”no” press_ave timelevels 1 eos_cs2_ave variable type REAL
 Group Names Variable Names Details con2prim_temps compact 0 press_old description Temporaries for the conservative to primitive conversion press_new dimensions 3 eos_dpdeps_temp distribution DEFAULT eos_dpdrho_temp group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL h_viscosity_temps compact 0 eos_c description Temporaries for H viscosity dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL

11.0.2 PUBLIC GROUPS
 Group Names Variable Names Details grhydro_eos_scalars compact 0 GRHydro_eos_handle description Handle number for EOS GRHydro_polytrope_handle dimensions 0 distribution CONSTANT group type SCALAR timelevels 1 variable type INT grhydro_minima compact 0 GRHydro_rho_min description Atmosphere values GRHydro_tau_min dimensions 0 distribution CONSTANT group type SCALAR timelevels 1 variable type REAL dens compact 0 dens description generalized particle number dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”Scalar” tensorweight=+1.0 jacobian=”inverse_jacobian” interpolator=”matter” timelevels 3 variable type REAL tau compact 0 tau description internal energy dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”Scalar” tensorweight=+1.0 jacobian=”inverse_jacobian” interpolator=”matter” timelevels 3 variable type REAL scon compact 0 scon description generalized momenta dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”D” tensorweight=+1.0 jacobian=”inverse_jacobian” interpolator=”matter” timelevels 3 vararray_size 3 variable type REAL bcons compact 0 Bcons description B-ﬁeld conservative variable dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”U” tensorparity=-1 tensorweight=+1.0 jacobian=”jacobian” interpolator=”matter” timelevels 3 vararray_size 3 variable type REAL
 Group Names Variable Names Details evec compact 0 Evec description Electric ﬁeld at edges dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”U” tensorweight=+1.0 jacobian=”jacobian” interpolator=”matter” timelevels 1 vararray_size 3 variable type REAL y_e_con compact 0 Y_e_con description Conserved electron fraction dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”Scalar” tensorweight=+1.0 jacobian=”inverse_jacobian” interpolator=”matter” timelevels 3 variable type REAL entropycons compact 0 entropycons description Conserved entropy density dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”Scalar” tensorweight=+1.0 jacobian=”inverse_jacobian” interpolator=”matter” timelevels 3 variable type REAL grhydro_tracers compact 0 tracer description Tracers dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”Scalar” timelevels 3 vararray_size number_of_tracers variable type REAL sdetg compact 0 sdetg description Sqrt of the determinant of the 3-metric dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” tensortypealias=”Scalar” tensorweight=+1.0 interpolator=”matter” checkpoint=”no” timelevels 1 variable type REAL psidc compact 0 psidc description Psi parameter for divergence cleaning dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”Scalar” tensorweight=+1.0 tensorparity=-1 jacobian=”inverse_jacobian” interpolator=”matter” timelevels 3 variable type REAL
 Group Names Variable Names Details densrhs compact 0 densrhs description Update term for dens dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL taurhs compact 0 taurhs description Update term for tau dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL srhs compact 0 srhs description Update term for s dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 vararray_size 3 variable type REAL bconsrhs compact 0 Bconsrhs description Update term for Bcons dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 vararray_size 3 variable type REAL avecrhs compact 0 Avecrhs description Update term for Avec dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 vararray_size 3 variable type REAL aphirhs compact 0 Aphirhs description Update term for Aphi dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL
 Group Names Variable Names Details psidcrhs compact 0 psidcrhs description Update term for psidc dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL entropyrhs compact 0 entropyrhs description Update term for entropycons dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL divb compact 0 divB description Magnetic ﬁeld constraint dimensions 3 distribution DEFAULT group type GF tags Prolongation=”Restrict” checkpoint=”no” tensorparity=-1 timelevels 1 variable type REAL bcom compact 0 bcom description bî: comoving contravariant magnetic ﬁeld 4-vector spatial components dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”U” tensorparity=-1 interpolator=”matter” timelevels 3 vararray_size 3 variable type REAL bcom0 compact 0 bcom0 description b0̂ component of the comoving contravariant magnetic ﬁeld 4-vector dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”Scalar” interpolator=”matter” timelevels 3 variable type REAL bcom_sq compact 0 bcom_sq description half of magnectic pressure: contraction of b_a bâ dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”Scalar” interpolator=”matter” timelevels 3 variable type REAL
 Group Names Variable Names Details lvel compact 0 lvel description local velocity vî dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”U” jacobian=”jacobian” interpolator=”matter” timelevels 3 vararray_size 3 variable type REAL lbvec compact 0 lBvec description local Magnetic ﬁeld components Bî dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”U” jacobian=”jacobian” tensorparity=-1 interpolator=”matter” timelevels 3 vararray_size 3 variable type REAL local_metric compact 0 gaa description local ADM metric g_ij gab dimensions 3 gac distribution DEFAULT gbb group type GF gbc tags Prolongation=”None” checkpoint=”no” gcc timelevels 3 variable type REAL local_extrinsic_curvature compact 0 kaa description local extrinsic curvature K_ij kab dimensions 3 kac distribution DEFAULT kbb group type GF kbc tags Prolongation=”None” checkpoint=”no” kcc timelevels 1 variable type REAL local_shift compact 0 betaa description local ADM shift ∖betaî betab dimensions 3 betac distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL grhydro_prim_bext compact 0 rhoplus description Primitive variables extended to the cell boundaries velxplus dimensions 3 velyplus distribution DEFAULT velzplus group type GF pressplus tags Prolongation=”None” checkpoint=”no” epsplus timelevels 1 w_lorentzplus variable type REAL
 Group Names Variable Names Details grhydro_scalars compact 0 ﬂux_direction description Which direction are we taking the ﬂuxes in and the oﬀsets xoﬀset dimensions 0 yoﬀset distribution CONSTANT zoﬀset group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT grhydro_atmosphere_mask compact 0 atmosphere_mask description Flags to say whether a point needs to be reset to the atmosphere dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” timelevels 1 variable type INT grhydro_atmosphere_mask_real compact 0 atmosphere_mask_real description Flags to say whether a point needs to be reset to the atmosphere. This is sync’ed (and possibly interpolated)! dimensions 3 distribution DEFAULT group type GF tags Prolongation=”sync” checkpoint=”no” timelevels 1 variable type REAL grhydro_atmosphere_descriptors compact 0 atmosphere_ﬁeld_descriptor dimensions 0 atmosphere_atmosp_descriptor distribution CONSTANT atmosphere_normal_descriptor group type SCALAR timelevels 1 variable type INT grhydro_cons_tracers compact 0 cons_tracer description The conserved tracer variable dimensions 3 distribution DEFAULT group type GF tags ProlongationParameter=”HydroBase::prolongation_type” tensortypealias=”Scalar” timelevels 3 vararray_size number_of_tracers variable type REAL grhydro_maxima_position compact 0 maxima_x description The position (coordinate values) of the maximum value of rho maxima_y dimensions 0 maxima_z distribution CONSTANT maximum_density group type SCALAR tags checkpoint=”no” timelevels 1 variable type REAL
 Group Names Variable Names Details maxrho_global compact 0 maxrho_global description store the global maximum of rho description for reﬁnment-grid steering dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type REAL grhydro_c2p_failed compact 0 GRHydro_C2P_failed description Mask that stores the points where C2P has failed dimensions 3 distribution DEFAULT group type GF tags Prolongation=”restrict” tensortypealias=”Scalar” checkpoint=”no” timelevels 1 variable type REAL grhydro_ﬂuxes compact 0 densﬂux description Fluxes for each conserved variable sxﬂux dimensions 3 syﬂux distribution DEFAULT szﬂux group type GF tauﬂux tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL grhydro_bﬂuxes compact 0 Bconsxﬂux description Fluxes for each B-ﬁeld variable Bconsyﬂux dimensions 3 Bconszﬂux distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL grhydro_psiﬂuxes compact 0 psidcﬂux description Fluxes for the divergence cleaning parameter dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL grhydro_entropyﬂuxes compact 0 entropyﬂux description Fluxes for the conserved entropy density dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL
 Group Names Variable Names Details grhydro_avecﬂuxes compact 0 Avecxﬂux description Fluxes for each Avec variable Avecyﬂux dimensions 3 Aveczﬂux distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL grhydro_aphiﬂuxes compact 0 Aphiﬂux description Fluxes for Aphi dimensions 3 distribution DEFAULT group type GF tags Prolongation=”None” checkpoint=”no” timelevels 1 variable type REAL evolve_y_e compact 0 evolve_Y_e description Are we evolving Y_e? Set in Paramcheck dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT evolve_temper compact 0 evolve_temper description Are we evolving temperature? Set in Paramcheck dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT evolve_entropy compact 0 evolve_entropy description Are we evolving entropy? Set in Paramcheck dimensions 0 distribution CONSTANT group type SCALAR tags checkpoint=”no” timelevels 1 variable type INT evolve_mhd