1 Commits

Author SHA1 Message Date
  Björn Steinbrink 21b2fcf8f0 Perf: Use struct-of-arrays instead of array-of-structs and precalculate s 8 months ago
2 changed files with 55 additions and 66 deletions
Split View
  1. +53
    -64
      src/barnes_hut.rs
  2. +2
    -2
      src/main.rs

+ 53
- 64
src/barnes_hut.rs View File

@ -7,34 +7,36 @@ const G: f64 = 0.00667408;
const THETA: f64 = 0.5;
pub struct Octree {
nodes: Vec<OctNode>
}
pub struct OctNode {
parent: isize, // -1 = no parent. TODO get rid of all the parent stuff - we never use it
children: [isize; 8], // TODO make this an Option<[usize; 8]>
mass: f64,
min_coord: Vector3D,
max_coord: Vector3D,
num_particles: usize,
avg_weighted_position: Vector3D, // Divide by mass for this to be accurate!
center_of_mass: Vector3D // Calculated once from avg_weighted_position
num_particles: Vec<usize>,
children: Vec<[isize; 8]>, // TODO make this an Option<[usize; 8]>
mass: Vec<f64>,
min_coord: Vec<Vector3D>,
max_coord: Vec<Vector3D>,
avg_weighted_position: Vec<Vector3D>, // Divide by mass for this to be accurate!
center_of_mass: Vec<Vector3D>, // Calculated once from avg_weighted_position
s: Vec<f64>,
}
impl Octree {
pub fn new(min: Vector3D, max: Vector3D) -> Self {
let mut nodes = Vec::new();
nodes.push(OctNode::new(min, max, -1)); // Create our root node
Octree { nodes: nodes }
Octree {
num_particles: vec![0],
children: vec![[-1; 8]],
mass: vec![0.0],
min_coord: vec![min],
max_coord: vec![max],
avg_weighted_position: vec![Vector3D::new(0.0, 0.0, 0.0)],
center_of_mass: Vec::new(),
s: vec![0.0],
}
}
pub fn construct_tree(&mut self, particles: Vec<RigidPoint>) {
self.add_particles_to_specific_node(particles, 0);
//calculate center of masses
for n in &mut self.nodes {
n.center_of_mass = n.avg_weighted_position / n.mass;
for (&mass, &awp) in self.mass.iter().zip(&self.avg_weighted_position) {
self.center_of_mass.push(awp / mass);
}
}
@ -46,39 +48,29 @@ impl Octree {
particle: &mut RigidPoint,
cur_node: usize,
delta: f64) {
let c_node = &self.nodes[cur_node];
if c_node.num_particles == 0 { return; }
let min_coord = c_node.min_coord;
let max_coord = c_node.max_coord;
if self.num_particles[cur_node] == 0 { return; }
let width = max_coord.x - min_coord.x;
let height = max_coord.y - min_coord.y;
let depth = max_coord.z - min_coord.z;
let s = if width > height { width } else { height };
let s = if depth > s { depth } else { s };
let s = self.s[cur_node];
let s_squared = s * s;
let dis_vec = self.nodes[cur_node].center_of_mass - particle.position;
let dis_vec = self.center_of_mass[cur_node] - particle.position;
let d_squared = dis_vec.dist_squared();
let theta_squared = THETA * THETA;
let ratio = s_squared / d_squared;
if c_node.num_particles == 1 || ratio < theta_squared {
if self.num_particles[cur_node] == 1 || ratio < theta_squared {
// println!("s/d = {}", ratio.sqrt());
//We can approximate
if d_squared < 0.1 { return; } // Worse, but higher-performance, method of smoothing
let other = &self.nodes[cur_node];
let dist = d_squared.sqrt();
let v = G * other.mass * delta / d_squared; // / (d_squared + 0.1).sqrt(); // Uncomment this to enable smoothing for short distances
let v = G * self.mass[cur_node] * delta / d_squared; // / (d_squared + 0.1).sqrt(); // Uncomment this to enable smoothing for short distances
particle.velocity += v * dis_vec / dist;
}
else {
//println!("A{}", self.nodes[cur_node].num_particles);
//Divide and conquer
for i in 0..8 {
self.barnes_hut_specific_node(particle, c_node.children[i] as usize, delta);
for &child in &self.children[cur_node] {
self.barnes_hut_specific_node(particle, child as usize, delta);
}
}
}
@ -89,24 +81,36 @@ impl Octree {
}
if particles.len() == 1 {
self.nodes[cur_node].mass += particles[0].mass;
self.nodes[cur_node].num_particles += 1;
self.nodes[cur_node].avg_weighted_position += particles[0].position * particles[0].mass;
self.mass[cur_node] += particles[0].mass;
self.num_particles[cur_node] += 1;
self.avg_weighted_position[cur_node] += particles[0].position * particles[0].mass;
return;
}
let center_point = (self.nodes[cur_node].min_coord + self.nodes[cur_node].max_coord) / 2.0;
let center_point = (self.min_coord[cur_node] + self.max_coord[cur_node]) / 2.0;
if self.nodes[cur_node].children[0] == -1 {
if self.children[cur_node][0] == -1 {
for i in 0..8 {
let idx = self.nodes.len();
let (min, max) = Octree::get_octant_bounding_box_from_id(self.nodes[cur_node].min_coord,
self.nodes[cur_node].max_coord,
let idx = self.num_particles.len();
let (min, max) = Octree::get_octant_bounding_box_from_id(self.min_coord[cur_node],
self.max_coord[cur_node],
center_point,
i);
self.nodes.push(OctNode::new(min, max, cur_node as isize));
self.nodes[cur_node].children[i] = idx as isize;
let width = max.x - min.x;
let height = max.y - min.y;
let depth = max.z - min.z;
let s = if width > height { width } else { height };
let s = if depth > s { depth } else { s };
self.num_particles.push(0);
self.children.push([-1; 8]);
self.mass.push(0.0);
self.min_coord.push(min);
self.max_coord.push(max);
self.avg_weighted_position.push(Vector3D::new(0.0, 0.0, 0.0));
self.s.push(s);
self.children[cur_node][i] = idx as isize;
}
}
else {
@ -118,9 +122,9 @@ impl Octree {
for particle in particles {
// Update ourself
self.nodes[cur_node].mass += particle.mass;
self.nodes[cur_node].num_particles += 1;
self.nodes[cur_node].avg_weighted_position += particle.position * particle.mass;
self.mass[cur_node] += particle.mass;
self.num_particles[cur_node] += 1;
self.avg_weighted_position[cur_node] += particle.position * particle.mass;
let octant_index = Octree::get_id_from_center(center_point, particle.position);
particle_octs[octant_index].push(particle);
@ -129,7 +133,7 @@ impl Octree {
//Recurse
for (i, particle_oct) in particle_octs.iter().enumerate() {
self.add_particles_to_specific_node(particle_oct.clone(),
self.nodes[cur_node].children[i] as usize);
self.children[cur_node][i] as usize);
}
}
@ -179,18 +183,3 @@ impl Octree {
(min_coord, max_coord)
}
}
impl OctNode {
fn new(min: Vector3D, max: Vector3D, parent: isize) -> Self {
OctNode {
parent: parent,
children: [-1, -1, -1, -1, -1, -1, -1, -1],
mass: 0.0,
min_coord: min,
max_coord: max,
num_particles: 0,
avg_weighted_position: Vector3D::new(0.0, 0.0, 0.0),
center_of_mass: Vector3D::new(0.0, 0.0, 0.0)
}
}
}

+ 2
- 2
src/main.rs View File

@ -245,8 +245,8 @@ fn calc_velocities(u: &mut Universe, delta: f64) {
println!("Barnes hut");
u.particles.par_iter_mut().for_each(|mut p| {
oct.barnes_hut(&mut p, delta);
u.particles.par_iter_mut().for_each(|p| {
oct.barnes_hut(p, delta);
});
println!("Done");


Loading…
Cancel
Save