SWIFT: A modern highly-parallel gravity and smoothed particle hydrodynamics solver for astrophysical and cosmological applications

TL;DR

SWIFT addresses the challenge of simulating coupled gravity and hydrodynamics for large-scale structure, galaxy formation, and planetary contexts on modern HPC systems. It advances a hydrodynamics-first, task-based architecture that coalesces SPH-based hydro with a Fast Multipole Method gravity solver (optionally coupled to a PM grid) and supports delta-f neutrino methods, multiple sub-grid galaxy formation models, and extensive I/O/analysis tooling. The code demonstrates strong weak and strong scaling to billions of particles and delivers a versatile public release with documentation, test problems, and analysis utilities, marking a significant step toward reproducible, large multi-physics simulations. By combining modular SPH flavours, adaptive time stepping, proxy-domain communications, and flexible cosmological integration, SWIFT enables high-fidelity, scalable simulations with broad scientific reach and practical impact for cosmology and galaxy evolution research.

Abstract

Numerical simulations have become one of the key tools used by theorists in all the fields of astrophysics and cosmology. The development of modern tools that target the largest existing computing systems and exploit state-of-the-art numerical methods and algorithms is thus crucial. In this paper, we introduce the fully open-source highly-parallel, versatile, and modular coupled hydrodynamics, gravity, cosmology, and galaxy-formation code SWIFT. The software package exploits hybrid shared- and distributed-memory task-based parallelism, asynchronous communications, and domain-decomposition algorithms based on balancing the workload, rather than the data, to efficiently exploit modern high-performance computing cluster architectures. Gravity is solved for using a fast-multipole-method, optionally coupled to a particle mesh solver in Fourier space to handle periodic volumes. For gas evolution, multiple modern flavours of Smoothed Particle Hydrodynamics are implemented. SWIFT also evolves neutrinos using a state-of-the-art particle-based method. Two complementary networks of sub-grid models for galaxy formation as well as extensions to simulate planetary physics are also released as part of the code. An extensive set of output options, including snapshots, light-cones, power spectra, and a coupling to structure finders are also included. We describe the overall code architecture, summarise the consistency and accuracy tests that were performed, and demonstrate the excellent weak-scaling performance of the code using a representative cosmological hydrodynamical problem with $\approx$$300$ billion particles. The code is released to the community alongside extensive documentation for both users and developers, a large selection of example test problems, and a suite of tools to aid in the analysis of large simulations run with SWIFT.

Figures (20)

• Figure 1: A selection of simulation results obtained with the Swift code, illustrating the huge range of problems that have already been targeted and the flexibility of the solver. The panels show: (a) a projection of the large-scale distribution of dark matter from a $10~{\rm Mpc}/h$ slice of the $(500~{\rm Mpc}/h)^3$ benchmark simulation of Schneider2016; (b) the temperature of the gas weighted by its velocity dispersion in a zoom-in simulation of a galaxy cluster using the Swift- Eagle galaxy formation model (§ \ref{['ssec:eagle']}) extracted from the runs of Altamura2022; (c) an idealised isolated galaxy from the Agora-suite Kim2016 simulated using the Gear model (§ \ref{['ssec:gear']}) rendered using pNbodypNbody; and (d) a snapshot extracted from a Moon-forming giant impact simulation of Kegerreis2022 using the planetary physics extension of the code (§ \ref{['ssec:planetary']}) and rendered using the Houdini software.
• Figure 2: The Verlet-list method. By constructing a mesh structure with cell sizes matching the search radius $H$ of particles, the neighbour-finding strategy is entirely set by the geometry of the cells and the list of potential candidates is thus exactly known. The particle in black only has potential neighbours in the cell where it resides or any of the 8 (26 in 3D) directly neighbouring cells (in grey). The smoothly varying nature of SPH leads to particles having similar $H$ in nearby regions, with this scale only varying slowly over the whole simulated domain.
• Figure 3: An example of interactions between regions of different densities, i.e. particles with different search radii. Particle $a$ will interact with the particles on the left and above using the smaller cells. It will interact with the particles on the right using the larger cell. The particle $b$ will only interact using the cells at the coarser level. Thanks to the nested grids, interactions happen at different levels in the hierarchy depending on the local search radius. Once the grid is constructed, all the possible interactions at the different levels are known without the need of a speculative tree-walk.
• Figure 4: A simplified graph of the tasks acting on a given cell for SPH and gravity during one time step in Swift. Dependencies are depicted as arrows and conflicts by dotted lines. Once the particles have been drifted to the current point in time, the first loop over neighbours can be run. The so-called "ghost" task serves mainly to reduce the number of dependencies between successive loops over the neighbours. Once the second loop has run, the time integration (§ \ref{['ssec:design:multi_dt']}) can be performed. In parallel to the SPH operations, the gravity tasks (condensed into a single one here for clarity) can be run as they act on different subsets of the data. To prevent different threads from over-writing each others' data, the various SPH loop tasks (1 self and 26 pairs) are prevented from running concurrently via our conflict mechanism. Additional loops over neighbours, used for instance in more advanced SPH implementations, in sub-grid models or for radiative transfer, can be added by repeating the same pattern. They can also be placed after the time integration tasks if they correspond to terms entering the equations in an operator splitting way.
• Figure 5: The execution of various tasks using 8 threads over the course of one time-step, extracted from a cosmological hydrodynamical simulation with $2\times128^3$ particles using only gravity and hydrodynamics on a shared-memory system. The different rows correspond to the different threads on the compute node. The work each thread performs is coloured to correspond to the task type it executes. Yellow, for instance, corresponds to a self-task performing gravity operation on a cell, whereas navy blue corresponds to a pair-task performing a 3 SPH loop over two cells. Note that some tasks displayed in the legend do not actually run in this example. For instance, no MPI-related send or recv tasks are executed here. We show them in the legend for consistency with Fig. \ref{['fig:design:task_graph_mpi']}. The long bands are actually a series of the same task acting on different cells one after the others. There are for instance 512 yellow tasks. As desired, the threads display essentially no idle time (white gaps) between operations and all end their work at very nearly the same time. In other words, the load balancing is near-perfect with no parallel performance loss. The small gap at the start corresponds to cost of deciding what tasks to activate for this step. Bands of a given colour can have different lengths, indicating that tasks can correspond to very different workloads depending on how many particles are present in the cell(s) on which they act. At a given point in time, different threads often process different task types, and hence solve a different set of equations. This is different from the traditional branch-and-bound parallelism approach where all threads perform the same action and have to wait until they have all completed it before moving to the next piece of physics.