A sparse octree gravitational N-body code that runs entirely on the GPU processor

Jeroen Bédorf; Evghenii Gaburov; Simon Portegies Zwart

Source-linked AI summary

A sparse octree gravitational N-body code that runs entirely on the GPU processor

Jeroen Bédorf, Evghenii Gaburov, Simon Portegies Zwart

arXiv:1106.1900v2Jun 9, 2011astro-ph.IMcs.DC

TL;DR

Large particle simulations need efficient hierarchical searches, but block time-step methods reduce available parallelism and shift costs toward construction, traversal, and integration. This paper moves octree construction, traversal, and time integration onto the GPU, achieving 2.8 million particles per second on a GTX480 while supporting portable many-core implementations.

Problem
For large particle sets, hierarchical octrees reduce search complexity, while block time-step simulations reduce parallelism and make non-interaction computations more significant.
Method
The authors develop GPU algorithms for octree construction and traversal, execute particle integration on the GPU, and implement them in a CUDA gravitational N-body code.
Results
2.8 million particles per second are processed at θ = 0.75 on a GTX480, with sorting, moving, tree construction, and multipole computation scaling linearly for N ≳ 10^6.
Takeaways & Limitations
The algorithms are portable to other massively parallel architectures and extend beyond gravitational forces to neighbor search, clump finding, fast multipole methods, and ray tracing.
Takeaways & Limitations
The calculations use single precision arithmetic, except for specified monopole, quadrupole, and force computations that use double precision.

Abstract

from arXiv · show

We present parallel algorithms for constructing and traversing sparse octrees on graphics processing units (GPUs). The algorithms are based on parallel-scan and sort methods. To test the performance and feasibility, we implemented them in CUDA in the form of a gravitational tree-code which completely runs on the GPU.(The code is publicly available at: http://castle.strw.leidenuniv.nl/software.html) The tree construction and traverse algorithms are portable to many-core devices which have support for CUDA or OpenCL programming languages. The gravitational tree-code outperforms tuned CPU code during the tree-construction and shows a performance improvement of more than a factor 20 overall, resulting in a processing rate of more than 2.8 million particles per second.

1. Introduction

Sparse octrees support efficient hierarchical computations for large, non-homogeneous particle distributions, but constructing and traversing them is difficult to parallelize on GPUs. The paper addresses this by moving tree construction, traversal, and integration onto the GPU.

Motivation: Sparse octrees recursively partition three-dimensional space and are suited to non-homogeneous particle distributions.The paper focuses on sparse rather than dense octrees, whose empty branches reflect the underlying particle distribution.
Motivation: For large particle sets, octrees reduce range-search complexity from O(N) to O(log N) per particle.Hierarchical methods become more efficient than brute-force methods for approximately N ≥ 10^5 particles.
Motivation: High-resolution simulations require GPUs or costly computing clusters because a single CPU is insufficient for approximately N ≥ 10^5 particles.GPUs offer an attractive alternative for these computationally demanding workloads.
Motivation: Block time-step simulations expose a GPU-parallelism bottleneck because only a fraction of particles may be updated at each step.Construction, traversal, and integration become relatively more important when fewer particles are active.
Contribution: The paper redesigns tree construction, traversal, and time integration for GPU execution, eliminating host-device transfers and supporting shared and block time steps.Although demonstrated in a gravitational N-body code, the algorithms are intended for related hierarchical methods.

2. Sparse octrees on GPUs

The GPU algorithms construct sparse octrees by sorting Morton keys and grouping particles level-by-level, then traverse them with breadth-first, stack-based parallel processing. Parallel scans, compacting, and GPU threads provide the principal sources of parallelism.

2.1. Tree construction: Morton keys preserve spatial locality while mapping three-dimensional particle positions into a sortable one-dimensional representation.Sorting particles by increasing Morton key improves tree construction and cache efficiency during traversal.
2.1. Tree construction: At each level, masked Morton keys group particles into cells, with cells containing at most Nleaf = 16 particles becoming leaves.Particles assigned to leaves are flagged and excluded from subsequent levels; larger cells remain nodes.
2.1. Tree construction: Tree cells are linked in parallel by locating each cell’s parent and first child, then updating the parent’s child counter.The linking step uses the sorted Morton-key array and connects parent-child relationships after construction.
2.2. Tree traverse: Breadth-first traversal processes groups of particles with GPU thread blocks rather than recursively traversing one particle at a time.The method uses groups of at most Ncrit particles and associates each group with a thread block.
2.2. Tree traverse: A stack-based traversal sends node children to the next-level stack, leaves to particle-cell interaction lists, and terminated cells to cell-group lists.Children and leaf particles are stored contiguously, so traversal needs only first-child or first-particle indices and counts.

3. Gravitational Tree-code

The gravitational application implements an entirely GPU-resident Barnes–Hut code, including integration, multipole computation, and tree traversal. Its accuracy-speed trade-off is controlled by the cell-opening criterion and θ.

3. Gravitational Tree-code: The implementation runs the gravitational Barnes–Hut tree code entirely on GPUs, including tree construction, traversal, integration, and cell-property computation.The implementation is based on the sparse-octree methods and is distinct from GPU codes that retain CPU-side execution.
3.1. Time integration: For clusters larger than approximately 10^5 particles, prediction-correction integration takes less than 1% of total execution time.This timing refers to the leapfrog predictor-corrector step.
3.2. Tree-cell properties: Multipole moments are computed from leaves upward, with each node calculated from the moments of cells one level below.The implementation uses monopole and quadrupole moments and computes node levels in parallel.
3.3. Cell opening criterion: The multipole acceptance criterion balances accuracy against execution time by deciding whether to accept a cell or open its children.Larger θ values open fewer cells and run faster, whereas smaller θ values traverse farther and improve accuracy.
3.3. Cell opening criterion: Using the minimal-distance MAC reduces acceleration error by 10%–50% for the same θ but makes computation almost 3 times slower.The higher cost results from accepting, or opening, more cells.

4. Performance and Accuracy

Bonsai evaluates GPU-based tree construction and traversal against tuned CPU and GPU implementations, measuring scaling, acceleration accuracy, and galaxy-merger behavior. The GPU algorithms show linear tree-construction scaling, substantial speedups over CPU code, competitive acceleration errors, and stable merger-simulation results.

4.1. Performance: Linear scaling appears for sorting, moving, tree construction, and multipole computation once N ≳ 10^6, whereas smaller runs underutilize parallel scans.More than 90% of wall-clock time is spent in tree traversal for θ = 0.75 and θ = 0.5.
4.1. Performance: Almost 30× separates the fastest GPU from the CPU for tree traversal, while CPU–GPU gaps across algorithms range from 2× to almost 30×.The comparison uses tuned CPU implementations and five GPUs spanning three generations.
4.2.1. Acceleration: Acceleration errors are slightly below Octgrav’s and comparable to CPU tree-code results, while remaining nearly independent of particle count.The comparison varies particle number and opening angle θ.
4.2.2. Galaxy merger: Merger-product density, cumulative-mass, and circular-velocity profiles are identical across codes apart from slight small-number-statistics differences.Black-hole separations are indistinguishable through the third pericenter passage at t = 300, after which numerical differences make them incomparable.
4.2.2. Galaxy merger: Bonsai exhibits more stable energy evolution than Octgrav and Partree, while reducing maximum energy error by roughly 2× when θ decreases from 0.75 to 0.5.The reported maximum occurs during the second pericenter passage at t ≈ 280.
4.2.2. Galaxy merger: Bonsai runs 1.6× faster than Octgrav at θ = 0.75 and 1.26× faster at θ = 0.5.The smaller speedup at θ = 0.5 is attributed to faster tree traversal in Octgrav.

5. Discussion and Conclusions

The paper presents GPU algorithms for constructing and traversing hierarchical data structures, integrating particle time integration to eliminate host–GPU transfers. Bonsai achieves high throughput, supports block-time-step simulations, and offers portability beyond gravitational N-body calculations.

5. Discussion and Conclusions: 2.8 million particles per second are processed with θ = 0.75 on one GTX480, while stable energy evolution supports routine GPU simulations.The current implementation uses one GPU, with memory limiting 5-million-particle simulations to approximately 1 gigabyte.
5. Discussion and Conclusions: Bonsai is more appropriate than Octgrav for large N > 10^6 and block-time-step simulations despite Octgrav’s approximately 10% faster tree traversal.Octgrav transfers the complete tree structure, particle array, and multipole moments to the GPU at every time step.
5. Discussion and Conclusions: GPU construction, traversal, and integration remove host communication and are particularly suitable for shared- and block-time-step simulations.For block time steps, shifting construction, prediction, and communication to the GPU addresses bottlenecks caused by fewer active particles.
5. Discussion and Conclusions: The algorithms are portable to other massively parallel architectures and applicable to neighbour search, clump finding, fast multipole methods, ray tracing, and smoothed particle hydrodynamics.OpenCL portability requires rewriting only the GPU functions, and the OpenCL version is described as work in progress.

Appendix A. Scan algorithms

The appendix introduces the all-prefix-sums operation as a foundation for the scan algorithms used in the paper. It defines scan over an array using a binary associative operator.

Appendix A. Scan algorithms: Parallel scan, or prefix sum, applies a binary associative operator ⊕ across an array to compute all prefix results.The definition begins with an array of n elements and an associative operator.

and returns the ordered set

The appendix illustrates prefix sums with addition and identifies their use as building blocks for several parallel-processing methods. These methods include compaction, splitting, sorting, and tree-related operations.

and returns the ordered set: With addition as ⊕, an all-prefix-sums operation produces the ordered cumulative values shown in the example array.The example’s output is [3 4 11 11 15 16 22 25].
and returns the ordered set: Prefix-sum algorithms support stream compaction, stream splitting, sorting, regular expressions, tree-depth determination, and histogram construction.The paper uses them as building blocks for GPU methods and refers readers to Blelloch for further descriptions.

Appendix A.1. Stream Compaction

Stream compaction removes invalid stream elements and places valid elements at the beginning of the output. The appendix presents it as an example of stream reduction.

Appendix A.1. Stream Compaction: Stream compaction removes invalid elements and packs valid elements at the start of the output stream.The paper also calls this operation stream reduction.

Appendix A.2. Split and Sort

Stream split preserves all input items while moving invalid elements behind valid ones, making it suitable for radix sort.

Stream split places invalid elements behind valid ones instead of discarding them, preserving item ordering for radix sort.The algorithm is applied once per integer bit from least to most significant.

Appendix A.3. Implementation

The implementation uses parallel scan-based primitives for GPU stream processing, including split, compaction, radix sort, and prefix-sum operations.

The implementation uses scan algorithms throughout the parallel octree construction and traversal pipelines.The authors select an efficient many-core scan method to avoid making scans an application bottleneck.
The stream-processing method counts valid elements, computes output offsets with parallel prefix-sum, and places valid elements at those offsets.
The chosen prefix-sum method is applied to operations in both tree construction and tree traversal.

Appendix B. Morton Key generation

Morton keys are generated directly on the GPU by converting positions to integers and dilating their bits into a locality-preserving key.

GPU Morton-key generation converts floating-point particle positions into integer positions before bit-based dilation.Direct computation avoids relatively inefficient lookup-table access by many parallel GPU threads.
The dilation primitive converts the first 10 bits of an integer into a 30-bit representation used to form Morton keys.
The implementation stores 60-bit particle keys in two 32-bit integers, supporting octrees with a maximum depth of 20 levels.The depth limit could be extended to 30 levels with 90-bit keys or by modifying construction at the deepest levels.