The GRHydro code: three-dimensional relativistic hydrodynamics

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

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

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

http://www.einsteintoolkit.org

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

GRHydro/test/GRHydro_test_shock.par

and will include the essential thorns

GRHydro EOS_Omni ADMBase ADMCoupling MoL

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 efficient) method can be chosen by setting ode_method to "genrk" which is a Shu-Osher type TVD Runge-Kutta evolution. Then the parameter mol_timesteps controls the number of intermediate steps and hence the order of accuracy. First to seventh order are currently supported.

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

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

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

2.1 Obtaining This Thorn

The public version of GRHydro can be found on the website http://www.einsteintoolkit.org.

2.2 Basic Usage

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

2.3 Special Behaviour

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

2.4 Interaction With Other Thorns

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

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

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

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

2.5 Support and Feedback

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

http://www.einsteintoolkit.org

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

3 Physical System

The equations of general relativistic hydrodynamics can be written in the flux conservative form \begin {equation} \label {eq:consform1} \partial _t {\bf q} + \partial _{x^i} {\bf f}^{(i)} ({\bf q}) = {\bf s} ({\bf q}), \end {equation} where \(\bf q\) is a set of conserved variables, \({\bf f}^{(i)} ({\bf q})\) the fluxes and \({\bf s} ({\bf q})\) the source terms. The 8 conserved variables are labeled \(D\), \(S_i\), \(\tau \), and \(\mathcal B^k\), where \(D\) is the generalized particle number density, \(S_i\) are the generalized momenta in each direction, \(\tau \) is an internal energy term, and \(\mathcal B^k\) is the densitized magnetic field. These conserved variables are composed from a set of primitive variables, which are \(\rho \), the rest-mass density, \(p\), the pressure, \(v^i\), the fluid 3-velocities, \(\epsilon \), the specific internal energy, \(B^k\) the magnetic field in the lab frame, and \(W\), the Lorentz factor, via the following relations  

\begin {eqnarray} \label {eq:prim2con} D &=& \sqrt {\gamma }W\rho \nonumber \\ S_i &=& \sqrt {\gamma } \left (\rho h^* W^{\,2} v_j-\alpha b^0b_j\right ) \nonumber \\ \tau &=& \sqrt {\gamma }\left (\rho h^* W^2 - p^*-(\alpha b^0)^2\right ) - D, \nonumber \\ \mathcal B^k &=& \sqrt {\gamma } B^k, \end {eqnarray}

where \(\gamma \) is the determinant of the spatial 3-metric \(\gamma _{ij}\) and \(h^* \equiv 1+\epsilon +\left (p + b^2\right )/\rho \), \(p^* = p + b^2/2\). \(b^\mu \) is the magnetic field in the fluid’s rest frame \(b^\mu = u_\nu ^*F^{\mu \nu }\) where \(^*F^{\mu \nu }\) is the dual of the Faraday tensor. It is related to \(B^k\) via

\begin {eqnarray} b^0&=&\frac {WB^kv_k}{\alpha }\,\,,\\ b^i&=&\frac {B^i}{W}+W(B^kv_k)\left (v^i-\frac {\beta ^i}{\alpha }\right )\,\,,\\ b^2&=&\frac {B^iB_i}{W^2}+(B^iv_i)^2\,\,. \end {eqnarray}

Only five of the primitive fluid variables are independent. The Lorentz factor is defined in terms of the velocities and the metric as \(W = (1-\gamma _{ij}v^i v^j)^{-1/2}\). Also one of the pressure, rest-mass density or specific internal energy terms is given in terms of the other two by an equation of state.

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

\begin {eqnarray} {\bf F}^i({\bf U}) &=& [D(\alpha v^i - \beta ^i), S_j(\alpha v^i - \beta ^i) + p\delta ^i_j, \tau (\alpha v^i - \beta ^i) + p v^i, \mathcal B^k (\alpha v^i - \beta ^i) - \mathcal B^i (\alpha v^k - \beta ^k)]^T\ . \end {eqnarray}

The source terms are

\begin {eqnarray} \label {source_terms} {\bf s}({\bf U}) = \Big [0, T^{\mu \nu }\big (\partial _\mu g_{\nu j} + \Gamma ^\delta _{\mu \nu } g_{\delta j}\big ), \alpha \big (T^{\mu 0}\partial _\mu \ln \alpha - T^{\mu \nu }\Gamma ^0_{\nu \mu } \big ), 0 \Big ]\ . \end {eqnarray}

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

The stress energy tensor is in turn given by

\begin {eqnarray} T^{\mu \nu } &=& \left ( \rho + \rho \epsilon + p + b^2 \right ) u^\mu u^\nu + \left (p + \frac {b^2}{2} \right ) g^{\mu \nu } - b^\mu b^\nu \\\nonumber &\equiv & \rho h^*u^\mu u^\nu + P^* g^{\mu \nu } - b^\mu b^\nu . \label {mhd-stress-energy-tensor} \end {eqnarray}

For more detail see the review of Font [9] and the GRHydro paper [21].

4 Numerical Implementation

TODO: describe MHD scheme in particular constrained transport and con2prim method used.

5 High Resolution Shock Capturing methods

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

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

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

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

\begin {eqnarray} \label {eq:molrhs1} \nonumber {\bf L}({\bf q}) & = & \int \!\!\!\! \int \!\!\!\! \int {\bf s} \,{\rm d}^3 x + \int _{x^2_{j-1/2}}^{x^2_{j+1/2}} \int _{x^3_{k-1/2}}^{x^3_{k+1/2}} {\bf f}^{(1)} ({\bf q} (x^1_{i-1/2}, y, z)) {\rm d} y \, {\rm d} z - \\ && \int _{x^2_{j-1/2}}^{x^2_{j+1/2}} \int _{x^3_{k-1/2}}^{x^3_{k+1/2}} {\bf f}^{(1)} ({\bf q} (x^1_{i+1/2}, y, z)) {\rm d} y \, {\rm d} z + \cdots \end {eqnarray}

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

In this implementation of GRHydro the right hand side operator \(\bf L\) will be simplified considerably by approximating the integrals by the midpoint rule to get \begin {equation} \label {eq:molrhs2} {\bf L}({\bf q}) = {\bf s}_{i,j,k} + {\bf f}^{(1)}_{i-1/2,j,k} - {\bf f}^{(1)}_{i+1/2,j,k} + \cdots \end {equation} where the dependence of the flux on \(\bf q\) and spatial position is implicit in the notation. Given this simplification, the calculation of the right hand side operator splits simply into the following two parts:

  1. Calculate the source terms \({\bf s}({\bf q}(x^1_i, x^2_j, x^3_k))\):

    Given that \(\bar {{\bf q}}\) is a second-order accurate approximation to \(\bf q\) at the midpoint of the cell, which is precisely \((x^1_i, x^2_j, x^3_k)\), for second-order accuracy it is sufficient to use \({\bf s}(\bar {{\bf q}}_{i,j,k})\).

  2. In each coordinate direction \(x^a\), calculate the intercell flux \({\bf f}^{(a)}({\bf q}_{i+1/2,j,k})\):

    From the initial data \(\bar {{\bf q}}\) given at time \(t^n\) we can reconstruct the data at the cell boundary, \(({\bf q}_{i+1/2,j,k})\) to any required accuracy in space (except in the vicinity of a shock, where only first-order accuracy is guaranteed). However this will only be zeroth-order accurate in time. To ensure first-order accuracy in time, we have to find \(({\bf q}_{i+1/2,j,k})(t)\) while retaining the high spatial order of accuracy. This requires two steps:

    1. Reconstruct the data \(\bf 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 \(({\bf q}_{i+1/2,j,k})\), which we call \({\bf q}_L\) and \({\bf q}_R\), where \({\bf q}_L\) is obtained from cell \(i\) (left cell) and \({\bf q}_R\) from cell \(i+1\) (right cell).

    2. The values \({\bf q}_{L,R}\) are used as initial data for the Riemann problem. This is the initial value problem given by the partial differential equation (??) with semi-infinite piecewise constant initial data \({\bf q}_{L,R}\). As the true function \(\bf q\) is probably not piecewise constant we will not get the exact solution of the general problem even if we solve the local Riemann problem exactly. However, it will be first-order accurate in time and retain the high order of accuracy in space which, as explained in the documentation to the MoL thorn, is sufficient to ensure that the method as a whole has a high order of accuracy.

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

6 Reconstruction

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

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

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

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

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

Switching our attention to conservation, we note that there is precisely one conservative first-order reconstruction method, \begin {equation} \label {eq:reconfirst} {\bf q}^{{\rm First}}(x) = \bar {{\bf q}}_i, \qquad x \in [ x_{i-1/2}, x_{i+1/2} ], \end {equation} and that any second-order conservative method can be written in terms of a slope or rather difference \(\Delta _i\) as \begin {equation} \label {eq:reconsecond} {\bf q}^{{\rm Second}}(x) = \bar {{\bf q}}_i + \frac {x - x_i}{x_{i+1/2} - x_{i-1/2}} \Delta _i, \qquad x \in [ x_{i-1/2}, x_{i+1/2} ]. \end {equation}

6.1 TVD Reconstruction

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

In the articles describing the GRAstro_Hydro code4 , this was described as an average of the two reconstructions, \begin {equation} {\bf q}^{{\rm TVD}}(x) = \phi ({\bf q}) {\bf q}^{{\rm Second}} + (1 - \phi ({\bf q})) {\bf q}^{{\rm First}}, \label {First_qTVD} \end {equation} where \(\phi \in [0,1]\) varies smoothly in some sense, and is zero near a shock and 1 in smooth regions. In Toro’s notation [20] (which we usually adopt here) the slope limiter function \(\phi \) (having the same attributes as above) directly multiplies the slope, giving \begin {equation} {\bf q}^{{\rm TVD}}(x) = \bar {{\bf q}}_i + \frac {x - x_i}{x_{i+1/2} - x_{i-1/2}} \phi ({\bf q}) \Delta _i, \qquad x \in [ x_{i-1/2}, x_{i+1/2} ]. \label {Toro_qTVD} \end {equation} Equations (??) and (??) are equivalent.

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

As an example, we show how TVD with the minmod limiter is implemented in the code. First, we define the minmod function: \begin {equation} \label {eq:tvdminmod} \mathrm {\mathbf {minmod}}(a,b) = \left \{ \begin {array}{c l} \text {min}(|a|,|b|) & \text {if } (a b > 0)\\ 0 & \text {otherwise}. \end {array}\right . \end {equation} For reconstructing \(\mathbf {q}\) we choose two differences \begin {equation} \begin {array}{lcl} \Delta _{\mathrm {upw}} & = & {\mathbf {q}}_{i} - {\mathbf {q}}_{i-1}\\ \Delta _{\mathrm {loc}} & = & {\mathbf {q}}_{i+1} - {\mathbf {q}}_{i}\,,\\ \end {array} \end {equation} and write \begin {equation} {\bf q}^{{\rm TVD,minmod}}(x) = \bar {{\bf q}}_i + \frac {x - x_i}{x_{i+1/2} - x_{i-1/2}} \mathbf {minmod}(\Delta _{\mathrm {upw}},\Delta _{\mathrm {loc}}), \qquad x \in [ x_{i-1/2}, x_{i+1/2} ]. \end {equation}

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 first step is to interpolate a quadratic polynomial to the cell boundary, \begin {equation} \label {eq:ppm1} 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 ), \end {equation} where \begin {equation} \label {eq:ppmdm1} \delta _m q_i = \left \{ \begin {array}{c l} \text {min}(|\delta q_i|, 2|q_{i+1} - q_i|, 2|q_i - q_{i-1}|) \text { sign}(\delta q_i) & \text {if } (q_{i+1} - q_i)(q_i - q_{i-1}) > 0, \\ 0 & \text {otherwise}. \end {array} \right ., \end {equation} and \begin {equation} \label {eq:ppmd1} \delta q_i = \frac {1}{2}(q_{i+1} - q_{i-1}). \end {equation} At this point we set both left and right states at the interface to be equal to this, \begin {equation} \label {eq:ppmset1} q_i^R = q_{i+1}^L = q_{1+1/2}. \end {equation}

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

\begin {eqnarray} \label {eq:ppmmonot} q_i^L = q_i^R = q_i & \text {if} & (q_i^R - q_i)(q_i - q_i^L) \leq 0 \\ q_i^L = 3 q_i - 2q_i^R & \text {if} & (q_i^R - q_i^L)\left ( q_i - \frac {1}{2} (q_i^L + q_i^R) \right ) > \frac {1}{6}(q_i^R - q_i^L)^2 \\ q_i^R = 3 q_i - 2q_i^L & \text {if} & (q_i^R - q_i^L)\left ( q_i - \frac {1}{2} (q_i^L + q_i^R) \right ) < -\frac {1}{6}(q_i^R - q_i^L)^2. \end {eqnarray}

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

\begin {eqnarray} \label {eq:ppmdetect} \rho _i^L & = & \rho _i^L (1-\eta ) + \left (\rho _{i-1} + \frac {1}{2} \delta _m \rho _{i-1} \right ) \eta \\ \rho _i^R & = & \rho _i^R (1-\eta ) + \left (\rho _{i+1} - \frac {1}{2} \delta _m \rho _{i+1} \right ) \eta \end {eqnarray}

where \(\eta \) is only applied if the discontinuity is mostly a contact (see [4] for the details) and is defined as \begin {equation} \label {eq:ppmeta} \eta = \text {max}(0, \text {min}(1, \eta _1 (\tilde {\eta } - \eta _2))), \end {equation} where \(\eta _1,\eta _2\) are positive constants and \begin {equation} \label {eq:ppmetatilde} \tilde {\eta } = \left \{ \begin {array}{c l} \frac {\rho _{i-2} - \rho _{i+2} + 4 \delta \rho _i}{12\delta \rho _i} & \text { if } \left \{ \begin {array}{l} \delta ^2\rho _{i+1}\delta ^2\rho _{i-1} < 0\\ (\rho _{i+1} - \rho _{i-1}) - \epsilon \text {min}(|\rho _{i+1}|,|\rho _{i-1}|) > 0\nonumber \end {array} \right . \\ 0 & \text {otherwise} \end {array} \right ., \end {equation} with \(\epsilon \) another positive constant and \begin {equation} \label {eq:ppmd2rho} \delta ^2\rho _i = \frac {\rho _{i+1} - 2\rho _i + \rho _{i-1}}{6\Delta x^2}. \end {equation}

Another step that is performed before monotonicity enforcement is to flatten the zone structure near shocks. This adds simple dissipation and is always in the code. In short, the reconstructions are again altered to \begin {equation} \label {eq:ppmflatten} q_i^{L,R} = \nu _i q_i^{L,R} + (1 - \nu _i) q_i, \end {equation} where \begin {equation} \label {eq:ppmflatten2} \nu _i = \left \{ \begin {array}{c l} {\rm max}(0, 1 - \text {max}(0, \omega _2 (\frac {p_{i+1} - p_{i-1}}{p_{i+2} - p_{i-2}} - \omega _1))) & \text { if } \omega _0 \text {min}(p_{i-1}, p_{i+1}) < |p_{i+1} - p_{i - 1}| \text { and } v^x_{i-1} - v^x_{i+1} > 0 \\ 1 & \text {otherwise} \end {array} \right . \end {equation} and \(\omega _0, \omega _1,\omega _2\) are constants.

The above flattening procedure is not the one in the original article of Colella and Woodward, but it has been adapted from it in order to have a stencil of three points. The original flattening procedure is also implemented in GRHydro. Instead of ??, it consists in the formula \begin {equation} \label {eq:ppmflatten-stencil4} q_i^{L,R} = \tilde \nu _i q_i^{L,R} + (1 - \tilde \nu _i) q_i, \end {equation} where

\begin {eqnarray} \tilde \nu _i &=& {\rm max}\Big (\nu _i,\nu _{i+{\rm sign}(p_{i-1}-p_{i+1})}\Big ) \end {eqnarray}

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

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 differences to measure the smoothness.

Let \(k\) be the order of the reconstruction. Suppose we are reconstructing the scalar function \(q\) in cell \(i\). We start with the cell \(I_i\). We then add to the stencil cell \(I_j\), where \(j = i \pm 1\), where we choose \(j\) to minimize the Newton divided differences

\begin {eqnarray} \label {enodd} q \left [ x_{i-1}, x_i \right ] & = & \frac {q_i - q_{i-1}}{x_{i+1/2} - x_{i-3/2}} \\ q \left [ x_i, x_{i+1} \right ] & = & \frac {q_{i+1} - q_i}{x_{i+3/2} - x_{i-1/2}}. \end {eqnarray}

We then recursively add more cells, minimizing the higher-order Newton divided differences \(q \left [ x_{i-k}, \dots , x_{i+j} \right ]\) defined by \begin {equation} \label {enodd2} 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}}. \end {equation} The reconstruction at the cell boundary is given by a standard \(k^{\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 \begin {equation} \label {enostencil1} S(i) = \left \{ I_{i-r}, \dots , I_{i+k-r-1} \right \}, \end {equation} 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 \begin {equation} \label {enoc1} q_{i+1/2} = \sum _{j=0}^{k-1} c_{rj} q_{i-r+j}, \qquad q_{i-1/2} = \sum _{j=0}^{k-1} c_{r-1,j} q_{i-r+j}. \end {equation} The constants \(c_{rj}\) are given by the rather complicated formula \begin {equation} \label {enoc2} c_{rj} = \left \{ \sum _{m=j+1}^k \frac { \sum _{l=0, l \neq m}^k \prod _{q=0, q \neq m, l}^k \left ( x_{i+1/2} - x_{i-r+q-1/2} \right ) }{ \prod _{l=0, l \neq m}^k \left ( x_{i-r+m-1/2} - x_{i-r+L-1/2} \right ) } \right \} \Delta x_{i-r+j}. \end {equation} This simplifies considerably if the grid is even. The coefficients for an even grid are given (up to seventh order) by [19].

7 Riemann Problems

Given the reconstructed data, we then solve a local Riemann problem in order to get the intercell flux. The Riemann problem is specified by an equation in flux conservative homogeneous form, \begin {equation} \label {eq:homconsform1} \partial _t {\bf q} + \partial _{x^i} {\bf f}^{(i)} ({\bf q}) = 0 \end {equation} with piecewise constant initial data \({\bf q}_{_{L,R}}\) separated by a discontinuity at \(x^{(1)}=0\). Flux terms for the other directions are given similarly. There is no intrinsic scale to this problem and so the solution must be a self similar solution with similarity variable \(\xi = x^{(1)}/t\). The solution is given in terms of waves which separate different states, with each state being constant. The waves are either shocks, across which all hydrodynamical variables change discontinuously, rarefactions, across which all the variables change continuously (the wave is not a single value of \(\xi \) for a rarefaction, but spreads across a finite range), or contact or tangential discontinuities, across which some but not all of the hydrodynamical variables change discontinuously and the rest are constant. The characteristics of the matter evolution converge and break at a shock, diverge at a rarefaction and are parallel at the other linear discontinuities.

The best references for solving the Riemann problem either exactly or approximately are [15], [20], [14]. Here, we start by giving a simple outline. We start by considering the \(N\) dimensional linear problem in one dimension given by \begin {equation} \label {lsrp} \partial _t {\bf q} + A \partial _{x} {\bf q} = 0 \ , \end {equation} where \(A\) is a \(N\times N\) matrix with constant coefficients. We define the eigenvalues \(\lambda ^j\) with associated right and left eigenvectors \({\bf r}^j,{\bf l}_j\), where the eigenvectors are normalized to each other (i.e., their dot product is \(\delta ^i_j\)). We shall assume that the eigenvectors span the space. The characteristic variables \({\bf w}_i\) are defined by \begin {equation} \label {charvar} {\bf w}_i = {\bf l}_i \cdot {\bf q}. \end {equation} Then equation (??) when written in terms of the characteristic variables becomes \begin {equation} \label {charvarrp} \partial _t {\bf w} + \Lambda \partial _x {\bf w} = 0, \end {equation} where \(\Lambda \) is the matrix containing the eigenvalues \(\lambda _i\) on the diagonals and zeros elsewhere. Hence each characteristic variable \({\bf 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 \(\bf q\) we order the variables in such a way that \(\lambda _1 \leq \dots \leq \lambda _N\). We also define the differences in the characteristic variables \(\Delta {\bf w}_i = ({\bf w}_i)_L - ({\bf w}_i)_R\) across the \(i^{{\rm th}}\) characteristic wave. These differences are single numbers (‘scalars’). We note that these differences can easily be found from the initial data using \begin {equation} \label {dw} \Delta {\bf w}_i = {\bf l}_i \cdot \left ( {\bf q}_L - {\bf q}_R \right ). \end {equation} As the change in the solution across each wave is precisely the difference in the associated characteristic variable, the solution of the Riemann problem in terms of characteristic variables can be written as either \begin {equation} \label {lsrpsol1} {\bf w}_i = ({\bf w}_i)_L + \sum _{j=1}^M \Delta {\bf w}_j {\bf e}^j \quad {\rm if}\ \lambda _M < \xi < \lambda _{M+1}, \end {equation} or \begin {equation} \label {lsrpsol2} {\bf w} = ({\bf w}_i)_R - \sum _{j=M+1}^N \Delta {\bf w}_j {\bf e}^j \quad {\rm if}\ \lambda _M < \xi < \lambda _{M+1}, \end {equation} or as the average \begin {equation} \label {lsrpsol3} {\bf w}_i = \frac {1}{2} \left ( ({\bf w}_i)_L + ({\bf w}_i)_R + \sum _{j=1}^M \Delta {\bf w}_j {\bf e}^j - \sum _{j=M+1}^N \Delta {\bf w}_j {\bf e}^j \right ) \quad {\rm if}\ \lambda _M < \xi < \lambda _{M+1}, \end {equation} where \({\bf e}^i\) is the column vector \(({\bf e}^i)_j = \delta ^i_j\).

Converting back to the original variables \(\bf q\) we have the solution \begin {equation} \label {lsrpsol4} {\bf q} = \frac {1}{2} \left ( {\bf q}_L + {\bf q}_R + \sum _{i=1}^M \Delta {\bf w}_i {\bf r}^i - \sum _{i=M+1}^N \Delta {\bf w}_i {\bf r}^i \right ) \quad {\rm if}\ \lambda _M < \xi < \lambda _{M+1}. \end {equation} In the case where we are only interested in the flux along the characteristic \(\xi = 0\) we can write the solution in the simple form \begin {equation} \label {lsrpsol6} {\bf f}({\bf q}) = \frac {1}{2} \left ( {\bf f}({\bf q}_L) + {\bf f}({\bf q}_R) - \sum _{i=1}^N | \lambda _i | \Delta {\bf w}_i {\bf r}^i \right ). \end {equation}

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

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

7.1 HLLE solver

The Harten-Lax-van Leer-Einfeldt (HLLE) solver of Einfeldt [7] is a simple two-wave approximation. We assume that the maximum and minimum wave speeds \(\xi _{\pm }\) are known. The solution of the Riemann problem is then given by requiring conservation to hold across the waves. The solution takes the form \begin {equation} \label {hlle1} {\bf q} = \left \{ \begin {array}[c]{r c l} {\bf q}_L & {\rm if} & \xi < \xi _- \\ {\bf q}_* & {\rm if} & \xi _- < \xi < \xi _+ \\ {\bf q}_R & {\rm if} & \xi > \xi _+, \end {array}\right . \end {equation} and the intermediate state \({\bf q}_*\) is given by \begin {equation} \label {hlle2} {\bf q}_* = \frac {\xi _+ {\bf q}_R - \xi _- {\bf q}_L - {\bf f}({\bf q}_R) + {\bf f}({\bf q}_L)}{\xi _+ - \xi _-}. \end {equation} If we just want the numerical flux along the boundary then this takes the form \begin {equation} \label {hlleflux} {\bf f}({\bf q}) = \frac {\widehat {\xi }_+{\bf f}({\bf q}_L) - \widehat {\xi }_-{\bf f}({\bf q}_R) + \widehat {\xi }_+ \widehat {\xi }_- ({\bf q}_R - {\bf q}_L)}{\widehat {\xi }_+ - \widehat {\xi }_-}, \end {equation} where \begin {equation} \label {hlle3} \widehat {\xi }_- = {\rm min}(0, \xi _-), \quad \widehat {\xi }_+ = {\rm max}(0, \xi _+). \end {equation}

Knowledge of the precise minimum and maximum characteristic velocities \(\xi _{\pm }\) requires knowing the solution of the Riemann problem. Instead, the characteristic velocities are usually found from the eigenvalues of the Jacobian matrix \(\partial {\bf f} / \partial {\bf 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 \([\xi _-, \xi _+]\).

If we set \(\alpha = {\rm max}(|\xi _-|, |\xi _+|)\) and replace the characteristic velocities \(\xi _{\pm }\) with \(\pm \alpha \), we find the Lax–Friedrichs flux (cf. also Tadmor’s semi-discrete scheme [13]) \begin {equation} \label {lfflux} {\bf f}({\bf q}) = \frac {1}{2} \left [ {\bf f}({\bf q}_L) + {\bf f}({\bf q}_R) + \alpha ({\bf q}_L -{\bf q}_R) \right ]. \end {equation} This is very diffusive, but also very stable.

7.2 Roe solver

The linearized solver of Roe [18] is probably the most popular approximate Riemann solver. The simplest interpretation is that the Jacobian \(\partial {\bf f} / \partial {\bf q}\) is linearized about some intermediate state. Then the conservation form reduces to the linear equation \begin {equation} \label {roe1} \partial _t {\bf q} + A \partial _x {\bf q} = 0, \end {equation} where \(A\) is a constant coefficient matrix. This is identical to equation (??) and so all the results of section 7 on linear systems hold. We reiterate that the standard form for the flux along the characteristic ray \(\xi =0\) is \begin {equation} \label {roe2} {\bf f}({\bf q}) = \frac {1}{2} \left ( {\bf f}({\bf q}_L) + {\bf f}({\bf q}_R) - \sum _{i=1}^N | \lambda _i | \Delta {\bf w}_i {\bf r}^i \right ). \end {equation}

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

  1. \(A({\bf q}_{{\rm Roe}}) \left ( {\bf q}_R - {\bf q}_L \right ) = {\bf f}({\bf q}_R) - {\bf f}({\bf q}_L)\),

  2. \(A({\bf q}_{{\rm Roe}})\) is diagonalizable with real eigenvalues,

  3. \(A({\bf q}_{{\rm Roe}}) \rightarrow \partial {\bf f} / \partial {\bf q}\) smoothly as \({\bf q}_{{\rm Roe}} \rightarrow {\bf q}\).

A true Roe average for relativistic hydrodynamics, i.e., an intermediate state that satisfies all these conditions, has been constructed by Eulderink [8]. However, frequently it is sufficient to use \begin {equation} \label {roe3} {\bf q}_{{\rm Roe}} = \frac {1}{2} \left ( {\bf q}_R + {\bf q}_L \right ), \end {equation} which satisfies only the last two conditions. For simplicity we have implemented this arithmetic average.

7.3 Marquina solver

Unlike all the other Riemann solvers introduced so far, the Marquina solver as outlined in [6] does not solve the Riemann problem completely. Instead, only the flux along the characteristic ray \(\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 fluid velocity is equal to the speed of sound of the fluid. In the context of Riemann problems, sonic points are found when the ray \(\xi =0\) is within a rarefaction wave.

Firstly define the left \({\bf l}({\bf q}_{L,R})\) and right \({\bf r}({\bf q}_{L,R})\) eigenvectors and the eigenvalues \(\lambda ({\bf q}_{L,R})\) of the Jacobian matrix \(\partial {\bf f} / \partial {\bf q}\) evaluated at the left and right states. Next define left and right characteristic variables \({\bf w}_{L,R}\) and fluxes \({\bf \phi }_{L,R}\) by \begin {equation} \label {marq1} ({\bf w}_i)_{L,R} = {\bf l}_i({\bf q}_{L,R}) \cdot {\bf q}_{L,R}, \quad ({\bf \phi }_i)_{L,R} = {\bf l}_i({\bf q}_{L,R}) \cdot {\bf f}({\bf q}_{L,R}). \end {equation}

Then the algorithm chooses the correct-sided characteristic flux if the eigenvalues \(\lambda _i({\bf q}_L)\), \(\lambda _i({\bf q}_R)\) have the same sign, and uses a Lax–Friedrichs type flux if they change sign. In full, the algorithm is given in figure 1.


\begin {equation} \label {marqalg} \begin {array}[l]{l} {\bf For}\ \, i = 1, \dots , N\ {\bf do} \\ \qquad \begin {array}[c]{l} {\bf If}\ \, \lambda _i({\bf q}_L) \lambda _i({\bf q}_R) > 0 \ \, {\bf then} \\ \qquad \qquad \begin {array}[c]{l} {\bf If}\ \, \lambda _i({\bf q}_L) > 0 \ \, {\bf then} \\ \qquad \qquad \qquad \begin {array}[c]{r c l} {\bf \phi }^i_+ & = & {\bf \phi }^i_L \\ {\bf \phi }^i_- & = & 0 \end {array} \\ {\bf else} \\ \qquad \qquad \qquad \begin {array}[c]{r c l} {\rm \phi }^i_+ & = & 0 \\ {\rm \phi }^i_- & = & {\bf \phi }^i_R \end {array} \\ {\bf end if} \end {array} \\ {\bf else} \\ \qquad \qquad \begin {array}[c]{r c l} \alpha ^i & = & {\rm max}(|\lambda _i({\bf q}_L), \lambda _i({\bf q}_R)|) \\ {\bf \phi }^i_+ & = & \frac {1}{2} \left ( {\bf \phi }^i_L + \alpha ^i {\bf w}^i_L \right ) \\ {\bf \phi }^i_- & = & \frac {1}{2} \left ( {\bf \phi }^i_R - \alpha ^i {\bf w}^i_R \right ) \end {array} \\ {\bf end if} \\ \end {array} \\ {\bf end do} \end {array} \end {equation}

Figure 1: The algorithm to calculate the Marquina flux.

Then the numerical flux is given by \begin {equation} \label {marqflux} {\bf f}({\bf q}) = \sum _{i=1}^N \left [ {\bf \phi }^i_+ {\bf r}^i ({\bf q}_L) + {\bf \phi }^i_- {\bf r}^i ({\bf q}_R) \right ]. \end {equation}

The above implementation is based on [1].

8 Other points in GRHydro

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

8.1 Source terms

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

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

For the following computations, we need the expression of some of the 4-Christoffel symbols \(\ {}^{(4)}\Gamma ^\rho _{\mu \nu }\) 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 \begin {equation} \label {eq:SourceADMg} \partial _t \gamma _{ij} = 2\left (- \alpha K_{ij} + \partial _{(i} \beta _{j)} - {}^{(3)}\Gamma ^k_{ij} \beta _k \right )\ . \end {equation} 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 {eqnarray} \label {compatible_derivative} \nabla _i\gamma ^{jk}=\partial _i\gamma ^{jk} + 2{}^{(3)}\Gamma ^j_{il}\gamma ^{lk} = 0 \ . \end {eqnarray}

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

\begin {eqnarray} \label {eq:Gamma000} {}^{(4)}\Gamma ^0_{00} = \frac {1}{2\alpha ^2}\Big [ -\partial _t\big (\beta _k\beta ^k\big )+2\alpha \partial _t\alpha + 2\beta ^i\partial _t\beta _i - \beta ^i\partial _i\big (\beta _k\beta ^k\big ) + 2\alpha \beta ^i\partial _i\alpha \Big ] \end {eqnarray}

and we expand the derivatives as

\begin {eqnarray} \label {de_t_beta2} \partial _t\big (\beta _k\beta ^k\big ) &=& \partial _t\big (\gamma _{jk}\beta ^j\beta ^k\big ) = 2\gamma _{jk}\beta ^j\partial _t\beta ^k + \beta ^j\beta ^k\partial _t\gamma _{jk} = \nonumber \\ &=& 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{}^{(3)}\Gamma ^i_{kj} \beta _i\beta ^j\beta ^k \end {eqnarray}

and

\begin {eqnarray} \label {de_i_beta2} \partial _i\big (\beta _k\beta ^k\big ) = \partial _i\big (\gamma ^{jk}\beta _j\beta _k\big ) = 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{}^{(3)}\Gamma ^j_{ik}\beta _j\beta ^k \ , \end {eqnarray}

where we have used (??) and (??), respectively. Inserting (??) and (??), equation (??) becomes

\begin {eqnarray} \label {eq:Gamma000_final} {}^{(4)}\Gamma ^0_{00} = \frac {1}{\alpha }\Big (\partial _t\alpha + \beta ^i\partial _i\alpha + K_{jk}\beta ^j\beta ^k \Big )\ . \end {eqnarray}

With the same strategy we then compute

\begin {eqnarray} \label {eq:SourceChr1a} {}^{(4)}\Gamma ^0_{i0} & = & - \frac {1}{2\alpha ^2} \Big [ \partial _i (\beta ^k \beta _k - \alpha ^2) - \beta ^j (\partial _i \beta _j - \partial _j \beta _i + \partial _t \gamma _{ij}) \Big ] = - \frac {1}{\alpha } \Big (\partial _i\alpha - \beta ^j K_{ij}\Big ) \end {eqnarray}

and

\begin {eqnarray} \label {eq:SourceChr0ij} {}^{(4)}\Gamma ^0_{ij} & = & - \frac {1}{2\alpha ^2} \Big [ \partial _i\beta _j + \partial _j\beta _i - \partial _t\gamma _{ij} - \beta ^k (\partial _i\gamma _{kj} + \partial _j\gamma _{ki} - \partial _k\gamma _{ij})\Big ] = - \frac {1}{\alpha } K_{ij}\ . \end {eqnarray}

Other more straightforward calculations give \begin {alignat} {3} \label {eq:SourceS3a} {}^{(4)}\Gamma _{00j} &=& {}^{(4)}\Gamma ^\nu _{0j}g_{\nu 0} & = \frac {1}{2} \partial _j \left ( \beta _k \beta ^k - \alpha ^2 \right ), \\ \nonumber \\ \label {eq:SourceS3b} {}^{(4)}\Gamma _{l0j} &=& {}^{(4)}\Gamma ^\nu _{lj}g_{\nu 0} & = \alpha K_{lj} + \partial _l\beta _j + \partial _j\beta _l - \beta _k{}^{(3)}\Gamma ^k_{lj}\ , \\ \nonumber \\ \label {eq:SourceS3c} {}^{(4)}\Gamma _{0lj} &=& {}^{(4)}\Gamma ^\nu _{0j}g_{\nu l} & = -\alpha K_{jl} + \partial _l\beta _j - \beta _k {}^{(3)}\Gamma ^k_{lj}\ , \\ \nonumber \\ \label {eq:SourceS3d} {}^{(4)}\Gamma _{lmj} &=& {}^{(4)}\Gamma ^\nu _{lj}g_{\nu m} & = {}^{(3)}\Gamma _{lmj}\ , \end {alignat}

where (??) was used to derive (??) and (??).

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_{\,k}\) are \begin {equation} \label {eq:SourceS1} \big ({\mathcal S}_{S_k}\big )_j = T^\mu _\nu \Gamma ^\nu _{\mu j} = T^{\mu \nu } \Gamma _{\mu \nu j}\ . \end {equation} After expanding the derivative in (??), the coefficient of the \(T^{\ 00}\) term in (??) becomes

\begin {eqnarray} \label {eq:SourceS4a} {}^{(4)}\Gamma _{00j} & = & \frac {1}{2} \beta ^l \beta ^m \partial _j \gamma _{lm} - \alpha \partial _j \alpha + \beta _m \partial _j \beta ^m. \end {eqnarray}

The coefficient of the \(T^{\,0i}\) term is

\begin {eqnarray} \label {eq:SourceS5a} {}^{(4)}\Gamma _{0ij} + {}^{(4)}\Gamma _{i0j} = \partial _j\beta _i = \beta ^l \partial _i \gamma _{jl} + \gamma _{il} \partial _j \beta ^l. \end {eqnarray}

The coefficient of the \(T^{\,lm}\) term is simply

\begin {eqnarray} \label {eq:SourceS6a} {}^{(3)}\Gamma _{lmj} = \frac {1}{2} \Big (\partial _j\gamma _{ml} + \partial _m\gamma _{jl} - \partial _l\gamma _{mj} \Big ). \end {eqnarray}

Finally, summing (??)–(??) we find

\begin {eqnarray} \label {eq:SourceS2a} \big ({\mathcal S}_{S_k}\big )_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^0_i\partial _j \beta ^i + \frac {1}{2} T^{lm} \partial _j \gamma _{lm} \ , \end {eqnarray}

which is the expression implemented in the code.

8.1.2 Source term for \(\tau \)

The source term for \(\tau \) is [cf. (??)] \begin {equation} \label {eq:SourceT1} {\mathcal S}_{\tau } = \alpha \left ( T^{\mu 0} \partial _{\mu } \alpha - \alpha T^{\mu \nu } {}^{(4)}\Gamma ^0_{\mu \nu }\right ). \end {equation} For clarity, again we consider separately the terms containing as a factor the different components of \(T^{\mu \nu }\). From (??) we find the coefficient of \(T^{\,00}\) to be

\begin {eqnarray} \label {eq:SourceT3a} \alpha \big (\partial _t \alpha -\alpha {}^{(4)}\Gamma ^0_{00}\big ) = -\alpha \big ( \beta ^i \partial _i \alpha + \beta ^k \beta ^l K_{kl}\big )\ . \end {eqnarray}

The coefficient of the term \(T^{\,0i}\) is given by

\begin {eqnarray} \label {eq:SourceT4a} \alpha \big (\partial _i \alpha - 2 \alpha {}^{(4)}\Gamma ^0_{i0}\big ) = 2 \alpha \beta ^j K_{ij} - \alpha \partial _i \alpha \end {eqnarray}

and, finally, the coefficient for \(T^{\,ij}\) is

\begin {eqnarray} \label {eq:SourceT5a} -\alpha ^2 {}^{(4)}\Gamma ^0_{ij} = \alpha K_{ij}\ . \end {eqnarray}

The final expression implemented in the code is thus

\begin {eqnarray} \label {eq:SourceT2a} {\mathcal S}_{\tau } = \alpha \big [ 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}\big ]\ . \end {eqnarray}

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 fluxes and sources we require the primitive variables \(\rho , v_i, P\). Conversion from primitive to conservative is given analytically by equation (??). Converting in the other direction is not possible in a closed form except in certain special circumstances.

There are a number of methods for converting from conservative to primitive variables; see [16]. Here we use a Newton-Raphson type iteration. If we are using a general equation of state such as an ideal gas, then we find a root of the pressure equation. Given an initial guess for the pressure \(\bar {P}\) we find the root of the function \begin {equation} \label {eq:pressure1} f = \bar {P} - P(\bar {\rho }, \bar {\epsilon }), \end {equation} where the approximate density and specific internal energy are given by

\begin {eqnarray} \label {eq:press1gives} \bar {\rho } & = & \frac {\tilde {D}}{\tilde {\tau } + \bar {P} + \tilde {D}} \sqrt { (\tilde {\tau } + \bar {P} + \tilde {D})^2 - S^2 }, \\ \bar {W} & = & \frac {\tilde {\tau } + \bar {P} + \tilde {D}}{\sqrt { (\tilde {\tau } + \bar {P} + \tilde {D})^2 - S^2 }}, \\ \bar {\epsilon } & = & \tilde {D}^{-1} \left ( \sqrt { (\tilde {\tau } + \bar {P} + \tilde {D})^2 - S^2 } - \bar {P} \bar {W} - \tilde {D} \right ). \end {eqnarray}

Here the conserved variables are all “undensitized”, e.g., \begin {equation} \label {eq:undens} \tilde {D} = \gamma ^{-1/2} D, \end {equation} where \(\gamma \) is the determinant of the 3-metric, and \(S^2\) is given by \begin {equation} \label {eq:s2} S^2 = \gamma ^{ij}\tilde {S}_i\tilde {S}_j. \end {equation}

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 \begin {equation} \label {eq:df} f' = 1 - \frac {\partial P}{\partial \rho }\frac {\partial \rho }{\partial P} - \frac {\partial P}{\partial \epsilon }\frac {\partial \epsilon }{\partial P}, \end {equation} where \(\frac {\partial P}{\partial \rho }\) and \(\frac {\partial P}{\partial \epsilon }\) given by calls to EOS_Base, and

\begin {eqnarray} \label {eq:df2} \frac {\partial \rho }{\partial P} & = & \frac {\tilde {D} S^2}{\sqrt {(\tilde {\tau } + \bar {P} + \tilde {D})^2 - S^2}(\tilde {\tau } + \bar {P} + \tilde {D})^2}, \\ \frac {\partial \epsilon }{\partial P} & = & \frac {\bar {P} S^2}{\rho \left ((\tilde {\tau } + \bar {P} + \tilde {D})^2 - S^2\right )(\tilde {\tau } + \bar {P} + \tilde {D})}. \\ \end {eqnarray}

For a polytropic type equation of state, the function is given by \begin {equation} \label {eq:polyf} f = \bar {\rho }\bar {W} - \tilde {D}, \end {equation} where \(\bar {\rho }\) is the variable solved for, the pressure, specific internal energy and enthalpy \(\bar {h}\) are set from the EOS and the Lorentz factor is found from \begin {equation} \label {eq:polyw} \bar {W} = \sqrt {1 + \frac {S^2}{(\tilde {D}\bar {h})^2}}. \end {equation} The derivative is given by \begin {equation} \label {eq:dpolyf} f' = \bar {W} - \frac {\bar {\rho }S^2 \bar {h}'}{\bar {W} \tilde {D}^2 \bar {h}^3}, \end {equation} where \begin {equation} \label {eq:dpolyenth} \bar {h}' = \bar {\rho }^{-1}\frac {\partial P}{\partial \rho }. \end {equation}

8.3 A note on the Roe and Marquina Riemann Solvers

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

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

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

8.4 The atmosphere

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

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

However GRHydro itself also knows about the atmosphere, of course. If the conserved variables \(D\) or \(\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\), \(\epsilon \), \(S_i\) and \(\tau \) are reset to be consistent with \(\rho \) (and \(D\)). Note that even though the polytropic equation of state gives us sufficient information to calculate a consistent value of \(D\), this is not done.

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

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

The reason why the routine that converts to the primitive variables does not ensure that \(D\) is consistent with the other variables is practical rather than accurate. If the value of the variables is set such that they all lie precisely on the atmosphere, then small errors (typically initially of the order of \(10^{-25}\) for a \(64^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 floor) they will result in visible secondary overtones in the oscillations of, e.g., the central density.

The parameters controlling the atmosphere are the following.

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

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

8.5 Advection of passive scalars (’tracers’)

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

\begin {equation} \label {eq:tracer} \partial _t { ( D X_k )} + \partial _{x^j} {\bf f}^{(j)} ({D X_k}) = 0\, , \end {equation} where \(D\) is the generalized particle number density as defined in Eq. (??). GRHydro currently supports any number of independent tracer variables. The following parameters have to be set to use the tracers:

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

8.5.1 Implementation and Limitations

9 History

The approximate time line is something like this:

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

9.1 Thorn Source Code

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

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

9.2 Thorn Documentation

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

9.3 Acknowledgements

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

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

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

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

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

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

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 efficient implementation of flux formulae in multidimensional relativistic hydrodynamical codes. Comput. Phys. Commun., 120:115–121, 1999.

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

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

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

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

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

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

[9]   J. A. Font. Numerical hydrodynamics in General Relativity. Living Rev. Relativity, 3, 2000. [Article in on-line journal], cited on 31/07/01, 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-diffusion equations. J. Comput. Phys., 160:241, 2000.

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

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

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

[21]    P. Mösta et al., ‘GRHydro: A new open source general-relativistic magnetohydrodynamics code for the Einstein Toolkit,” Class. Quant. Grav. 31, 015005 (2014) doi:10.1088/0264-9381/31/1/015005 [arXiv:1304.5544 [gr-qc]].

10 Parameters




constrain_to_1d
Scope: private BOOLEAN



Description: Set fluid 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 fix the carbuncle instability seen at strong shocks



Default: no






atmo_falloff_power
Scope: restricted REAL



Description: The power at which the atmosphere level falls off as (atmo_falloof_radius/r)**atmo_falloff_power



Range Default: 0.0
0:*
Anything positive






atmo_falloff_radius
Scope: restricted REAL



Description: The radius for which the atmosphere starts to falloff as (atmo_falloff_radius/r)**atmo_falloff_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
flat
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 - offset)))”






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_offset
Scope: restricted REAL



Description: The offset 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_bfield
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 fluid 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::DiffRho?



Range Default: First diff
First diff
Standard first differences
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 buffer 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_reflevel
Scope: restricted INT



Description: Start warning on given refinement 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_fix
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 finding 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 specific 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_cutoff
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_reflevel
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 reflevel number) or -1 (off)






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 simplification 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 finite volume method”
Flux Split FD
Finite difference using Lax-Friedrichs flux splitting






min_tracer
Scope: restricted REAL



Description: The floor 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






myfloor
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 off 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 flattening



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_flatten
Scope: restricted KEYWORD



Description: Which flattening procedure should the PPM solver use?



Range Default: stencil_3
stencil_3
our flattening procedure, which requires only stencil 3 and which may work
stencil_4
original C&W flattening 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 modified (enhanced) ppm scheme



Range Default: (none)
0:1
0 (off, default) or 1 (on)






ppm_mppm_debug_eigenvalues
Scope: restricted INT



Description: write eigenvalues into debug grid variables



Range Default: (none)
0:1
0 (off, default) or 1 (on)






ppm_omega1
Scope: restricted REAL



Description: Omega1 for PPM zone flattening



Range Default: 0.75
:
Anything goes. Default is from Colella & Woodward






ppm_omega2
Scope: restricted REAL



Description: Omega2 for PPM zone flattening



Range Default: 10.0
:
Anything goes. Default is from Colella & Woodward






ppm_omega_tracer
Scope: restricted REAL



Description: Omega for tracer PPM zone flattening



Range Default: 0.50
:
Anything goes. Default is from Plewa & Mueller






ppm_small
Scope: restricted REAL



Description: A floor 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 off (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 off (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 flux
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 differencing of the source terms



Range Default: 2
2
2nd order finite differencing
4
4th order finite differencing






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 off
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 off
0:*
damping radius at which we start to damp






track_divb
Scope: restricted BOOLEAN



Description: Track the magnetic field constraint violations



Default: no






transport_constraints
Scope: restricted BOOLEAN



Description: Use constraint transport for magnetic field



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 floor 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_fluxes
Scope: restricted BOOLEAN



Description: Weight the flux terms by the cell surface areas



Default: no






weno_adaptive_epsilon
Scope: restricted BOOLEAN



Description: use modified 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






use_mask
Scope: shared from SPACEMASKBOOLEAN



11 Interfaces

General

Implements:

grhydro

Inherits:

admbase

boundary

spacemask

tmunubase

hydrobase

Grid Variables

11.0.1 PRIVATE GROUPS





  Group Names     Variable Names   Details    




inlastmolpoststep InLastMoLPostStep compact 0
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 execute_MoL_Step compact 0
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 execute_MoL_PostStep compact 0
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 whichpsidcspeed compact 0
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




flux_splitting compact 0
densfplus description Fluxes for use in the flux 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 flux splitting
fs_alpha2 dimensions 0
fs_alpha3 distribution CONSTANT
fs_alpha4 group type SCALAR
fs_alpha5 timelevels 1
variable type REAL




y_e_plus Y_e_plus compact 0
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 Y_e_minus compact 0
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 tempplus compact 0
description Plus state for the temperature
dimensions 3
distribution DEFAULT
group type GF
tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




tempminus tempminus compact 0
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_flux compact 0
cons_tracerflux 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_flux_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 flux 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 GRHydro_maxima_iteration compact 0
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




diffrho DiffRho compact 0
description The first difference 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 dens compact 0
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 tau compact 0
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 scon compact 0
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 Bcons compact 0
description B-field 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 Evec compact 0
description Electric field 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 Y_e_con compact 0
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 entropycons compact 0
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 sdetg compact 0
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 psidc compact 0
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 densrhs compact 0
description Update term for dens
dimensions 3
distribution DEFAULT
group type GF
tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




taurhs taurhs compact 0
description Update term for tau
dimensions 3
distribution DEFAULT
group type GF
tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




srhs srhs compact 0
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 Bconsrhs compact 0
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 Avecrhs compact 0
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 Aphirhs compact 0
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 psidcrhs compact 0
description Update term for psidc
dimensions 3
distribution DEFAULT
group type GF
tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




entropyrhs entropyrhs compact 0
description Update term for entropycons
dimensions 3
distribution DEFAULT
group type GF
tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




divb divB compact 0
description Magnetic field constraint
dimensions 3
distribution DEFAULT
group type GF
tags Prolongation=”Restrict” checkpoint=”no” tensorparity=-1
timelevels 1
variable type REAL




bcom bcom compact 0
description bî: comoving contravariant magnetic field 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 bcom0 compact 0
description bˆ0  component of the comoving contravariant magnetic field 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 bcom_sq compact 0
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 lvel compact 0
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 lBvec compact 0
description local Magnetic field 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
flux_direction description Which direction are we taking the fluxes in and the offsets
xoffset dimensions 0
yoffset distribution CONSTANT
zoffset 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_field_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 maxrho_global compact 0
description store the global maximum of rho
  description for refinment-grid steering
dimensions 0
distribution CONSTANT
group type SCALAR
tags checkpoint=”no”
timelevels 1
variable type REAL




grhydro_c2p_failed GRHydro_C2P_failed compact 0
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_fluxes compact 0
densflux description Fluxes for each conserved variable
sxflux dimensions 3
syflux distribution DEFAULT
szflux group type GF
tauflux tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




grhydro_bfluxes compact 0
Bconsxflux description Fluxes for each B-field variable
Bconsyflux dimensions 3
Bconszflux distribution DEFAULT
group type GF
tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




grhydro_psifluxes compact 0
psidcflux 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_entropyfluxes compact 0
entropyflux 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_avecfluxes compact 0
Avecxflux description Fluxes for each Avec variable
Avecyflux dimensions 3
Aveczflux distribution DEFAULT
group type GF
tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




grhydro_aphifluxes compact 0
Aphiflux description Fluxes for Aphi
dimensions 3
distribution DEFAULT
group type GF
tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




evolve_y_e evolve_Y_e compact 0
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 evolve_temper compact 0
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 evolve_entropy compact 0
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 evolve_MHD compact 0
description Are we doing MHD? Set in ParamCheck
dimensions 0
distribution CONSTANT
group type SCALAR
tags checkpoint=”no”
timelevels 1
variable type INT








  Group Names     Variable Names   Details    




evolve_lorenz_gge evolve_Lorenz_gge compact 0
description Are we evolving the Lorenz gauge?
dimensions 0
distribution CONSTANT
group type SCALAR
tags checkpoint=”no”
timelevels 1
variable type INT




grhydro_reflevel GRHydro_reflevel compact 0
description Refinement level GRHydro is working on right now
dimensions 0
distribution CONSTANT
group type SCALAR
tags checkpoint=”no”
timelevels 1
variable type INT




y_e_con_rhs Y_e_con_rhs compact 0
description RHS for the electron fraction
dimensions 3
distribution DEFAULT
group type GF
tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




y_e_con_flux Y_e_con_flux compact 0
description Flux for the electron fraction
dimensions 3
distribution DEFAULT
group type GF
tags Prolongation=”None” checkpoint=”no”
timelevels 1
variable type REAL




Uses header:

SpaceMask.h

carpet.hh

Provides:

SpatialDet to

UpperMet to

Con2PrimPoly to

Con2PrimGenM to

Con2PrimGenMee to

Con2PrimGen to

Con2PrimPolyM to

Prim2ConGen to

Prim2ConPoly to

Prim2ConGenM to

Prim2ConGenM_hot to

Prim2ConPolyM to

12 Schedule

This section lists all the variables which are assigned storage by thorn EinsteinEvolve/GRHydro. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Storage

 

Always: Conditional:
execute_MoL_Step execute_MoL_PostStep
evolve_MHD HydroBase::temperature[timelevels]
evolve_Y_e HydroBase::entropy[timelevels]
evolve_temper HydroBase::Bvec[timelevels]
GRHydro_reflevel GRHydro::Bcons[timelevels]
InLastMoLPostStep Bconsrhs
sdetg psidc[timelevels]
densrhs HydroBase::Bvec[timelevels]
taurhs HydroBase::Avec[timelevels]
srhs HydroBase::Aphi[timelevels]
GRHydro_C2P_failed[1] dens[timelevels]
  Bconsrhs
  psidcrhs
  tau[timelevels]
  whichpsidcspeed
  divB
  GRHydro::bcom[timelevels]
  GRHydro::bcom0[timelevels]
  GRHydro::bcom_sq[timelevels]
  Avecrhs
  evolve_Lorenz_gge
  Aphirhs
  divB
  GRHydro::bcom[timelevels]
  scon[timelevels]
  GRHydro::bcom0[timelevels]
  GRHydro::bcom_sq[timelevels]
  HydroBase::entropy[timelevels]
  GRHydro::entropycons[timelevels]
  entropyrhs
  evolve_MHD
  evolve_Y_e
  evolve_temper
  evolve_entropy
  GRHydro_reflevel
  particles[timelevels]
  InLastMoLPostStep
  densrhs
  taurhs
  srhs
  GRHydro_eos_scalars
  GRHydro_minima
  GRHydro_scalars
  particle_rhs
  particle_arrays
  GRHydro_tracers[timelevels]
  Y_e_con[timelevels]
  GRHydro_cons_tracers[timelevels]
  GRHydro_tracer_rhs
  Evec
  ADMBase::metric[timelevels] ADMBase::curv[timelevels]
  ADMBase::lapse[timelevels]
  ADMBase::shift[timelevels]
  GRHydro_coords[timelevels]
  GRHydro_coords_rhs
  lvel[timelevels]
  lBvec[timelevels]
  Y_e_con_rhs Y_e_con_flux Y_e_plus Y_e_minus
  local_metric[timelevels]
  local_extrinsic_curvature
  local_shift
  H_viscosity_temps
  GRHydro_atmosphere_mask
  GRHydro_atmosphere_mask_real
  GRHydro_atmosphere_descriptors
  fs_alpha
  GRHydro_tracers[timelevels]
  GRHydro_cons_tracers[timelevels]
  HydroBase::Y_e[timelevels]
  GRHydro_tracer_rhs
  GRHydro_trivial_rp_gf_group
  GRHydro_mppm_eigenvalues
  tempplus tempminus
   

Scheduled Functions

CCTK_BASEGRID

  grhydro_reset_execution_flags

  initially set execution flags to ’yeah, execute’!

 

  Language: c
  Options: global
  Type: function

MoL_Step (conditional)

  grhydro_set_execution_flags

  check if we need to execute rhs / post-step calculation

 

  After: mol_decrementcounter
  Before: mol_poststepmodify
  Language: c
  Options: level
  Type: function

GRHydroRHS (conditional)

  grhydro_evolvecoords

  evolve the coordinates for the comoving shift

 

  Language: fortran
  Type: function

HydroBase_Prim2ConInitial (conditional)

  primitive2conservativecells

  convert initial data given in primive variables to conserved variables

 

  Language: fortran
  Type: function

HydroBase_Con2Prim (conditional)

  conservative2primitivepolytypem

  convert back to primitive variables (polytype) - mhd version

 

  If: grhydro::execute_mol_poststep
  Language: fortran
  Type: function

HydroBase_Prim2ConInitial (conditional)

  primitive2conservativepolycellsm

  convert initial data given in primive variables to conserved variables - mhd version

 

  Language: fortran
  Type: function

HydroBase_Con2Prim (conditional)

  conservative2primitivepolytypeam

  convert back to primitive variables (polytype) - mhd with avec version

 

  If: grhydro::execute_mol_poststep
  Language: fortran
  Type: function

HydroBase_Prim2ConInitial (conditional)

  primitive2conservativepolycellsam

  convert initial data given in primive variables to conserved variables - mhd with avec version

 

  Language: fortran
  Type: function

HydroBase_Con2Prim (conditional)

  conservative2primitivepolytype

  convert back to primitive variables (polytype)

 

  If: grhydro::execute_mol_poststep
  Language: fortran
  Type: function

HydroBase_Prim2ConInitial (conditional)

  primitive2conservativepolycells

  convert initial data given in primive variables to conserved variables

 

  Language: fortran
  Type: function

HydroBase_Boundaries (conditional)

  do_grhydro_boundaries

  grhydro boundary conditions group

 

  Type: group

HydroBase_PostStep (conditional)

  grhydro_atmospheremaskboundaries

  apply boundary conditions to primitives

 

  Before: hydrobase_boundaries
    grhydro_primitiveinitialguessesboundaries
  Type: group

GRHydro_AtmosphereMaskBoundaries (conditional)

  grhydro_selectatmospheremaskboundaries

  select atmosphere mask for boundary conditions

 

  Language: fortran
  Options: level
  Sync: grhydro_atmosphere_mask_real
  Type: function

HydroBase_PostStep (conditional)

  grhydro_comovingshift

  comoving shift

 

  After: hydrobase_con2prim
  Language: fortran
  Type: function

GRHydro_AtmosphereMaskBoundaries (conditional)

  applybcs

  apply boundary conditions to real-valued atmosphere mask

 

  After: grhydro_selectatmospheremaskboundaries
  Type: group

HydroBase_Select_Boundaries (conditional)

  grhydro_boundaries

  select grhydro boundary conditions

 

  If: grhydro::execute_mol_poststep
  Language: fortran
  Options: level
  Sync: dens
    tau
    scon
    hydrobase::w_lorentz
    hydrobase::rho
    hydrobase::press
    hydrobase::eps
    hydrobase::vel
    bcons
    entropycons
    hydrobase::bvec
    psidc
    grhydro_cons_tracers
    grhydro_tracers
    hydrobase::temperature
    hydrobase::entropy
    hydrobase::y_e
    y_e_con
    lvel
    lbvec
  Type: function

HydroBase_Select_Boundaries (conditional)

  grhydro_boundaries

  select grhydro boundary conditions

 

  If: grhydro::execute_mol_poststep
  Language: fortran
  Options: level
  Sync: dens
    tau
    scon
    bcons
    entropycons
    psidc
    grhydro_cons_tracers
    y_e_con
  Type: function

CCTK_POSTREGRID (conditional)

  grhydro_primitiveboundaries

  apply boundary conditions to all primitives

 

  Before: mol_poststep
  Type: group

CCTK_POSTREGRIDINITIAL (conditional)

  grhydro_primitiveboundaries

  apply boundary conditions to all primitives

 

  Before: mol_poststep
  Type: group

HydroBase_PostStep (conditional)

  grhydro_primitiveinitialguessesboundaries

  apply boundary conditions to those primitives used as initial guesses

 

  Before: hydrobase_boundaries
  If: grhydro::inlastmolpoststep
  Type: group

GRHydro_PrimitiveInitialGuessesBoundaries (conditional)

  grhydro_selectprimitiveinitialguessesboundaries

  select initial guess primitive variables for boudary conditions

 

  Language: fortran
  Options: level
  Sync: hydrobase::press
    hydrobase::rho
    hydrobase::eps
    hydrobase::temperature
    lvel
    lbvec
  Type: function

GRHydro_PrimitiveBoundaries (conditional)

  grhydro_selectprimitiveboundaries

  select primitive variables for boundary conditions

 

  Language: fortran
  Options: level
  Sync: hydrobase::press
    hydrobase::entropy
    hydrobase::y_e
    grhydro_tracers
    hydrobase::rho
    hydrobase::eps
    hydrobase::temperature
    lvel
    lbvec
  Type: function

GRHydro_PrimitiveInitialGuessesBoundaries (conditional)

  grhydro_selectprimitiveinitialguessesboundaries

  select initial guess primitive variables for boudary conditions

 

  Language: fortran
  Options: level
  Sync: hydrobase::press
    hydrobase::rho
    hydrobase::eps
    hydrobase::vel
    hydrobase::bvec
    hydrobase::temperature
  Type: function

GRHydro_PrimitiveBoundaries (conditional)

  grhydro_selectprimitiveboundaries

  select primitive variables for boundary conditions

 

  Language: fortran
  Options: level
  Sync: hydrobase::press
    hydrobase::entropy
    hydrobase::y_e
    grhydro_tracers
    hydrobase::rho
    hydrobase::eps
    hydrobase::vel
    hydrobase::bvec
    hydrobase::temperature
  Type: function

CCTK_STARTUP (conditional)

  grhydro_startup

  startup banner

 

  Language: fortran
  Type: function

GRHydro_PrimitiveInitialGuessesBoundaries (conditional)

  applybcs

  apply boundary conditions to initial guess primitive variables

 

  After: grhydro_selectprimitiveinitialguessesboundaries
  Type: group

GRHydro_PrimitiveBoundaries (conditional)

  applybcs

  apply boundary conditions to all primitive variables

 

  After: grhydro_selectprimitiveboundaries
  Type: group

MoL_PostStep (conditional)

  grhydro_setlastmolpoststep

  set grid scalar inlastmolpoststep if this is the last mol poststep call

 

  Language: c
  Options: level
  Type: function

MoL_Step (conditional)

  grhydro_clearlastmolpoststep

  reset inlastmolpoststep to zero

 

  After: mol_poststep
  Language: c
  Options: level
  Type: function

CCTK_WRAGH (conditional)

  grhydro_clearlastmolpoststep

  initialize inlastmolpoststep to zero

 

  Language: c
  Options: global-early
  Type: function

HydroBase_PostStep (conditional)

  grhydropostsyncatmospheremask

  set integer atmosphere mask from synchronized real atmosphere mask

 

  After: grhydro_atmospheremaskboundaries
  Language: fortran
  Type: function

HydroBase_PostStep (conditional)

  grhydro_atmosphereresetm

  reset the atmosphere - mhd version

 

  After: grhydropostsyncatmospheremask
  Before: hydrobase_boundaries
    grhydro_primitiveinitialguessesboundaries
  If: grhydro::inlastmolpoststep
  Language: fortran
  Type: function

HydroBase_PostStep (conditional)

  grhydro_atmosphereresetam

  reset the atmosphere - mhd with avec version

 

  After: grhydropostsyncatmospheremask
  Before: hydrobase_boundaries
    grhydro_primitiveinitialguessesboundaries
  If: grhydro::inlastmolpoststep
  Language: fortran
  Type: function

HydroBase_PostStep (conditional)

  grhydro_atmospherereset

  reset the atmosphere

 

  After: grhydropostsyncatmospheremask
  Before: hydrobase_boundaries
    grhydro_primitiveinitialguessesboundaries
  If: grhydro::inlastmolpoststep
  Language: fortran
  Type: function

CCTK_ANALYSIS (conditional)

  grhydro_check_rho_minimum

  check whether somewhere rho(i,j,k) < grhydro_rho_min and produce a warning

 

  Language: fortran
  Triggers: hydrobase::rho
    hydrobase::press
    hydrobase::eps
  Type: function

CCTK_STARTUP (conditional)

  grhydro_registermask

  register the hydro masks

 

  Language: c
  Type: function

CCTK_ANALYSIS (conditional)

  grhydro_refinementlevel

  calculate current refinement level

 

  Before: grhydro_check_rho_minimum
  Language: fortran
  Triggers: hydrobase::rho
    hydrobase::press
    hydrobase::eps
  Type: function

CCTK_BASEGRID

  reset_grhydro_c2p_failed

  initialise the mask function that contains the points where c2p has failed (at basegrid)

 

  Language: fortran
  Type: function

CCTK_PRESTEP

  reset_grhydro_c2p_failed

  reset the mask function that contains the points where c2p has failed (at prestep)

 

  Language: fortran
  Type: function

CCTK_EVOL

  sync_grhydro_c2p_failed

  syncronise the mask function that contains the points where c2p has failed

 

  After: mol_evolution
  Language: fortran
  Sync: grhydro_c2p_failed
  Type: function

CCTK_POSTSTEP

  check_grhydro_c2p_failed

  check the mask function that contains the points where c2p has failed and report an error in case a failure is found

 

  After: grhydro_refinementlevel
  Language: fortran
  Type: function

AddToTmunu (conditional)

  grhydro_tmunum

  compute the energy-momentum tensor - mhd version

 

  Language: fortran
  Type: function

HydroBase_PostStep (conditional)

  grhydro_calcbcom

  compute comoving magnetic field, pressure, etc...

 

  After: grhydrotransformprimtoglobalbasis
  Language: fortran
  Type: function

AddToTmunu (conditional)

  grhydro_tmunu

  compute the energy-momentum tensor

 

  Language: fortran
  Type: function

CCTK_POSTPOSTINITIAL (conditional)

  settmunu

  calculate the stress-energy tensor

 

  After: con2prim
  Before: admconstraintsgroup
  Type: group

MoL_PseudoEvolution (conditional)

  grhydroanalysis

  calculate analysis quantities

 

  Type: group

CCTK_PARAMCHECK (conditional)

  grhydro_paramcheck

  check parameters

 

  Language: fortran
  Type: function

GRHydroAnalysis (conditional)

  grhydro_calcdivb

  calculate divb

 

  After: grhydro_analysis_init
  Language: fortran
  Type: function