Implement curved spacetime hydrodynamics

[JonathanGorard]

The special relativistic 5-moment equations currently implemented in Gkeyll permit a straightforward extension to curved spacetimes (i.e. to fully general relativistic hydrodynamics in the test fluid limit) by means of the fully-conservative, strongly-hyperbolic "Valencia" formulation of Banyuls, Font, Ibáñez, Martí and Miralles. The only differences between these curved spacetime hydrodynamics equations and the flat spacetime case are a gauge dependence (i.e. an appearance of the lapse function and the shift vector) within the flux vector, and the appearance of source terms proportional to certain derivatives of the Christoffel symbols. Both of these differences can be absorbed into a small modification of the finite-volume boundary extrapolation step, such that the equations that one is solving at the cell centers remain purely special relativistic.

Subsonic, steady-state accretion of a stiff relativistic fluid (i.e. an ultrarelativistic fluid with Gamma = 2, i.e. where P = rho) onto a Schwarzschild or Kerr black hole, quite surprisingly, admits an analytic solution due to Petrich, Shapiro and Teukolsky, against which one can therefore make direct numerical comparisons for the sake of unit and regression testing. Having a robust and well-validated general relativistic hydrodynamics model (at least in the test fluid limit) implemented within Gkeyll is a crucial initial step towards being able to do full GRMHD and even kinetic simulations in curved (static) spacetimes. With this in mind, the implementation steps that I currently have in mind are:

1. Increase the robustness of the current Gkeyll implementation of the Eulderink-Mellema iterative scheme for doing conservative -> primitive variable reconstruction in SRHD, especially in the presence of strongly relativistic shocks, low densities and low pressures. This can be achieved by first modifying the formula for extracting the Lorentz factor W by pulling the (C_0 * xi) term outside of the brackets, since this term becomes arbitrarily large for strongly relativistic shocks, and the term inside the brackets should be kept at order unity if possible. From there, one can enforce a minimum value of the fluid density, and then recalculate the pressure from the (potentially modified value of the) specific relativistic enthalpy. This scheme can be validated against both the mildly and strongly relativistic blast wave tests of Del Zanna and Bucciatini, since the latter problem in particular is an extremely effective and demanding test for assessing the robustness of the primitive reconstruction scheme. Riemann problems for the one-dimensional relativistic hydrodynamics equations admit an analytic solution due to Pons, Martí and Müller, against which we can write appropriate unit and regression test comparisons (along with the 2D quadrants test, which currently fails). Though not strictly a Riemann problem, we can also perform comparisons against the perturbed density test, since this will be crucial for assessing the extent to which we are able to resolve fine, oscillatory features.
2. Implement an abstract context struct for storing relevant spacetime data (for the case of black hole spacetimes, this would likely consist of parameters for mass, spin, position of the black hole, and possibly linear momentum of black hole if we intend to support non-comoving frames), which is then passed to all wv_gr_euler functions.
3. Add helper functions for computing relevant spacetime quantities on the basis of this context struct; in particular, we will need to add functions for calculating the components of the spatial and spacetime metrics (plus their determinants, and their inverses, for raising and lowering indices), the lapse function (and its derivatives), the shift vector (and its derivatives), and the components of the extrinsic curvature tensor. For the case of black hole spacetimes in Kerr-Schild coordinates, these quantities all admit known analytic solutions, which we can implement directly. Ideally, we should aim to implement this in a manner that is sufficiently general as to admit most or all standard coordinatization schemes for black hole spacetimes (e.g. Boyer-Lindquist/oblate spheroidal coordinates, general horizon-adapted coordinates, Cartesian Kerr-Schild coordinates, spherical Kerr-Schild coordinates, etc.)
4. Implement excision boundary conditions at the event horizon, in the case that one chooses not to use horizon-adapted coordinates. This can be achieved either by using the existing (provisional) multi-block infrastructure to excise a given region of the domain directly in computational coordinates, and then impose physical boundary conditions at the excision point, or by introducing a fictitious scalar field that delineates the interior and exterior spacetime regions (and exploiting the sign change in this field to allow the solver to detect the event horizon efficiently, similar to a level set method). I am currently undecided regarding which of these two approaches is preferable.
5. Implement appropriate modifications of the Roe, HLL and HLLC solvers for the curved spacetime Euler case (incorporating gauge and metric dependence in the boundary extrapolation step for the Riemann solver), along with wave speed estimates based on the known eigendecomposition of the full GRHD system due to Banyuls et al.
6. Perform comparisons against analytic spherical Bondi-Hoyle and non-spherical Hoyle-Lyttleton accretion problems onto both Schwarzschild and Kerr black holes, comparing against the analytic solutions of Petrich et al. in the stiff relativistic case, and the numerical solutions of Font, Ibáñez and Papadopoulos in the more general ultrarelativistic case in axisymmetry.

[ammarhakim]

This looks like an excellent plan, specially given the work we will pursue this summer on the astro front. If we have the hyrdo, Maxwell and implicit sources (Lorentz + currents) completed this summer then we will be set to do some really cool, and totally unique, simulations of the plasma around compact objects!