Created
March 12, 2026 19:54
-
-
Save dineshadepu/031013c29975406e81a63c0c57368f9b to your computer and use it in GitHub Desktop.
Parallel iteration using Rayon
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| use rayon::prelude::*; | |
| pub struct Particles { | |
| /// Total number of particles in the simulation. | |
| pub n: usize, | |
| pub x: Vec<f64>, | |
| pub u: Vec<f64>, | |
| pub f: Vec<f64>, | |
| /// Mass of particles. | |
| pub m: Vec<f64>, | |
| /* ---------------- Material / physics properties ---------------- */ | |
| /// Density (mainly for SPH / continuum methods). | |
| pub rad: Vec<f64>, | |
| pub tng_ids: Vec<u32>, | |
| } | |
| impl Particles { | |
| /// Create a new particle container with `n` particles | |
| /// and `n_rb` rigid bodies. | |
| /// | |
| /// All fields are initialized to zero. | |
| /// User is expected to fill positions, masses, kinds, etc. | |
| pub fn new(n: usize, max_tng_cnts: usize) -> Self { | |
| Self { | |
| n, | |
| x: vec![0.0; n], | |
| u: vec![0.0; n], | |
| f: vec![0.0; n], | |
| m: vec![0.0; n], | |
| rad: vec![0.0; n], | |
| tng_ids: vec![u32::MAX; n * max_tng_cnts], | |
| } | |
| } | |
| } | |
| pub trait ForceModel { | |
| fn compute(&mut self, p: &mut Particles); | |
| } | |
| pub struct DEMForce { | |
| pub kn: f64, | |
| } | |
| impl DEMForce { | |
| pub fn new(kn: f64) -> Self { | |
| Self { kn } | |
| } | |
| } | |
| impl ForceModel for DEMForce { | |
| fn compute(&mut self, p: &mut Particles) { | |
| let kn = self.kn; | |
| let x = &p.x; | |
| let u = &p.u; | |
| let m = &p.m; | |
| let rad = &p.rad; | |
| let force = &mut p.f; | |
| let tng_ids = &mut p.tng_ids; | |
| let max_contacts = 8; | |
| force | |
| .par_iter_mut() | |
| .zip(p.tng_ids.par_chunks_mut(max_contacts)) | |
| .enumerate() | |
| .for_each(|(i, (fi, tng_ids_i))| { | |
| let xi = x[i]; | |
| let ri = rad[i]; | |
| let mut fx = 0.0; | |
| let mut cnt = 0; | |
| for j in 0..x.len() { | |
| let dx = xi - x[j]; | |
| let r2 = dx * dx; | |
| if r2 < 1e-24 { | |
| continue; | |
| } | |
| let rij = r2.sqrt(); | |
| let overlap = ri + rad[j] - rij; | |
| if overlap <= 0.0 { | |
| continue; | |
| } | |
| fx += kn * overlap * dx; | |
| if cnt < max_contacts { | |
| tng_ids_i[cnt] = j as u32; | |
| cnt += 1; | |
| } | |
| } | |
| *fi += fx; | |
| }); | |
| } | |
| } | |
| fn main() {} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment