⚛️ general relativity

GRACE: An Open-Source Framework for GPU-Accelerated Numerical Relativity

This paper introduces GRACE, a new open-source, GPU-accelerated numerical relativity framework built on Kokkos and p4est that couples ideal GRMHD with the Z4c formulation of Einstein's equations, demonstrating high performance and accuracy through extensive validation tests and binary neutron star merger simulations.

Original authors: Carlo Musolino, Christian Ecker, Konrad Topolski, Marie Cassing, Keneth Miler, Harry Ho-Yin Ng, Khalil Pierre, Elias R. Most, Luciano Rezzolla

Published 2026-07-14

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Carlo Musolino, Christian Ecker, Konrad Topolski, Marie Cassing, Keneth Miler, Harry Ho-Yin Ng, Khalil Pierre, Elias R. Most, Luciano Rezzolla

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the universe as a giant, chaotic dance floor where black holes spin, neutron stars collide, and magnetic fields twist like spaghetti. For decades, scientists have tried to simulate this dance on computers to understand how gravity and matter behave in the most extreme conditions. But there's a problem: the computers they use are often too slow, or the software they wrote is built for old-school processors that just can't keep up with the speed of modern graphics cards (GPUs).

Enter GRACE (General Relativistic Astrophysics Code for Exascale). Think of GRACE as a brand-new, high-speed race car built from scratch specifically to run on the latest, most powerful computer hardware available today. It's not just a patch-up job; it's a complete rebuild designed to be "portable," meaning it can drive just as fast on a computer with a graphics card from NVIDIA, AMD, or even Intel, without needing a new engine for each one.

The Big Idea: A Digital Time Machine

The main finding of this paper is that GRACE works. The authors built this framework to solve the complex equations that describe how space and time warp (Einstein's equations) while simultaneously tracking how super-dense, magnetized fluids (like those inside neutron stars) move. They tested it against a whole menu of standard "driving tests" to make sure it doesn't crash.

They simulated everything from simple shockwaves in a flat, empty room to the violent crash of two neutron stars. In every case, GRACE produced results that matched other trusted codes and theoretical predictions. For instance, when they simulated two neutron stars spiraling into each other, the "gravitational waves" (the ripples in space-time) GRACE predicted looked almost identical to those from a different, well-known code called FIL. This suggests that GRACE is a reliable new tool for exploring the universe.

What It's NOT (And What It Rules Out)

It's important to know what GRACE doesn't do, or at least what it hasn't proven yet. The paper explicitly states that while they tested the code with magnetic fields, the simulations of neutron star mergers were not fine-grained enough to perfectly capture the tiny, chaotic "dynamo" effects that happen right at the moment of collision.

Think of it like trying to film a hurricane with a camera that has a slightly blurry lens. You can see the storm is there, you can see it spinning, and you can see the wind is strong, but you can't quite make out the individual raindrops or the tiny swirls inside the clouds. The authors are clear: their current simulations show that magnetic fields get stronger after the crash, but the exact amount of amplification isn't "converged" yet. They need even more powerful computers to see those tiny details clearly. They also ruled out the idea that their code is a magic bullet for everything right now; it doesn't yet include the complex physics of neutrinos (ghostly particles) or the heat exchange of weak nuclear forces, though they plan to add those later.

How Sure Are They?

The authors are very confident about the code's ability to run fast and handle the math correctly, but they are careful to call their results "simulations" rather than absolute truths about the universe.

The Speed Test: They measured how fast GRACE runs on different machines. On a super-fast AMD chip, it can update about 24.8 million tiny grid cells per second. On an NVIDIA A100, it does about 10.5 million. On a standard Intel CPU, it's much slower at 0.33 million. This proves the code is built to exploit the speed of modern graphics cards.
The Accuracy Test: When they compared GRACE to the FIL code for a binary neutron star merger, the difference in the timing of the crash was less than 0.01 milliseconds. The difference in the phase of the gravitational waves was tiny, consistent with the natural "pixelation" errors of any computer simulation. This suggests the two codes are agreeing on the physics, not just getting lucky.
The Limits: They admit that for the most complex test (the magnetized merger), the magnetic field amplification they saw (about 10 times stronger) is likely an underestimate because their grid wasn't fine enough. They suggest that to see the real, full effect, they would need a grid spacing of about 10–20 meters, which is much finer than the 185 meters they used in their best run.

The "Forest" and the "Ghost Zones"

To make this work, the team used some clever tricks. Imagine the simulation grid as a forest of trees. Where the action is happening (like near the colliding stars), the trees are tiny and dense to catch every detail. Far away, where nothing is happening, the trees are huge and sparse to save space. This is called "Adaptive Mesh Refinement."

But here's the tricky part: when a tiny tree meets a big tree, the data has to pass between them without getting lost. The authors built a system of "ghost zones"—imaginary buffers that let the different-sized grids talk to each other smoothly. They also used a special method called "constrained transport" to make sure the magnetic field lines never break or get tangled in a way that violates the laws of physics (specifically, that magnetic monopoles don't exist).

The Verdict

GRACE is a successful, open-source framework that proves you can run complex, universe-crushing simulations on modern, heterogeneous computer clusters. It's not a finished product that solves every mystery of the cosmos; rather, it's a highly capable new engine that has passed its driving test. The authors have shown it can handle the math, scale up to hundreds of processors, and produce results that agree with other top-tier codes.

They are now ready to take this engine for a longer drive, adding in the missing pieces like neutrino physics and better magnetic field resolution, but for now, they have successfully built the car and proven it can drive on the track.

Technical Summary: GRACE – An Open-Source Framework for GPU-Accelerated Numerical Relativity

Problem and Motivation

Numerical relativity (NR) has matured to a point where it provides quantitative predictions for relativistic astrophysical systems, such as binary neutron star (BNS) mergers and black hole dynamics. However, the field increasingly relies on high-performance computing (HPC) facilities featuring heterogeneous architectures, where Graphics Processing Units (GPUs) supplement traditional CPUs. While many NR codes have been adapted for GPUs, this transition presents significant challenges. Algorithms well-suited for CPUs often perform poorly on GPU architectures due to differences in memory access patterns and parallelism requirements. Furthermore, existing frameworks often lack a unified approach to portability across diverse hardware (e.g., NVIDIA, AMD, Intel) while maintaining high efficiency. There is a need for a framework built from the ground up to leverage modern heterogeneous HPC systems efficiently, utilizing open-source libraries to ensure maintainability and reproducibility.

Methodology

The authors present GRACE (General Relativistic Astrophysics Code for Exascale), a new NR framework designed for efficiency and portability on heterogeneous platforms.

Core Architecture and Portability

Kokkos: GRACE is built exclusively on open-source libraries, utilizing Kokkos as a minimal abstraction layer. This allows the code to target multiple backends (CUDA, HIP, OpenMP) from a single source tree, ensuring performance portability across CPU and GPU architectures.
Grid Management: The numerical grid is constructed as a forest of octrees using the p4est library. p4est runs exclusively on the host (CPU) to manage connectivity, neighbor lookups, tree balancing, and load balancing. All data-bearing operations on grid cells are executed on the device (GPU) using custom-written kernels.
Adaptive Mesh Refinement (AMR): GRACE supports both Fixed Mesh Refinement (FMR) and AMR. The grid maintains a 2:1 balance between neighbors. Refinement and coarsening are handled via "prolongation" and "restriction" operations performed on the GPU.

Physical Formulations

Einstein Equations: The code evolves the Einstein equations using the Z4c formulation, which includes constraint-damping terms to control numerical errors. The system uses the "puncture gauge" with 1+log slicing for the lapse and Gamma-driver equations for the shift.
GRMHD: The framework solves the equations of Ideal General-Relativistic Magnetohydrodynamics (GRMHD) for plasmas with infinite conductivity.
- Conservation Laws: The system uses a conservative formulation with high-resolution shock-capturing (HRSC) methods.
- Riemann Solvers: It employs approximate Riemann solvers, specifically Harten–Lax–van Leer–Einfeldt (HLLE) or Local Lax–Friedrichs (LLF).
- Reconstruction: The code utilizes a hierarchy of reconstruction operators, with fifth-order Weighted Essentially Non-Oscillatory (WENO-Z) schemes as the default.
- Constrained Transport (CT): To maintain the divergence-free condition of the magnetic field ( $\nabla \cdot \mathbf{B} = 0$ ) to machine precision, GRACE employs a constrained transport approach. The magnetic field is stored at cell faces, and electromotive forces (EMFs) are computed at cell edges. The implementation uses the "CT-contact" scheme, which averages EMFs from surrounding faces and corrects them with upwinded derivative terms.
Equation of State (EOS): GRACE supports four EOS families: ideal gas, piecewise polytropic, hybrid (cold + thermal), and fully tabulated 3D EOSs (CompOSE format).

Numerical Implementation Details

Time Integration: The system uses explicit Runge-Kutta (RK) schemes, with a four-stage, fourth-order low-storage RK method (3S*) as the default.
Conservative-to-Primitive (C2P): The code implements a robust inversion scheme (based on [92]) to recover primitive variables from conserved ones, with fallback mechanisms (entropy-based inversion) for ill-conditioned states.
Boundary Conditions: Outer boundaries use Sommerfeld-type radiative conditions. Ghost zones are updated using a task-based execution model represented as a Directed Acyclic Graph (DAG) to optimize memory-bound and compute-bound operations on GPUs.
Flux Correction: An a posteriori first-order flux correction (FOFC) is implemented to prevent unphysical states in regions of strong gradients or low density.

Key Contributions

Framework Development: The creation of GRACE, a novel NR framework built from scratch for heterogeneous computing, leveraging Kokkos for portability and p4est for AMR.
GPU Optimization: Implementation of a task-based ghost-zone update system and custom kernels that allow efficient execution on GPUs, including specific handling of constrained transport and magnetic field divergence preservation across refinement boundaries.
Validation Suite: A comprehensive validation of the code against a wide range of standard tests, including:
- Magnetized shock tubes and magnetic rotors in flat spacetime.
- Bondi accretion onto a Schwarzschild black hole.
- Neutron-star oscillations in both fixed (Cowling) and dynamical spacetimes.
- Ringdown of perturbed spinning punctures.
- Binary black hole inspiral and merger.
Binary Neutron Star Simulations: Evolution of two distinct BNS merger scenarios:
- An equal-mass, unmagnetized system with an ideal-gas EOS.
- An unequal-mass ( $q=0.8$ ), magnetized system with a finite-temperature, tabulated SFHo EOS.
Cross-Code Comparison: A detailed comparison of the magnetized BNS merger results with the FIL (Frankfurt/IllinoisGRMHD) code, demonstrating agreement in inspiral dynamics within the limits of discretization errors.

Results

Accuracy: GRACE successfully reproduces exact solutions for shock tubes and quasi-analytic solutions for Bondi accretion. It accurately captures neutron-star oscillation spectra (fundamental and overtone modes) and binary black hole ringdown waveforms.
Convergence: The code exhibits expected convergence orders. For binary black hole mergers, spatial discretization errors (sixth-order) dominate the error budget. For BNS inspirals, the convergence order is observed to be approximately $p \approx 2.8$ , lying between the second-order accuracy of the finite-volume matter sector and the fourth-order accuracy of the spacetime sector.
Conservation: The code conserves rest-mass to the level of accumulated round-off error ( $\lesssim 10^{-10}$ ) during the inspiral phase. The constrained transport scheme maintains the discrete divergence-free magnetic field constraint to machine precision ( $\sim 10^{-19}$ ) across refinement boundaries and MPI partitions.
Performance:
- Throughput: On an AMD MI300A APU, GRACE achieves $\sim 2.48 \times 10^7$ zone-cycles per second (fixed spacetime) and $\sim 1.68 \times 10^7$ (coupled GRMHD-Z4c). On an NVIDIA A100, these values are $\sim 1.05 \times 10^7$ and $\sim 6.7 \times 10^6$ , respectively.
- Scaling: Strong scaling tests show efficiencies of $\sim 44\%$ (uniform grid) and $\sim 54\%$ (FMR grid) when increasing device counts by factors of 16 and 8. Weak scaling tests demonstrate efficiencies of $\sim 90\%$ (uniform) and $\sim 98\%$ (FMR) up to 128 and 256 devices, respectively.
- Production Runs: Simulations of BNS mergers on the Viper cluster (AMD MI300A) achieved sustained evolution rates of $\sim 300$ $M_\odot$ /h, with computational costs scaling roughly as $h^{-4}$ , consistent with theoretical expectations for 3+1 dimensional evolution.

Significance and Claims

The paper claims that GRACE represents a significant step forward in making numerical relativity accessible and efficient on modern heterogeneous HPC platforms. By being built entirely on open-source libraries (Kokkos, p4est), it ensures long-term maintainability and portability. The authors emphasize that GRACE is not merely a port of existing codes but a new framework designed to overcome the specific algorithmic challenges of GPU architectures, particularly regarding memory access and parallelism.

The validation against the FIL code and the successful reproduction of complex phenomena (such as magnetized BNS mergers with tabulated EOSs) serve as evidence of the code's correctness and robustness. The authors modestly note that while the current resolution is sufficient for validating the code's dynamics and convergence, capturing small-scale dynamo effects in post-merger magnetic field amplification requires resolutions beyond the scope of this validation paper. The release of GRACE alongside GRACEpy (a Python environment for post-processing) is presented as a contribution to the community's ability to perform and analyze high-fidelity NR simulations.