Particular is a crate providing a simple way to compute N-body interaction of particles in Rust.
Please note that this branch is for development purposes and may not represent the latest stable
release of the library. For the most recent stable version, refer to the
latest branch.
The main goal of this crate is to provide users with a simple API to set up N-body simulations that can easily be integrated into existing game and physics engines. Thus it does not include anything related to numerical integration or other similar tools and instead only focuses on the calculations of the interactions between particles.
In order to do this, particular provides multiple algorithms to compute the
interactions between particles, allowing users to choose the one that best fits their needs in
terms of performance and accuracy.
There are currently 2 algorithms: Brute-force and Barnes-Hut.
Generally speaking, the Brute-force algorithm is more accurate, but slower. The Barnes-Hut
algorithm allows trading accuracy for speed by increasing the theta parameter.
You can see more about their relative performance here.
Particular uses rayon for parallelization and
wgpu for GPU computation.
Enable the respective parallel and gpu features to access the relevant algorithms.
At its core, Particular is a simple set of traits designed to simplify the implementation of
different algorithms to compute interactions between particles. These interactions can be of any
nature, but a common example is the gravitational interaction between particles, which is
provided by the gravity module and the main focus of the library.
This module provides implementations using popular vector math libraries such as glam and
nalgebra. You you will need to enable the corresponding features to use them.
The Position and Mass traits provide methods to retrieve the position, mass and
gravitational parameter of a particle. Implementing these traits for a type will result in
multiple blanket implementations that will allow it to be used with all available algorithms for
computing different interactions.
When the type has fields named position and mass or mu, you can derive these traits. An
optional attribute #[G = ...] defining the gravitational constant can be added on the field.
If the attribute is missing, the value of the gravitational constant is 1.0.
use particular::prelude::*;
use glam::DVec3;
#[derive(Position, Mass)]
struct Planet {
velocity: DVec3,
position: DVec3,
#[G = 6.67430e-11]
mass: f64,
}If you can't implement these traits, you can use the fact that it is implemented for tuples of a position and a mass or gravitational parameter instead of creating an intermediate type.
use particular::prelude::*;
use glam::DVec3;
let particle = (DVec3::new(-1.0, 0.0, 1.0), 5.0);
assert_eq!(particle.position(), DVec3::new(-1.0, 0.0, 1.0));
assert_eq!(particle.mass(), 5.0);
assert_eq!(particle.mu(), 5.0);In order to compute the gravitational interaction between particles, you can use the various
methods implemented on storages of particles.
Note that most algorithms return iterators that borrow the particles, and as such, rust's
borrowing rules will prevent you from mutating these particles as you iterate the computed
interactions, even if the interaction does not use the fields you are mutating. You can either
immediately collect the interactions or use separate collections for the properties you need to
mutate.
use particular::prelude::*;
use particular::gravity::newtonian::Acceleration;
let accelerations: Vec<_> = planets.brute_force(Acceleration::checked()).collect::<Vec<_>>();
for (planet, acceleration) in planets.iter_mut().zip(accelerations) {
planet.velocity += acceleration * DT;
planet.position += planet.velocity * DT;
}You can alternatively use indices, but note that this requires making sure the positions don't change as the interaction is computed or you will get unexpected behaviour.
use particular::prelude::*;
use particular::gravity::newtonian::AccelerationSoftened;
for i in 0..planets.len() {
let between = Between(&planets[i], planets.as_slice());
let acceleration = between.brute_force(AccelerationSoftened::checked(1e-3));
planets[i].velocity += acceleration * DT;
}
for planet in planets.iter_mut() {
planet.position += planet.velocity * DT;
}Storages and built-in Interaction implementations
Storages are containers that make it easy to apply certain optimisation or algorithms on collections of particles when computing their gravitational interaction.
Between is a tuple struct of two objects where the first one is conventionally the affected
object and the second one is the affecting object. These objects can be particles or storages
of particles. Between is the underlying type used to implement most algorithms allowing
trivial implementation of other storage types.
The Reordered storage defines a slice of particles and stores a copy of them in an
Ordered storage. These two storages make it easy for algorithms to skip particles that do
not affect other particles when computing interactions. A RootedOrthtree is a tree structure
(equivalent to a quadtree in 2D, an octree in 3D etc.) that stores particles to be used in
Barnes-Hut algorithms.
use particular::prelude::*;
use particular::gravity::newtonian::{Acceleration, TreeData};
use glam::Vec3;
// Particles here are simple gravitational fields.
let particles: Vec<GravitationalField<Vec3, f32>> = vec![
// ...
];
// Function to determine if a field will affect particles.
fn is_massive(p: &GravitationalField<Vec3, f32>) -> bool {
p.m != 0.0
}
// Create an `Ordered` storage to split massive and massless particles.
let ordered = Ordered::new(&particles, is_massive);
// We can build a `RootedOrthtree` from just the massive particles. Note that this is done
// automatically when computing over an [`Ordered`] or [`Reordered`] storage for the Barnes-Hut
// algorithm.
let tree = RootedOrthtree::new(
ordered.affecting(),
|p| p.position().to_array(),
|slice| GravitationalField::centre_of_mass(slice.iter().cloned()),
);
// Do something with the tree.
for (node, data) in std::iter::zip(&tree.get().nodes, &tree.get().data) {
// ...
}
// Compute the acceleration the massive particles in a tree exert on the massless particles.
let between = Between(ordered.non_affecting(), &tree);
let accelerations = between.barnes_hut(0.5, Acceleration::checked()).collect::<Vec<_>>();Custom Interaction implementations
Implementing Interaction for YourInteraction allows
it to be used with the CPU brute-force algorithms. Other algorithms may require you to implement
other traits to be used, namely SimdInteraction, ReduceSimdInteraction,
BarnesHutInteraction and TreeInteraction. Refer to the documentation of the specific
algorithms for more information.
use particular::prelude::*;
use glam::DVec3;
#[derive(Clone, Copy)]
struct Body {
position: DVec3,
mass: f64,
}
#[derive(Clone)]
struct GravitationalForce(pub f64);
impl Interaction<Between<&Body, &Body>> for GravitationalForce {
type Output = DVec3;
fn compute(&mut self, Between(affected, affecting): Between<&Body, &Body>) -> DVec3 {
if affected.position == affecting.position {
return DVec3::ZERO;
}
let r = affecting.position - affected.position;
let l = r.length_squared();
let f = self.0 * affected.mass * affecting.mass / (l * l.sqrt());
r * f
}
}
// Distance: AU, Mass: M☉ (solar mass)
const G: f64 = 4.0 * std::f64::consts::PI * std::f64::consts::PI;
let sun = Body::new(DVec3::ZERO, 1.0);
let earth = Body::new(DVec3::new(1.0, 0.0, 0.0), 3.0027e-6);
let jupiter = Body::new(DVec3::new(5.2, 0.0, 0.0), 0.000954588);
let solar_system = [sun, earth, jupiter];
let mut forces = solar_system.brute_force(GravitationalForce(G));
let sun_earth = GravitationalForce(G).compute(Between(&sun, &earth));
let sun_jupiter = GravitationalForce(G).compute(Between(&sun, &jupiter));
let earth_jupiter = GravitationalForce(G).compute(Between(&earth, &jupiter));
assert_eq!(forces.next(), Some(sun_earth + sun_jupiter));
assert_eq!(forces.next(), Some(-sun_earth + earth_jupiter));
assert_eq!(forces.next(), Some(-sun_jupiter - earth_jupiter));This project is licensed under either of Apache License, Version 2.0 or MIT license, at your option.
Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in this project by you, as defined in the Apache 2.0 license, shall be dual licensed as above, without any additional terms or conditions.