GRaM-X: Dynamical-Spacetime GRMHD on GPUs

• GRaM-X is a GPU-accelerated GRMHD code that couples dynamical spacetime evolution with adaptive mesh refinement within the Einstein Toolkit.
• It employs a mixed centering strategy with high-order TVD/WENO reconstructions and the HLLE Riemann solver to accurately solve coupled Einstein and MHD equations.
• Validated through shock tubes, cylindrical explosions, and 3D TOV star benchmarks, GRaM-X achieves 40%-50% weak-scaling efficiency on exascale GPU systems.

GRaM-X, expanded as General Relativistic accelerated Magnetohydrodynamics on AMReX, is a GPU-accelerated dynamical-spacetime general relativistic magnetohydrodynamics (GRMHD) code developed within the Einstein Toolkit to target GPU-based exascale systems. It extends the Einstein Toolkit GRMHD capability by combining dynamical spacetime evolution, 3D adaptive mesh refinement (AMR) on GPUs via CarpetX and AMReX, the Z4c formalism for general relativity, and the Valencia formulation for GRMHD. The code supports both analytic and tabulated equations of state, implements TVD and WENO reconstruction together with the HLLE Riemann solver, and is validated on standard static- and dynamical-spacetime benchmarks. The reported weak-scaling result is approximately $40\%-50\%$ efficiency on 2304 Summit nodes, corresponding to 13824 NVIDIA V100 GPUs, relative to single-node performance (Shankar et al., 2022).

1. Definition and position within the Einstein Toolkit

GRaM-X is presented as a new GPU-accelerated dynamical-spacetime GRMHD code for the Einstein Toolkit. Its stated purpose is to enable dynamical-spacetime GRMHD simulations on GPU-based exascale systems, a setting motivated by applications such as compact-object mergers, magnetized neutron stars, and core-collapse supernovae. Within the Einstein Toolkit software stack, it is built as a Cactus-based capability layered on the new CarpetX driver, which itself leverages AMReX as the underlying AMR framework (Shankar et al., 2022).

The code is described as a complete redesign of GRHydro for GPU-era computing rather than as an incremental port of earlier CPU-centric infrastructure. The distinction emphasized in the literature is not that GRaM-X is the first GPU GRMHD code in general, but that it combines dynamical spacetime, GRMHD, 3D AMR, and Einstein Toolkit integration in a GPU-oriented implementation. This places it at the intersection of numerical relativity, relativistic astrophysics, and exascale heterogeneous computing (Shankar et al., 2022).

In architectural terms, the relevant layers are Cactus as the multiphysics framework, CarpetX as the AMR and driver layer, and AMReX as the low-level block-structured AMR and heterogeneous runtime substrate. That organization makes GRaM-X part of the broader Einstein Toolkit ecosystem rather than a standalone monolithic code. A plausible implication is that its design prioritizes interoperability with existing Einstein Toolkit components alongside accelerator performance.

2. Software architecture, AMR model, and heterogeneous execution

GRaM-X uses block-structured Berger–Oliger AMR through CarpetX. A coarse rectangular mesh is overlaid by rectangular refined grids organized into refinement levels with factor-of-two spacing ratios, and proper nesting is required so that level $L+1$ lies entirely within level $L$. Refined grids carry ghost zones, which are filled either by copying from neighboring blocks or MPI ranks or by prolongation from the next coarser level; restriction updates coarse data where fine grids overlay coarse regions (Shankar et al., 2022).

A notable implementation choice is the mixed centering of fields. Vertex-centered variables are used for spacetime geometry, while conserved fluid quantities are cell-centered. This is coupled to higher-order generic templated stencils in CarpetX for inter-mesh transport, so that both vertex-centered geometry and cell-centered conservative data are handled consistently across AMR boundaries. For hydrodynamics, face-centered fluxes are synchronized conservatively across blocks and coarse-fine interfaces. Because the code currently uses a global time step across all AMR levels, refluxing is unnecessary (Shankar et al., 2022).

The parallelization strategy is explicitly three-level: OpenMP shared-memory parallelism on CPUs, GPU kernel parallelism through AMReX accelerator facilities, and MPI across nodes. AMReX GPU kernel launches are used via lambda functions, and the software stack is stated to support CUDA, HIP/ROC, and DPC++/SYCL as applicable. On CPUs, CarpetX and GRaM-X use OpenMP and NSIMD for vectorization. The design goal on GPUs is that all data remain resident in accelerator memory except during I/O, with even simple reductions performed on device because device-host transfers are too expensive (Shankar et al., 2022).

The mesh is decomposed into rectangular blocks, which are the units of scheduling and distribution. A representative block size given in the implementation discussion is $8\times 8\times 8$ cells. On CPUs, tiles are chosen with cache behavior in mind; for example, $1000\times 4\times 4$ is cited as a typical efficient OpenMP tile shape. A central exascale-oriented design decision is the absence, at this stage, of time subcycling: all refinement levels use the same time step. The stated rationale is that this simplifies synchronization and improves parallel scalability because all levels can evolve simultaneously rather than serializing level advances (Shankar et al., 2022).

3. Physical formulation: Z4c spacetime evolution and Valencia GRMHD

GRaM-X evolves both spacetime and magnetized relativistic fluid. For spacetime, it uses the Z4c formalism for Einstein’s equations, together with $1+\log$ foliation for the lapse and a $\Gamma$-driver condition for the shift. The parameters reported for the gauge and damping system are

$$\kappa_1 = 0.02, \qquad \kappa_2 = 0.0, \qquad \mu_L=\frac{2}{\alpha}, \qquad \mu_S=1, \qquad \eta=2.$$

Although Z4c introduces additional dynamical fields, other Einstein Toolkit thorns interface through the standard ADM variables $\gamma_{ij}$, $K_{ij}$, $L+1$0, $L+1$1, $L+1$2, and $L+1$3 (Shankar et al., 2022).

For the matter sector, GRaM-X adopts the Valencia formulation of ideal GRMHD. The governing conservation laws are written as

$L+1$4

with ideal MHD imposed through

$L+1$5

The stress-energy tensor is decomposed into hydrodynamic and electromagnetic parts, and the total system is expressed in flux-conservative form. The primitive variables are

$L+1$6

while the conserved variables are

$L+1$7

with

$L+1$8

$L+1$9

$L$0

$L$1

Here $L$2, and $L$3 (Shankar et al., 2022).

A distinctive numerical consequence of the mixed centering strategy is that $L$4 is stored at vertices whereas fluid variables are cell-centered cell averages. To construct $L$5, fluid variables must therefore be interpolated from cell centers to vertices using fourth-order symmetric 3D interpolation. This interpolation step is not an incidental detail; it materially affects the accuracy of the coupled spacetime-matter evolution, particularly in dynamical-spacetime tests (Shankar et al., 2022).

4. Numerical methods, equations of state, and magnetic-field treatment

GRaM-X uses the Method of Lines for all evolved fields. A typical right-hand-side evaluation proceeds by recovering primitive variables from conserved variables, reconstructing primitive states to faces, solving Riemann problems, and then forming flux divergences and source terms. The primary conservative-to-primitive inversion is a 3D Newton–Raphson method solving for $L$6. A Newman–Hamlin method serves as a fallback, and if both fail, bisection in temperature is used as a guaranteed-convergence but expensive final fallback. The rationale given for preferring the 3D Newton–Raphson method is that Newman’s method requires about 20 times more equation-of-state calls, which is particularly costly for tabulated EOS (Shankar et al., 2022).

The reconstruction options implemented are TVD with limiter choices including minmod, and WENO5. TVD is described as second-order accurate in smooth monotonic regions and first order near shocks or extrema, whereas WENO is fifth-order accurate in smooth monotonic flows and intended for production simulations. The Riemann solver implemented in the reported version is HLLE. Time integration uses RK4 in all validation problems except benchmark setup A, for which RK2 is employed (Shankar et al., 2022).

GRaM-X supports both analytic and tabulated equations of state. The analytic options explicitly mentioned are the polytropic EOS,

$L$7

and the ideal-gas or $L$8-law EOS,

$L$9

For realistic nuclear matter, the tabulated EOS support includes tables containing 19 fluid variables such as pressure $8\times 8\times 8$0, specific enthalpy $8\times 8\times 8$1, and speed of sound $8\times 8\times 8$2 as functions of $8\times 8\times 8$3. These tables are loaded into unified memory so that both CPU and GPU can access them (Shankar et al., 2022).

The treatment of the magnetic divergence constraint is problem dependent in the published study. The 1D shock-tube tests are run without constrained transport, while the 2D cylindrical explosion, 2D magnetic rotor, and 3D TOV star tests use constrained transport to preserve the divergence-free magnetic-field condition. The paper notes the advantages of placing the magnetic vector potential $8\times 8\times 8$4 at edges for maintaining $8\times 8\times 8$5 exactly up to floating-point precision, but it does not give a full formal description of the constrained-transport algorithm used in GRaM-X (Shankar et al., 2022).

5. Verification, benchmark problems, and measured performance

The published validation suite spans static spacetime, special-relativistic MHD, and full dynamical-spacetime evolution. On static backgrounds, GRaM-X is tested with the five standard Balsara shock tubes, a 2D cylindrical explosion, and a 2D magnetic rotor. On dynamical spacetimes, the principal benchmark is the oscillation of a 3D magnetized TOV star (Shankar et al., 2022).

Test class	Configuration	Reported outcome
1D MHD shock tubes	Balsara 1–5	Very good agreement with analytic solutions
2D cylindrical explosion	Flat spacetime, constrained transport	Very good agreement with GRHydro and Komissarov
2D magnetic rotor	Flat spacetime, constrained transport	Matches prior published rotor solutions very well
3D TOV star oscillations	Dynamical spacetime, Z4c	Convergence of central density and good mode-frequency agreement

In the Balsara shock tubes, the domain is $8\times 8\times 8$6 with 1600 points and $8\times 8\times 8$7, using RK4, CFL $8\times 8\times 8$8, TVD reconstruction with minmod, HLLE, and von Neumann boundaries. The code is reported to reproduce the elementary-wave structure accurately, including the more challenging features of the relativistic cases. The Balsara 3 test is identified as the most difficult, and the reported undershoot with TVD at $8\times 8\times 8$9 improves substantially when WENO is used (Shankar et al., 2022).

The 2D cylindrical explosion is evolved on $1000\times 4\times 4$0 with a $1000\times 4\times 4$1 grid, $1000\times 4\times 4$2, a uniform field $1000\times 4\times 4$3, constrained transport, TVD plus minmod, HLLE, RK4, and CFL $1000\times 4\times 4$4, up to $1000\times 4\times 4$5. The paper reports very good agreement with GRHydro and Komissarov, and attributes minor discrepancies to different centering conventions. The 2D magnetic rotor, on $1000\times 4\times 4$6 with a $1000\times 4\times 4$7 grid and an initial maximum edge velocity $1000\times 4\times 4$8, is likewise reported to match published reference solutions very well (Shankar et al., 2022).

The 3D TOV-star test is the key dynamical-spacetime validation. The equilibrium star is constructed from the Tolman–Oppenheimer–Volkoff solution with $1000\times 4\times 4$9, $1+\log$0, and central density $1+\log$1, yielding $1+\log$2 and $1+\log$3. A poloidal magnetic field is initialized through a vector potential, and four AMR resolutions are used with finest spacings $1+\log$4, $1+\log$5, $1+\log$6, and $1+\log$7. The reported convergence order of the central density is approximately $1+\log$8, and the Fourier spectrum agrees well with perturbation theory for the $1+\log$9, $\Gamma$0, and $\Gamma$1 modes, while $\Gamma$2 is slightly shifted at the highest tested resolution (Shankar et al., 2022).

Performance measurements are reported primarily on OLCF’s Summit. Single-GPU throughput on Summit V100 hardware reaches $\Gamma$3 to $\Gamma$4 zone-cycles/s for static-spacetime setup A, $\Gamma$5 to $\Gamma$6 zone-cycles/s for static-spacetime setup B with tabulated EOS and RK4, and $\Gamma$7 zone-cycles/s for dynamical-spacetime setup C. The headline weak-scaling result is approximately $\Gamma$8 efficiency on 2304 nodes, i.e. 13824 NVIDIA V100 GPUs, relative to one node. For setup C, the reported efficiency remains roughly flat from 36 nodes to 2304 nodes, which the paper interprets as stable large-scale scaling once inter-node communication becomes dominant. The principal bottleneck identified is ghost-zone communication bandwidth rather than latency (Shankar et al., 2022).

6. Later extensions, comparison points, and documented limitations

A later paper extends GRaM-X from ideal GRMHD to full resistive GRMHD in dynamical spacetime. The resistive module evolves the full resistive GRMHD equations using Z4c, treats stiff source terms with a second-order IMEX-SSP2 $\Gamma$9 scheme, employs a one-dimensional conservative-to-primitive recovery, and uses HLLE with TVD and WENO reconstruction. Its validation suite includes 1D shock tubes, current sheets, Alfvén waves, 2D cylindrical explosions, and 3D TOV stars. The stated scientific motivation is that ideal GRMHD neglects resistive effects needed for magnetic reconnection, dissipation, magnetospheric dynamics, and transitions between high-conductivity interior plasma and low-conductivity or electrovacuum exterior regions (Azizi et al., 21 Oct 2025).

Within the same Einstein Toolkit, CarpetX, and AMReX design space, AsterX is presented as a closely related GPU-accelerated GRMHD code for dynamical spacetimes. It also uses Valencia GRMHD and Z4c, but is distinguished by a staggered vector-potential magnetic-field treatment, generalized Lorenz gauge, and subcycling in time. Its reported weak-scaling efficiency is about $$\kappa_1 = 0.02, \qquad \kappa_2 = 0.0, \qquad \mu_L=\frac{2}{\alpha}, \qquad \mu_S=1, \qquad \eta=2.$$0 on 4096 Frontier nodes compared to 8-node performance, and the paper explicitly positions AsterX as operating in a similar algorithmic and performance regime to GRaM-X (Kalinani et al., 2024). This suggests that GRaM-X has become part of a broader family of Einstein Toolkit–native exascale GRMHD codes rather than an isolated development.

The documented limitations of the 2022 GRaM-X release are specific. No AMR performance benchmarks are included in that paper; all reported scaling measurements are unigrid. Time subcycling is not yet implemented. The spacetime solver is not yet fully GPU-optimized and exhibits register spills. HLLE is the only Riemann solver presently implemented in that version, with Roe and Marquina identified as future additions. The mixed centering of hydro and spacetime variables makes $$\kappa_1 = 0.02, \qquad \kappa_2 = 0.0, \qquad \mu_L=\frac{2}{\alpha}, \qquad \mu_S=1, \qquad \eta=2.$$1 approximate and increases numerical error in coupled evolutions such as TOV oscillations. The paper also does not mention neutrino transport in the current version, although testing and implementation of an M1 moment-based neutrino transport scheme is explicitly reported as future work (Shankar et al., 2022).

In that sense, GRaM-X is best understood as an exascale-oriented Einstein Toolkit GRMHD capability whose defining features are GPU execution, dynamical spacetime, 3D AMR compatibility through CarpetX and AMReX, and support for both analytic and tabulated EOS. Its later resistive extension broadens the mathematical model from ideal to resistive GRMHD, while related CarpetX-based codes such as AsterX provide an external comparison point for formulation choices, AMR strategy, and large-scale scaling behavior (Azizi et al., 21 Oct 2025).