--- title: 'CRK-HACC: GPU Cosmology & Hydro Framework' url: https://www.emergentmind.com/topics/crk-hacc type: topic --- # CRK-HACC: GPU Cosmology & Hydro Framework 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 [2202.02840]. 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 [2310.16122], [2510.03557], [2511.21921]. ## 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 [2202.02840]. 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 [2202.02840]. The stated motivation is accurate modeling of baryonic effects for the generation of precise synthetic sky predictions for upcoming observational surveys [2202.02840]. 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 [2310.16122]. That study further states that production runs such as Borg Cube and Farpoint typically evolve \(2\times N^3\) particles, with equal numbers of dark matter and gas, from redshift \(z\approx 200\) to \(z\approx 0\), coupling GPU-accelerated short-range gravity and hydrodynamics kernels with a long-range particle-mesh Poisson solver [2310.16122]. 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 [2510.03557], and as a framework extended with radiative cooling, star formation, stellar evolution, and AGN feedback to model baryonic effects self-consistently at survey scale [2511.21921]. 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 [2202.02840], [2310.16122]. In the SYCL implementation study, the reproducing kernel is written as \[ W^R_{ij} = A_i\,W_{ij} + \mathbf{B}_i\cdot\mathbf{r}_{ij}\,W_{ij}, \] with the scalar \(A_i\) and vector \(\mathbf{B}_i\) chosen so that zeroth- and first-order consistency conditions hold: \[ \sum_j m_j\bigl(A_i\,W_{ij} + \mathbf{B}_i\!\cdot\!\mathbf{r}_{ij}\,W_{ij}\bigr)=1, \] \[ \sum_j m_j\,\mathbf{r}_{ij}\,\bigl(A_i\,W_{ij} + \mathbf{B}_i\!\cdot\!\mathbf{r}_{ij}\,W_{ij}\bigr)=\mathbf{0} \] [2310.16122]. 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 [2310.16122]. The exascale paper expresses the same idea in an alternative notation: \[ W^{\mathrm{CRK}}(r_{ij},h_i) = \Psi_i(\mathbf{r}_{ij})\,W(r_{ij},h_i), \] where \(\Psi_i(\mathbf{r})\) is a polynomial prefactor chosen so that \[ \sum_j m_j\,\mathbf{r}_{ij}\,W^{\mathrm{CRK}}(r_{ij},h_i)=\mathbf{0}, \qquad \sum_j m_j\,W^{\mathrm{CRK}}(r_{ij},h_i)=\rho_i \] [2510.03557]. The associated density estimate is \[ \rho_i=\sum_j m_j\,W^{\mathrm{CRK}}(r_{ij},h_i), \] with momentum and energy equations given by \[ \frac{d\mathbf{v}_i}{dt} = -\sum_j m_j \Bigl( \frac{P_i}{\rho_i^2} + \frac{P_j}{\rho_j^2} + \Pi_{ij} \Bigr)\nabla_i W^{\mathrm{CRK}}(r_{ij},h_{ij}), \] \[ \frac{du_i}{dt} = \frac{1}{2}\sum_j m_j \Bigl( \frac{P_i}{\rho_i^2} + \frac{P_j}{\rho_j^2} + \Pi_{ij} \Bigr) (\mathbf{v}_i-\mathbf{v}_j)\cdot\nabla_i W^{\mathrm{CRK}}(r_{ij},h_{ij}) \] [2510.03557]. The galaxy-formation paper gives a more implementation-specific form using the Wendland \(C^4\) kernel, \[ K_{ij}=A_i\bigl[1+B_i\cdot x_{ij}\bigr]\,W(|x_{ij}|,h_i), \] and writes the momentum equation in compatible-energy form as \[ \frac{d\mathbf{v}_i}{dt} = -\sum_j m_j \Bigl[ \frac{P_i}{\rho_i^2}F_{ij}(h_i)+\frac{P_j}{\rho_j^2}F_{ji}(h_j) \Bigr] +\mathbf{a}_i^{\rm visc}, \] with \[ P_i=(\gamma-1)\rho_i u_i,\qquad \rho_i=\sum_j m_j K_{ij}, \] and the energy equation \[ \frac{du_i}{dt} = \frac{1}{2}\sum_j m_j \Bigl[ \frac{P_i}{\rho_i^2}F_{ij}(h_i)+\frac{P_j}{\rho_j^2}F_{ji}(h_j) \Bigr] (\mathbf{v}_i-\mathbf{v}_j)\cdot\hat x_{ij} +\Lambda_i(\rho_i,u_i,Z_i) \] [2511.21921]. 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 \(\rho_i\), with individual particle masses modified only by subgrid processes such as winds, enrichment, and accretion [2511.21921]. 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 [2202.02840], [2310.16122], [2511.21921]. 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 \[ \Phi(\mathbf{x})=\Phi^{\rm LR}(\mathbf{x})+\Phi^{\rm SR}(\mathbf{x}), \] with the long-range component computed globally on a uniform FFT grid of size \(N_{\rm grid}=12\,600^3\) 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 [2510.03557]. The same account states that a high-order filter \(G(k)\) 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 [2510.03557]. 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 [2510.03557]. 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 \(\Delta t_i\) [2511.21921]. That paper also details operator coupling for subgrid physics by first-order Strang splitting at each level, schematically \[ U_H(\Delta t_H)\to U_H(\Delta t_H)\circ U_H^{\rm cool}(\Delta t_H)\circ U_H^{\rm SF+wind}(\Delta t_H)\circ\dots \] [2511.21921]. To avoid unphysical smoothing-length oscillations under deep subcycling, it applies a combined smoothing-length derivative, \[ \dot h_i = \frac{h_i^{\rm T}-h_i} {\Delta t_{\rm a}+\tfrac{1}{3}\,h_i\,(\nabla\!\cdot\!v)_i}, \] where \(h_i^T\) is the neighbor-count-corrected length and both active and passive particles update \(h_i\) smoothly [2511.21921]. It further introduces a divergence-based regularization that disables the reproducing-kernel correction when \[ \exp[|\nabla\!\cdot\!v_i|\Delta t] > 8, \] reverting to standard SPH to prevent pathological neighbor configurations [2511.21921]. 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 [2202.02840]. 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 [2511.21921]. 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 [2511.21921]. 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 \(k\)-d tree subdivides particles into leaves of \(\mathcal{O}(10^2)\) particles [2510.03557]. At each global PM step, the chaining-mesh and tree build is \(\mathcal{O}(N)\) on the host and approximately \(1.7\%\) of total time-to-solution, while GPU kernels compute leaf-to-leaf interactions for pairwise short-range gravity and hydrodynamic sums [2510.03557]. 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 \(\phi_{ij}\) without redundant memory loads; this reduces register pressure, minimizes global loads, and localizes atomics to per-leaf reductions [2510.03557]. The 2023 performance-portability study concentrates on the migration from CUDA to SYCL for GPUs from AMD, Intel, and NVIDIA [2310.16122]. 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 \(85\%\) of GPU time [2310.16122]. 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 [2310.16122]. 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 \(S=32\) on NVIDIA A100, \(S=64\) on AMD MI250X, and ultimately \(S=16\) on Intel GPUs to balance register pressure against occupancy [2310.16122]. The study further developed four variants of a critical half-warp leaf-interaction pattern: Select, Memory, Broadcast, and vISA [2310.16122]. It then evaluates portability using the harmonic-mean efficiency metric \[ \Phi(a,p,H)= \begin{cases} \dfrac{|H|}{\sum_{i\in H}1/e_i(a,p)} & \text{if }e_i(a,p)\neq 0\ \forall i\in H,\\[6pt] 0 & \text{otherwise}, \end{cases} \] with \[ e_i(a,p)=\frac{T_{\mathrm{best},i}}{T_i(a,p)}, \] and code divergence using the Jaccard distance \[ d_{i,j}(a,p)=1-\frac{|c_i\cap c_j|}{|c_i\cup c_j|}, \qquad \mathrm{CD}(a,p,H)=\binom{|H|}{2}^{-1}\sum_{in_{\rm sfr}=0.13\ {\rm cm}^{-3}\) and \(u\leq u_{\rm ISM}\), with \[ \dot M_* = x_c\,\frac{m}{t_{\rm sfr}}, \qquad t_{\rm sfr}(\rho)\propto n_{\rm H}^{-1/2}, \] integrated probabilistically through \(p_{\rm star}=1-\exp[-(\dot m/m)\Delta t]\) [2511.21921]. The model also includes metallic floors \(Z_{\rm ISM,min}=0.25 Z_\odot\) and \(Y_{\rm ISM,min}\) to mimic early enrichment at coarse resolution [2511.21921]. Galactic winds are kinetic and decoupled, launched at rate \[ \dot M_{\rm w}=\eta_{\rm w}\,\dot M_*, \qquad \eta_{\rm w}=\frac{2E_{\rm w}}{v_{\rm w}^2}, \qquad E_{\rm w}=e_{\rm w}E_0, \] with \(E_0=1.73\times 10^{49}\ {\rm erg}\ M_\odot^{-1}\) and calibrated \(e_{\rm w}=0.5\) [2511.21921]. Wind velocity scales with local dark-matter velocity dispersion, \[ v_{\rm w}=\kappa_{\rm w}\,\sigma_{\rm DM}\,\Bigl[\frac{H_0}{H(z)}\Bigr]^{1/3}, \qquad \kappa_{\rm w}=3.0 \] [2511.21921]. 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 [2511.21921]. 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 \(m_{\rm BH}\sim 2\times 10^8 M_\odot\) and internal mass \(M_{{\rm BH},0}=8\times 10^5 h^{-1}M_\odot\) [2511.21921]. Gas accretion follows Bondi–Hoyle, \[ \dot M_{\rm Bondi} = 4\pi G^2 M_{\rm BH}^2 \frac{\rho}{(c_s^2+v_{\rm rel}^2)^{3/2}}, \qquad \dot M_{\rm BH}=\min(\dot M_{\rm Bondi},\dot M_{\rm Edd}), \] with radiative efficiency \(\epsilon_r=0.2\) and fixed \(v_{\rm rel}=8\ {\rm km\ s}^{-1}\), and the feedback prescription uses high-accretion thermal and low-accretion kinetic modes [2511.21921]. Calibration was carried out at baryon mass \(m_b\sim 2\times 10^8 h^{-1}M_\odot\) in an \(L=128 h^{-1}\) Mpc box via a \(64\)-point Latin hypercube over five parameters: \(\kappa_w\), \(e_w\), \(M_{\rm seed}\), \(\epsilon_{\rm kin}\), and \(v_{\rm jet}\) [2511.21921]. The fiducial values are \(\kappa_w=3.0\), \(e_w=0.5\), \(M_{\rm seed}=8\times 10^5 h^{-1}M_\odot\), \(\epsilon_{\rm kin}=1.3\), and \(v_{\rm jet}=5.1\ {\rm km\ s}^{-1}\) [2511.21921]. Calibration targets were galaxy stellar mass functions at \(z=0,1,2\) and low-redshift cluster gas-density profiles [2511.21921]. 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 [2510.03557]. Its reported scientific outputs include approximately \(570{,}000\) 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 \(z\approx 9\to 0\) for weak lensing, BAO, RSD, tSZ, and kSZ forecasts [2510.03557]. The comparison study reports that CRK-HACC matches observational data and UniverseMachine in the GSMF at \(z=0,1,2\), closely tracks other simulations in cosmic star-formation-rate density and cosmic stellar-mass density at \(z\lesssim 3\), reproduces local sSFR and quenched-fraction trends with systematic offsets at high \(z\) 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 [2511.21921]. 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 [2202.02840]; the 2023 work addresses portability across AMD, Intel, and NVIDIA GPUs [2310.16122]; the 2025 exascale report demonstrates a four-trillion-particle full-sky run with integrated analysis and I/O [2510.03557]; and the 2025 galaxy-formation paper adds calibrated subgrid physics for survey-scale baryonic modeling [2511.21921]. 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. Source: https://www.emergentmind.com/topics/crk-hacc