PaperlayerView exact revision on arXiv

Source-linked AI summary

NBODY6++GPU: Ready for the gravitational million-body problem

Long Wang, Rainer Spurzem, Sverre Aarseth, Keigo Nitadori, Peter Berczik, M. B. N. Kouwenhoven, Thorsten Naab

arXiv:1504.03687v2astro-ph.IMastro-ph.SR

TL;DR

Dense and collisional star clusters require direct pairwise gravitational interactions, motivating larger and more realistic simulations. The paper presents NBODY6++GPU, combining MPI with GPU and AVX/SSE acceleration, and demonstrates million-particle globular-cluster simulations with stellar evolution and primordial binaries.

  • Problem

    Dense and collisional star clusters require pairwise gravitational interactions between many or all stars to integrate their motions accurately.

  • Method

    NBODY6++GPU combines the MPI-parallelized NBODY6++ code with GPU and AVX/SSE acceleration on each node.

  • Results

    400-2000 times faster than NBODY6, NBODY6++GPU enables million-particle globular-cluster simulations with 5% primordial binaries at about one hour per half-mass crossing time.

  • Takeaways & Limitations

    A million-body cluster can be simulated for about 20 crossing times in one day on 320 CPU cores with 32 GPUs.

  • Takeaways & Limitations

    NBODY6++GPU is not suitable for clusters that reserve most CPU cores for pure CPU tasks because it relies heavily on CPU calculations and data movement.

Abstract

from arXiv · show

Accurate direct $N$-body simulations help to obtain detailed information about the dynamical evolution of star clusters. They also enable comparisons with analytical models and Fokker-Planck or Monte-Carlo methods. NBODY6 is a well-known direct $N$-body code for star clusters, and NBODY6++ is the extended version designed for large particle number simulations by supercomputers. We present NBODY6++GPU, an optimized version of NBODY6++ with hybrid parallelization methods (MPI, GPU, OpenMP, and AVX/SSE) to accelerate large direct $N$-body simulations, and in particular to solve the million-body problem. We discuss the new features of the NBODY6++GPU code, benchmarks, as well as the first results from a simulation of a realistic globular cluster initially containing a million particles. For million-body simulations, NBODY6++GPU is $400-2000$ times faster than NBODY6 with 320 CPU cores and 32 NVIDIA K20X GPUs. With this computing cluster specification, the simulations of million-body globular clusters including $5\%$ primordial binaries require about an hour per half-mass crossing time.

1 INTRODUCTION

Direct N-body simulations model dense, collisional star clusters accurately but become prohibitively expensive as particle numbers grow. NBODY6++ extends NBODY6 for parallel supercomputers, motivating further GPU- and multi-node acceleration.

  • Motivation: Dense, collisional clusters require pairwise gravitational interactions among many or all stars for correct stellar-motion integration.NBODY6 is designed for this regime, where close encounters and stellar collisions may occur.
  • Computational challenge: O(N^7/3) is the computational cost per crossing time for homogeneous systems using individual time steps.The number of steps per particle scales roughly as N^1/3 under the stated Hermite-scheme assumptions.
  • Computational challenge: O(N^10/3 / ln N) is the cost per half-mass relaxation time for homogeneous systems.This scaling follows from the relation between half-mass relaxation time and crossing time.
  • Prior acceleration: NBODY6++ extended NBODY6 with MPI parallelization for general parallel supercomputers, while GPU and SSE/AVX methods enabled larger simulations.Earlier GPU-based work made simulations with N ∼ 10^5 possible on a single workstation.
  • Prior acceleration: Previous direct simulations reached 200,000 and 262,500 stars, establishing the scale immediately preceding million-body modeling.The cited studies included thousands of primordial binaries in globular-cluster simulations.

2 THE FEATURES OF NBODY6/6++

NBODY6 combines fourth-order Hermite integration, hierarchical block time steps, and the Ahmad-Cohen neighbor scheme to reduce integration and force-calculation overheads. These efficiencies become constrained by parallel implementation complexity and neighbor-related overheads.

  • Integration: NBODY6 uses fourth-order Hermite integration with hierarchical block time steps that integrate active particles together.Time steps are quantized, usually as integer powers of 0.5.
  • Force calculation: The Ahmad-Cohen scheme separates regular long-range forces from irregular neighbor forces using per-particle neighbor lists.Regular forces are evaluated at larger time steps, while neighbor forces are accumulated during irregular steps.
  • Force calculation: N^1/4 is the approximate speed gain from the Ahmad-Cohen scheme in sequential computing.The scheme's benefit is associated with separating regular and irregular force calculations.
  • Parallel limitations: Parallel computing limits the Ahmad-Cohen gain because implementing the algorithm is complex.The block-step reduction in particle-prediction overhead is also strongly limited in the neighbor scheme.
  • Integration: Direct N-body integration requires careful treatment because energy conservation accuracy is important for self-gravitating systems.The cited three-body analysis found a requirement of better than 1/10th total-energy conservation, though its extension to large-N systems is uncertain.

3.1 MPI parallelization of NBODY6++

NBODY6++ uses MPI with a copy algorithm that distributes active particles across processors while replicating the complete particle dataset on each processor. A ring algorithm reduces memory use but requires additional neighbor-data communication.

  • Copy algorithm: MPI NBODY6++ parallelizes both regular and irregular forces by assigning different subsets of active particles to processors.Each MPI processor retains the complete particle dataset under the copy algorithm.
  • Ring algorithm: The ring algorithm splits the particle dataset across MPI processors and reduces memory cost per process.Its trade-off is extra communication for neighbor particles absent from the local process.
  • Algorithm trade-off: The copy algorithm avoids extra neighbor-particle communication during irregular-force calculations but imposes a particle-number limit through memory use.The passage presents this as the principal trade-off between the two MPI strategies.

3.2 Basic NBODY6-GPU implementation

NBODY6-GPU combines OpenMP, CUDA, and AVX/SSE acceleration for workstation-scale simulations, assigning regular-force work to GPUs and irregular work to CPU vector parallelism. Its single-node design limits the particle number feasible within a reasonable simulation time.

  • Architecture: NBODY6-GPU uses OpenMP, CUDA, and AVX/SSE parallel methods to accelerate workstation and desktop computations.The implementation is shared-memory and does not provide a massively parallel MPI version.
  • Architecture: Single-node operation limits NBODY6-GPU to particle numbers feasible within a reasonable simulation time.The absence of massively parallel MPI implementation is the stated constraint.
  • GPU workload: O(N) regular-force calculation and O(N^2) potential-energy calculation are assigned to GPU acceleration.The GPU also efficiently accumulates neighbor lists during regular-force calculation.
  • CPU workload: AVX/SSE and OpenMP accelerate neighbor prediction and irregular-force calculation because the AC neighbor scheme is difficult to parallelize on GPUs.AVX/SSE supports vector calculation in the CPU cache.

3.3 Code improvements in NBODY6++GPU

NBODY6++GPU combines GPU, AVX/SSE, MPI, and OpenMP acceleration while redesigning initialization, prediction consistency, scheduling, and force correction for million-particle simulations.

  • Acceleration: GPU acceleration removes long-range regular-force calculation as the dominant computational cost.AVX/SSE then accelerates prediction and irregular neighbor-force calculations.
  • Time-step scheduling: Active-particle selection can cost O(N) per block time step, exceeding irregular-force work when N ≫ Ni⟨Nb⟩.A sorting-list algorithm is introduced to address the expense of searching all particles, including at deep block levels.
  • Initialization: MPI, GPU, and OpenMP parallelization improves initialization, particularly for large particle numbers and many primordial binaries.The revised initialization targets force, neighbor-list, and primordial KS-binary setup.
  • Initialization: The initialization handles masses, phase-space variables, N-body-unit scaling, forces, neighbor lists, time steps, and primordial KS binaries.Its four stages include parameter generation, scaling, force-related setup, and KS-binary initialization.
  • Prediction consistency: Parallel particle copies can produce inconsistent predictions and divergent time steps across MPI processors.The implementation predicts every particle to the current time before nonparallel stellar-evolution, KS, and hierarchical-regularization operations.
  • Stellar evolution: Neighbor-force correction remains costly during frequent stellar mass loss because corrections with cost O(N) per particle cannot be avoided.OpenMP accelerates but does not eliminate this performance loss.

3.4 Hybrid MPI parallelization

NBODY6++GPU uses a hybrid multi-node design in which MPI distributes nodes, OpenMP manages irregular-force threads and GPUs, and GPUs hold comparable particle-data subsets for regular calculations.

  • Hybrid parallelization: Each computing node runs one MPI process, which opens multiple OpenMP threads for irregular-force calculations.OpenMP threads control the GPUs within each node.
  • Hybrid parallelization: Each GPU receives a similar particle-dataset size for regular-force and potential-energy calculations.

4 PERFORMANCE TEST

NBODY6++GPU combines GPU, AVX/SSE, OpenMP, and MPI parallelization to accelerate regular and irregular integration across nodes. Tests show strong performance for million-particle simulations, while communication and KS calculations limit scaling in some binary-rich cases.

  • 4.1 Pure MPI and hybrid MPI: 33 times faster regular force integration is achieved with GPUs than with pure MPI.The comparison uses the same four-node CPU configuration, with GPUs added for the hybrid test.
  • 4.1 Pure MPI and hybrid MPI: 3 times faster irregular integration is achieved using AVX/SSE with OpenMP, while OpenMP reduces MPI communication cost by 5−10 times.The hybrid arrangement performs irregular-force work through threaded CPU vectorization and lowers communication frequency.
  • 4.2 Scaling with different particle numbers and processors: 400−2000 times speed-up over sequential NBODY6 is reached for one million particles on 16 nodes.This configuration uses 320 CPU cores and 86016 GPU cores; total time is about 800 s for one million particles.
  • 4.2 Scaling with different particle numbers and processors: Regular and irregular integration times are close to the ideal parallel limit for million-body simulations, but smaller systems depart from ideal scaling.For particle numbers below 10^5, internal memory access and modification dominate because each node performs fewer operations.
  • 4.3 Time fraction for different parts: Half of the calculation time is spent in KS procedures for 1024k particles with 5% binaries when Nnode ≥8, making KS the bottleneck.Without binaries, MPI communication and data movement consume about half the total time at 1024k particles on 16 nodes and 128k particles on 8 nodes.
  • 4.4 Sorting list algorithm for selecting active particles: 5 times faster sorting is obtained with the sorting list algorithm than with the temporary list algorithm for a 1024k-star model with 5% primordial binaries.The test uses a King sphere with W0 = 6 on 8 nodes.

5 APPLICATION

The million-particle application uses an NGC 4372-like globular-cluster model to examine computing performance and parameter adjustments in NBODY6++GPU. Adjusting the KS regularization criteria reduced computing time after initially unsuitable settings, while further models remained in progress.

  • 5 APPLICATION: The application section identifies million-particle globular-cluster modeling as the main target of NBODY6++GPU.The supplied figure captions also distinguish performance tests by node count, particle number, and primordial-binary fraction.
  • 5 APPLICATION: The simulation models a globular cluster with 1M stars, 5% primordial binaries, a 7.5 pc initial half-mass radius, and an 89.2 pc tidal radius.One N-body time unit corresponds to 0.622 Myr.
  • 5 APPLICATION: 3000 s per N-body time unit initially rose when small-time-step particles formed, then fell to about 1500 s after adjustments.The particle count decreased only slightly during 4500 N-body time units, while computing speed later increased.
  • 5 APPLICATION: The initial KS auto-adjustment produced Rcl = 1.4×10^-6 and ∆tcl = 6.8×10^-8 N-body units, which were too small to regularize many wide binaries.The simulation instead used Rcl = 5.0×10^-6 and ∆tcl ≤2.0×10^-7.
  • 5 APPLICATION: Enlarging the criteria to Rcl = 1.0×10^-5 and ∆tcl = 5.0×10^-7 sped up computing after 2800 time units.The original criterion was designed for N = 10^2–10^3 and was unsuitable for the million-body model’s outer region.
  • 5 APPLICATION: Several additional models were still in progress and reserved for a future publication.The paper therefore reports only an initial application set here.

6 DISCUSSION

The discussion identifies scaling and deployment limits arising from communication, KS processing, memory, coding effort, and cluster scheduling. It also outlines future improvements and documentation needs for broader use.

  • 6 DISCUSSION: 86,000 GPU cores and 320 CPU cores mark the approximate saturation point of the reported performance scaling.The authors attribute this comparison partly to NBODY6 and NBODY6++ requiring fewer operations than standard Hermite codes.
  • 6 DISCUSSION: Data movement and MPI communication become bottlenecks at large node counts, while KS integration dominates when many primordial binaries are present.The calculation-to-communication ratio differs across regular force, irregular force, and KS perturbation calculations.
  • 6 DISCUSSION: Table 2 estimates calculation and communication costs for the principal force and KS perturbation components.The table supports the discussion’s analysis of which components scale well and which become communication-limited.
  • 6 DISCUSSION: KS parallelization benefits from separating unperturbed and perturbed binaries because unperturbed systems require few operations and typically have no perturbers.The discussion also notes a separate effort to parallelize KS with block time steps.
  • 6 DISCUSSION: KS initialization and termination can be costly for wide binaries that frequently switch between KS and Hermite solutions.Proposed improvements include reusing existing values and reverse neighbor lists, but the latter requires substantial memory and coding effort.
  • 6 DISCUSSION: NBODY6++GPU is unsuitable for clusters that reserve most CPU cores for non-GPU tasks because it relies heavily on CPU computation and data movement.Shared-node contention can make MPI barrier time approach half of total computing time.
  • 6 DISCUSSION: The paper provides documentation because the long-developed NBODY6 and NBODY6++ codes have become difficult for beginners.The documentation includes detailed input and output descriptions, with further implementation details planned.
  • 6 DISCUSSION: Future hardware and code improvements may enable larger particle numbers, but communication and data management remain key constraints near central black holes.The authors suggest hybrid codes for Exaflop/s-scale regimes rather than full N^2 force computation everywhere.

7 CONCLUSIONS

NBODY6++GPU addresses the longstanding challenge of direct simulations of large globular clusters by combining multiple parallelization and acceleration methods. The code reaches million-particle simulations with realistic physics and substantially reduces computational time.

  • Direct simulations of realistic globular clusters had remained difficult because hardware and software limits constrained large-particle-number calculations.Earlier work argued that breakthroughs in parallel computing were needed to overcome this challenge.
  • NBODY6++GPU combines MPI, GPU, and AVX/SSE acceleration with NBODY6++ for massively parallel multi-node simulations.The code integrates GPU and vector-instruction acceleration on each node.
  • 400-2000 times faster than NBODY6 is achieved for non-binary cases on 320 CPU cores and 32 NVIDIA K20X GPUs.GPU acceleration makes long-range force calculations about 33 times faster, while AVX/SSE accelerates prediction and neighbor forces by a factor of 3.
  • One million particles can be simulated in a realistic globular cluster with stellar evolution and 5% primordial binaries for several Gyr.The reported computational cost is about one hour per half-mass crossing time.
  • About 20 crossing times per day are possible for a million-body cluster on 320 cores with 32 GPUs.The authors identify communication, bandwidth, latency, and data-management improvements as routes toward larger simulations.

APPENDIX A: DIFFERENCES BETWEEN NBODY6++ AND NBODY6

The appendix notes that NBODY6++GPU differs from NBODY6 and directs readers to the code manual for the detailed comparison.

  • NBODY6++GPU has several important differences from NBODY6, with further details provided in Table A and the code manual.
Loading 1504.03687v2…