cellular_raza_core/backend/chili/
update_mechanics.rs

1#[cfg(feature = "tracing")]
2use tracing::instrument;
3
4use super::{
5    AdamsBashforth, CellBox, Communicator, MechanicsAdamsBashforthSolver, SimulationError,
6    SubDomainBox, SubDomainPlainIndex, UpdateMechanics, UpdateNeighborSensing, Voxel,
7    VoxelPlainIndex,
8};
9use cellular_raza_concepts::*;
10
11/// Send about the position of cells between threads.
12///
13/// This type is used during the update steps for cellular mechanics
14/// [update_mechanics_interaction_step_1](super::datastructures::SubDomainBox::update_mechanics_interaction_step_1).
15/// The response to [PosInformation] is the [ForceInformation] type.
16/// Upon requesting the acting force, by providing the information stored in this struct,
17/// the requester obtains the needed information about acting forces.
18/// See also the [cellular_raza_concepts::Interaction] trait.
19pub struct PosInformation<Pos, Vel, Inf> {
20    /// Current position
21    pub pos: Pos,
22    /// Current velocity
23    pub vel: Vel,
24    /// Information shared between cells
25    pub info: Inf,
26    /// Index of cell in stored vector
27    ///
28    /// When returning information, this property is needed in order
29    /// to get the correct cell in the vector of cells and update its properties.
30    pub cell_index_in_vector: usize,
31    /// Voxel index of the sending cell.
32    /// Information should be returned to this voxel.
33    pub index_sender: VoxelPlainIndex,
34    /// Voxel index of the voxel from which information is requested.
35    /// This index is irrelevant after the initial query has been sent.
36    pub index_receiver: VoxelPlainIndex,
37}
38
39/// Return type to the requested [PosInformation].
40///
41/// This type is returned after performing all necessary force calculations in
42/// [update_mechanics_interaction_step_2](super::datastructures::SubDomainBox::update_mechanics_interaction_step_2).
43/// The received information is then used in combination with the already present information
44/// to update the position and velocity of cells in
45/// [update_mechanics_interaction_step_3](super::datastructures::SubDomainBox::update_mechanics_interaction_step_3).
46pub struct ForceInformation<For> {
47    /// Overall force acting on cell.
48    ///
49    /// This force is already combined in the sense that multiple forces may be added together.
50    pub force: For,
51    /// Index of cell in stored vector
52    ///
53    /// This property works in tandem with [Self::index_sender] in order to send
54    /// the calculated information to the correct cell and update its properties.
55    pub cell_index_in_vector: usize,
56    /// The voxel index where information is returned to
57    pub index_sender: VoxelPlainIndex,
58}
59
60/// Send cell and its AuxStorage between threads.
61pub struct SendCell<Cel, Aux>(pub VoxelPlainIndex, pub Cel, pub Aux);
62
63impl<C, A> Voxel<C, A> {
64    #[cfg_attr(feature = "tracing", instrument(skip_all))]
65    pub(crate) fn calculate_force_between_cells_external<
66        Pos,
67        Vel,
68        For,
69        Inf,
70        Float,
71        const N: usize,
72    >(
73        &mut self,
74        ext_pos: &Pos,
75        ext_vel: &Vel,
76        ext_inf: &Inf,
77        neighbor_sensing_func: impl Fn(&mut A, &C, &Pos, &Pos, &Inf) -> Result<(), CalcError>,
78    ) -> Result<Option<For>, CalcError>
79    where
80        C: cellular_raza_concepts::Interaction<Pos, Vel, For, Inf>
81            + cellular_raza_concepts::Position<Pos>
82            + cellular_raza_concepts::Velocity<Vel>,
83        A: UpdateMechanics<Pos, Vel, For, N>,
84        Float: num::Float,
85        For: Xapy<Float> + core::ops::AddAssign,
86    {
87        use core::borrow::BorrowMut;
88        let one_half = Float::one() / (Float::one() + Float::one());
89        let mut force = None;
90        for (cell, aux_storage) in self.cells.iter_mut() {
91            let (f1, f2) = cell.calculate_force_between(
92                &cell.pos(),
93                &cell.velocity(),
94                &ext_pos,
95                &ext_vel,
96                &ext_inf,
97            )?;
98            aux_storage.add_force(f1.xa(one_half));
99            if let Some(f) = force.borrow_mut() {
100                *f = f2.xapy(one_half, &*f);
101            } else {
102                force = Some(f2.xa(one_half));
103            }
104
105            neighbor_sensing_func(aux_storage, &cell, &cell.pos(), &ext_pos, &ext_inf)?;
106        }
107        Ok(force)
108    }
109}
110
111/// Empty function which is used when no neighbor sensing is activated.
112/// This function is called by the [cellular_raza_core_proc_macro::run_simulation] macro.
113pub fn neighbor_sensing_single_empty<Pos, Inf, A, C>(
114    _aux_storage: &mut A,
115    _cell: &C,
116    _p1: &Pos,
117    _p2: &Pos,
118    _i2: &Inf,
119) -> Result<(), CalcError> {
120    Ok(())
121}
122
123/// Uses [cellular_raza_concepts::NeighborSensing] trait to sense neighbors
124/// This function is called by the [cellular_raza_core_proc_macro::run_simulation] macro.
125pub fn neighbor_sensing_single<Pos, Acc, Inf, A, C>(
126    aux_storage: &mut A,
127    cell: &C,
128    p1: &Pos,
129    p2: &Pos,
130    i2: &Inf,
131) -> Result<(), CalcError>
132where
133    A: UpdateNeighborSensing<Acc>,
134    C: NeighborSensing<Pos, Acc, Inf>,
135{
136    let mut acc = aux_storage.get_accumulator();
137    cell.accumulate_information(&p1, &p2, &i2, &mut acc)?;
138    Ok(())
139}
140
141impl<I, S, C, A, Com, Sy> SubDomainBox<I, S, C, A, Com, Sy>
142where
143    S: SubDomain,
144{
145    /// Update cells position and velocity
146    ///
147    /// We assume that cells implement the
148    /// [Mechanics](cellular_raza_concepts::Mechanics) and
149    /// [Interaction](cellular_raza_concepts::Interaction) traits.
150    /// Then, threads will exchange information in the [PosInformation] format
151    /// to calculate the forces acting on the cells.
152    #[cfg_attr(feature = "tracing", instrument(skip_all))]
153    pub fn update_mechanics_interaction_step_1<Pos, Vel, For, Float, Inf, const N: usize>(
154        &mut self,
155        neighbor_sensing_func: impl Fn(&mut A, &C, &Pos, &Pos, &Inf) -> Result<(), CalcError>,
156    ) -> Result<(), SimulationError>
157    where
158        Pos: Clone,
159        Vel: Clone,
160        Inf: Clone,
161        C: cellular_raza_concepts::Position<Pos>,
162        C: cellular_raza_concepts::Velocity<Vel>,
163        C: cellular_raza_concepts::Mechanics<Pos, Vel, For, Float>,
164        C: cellular_raza_concepts::Interaction<Pos, Vel, For, Inf>,
165        A: UpdateMechanics<Pos, Vel, For, N>,
166        For: Xapy<Float> + core::ops::AddAssign,
167        Float: num::Float + core::ops::AddAssign,
168        <S as SubDomain>::VoxelIndex: Ord,
169        S: SubDomainMechanics<Pos, Vel>,
170        Com: Communicator<SubDomainPlainIndex, PosInformation<Pos, Vel, Inf>>,
171    {
172        let one_half: Float = Float::one() / (Float::one() + Float::one());
173        let voxel_indices: Vec<_> = self.voxels.keys().map(|k| *k).collect();
174
175        for voxel_index in voxel_indices {
176            let n_cells = self.voxels[&voxel_index].cells.len();
177            // Calculate Interactions internally
178            for n in 0..n_cells {
179                let vox = self.voxels.get_mut(&voxel_index).unwrap();
180                let p1 = vox.cells[n].0.pos();
181                let v1 = vox.cells[n].0.velocity();
182                let i1 = vox.cells[n].0.get_interaction_information();
183
184                let mut cells_mut = vox.cells.iter_mut();
185                let (c1, aux1) = cells_mut.nth(n).unwrap();
186                for _ in n + 1..n_cells {
187                    let (c2, aux2) = cells_mut.nth(0).unwrap();
188
189                    let p2 = c2.pos();
190                    let v2 = c2.velocity();
191                    let i2 = c2.get_interaction_information();
192
193                    let (force1, force2) = c1.calculate_force_between(&p1, &v1, &p2, &v2, &i2)?;
194                    aux1.add_force(force1.xa(one_half));
195                    aux2.add_force(force2.xa(one_half));
196
197                    let (force2, force1) = c2.calculate_force_between(&p2, &v2, &p1, &v1, &i1)?;
198                    aux1.add_force(force1.xa(one_half));
199                    aux2.add_force(force2.xa(one_half));
200
201                    neighbor_sensing_func(aux1, &c1, &p1, &p2, &i2)?;
202                    neighbor_sensing_func(aux2, &c2, &p2, &p1, &i1)?;
203                }
204
205                let neighbors = vox.neighbors.clone();
206                let mut force = None;
207                for neighbor_index in &neighbors {
208                    match self.voxels.get_mut(&neighbor_index) {
209                        Some(vox) => {
210                            if let Some(f) = vox.calculate_force_between_cells_external(
211                                &p1,
212                                &v1,
213                                &i1,
214                                &neighbor_sensing_func,
215                            )? {
216                                match &mut force {
217                                    Some(f2) => *f2 = f.xapy(Float::one(), &f2),
218                                    f2 @ None => *f2 = Some(f),
219                                }
220                            }
221                        }
222                        None => self.communicator.send(
223                            &self.plain_index_to_subdomain[&neighbor_index],
224                            PosInformation {
225                                index_sender: voxel_index,
226                                index_receiver: neighbor_index.clone(),
227                                pos: p1.clone(),
228                                vel: v1.clone(),
229                                info: i1.clone(),
230                                cell_index_in_vector: n,
231                            },
232                        )?,
233                    }
234                }
235            }
236        }
237        Ok(())
238    }
239}
240
241impl<I, S, C, A, Com, Sy> SubDomainBox<I, S, C, A, Com, Sy>
242where
243    S: SubDomain,
244{
245    /// Calculates the custom [force](SubDomainForce) of
246    /// the domain on the cells.
247    #[cfg_attr(feature = "tracing", instrument(skip_all))]
248    pub fn calculate_custom_domain_force<Pos, Vel, For, Inf, const N: usize>(
249        &mut self,
250    ) -> Result<(), SimulationError>
251    where
252        Pos: Clone,
253        Vel: Clone,
254        C: cellular_raza_concepts::Position<Pos>,
255        C: cellular_raza_concepts::Velocity<Vel>,
256        C: cellular_raza_concepts::InteractionInformation<Inf>,
257        A: UpdateMechanics<Pos, Vel, For, N>,
258        S: SubDomainForce<Pos, Vel, For, Inf>,
259    {
260        for (cell, aux) in self
261            .voxels
262            .iter_mut()
263            .map(|(_, vox)| vox.cells.iter_mut())
264            .flatten()
265        {
266            let f = self.subdomain.calculate_custom_force(
267                &cell.pos(),
268                &cell.velocity(),
269                &cell.get_interaction_information(),
270            )?;
271            aux.add_force(f);
272        }
273        Ok(())
274    }
275
276    /// Update cells position and velocity
277    ///
278    /// We assume that cells implement the
279    /// [Mechanics](cellular_raza_concepts::Mechanics) and
280    /// [Interaction](cellular_raza_concepts::Interaction) traits.
281    /// Then, threads will use the previously exchanged [PosInformation] to calculate forces
282    /// and send back information about the acting force in the [ForceInformation] format.
283    /// In addition, this method also applies the inverse force to local cells.
284    #[cfg_attr(feature = "tracing", instrument(skip_all))]
285    pub fn update_mechanics_interaction_step_2<Pos, Vel, For, Float, Inf, const N: usize>(
286        &mut self,
287        determinism: bool,
288        neighbor_sensing_func: impl Fn(&mut A, &C, &Pos, &Pos, &Inf) -> Result<(), CalcError>,
289    ) -> Result<(), SimulationError>
290    where
291        For: Xapy<Float>,
292        A: UpdateMechanics<Pos, Vel, For, N>,
293        Float: num::Float,
294        Pos: Clone,
295        Vel: Clone,
296        For: Clone + core::ops::AddAssign,
297        C: cellular_raza_concepts::Position<Pos>,
298        C: cellular_raza_concepts::Velocity<Vel>,
299        C: cellular_raza_concepts::Mechanics<Pos, Vel, For, Float>,
300        C: cellular_raza_concepts::Interaction<Pos, Vel, For, Inf>,
301        Com: Communicator<SubDomainPlainIndex, PosInformation<Pos, Vel, Inf>>,
302        Com: Communicator<SubDomainPlainIndex, ForceInformation<For>>,
303    {
304        // Receive PositionInformation and send back ForceInformation
305        let mut received_infos = <Com as Communicator<
306            SubDomainPlainIndex,
307            PosInformation<Pos, Vel, Inf>,
308        >>::receive(&mut self.communicator);
309        if determinism {
310            received_infos.sort_by_key(|pos_info| pos_info.index_sender);
311        }
312        for pos_info in received_infos.iter() {
313            let vox = self.voxels.get_mut(&pos_info.index_receiver).ok_or(
314                cellular_raza_concepts::IndexError(format!(
315                    "EngineError: Voxel with index {:?} of PosInformation can not be\
316                            found in this thread.",
317                    pos_info.index_receiver
318                )),
319            )?;
320            // Calculate force from cells in voxel
321            if let Some(force) = vox.calculate_force_between_cells_external(
322                &pos_info.pos,
323                &pos_info.vel,
324                &pos_info.info,
325                &neighbor_sensing_func,
326            )? {
327                // Send back force information
328                // let thread_index = self.plain_index_to_subdomain[&pos_info.index_sender];
329                self.communicator.send(
330                    &self.plain_index_to_subdomain[&pos_info.index_sender],
331                    ForceInformation {
332                        force,
333                        cell_index_in_vector: pos_info.cell_index_in_vector,
334                        index_sender: pos_info.index_sender,
335                    },
336                )?;
337            }
338        }
339        Ok(())
340    }
341
342    /// Receive all calculated forces and include them for later update steps.
343    #[cfg_attr(feature = "tracing", instrument(skip(self)))]
344    pub fn update_mechanics_interaction_step_3<Pos, Vel, For, const N: usize>(
345        &mut self,
346        determinism: bool,
347    ) -> Result<(), SimulationError>
348    where
349        A: UpdateMechanics<Pos, Vel, For, N>,
350        Com: Communicator<SubDomainPlainIndex, ForceInformation<For>>,
351    {
352        // Update position and velocity of all cells with new information
353        let mut received_infos = <Com as Communicator<
354            SubDomainPlainIndex,
355            ForceInformation<For>,
356        >>::receive(&mut self.communicator);
357        if determinism {
358            received_infos.sort_by_key(|force_info| force_info.index_sender);
359        }
360        for obt_forces in received_infos {
361            let error_1 = format!(
362                "EngineError: Sender with plain index {:?} was ended up in location\
363                where index is not present anymore",
364                obt_forces.index_sender
365            );
366            let vox = self
367                .voxels
368                .get_mut(&obt_forces.index_sender)
369                .ok_or(cellular_raza_concepts::IndexError(error_1))?;
370            let error_2 = format!(
371                "\
372                EngineError: Force Information with sender index {:?} and \
373                cell at vector position {} could not be matched",
374                obt_forces.index_sender, obt_forces.cell_index_in_vector
375            );
376            match vox.cells.get_mut(obt_forces.cell_index_in_vector) {
377                Some((_, aux_storage)) => Ok(aux_storage.add_force(obt_forces.force)),
378                None => Err(cellular_raza_concepts::IndexError(error_2)),
379            }?;
380        }
381        Ok(())
382    }
383
384    /// Applies boundary conditions to cells. For the future, we hope to be using previous and
385    /// current position of cells rather than the cell itself.
386    #[cfg_attr(feature = "tracing", instrument(skip_all))]
387    pub fn apply_boundary<Pos, Vel>(&mut self) -> Result<(), BoundaryError>
388    where
389        C: cellular_raza_concepts::Position<Pos>,
390        C: cellular_raza_concepts::Velocity<Vel>,
391        S: SubDomainMechanics<Pos, Vel>,
392    {
393        for (cell, _) in self
394            .voxels
395            .iter_mut()
396            .map(|(_, voxel)| voxel.cells.iter_mut())
397            .flatten()
398        {
399            let mut pos = cell.pos();
400            let mut vel = cell.velocity();
401            self.subdomain.apply_boundary(&mut pos, &mut vel)?;
402            cell.set_pos(&pos);
403            cell.set_velocity(&vel);
404        }
405        Ok(())
406    }
407
408    /// Sort new cells into respective voxels
409    ///
410    /// This step determines if a cell is still in its correct location
411    /// after its position has changed. This can be due to the
412    /// [SubDomainBox::update_mechanics_interaction_step_3] method or due to other effects such as
413    /// cell-division by the [cellular_raza_concepts::Cycle] trait.
414    ///
415    /// If the cell is not in the correct voxel, we either directly insert this cell into the voxel
416    /// or send it to the other [SubDomainBox] to take care of this.
417    #[cfg_attr(feature = "tracing", instrument(skip_all))]
418    pub fn sort_cells_in_voxels_step_1(&mut self) -> Result<(), SimulationError>
419    where
420        Com: Communicator<SubDomainPlainIndex, SendCell<CellBox<C>, A>>,
421        S: SortCells<C, VoxelIndex = <S as SubDomain>::VoxelIndex>,
422        <S as SubDomain>::VoxelIndex: Eq + core::hash::Hash + Ord,
423    {
424        // Store all cells which need to find a new home in this variable
425        let mut find_new_home_cells = Vec::<_>::new();
426
427        for (voxel_index, vox) in self.voxels.iter_mut() {
428            // Drain every cell which is currently not in the correct voxel
429            // TODO use drain_filter when stabilized
430            let mut errors = Vec::new();
431            let (new_voxel_cells, old_voxel_cells): (Vec<_>, Vec<_>) =
432                vox.cells.drain(..).partition(|(c, _)| {
433                    let cell_index = self.subdomain.get_voxel_index_of(&c);
434                    match cell_index {
435                        Ok(index) => {
436                            let plain_index = self.voxel_index_to_plain_index.get(&index);
437                            match plain_index {
438                                Some(i) => i != voxel_index,
439                                None => {
440                                    errors.push(SimulationError::from(IndexError(format!(
441                                        "Could not find matching plain index for given voxel index"
442                                    ))));
443                                    false
444                                }
445                            }
446                            // &plain_index != voxel_index
447                        }
448                        Err(e) => {
449                            errors.push(e.into());
450                            false
451                        }
452                    }
453                });
454            find_new_home_cells.extend(new_voxel_cells.into_iter().map(|c| (*voxel_index, c)));
455            vox.cells = old_voxel_cells;
456        }
457
458        // Send cells to other multivoxelcontainer or keep them here
459        for (voxel_index, (cell, aux_storage)) in find_new_home_cells {
460            let ind = self.subdomain.get_voxel_index_of(&cell)?;
461            let cell_index = self.voxel_index_to_plain_index[&ind];
462            match self.voxels.get_mut(&cell_index) {
463                // If new voxel is in current multivoxelcontainer then save them there
464                Some(vox) => {
465                    vox.cells.push((cell, aux_storage));
466                    Ok(())
467                }
468                // Otherwise send them to the correct other multivoxelcontainer
469                None => <Com as Communicator<SubDomainPlainIndex, SendCell<CellBox<C>, A>>>::send(
470                    &mut self.communicator,
471                    &self.plain_index_to_subdomain[&cell_index],
472                    SendCell(voxel_index, cell, aux_storage),
473                ),
474            }?;
475        }
476        Ok(())
477    }
478
479    /// Sort new cells into respective voxels
480    ///
481    /// After having sent cells to the new [SubDomainBox] in the
482    /// [SubDomainBox::sort_cells_in_voxels_step_1] method, we receive these new cells and insert
483    /// them into their respective voxels.
484    #[cfg_attr(feature = "tracing", instrument(skip_all))]
485    pub fn sort_cells_in_voxels_step_2(&mut self, determinism: bool) -> Result<(), SimulationError>
486    where
487        Com: Communicator<SubDomainPlainIndex, SendCell<CellBox<C>, A>>,
488        <S as SubDomain>::VoxelIndex: Eq + core::hash::Hash + Ord,
489        S: SortCells<C, VoxelIndex = <S as SubDomain>::VoxelIndex>,
490    {
491        // Now receive new cells and insert them
492        let mut received_cells = <Com as Communicator<
493            SubDomainPlainIndex,
494            SendCell<CellBox<C>, A>,
495        >>::receive(&mut self.communicator);
496        if determinism {
497            received_cells.sort_by_key(|send_cell| send_cell.0);
498        }
499        for sent_cell in received_cells {
500            let SendCell(_, cell, aux_storage) = sent_cell;
501            let index =
502                self.voxel_index_to_plain_index[&self.subdomain.get_voxel_index_of(&cell)?];
503
504            match self.voxels.get_mut(&index) {
505                Some(vox) => Ok(vox.cells.push((cell, aux_storage))),
506                None => Err(cellular_raza_concepts::IndexError(format!(
507                    "Cell with index {:?} was sent to subdomain which does not hold this index",
508                    index
509                ))),
510            }?;
511        }
512        Ok(())
513    }
514}
515
516/// We assume that cells implement the
517/// [Mechanics](cellular_raza_concepts::Mechanics) and
518/// [Interaction](cellular_raza_concepts::Interaction) traits.
519/// In this last step, all [ForceInformation] are gathered and used to update the
520/// cells positions and velocities.
521///
522/// For the future, we hope to provide an abstracted method to use any of our implemented
523/// solvers.
524/// The solver currently limits the number of saved previous increments in the
525/// [UpdateMechanics] trait.
526///
527/// Currently, we employ the [mechanics_adams_bashforth_3](super::mechanics_adams_bashforth_3)
528/// solver.
529#[allow(private_bounds)]
530pub fn local_mechanics_update<
531    C,
532    A,
533    Pos,
534    Vel,
535    For,
536    #[cfg(feature = "tracing")] Float: core::fmt::Debug,
537    #[cfg(not(feature = "tracing"))] Float,
538    const N: usize,
539>(
540    cell: &mut C,
541    aux_storage: &mut A,
542    dt: Float,
543    rng: &mut rand_chacha::ChaCha8Rng,
544) -> Result<(), SimulationError>
545where
546    A: UpdateMechanics<Pos, Vel, For, N>,
547    C: cellular_raza_concepts::Mechanics<Pos, Vel, For, Float> + Clone,
548    C: cellular_raza_concepts::Position<Pos>,
549    C: cellular_raza_concepts::Velocity<Vel>,
550    Float: num::Float + Copy + num::FromPrimitive,
551    Pos: Xapy<Float> + Clone,
552    Vel: Xapy<Float> + Clone,
553    Vel: Clone,
554    MechanicsAdamsBashforthSolver<N>: AdamsBashforth<N>,
555{
556    use super::solvers::{AdamsBashforth, MechanicsAdamsBashforthSolver};
557    <MechanicsAdamsBashforthSolver<N> as AdamsBashforth<N>>::update(cell, aux_storage, dt, rng)?;
558    Ok(())
559}
560
561/// Perform the [Interaction::react_to_neighbors] function and clear current neighbors.
562pub fn local_interaction_react_to_neighbors<C, A, Pos, Acc, Inf, Float>(
563    cell: &mut C,
564    aux_storage: &mut A,
565    _dt: Float,
566    _rng: &mut rand_chacha::ChaCha8Rng,
567) -> Result<(), cellular_raza_concepts::CalcError>
568where
569    C: cellular_raza_concepts::NeighborSensing<Pos, Acc, Inf>,
570    C: cellular_raza_concepts::Position<Pos>,
571    A: UpdateNeighborSensing<Acc>,
572{
573    let mut acc = aux_storage.get_accumulator();
574    cell.react_to_neighbors(&acc)?;
575    <C as NeighborSensing<Pos, Acc, Inf>>::clear_accumulator(&mut acc);
576    Ok(())
577}
cellular_raza_core/backend/chili/update_mechanics.rs

cellular_raza_core/backend/chili/
update_mechanics.rs
