I'm trying to parrallelize some particle-in-cell code, and part of the design involves mutating the contents of vectors by index at some points. I've done this by creating mutable pointers as in this SO question, where each time I'm using this, either the indices are unique or I'm ensuring there isn't concurrent access to a single index
I realize in the code below, there is an alternative way of designing the iterators that can avoid using the mutable pointers, so that I'm iterating over a tuple of position and keepers, but there is another function where for the life of me I can't think of a way to avoid it without hitting some major performance penalty.
this is the structure that these functions are implemented for.
pub struct Particles {
//receiver: ThingBuf,
pub velocities: Vec<na::Vector3<f64>>,
pub positions: Vec<na::Vector3<f64>>,
pub types: Vec<ParticleType>,
pub parent_cells: Vec<usize>,
keepers: Vec<bool>,
}
this is the pointer I've created:
#[derive(Copy, Clone)]
pub(crate) struct BoolPointer(pub(crate) *mut bool);
unsafe impl Send for BoolPointer {}
unsafe impl Sync for BoolPointer {}
this function should be modifying the keepers vec:
///cell membership goes like this:
/// 1. abs(x) starting from the origin,
/// - 0 and negative numbers are even
/// - positive numbers are odd
/// 2. y just goes up from the bottom,
/// 3. z from front to back, so
/// if you have a 2 x 4 x 3 grid it's members would be
/// 6 14 22 7 15 24
/// 4 12 20 5 13 21
/// 2 10 18 3 11 19
/// 0 8 16 1 9 17
pub fn index(
&mut self,
sample_sender: &SampleUpdateSender,
num_x: usize,
num_y: usize,
num_z: usize,
) {
let num_cells = num_x * num_y * num_z;
//assuming number of cells must be even in the x direction
let half_x = num_x / 2;
let dy: f64 = 2. / num_y as f64;
let dz: f64 = 2. / num_y as f64;
let Particles {
velocities,
positions,
types: _,
parent_cells,
keepers,
} = self;
let particle_count = velocities.len();
let cell_ptr = ParentCellPointer(parent_cells.as_mut_ptr());
let keep_ptr = BoolPointer(keepers.as_mut_ptr());
(0..particle_count)
.into_par_iter()
.chunks(4)
.map(move |chunk| {
chunk //iterate syncronously for chunk size
.iter()
.map(|i| {
let y_offset =
((positions[*i].y + 1.0 / dy).floor() as usize).min(num_y - 1) * 2;
let z_offset = ((positions[*i].z + 1.0 * dz).floor() as usize)
.min(num_z - 1)
* (num_y * 2);
//figure out where a particle belongs based of it's location along the x-axis
let cell_membership: usize = match positions[*i] {
//scenario one: it's left of 0 or 0
p if (-1.0..=0.0).contains(&p.x) => {
let x_offset =
(((p.x).abs() * half_x as f64).floor() as usize).min(half_x);
x_offset + y_offset + z_offset
}
//scenario two: it's right of zero or one
p if (0.0..=1.).contains(&p.x) => {
let x_offset =
(((p.x).abs() * half_x as f64).floor() as usize).min(half_x);
x_offset + y_offset + z_offset + 1
}
//scenario three: it's no longer relevant
_ => num_cells,
};
let parent = unsafe { &mut *{ cell_ptr }.0.add(*i) };
*parent = cell_membership;
if cell_membership == num_cells {
println!("yeet {particle_count}");
let keep = unsafe { &mut *{ keep_ptr }.0.add(*i) };
*keep = false;
}
(cell_membership, *i)
})
.collect::<Vec<(usize, usize)>>()
}) //send each chunk of indices to sample update function
.for_each(|chunk| sample_sender.send(chunk).unwrap());
}
however, in this function (called shortly after) when printing out keepers, it turns out they are what they were instantiated to: true, so nothing was actually changed by the index function
/// any particles outside of the bounds of the grid need to be discarded as they cannot be indexed.
/// Since this can happen only at the x=-1 or x=1 boundaries, we only need to check the x coordinate
pub(crate) fn filter_out_of_scope(&mut self, num_cells: usize) {
// self.vec.sort_by_key(|p| (p.parent_cell));
// if let Some(cutoff) = particles.iter().rposition(|p| p.parent_cell < num_cells) {
// particles.truncate(cutoff + 1);
// }
let Particles {
velocities,
positions,
types,
parent_cells,
keepers,
} = self;
println!("Starting filter");
println!("keepers: {:?}", keepers);
rayon::join(
|| {
rayon::join(
|| {
let mut keep = keepers.iter();
velocities.retain(|_| *keep.next().unwrap());
},
|| {
let mut keep = keepers.iter();
positions.retain(|_| *keep.next().unwrap());
},
)
},
|| {
rayon::join(
|| {
let mut keep = keepers.iter();
parent_cells.retain(|_| *keep.next().unwrap());
},
|| {
let mut keep = keepers.iter();
types.retain(|_| *keep.next().unwrap());
},
)
},
);
keepers.clear();
keepers.resize(velocities.len(), true);
}
why isn't the value pointed to by the mutable pointer changing?
unrelated but the way I'm using the retain function was copied almost verbatim from the rust docs for retain, so why is it panicking on unwrap none option (the iterator being consumed)?