Particular
Particular is a crate providing a simple way to simulate N-body gravitational interaction of particles in Rust.
Change log
Goals
The main goal of this crate is to provide users with a simple API to set up N-body gravitational simulations that can easily be integrated into existing game and physics engines. Thus it does not concern itself with numerical integration or other similar tools and instead only focuses on the acceleration calculations.
Particular is also built with performance in mind and provides multiple ways of computing the acceleration between particles.
Computation algorithms
There are currently 2 algorithms used by the available compute methods: 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 available compute methods.
Using Particular
Particular consists of two "modules", one that takes care of the abstraction of the computation of the gravitational forces between bodies for different floating-point types and dimensions, and one that facilitates usage of that abstraction for user-defined andnon-user-defined types. For most simple use cases, the latter is all that you need to know about.
Simple usage
The Particle trait provides the main abstraction layer between the internal representation
of the position and mass of an object in N-dimensional space and external types by defining
methods to retrieve a position and a gravitational parameter.
These methods respectively return an array of scalars and a scalar, which are converted using
the point_mass method to interface with the underlying algorithm implementations.
Implementing the Particle trait
When possible, it can be useful to implement Particle on a type.
Deriving
Used when the type has fields named position and mu:
#[derive(Particle)]
#[dim(3)]
struct Body {
position: Vec3,
mu: f32,
// ...
}
Manual implementation
Used when the type does not directly provide a position and a gravitational parameter.
struct Body {
position: Vec3,
mass: f32,
// ...
}
impl Particle for Body {
type Array = [f32; 3];
fn position(&self) -> [f32; 3] {
self.position.into()
}
fn mu(&self) -> f32 {
self.mass * G
}
}
If you can't implement Particle on a type, you can use the fact that it is implemented for
tuples of an array and its scalar type instead of creating an intermediate type.
let particle = ([1.0, 1.0, 0.0], 5.0);
assert_eq!(particle.position(), [1.0, 1.0, 0.0]);
assert_eq!(particle.mu(), 5.0);
Computing and using the gravitational acceleration
In order to compute the accelerations of your particles, you can use the accelerations method
on iterators, passing in a mutable reference to a ComputeMethod of your choice. It returns
the acceleration of each iterated item, preserving the original order.
Because it collects the mapped particles in a ParticleReordered in order to optimise the
computation of forces of massless particles, this method call results in one additional
allocation. See the advanced usage section for information on how to opt out.
When the iterated type implements Particle
for (acceleration, body) in bodies.iter().accelerations(&mut cm).zip(&mut bodies) {
body.velocity += Vec3::from(acceleration) * DT;
body.position += body.velocity * DT;
}
When the iterated type doesn't implement Particle
// Items are a tuple of a velocity, a position and a mass.
// We map them to a tuple of the positions as an array and the mu,
// since this implements `Particle`.
let accelerations = items
.iter()
.map(|(_, position, mass)| (*position.as_array(), *mass * G))
.accelerations(&mut cm);
for (acceleration, (velocity, position, _)) in accelerations.zip(&mut items) {
*velocity += Vec3::from(acceleration) * DT;
*position += *velocity * DT;
}
Advanced usage
In some instances the iterator abstraction provided by particular might not be flexible enough. For example, you might need to access the tree built from the particles for the Barnes-Hut algorithm, want to compute the gravitational forces between two distinct collections of particles, or both at the same time.
The PointMass type
The underlying type used in storages is the PointMass, a simple representation in
N-dimensional space of a position and a gravitational parameter. Instead of going through a
ComputeMethod, you can directly use the different generic methods available to compute the
gravitational forces between PointMasses, with variants optimised for scalar and simd types.
Example
use particular::math::Vec2;
use storage::PointMass;
let p1 = PointMass::new(Vec2::new(0.0, 1.0), 1.0);
let p2 = PointMass::new(Vec2::new(0.0, 0.0), 1.0);
let softening = 0.0;
assert_eq!(p1.force_scalar::<false>(p2.position, p2.mass, softening), Vec2::new(0.0, -1.0));
Storages and built-in ComputeMethod implementations
Storages are containers that make it easy to apply certain optimisation or algorithms on collections of particles when computing their gravitational acceleration.
The ParticleSystem storage defines an affected slice of particles and a massive storage,
allowing algorithms to compute gravitational forces the particles in the massive storage exert
on the affected particles. It is used to implement most compute methods, and blanket
implementations with the other storages allow a ComputeMethod implemented with
ParticleSliceSystem or ParticleTreeSystem to also be implemented with the other
storages.
The ParticleReordered similarly defines a slice of particles, but stores a copy of them in a
ParticleOrdered. These two storages make it easy for algorithms to skip particles with no
mass when computing the gravitational forces of particles.
Example
use particular::math::Vec3;
let particles = vec![
// ...
];
// Create a `ParticleOrdered` to split massive and massless particles.
let ordered = ParticleOrdered::from(&*particles);
// Build a `ParticleTree` from the massive particles.
let tree = ParticleTree::from(ordered.massive());
// Do something with the tree.
for (node, data) in std::iter::zip(&tree.get().nodes, &tree.get().data) {
// ...
}
let bh = &mut sequential::BarnesHut { theta: 0.5 };
// The implementation computes the acceleration exerted on the particles in
// the `affected` slice.
// As such, this only computes the acceleration of the massless particles.
let accelerations = bh.compute(ParticleSystem {
affected: ordered.massless(),
massive: &tree,
});
Custom ComputeMethod implementations
In order to work with the highest number of cases, built-in compute method implementations may
not be the most appropriate or optimised for your specific use case. You can implement the
ComputeMethod trait on your own type to satisfy your specific requirements but also if you
want to implement other algorithms.
Example
use particular::math::Vec3;
struct MyComputeMethod;
impl ComputeMethod<ParticleReordered<'_, Vec3, f32>> for MyComputeMethod {
type Output = Vec<Vec3>;
#[inline]
fn compute(&mut self, storage: ParticleReordered<Vec3, f32>) -> Self::Output {
// Only return the accelerations of the massless particles.
sequential::BruteForceScalar.compute(ParticleSystem {
affected: storage.massless(),
massive: storage.massive(),
})
}
}
License
This project is licensed under either of Apache License, Version 2.0 or MIT license, at your option.
Contribution
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.