Gravity Simulator

Newtonian N-body gravity simulator accelerated with C library

This is a student project developed for learning purpose. Other packages such as REBOUND are recommended for more functionality.

Features:

• Ten integrators including WHFast and IAS15
• Barnes-Hut algorithm
• GravitySimulator API (prototype)
• CUDA implementation for the pairwise acceleration function
• Multiple sample projects

TODOs (if I have time):

• Barnes-Hut algorithm with CUDA

Documentation

• Quick Start
  • Python version
  • Installation
  • Some notes
• Running the program in terminal
• Sample projects
• GravitySimulator API
• Default systems
• Integrators
  • Simple methods
  • Embedded Runge-Kutta methods
  • IAS15
  • WHFast
• Making changes to saved systems
• Saving the results
• Output animations in .gif
• Compensated summation
• Feedback and Bugs
• Data Sources
• References
• Acknowledgement

Quick Start

Python version

This program requires Python version 3.10 or higher.

Installation

Download the source files, or clone this repository by running the following command in terminal:

git clone https://github.com/alvinng4/Gravity-Simulator

Install the required packages by

pip install .

If the installation is not successful, install the following packages manually:

matplotlib==3.8.3
numpy==1.26.4
rich==13.7.1

If the program failed to load the C library, you may choose the NumPy option (slow), or choose to recompile the C library. Simply go to the src folder and run

make [CC=gcc] [USE_CUDA=1]

Then, move the c_lib.dylib, c_lib.dll or c_lib.so file into the gravity_sim folder (One of them will be generated depending on your operation system).

mv c_lib.dylib ../gravity_sim
mv c_lib.dll   ../gravity_sim
mv c_lib.so    ../gravity_sim

If you wish to use CUDA GPU acceleration, you will also need to recompile the library with the USE_CUDA=1 flag, which requires nvcc to be installed in your system.

Some notes

• This project offers two user-interface: API and CLI
• The default unit for this project is solar masses, AU and days, with $G = 0.00029591220828411956 \text{ M}_\odot^{-1} \text{ AU}^3 \text{ day}^{-2}$. It is possible to change this value in the API by changing system.G.
• Animations, simulation results, etc. will be stored to gravity_sim/result by default, unless a file path is specified.
• Complex animations like the asteroid belt cannot be done solely with the API functions. Sample scripts are provided in this repository (See Sample projects)
• Check the examples folder for API tutorial and sample projects
• For WHFast, features including CUDA acceleration, compensated_summation, barnes-hut and resume_simulation are not available due to implementation difficulties and bugs.

Running the program in terminal

Once you have downloaded the source files, navigate to the source directory in terminal and run

python gravity_sim [-n|--numpy]

-n, --numpy: run the program with NumPy instead of C library

Sample projects

Some projects are done with the API. The scripts are stored at the examples folder.

Simulating the solar system for 1 million years

Asteroid belt animation

Formation of Kirkwood gap

Videos:

• https://www.youtube.com/watch?v=AEyjIF-8zT0
• https://www.youtube.com/watch?v=jHLLr7ACvDQ

Galaxy collision

Note: This simulation is for demonstration only. The initial conditions are very different from actual galaxies.

Videos:

• https://www.youtube.com/watch?v=g45o0Vf4zDk
• https://www.youtube.com/watch?v=OdfgFYLYcVw

GravitySimulator API

You may import the GravitySimulator API from gravity_sim to perform gravity simulation. See examples/tutorial.ipynb or Sample projects for some example usage. If your computer cannot render jupyter notebook (files end with .ipynb), you can view them on github.

from gravity_sim import GravitySimulator

grav_sim = GravitySimulator()
grav_sim.integration_mode = "c_lib"  # "numpy" is also available

system = grav_sim.create_system()
grav_sim.set_current_system(system)
system.load("solar_system")

grav_sim.launch_simulation(
    integrator="ias15",
    tf=grav_sim.years_to_days(1000.0),
    tolerance=1e-9,
    store_every_n=50,
)

grav_sim.plot_2d_trajectory()
grav_sim.plot_rel_energy()
grav_sim.save_results()

launch_simulation

launch_simulation() is the main method for launching the simulation.

Parameter	Type	Default	Description
integrator	str	Required	Name of the integrator
tf	float	Required	Simulation time (days)
dt	float	None	Time step (days)
tolerance	float	None	Tolerance for adaptive step integrators
store_every_n	int	1	Store results every n steps
acceleration_method	str	pairwise	Method for calculating accelerations
storing_method	str	default	Method for storing simulation results
flush_path	str	None	Path to flush intermediate results.
no_progress_bar	bool	False	If True, disables the progress bar.
no_print	bool	False	If True, disables some printing to console
softening_length	float	0.0	Softening length for acceleration calculation
**kwargs	dict	-	Additional keyword arguments.

integrators

euler, euler_cromer, rk4, leapfrog, rkf45, dopri, dverk, rkf78, ias15, whfast

acceleration_method

• pairwise
  • Brute force pairwise calculations for gravitational acceleration
  • Time complexity: $O(N^2)$
  • Variants: pairwise_cuda, pairwise_float_cuda, pairwise_float_comp_sum_cuda (Not available for WHFast)
    • Parallelized calculation with CUDA GPU acceleration.
    • pairwise_float_cuda uses single precision to provide furthur speed up with a loss of precision
    • pairwise_float_comp_sum_cuda uses compensated summation to reduce precision loss (See Compensated summation)
• massless
  • Similar to pairwise, but seperate the calculations for massive and massless particles
  • Time complexity: $O(M^2 + MN)$, where $M$ and $N$ are the number of massive and massless particles respectively
• barnes-hut
  • Calculate gravitational acceleration with barnes-hut algorithm
  • Time complexity: $O(N \log{N})$
  • **kwargs: barnes_hut_theta
    • Threshold for Barnes-hut algorithm, default = 0.5

storing_method

• default
  • Store solutions directly into memory
• flush
  • Flush intermediate results into a csv file in gravity_sim/results to reduce memory pressure.
• no_store
  • To not store any result, typically used for benchmarking.

Default systems

Some systems are available by default.

System	Description
circular_binary_orbit	A circular orbit formed by two stars
eccentric_binary_orbit	An eccentric orbit formed by two stars
3d_helix	An upward helix consists of three stars
sun_earth_moon	The Sun, Earth, and Moon system
figure-8	A "figure-8" orbit involving three stars
pyth-3-body	Three stars arranged in a triangle with length ratios of 3, 4, and 5. It is a highly chaotic orbit with close encounters that can be used to test the difference between fixed and variable step size integrators.
solar_system	Solar System with the Sun and the planets
solar_system_plus	solar_system with the inclusion of Pluto, Ceres, and Vesta

Integrators

Simple methods

Below are four simple fixed step size methods to simulate the system with a given step size $\text{d}t$.

Simple methods
Euler
Euler Cromer
Fourth Order Runge-Kutta (RK4)
Leapfrog

Embedded Runge-Kutta methods

Embedded RK methods are adaptive methods that decides the step size automatically based on the estimated error. It can resolve close encounters but fail to conserve energy over long time scele.

Embdedded Runge-Kutta methods	Recommended tolerance*
Runge–Kutta–Fehlberg 4(5)	$10^{-8}$ to $10^{-14}$
Dormand–Prince method (DOPRI) 5(4)	$10^{-8}$ to $10^{-14}$
Verner's method (DVERK) 6(5)	$10^{-8}$ to $10^{-14}$
Runge–Kutta–Fehlberg 7(8)	$10^{-4}$ to $10^{-8}$

*For reference only

IAS15

IAS15 (Implicit integrator with Adaptive time Stepping, 15th order) is a highly optimized integrator with extremely high accuracy. It is the default method for this project.

The recommended tolerance* is $10^{-9}$. Since the integrator is 15th order, changing the tolerance results in little improvement in performance, but a huge penalty in accuracy. Therefore, it is not recommended to change this tolerance.

*For reference only

WHFast

WHFast is a second order symplectic method with fixed step size, which conserves energy over long integration period. This integrator cannot resolve close encounter.

**kwargs for WHFast

Argument	Description	Default Value
kepler_tol	Tolerance in solving the Kepler's equation	$10^{-12}$
kepler_max_iter	Maximum number of iterations in solving Kepler's equation	500
kepler_auto_remove	Integer flag to indicate whether to remove objects that failed to converge in Kepler's equation	False
kepler_auto_remove_tol	Tolerance for removing objects that failed to converge in Kepler's equation	$10^{-8}$

Warning

When using WHFast, the order of adding objects matters. Since WHFast use Jacobi coordinate, we must add the inner object first, followed by outer objects relative to the central star. For convenience, you may also add the objects in any order, then call system.sort_by_distance(primary_object_name) or system.sort_by_distance(primary_object_index)

Making changes to saved systems

If you wish to make any changes to saved systems, you can access the file at

gravity_simulator/gravity_sim/customized_systems.csv

The data follow the format

Name, Gravitational constant, Number of objects, m1, ..., x1, y1, z1, ..., vx1, vy1, vz1, ...

You may also load the system in API and then save the system after making the changes.

Saving the results

If you saved the results, the numerical data will be stored in the following folder:

Gravity-Simulator/gravity_sim/results

The file will starts with the metadata which starts with #. Missing information will be saved as None. More rows may be added in the future.

Below is an example:

# Data saved on (YYYY-MM-DD): 2024-07-27
# System Name: solar_system
# Integrator: IAS15
# Number of objects: 9
# Gravitational constant: 0.00029591220828411956
# Simulation time (days): 73048.4378
# dt (days): None
# Tolerance: 1e-09
# Data size: 64596
# Store every nth point: 1
# Run time (s): 1.4667500000214204
# masses: 1.0 1.6601208254589484e-07 2.447838287796944e-06 3.0034896154649684e-06 3.2271560829322774e-07 0.0009547919099414248 0.00028588567002459455 4.36624961322212e-05 5.151383772654274e-05

Then, the actual data will be saved in the default unit (solar masses, AU and days), and follow this format:

time, dt, total energy, x1, y1, z1, ... vx1, vy1, vz1, ...

The saved data file can be read by the program. Even if the metadata is corrupted or missing, the program can still read the data file, although some information could be missing.

Output animations in .gif

You may output the simple trajectory animations in 2D / 3D in .gif. The output file would be stored in gravity_sim/result.

Available parameters:

Parameter	Description
FPS	Frames per second
Animation length	Desired length of the output gif
plot_every_nth_point	File name without extension
dpi	The resolution of the animation (dots per inch)
is_dynamic_axes	Rescale the axes limit dynamically
axes_lim	An array for the axes limits
is_maintain_fixed_dt	Attempt to maintain fixed step size with variable time step data
traj_len	Length of the trajectory in number of data points

Some parameters are not listed here due to limit of space.

Compensated summation

A method known as compensated summation [1], [4] is implemented for all integrators EXCEPT WHFast:

When we advance our system by $\text{d}t$, we have

$x_{n+1} = x_n + \delta x$

Since $\delta x$ is very small compared to $x_n$, many digits of precision will be lost. By compensated summation, we keep track of the losing digits using another variable, which allows us to effectively eliminates round off error with very little cost.

However, for WHFast, the improvement is little but takes 10% longer run time. Therefore, it is excluded from this method.

Feedback and Bugs

If you find any bugs or want to leave some feedback, please feel free to let me know by opening an issue or sending an email to [email protected].

Data Sources

The solar system positions and velocities data at 1/Jan/2024 are collected from the Horizons System [2]. Gravitational constant, and masses of the solar system objects are calculated using the data from R.S. Park et. al. [3].

References

1. E. Hairer, C. Lubich, and G. Wanner, "Reducing Rounding Errors" in Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Springer, 2004, pp. 272-274.
2. Horizons System, Jet Propulsion Laboratory, https://ssd.jpl.nasa.gov/horizons/
3. R. S. Park, et al., 2021, “The JPL Planetary and Lunar Ephemerides DE440 and DE441”, https://ssd.jpl.nasa.gov/doc/Park.2021.AJ.DE440.pdf, Astronomical Journal, 161:105.
4. H. Rein, and D. S. Spiegel, 2014, "IAS15: A fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits", Monthly Notices of the Royal Astronomical Society 446: 1424–1437.

Acknowledgement

The integrators in this project were developed with great assistance from the following book:

• J. Roa, et al. Moving Planets Around: An Introduction to N-Body Simulations Applied to Exoplanetary Systems, MIT Press, 2020