PaperlayerView exact revision on arXiv

Source-linked AI summary

Simulating cosmic structure formation with the GADGET-4 code

Volker Springel, Rüdiger Pakmor, Oliver Zier, Martin Reinecke

arXiv:2010.03567v3astro-ph.IMastro-ph.CO

TL;DR

Cosmic-structure simulations need codes that remain accurate, scalable, and adaptable as problem sizes and hardware capabilities grow. This paper presents and evaluates GADGET-4, a publicly released, multipurpose update with new gravity, integration, hydrodynamics, parallelization, and analysis capabilities. The authors report robust converged test solutions and suitability for very large and high-dynamic-range simulations, while noting specific scope boundaries and test-distribution caveats.

  • Problem

    Cosmic-structure simulations require more accurate, larger, and higher-dynamic-range calculations, but these demands need more scalable and flexible numerical methods.

  • Method

    The paper describes and evaluates GADGET-4, combining alternative gravity solvers, hierarchical integration, hydrodynamics, parallelization, and on-the-fly structure-processing tools.

  • Results

    GADGET-4 provides converged, timestep- and force-accuracy-robust solutions for public test problems, with FMM offering momentum conservation and potential speed advantages in some situations.

  • Takeaways & Limitations

    The public release provides a reference code for accurate cosmic-structure simulations spanning very large problem sizes and very high dynamic ranges.

  • Takeaways & Limitations

    Stretched boxes support group finding and hydrodynamics, but secondary high-resolution PM meshes are not supported.

Abstract

from arXiv · show

Numerical methods have become a powerful tool for research in astrophysics, but their utility depends critically on the availability of suitable simulation codes. This calls for continuous efforts in code development, which is necessitated also by the rapidly evolving technology underlying today's computing hardware. Here we discuss recent methodological progress in the GADGET code, which has been widely applied in cosmic structure formation over the past two decades. The new version offers improvements in force accuracy, in time-stepping, in adaptivity to a large dynamic range in timescales, in computational efficiency, and in parallel scalability through a special MPI/shared-memory parallelization and communication strategy, and a more-sophisticated domain decomposition algorithm. A manifestly momentum conserving fast multipole method (FMM) can be employed as an alternative to the one-sided TreePM gravity solver introduced in earlier versions. Two different flavours of smoothed particle hydrodynamics, a classic entropy-conserving formulation and a pressure-based approach, are supported for dealing with gaseous flows. The code is able to cope with very large problem sizes, thus allowing accurate predictions for cosmic structure formation in support of future precision tests of cosmology, and at the same time is well adapted to high dynamic range zoom-calculations with extreme variability of the particle number density in the simulated volume. The GADGET-4 code is publicly released to the community and contains infrastructure for on-the-fly group and substructure finding and tracking, as well as merger tree building, a simple model for radiative cooling and star formation, a high dynamic range power spectrum estimator, and an initial conditions generator based on second-order Lagrangian perturbation theory.

1 INTRODUCTION

GADGET-4 updates a widely used cosmic-structure simulation code to address demands for larger, more accurate, scalable, and flexible calculations. The paper comprehensively describes and evaluates its numerical methods, implementation, and processing infrastructure.

  • Motivation: Growing scientific needs for larger statistical power, higher resolution, and broader dynamic ranges motivate the GADGET-4 update.These demands require more scalable integration schemes for multi-physics and multi-scale calculations.
  • Design goals: GADGET-4 targets improved scalability, larger simulations, high dynamic-range timestepping, force accuracy, and a more modular code architecture.The design also aims to help scientific users develop extensions more easily.
  • Design goals: The code prioritizes flexible multi-purpose use over optimization for a narrow application type.New methods include workload-balancing domain decomposition, hierarchical time integration, optional FMM, and multiple SPH formulations.
  • Infrastructure: GADGET-4 supports hybrid MPI/shared-memory parallelization and on-the-fly FOF, SUBFIND, merger-tree, and substructure-processing functionality.SUBFIND-HBT incorporates past subhalo information for robust and computationally efficient tracking.
  • Scope: The paper provides a comprehensive technical description and evaluation intended as a reference for GADGET-4 simulation methodology.The code had already seen initial applications in the literature.

2 GRAVITY CALCULATION

GADGET-4 computes gravitational interactions using softened, periodic potentials and several hierarchical force solvers. Its treatment accounts for periodic-image corrections, conservative tree opening, and accuracy limitations of interpolation-based corrections.

  • Softening: The adopted softening becomes significant near r ∼1.0ϵ0, although it vanishes completely only beyond the kernel’s 2.8ϵ0 outer edge.Thus ϵ0, rather than 2.8ϵ0, characterizes the effective softening scale.
  • Periodic potential: Periodic gravity is computed from the softened nearest-image interaction plus a correction representing the infinite periodic grid.The nearest-image interaction alone receives softening, while the correction remains unsoftened.
  • Periodic potential: 10^-5 relative force error already arises near 1 percent of the box size when the nearest-image Newtonian force replaces the full periodic force.The discrepancy grows toward order unity near half the box size and becomes non-isotropic.
  • Tree accuracy: Tree multipole expansions are rejected inside a geometric exclusion zone to protect against pathological force errors and Taylor-expansion breakdown.The exclusion test supplements the opening criteria used during top-down tree walks.
  • Periodic corrections: The new Taylor-expansion treatment of Ewald corrections reaches relative derivative-tensor accuracy of approximately 10^-10, unlike the lookup-table approach limited near 10^-3.The Taylor approach uses analytically differentiated Ewald expressions tabulated to higher order.

2.3 TreePM approach

TreePM splits periodic gravity into short-range tree forces and long-range Fourier forces. The split uses Ewald-inspired filtering, with accuracy and cost controlled by mesh and lookup-table choices.

  • Method: TreePM obtains the total potential as φ(x) = φsr(x) + φlr(x), evaluating short-range interactions with a tree and long-range interactions with Fourier methods.The method is a specialized Ewald decomposition using a small real-space split scale.
  • Short-range forces: A 256-entry one-dimensional lookup table supports maximum relative force errors as low as 10^-5, while 48 entries suffice for 10^-3 accuracy.Finer tables reduce interpolation errors but lower processor-cache efficiency.
  • Long-range forces: The long-range component is computed by particle-mesh density assignment, FFT transformation, Green’s-function multiplication, and inverse transformation.Forces are then obtained on the same mesh.
  • Long-range forces: The exponential Fourier cutoff provides a smooth scale beyond which long-range modes decline rapidly, supporting accurate mesh evaluation when the Nyquist frequency is sufficient.The mesh parameter Asmth controls the transition scale and the allocation of Fourier modes.
  • Accuracy-cost trade-off: Increasing Asmth improves long-range accuracy but moves the short-range cutoff to larger scales, while the fastest setting uses Asmth = 1.25.The Asmth trade-off makes optimum performance problem-dependent.

2.4 Fast multipole method

GADGET-4 adds a Cartesian FMM as an alternative to the one-sided tree, using symmetric source-sink expansions to conserve momentum and support multiple boundary and softening treatments. Its efficiency depends on the timestep workload.

  • Method: The Cartesian FMM is introduced as an alternative to the one-sided tree to improve hierarchical multipole calculations when high force accuracy is required.GADGET-4 extends the approach to distributed-memory parallelization, multiple softenings, periodic boundaries, and FMM-PM coupling.
  • Momentum conservation: Symmetric source-sink expansions yield manifest momentum conservation, with the vector sum of force errors equal to zero to machine precision.The ordinary tree algorithm does not generally provide this property and can accumulate substantial total-momentum errors.
  • Implementation: The method supports expansion orders p = 1 through p = 5, while some highest-order moments can be omitted because they contribute only constants to the potential.The omitted terms therefore do not affect particle forces.
  • Algorithm: FMM constructs node multipoles recursively and evaluates expansions for interacting node pairs through a dual tree walk.The interaction uses expansion centers at the nodes’ centers of mass, where dipole moments vanish.
  • Performance boundary: FMM is poorly matched to partial force calculations because computing forces for 5% of sinks still requires effort close to a full calculation.This disadvantage disappears with hierarchical time integration, where source and sink sets remain equal.

2.5 The FMM-PM approach

GADGET-4 extends TreePM to an FMM-PM formulation by retaining the mesh treatment for long-range forces and replacing the short-range tree calculation with FMM. The approach uses a modified Green’s function and cutoff-based interaction pruning while preserving the stated force-accuracy framework.

  • FMM-PM formulation: FMM-PM retains the PM calculation and replaces the TreePM short-range force evaluation with the fast multipole method.The formulation changes the real-space Green’s function from the Newtonian form to the corresponding short-range force kernel.
  • Interaction control: Interactions are omitted when the nearest sides of two nodes lie beyond the cutoff distance rcut, preventing further refinement of those interactions.The opening criteria remain unchanged, and relative force errors are compared with total gravitational accelerations including PM contributions.
  • Mesh placement: The secondary PM mesh is positioned by measuring the high-resolution region’s diameter, selecting the next larger power-of-two subdivision, and shrink-wrapping touched cells into a cubical footprint.This construction aligns the secondary mesh with a level of the global oct-tree.

2.6 Accelerating short-range force calculations through a secondary mesh

GADGET-4 accelerates extreme zoom simulations with a secondary high-resolution PM mesh whose oct-tree-aligned region enables a refined, momentum-conserving force split. The method reduces high-resolution tree work but introduces additional PM computation and has parameter-dependent benefits.

  • Motivation: A single full-box PM grid provides limited speed-up in zoom simulations because high-resolution particles can occupy only a few coarse-grid cells.The problem arises from the extreme particle-density variation within the simulation domain.
  • Region construction: GADGET-4 enlarges the particle-defined high-resolution region to align with oct-tree node boundaries, producing a cubical region for a clean force split.Nodes at the corresponding tree level and below either lie fully inside or fully outside the high-resolution region.
  • Secondary mesh: The secondary mesh computes intermediate-scale forces for mass inside the enlarged region using a zero-padded FFT and a difference of coarse- and fine-mesh short-range kernels.Zero padding realizes vacuum boundary conditions for this intermediate calculation.
  • Force split: High-resolution particles combine background-PM and intermediate-mesh forces with tree or FMM forces using rHR cut for internal partners and rLR cut for external partners.Particles outside the high-resolution node retain the ordinary background PM plus extended-cutoff short-range calculation.
  • Conservation: The refined procedure preserves manifest force antisymmetry when FMM is used because the PM algorithm also supplies antisymmetric forces.This addresses the asymmetry present in the earlier approach.
  • Trade-offs: The secondary mesh can shorten high-resolution tree walks, but its extra zero-padded PM calculation may be more expensive at a given grid resolution.Its benefit is therefore conditional rather than guaranteed; the stated criterion is only a rough indication of when it may help.
  • Gravitational softening: For particle pairs with individual softenings, GADGET-4 uses the larger softening length, preserving force antisymmetry and ensuring each particle’s interaction is softened at least by its assigned value.Tree nodes store their largest constituent softening to guide multipole treatment.

2.8 Periodicity only in two dimensions

GADGET-4 supports gravity with periodicity in only two dimensions and extends its algorithms to stretched periodic boxes. These capabilities broaden boundary-condition support, but stretched boxes impose geometric and mesh constraints and do not support secondary high-resolution meshes.

  • Mixed periodicity: GADGET-4 allows gravity to be periodic in two spatial dimensions while remaining non-periodic in the third, targeting stratified or shearing-box simulations.The implementation derives and evaluates the corresponding mixed-boundary potential through an Ewald-based correction.
  • Force calculation: The mixed-boundary implementation treats the nearest periodic images in the two periodic directions with the ordinary kernel and adds a Taylor-expanded correction from a lookup table.Tree, TreePM, and FMM-PM methods are supported for these boundary conditions.
  • Stretched boxes: GADGET-4 supports stretched periodic domains for one-sided Tree, FMM, TreePM, and FMM-PM gravity calculations.The domain may be stretched by different factors along each spatial dimension.
  • Constraints: Stretched-box calculations constrain stretch factors so fully filled tree nodes remain cubical and each dimension contains an integer number of PM cells.Group finding and hydrodynamics remain operational, but secondary high-resolution PM meshes are unavailable.
  • Coordinate representation: Integer coordinate storage provides uniform positional resolution and supports branch-free periodic mapping, fast oct-tree construction, and lossless Peano-Hilbert key generation.GADGET-4 internally uses 32-bit, 64-bit, or 128-bit integers for coordinates before converting differences for force calculations.

3 FORCE ACCURACY TESTS

The tests show that FMM force errors depend strongly on node-size symmetry and lack the cancellation boost seen in Tree calculations. Relative opening criteria improve efficiency and provide accuracy-dependent choices of expansion order, while high-redshift and PM-coupled regimes impose important accuracy limits.

  • FMM force accuracy between two interacting nodes: Equal-sized interacting nodes produce symmetric and lowest overall FMM force errors, whereas asymmetric nodes give the smaller node smaller relative errors.This motivates opening the larger node when necessary so interacting nodes have equal sizes.
  • FMM force accuracy between two interacting nodes: For equal expansion order and opening angle, Tree errors are lower than FMM errors, although both scale with the same power-law order.FMM incurs an additional sink-side expansion, producing an approximately constant error offset of ≃4 relative to Tree in the tested configuration.
  • Force errors for an isolated particle distribution: Tree errors gain an accuracy boost from cancellation of uncorrelated partial-force errors, whereas FMM sink-side errors remain correlated and show little such cancellation.The Tree scaling steepens by about ∼1.0 for p = 2 and slightly less for p = 3 in the tested particle distribution.
  • Force errors for an isolated particle distribution: The relative opening criterion is more efficient than the geometric criterion, and the fastest algorithm depends on the target accuracy: FMM-O5 below 10−4, FMM-O4 from 10−4 to 10−3, and FMM-O3 from 10−3 to 10−2.FMM-O2 becomes inefficient at high accuracy, while higher expansion orders become worthwhile when stricter accuracy is required.
  • Cosmological force-accuracy tests: At high redshift, force errors for a given opening angle are larger, and Tree-O1 and FMM-O1 without a PM mesh are essentially unusable.The high-redshift behavior reflects the small peculiar potential and peculiar forces in the nearly unperturbed density field.

4 TIME INTEGRATION OF COLLISIONLESS PARTICLES

GADGET-4 combines leapfrog-based cosmological integration with hierarchical timestepping and multiple gravity solvers to improve accuracy, momentum conservation, and efficiency across large timescale ranges. Performance depends strongly on timestep hierarchy depth, solver choice, accuracy demands, and computing environment.

  • Integration scheme: GADGET-4 advances collisionless particles with a kick-drift-kick leapfrog scheme using logarithmic scale factor time τ = ln a.A constant step in τ corresponds to a fixed fraction of the current Hubble time.
  • Integration scheme: The unmerged half-step kicks preserve second-order accuracy for output velocities and simplify nesting hydrodynamics within gravity steps.Only one force calculation per timestep is required because the final force can be reused for the next step’s first kick.
  • Hierarchical time integration: Hierarchical integration removes asymmetric active-particle coupling and restricts tree construction to particles participating in the current force calculation.The resulting local operators have equal source and sink sets, giving manifest momentum conservation.
  • Hierarchical time integration: For very deep timestep hierarchies, hierarchical integration can outperform ordinary integration by avoiding full-system tree construction and passive-particle drift overhead.Its relative cost depends on timestep-bin occupancy; the hierarchy may also improve accuracy at equal timestep size.
  • Gravity solver choices: GADGET-4 supports time-unsplit Tree-PM/FMM-PM calculations that choose between PM-accelerated and pure Tree/FMM forces on the fly.The most efficient configuration remains problem-dependent, and explicit tests are needed for scientifically relevant setups.

5 HYDRODYNAMICAL DISCRETIZATION

GADGET-4 modernizes smoothed particle hydrodynamics while retaining its adaptive, robust particle-based formulation. It supports entropy-conserving and pressure-based variants, improved kernels and viscosity treatment, and basic cooling and star-formation modules.

  • SPH motivation: SPH provides automatic spatial adaptivity, robustness, and effective treatment of nearly empty space, but standard formulations can be noisy and suppress fluid instabilities.These limitations motivate ongoing formulation improvements.
  • SPH formulations: GADGET-4 supports both conservative entropy-based vanilla SPH and pressure-based P-SPH, which behaves better at contact discontinuities.The pressure-based formulation can improve interface behavior at the cost of somewhat higher noise elsewhere.
  • SPH formulations: The entropy-conserving formulation preserves energy and entropy but introduces spurious surface tension that suppresses or slows instability growth at phase boundaries.P-SPH addresses this issue through a pressure-based variational formulation.
  • Kernel choices: Wendland kernels are included as alternatives because they prevent clumping instability and permit larger neighbour numbers for reducing interpolation errors.Recommended neighbour counts range from 92 for C2 to 356 for C6.
  • Artificial viscosity: The artificial-viscosity treatment combines bulk and von Neumann-Richtmyer viscosity, with typical α values of approximately 0.5–1.0 and β = 2α.The formulation compensates viscous work through entropy generation to conserve total energy.
  • Additional physics: GADGET-4 includes basic radiative cooling for hydrogen and helium in collisional ionization equilibrium and star formation based on a coarse-grained interstellar-medium model.These modules are intended to help users begin galaxy-formation simulations rather than provide comprehensive subgrid physics.

6 PARALLELIZATION STRATEGY

GADGET-4 uses spatial domain decomposition, hybrid MPI/shared-memory parallelism, and buffered tree communication to improve scalability while preserving computational semantics as far as possible. The design targets large simulations and many-core systems, but some strategies introduce trade-offs in determinism, portability, or speed.

  • 6.1 Domain decomposition: GADGET-4 subdivides the simulated volume into disjoint regions mapped to MPI ranks to balance computational and memory loads.The decomposition is designed so that domain boundaries do not alter tree geometry, force errors, or the multipole expansions seen by particles.
  • 6.2 Hybrid parallelization: Many-core systems make decomposition, load balancing, and MPI scalability increasingly difficult as a fixed problem is distributed over more ranks.These challenges motivate shared-memory communication within nodes and more scalable parallel strategies.
  • 6.2 Hybrid parallelization: MPI-3 shared-memory allocation lets ranks on the same node replace intra-node message passing with direct memory reads and writes.GADGET-4 assigns a communication rank to fetch remote-node tree data, buffers requests into larger packages, and shares imported data among local ranks.
  • 6.3 Data layout and communication scheme in parallel tree walks: Buffered tree communication can break binary floating-point determinism because imported nodes may arrive in different orders and change accumulation order.The approach nevertheless reduces repeated tree-walk work by retaining starting-node indices for partial walks on foreign processors.
  • 6.4 Vectorization: The AVX-based SPH kernel achieved only a 5-10% speed-up over the ordinary serial routine on Intel Xeon Gold 6138 CPUs.The implementation processes four neighbour particles in parallel using 256-bit AVX instructions, while the serial routine remains the default.

7 ON THE FLY ANALYSIS AND OTHER FEATURES

GADGET-4 integrates scalable on-the-fly group, substructure, merger-tree, and light-cone analysis tools for large simulations and high-resolution zoom calculations. Its SUBFIND-HBT variant uses temporal membership information to identify and track substructures, while additional outputs reduce storage demands and support downstream mock-observation applications.

  • Group and substructure finding: GADGET-4 includes FOF, SUBFIND, and SUBFIND-HBT finders capable of processing very large simulations and high-resolution zoom calculations.SUBFIND-HBT uses past group and subhalo membership information, unlike the positional-space SUBFIND workflow.
  • Group and substructure finding: SUBFIND identifies bound substructures by locating density peaks, growing them to saddle points, and applying gravitational unbinding.SUBFIND-HBT can optionally use previous subhalo membership when assigning substructure and background peaks.
  • Temporal substructure tracking: SUBFIND-HBT processes snapshots in temporal order, uses prior memberships to form candidates, and unbinds non-background candidates at their new phase-space coordinates.By construction, the method can recover at most the mass a substructure had in the previous output before subsequent stripping.
  • Temporal substructure tracking: SUBFIND-HBT is simpler and computationally cheaper than SUBFIND because it avoids adaptive density estimation, serial saddle-point processing, and many costly unbinding operations.The supplied passage identifies these avoided steps as the source of the computational difference.
  • Merger trees: Merger trees are built on the fly from group and subhalo catalogues linked through retained particle-membership information, without storing full particle data.This enables finely time-resolved merger trees while keeping catalogue data volumes small.
  • Light-cone and auxiliary outputs: GADGET-4 can generate light-cone outputs, identify halos directly on the lightcone, and save halos instead of particles to substantially reduce storage volume.The resulting outputs support applications such as HOD mocks and gravitational lensing without discontinuities from patched time slices.

8 CODE VALIDATION AND CONVERGENCE TESTS

The validation tests show that GADGET-4 can achieve converged, robust results with sufficiently accurate force calculations and time integration, while remaining sensitive to under-resolved settings and cross-code systematics.

  • Force accuracy: Force errors can be pushed below 10^-2 without large outliers, supporting the in-situ validation of parallel force computations.Random translations decorrelate force errors in time, and the authors expect the resulting errors to be sufficiently small.
  • Time integration: A timestep of Δln(a) = 5.6 × 10^-4 or smaller converges the power spectrum to better than 1% down to the softening scale.This requires at least 8192 steps from z = 99 to z = 0 in the tested 25 h−1Mpc box.
  • Time integration: Using only 64 steps can limit fundamental-mode deviations below 1%, but is inadequate in mildly non-linear and highly non-linear regimes.The required timestep therefore depends strongly on the physical scale and degree of non-linearity being studied.
  • Time integration: Density profiles remain indistinguishable from finely timed runs for Δln(a) ≤ 5.6 × 10^-4, whereas twice larger steps systematically reduce density near 10 h−1kpc.Insufficient time-integration accuracy generically lowers central halo density.
  • Robustness: With sufficiently accurate force and time integration, results are consistent across the supported force algorithms and time-stepping schemes, though this offers no guarantee against poor settings.Force-accuracy reductions remain relatively stable when randomizations are applied.
  • Cross-code comparison: Across codes, power spectra agree reasonably for k < 10 h Mpc−1 but diverge strongly in the highly non-linear regime, while halo profiles differ by up to 30% near the softening scale.The authors associate these systematic differences with factors including force softening and time-stepping accuracy.
  • Hydrodynamics: The Wendland C6 kernel with 66 neighbours achieves L1 ∝N^-0.9, close to the theoretical first-order optimum for the tested discontinuous flow.Convergence depends mainly on the kernel and indirectly on neighbour number; the two SPH formulations have similar convergence behaviour and runtime effects.

9 EXEMPLARY APPLICATIONS

GADGET-4 is demonstrated on lightcone catalogues, cluster hydrodynamics, galaxy mergers, and mixed-boundary self-gravitating sheets. These examples test observational outputs, SPH behaviour, subgrid physics, and gravity-solver accuracy.

  • Lightcone output: A 1 h−1Gpc ΛCDM box produces full-sky halo catalogues from z = 0.4 to z = 0 while outputting particles only in a 5 h−1Mpc-thick lightcone disk.The initial conditions are generated at z = 99, and group catalogues are produced without saving the full particle lightcone.
  • Lightcone output: On-the-fly lightcone group finding supports observable halo mass functions representing absolute sky counts above mass thresholds to specified look-back redshifts.The catalogues correspond to objects that would be observationally accessible under ideal capabilities.
  • Santa Barbara cluster: The Santa Barbara cluster shows resolution-consistent dark-matter and gas profiles, while central entropy remains method-dependent across hydrodynamical codes.The pressure-based and density-based SPH formulations give very similar results in this test.
  • Limitations: The Santa Barbara validation remains limited because the correct central entropy profile is not fully established and SPH can suppress turbulent mixing.The suppression is slightly more pronounced in the density-entropy formulation, while some mesh-code entropy cores may be excessively high.
  • Galaxy merger: The merger sequence shows tidal arms, gas inflow and outflow, dynamical-friction-driven re-encounter, disk destruction, a moderate starburst, and shock-heated hot gas.The calculation uses hierarchical time integration, with results essentially invariant to the integration settings tested.
  • Galaxy merger: The merger consumes gas on a timescale of about 1 Gyr, with star formation enhanced at first encounter and rapidly altered by the later encounter.The star-formation history is modulated as the galaxies interact, linger near turn-around radii, and eventually merge.
  • Mixed-boundary gravity: With periodicity in x and y but open z boundaries, perturbations dissipate in the thick sheet and grow into spherical fragments in the thin sheet.The contrasting evolution matches the expected stability of the two configurations.
  • Mixed-boundary gravity: The thick sheet reaches Spitzer’s analytic isothermal-sheet profile, indirectly validating the gravity solver under mixed boundary conditions.The comparison is made between the final SPH particle densities and the analytic solution.

10 CODE PERFORMANCE AND SCALABILITY

GADGET-4 targets flexible performance across uniform, hydrodynamic, cosmological, and extreme zoom workloads. Scaling is generally strong, but efficiency declines when per-core work becomes too small or domain decomposition dominates.

  • Design goal: The code is designed as a flexible multi-purpose system prioritizing reliability and adaptability over specialization for a single hardware or problem class.The performance tests therefore include realistic cosmological and extreme zoom simulations rather than only uniform benchmark loads.
  • SPH scaling: Strong scalability is generally good for the SPH blob test until the load per core falls to a few thousand particles.At that point, domain decomposition costs become large relative to the work distributed; larger problems extend the scalable core count.
  • Weak scaling: Weak scaling reaches about 66% parallel efficiency on 10^5 cores relative to one core, corresponding to a speed-up by a factor of 66000.The test keeps the particle load per core fixed while increasing the particle, mesh, and partition sizes together.
  • Cosmological scaling: Homogeneous cosmological runs test strong scaling from one core to 10^4 cores while the particle load per core decreases from over 134 million to little more than 10000.The tests cover both high-redshift initial conditions and the highly clustered z = 0 state.
  • Force-solver scaling: TreePM and FMM-PM variants show similar scalability, although absolute and relative speeds depend on expansion order, force accuracy, and other settings.The production-oriented tests avoid extremely small opening angles that could improve formal scaling while worsening total runtime.
  • Hydrodynamic scaling: The Santa Barbara hydrodynamical simulation scales well through 1280 cores, with further speed-up beginning to degrade as more nodes are added.The test evolves the clustered z ∼ 0 state across full timestep-hierarchy cycles.
  • Zoom scaling: The demanding Milky Way zoom achieves satisfactory strong-scaling efficiency from 20 to 1280 cores, although Tree-PM variants scale somewhat better than the tested FMM-PM-O3 setup.The FMM-PM-O3 test uses hierarchical time integration and is included to broaden the performance coverage, not because it scales best.

11 PUBLIC CODE RELEASE

GADGET-4 is publicly released through a collaborative git-based platform intended to improve maintainability, contribution, and reproducibility. The code also adds safeguards and architectural changes for memory management and modularity.

  • Motivation: The release responds to fragmented GADGET-3 development, which complicated validation and reproducibility while duplicating effort.Multiple independently developed versions created a difficult-to-navigate code landscape.
  • Collaborative release: The git-based public release preserves commit history and supports moderated patches, pull requests, branching, and sharing of private developments.This replaces static source archives with a collaborative development workflow.
  • Code maintenance: New compile-time features and run-time parameters must be documented or the code refuses to compile.A build-time control script is used to guard against undocumented additions.
  • Memory management: GADGET-4’s memory system prevents leaks and fragmentation, tracks allocation responsibility, and protects nodes from swapping.It also reports precise memory usage and allocation sources.
  • Architecture: Moving to C++ organizes the code into purpose-specific classes and uses inheritance and templating to improve modularity and reduce duplication.The transition modifies older C code rather than redesigning the system entirely.

12 DISCUSSION AND CONCLUSIONS

GADGET-4 is presented as a publicly released reference code for accurate, versatile, robust, scalable cosmic-structure simulations. Its updates span gravity, time integration, load balancing, hydrodynamics, large-scale communication, and analysis infrastructure, with tests showing converged solutions robust to timestep and force-accuracy changes.

  • Overall goals: GADGET-4 targets accuracy, versatility, robustness, scalability, very large problem sizes, and high-dynamic-range simulations.The authors position it as a reference code for numerical cosmology.
  • Gravity and integration: FMM provides a momentum-conserving alternative to tree gravity, with higher-order and hybrid TreePM or FMM-PM modes; momentum conservation can persist under adaptive local timestepping.The latter property follows from combining FMM with hierarchical gravity integration.
  • Scalability: The revised domain decomposition balances multiple constraints while avoiding memory imbalance, although decomposition cost can limit strong scalability for small problems.A future hierarchical treatment of shared-memory nodes is expected to reduce this cost further.
  • Hydrodynamics: The code supports improved vanilla SPH and pressure-based SPH, hybrid MPI/shared-memory communication, and separate hydrodynamic and gravitational timesteps.Hydrodynamics can subsample gravitational dynamics, and simple cooling and star-formation prescriptions are included.
  • Scalability: Scalability barriers are removed through flexible column-based distributed FFTs, arbitrarily stretched boxes, restricted gravity periodicity, and support beyond former billion-object and message-size limits.The code divides oversized communications across several calls when needed.
  • Analysis and usability: Built-in merger trees, lightcone outputs, lightcone halo identification, power-spectrum measurement, performance metrics, and HDF5 output improve analysis and usability.These facilities support rapid diagnosis and analysis of simulation results.
  • Validation: Test applications produce converged solutions for public problems that remain robust when timestepping and force accuracy change.The authors also report independently demonstrated correct and unbiased forces.

APPENDIX A: MULTIPOLE EXTENSIONS AT DIFFERENT ORDER

Appendix A specifies the multipole extensions used in the FMM potential approximation at different expansion orders. It lists higher-order kernel derivatives and Newtonian forms for the interaction coefficients.

  • Potential approximation: The appendix introduces expressions for the FMM potential approximation at different multipole expansion orders p.These expressions define the successive terms used by the approximation.
  • Expansion orders: The displayed expansions add progressively higher multipole terms involving Q2·D2 through Q5·D5 and powers of the displacement vector a.The sequence includes coefficients such as 1/6, 1/24, and 1/120 for higher-order contributions.
  • Kernel specification: Explicit derivatives of g(r) are supplied for the interaction kernel used in equations (29)–(34).The appendix thereby makes the derivative inputs to the multipole expressions explicit.
  • Newtonian case: For Newtonian gravity, the derivatives are given as alternating signed inverse powers from g0 = 1/r through g7 = −135135/r15.The sequence progresses through odd inverse powers of r.
  • Kernel derivatives: Expressions for D6 and D7 provide higher-order derivatives of the interaction kernel beyond the derivatives in equations (29)–(34).The displayed index contractions include extended Kronecker-delta combinations.
Loading 2010.03567v3…