Source-linked AI summary
REBOUND: An open-source multi-purpose N-body code for collisional dynamics
Hanno Rein, Shang-Fei Liu
TL;DR
REBOUND addresses the need for an open-source, multi-purpose N-body code for collisional dynamics and related systems. It combines modular integrators, gravity and collision algorithms, boundary conditions, and parallel implementations, then evaluates their accuracy and scaling. The authors report efficient operation from desktop machines to large clusters, with particle capacity practically limited by available CPUs in weak scaling tests.
Problem
REBOUND targets collisional N-body problems, including planetary rings, for which the authors report no publicly available code capable of solving the described cases.
Method
The paper develops a modular open-source N-body code with symplectic integrators, tree and plane-sweep collision or gravity modules, configurable boundaries, and MPI/OpenMP parallelization.
Results
The octree self-gravity and collision modules run efficiently on multi-core desktops and large clusters, with weak scaling showing no practical particle limit beyond available CPUs.
Takeaways & Limitations
REBOUND provides a customizable framework for collisional and collision-less N-body simulations across varied astrophysical settings and computing platforms.
Takeaways & Limitations
Collision detection is approximate and can miss collisions or report collisions that do not occur.
Abstract
from arXiv · showhide
REBOUND is a new multi-purpose N-body code which is freely available under an open-source license. It was designed for collisional dynamics such as planetary rings but can also solve the classical N-body problem. It is highly modular and can be customized easily to work on a wide variety of different problems in astrophysics and beyond. REBOUND comes with three symplectic integrators: leap-frog, the symplectic epicycle integrator (SEI) and a Wisdom-Holman mapping (WH). It supports open, periodic and shearing-sheet boundary conditions. REBOUND can use a Barnes-Hut tree to calculate both self-gravity and collisions. These modules are fully parallelized with MPI as well as OpenMP. The former makes use of a static domain decomposition and a distributed essential tree. Two new collision detection modules based on a plane-sweep algorithm are also implemented. The performance of the plane-sweep algorithm is superior to a tree code for simulations in which one dimension is much longer than the other two and in simulations which are quasi-two dimensional with less than one million particles. In this work, we discuss the different algorithms implemented in REBOUND, the philosophy behind the code's structure as well as implementation specific details of the different modules. We present results of accuracy and scaling tests which show that the code can run efficiently on both desktop machines and large computing clusters.
1. Introduction
REBOUND is an open-source collisional N-body code intended for planetary rings and other astrophysical or interdisciplinary problems. It also supports statistical collision-frequency estimates and classical collision-less N-body simulations.
- Motivation: REBOUND addresses collisional dynamics problems for which the authors state no publicly available code was capable of solving the described cases.The code was released under the GPLv31 open-source license.
- Applications: Collisional N-body simulations model systems such as planetary rings, where collision and orbital time-scales can be comparable and self-gravity matters.Planetary-ring simulations are usually performed using the shearing-sheet approximation.
- Applications: Collisions provide the dissipative mechanism for forming planetesimals from gravitationally bound swarms of boulders.
- Applications: REBOUND can approximate collision frequencies with super-particles when individual collisions need not be modeled exactly, including transitional and debris discs.
- Applications: The code also handles collision-less classical N-body problems involving test particles and massive particles.A symplectic mixed-variable integrator can follow both particle types.
2. Overview of the code structure
REBOUND combines portable C99 implementation with modular problem construction, configurable domain and boundary handling, and parallelization support. Users select and adapt modules through problem-specific makefiles and symbolic links.
- Implementation: REBOUND is written in ISO C99, runs on POSIX platforms, and requires no external libraries in its simplest form.OpenGL, GLUT, and LIBPNG support optional visualization and screenshot features.
- Parallelization: OpenMP supports shared-memory parallelization, MPI supports distributed-memory systems, and both can be combined in hybrid execution.
- Modularity: The code is organized around user-selectable gravity, collision, boundary, and integration modules that can be replaced or newly implemented with minimal effort.
- Modularity: Problem makefiles select modules through symbolic links, allowing different projects to use different configurations without changing the source directory.New modules can be copied into a problem directory, modified, and linked through its makefile.
- Domain and boundaries: The computational domain uses cubic boxes that can form elongated regions and support static MPI domain decomposition.Tree-based gravity or collision detection uses one tree structure per root box.
- Domain and boundaries: Open, periodic, and shear-periodic boundary conditions support particle removal, ghost-box periodicity, and shearing-sheet simulations, respectively.
3. Integrators
REBOUND provides second-order integrators within a Drift-Kick-Drift framework, including leap-frog, Wisdom-Holman, and the symplectic epicycle integrator. Their formal symplectic properties and practical performance depend on the force model, time-step, and problem geometry.
- General framework: REBOUND’s integrators are second-order accurate and symplectic, but approximated self-gravity, collisions, or velocity-dependent forces formally break symplecticity.
- General framework: All integrators use a Drift-Kick-Drift scheme whose Hamiltonian split determines the implementation of the three sub-steps.The split commonly separates kinetic and potential Hamiltonians.
- Time-stepping: REBOUND uses one shared time-step, which is fixed by default but can be made variable through implementation changes.Variable time-steps generally break the integrator’s symplectic nature.
- Time-stepping: Users must choose a sufficiently small time-step because multiple collisions or strongly curved trajectories can affect results.The paper recommends checking numerical convergence with different time-steps when studying new parameter regimes.
- Leap-frog: Leap-frog is a second-order symplectic integrator for non-rotating frames that advances positions and velocities through three simple Euler sub-steps.
- Wisdom-Holman Mapping: The Wisdom-Holman mapping solves two-body Keplerian motion accurately during drift and is suited to central 1/r potentials with small perturbations.It is slower than leap-frog because Kepler’s equation is solved iteratively for every particle at every time-step.
- Symplectic Epicycle Integrator: The symplectic epicycle integrator applies Hill’s approximation in a rotating frame and has similar Kepler-potential properties to Wisdom-Holman mapping.Its error after one time-step scales as the third power of the time-step times the gravitational-to-Coriolis force ratio.
- Symplectic Epicycle Integrator: SEI is as fast as a standard non-symplectic integrator and allows the radial and vertical epicyclic frequencies to be specified independently.This can enhance mid-plane particle density in planetary-ring simulations to represent self-gravity effects.
4. Gravity
REBOUND provides direct and Barnes–Hut gravity modules, with octree approximations controlled by an opening-angle criterion and parallelized through static domain decomposition and essential trees.
- Direct summation computes gravitational accelerations with O(N · Nactive) performance and is efficient mainly when Nactive ≲102.
- Barnes–Hut reduces computation from O(N2) to O(N log N) by approximating distant particle groups using their total mass and center of mass.
- The octree is built from cubic cells, recursively subdivided into up to eight sub-cells, with particles stored in leaf nodes.
- For force evaluation, cells are approximated when θ = w/R is below θcrit; otherwise, sub-cells are opened to balance accuracy and speed.
- Quadrupole expansions improve accuracy at fixed θcrit but increase computational cost.
- MPI parallelization uses static domain decomposition and essential trees containing remote cells that may be accessed during force calculations.
- The parallelization is efficient for homogeneous particle distributions with more than a few thousand particles per node; otherwise, communication dominates runtime.
5. Collisions
REBOUND provides multiple collision-detection modules, including direct, octree-based, and plane-sweep searches, with trade-offs in accuracy, complexity, and geometric applicability. The plane-sweep method is especially efficient for low-dimensional systems where the swept-list size remains small.
- Collision detection: REBOUND's collision modules may miss collisions or detect false collisions because they approximate trajectories or require particles to overlap.Straight-line trajectory approximations affect the sweep methods, while direct and tree methods detect collisions only for overlapping particles.
- Collision model: Collisions use a free-slip hard-sphere model that conserves momentum and energy while allowing constant or velocity-dependent normal restitution.Particle spin is not supported.
- Direct nearest-neighbor search: The direct nearest-neighbor search checks every particle pair once per time-step and therefore scales as O(N^2), making it suitable mainly for small systems.For most cases, the octree search is faster.
- Octree: The octree collision search identifies overlapping particles through cell-opening criteria and scales as O(N log(N)) for all N particles.It can be parallelized with OpenMP and MPI and can share a tree structure with gravity calculations.
- Plane-sweep algorithm: The plane-sweep modules move a conceptual plane along x or azimuthal angle and track intersected trajectories to identify candidate collisions.The x-sweep also applies to shearing sheets, while the φ-sweep targets narrow rings in global simulations.
- Plane-sweep algorithm: One plane-sweep time-step costs O(N · NSWEEPL) and is efficient only when NSWEEPL ≲ log(N); experiments find it can outperform octree searches by many orders of magnitude when NSWEEPL is small.This favors low-dimensional systems with elongated or quasi-two-dimensional geometry.
6. Test problems
REBOUND’s tests verify force accuracy, conservation, long-term Solar System integration, and planetary-ring viscosity calculations. The results generally agree with exact calculations or previous studies, while identifying a collision-treatment discrepancy for permanent particle clumps.
- Test problems: The test suite evaluates tree-code accuracy, conservation of energy and momentum, long-term Wisdom-Holman integration, and planetary-ring viscosity.These tests cover self-gravity, collisions, and symplectic integration across several physical regimes.
- Force accuracy: One and two orders of magnitude: quadrupole force errors are lower than monopole errors at θcrit ∼0.5 and θcrit ∼0.1, respectively.The quadrupole expansion costs more CPU time at fixed θcrit but is faster for θcrit ≲1 at fixed accuracy.
- Conservation tests: Perfectly elastic collisions conserve total momentum and energy up to machine precision for every collision-detection algorithm.The test uses a non-rotating periodic box with non-gravitating collisional particles.
- Long-term integration: 200 Myr integrations reproduce Pluto’s two libration frequencies, 3.8 Myr and 34 Myr, in agreement with Applegate et al. (1986).The integration used four massive planets, Pluto as a test particle, and a 40-day timestep; runtime was under two hours on one CPU.
- Planetary-ring viscosity: The planetary-ring viscosity results are in good agreement with previous results, but the collisional component is slightly higher when permanent particle clumps form.The discrepancy is attributed to different collision treatment and ambiguity in defining collisional viscosity when particles remain in contact.
7. Scaling
Scaling tests show efficient MPI and OpenMP parallelization, with performance depending on particles per processor and problem geometry. Collision algorithms have distinct scaling regimes, making algorithm choice important for the simulation configuration.
- Strong scaling: MPI strong scaling is linear at small processor counts but degrades below approximately 2000 particles per processor because essential-tree communication increases.The tests used 12.5k, 50k, 200k, and 800k particles.
- Strong scaling: One million particles can fully utilize 64 processor cores in the MPI scaling test.This result was obtained with the shearing-sheet gravity and collision tree modules.
- Weak scaling: Weak-scaling runtime follows O(Npp log(Npp k)) for k processors, while gravity runtime and communication increase logarithmically with total particle number.Collision communication per processor does not increase because only direct neighboring nodes are evaluated.
- OpenMP and MPI: OpenMP outperforms MPI when particles per node are few and communication dominates, whereas MPI is faster for very large simulations with Npp ≳5000.OpenMP can accelerate communication-bound MPI runs and multi-core desktop simulations.
- Collision detection: Collision-detection scaling depends strongly on geometry: plane-sweep scales as O(N^1.5) at fixed aspect ratio and O(N) in elongated domains, while tree scaling can reach O(N^2).The tree method scales as O(N log(N)) for constant-aspect-ratio tests but becomes effectively direct in elongated boxes.
- Collision detection: The plane-sweep module can be faster by many orders of magnitude in non-self-gravitating simulations whose size varies mainly along one dimension.The collision algorithm should therefore be selected according to the problem geometry.
8. Summary
REBOUND is an open-source, modular N-body code designed for collisional dynamics and adaptable to other applications. Its parallelized gravity and collision modules, together with plane-sweep collision detection, support efficient simulations across desktop machines and large clusters.
- REBOUND is a freely available, open-source multi-purpose N-body code for collisional dynamics.
- Its modular structure lets users choose or implement numerical integrators, boundary conditions, gravity solvers, and collision algorithms.
- MPI and OpenMP parallelize the octree self-gravity and collision-detection modules, which run efficiently on multi-core desktops and large clusters.
- Plane-sweep collision detection scales linearly in effectively low-dimensional systems and outperforms tree searches for elongated or narrow domains.It is also superior for quasi-two-dimensional simulations with fewer than about one million particles.
- Three symplectic integrators support rotating and non-rotating frames, including Wisdom-Holman mapping for exact long-term orbit integrations.
- The implemented methods are expected to apply beyond astrophysics, including molecular dynamics and granular flows.