IllinoisGRMHD: An Open-Source, User-Friendly GRMHD Code for Dynamical Spacetimes

• Without detailed electromagnetic-counterpart models, interpreting observed counterparts may be limited; constructing such models requires resolving many length and time scales and coupling nonlinear Einstein and GRMHD equations.Source: §1 · p. 2
• OrigGRMHD’s vector-potential AMR formulation preserves the no-monopoles constraint across refinement boundaries, supporting modeling of multi-scale magnetized flows.Source: §1 · p. 2
• The community could benefit from consolidating development around a shared dynamical-spacetime GRMHD code, while GRHydro’s AMR monopole issue and reliance on black-hole excision have limited its broader use.Source: §1 · pp. 3–4
• Despite OrigGRMHD’s track record, its limited documentation, modularity, coding consistency, and maintenance hindered community adoption, prompting a ground-up rewrite guided by four design principles.Source: §1 · p. 3
• Validation tests show IllinoisGRMHD reproduces OrigGRMHD, matches or improves its performance and scalability relative to OrigGRMHD and GRHydro, and agrees quantitatively with GRHydro in dynamical-spacetime TOV-star evolutions.Source: §1 · p. 4

2. Basic Equations

• IllinoisGRMHD couples Einstein gravity to ideal-MHD matter, assuming a perfect-fluid stress-energy tensor and infinite conductivity, with GRMHD quantities sourcing Einstein’s equations.Source: §2 · p. 5
• For numerical evolution, the spacetime is split into ADM 3+1 form and evolved with BSSN, a formulation chosen for stable, long-term, accurate integration.Source: §2 · pp. 5–6
• The GRMHD equations use conservative form with high-resolution shock capturing, conserving rest mass to roundoff without source terms and satisfying shock jump conditions.Source: §2 · pp. 6–7
• Directly discretizing the induction equation can violate the divergence-free magnetic-field constraint, especially on AMR grids; evolving the vector potential preserves it to roundoff error.Source: §2 · p. 7
• The hybrid equation of state combines cold and thermal components, with Γ_th setting kinetic-to-thermal conversion at shocks; the conservative-to-primitive solver supports only one polytropic index.Source: §2 · p. 8

2.1. Outer Boundary Conditions

• At outer faces, density, pressure, and velocity use copied zero-gradient values, with incoming velocity set to zero to impose outflow conditions.Source: §2.1 · p. 9

3. Basic Algorithms

• IllinoisGRMHD removes obsolete and nonrobust functionality, standardizes core algorithms, and replaces duplicated code with a compact library of optimized functions.The retained algorithms are those used in the Illinois group’s latest GRMHD publications.Source: §3 · p. 9
• Conservative-to-primitive recovery uses a Newton-Raphson root-finder to compute primitive variables from updated conservatives, with consistency checks before and after the solver.This recovery restores primitive variables needed for the next iteration or timestep.Source: §3 · pp. 10–11

3.1. Evolution of ρ∗, ˜Si, and ˜τ

• The conservative variables ρ∗, ˜Si, and ˜τ are evolved with a second-order finite-volume HRSC scheme using PPM reconstruction and an HLL approximate Riemann solver.The approximate characteristic speeds can overestimate the maximum by a factor ≤2, increasing dissipation that may smear physical features but helps stabilize evolutions; RK4 is typically chosen to minimize total error in spacetime and fluid evolution.Source: §3.1 · pp. 11–13

3.2. Vector-Potential-Based Constrained Transport Scheme

• Vector-potential constrained transport maintains ∇·B = 0 and yields identical output to standard staggered constrained transport on uniform grids, while allowing any AMR-boundary interpolation scheme if B is computed from interpolated A.Standard constrained transport on AMR requires specially devised, implementation-tuned prolongation and restriction to maintain the divergence constraint at refinement boundaries.Source: §3.2 · p. 13

3.3. Evolution of the densitized EM scalar potential [√γΦ]

• For the densitized electromagnetic scalar potential, IllinoisGRMHD interpolates right-hand-side quantities to consistent staggered points before taking derivatives, rather than applying reconstruction.The shift-advection term is handled by first interpolating the shift to the required staggered point.Source: §3.3 · p. 16
• For the densitized EM scalar-potential derivative, advection uses the metric thorn’s upwinding strategy; Term (3) requires interpolating only α using Term (1)’s strategy.The derivative is evaluated with second-order operators, with corresponding derivatives in the y and z directions following the x-direction procedure.Source: §3.3 · pp. 16–17

3.4. Conservatives-to-Primitives Solver

• IllinoisGRMHD uses a two-dimensional Newton-Raphson solver to recover primitive variables, validating conservatives first and setting pressure to Pcold if the solver rarely fails.Failures occur almost always in a low-density atmosphere or inside a black hole.Source: §3.4 · p. 17
• After inversion, superluminal speeds are reduced to W=10, pressure is floored at 0.9Pcold, and capped at 100Pcold for ρ0<100ρatm or deep inside horizons.Raising the pressure floor to Pcold exactly causes large central density drifts in TOV star evolutions.Source: §3.4 · pp. 17–18

3.5. Outer Boundary Conditions for Ai, [√γΦ], P, ρ0, and vi

• Each RK4 iteration uses three ghost zones; vector potentials receive linear-extrapolation boundaries, while hydrodynamic variables use zero-derivative outflow boundaries.With xy-plane symmetry, the negative-z outer condition is omitted so Cactus/Carpet can impose reflection symmetry.Source: §3.5 · p. 18

3.6. Linkage of IllinoisGRMHD to the Rest of the Einstein Toolkit

• Through TmunuBase, IllinoisGRMHD is compatible with Einstein Toolkit spacetime formulations coupled to ADMBase, including BSSN and conformal and covariant Z4 [1].The paper’s dynamical-spacetime evolutions couple IllinoisGRMHD to the McLachlan BSSN thorn [15].Source: §3.6 · pp. 18–19

4.1. IllinoisGRMHD and OrigGRMHD: Roundoff-Error Agreement in Evolving Magnetized TOV Star

• On AMR grids, IllinoisGRMHD and OrigGRMHD track the same roundoff-error curve when evolving a weakly magnetized TOV star.The test compares their agreement with the expected roundoff error from a 15th-significant-digit perturbation.Source: §4.1 · p. 20, Figure 1 · p. 21

4.2. IllinoisGRMHD and OrigGRMHD: Roundoff-Error Agreement in Evolving Random Initial Data

• The random-data test uses discontinuous, stochastic GRMHD flows and spacetime fields, while holding the spacetime fixed rather than evolving Einstein's equations.The setup was designed as a numerical stress test and is not intended to satisfy Einstein's equations.Source: §4.2 · p. 20
• In that harsh random-data test, IllinoisGRMHD and OrigGRMHD agree within expected roundoff error for at least 25 light-crossing times.The OrigGRMHD perturbation baseline plateaus at about 13 digits of agreement.Source: §4.2 · p. 21, Figure 2 · p. 22

4.3. IllinoisGRMHD and GRHydro: Unmagnetized TOV Star Convergence Tests

• The GRHydro comparison uses identical initial data, grids, spacetime evolution, reconstruction, and Riemann-solver choices, making it a proxy comparison with OrigGRMHD.Both codes use the same basic algorithms in this test but differ in implementation.Source: §4 · p. 19, §4.3 · p. 23
• Both codes' central-density drift and Hamiltonian constraint violations converge to zero at roughly second order, while GRHydro's constraint-convergence order fluctuates more.At the highest resolution, GRHydro has roughly 8% higher Hamiltonian constraint violation than IllinoisGRMHD.Source: §4.3 · p. 23, Figure 3 · p. 24
• For the unmagnetized TOV test, IllinoisGRMHD's absolute central-density drift is higher, but the fitted drift rates are within about one standard deviation.The fits use data from 5 to 55 dynamical timescales; doubling the outer boundary had no qualitative effect.Source: §4.3 · pp. 23–24, Figure 4 · p. 25

5. Performance Benchmarks

• IllinoisGRMHD's staggered vector-potential formulation costs about 60% more for flux reconstructions, but prevents monopoles from being generated at AMR boundaries.The cited comparison says GRHydro's AMR scheme cannot guarantee this constraint.Source: §5 · pp. 25–26
• IllinoisGRMHD outperforms OrigGRMHD by 1.3× at 32 cores and 1.6× at 256 and 2,048 cores.Independent tests attribute about 10–20% of the performance difference to OrigGRMHD's older Cactus/Carpet infrastructure.Source: §5 · p. 26
• The benchmarks project a 1.6–1.8× computational-resource boost for groups adopting IllinoisGRMHD from standard GRHydro or OrigGRMHD.Independent tests indicate a performance gap of about 2× for fixed-spacetime-background GRMHD simulations.Source: §5 · p. 26
• Across 32–2,048 cores, IllinoisGRMHD delivers 1.7–1.8× standard GRHydro's per-core throughput, and is about 16% faster than GRHydro-experimental at 2,048 cores.At 32 cores, IllinoisGRMHD and GRHydro-experimental are within measurement error.Source: §5 · p. 26, Figure 5 · p. 27
• For a typical production run, IllinoisGRMHD's OpenMP/MPI hybrid is slightly more than 40% faster than pure MPI.OpenMP reduces network load by allowing fewer MPI processes and larger grids per process.Source: §5 · p. 27

6. Conclusions and Future Work

• IllinoisGRMHD was rewritten to improve user-friendliness and modularity, with documentation, basic programming requirements, and low-level CFD routines separated into an extensible function library.These changes were intended to facilitate community adoption.Source: §6 · p. 29
• IllinoisGRMHD reproduces OrigGRMHD results to roundoff-level precision for magnetized neutron stars and discontinuous random initial data.It was designed as a drop-in replacement for the original code.Source: §6 · p. 29
• IllinoisGRMHD and GRHydro both show approximate second-order convergence, but differ in constraint violations and central-density drift during unmagnetized TOV-star evolution.GRHydro has about 8% higher Hamiltonian constraint violations, while its absolute central-density drift is lower; the drift rates match after initial settling over a few dynamical timescales.Source: §6 · p. 29
• IllinoisGRMHD is 1.7–1.8 times faster than standard GRHydro and about 16% faster than GRHydro-experimental at 1,024 to 2,048 cores.The standard-GRHydro comparison is for production-size, AMR-enabled GRMHD runs on Stampede; the experimental version matches IllinoisGRMHD at best at small core counts.Source: §6 · pp. 29–30
• IllinoisGRMHD enforces the no-monopole constraint on multi-scale AMR flows and can stably evolve flows into black-hole horizons over long timescales without excision.The authors state that GRHydro cannot guarantee the no-monopole constraint at AMR boundaries and requires excision for black-hole flows.Source: §6 · pp. 29–30
• IllinoisGRMHD outperforms OrigGRMHD by a factor of 1.3–1.6.The authors state that adopting IllinoisGRMHD would increase computational resources available for dynamical-spacetime GRMHD runs.Source: §6 · p. 30
• Future work includes comparing algorithmic implementations and extending IllinoisGRMHD with additional solvers, equations of state, and symmetry conditions.The current version supports only single gamma-law equations of state and either no symmetry or bitant symmetry across the xy-plane.Source: §6 · p. 30

Appendix A. Physical and Computational Setup for Unmagnetized TOV Star Convergence Tests

• The initial data describe a cold, degenerate TOV neutron star with central density 0.129285309, K = 1, Γ = 2, and compaction MNS/RNS ≈0.1467.The model uses P = Kρ^Γ and is generated with the Einstein Toolkit's built-in TOVSolver thorn.Source: Appendix A · p. 31
• IllinoisGRMHD and GRHydro are coupled to McLachlan BSSN and use matched evolution algorithms, with HLLE and PPM selected for GRHydro.The codes differ in implementation details, including reconstructing internal energy in GRHydro rather than pressure.Source: Appendix A · p. 31
• The tests use four nested refinement cubes and resolve the neutron-star diameter with approximately 51, 64, or 80 gridpoints in low-, medium-, and high-resolution runs.The finest-level grid spacing is ∆x = {0.03906, 0.03125, 0.025}RNS; the atmosphere floor is 10^-7 times the initial central density.Source: Appendix A · p. 31

Appendix B. Physical and Computational Setup for Magnetized TOV Star Code Validation Tests

• The magnetized initial data use the same stable neutron star as Appendix A, seeded at t = 0 with a weak, purely poloidal magnetic field.The maximum initial magnetic-to-gas pressure ratio β^-1 is 0.83×10^-3, with Pcut set to 4% of maximum pressure.Source: Appendix B · p. 32
• IllinoisGRMHD and OrigGRMHD are coupled to the Illinois group's closed-source BSSN dynamical-spacetime module for these tests.The initial star density, pressure, and metric profiles come from the code of and agree to machine precision with the unmagnetized-test TOVSolver data.Source: Appendix B · p. 32
• The AMR setup resolves the neutron-star diameter with 26 gridpoints, intentionally using low resolution and a close outer boundary to magnify truncation-error and boundary-condition effects.The finest-grid spacing is ∆x ≈ 0.078RNS; the atmosphere floor is 10^-7 times the initial central density.Source: Appendix B · p. 32

Appendix C. Random Initial Data

• Each gridpoint receives a double-precision random number ξ ∈ (−0.1, 0.1) from a coordinate-based seed, making data reproducible across processor decompositions on identical grids.The module uses the gridpoint coordinate to generate the seed.Source: Appendix C · p. 32
• The initial prescription sets spacetime and GRMHD quantities from a stochastic perturbation atop a smoothly oscillating function, with remaining quantities computed from the prescribed values.The setup enforces a Lorentzian spacetime metric, a positive-definite three-metric, and the BSSN unit-determinant constraint by rescaling the spatial metric components.Source: Appendix C · p. 33
• The random-data tests use one refinement level and disable spacetime evolution, along with prolongation and restriction operations on spacetime variables.The coarser grid has spacing 1.0 and half-side-length 10; the finer grid has spacing 0.5 and half-side-length 2.Source: Appendix C · p. 33
• During evolution, the Lorentz factor limit W = 10 is exceeded and enforced about 400 times despite initially independent random values at neighboring gridpoints.The initial magnetic-to-gas pressure ratio ranges from approximately 10^-3 to about 10, while density fluctuations span about four orders of magnitude.Source: Appendix C · pp. 33–34