CRK-HACC: GPU Cosmology & Hydro Framework
- CRK-HACC is a GPU-accelerated cosmological simulation framework that couples HACC's N-body gravity solver with CRKSPH to capture baryonic gas dynamics.
- It integrates separation-of-scale gravity solvers, mesh-free hydrodynamics, in situ analysis, and multi-tiered I/O to efficiently produce survey-scale synthetic skies at exascale.
- The framework incorporates calibrated subgrid physics for radiative cooling, star formation, and AGN feedback, addressing observational benchmarks in large-scale structure studies.
CRK-HACC is a GPU-accelerated, particle-based cosmological simulation framework for hydrodynamics in large-scale structure formation. It extends the Hardware/Hybrid Accelerated Cosmology Code (HACC) by coupling the HACC gravitational N-body solver to Conservative Reproducing-Kernel Smoothed Particle Hydrodynamics, denoted CRKSPH or CRK-SPH, in order to resolve gas hydrodynamics alongside dark matter and thereby model baryonic effects in cosmology simulations (Frontiere et al., 2022). Across its published descriptions, CRK-HACC is presented as a codesigned framework for modern GPU-accelerated and exascale supercomputers, combining separation-of-scale gravity solvers, mesh-free hydrodynamics, in situ analysis, and multi-tiered I/O for survey-scale synthetic-sky production (Rangel et al., 2023, Frontiere et al., 3 Oct 2025, Frontiere et al., 26 Nov 2025).
1. Origin, scientific objective, and relation to HACC
CRK-HACC was introduced as an extension of HACC to resolve gas hydrodynamics in simulations of the universe’s large-scale structure (Frontiere et al., 2022). The 2022 formulation emphasizes that the new framework couples the HACC gravitational N-body solver with a modern SPH approach called CRKSPH, whose defining property is the use of smoothing functions that exactly interpolate linear fields while manifestly preserving conservation laws for momentum, mass, and energy (Frontiere et al., 2022). The stated motivation is accurate modeling of baryonic effects for the generation of precise synthetic sky predictions for upcoming observational surveys (Frontiere et al., 2022).
Later work situates CRK-HACC within a broader production workflow. In the performance-portability study, HACC is described as a particle-based N-body cosmology code that has run on DOE leadership machines for over a decade, simulating the growth of structure under gravity, while CRK-HACC extends it by adding baryonic hydrodynamic physics through CRKSPH (Rangel et al., 2023). That study further states that production runs such as Borg Cube and Farpoint typically evolve particles, with equal numbers of dark matter and gas, from redshift to , coupling GPU-accelerated short-range gravity and hydrodynamics kernels with a long-range particle-mesh Poisson solver (Rangel et al., 2023).
In the exascale and galaxy-formation descriptions, CRK-HACC is further characterized as a cosmological hydrodynamics code built for the extreme scalability requirements set by modern cosmological surveys (Frontiere et al., 3 Oct 2025), and as a framework extended with radiative cooling, star formation, stellar evolution, and AGN feedback to model baryonic effects self-consistently at survey scale (Frontiere et al., 26 Nov 2025). Taken together, these descriptions identify CRK-HACC not merely as a hydrodynamics add-on to HACC, but as a full cosmological hydro-gravity framework whose primary scientific role is end-to-end modeling of structure formation with baryons.
2. Hydrodynamic formulation: CRKSPH and conservative reproducing kernels
The hydrodynamic core of CRK-HACC is CRKSPH, a Conservative Reproducing-Kernel Smoothed Particle Hydrodynamics method (Frontiere et al., 2022, Rangel et al., 2023). In the SYCL implementation study, the reproducing kernel is written as
with the scalar and vector chosen so that zeroth- and first-order consistency conditions hold:
(Rangel et al., 2023). The same study states that these corrected kernels allow conservative discretization of fields, and that the resulting equations conserve mass, momentum, and energy to machine precision aside from time-integration error (Rangel et al., 2023).
The exascale paper expresses the same idea in an alternative notation: where is a polynomial prefactor chosen so that
0
(Frontiere et al., 3 Oct 2025). The associated density estimate is
1
with momentum and energy equations given by
2
3
(Frontiere et al., 3 Oct 2025).
The galaxy-formation paper gives a more implementation-specific form using the Wendland 4 kernel,
5
and writes the momentum equation in compatible-energy form as
6
with
7
and the energy equation
8
(Frontiere et al., 26 Nov 2025). In that account, artificial viscosity uses a limited Monaghan-type form to capture shocks without smearing smooth flows, and mass conservation follows directly from the SPH summation for 9, with individual particle masses modified only by subgrid processes such as winds, enrichment, and accretion (Frontiere et al., 26 Nov 2025).
A recurrent misconception is to assimilate CRK-HACC to conventional SPH without qualification. The published formulations instead stress linear reproduction, corrected kernels, and manifest conservation as the numerical signature of CRKSPH (Frontiere et al., 2022, Rangel et al., 2023, Frontiere et al., 26 Nov 2025). A plausible implication is that the framework is intended to address known consistency limitations of standard SPH while preserving the conservative structure that makes particle hydrodynamics attractive in cosmological settings.
3. Gravity–hydrodynamics coupling, timestepping, and stabilization
CRK-HACC couples hydrodynamics to gravity through a separation-of-scales strategy. The exascale description decomposes the gravitational potential as
0
with the long-range component computed globally on a uniform FFT grid of size 1 via a spectrally filtered particle-mesh solver using SWFFT, and the short-range component evaluated locally via a GPU-resident tree that includes both dark-matter pairwise forces and hydrodynamic accelerations from CRK-SPH (Frontiere et al., 3 Oct 2025). The same account states that a high-order filter 2 suppresses aliasing and ensures smooth hand-off to the short-range solver, and that the two-scale approach permits coarse PM steps in FP64 and fine local interactions in FP32 on the GPU with negligible loss of accuracy (Frontiere et al., 3 Oct 2025).
The time-integration strategy is hierarchical. The exascale paper states that the CRK-SPH equations are integrated with a hierarchical split-timestep scheme following Saitoh and Makino, grouping particles into bins by local CFL and feedback constraints so that only active particles advance at each subcycle (Frontiere et al., 3 Oct 2025). The galaxy-formation paper describes a reversible Kick-Drift-Kick symplectic integrator with hierarchical, power-of-two subcycling for hydrodynamic timesteps, in which long-range PM forces are updated at coarse PM steps in scale factor while short-range gravity and hydrodynamics subcycle according to local CFL-limited timesteps 3 (Frontiere et al., 26 Nov 2025).
That paper also details operator coupling for subgrid physics by first-order Strang splitting at each level, schematically
4
(Frontiere et al., 26 Nov 2025). To avoid unphysical smoothing-length oscillations under deep subcycling, it applies a combined smoothing-length derivative,
5
where 6 is the neighbor-count-corrected length and both active and passive particles update 7 smoothly (Frontiere et al., 26 Nov 2025). It further introduces a divergence-based regularization that disables the reproducing-kernel correction when
8
reverting to standard SPH to prevent pathological neighbor configurations (Frontiere et al., 26 Nov 2025).
These details clarify that CRK-HACC is neither a pure gravity code with post-processed baryons nor a hydro code with an external gravity module. Its published design is an explicitly coupled gravity–hydro system in which long-range PM evolution, local tree interactions, SPH hydrodynamics, and subgrid operators are orchestrated within a multirate timestep hierarchy.
4. GPU-resident implementation and performance portability
CRK-HACC inherits the codesign strategies of the HACC solver and is built to run on modern GPU-accelerated supercomputers (Frontiere et al., 2022). The galaxy-formation description states that the code structure uses a hybrid Tree-PM gravity solver with long-range forces via a 3D mesh and FFT through SWFFT and short-range forces via a GPU-resident Barnes–Hut tree with highly optimized neighbor search (Frontiere et al., 26 Nov 2025). Hydrodynamic and subgrid kernels are implemented fully on GPU, particle data is stored in structure-of-arrays format for coalesced memory access, and particle overloading duplicates boundary regions to neighboring MPI ranks so that hydrodynamic interactions proceed without synchronous communication until the next PM step (Frontiere et al., 26 Nov 2025).
The exascale paper gives a more detailed algorithmic account of the local solver. Each MPI rank overloads its domain by a fixed chaining-mesh buffer four PM-cells wide, and within each bin a shallow 9-d tree subdivides particles into leaves of 0 particles (Frontiere et al., 3 Oct 2025). At each global PM step, the chaining-mesh and tree build is 1 on the host and approximately 2 of total time-to-solution, while GPU kernels compute leaf-to-leaf interactions for pairwise short-range gravity and hydrodynamic sums (Frontiere et al., 3 Oct 2025). The same paper identifies a warp-splitting optimization in which half-warps load “i” and “j” states once, use register-shuffle instructions to exchange partials, and accumulate 3 without redundant memory loads; this reduces register pressure, minimizes global loads, and localizes atomics to per-leaf reductions (Frontiere et al., 3 Oct 2025).
The 2023 performance-portability study concentrates on the migration from CUDA to SYCL for GPUs from AMD, Intel, and NVIDIA (Rangel et al., 2023). It reports that CRK-HACC’s short-range gravity and hydrodynamics kernels were originally hand-tuned CUDA with approximately 30 kSLOC of device code, and that each timestep invokes roughly five hot kernels—Geometry, Corrections, Extras, Acceleration, and Energy—that account for at least 4 of GPU time (Rangel et al., 2023). The migration employed SYCLomatic for CUDA-to-SYCL translation, supplemented by a small Clang LibTooling pass that converted migrated kernels into function objects so they could be passed directly to parallel_for without lambdas (Rangel et al., 2023).
Subgroup management is a central issue in that study. Because CUDA uses warps of size 32, whereas SYCL requires explicit subgroup specification for deterministic behavior, the implementation chose subgroup size 5 on NVIDIA A100, 6 on AMD MI250X, and ultimately 7 on Intel GPUs to balance register pressure against occupancy (Rangel et al., 2023). The study further developed four variants of a critical half-warp leaf-interaction pattern: Select, Memory, Broadcast, and vISA (Rangel et al., 2023). It then evaluates portability using the harmonic-mean efficiency metric
8
with
9
and code divergence using the Jaccard distance
0
Measured on the GPU-only time for five timesteps of a 1-particle adiabatic test problem on one node with eight MPI ranks, the SYCL version of CRK-HACC achieves a performance portability of 2 with a code divergence of almost 3 (Rangel et al., 2023). More specifically, the mixed SYCL strategy Select+vISA yields 4, while Memory-only yields 5, Select-only yields 6, and the cross-language CUDA+HIP+SYCL combination yields 7 (Rangel et al., 2023). The same work reports code convergence of approximately 8 for the mixed-SYCL Select+vISA approach, compared with approximately 9 for a full CUDA+HIP+SYCL approach (Rangel et al., 2023). This directly addresses another common misconception: single-source portability in CRK-HACC does not mean zero specialization. The published result is instead that small, targeted specializations can greatly improve performance portability without significantly impacting programmer productivity (Rangel et al., 2023).
5. Exascale realization: Frontier-E, in situ analysis, and multi-tiered I/O
The largest published CRK-HACC deployment is the Frontier-E full-sky simulation (Frontiere et al., 3 Oct 2025). That run used a 0 Gpc comoving box, 1 total particles, and equal baryon and dark-matter tracers, with a PM grid of 2 and Planck-like 3CDM cosmology 4 (Frontiere et al., 3 Oct 2025). The published performance figures are 5 PFLOPs peak, 6 PFLOPs sustained, and throughput of 7 particles s8 on 9 Frontier nodes, corresponding to 0 GPU dies (Frontiere et al., 3 Oct 2025).
The same paper reports weak scaling of 1 efficiency from 2 to 3 nodes and strong scaling of 4 efficiency on a fixed 5 grid, with scaling laws described as essentially ideal, time proportional to 6 in weak scaling (Frontiere et al., 3 Oct 2025). It states that the overall local complexity per PM step is 7, and that GPU utilization was sustained at 8–9 across 0 Frontier nodes for the local tree solver (Frontiere et al., 3 Oct 2025).
CRK-HACC’s exascale design includes a fully GPU-resident in situ analysis pipeline because performing halo finding and clustering after the fact at more than 1 PB of raw data is deemed infeasible (Frontiere et al., 3 Oct 2025). Embedded analyses include Friends-of-Friends and DBSCAN halo finders via the ArborX library, mock-survey calculators such as light-cone assembly and SZ and X-ray synthetic maps, and baryon/dark-matter field statistics including power spectra and PDFs (Frontiere et al., 3 Oct 2025). Executed immediately after each PM step on the device, this analysis consumes 2 of time-to-solution, whereas short-range forces consume 3 (Frontiere et al., 3 Oct 2025). The galaxy-formation account adds that in situ analysis also includes a DBSCAN galaxy finder and FOF/SO halos via ArborX on GPU, enabling on-the-fly AGN seeding and galaxy property catalogs (Frontiere et al., 26 Nov 2025).
For data management, CRK-HACC uses a decentralized, node-local staging strategy with synchronous writes to local NVMe SSD, asynchronous background bleed of complete files to the Orion Lustre PFS via OS move, and rolling retention and purge of old checkpoints (Frontiere et al., 3 Oct 2025). The exascale paper reports checkpoint volumes of approximately 4–5 TB per full PM step, more than 6 PB total checkpoint data, and 7 PB of scientific outputs (Frontiere et al., 3 Oct 2025). It further reports an effective aggregated write bandwidth of 8 TB/s to Orion, exceeding its 9 TB/s peak, while never stalling the solver, with only 0 of total time-to-solution spent in I/O (Frontiere et al., 3 Oct 2025).
The exascale literature therefore presents CRK-HACC as a system-level framework in which solver design, analysis, and I/O are co-optimized. This suggests that the code’s significance for cosmological surveys lies not only in numerical hydro accuracy but also in the practical ability to produce survey-scale outputs and diagnostics within operational runtime and storage constraints.
6. Astrophysical subgrid modeling, calibration, and scientific outputs
The 2025 galaxy-formation extension augments CRK-HACC with a suite of subgrid models for radiative cooling, star formation, stellar evolution, and AGN feedback (Frontiere et al., 26 Nov 2025). For radiative cooling and heating, gas is treated as optically thin and in ionization equilibrium with a uniform UV/X-ray background, attenuated in self-shielded regions according to
1
while metal-line and primordial cooling rates are tabulated with CLOUDY on a 5-D grid and applied using Townsend’s exact integration method,
2
(Frontiere et al., 26 Nov 2025).
The star-formation model uses a two-phase ISM following Springel and Hernquist, with cold clouds at 3 K embedded in a hot phase 4 K and effective pressure
5
with 6 and 7 K (Frontiere et al., 26 Nov 2025). Gas forms stars stochastically when 8 and 9, with
0
integrated probabilistically through 1 (Frontiere et al., 26 Nov 2025). The model also includes metallic floors 2 and 3 to mimic early enrichment at coarse resolution (Frontiere et al., 26 Nov 2025).
Galactic winds are kinetic and decoupled, launched at rate
4
with 5 and calibrated 6 (Frontiere et al., 26 Nov 2025). Wind velocity scales with local dark-matter velocity dispersion,
7
(Frontiere et al., 26 Nov 2025).
Chemical enrichment is modeled through star particles representing single stellar populations, with mass loss from Type Ia and core-collapse supernovae and OB/AGB winds obtained by analytically integrating fits to FIRE-3 rates (Frontiere et al., 26 Nov 2025). For black-hole growth and AGN feedback, black holes are seeded in galaxies identified in situ once they exceed the seed-mass criterion, with particle mass 8 and internal mass 9 (Frontiere et al., 26 Nov 2025). Gas accretion follows Bondi–Hoyle,
00
with radiative efficiency 01 and fixed 02, and the feedback prescription uses high-accretion thermal and low-accretion kinetic modes (Frontiere et al., 26 Nov 2025).
Calibration was carried out at baryon mass 03 in an 04 Mpc box via a 05-point Latin hypercube over five parameters: 06, 07, 08, 09, and 10 (Frontiere et al., 26 Nov 2025). The fiducial values are 11, 12, 13, 14, and 15 (Frontiere et al., 26 Nov 2025). Calibration targets were galaxy stellar mass functions at 16 and low-redshift cluster gas-density profiles (Frontiere et al., 26 Nov 2025).
The exascale paper states that Frontier-E includes radiative and metal-line cooling through the GRACKLE library, star formation and supernova feedback, stellar chemical enrichment, and AGN feedback, calibrated on Perlmutter mid-scale runs to reproduce galactic stellar mass functions, cluster gas fractions, and IGM thermodynamics (Frontiere et al., 3 Oct 2025). Its reported scientific outputs include approximately 17 resolved galaxy clusters, continuous light-cone outputs for DESI, Euclid, LSST, Roman, and SPHEREx, synthetic SZ, X-ray, optical, IR, and radio maps computed in situ, and full redshift coverage 18 for weak lensing, BAO, RSD, tSZ, and kSZ forecasts (Frontiere et al., 3 Oct 2025).
The comparison study reports that CRK-HACC matches observational data and UniverseMachine in the GSMF at 19, closely tracks other simulations in cosmic star-formation-rate density and cosmic stellar-mass density at 20, reproduces local sSFR and quenched-fraction trends with systematic offsets at high 21 comparable to other codes, follows abundance-matching determinations in the stellar-mass–halo-mass relation, reproduces the SDSS stellar mass–metallicity relation, matches the Illustris-TNG black-hole-mass–stellar-mass slope, and yields halo gas fractions near cluster data when no hydrostatic mass correction is applied (Frontiere et al., 26 Nov 2025). These are validation statements rather than a claim of exact observational agreement in every diagnostic; the published framing is explicitly comparative and benchmark-driven.
A final misconception is that CRK-HACC’s recent development concerns only hardware portability. The combined literature shows a broader trajectory: the 2022 introduction establishes the hydro extension to HACC (Frontiere et al., 2022); the 2023 work addresses portability across AMD, Intel, and NVIDIA GPUs (Rangel et al., 2023); the 2025 exascale report demonstrates a four-trillion-particle full-sky run with integrated analysis and I/O (Frontiere et al., 3 Oct 2025); and the 2025 galaxy-formation paper adds calibrated subgrid physics for survey-scale baryonic modeling (Frontiere et al., 26 Nov 2025). The framework’s defining characteristic is therefore the joint treatment of numerical hydrodynamics, gravity coupling, exascale implementation, and astrophysical modeling within a single production cosmology code.