Source-linked AI summary
High Performance Direct Gravitational N-body Simulations on Graphics Processing Units -- II: An implementation in CUDA
Robert G. Belleman, Jeroen Bedorf, Simon Portegies Zwart
TL;DR
The paper examines whether a commercial GPU can accelerate gravitational N-body force calculations as an alternative to specialized GRAPE hardware. It implements CUDA for direct and Barnes-Hut simulations, finding competitive performance for sufficiently large systems while retaining precision and data-transfer limitations. GPUs therefore offer an attractive alternative to GRAPE-6Af within the tested scope.
Problem
The paper investigates whether programmable commercial GPUs can provide an alternative to specialized GRAPE hardware for gravitational N-body simulations.
Method
The authors implement CUDA force evaluation on an NVIDIA GeForce 8800GTX and test two direct Hermite-integrator codes plus a Barnes-Hut treecode, with motion integration on the host CPU.
Results
For N > 512 with softening, kirin performs comparably to GRAPE-6Af; for treecodes it outperforms GRAPE across all N, while energy errors remain of order |∆E|/E ∼10^-7.
Takeaways & Limitations
Modern GPUs offer an attractive alternative to GRAPE-6Af, particularly for large particle counts and methods compatible with their parallel performance and precision.
Takeaways & Limitations
GPU performance is reduced by CPU-GPU memory-copy overhead and limited single precision; without softening or with very small time steps, GRAPE-6Af can be faster.
Abstract
from arXiv · showhide
We present the results of gravitational direct $N$-body simulations using the Graphics Processing Unit (GPU) on a commercial NVIDIA GeForce 8800GTX designed for gaming computers. The force evaluation of the $N$-body problem is implemented in ``Compute Unified Device Architecture'' (CUDA) using the GPU to speed-up the calculations. We tested the implementation on three different $N$-body codes: two direct $N$-body integration codes, using the 4th order predictor-corrector Hermite integrator with block time-steps, and one Barnes-Hut treecode, which uses a 2nd order leapfrog integration scheme. The integration of the equations of motions for all codes is performed on the host CPU. We find that for $N > 512$ particles the GPU outperforms the GRAPE-6Af, if some softening in the force calculation is accepted. Without softening and for very small integration time steps the GRAPE still outperforms the GPU. We conclude that modern GPUs offer an attractive alternative to GRAPE-6Af special purpose hardware. Using the same time-step criterion, the total energy of the $N$-body system was conserved better than to one in $10^6$ on the GPU, only about an order of magnitude worse than obtained with GRAPE-6Af. For $N \apgt 10^5$ the 8800GTX outperforms the host CPU by a factor of about 100 and runs at about the same speed as the GRAPE-6Af.
1 Introduction
GPUs evolved into highly parallel programmable processors with substantial raw performance, but programming complexity and limited precision constrained their scientific use. This paper presents a CUDA implementation for gravitational N-body simulations using direct integration and a Barnes-Hut treecode.
- GPU motivation: 350 GFLOP/s on the GeForce 8800GTX illustrates the raw performance advantage of GPUs over conventional CPUs.The hardware is optimized for vectorizable parallel operations and streaming large data volumes.
- GPU motivation: GPU programming progressed from hardware-specific assembly to higher-level shading languages that let programmers focus on the computational problem.Cg and GLSL reduced the need for detailed hardware-specific programming knowledge.
- Scientific applications: Scientific GPU computing was demonstrated for PDE solvers, ray tracing, image segmentation, and gravitational simulations.Early implementations often represented scientific data as colored pixels stored in textures.
- GPU limitations: Single-precision GPU arithmetic limited applications requiring higher numerical precision.CUDA release notes anticipated GPUs supporting 64-bit double-precision arithmetic in late 2007.
- Paper contribution: The paper applies a CUDA implementation to gravitational N-body systems using both direct integration and a Barnes-Hut treecode.It extends the authors’ earlier Cg implementation, which outperformed CPUs by about an order of magnitude for N >∼10^4.
2 Background
Direct gravitational N-body simulation requires an expensive, parallelizable force calculation, motivating specialized GRAPE hardware and programmable GPU alternatives. CUDA exposes GPU parallelism and memory hierarchies for implementing these computations, while limited precision favors some methods over others.
- N-body problem: The total force on each particle is obtained by summing pairwise forces from all N particles.The force equation uses the Newton constant, particle masses, and particle positions.
- Specialized hardware: 64 TFLOP/s peak speed made GRAPE6 a powerful special-purpose N-body computer, while the GRAPE-6Af provided 123 GFLOP/s on a PCI card.The GRAPE-6Af used four GRAPE6 chips and supported up to 131072 particles.
- GPU alternatives: GPUs provide an alternative coprocessor because many processing units perform the same operations on different input streams using SIMD.Earlier GPU implementations used shared time steps or softening, while paper I used block time steps for force, jerk, and potential calculations.
- CUDA model: CUDA exposes the GPU as a parallel streaming processor and provides more flexible hardware mapping than earlier interfaces such as Cg.CUDA programs contain GPU kernels executed by multiple threads organized into bundles that communicate through shared memory.
- CUDA model: GPU memory is hierarchical: larger memories are more flexible but slower, while smaller memories are more restrictive but faster.Efficient CUDA mapping therefore requires storing frequently used data appropriately within the hierarchy.
- Numerical precision: 32-bit GPU precision hinders high-precision direct integrations but is suitable for lower-precision methods such as Barnes-Hut treecodes.The implementation was tested with both direct integration and a treecode.
3 Implementation
The implementation divides N-body integration between the CPU and GPU: the GPU computes force-related quantities, while the CPU performs prediction and correction. CUDA kernels use block time-steps, shared-memory tiling, and a GRAPE-compatible library interface to improve utilization and portability.
- 3.1 Decomposition over CPU and GPU: The Hermite scheme uses block time-steps, with prediction and correction on the CPU and direct-summation acceleration, jerk, and potential on the GPU.Only particles requiring updates are integrated in each block.
- 3.1 Decomposition over CPU and GPU: CPU–GPU memory copies occur at every block time-step, creating greater overhead than GRAPE-6Af because prediction remains on the CPU.The exchanged data includes particle mass, position, and velocity as inputs, and acceleration, jerk, and potential as outputs.
- 3.1 Decomposition over CPU and GPU: Each GPU thread processes one particle while bundles repeatedly cache particle data into shared memory and accumulate results in registers.The kernel processes all N particles in shared-memory tiles, removes self-interaction from the potential, and writes outputs to global memory.
- 3.1 Decomposition over CPU and GPU: The GPU performs N^2 calculations rather than 1/2N(N −1), avoiding costly internal communication and synchronization.The extra arithmetic is accepted because reducing the work would lower GPU performance through communication overhead.
- 3.1 Decomposition over CPU and GPU: Shared memory reduces global-memory latency by reusing cached particle data, while 128-thread bundles calculate 128 particles in parallel.Global memory is slow and shared memory is fast but limited, so particles are precached and processed in consecutive bursts.
- 3.2 Optimizing GPU utilization: The implementation uses single-precision values and splits small blocks below 4096 particles into sequential parts to improve GPU utilization.A GRAPE6-compatible library lets existing GRAPE-linked applications use kirin with minimal changes; branch checks can reduce performance by roughly 10% each.
4 Results
The CUDA implementation is evaluated against GRAPE-6Af, the earlier Cg implementation, and a host-only CPU across direct and treecode N-body workloads. Results report performance, peak GPU throughput, simulation energy error, and the effects of softening and algorithm-specific optimization.
- Direct N-body integration: The experiments compare CUDA kirin with GRAPE-6Af, Cg on the GeForce 8800GTX, and a host-only Intel Xeon implementation.Direct simulations run over 0.5 N-body time units, with timings measured from t = 0.25 to t = 0.5; the baseline softening is 1/256.
- Direct N-body integration: 340 GFLOP/s peak force-calculation performance nearly reaches the 8800GTX theoretical peak of 346 GFLOP/s.These measurements use force calculation only; performance is computed with k = 38 in equation 3.
- Accuracy measurements: Energy-error measurements compare implementations over 0.5 N-body time units with ǫ = 1/256 using the same input parameters as the performance tests.The relative error is defined by comparing the system energy at the start and end of the simulation.
- N-body integration using the treecode: The treecode uses a GPU adaptation of a GRAPE-designed Barnes-Hut implementation, with a second version omitting jerk and nearest-neighbour calculations.Because the treecode integrates using acceleration only, omitting these calculations produces a performance gain of a factor of two.
5 Discussion
The CUDA-based kirin library is competitive with GRAPE-6Af for sufficiently large systems, especially with softening, while performance and accuracy depend on integration settings and GPU utilization. Treecode execution benefits from avoiding GRAPE memory-transfer limits, although GPU precision and hardware-specific tuning remain constraints.
- GPU adoption offers a lower-cost, programmable alternative to GRAPE, but single precision remains problematic and the comparison uses only the smallest one-module GRAPE-6Af.
- N <~512 particles: GRAPE-6Af remains faster by about a factor of two with block time-steps.
- N > 512 particles: CUDA kirin reaches comparable speed to GRAPE-6Af because most block time-steps use the GPU at full capacity.
- GPU performance depends on GPU-specific thread and bundle counts constrained by available registers, making generally optimal configuration difficult.
- For N > 512 and ǫ = 1/256, kirin performs comparably to GRAPE-6Af; without softening, smaller integration steps reduce kirin performance below GRAPE-6Af.
- Increasing integrator accuracy by twofold doubles GRAPE computation time and reduces energy error 24-fold, but increases GPU computation time about tenfold while barely reducing its energy error.
- The treecode GPU is an order of magnitude faster than the CPU and outperforms GRAPE for all N because it avoids memory-transfer limitations; GPU and GRAPE energy errors are comparable.
- Across N = 256 to 65536, GPU and GRAPE relative energy errors are of order |∆E|/E ~10^-7, while reduced time steps further improve GRAPE accuracy but not GPU accuracy.
A kirin library functions
The kirin library provides GPU counterparts to GRAPE6 functions, enabling existing GRAPE6-based code to use the GPU with recompilation and relinking. Its functions cover connection management, configuration, particle transfer, prediction, and force-calculation control.
- All GRAPE6 library functions have equivalent GPU implementations, so existing code only needs recompilation and relinking to use them.
- The open function connects to the GPU and initializes local buffers, while close releases local and GPU memory.
- The npipes function configures the GPU pipeline count, which can affect performance and is not fixed as on GRAPE.
- The particle function buffers particle data locally before sending it to the GPU after prediction.
- The ti function sets the next time step, starts host-side prediction, and sends predicted particles to the GPU.
- The firsthalf and lasthalf functions respectively start force calculations and retrieve their results for specified particles.