Source-linked AI summary
Accelerating NBODY6 with Graphics Processing Units
Keigo Nitadori, Sverre J. Aarseth
TL;DR
NBODY6’s direct-force cost limits larger N-body simulations. The paper accelerates its regular-force calculation with GPUs while handling irregular forces on the CPU, yielding a balanced NBODY6-GPU implementation for difficult simulations and particle numbers up to 2 × 10^5 on dual-GPU systems.
Problem
The N^2 cost of direct force calculations has hindered extending the particle number in NBODY6 simulations.
Method
The paper places regular-force evaluation on GPUs and uses SSE/AVX and OpenMP CPU procedures for irregular forces, with additional prediction and regularization strategies.
Results
The current NBODY6-GPU formulation performs well under difficult conditions and is balanced for simulations with 10^4-2 × 10^5 particles on dual-GPU systems.
Takeaways & Limitations
The regular component of the Ahmad–Cohen neighbour scheme is well suited to GPU hardware, while irregular-force acceleration is better handled by optimized CPU code.
Takeaways & Limitations
Strict reproducibility requires omitting some non-thread-safe parallel procedures, at a loss of efficiency, and GPU execution is not beneficial for the irregular force in the current implementation.
Abstract
from arXiv · showhide
We describe the use of Graphics Processing Units (GPUs) for speeding up the code NBODY6 which is widely used for direct $N$-body simulations. Over the years, the $N^2$ nature of the direct force calculation has proved a barrier for extending the particle number. Following an early introduction of force polynomials and individual time-steps, the calculation cost was first reduced by the introduction of a neighbour scheme. After a decade of GRAPE computers which speeded up the force calculation further, we are now in the era of GPUs where relatively small hardware systems are highly cost-effective. A significant gain in efficiency is achieved by employing the GPU to obtain the so-called regular force which typically involves some 99 percent of the particles, while the remaining local forces are evaluated on the host. However, the latter operation is performed up to 20 times more frequently and may still account for a significant cost. This effort is reduced by parallel SSE/AVX procedures where each interaction term is calculated using mainly single precision. We also discuss further strategies connected with coordinate and velocity prediction required by the integration scheme. This leaves hard binaries and multiple close encounters which are treated by several regularization methods. The present nbody6-GPU code is well balanced for simulations in the particle range $10^4-2 \times 10^5$ for a dual GPU system attached to a standard PC.
1 INTRODUCTION
Direct N-body calculations have long been limited by their O(N^2) force cost. NBODY6 advances this problem through neighbour schemes, specialized hardware, and GPU-oriented treatment of regular and irregular forces.
- The Ahmad–Cohen neighbour scheme split the total force into slowly changing regular and shorter-timescale irregular components.
- GRAPE-type computers later accelerated force calculations, followed by increasingly available GPUs and CUDA programming.
- The paper focuses mainly on stand-alone NBODY6–GPU systems using one or two GPUs.
- NBODY6++ is noted as a parallel alternative capable of larger particle numbers with several hardware types and intended GPU support.
- The implementation accelerates regular forces on GPUs while evaluating irregular forces with parallel SSE and OpenMP because GPU overheads are too large for that component.
2 BASIC NBODY6 CODE
NBODY6 combines Hermite integration, block-steps, neighbour selection, and regularization to integrate systems with individual time-steps and close encounters. The implementation also omits optional stellar-evolution functionality for simplicity in GPU development.
- NBODY6 replaced the fourth-order Adams method with an equivalent Hermite formulation using two force polynomials.
- Hermite integration with hierarchical block-steps reduces coordinate and velocity prediction overheads and supports parallel procedures.
- Neighbour radii are adjusted using local density contrast, while neighbour changes avoid derivative-correction overhead when results are requested at commensurate times.
- Two-body and multiple close encounters are handled with KS and chain regularization, including high-order Hermite integration for KS equations.
- Synthetic stellar evolution is omitted from the GPU-development discussion because it is optional and adds relatively little host CPU time.
3 NEW IMPLEMENTATIONS
The new implementation assigns regular-force calculation and neighbour-list construction to CUDA GPUs while handling irregular forces on CPUs with SSE/AVX and OpenMP. It also adapts neighbour selection, prediction, precision, and force consistency to the block-step integration scheme.
- Software design: NBODY 6-GPU uses CUDA with multiple-GPU support for regular forces and CPU SSE/AVX plus OpenMP for irregular forces.The GPU library is GPUNB, while GPUIRR uses the CPU because the irregular calculation is fine-grained.
- Velocity neighbour criterion: The velocity-based neighbour criterion permits larger regular time-steps while adding high-velocity particles before close approach.It increases the cost of each regular pairwise interaction but provides a safety condition for rapidly moving particles.
- Block-step procedure: Each block-step obtains the next integration time, forms active regular and irregular lists, predicts particles, evaluates forces, corrects particles, and updates neighbours.The regular corrector also adjusts force polynomials when the neighbour list changes.
- The GPUNB library: GPUNB evaluates NiNj interactions in one call, skips neighbour contributions during accumulation, and records neighbour indices for all active particles in parallel.GPU thread-level parallelism covers multiple i-particles, while reductions and list gathering limit host data transfer.
- Precision and force evaluation: The implementation uses mainly single precision for regular forces, but two-float coordinates for irregular close encounters and retains sufficient precision for force summation.Regular-force close interactions are excluded, whereas irregular calculations may require higher coordinate precision.
- Prediction: Prediction strategy switches between all-particle and necessary-particle prediction because neighbour-list sorting and random memory access can make selective prediction expensive.A turn-around point determines which strategy is used for a block-step.
4 PERFORMANCE
Performance tests show that GPU acceleration substantially reduces regular-force costs, while irregular-force, host-side, reproducibility, and advanced-evolution considerations remain important. The implementation achieves strong scaling and long-term energy stability across tested particle numbers, but has stated scope and performance caveats.
- System performance: 56 times faster at N = 24 000 than the corresponding original NBODY 6 time on System A, with Tcomp ∝N 1.8 asymptotically.The comparison concerns wall-clock timing; regular-force performance on System A exceeded 1 TFLOPS.
- GPU scaling: 10-fold regular-force gain from the tuned NBODY 6–SSE version to one GPU, doubled again when using two GPUs.Total simulation time scales less well because other components remain costly.
- GPU scaling: 45 to 23 sec for the regular force at N = 64 k with two GPUs, while total time decreases from 99 to 77 sec.The authors suggest running two separate simulations simultaneously on a dual-GPU machine when N is not too large.
- Reproducibility and limitations: Strict reproducibility requires omitting two non-thread-safe OpenMP treatments, reducing efficiency; timings instead used full optimization.The early stage was considered essentially reproducible despite the optimized timing configuration.
- Accuracy and stability: 1.5 × 10^-7 relative energy error for tested N = 128 and 256 k, compared with 4 × 10^-7 for the original code at N = 8 and 16 k.The GPU potential-energy evaluation has an intrinsic relative error of 1 × 10^-8.
- Accuracy and stability: −1 × 10^-4 accumulated total-energy change at minimum core radius in an N = 32 k equal-mass core-collapse study.A comparable drift was observed beyond core collapse for N = 16 k, supporting the reported long-term stability of the Hermite scheme.
5 CONCLUSIONS
The implementations divide regular and irregular force calculations across GPUs and host-side parallel procedures, achieving significant gains and comparable wall-clock times across particle numbers. Multiple GPUs scale the regular-force calculation, while AVX further accelerates the irregular component.
- Force calculation: The regular force was implemented on GPUs with GPUNB, while the irregular force uses the separate GPUIRR library because GPU overheads were too large for that component.The two force components therefore require different acceleration strategies.
- Force calculation: SSE and OpenMP reduced irregular-force wall-clock time until it became comparable with the regular-force calculation across a range of particle numbers.
- Hardware scaling: AVX support added further speed-up to the irregular-force calculation.
- Hardware scaling: Multiple GPUs benefit the regular-force calculation, which scales well and may gain further speed-up with four GPUs and an octocore CPU.
- Scope: The current NBODY 6-GPU formulation performs well under a variety of difficult conditions, although large N-body simulations remain expensive.
APPENDIX A: GLOSSARY
The glossary defines the software, hardware, and parallel-computing terminology used to describe GPU-accelerated force calculations. It also distinguishes particles and loop-level forms of parallelism.
- Software: GCC denotes the GNU Compiler Collection, including gcc, g++, gfortran, and other languages.
- Software: API denotes an Application Programming Interface providing definitions or prototypes for FORTRAN subroutines and C functions.
- Hardware and instructions: SSE and AVX are x86 vector instruction sets supporting four single-precision values on 128-bit registers and eight single-precision or four double-precision values on 256-bit registers, respectively.
- GPU computing: CUDA is a framework for general-purpose computing on NVIDIA GPUs, including language, compiler, runtime library, and device driver.
- Parallel computing: OpenMP uses language directives to utilize multiple processors in shared memory, while SIMD applies one operation to multiple data values.
- Particle and loop roles: An i-particle feels gravitational force, a j-particle supplies it, and i- and j-parallelism parallelize the outer and inner interaction loops, respectively.
APPENDIX B: CUDA CODES
The CUDA appendix presents the innermost regular-force kernel and its supporting data structures. Each GPU thread handles one i-particle, accumulates interactions from j-particles, and includes additional velocity-criterion operations.
- Code listings: Listings B1 and B2 provide CUDA C++ definitions of data structures and the innermost gravity-calculation kernel.
- Kernel structure: Each CUDA thread stores one i-particle in registers, evaluates its force from particle j, and accumulates the result.
- Kernel structure: The kernel uses a thread block with NNB_PER_BLOCK set to 256, a power-of-two number of neighbours per block.
- Force calculation: The gravity calculation forms squared distance and position-velocity products using single-precision variables r2 and rv.
- Velocity criterion: The velocity criterion adds 5 multiplications, 6 additions, and 1 minimum operation to the kernel.
APPENDIX C: FORTRAN CODES
The FORTRAN appendix traces one block-step cycle through prediction, regular-force evaluation, irregular-force correction, neighbour-list updates, and selection of the next block time. It shows how GPU libraries and OpenMP procedures are coordinated during integration.
- Block-step control: The block-step cycle begins by identifying the next block time and forming candidate lists for new irregular and regular forces.
- Irregular-force processing: The program predicts active or all particles, evaluates new irregular forces and derivatives, and chooses between sequential and parallel irregular integration.
- Irregular-force processing: For small active sets, irregular steps are corrected sequentially; otherwise, irregular correction is performed in parallel with OpenMP.
- Regular-force processing: The regular-force stage predicts particles on the host, sends particles to the GPU, and evaluates forces, derivatives, and neighbour lists in blocks of NIMAX=1024.
- Neighbour-list handling: After GPU evaluation, neighbour lists are copied back and registered with the irregular-force library for each particle.
- Correction and update: The cycle sends corrected active particles to GPUIRR, performs Hermite correction, updates changed neighbour lists, and determines the next block time.