GRMHD in dynamical spacetimes

Major changes

This long-delayed merge request adds a new module for GRMHD in dynamical spacetimes. The number of changes are quite extensive, but they are mostly independent of functionality elsewhere in the code. Here is a brief summary of the most important changes:

• The dyngr folder contains code for DynGR, the new GRMHD solver based on the Valencia formulation so that it is suitable for dynamical spacetimes. The solver tries to piggyback off of the existing MHD routines where possible, but it implements its own Riemann solvers and source terms, and the EOS-agnostic implementation introduces some subtle complications that make it incompatible with the MHD flux and FOFC methods.
• The CornerE function now contains extra branches to handle the Valencia solver.
• The tmunu folder contains code for the new Tmunu object for storing the stress-energy tensor in its 3+1 decomposed form (E, S_i, S_{ij}) for use in the Z4c RHS. This is slightly different from GR-Athena++, which stores an object mat inside Z4c.
• The numerical-relativity folder is for the NumericalRelativity object, which acts as a task orchestrator for DynGR and Z4c. Rather than creating their own TaskList objects directly, DynGR and Z4c both submit tasks to TaskQueue objects in NumericalRelativity, each of which has a set of required dependencies and a set of optional dependencies which are conditionally enabled depending on the other physics modules in the code (e.g., Z4c_CalcRHS has a required dependency on Z4c_CopyU and an optional dependency on MHD_SetTmunu which is only enabled if MHD is also active). NumericalRelativity then loops over each TaskQueue and attempts to build a TaskList from it which satisfies all the correct dependencies. This infrastructure seems excessive at the moment, but it will greatly simplify things as we add additional physics.
• The eos/primitive-solver folder includes PrimitiveSolver, a library for handling EOS calls and the conserved-to-primitive inversion for DynGR. This is necessary because the existing EOS infrastructure in AthenaK assumes a gamma-law equation of state, and we need support for generic equations of state. Right now there is robust support for an ideal gas, but a hybrid piecewise-polytropic EOS has been lightly tested, and there is preliminary support for a tabulated finite-temperature EOS.
• Several problem generators have been modified to support DynGR. Most of this work involves assigning pressure as a primitive variable instead of internal energy density (the former tends to be better behaved for reconstruction when using tabulated equations of state) and instantiating the ADM metric, but a few also have new features for testing DynGR.
• Several new problem generators have been added, such as dyngr_tov, lorene_bns, elliptica_tov, etc.
• The Coordinates class has been modified to support updating the excision mask dynamically during runtime according to a user-specified excision scheme in the parameter file via coords/excision_scheme. The default scheme, fixed, maintains the original behavior of the code, so existing tests requiring excision should not be affected by this change. The other currently implemented scheme, lapse, is an ad-hoc scheme that applies excision everywhere the lapse falls below a specified cutoff value.
• The utils folder now includes code for a new class, TableReader, which is currently used by the tabulated EOS implementation. However, the structure of the tables (ASCII header + binary data) is generic enough that there's no reason it can't be used elsewhere, too.

Unfinished tasks

There are a number of tasks still left to do. However, the code has introduced enough new features that I figured it would be better to merge back sooner rather than later. Here's what is still left:

• There are no automated CI tests for DynGR. At the bare minimum we need a test that probes the GRMHD alone and a test that probes the GRMHD+Z4c, but it would also be useful to have tests for the tabulated EOS.
• DynGR uses pressure as a primitive variable, but it is still labeled eint in output files. There needs to be logic to handle this properly. Perhaps it can be done when the output infrastructure is refactored into something more maintainable.
• The code currently passes the CI pipeline, but a handful of pgens (blast, kh, gr_torus, etc.) not included in the built-in tests need to be validated more rigorously to ensure that changes for DynGR didn't cause any unexpected side effects in other cases.
• It may be worth separating out the TOV solver into utils. There are a variety of TOV pgens floating around, both within and without the repository, and there's a lot of unnecessary code duplication going on.

Validation

Here are a few tests showing that the code works as expected:

• MHD shock tubes (Balsara 1 and 4)

• Power spectrum of TOV in Cowling approximation with linear perturbation (vertical lines show predicted oscillation frequencies from perturbation theory)

• Power spectrum of TOV in full GR (vertical lines show predicted oscillation frequencies from perturbation theory)

• Unmagnetized binary neutron star merger

The BNS test achieves approximate second-order convergence of the waveform phase during the inspiral as expected. The merger time differs from GR-Athena++ by \Delta t \approx 0.05~\mathrm{ms}. The results do diverge somewhat in the post-merger phase, particularly in predicted collapse time, but this a combination of insufficient resolution and the extreme sensitivity of collapse to initial conditions.