Source-linked AI summary
HACC: Simulating Sky Surveys on State-of-the-Art Supercomputing Architectures
Salman Habib, Adrian Pope, Hal Finkel, Nicholas Frontiere, Katrin Heitmann, David Daniel, Patricia Fasel, Vitali Morozov, George Zagaris, Tom Peterka, Venkatram Vishwanath, Zarija Lukic, Saba Sehrish, Wei-keng Liao
TL;DR
Large-scale sky surveys generate data for studying dark energy, dark matter, primordial fluctuations, and neutrino properties, creating a need for large cosmological simulations. HACC provides a portable simulation framework with in situ analysis and architecture-specific optimizations, supporting applications from synthetic catalogs to simulation databases while achieving strong scaling results on leadership systems.
Problem
Sky surveys produce unprecedented data for investigating dark energy, dark matter, primordial fluctuations, neutrino properties, and astrophysical structure formation.
Method
HACC combines architecture-portable cosmological N-body simulation with optimized force calculations and in situ analysis tools for reducing and interpreting very large outputs.
Results
HACC supports large simulation campaigns for cosmological emulators, covariance matrices, survey optimization, synthetic catalogs, and MCMC-based parameter estimation, including a ∼3.6 trillion-particle simulation.
Takeaways & Limitations
HACC is designed to support precision cosmological predictions, survey strategy optimization, pipeline validation, and broad exploitation of large simulation databases.
Takeaways & Limitations
HACC’s ongoing development still includes future-architecture algorithms and additional physics and analysis capabilities such as gas physics, feedback, and more sophisticated galaxy modeling.
Abstract
from arXiv · showhide
Current and future surveys of large-scale cosmic structure are associated with a massive and complex datastream to study, characterize, and ultimately understand the physics behind the two major components of the 'Dark Universe', dark energy and dark matter. In addition, the surveys also probe primordial perturbations and carry out fundamental measurements, such as determining the sum of neutrino masses. Large-scale simulations of structure formation in the Universe play a critical role in the interpretation of the data and extraction of the physics of interest. Just as survey instruments continue to grow in size and complexity, so do the supercomputers that enable these simulations. Here we report on HACC (Hardware/Hybrid Accelerated Cosmology Code), a recently developed and evolving cosmology N-body code framework, designed to run efficiently on diverse computing architectures and to scale to millions of cores and beyond. HACC can run on all current supercomputer architectures and supports a variety of programming models and algorithms. It has been demonstrated at scale on Cell- and GPU-accelerated systems, standard multi-core node clusters, and Blue Gene systems. HACC's design allows for ease of portability, and at the same time, high levels of sustained performance on the fastest supercomputers available. We present a description of the design philosophy of HACC, the underlying algorithms and code structure, and outline implementation details for several specific architectures. We show selected accuracy and performance results from some of the largest high resolution cosmological simulations so far performed, including benchmarks evolving more than 3.6 trillion particles.
1. Introduction
Large sky surveys create demanding cosmological and astrophysical data-analysis needs, requiring large, accurate simulations and codes that exploit evolving supercomputers. HACC addresses these needs through portable, architecture-flexible algorithms, in situ analysis, and demonstrated trillion-particle-scale performance.
- Sky surveys probe dark energy, dark matter, primordial fluctuations, neutrino properties, and astrophysical structure formation.
- Large-scale simulations provide predictions, systematic studies, synthetic catalogs, and other modeling needed to extract information from surveys.
- HACC was designed for portability across branching supercomputer architectures while supporting multiple programming paradigms and sustained performance.
- HACC combines architecture-independent long-range computation with tunable implementations and supports in situ tools for reducing petabyte-scale outputs.
2. General Features of the HACC Code Framework
Cosmological N-body simulations replace a difficult high-dimensional Vlasov-Poisson calculation with particle dynamics and hybrid force solvers. HACC combines an architecture-independent spectral particle-mesh layer with architecture-tunable short-range methods.
- The Vlasov-Poisson equation is computationally difficult because of high dimensionality and nonlinear multistreaming, motivating N-body particle methods.
- Trillions of particles and large spatial dynamic ranges rule out brute-force particle-particle force calculations.
- HACC splits gravity into an architecture-independent spectral particle-mesh long/medium-range solver and an architecture-tunable particle-based short-range solver.
- Short-range implementations use direct particle-particle, tree, or combined PP-TreePM methods according to the target architecture.
- The long-range solver remains unchanged across architectures, while short-range solvers and subcycled time stepping receive architecture-specific implementations.
2.2. Particle Overloading
HACC uses particle overloading to buffer domain boundaries and reduce communication during force calculations. The approach improves locality and scaling but incurs memory and controllable short-range-force costs.
- Particle overloading extends each domain with a 3–10 Mpc passive-particle zone around its nominal boundary.
- The particle cache avoids communication for several time steps and keeps particle-mesh deposition and force interpolation local to each node.
- Short-range solvers require no communication infrastructure under overloading, so weak scaling depends on the long-range solver.
- The overloaded particle-mesh solver is formally exact, but short-range forces near passive-domain edges introduce controllable errors.
- Overloading costs memory efficiency and introduces boundary-force error, mitigated by balancing zone thickness with cache-refresh frequency.
2.3. The Long-Range Force
HACC computes long-range forces with a filtered spectral particle-mesh solver and a scalable two-dimensional pencil-decomposed FFT. Its filtering and spectral differentiation control force anisotropy while its data redistribution supports large-scale parallel execution.
- At separations above the force-matching scale, the filtered PM force approximates the Newtonian limit; Figure 3 tunes matching at three grid cells.The figure shows orientation-averaged force and its 1-σ anisotropy deviation.
- HACC combines CIC deposition, a Gausssinc spectral filter, fourth-order Super-Lanczos differentiation, and a sixth-order influence function for the long-range solver.The filter reduces CIC anisotropy noise and enables tunable force matching.
- The long-range force is computed from large FFTs, with the influence function representing the inverse discrete Laplacian.Each field-gradient component requires an independent FFT, but this is a small fraction of total force cost.
- HACC redistributes data between three-dimensional spatial domains and two-dimensional FFT pencils, overlapping transposes and FFTs to reduce communication hotspots.The data returns to a three-dimensional decomposition in the spectral domain.
- The parallel FFT has run on up to ∼1.5 million MPI ranks while communication efficiency remains unchanged as grid size increases on fixed processor counts.This behavior validates the implementation’s scaling design.
2.4. The Short-Range Force
HACC computes short-range forces by correcting the filtered PM force with a particle-based solver whose implementation is tuned for architecture and locality. RCB tree design and particle grouping improve both performance and, in clustered regions, accuracy.
- Below the force-matching scale rs, HACC subtracts the residual filtered grid force from the exact Newtonian force to obtain the short-range force.Above rs, only the long-range force is needed.
- A short matching scale enables a fifth-order polynomial force approximation that can be vectorized for high performance.This avoids the communication penalty of look-ups while simplifying force evaluation.
- HACC uses tree solvers on non-accelerated systems, while accelerated systems can use either tree or P3M methods.The production tree implementation pairs an RCB tree with a tuned polynomial force kernel.
- The RCB tree’s spatial locality lowers cache-miss rates because each leaf’s short-range calculation uses nearby leaf nodes.Particles are partitioned into disjoint memory buffers during recursive subdivision.
- RCB leaf nodes group particles to shift work from slow tree walking into the force kernel, with tens or hundreds of particles per leaf possible before the crossover.The benefit arises when reduced walk time exceeds the added kernel work.
- Keeping over 100 particles per leaf node can improve force accuracy by several orders of magnitude in highly clustered regions.More nearby particles are then handled exactly.
2.5. Time-Stepping
HACC advances particles with second-order symplectic maps, splitting long- and short-range dynamics when needed. Multi-level timestepping exploits clustering to accelerate evolution while preserving accuracy with small effects.
- HACC’s PM integrator composes stream and kick maps in symmetric SKS or KSK second-order symplectic steps.Stream updates positions at fixed velocity, while kick updates velocity using the force.
- For combined short- and long-range forces, HACC subcycles the short-range map because the long-range force varies relatively slowly.The Hamiltonian is split into Hsr and Hlr components.
- Typical problems use 3-10 short time steps for each long time step, depending on accuracy requirements.
- Multi-level timestepping uses density information from short-range solvers to provide speed-ups by a factor of 2-3 with only small effects on accuracy.The approach is motivated by highly clustered late-time particle distributions.
2.6. Code Units
HACC uses comoving coordinates and dimensionless internal variables, with the scale factor serving as the nominal evolution variable. Code units are tied to the long-range solver’s grid spacing and cosmological parameters.
- HACC represents positions and velocities in comoving coordinates and stores all variables internally in dimensionless form.
- The fundamental scaling length x0 is one long-range PM grid-cell length, L/(Ng −1), where L is box size and Ng is the one-dimensional grid count.
- HACC evolves using y = a^α with nominal α = 1 and performs time-stepping through the transformed derivative d/d˜t = (αyH/H0)d/dy.
2.7. Memory Management
HACC addresses memory fragmentation and variable particle-memory demands through optimized allocation and monitoring strategies designed for simulations operating near machine memory limits.
- Memory constraints: Large particle counts push simulations close to machine memory limits, making memory fragmentation a serious problem encountered on Blue Gene/Q systems.Fragmentation can leave insufficient contiguous memory despite substantial total free memory.
- Bigchunk allocator: Bigchunk mitigates fragmentation by allocating a large memory pool and distributing it among HACC subsystems.It initially tracks memory use by time-step phase, then allocates based on the maximum observed phase usage.
- Memory variation: Particle-memory requirements vary across MPI ranks because structure formation produces different particle counts in overloaded spatial sub-volumes.HACC records extrema of particle-memory usage and can adjust extra memory when restarting if initial estimates are insufficient.
- Memory-efficient reordering: An integer permutation array and scratch space enable out-of-place reordering one particle attribute at a time without additional allocation.The scratch array is sized for any single particle attribute.
2.8. I /O Strategy
HACC provides multiple I/O strategies because performance and scalability depend on use case, architecture, and data size. Its formats support scalable output, integrity validation, spatial queries, and interoperability with analysis tools.
- I/O strategies: HACC uses three main I/O approaches because no single strategy performs well across all use cases, architectures, and data sizes.The approaches include one file per process, many processes per file, and a single parallel-netCDF file.
- One file per process: One file per process provides the best write bandwidth but is practical only below roughly 10,000 processes and creates burdensome file management at larger scales.Reading with a different process count also requires redistribution of all data in memory.
- GenericIO: GenericIO reduces file-system contention by assigning processes disjoint file regions and reduces file counts by writing one file per I/O node.It supports MPI I/O in collective or non-collective mode and non-collective POSIX I/O.
- Data integrity: 64-bit CRC codes validate each variable for each rank and detect corruption during storage, transfer, or transmission through storage subsystems.A single-variable error was detected during a run producing 100 TB of checkpoint files.
- Performance: GenericIO achieves about two-thirds of peak performance under production conditions on IBM Blue Gene/Q.Tests recorded performance very close to the peak achievable before production-condition contention was considered.
- Parallel netCDF: Parallel netCDF organizes and indexes particles by spatial blocks, while supporting collective redistribution, block retrieval, and bounding-box queries.The netCDF format is integrated with many visualization and analysis tools and has an established user community.
- Parallel netCDF performance: 56% of ideal peak bandwidth is achieved using one shared file with collective I/O aggregation and indexed high-level formatting.The comparison uses an ideal 128-OST peak of approximately 26.7 GiB/s, while the tested strategy also stores index data.
3. Short-Range Force: Architecture Specific Implementations
HACC combines architecture-specific short-range solvers with adaptive multilevel time stepping and parallel tree strategies to sustain efficient, accurate force calculations across supercomputers.
- Architecture-specific implementations: HACC selects short-range solvers according to nodal architecture, including RCB trees for Blue Gene/Q and tuned kernels for accelerated systems.The implementation section emphasizes architecture-dependent algorithms and optimizations.
- Multilevel time stepping: A local density estimate from RCB-tree leaves assigns finer short-range time steps to higher-density particle regions.The kick count is controlled by an adjustable linear scaling parameter and a user-provided maximum level.
- Multilevel time stepping: A factor-of-two speed-up is achieved with multilevel stepping while retaining better-than-0.2% power-spectrum accuracy to kNy in the test case.At kNy/2, the reported error is approximately 0.05%.
- Accuracy control: Halo substructure varies only mildly as adaptive time-stepping parameters change, apart from differences attributed to FOF linking noise.The tests examine the largest and second-largest halos using tessellation-based density fields and spherical profiles.
- RCB-tree optimization: Multiple RCB trees built independently within fixed spatial bins improve force-computation performance and increase thread-level parallelism.Nearest-neighbor tree searches and distributed leaf-node/tree interaction pairs also support load balancing.
4. Code Verification and Testing
HACC verification tests show close agreement across implementations, reference codes, simulation scales, halo statistics, and matter power spectra.
- Code comparison: HACC tests compare multiple implementations with Gadget-2 using matter power spectra, halo structure, and halo statistics.The comparisons include PPTreePM, Cell-accelerated, and GPU-accelerated HACC versions.
- Code comparison: Better than 0.5% agreement across the full k-range is obtained among HACC versions and Gadget-2, while TreePM variants agree better than 0.1%.The comparison extends to the particle Nyquist wavenumber.
- Halo validation: Visual comparisons find that nearly all small halos and many fine details in larger halos are shared across HACC implementations and Gadget-2.The analysis includes halos with at least 100 particles and identifies missing dark-blue halos as below-threshold objects.
- Halo validation: Halo mass profiles agree very well beyond the small-radius binning shot-noise regime, with additional statistics reported in Table 3.The comparison covers the three largest halos and multiple algorithms.
- Large-scale simulations: Fractions-of-a-percent agreement is found between the Outer Rim and Q Continuum matter power spectra, and agreement with the Coyote emulator is at the 2% level.The emulator comparison is at its stated expected accuracy.
- Large-scale simulations: The halo mass function agrees across three same-cosmology simulations spanning approximately 30 billion to over one trillion particles and 1–4 Gpc boxes.The simulations use different short-range force implementations and resolutions.
5. In situ Analysis Tools
HACC includes in situ and stand-alone analysis tools that process simulation data during execution or afterward, reducing dependence on full raw-output storage.
- Framework: CosmoTools provides in situ and stand-alone analysis, with outputs available at run time or during post-processing.The framework can connect a ParaView server through Catalyst while the simulation runs.
- Analysis tools: The framework includes halo finders, merger-tree constructors, and parallel Voronoi tessellation for continuous density-field representations.The Voronoi mesh has cell volumes inversely proportional to particle spacing.
6. Selected Performance Results
HACC shows near-ideal weak scaling across Blue Gene/Q systems and Titan, while strong scaling remains effective over substantial core and node ranges.
- Weak scaling: The weak-scaling tests keep particles and physical volume per node or core fixed, so wall-clock time should remain constant as the system grows.Under this setup, time to solution per particle per step decreases while absolute TFlops per second rises.
- Weak scaling: 13.94 PFlops and 90% parallel efficiency were achieved at 1,572,864 cores on Sequoia with ideal weak scaling.The test used 2 million particles per core in an approximately (100 h^-1Mpc)^3 volume per core.
- Weak scaling: HACC weak-scales essentially perfectly to Titan’s full machine with 32 million particles per node in a fixed (256 h^-1Mpc)^3 nodal volume.The GPU version used one MPI rank per node.
- Strong scaling: Strong scaling remained effective from 512 to 16,384 Blue Gene/Q cores and from 32 to 8,192 Titan nodes, degrading only beyond specified points.Blue Gene/Q performance degraded somewhat at 16,384 cores, while Titan degradation began after 2,048 nodes.
- Strong scaling: HACC remains effective at approximately 100 MB per core, demonstrating operation under memory constraints comparable to next-generation many-core systems.The paper identifies reduced effective memory per core as the significance of the strong-scaling tests.
7. Conclusion and Outlook
HACC is positioned as a framework for synthetic catalogs, precision cosmological predictions, and large simulation campaigns supporting sky-survey analysis. Its future development targets both evolving architectures and expanded physical and analysis capabilities.
- Framework role: HACC is designed to produce synthetic catalogs and run large campaigns for precision predictions of cosmological observables.These uses support interpretation of survey data, survey-strategy optimization, and validation of data and analysis pipelines.
- Outlook: Future HACC development emphasizes architecture-oriented algorithms and optimization, including power management, fault tolerance, NVRAM, and alternative programming models.This is one of two stated directions for continued evolution.
- Outlook: A second development direction adds physics and analysis capabilities, including gas physics, feedback processes, more sophisticated semi-analytic galaxy modeling, and validation.These additions target both simulation modeling and downstream analysis.
- Framework role: HACC campaigns support cosmological emulators, covariance-matrix determination, survey-design and pipeline testing, and MCMC parameter estimation across multiple probes.These applications span prediction, uncertainty characterization, synthetic-catalog testing, and parameter inference.