cuHARM: GPU-Accelerated GRMHD Code

• cuHARM is a GPU-accelerated GRMHD code built in CUDA-C with OpenMP for multi-GPU execution, enabling efficient 3D simulations of black-hole accretion and jet formation.
• It employs a conservative finite-volume formulation with flux-CT and HARM-style reconstruction to maintain stability and accuracy in Kerr-spacetime simulations.
• Extensions include systematic MAD spin surveys and a fully covariant, GPU-accelerated radiative module to capture anisotropic radiation transport in challenging regimes.

Searching arXiv for recent cuHARM papers and related code/application papers. cuHARM is a GPU-accelerated general-relativistic magnetohydrodynamics code in the HARM family, written in CUDA-C and using OpenMP for multi-GPU execution on a single workstation. In its original form, it was introduced as a GR-MHD solver for Kerr-spacetime accretion and jet problems, with emphasis on 3D simulations of ADAF disks and SANE states (Bégué et al., 2022). Subsequent work employed cuHARM for broad spin surveys of magnetically arrested disks, characterizing angular-momentum transport, jet morphology, magnetic-pressure structure, and horizon diagnostics across retrograde, non-rotating, and prograde black holes (Zhang et al., 2023). A later extension added a fully covariant, GPU-accelerated radiation module that directly solves the time-dependent radiative transfer equation on a geodesic angular grid and couples the resulting radiation field back to the GRMHD evolution, thereby promoting cuHARM into a general relativistic radiation magnetohydrodynamics platform (Wallace et al., 21 Aug 2025).

1. Code identity and lineage

cuHARM was introduced as “a new GPU-accelerated general-relativistic magnetohydrodynamic (GR-MHD) code” based on HARM and developed with “many major modifications and improvements” relative to a literal port (Bégué et al., 2022). In later application work it is described as a “GPU-accelerated GRMHD code” developed “on HARMPI,” using “openMP and cuda,” and “thoroughly optimized for maximal harness of the power available in NVidia GPUs” (Zhang et al., 2023).

The defining architectural claim is that the expensive numerical workload is executed on GPUs rather than being CPU-dominated. The 2022 code paper states that cuHARM is written in CUDA-C and uses OpenMP to manage multi-GPU execution on a single node, with CPU involvement mainly for output and GPU-to-GPU data transfer (Bégué et al., 2022). The 2023 MAD study reiterates that “all calculations of cuHARM are accelerated by GPU (only the data transfer and exports are powered by CPU)” (Zhang et al., 2023). This implementation choice was motivated by the computational cost of global 3D GRMHD simulations of black-hole accretion and jet formation.

The code’s scientific niche is the numerical study of accretion disks, jet launching, and magnetized plasma dynamics in fixed Kerr geometry. The original application targeted ADAF disks in the SANE regime (Bégué et al., 2022); the later spin survey focused on MAD accretion flows across negative, zero, and positive black-hole spin (Zhang et al., 2023). The radiative extension broadened the intended scope to radiatively important black-hole accretion states, including slim disks and super-Eddington or moderately radiative systems (Wallace et al., 21 Aug 2025).

2. GRMHD formulation and numerical scheme

In its non-radiative form, cuHARM solves the equations of ideal, non-resistive GRMHD in fixed spacetime. The governing equations are mass conservation,

$\nabla_\mu \left ( \rho u^\mu \right ) = 0,$

energy-momentum conservation,

$\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$

and the homogeneous Maxwell equation,

$\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$

(Bégué et al., 2022). Under ideal MHD,

$u_\mu F^{\mu \nu} = 0,$

with magnetic four-vector

$b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$

and stress-energy tensor

$$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$

(Bégué et al., 2022).

Following HARM, cuHARM adopts a conservative finite-volume formulation. The evolved conserved variables are

$$U = \sqrt{-g}(\rho u^t, T^t_{~t}, T^{t}_{~i}, B^i),$$

while primitive variables are

$P = (\rho, u, \tilde u^i, B^i)$

(Bégué et al., 2022). Instead of evolving $u^i$ directly, the code uses modified velocities

$$\tilde u^i \equiv u^i + \frac{\Gamma \beta^i}{\alpha},$$

a HARM-style choice intended to improve numerical stability and maintain timelike reconstructed four-velocities (Bégué et al., 2022).

The discretization is a Godunov-type finite-volume method on a structured 3D grid. Primitive reconstruction uses piecewise linear reconstruction with the monotonized central limiter (Bégué et al., 2022). The default numerical flux in the 2022 science runs is Lax–Friedrichs,

$\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$0

although HLL fluxes were also implemented but not used in that study because Lax–Friedrichs produced fewer failures for the chosen setups (Bégué et al., 2022). The characteristic wave-speed estimate follows the HARM approximation rather than the full quartic relativistic MHD dispersion relation (Bégué et al., 2022).

For divergence control, cuHARM uses constrained transport, specifically the flux-CT method of Tóth as adopted in HARM (Bégué et al., 2022). The paper explicitly notes the known drawbacks of flux-CT—large stencil, improper 1D limit, and lack of upwinding—while retaining it for simplicity (Bégué et al., 2022). Time integration is second order, using a predictor-corrector or midpoint update with CFL parameter $\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$1 (Bégué et al., 2022).

Primitive-variable recovery is handled through a hierarchy of inversion strategies: a 2D inversion solving for $\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$2 and $\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$3, a 1D inversion solving for $\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$4, and an entropy-based 1D inversion used either as fallback or directly in highly magnetized regions (Bégué et al., 2022). The code evolves entropy,

$\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$5

to stabilize problematic cells and provide backup when energy-based primitive recovery fails (Bégué et al., 2022).

3. Coordinates, boundaries, floors, and execution model

cuHARM is designed for Kerr black-hole problems in horizon-penetrating coordinates. The 2022 code paper uses the modified Kerr–Schild metric and a coordinate transformation

$\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$6

$\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$7

with equatorial refinement and polar “cylindrification” to relax the timestep constraint (Bégué et al., 2022). The later MAD application uses physical horizon-penetrating spherical Kerr-Schild coordinates and numerical modified Kerr-Schild coordinates, explicitly tying the mapping to the cuHARM implementation (Zhang et al., 2023).

Boundary conditions are standard for global accretion-torus calculations: inflow at the inner radial boundary, outflow at the outer radial boundary, periodicity in $\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$8, and reflective polar boundaries in $\nabla_\mu \left ( T^{\mu\nu}\right ) = 0,$9 (Bégué et al., 2022, Zhang et al., 2023). Near the poles, additional damping or zeroing measures are used for stability (Bégué et al., 2022, Wallace et al., 21 Aug 2025).

Like other GRMHD codes evolving magnetized funnels and tenuous coronae, cuHARM relies on floors and ceilings. In the 2022 paper the density and internal-energy floors are

$\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$0

$\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$1

with magnetization limits

$\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$2

and Lorentz-factor cap

$\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$3

(Bégué et al., 2022). The MAD study uses the same flooring model as the code paper and remarks that weak polar channels in some plots are floor artifacts rather than physical structures (Zhang et al., 2023).

The execution model is single-node multi-GPU. Each GPU is assigned to one OpenMP thread, and the domain is decomposed along the $\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$4 direction so that each device advances a predetermined slice of the global domain (Bégué et al., 2022). The code is not described as multi-node or MPI-enabled in the 2022 paper (Bégué et al., 2022). This suggests a design optimized for GPU-rich workstations or servers rather than distributed CPU clusters.

4. Computational performance and hardware profile

The original code paper frames cuHARM primarily as an accessibility-enabling platform for 3D GRMHD. A headline result is that a $\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$5 simulation is “well within the reach of an Nvidia DGX-V100 server,” with the computation being “a factor about 10 times faster if only the CPU was used,” i.e. about tenfold speedup relative to CPU-only execution (Bégué et al., 2022). The same paper states that 3D simulations around $\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$6 take about 72 hours to reach

$\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$7

on a DGX system with $\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$8V100 GPUs, and that the runs presented there require from 2 to 10 days depending on setup and resolution (Bégué et al., 2022).

The performance discussion is unusually explicit about optimization limits. A naive implementation achieved roughly $\nabla_\mu \left ( ^* F^{\mu \nu}\right ) = 0$9 compute throughput and $u_\mu F^{\mu \nu} = 0,$0 memory throughput on an RTX-class gaming GPU, but performance on V100 hardware was much worse because of global-memory bottlenecks (Bégué et al., 2022). The implemented CUDA optimizations included heavy use of shared memory, warp-level primitives, and organization of calculations along the azimuthal direction so that metric data—independent of $u_\mu F^{\mu \nu} = 0,$1 in Kerr—can be reused more efficiently (Bégué et al., 2022). These optimizations reduced execution time by about a factor of $u_\mu F^{\mu \nu} = 0,$2 relative to the naive version, although compute and memory throughput on V100 remained around $u_\mu F^{\mu \nu} = 0,$3 (Bégué et al., 2022).

The 2023 MAD study reports that cuHARM had by then reached “more than 50\% computation efficiency on NVidia A100 cards,” although the production runs in that paper were performed on an Nvidia DGX-V100 server with 8 V100 GPUs (Zhang et al., 2023). The paper also describes the simulations as feasible on a “single multi-GPU workstation” (Zhang et al., 2023). This suggests that the code’s optimization trajectory continued after the original 2022 performance characterization, though no full scaling study is provided there.

The 2025 radiation paper indicates that the radiation module is fully GPU-accelerated but not yet as optimized as the mature MHD kernels (Wallace et al., 21 Aug 2025). It describes two current optimizations—exploiting azimuthal independence of metric and tetrad data in Kerr-Schild setups, and reduction kernels for radiation stress-energy computation—while stating that much larger speedups should still be possible through shared memory, more arithmetic per memory access, and overlap of communication with interior computation (Wallace et al., 21 Aug 2025). In memory terms, the gray radiative accretion application already requires about 40 GB RAM, and the paper estimates that full frequency dependence over $u_\mu F^{\mu \nu} = 0,$4 decades with $u_\mu F^{\mu \nu} = 0,$5 bins per decade would expand memory footprint by a factor of $u_\mu F^{\mu \nu} = 0,$6, to order 2 TB (Wallace et al., 21 Aug 2025).

5. Scientific applications: SANE disks and MAD spin surveys

The first cuHARM paper applies the code to 3D ADAF disks in the SANE state around a Kerr black hole with

$u_\mu F^{\mu \nu} = 0,$7

(Bégué et al., 2022). Nine runs are presented with Fishbone–Moncrief tori, weak poloidal seed fields defined through

$u_\mu F^{\mu \nu} = 0,$8

random 4% pressure perturbations to seed MRI, and resolutions up to $u_\mu F^{\mu \nu} = 0,$9 (Bégué et al., 2022). The paper’s three headline findings are that increasing the magnetic field while remaining in the SANE state does not significantly affect the mass accretion rate, that simultaneous increase of disk size and magnetic field can destroy the jet once the magnetic flux through the horizon decreases below a certain limit, and that jet structure is only a weak function of adiabatic index, with relativistic gas tending to produce a wider jet (Bégué et al., 2022).

A key diagnostic is the normalized horizon magnetic flux,

$b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$0

which remains of order $b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$1 in the SANE runs and therefore well below classical MAD levels (Bégué et al., 2022). The paper reports that sustained jets disappear when $b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$2 becomes too small, with the transition occurring around

$b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$3

in the particular set of runs shown (Bégué et al., 2022). It also finds an empirical relation

$b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$4

indicating that the quasi-steady accretion rate correlates mainly with initial disk mass rather than initial magnetization or adiabatic index in these SANE models (Bégué et al., 2022).

The 2023 paper shifts to the MAD regime and uses 14 global 3D simulations of ADAF accretion “carried with cuHARM” across a wide spin sequence,

$b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$5

(Zhang et al., 2023). The simulations begin from large Fishbone–Moncrief tori with

$b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$6

adiabatic index

$b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$7

and single-loop poloidal magnetic fields normalized by $b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$8, with additional retrograde tests at weaker field strengths (Zhang et al., 2023). Most runs use resolution

$b^\mu \equiv ~^*F^{\mu\nu} u_\nu,$9

with higher-resolution checks at

$$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$0

(Zhang et al., 2023).

That study identifies several spin-dependent properties of the MAD state. The abstract states that jets transport a significant amount of angular momentum to infinity in the form of Maxwell stresses; for high positive spin, the rate of angular-momentum transport is about five times larger than for negative spin, and it is nearly absent for a non-rotating black hole (Zhang et al., 2023). The paper also finds that the mass accretion rate and the MAD parameter, both calculated at the horizon, are not correlated, but their time derivatives are anti-correlated for every spin (Zhang et al., 2023). For zero spin, the toroidal magnetic-pressure contribution is negligible near the horizon compared with the radial component, whereas for fast spinning black holes it becomes comparable and even dominant for high positive spin (Zhang et al., 2023). Negative-spin jets are narrower and more fluctuating than their positive-spin counterparts, while the non-rotating case yields the widest and least variable weak jet (Zhang et al., 2023).

These studies establish cuHARM not merely as a code-description artifact but as a production tool for long-duration 3D GRMHD parameter surveys. This suggests that cuHARM’s practical value lies as much in enabling systematic comparative studies as in individual showcase simulations.

6. Radiation extension and direct covariant transfer

The 2025 extension adds “fully covariant treatment of radiation transport and the subsequent radiation backreaction on the dynamics of the fluid,” transforming cuHARM into a general relativistic radiation magnetohydrodynamics code (Wallace et al., 21 Aug 2025). The conceptual novelty is that the radiation field is not represented by low-order moments plus a closure; instead, cuHARM directly evolves the specific intensity on a discrete angular grid at every spatial cell (Wallace et al., 21 Aug 2025).

The radiation-augmented fluid equations are

$$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$1

$$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$2

$$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$3

where $$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$4 is the radiation four-force density (Wallace et al., 21 Aug 2025). The transfer scheme uses three frames: the coordinate frame, a locally Minkowskian tetrad frame, and the comoving fluid frame (Wallace et al., 21 Aug 2025). Photon direction in the tetrad frame is represented by

$$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$5

and the covariant time-dependent radiative transfer equation is written in flux-conservative form as

$$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$6

(Wallace et al., 21 Aug 2025).

The present implementation is gray. The angular grid is geodesic, constructed from refinements of an icosahedron projected onto the sphere and converted to its Voronoi dual. Generation $$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$7 has

$$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$8

cells, giving 12, 42, 162, 642, 2562, and 10242 directions for G0 through G5 (Wallace et al., 21 Aug 2025). Spatial reconstruction uses the same piecewise linear interpolation as the MHD solver, while angular reconstruction uses minimum-angle plane reconstruction (Wallace et al., 21 Aug 2025).

A central design decision is the use of a deterministic grid-based transport solver rather than M1 or Monte Carlo. The paper states that M1 is correct in optically thin and optically thick asymptotic limits but fails for multi-beam or otherwise non-isotropic fields in the intermediate regime, whereas Monte Carlo handles angular structure naturally in thin media but becomes expensive at high optical depth because of repeated scattering (Wallace et al., 21 Aug 2025). cuHARM’s angularly resolved transport is intended precisely for anisotropic radiation in regions of intermediate optical depth, such as disk–funnel transition regions.

The radiation–matter coupling is operator split into an explicit transport step and an implicit local interaction step in the fluid frame (Wallace et al., 21 Aug 2025). Thermal emissivity is

$$T^{\mu\nu} = (h + b^2) u^\mu u^\nu + \left (p_g + \frac{b^2}{2} \right )g^{\mu \nu} - b^\mu b^\nu$$9

scattering emissivity is

$$U = \sqrt{-g}(\rho u^t, T^t_{~t}, T^{t}_{~i}, B^i),$$0

and the interaction equation in the comoving frame becomes

$$U = \sqrt{-g}(\rho u^t, T^t_{~t}, T^{t}_{~i}, B^i),$$1

(Wallace et al., 21 Aug 2025). The paper develops three temperature-coupling methods because explicit-temperature interaction proved unstable for astrophysical opacities with strong temperature dependence (Wallace et al., 21 Aug 2025). Radiation backreaction is inserted into the GRMHD update via a conservative estimate of the radiation four-force derived from the time change of the radiation stress-energy tensor rather than direct angle-integration of emissivity and absorptivity (Wallace et al., 21 Aug 2025).

The verification suite is extensive and includes spherical light expansion, inverse-square decay, colliding beams, a 2D hohlraum, tetrad-choice comparisons, snake coordinates, photon circular orbits around Kerr black holes, static radiative diffusion, diffusion plus advection, one-cell radiation–matter equilibration, and 1D radiative shocks (Wallace et al., 21 Aug 2025). The black-hole demonstration problem evolves a radiative Fishbone–Moncrief torus with

$$U = \sqrt{-g}(\rho u^t, T^t_{~t}, T^{t}_{~i}, B^i),$$2

resolution $$U = \sqrt{-g}(\rho u^t, T^t_{~t}, T^{t}_{~i}, B^i),$$3, G2 angular grid, and runtime to $$U = \sqrt{-g}(\rho u^t, T^t_{~t}, T^{t}_{~i}, B^i),$$4 (Wallace et al., 21 Aug 2025). Compared to a non-radiative SANE torus, the radiative disk becomes thinner and denser, with mean accretion rate

$$U = \sqrt{-g}(\rho u^t, T^t_{~t}, T^{t}_{~i}, B^i),$$5

and disk midplane cooling to $$U = \sqrt{-g}(\rho u^t, T^t_{~t}, T^{t}_{~i}, B^i),$$6 K (Wallace et al., 21 Aug 2025).

7. Scope, limitations, and significance

Across its documented development, cuHARM has moved from a GPU-native ideal-GRMHD code for black-hole accretion and jets (Bégué et al., 2022) to a tool used for systematic 3D MAD spin surveys (Zhang et al., 2023), and then to a directly radiative GRMHD framework with angularly resolved transport (Wallace et al., 21 Aug 2025). Several strengths recur across these stages: a conservative finite-volume GRMHD core, horizon-penetrating Kerr coordinates, practical multi-GPU execution on a single server, and explicit attention to diagnostics relevant for accretion and jet physics.

The published limitations are equally clear. In the 2022 formulation, cuHARM is single-node rather than multi-node, uses PLM rather than PPM, retains flux-CT with its known drawbacks, and omits radiation, resistivity, electron thermodynamics, and dynamical spacetime (Bégué et al., 2022). The MAD application paper does not restate algorithmic internals such as reconstruction order, Riemann solver, constrained-transport details, or primitive recovery method, and explicitly points readers to the earlier code paper for those elements (Zhang et al., 2023). The 2025 radiation module remains gray, assumes isotropic and elastic scattering in the fluid frame, omits a full frequency-dependent Compton kernel, and imposes high memory and compute costs because direct angular transport effectively raises the problem dimensionality from 3D space to 5D in the gray case (Wallace et al., 21 Aug 2025).

A common misconception would be to equate cuHARM solely with its original 2022 code-paper implementation. The later literature shows that cuHARM denotes both a production GRMHD code used for extensive MAD studies and, after 2025, a broader radiation-GRMHD platform (Zhang et al., 2023, Wallace et al., 21 Aug 2025). Another plausible misconception is that its GPU acceleration automatically implies fully saturated hardware efficiency; the papers are explicit that performance gains are substantial but that the code remains bandwidth-limited in important kernels and that the radiation side is not yet fully optimized (Bégué et al., 2022, Wallace et al., 21 Aug 2025).

Within computational relativistic astrophysics, cuHARM occupies a specific methodological position. It inherits HARM-style conservative GRMHD numerics and adapts them to GPU-centric hardware (Bégué et al., 2022). It has been used to recover and extend MAD phenomenology across spin, especially in angular-momentum transport and magnetic-pressure decomposition (Zhang et al., 2023). Its radiative extension chooses direct angular transport over closure-based moment schemes, prioritizing anisotropic and intermediate-optical-depth fidelity at the cost of substantial memory and arithmetic expense (Wallace et al., 21 Aug 2025). This suggests a codebase oriented toward technically demanding, physics-rich accretion problems where single-node GPU acceleration can substitute for larger CPU-centric infrastructure.