cellular_raza_core/backend/cpu_os_threads/
domain_decomposition.rs

1use super::{Agent, ForceBound, InteractionInf, PositionBound, VelocityBound};
2use cellular_raza_concepts::domain_old::*;
3use cellular_raza_concepts::reactions_old::*;
4use cellular_raza_concepts::*;
5
6use super::errors::*;
7use crate::storage::StorageManager;
8
9use std::collections::{BTreeMap, HashMap};
10
11use std::ops::AddAssign;
12use std::ops::{Add, Mul};
13
14use crossbeam_channel::{Receiver, SendError, Sender};
15use hurdles::Barrier;
16
17use num::Zero;
18use serde::{Deserialize, Serialize};
19
20use rand_chacha::ChaCha8Rng;
21
22#[cfg(feature = "tracing")]
23use tracing::instrument;
24
25/// Wrapper around the user-defined simulation [Domain]
26#[derive(Clone, Serialize, Deserialize)]
27pub struct DomainBox<Dom> {
28    /// Raw simulation [Domain]
29    pub domain_raw: Dom,
30}
31
32impl<Dom> From<Dom> for DomainBox<Dom> {
33    fn from(domain: Dom) -> DomainBox<Dom> {
34        DomainBox { domain_raw: domain }
35    }
36}
37
38impl<Cel, Ind, Vox, Dom> cellular_raza_concepts::domain_old::Domain<CellBox<Cel>, Ind, Vox>
39    for DomainBox<Dom>
40where
41    Dom: cellular_raza_concepts::domain_old::Domain<Cel, Ind, Vox>,
42    Vox: Send + Sync,
43    Cel: Send + Sync,
44{
45    fn apply_boundary(&self, cbox: &mut CellBox<Cel>) -> Result<(), BoundaryError> {
46        self.domain_raw.apply_boundary(&mut cbox.cell)
47    }
48
49    fn get_neighbor_voxel_indices(&self, index: &Ind) -> Vec<Ind> {
50        self.domain_raw.get_neighbor_voxel_indices(index)
51    }
52
53    fn get_voxel_index(&self, cbox: &CellBox<Cel>) -> Ind {
54        self.domain_raw.get_voxel_index(&cbox.cell)
55    }
56
57    fn get_all_indices(&self) -> Vec<Ind> {
58        self.domain_raw.get_all_indices()
59    }
60
61    fn generate_contiguous_multi_voxel_regions(
62        &self,
63        n_regions: usize,
64    ) -> Result<Vec<Vec<(Ind, Vox)>>, CalcError> {
65        self.domain_raw
66            .generate_contiguous_multi_voxel_regions(n_regions)
67    }
68}
69
70/// This is a purely implementational detail and should not be of any concern to the end user.
71pub type PlainIndex = usize;
72
73pub(crate) struct IndexBoundaryInformation<Ind> {
74    pub index_original_sender: PlainIndex,
75    pub index_original_sender_raw: Ind,
76    pub index_original_receiver: PlainIndex,
77}
78
79pub(crate) struct ConcentrationBoundaryInformation<ConcVecExtracellular, Ind> {
80    pub index_original_sender: PlainIndex,
81    pub concentration_boundary: BoundaryCondition<ConcVecExtracellular>,
82    pub index_original_receiver_raw: Ind,
83}
84
85pub(crate) struct PosInformation<Pos, Vel, Inf> {
86    pub pos: Pos,
87    pub vel: Vel,
88    pub info: Inf,
89    pub count: usize,
90    pub index_sender: PlainIndex,
91    pub index_receiver: PlainIndex,
92}
93
94pub(crate) struct ForceInformation<For> {
95    pub force: For,
96    pub count: usize,
97    pub index_sender: PlainIndex,
98}
99
100/// Wrapper for a [Voxel] struct that includes more information and cells.
101#[derive(Serialize, Deserialize, Clone)]
102pub struct VoxelBox<
103    Ind,
104    Pos,
105    Vel,
106    For,
107    Vox,
108    Cel,
109    ConcVecExtracellular,
110    ConcBoundaryExtracellular,
111    ConcVecIntracellular,
112> {
113    pub(crate) plain_index: PlainIndex,
114    pub(crate) index: Ind,
115    pub(crate) voxel: Vox,
116    pub(crate) neighbors: Vec<PlainIndex>,
117    pub(crate) cells: Vec<(
118        CellBox<Cel>,
119        AuxiliaryCellPropertyStorage<Pos, Vel, For, ConcVecIntracellular>,
120    )>,
121    pub(crate) new_cells: Vec<(Cel, Option<CellIdentifier>)>,
122    pub(crate) id_counter: u64,
123    pub(crate) rng: ChaCha8Rng,
124    pub(crate) extracellular_concentration_increments: Vec<(Pos, ConcVecExtracellular)>,
125    pub(crate) concentration_boundaries: Vec<(Ind, BoundaryCondition<ConcBoundaryExtracellular>)>,
126}
127
128impl<
129    Ind,
130    Vox,
131    Cel,
132    Pos,
133    Vel,
134    For,
135    ConcVecExtracellular,
136    ConcBoundaryExtracellular,
137    ConcVecIntracellular,
138> Id
139    for VoxelBox<
140        Ind,
141        Pos,
142        Vel,
143        For,
144        Vox,
145        Cel,
146        ConcVecExtracellular,
147        ConcBoundaryExtracellular,
148        ConcVecIntracellular,
149    >
150{
151    type Identifier = PlainIndex;
152
153    fn get_id(&self) -> PlainIndex {
154        self.plain_index
155    }
156
157    fn ref_id(&self) -> &Self::Identifier {
158        &self.plain_index
159    }
160}
161
162/// Contains auxiliary information to update cellular properties
163#[derive(Serialize, Deserialize, Clone)]
164pub struct AuxiliaryCellPropertyStorage<Pos, Vel, For, ConcVecIntracellular> {
165    force: For,
166    intracellular_concentration_increment: ConcVecIntracellular,
167    pub(crate) cycle_events: Vec<CycleEvent>,
168    neighbor_count: usize,
169
170    inc_pos_back_1: Option<Pos>,
171    inc_pos_back_2: Option<Pos>,
172    inc_vel_back_1: Option<Vel>,
173    inc_vel_back_2: Option<Vel>,
174}
175
176impl<Pos, Vel, For, ConcVecIntracellular> Default
177    for AuxiliaryCellPropertyStorage<Pos, Vel, For, ConcVecIntracellular>
178where
179    For: Zero,
180    ConcVecIntracellular: Zero,
181{
182    fn default() -> AuxiliaryCellPropertyStorage<Pos, Vel, For, ConcVecIntracellular> {
183        AuxiliaryCellPropertyStorage {
184            force: For::zero(),
185            intracellular_concentration_increment: ConcVecIntracellular::zero(),
186            cycle_events: Vec::new(),
187            neighbor_count: 0,
188
189            inc_pos_back_1: None,
190            inc_pos_back_2: None,
191            inc_vel_back_1: None,
192            inc_vel_back_2: None,
193        }
194    }
195}
196
197impl<
198    Ind,
199    Vox,
200    Cel,
201    Pos,
202    Vel,
203    For,
204    ConcVecExtracellular,
205    ConcBoundaryExtracellular,
206    ConcVecIntracellular,
207>
208    VoxelBox<
209        Ind,
210        Pos,
211        Vel,
212        For,
213        Vox,
214        Cel,
215        ConcVecExtracellular,
216        ConcBoundaryExtracellular,
217        ConcVecIntracellular,
218    >
219where
220    Ind: Clone,
221    For: num::Zero,
222    ConcVecIntracellular: Zero,
223{
224    pub(crate) fn new(
225        plain_index: PlainIndex,
226        index: Ind,
227        voxel: Vox,
228        neighbors: Vec<PlainIndex>,
229        cells: Vec<CellBox<Cel>>,
230        rng_seed: u64,
231    ) -> VoxelBox<
232        Ind,
233        Pos,
234        Vel,
235        For,
236        Vox,
237        Cel,
238        ConcVecExtracellular,
239        ConcBoundaryExtracellular,
240        ConcVecIntracellular,
241    > {
242        use rand::SeedableRng;
243        let n_cells = cells.len() as u64;
244        VoxelBox {
245            plain_index,
246            index,
247            voxel,
248            neighbors,
249            cells: cells
250                .into_iter()
251                .map(|cell| (cell, AuxiliaryCellPropertyStorage::default()))
252                .collect(),
253            new_cells: Vec::new(),
254            id_counter: n_cells,
255            rng: ChaCha8Rng::seed_from_u64(rng_seed),
256            extracellular_concentration_increments: Vec::new(),
257            concentration_boundaries: Vec::new(),
258        }
259    }
260}
261
262impl<
263    Ind,
264    Vox,
265    Cel,
266    Pos,
267    Vel,
268    For,
269    ConcVecExtracellular,
270    ConcBoundaryExtracellular,
271    ConcVecIntracellular,
272>
273    VoxelBox<
274        Ind,
275        Pos,
276        Vel,
277        For,
278        Vox,
279        Cel,
280        ConcVecExtracellular,
281        ConcBoundaryExtracellular,
282        ConcVecIntracellular,
283    >
284where
285    Ind: Clone,
286    Pos: Serialize + for<'a> Deserialize<'a>,
287    Vel: Serialize + for<'a> Deserialize<'a>,
288    Cel: Serialize + for<'a> Deserialize<'a>,
289{
290    #[cfg_attr(feature = "tracing", instrument(skip_all))]
291    fn calculate_custom_force_on_cells(&mut self) -> Result<(), CalcError>
292    where
293        Vox: Voxel<Ind, Pos, Vel, For>,
294        Ind: Index,
295        Pos: PositionBound,
296        Vel: VelocityBound,
297        For: ForceBound,
298        Cel: Mechanics<Pos, Vel, For>,
299        Cel: Position<Pos>,
300        Cel: Velocity<Vel>,
301    {
302        for (cell, aux_storage) in self.cells.iter_mut() {
303            match self
304                .voxel
305                .custom_force_on_cell(&cell.pos(), &cell.velocity())
306            {
307                Some(Ok(force)) => Ok(aux_storage.force += force),
308                Some(Err(e)) => Err(e),
309                None => Ok(()),
310            }?;
311        }
312        Ok(())
313    }
314
315    /// This code relies on the following snippet to work
316    /// ```
317    /// let n_elements: usize = 6;
318    ///
319    /// let mut v: Vec<_> = (0..n_elements).collect();
320    ///
321    /// for n in 0..v.len() {
322    ///     for m in n+1..v.len() {
323    ///         let mut elements_mut = v.iter_mut();
324    ///
325    ///         let vn = elements_mut.nth(n).unwrap();
326    ///         let vm = elements_mut.nth(m-n-1).unwrap();
327    ///
328    ///         println!("Cells {} and {} interact", vn, vm);
329    ///         assert_ne!(vn, vm);
330    ///         assert_eq!(m>n,true);
331    ///     }
332    /// }
333    /// ```
334    #[cfg_attr(feature = "tracing", instrument(skip_all))]
335    fn calculate_force_between_cells_internally<Inf>(&mut self) -> Result<(), CalcError>
336    where
337        Vox: Voxel<Ind, Pos, Vel, For>,
338        Ind: Index,
339        Pos: PositionBound,
340        Vel: VelocityBound,
341        For: ForceBound,
342        Cel: Interaction<Pos, Vel, For, Inf> + Mechanics<Pos, Vel, For> + Clone,
343        Cel: Position<Pos>,
344        Cel: Velocity<Vel>,
345    {
346        use cellular_raza_concepts::InteractionInformation;
347
348        for n in 0..self.cells.len() {
349            for m in n + 1..self.cells.len() {
350                let mut cells_mut = self.cells.iter_mut();
351                let (c1, aux1) = cells_mut.nth(n).unwrap();
352                let (c2, aux2) = cells_mut.nth(m - n - 1).unwrap();
353
354                let p1 = c1.pos();
355                let v1 = c1.velocity();
356
357                let p2 = c2.pos();
358                let v2 = c2.velocity();
359                let i2 = c2.get_interaction_information();
360
361                let (force1, force2) = c1.calculate_force_between(&p1, &v1, &p2, &v2, &i2)?;
362                aux1.force += force1 * 0.5;
363                aux2.force += force2 * 0.5;
364            }
365        }
366        Ok(())
367    }
368
369    #[cfg_attr(feature = "tracing", instrument(skip_all))]
370    fn calculate_force_between_cells_external<Inf>(
371        &mut self,
372        ext_pos: &Pos,
373        ext_vel: &Vel,
374        ext_inf: &Inf,
375    ) -> Result<For, CalcError>
376    where
377        Vox: Voxel<Ind, Pos, Vel, For>,
378        Ind: Index,
379        Pos: PositionBound,
380        Vel: VelocityBound,
381        For: ForceBound,
382        Cel: Interaction<Pos, Vel, For, Inf> + Mechanics<Pos, Vel, For>,
383        Cel: Position<Pos>,
384        Cel: Velocity<Vel>,
385    {
386        let mut force = For::zero();
387        for (cell, aux_storage) in self.cells.iter_mut() {
388            let (f1, f2) = cell.calculate_force_between(
389                &cell.pos(),
390                &cell.velocity(),
391                &ext_pos,
392                &ext_vel,
393                &ext_inf,
394            )?;
395            aux_storage.force += f1;
396            force += f2;
397        }
398        Ok(force)
399    }
400
401    #[cfg_attr(feature = "tracing", instrument(skip_all))]
402    fn update_cell_cycle(&mut self, dt: &f64) -> Result<(), SimulationError>
403    where
404        Cel: Cycle<Cel>,
405        AuxiliaryCellPropertyStorage<Pos, Vel, For, ConcVecIntracellular>: Default,
406    {
407        // Update the cell individual cells
408        self.cells
409            .iter_mut()
410            .map(|(cbox, aux_storage)| {
411                // Check for cycle events and do update if necessary
412                let mut remaining_events = Vec::new();
413                for event in aux_storage.cycle_events.drain(..) {
414                    match event {
415                        CycleEvent::Division => {
416                            let new_cell = Cel::divide(&mut self.rng, &mut cbox.cell)?;
417                            self.new_cells.push((new_cell, Some(cbox.get_id())));
418                        }
419                        CycleEvent::Remove => remaining_events.push(event),
420                        CycleEvent::PhasedDeath => {
421                            remaining_events.push(event);
422                        }
423                    };
424                }
425                aux_storage.cycle_events = remaining_events;
426                // Update the cell cycle
427                if aux_storage.cycle_events.contains(&CycleEvent::PhasedDeath) {
428                    match Cel::update_conditional_phased_death(&mut self.rng, dt, &mut cbox.cell)? {
429                        true => aux_storage.cycle_events.push(CycleEvent::Remove),
430                        false => (),
431                    }
432                } else {
433                    match Cel::update_cycle(&mut self.rng, dt, &mut cbox.cell) {
434                        Some(event) => aux_storage.cycle_events.push(event),
435                        None => (),
436                    }
437                }
438                Ok(())
439            })
440            .collect::<Result<(), SimulationError>>()?;
441
442        // Remove cells which are flagged for death
443        self.cells
444            .retain(|(_, aux_storage)| !aux_storage.cycle_events.contains(&CycleEvent::Remove));
445
446        // Include new cells
447        self.cells
448            .extend(self.new_cells.drain(..).map(|(cell, parent_id)| {
449                self.id_counter += 1;
450                (
451                    CellBox::new(
452                        VoxelPlainIndex(self.plain_index),
453                        self.id_counter,
454                        cell,
455                        parent_id,
456                    ),
457                    AuxiliaryCellPropertyStorage::default(),
458                )
459            }));
460        Ok(())
461    }
462}
463
464// This object has multiple voxels and runs on a single thread.
465// It can communicate with other containers via channels.
466/// Subdomain which consists of multiple [Voxels](Voxel)
467pub struct MultiVoxelContainer<
468    Ind,
469    Pos,
470    Vel,
471    For,
472    Inf,
473    Vox,
474    Dom,
475    Cel,
476    ConcVecExtracellular = (),
477    ConcBoundaryExtracellular = (),
478    ConcVecIntracellular = (),
479> where
480    Pos: Serialize + for<'a> Deserialize<'a>,
481    For: Serialize + for<'a> Deserialize<'a>,
482    Vel: Serialize + for<'a> Deserialize<'a>,
483    Cel: Serialize + for<'a> Deserialize<'a>,
484    Dom: Serialize + for<'a> Deserialize<'a>,
485    ConcVecExtracellular: Serialize + for<'a> Deserialize<'a> + 'static,
486    ConcBoundaryExtracellular: Serialize + for<'a> Deserialize<'a>,
487    ConcVecIntracellular: Serialize + for<'a> Deserialize<'a>,
488{
489    pub(crate) voxels: BTreeMap<
490        PlainIndex,
491        VoxelBox<
492            Ind,
493            Pos,
494            Vel,
495            For,
496            Vox,
497            Cel,
498            ConcVecExtracellular,
499            ConcBoundaryExtracellular,
500            ConcVecIntracellular,
501        >,
502    >,
503
504    // TODO
505    // Maybe we need to implement this somewhere else since
506    // it is currently not simple to change this variable on the fly.
507    // However, maybe we should be thinking about specifying an interface to use this function
508    // Something like:
509    // fn update_domain(&mut self, domain: Domain) -> Result<(), BoundaryError>
510    // And then automatically have the ability to change cell positions if the domain shrinks/grows
511    // for example but then we might also want to change the number of voxels and redistribute cells
512    // accordingly
513    // This needs much more though!
514    pub(crate) domain: DomainBox<Dom>,
515    pub(crate) index_to_plain_index: BTreeMap<Ind, PlainIndex>,
516    pub(crate) plain_index_to_thread: BTreeMap<PlainIndex, usize>,
517    pub(crate) index_to_thread: BTreeMap<Ind, usize>,
518
519    pub(crate) senders_cell: HashMap<
520        usize,
521        Sender<(
522            CellBox<Cel>,
523            AuxiliaryCellPropertyStorage<Pos, Vel, For, ConcVecIntracellular>,
524        )>,
525    >,
526    pub(crate) senders_pos: HashMap<usize, Sender<PosInformation<Pos, Vel, Inf>>>,
527    pub(crate) senders_force: HashMap<usize, Sender<ForceInformation<For>>>,
528
529    pub(crate) senders_boundary_index: HashMap<usize, Sender<IndexBoundaryInformation<Ind>>>,
530    pub(crate) senders_boundary_concentrations:
531        HashMap<usize, Sender<ConcentrationBoundaryInformation<ConcBoundaryExtracellular, Ind>>>,
532
533    // Same for receiving
534    pub(crate) receiver_cell: Receiver<(
535        CellBox<Cel>,
536        AuxiliaryCellPropertyStorage<Pos, Vel, For, ConcVecIntracellular>,
537    )>,
538    pub(crate) receiver_pos: Receiver<PosInformation<Pos, Vel, Inf>>,
539    pub(crate) receiver_force: Receiver<ForceInformation<For>>,
540
541    pub(crate) receiver_index: Receiver<IndexBoundaryInformation<Ind>>,
542    pub(crate) receiver_concentrations:
543        Receiver<ConcentrationBoundaryInformation<ConcBoundaryExtracellular, Ind>>,
544
545    // Global barrier to synchronize threads and make sure every information is sent before further
546    // processing
547    pub(crate) barrier: Barrier,
548
549    pub(crate) storage_cells: StorageManager<CellIdentifier, CellBox<Cel>>,
550    pub(crate) storage_voxels: StorageManager<
551        PlainIndex,
552        VoxelBox<
553            Ind,
554            Pos,
555            Vel,
556            For,
557            Vox,
558            Cel,
559            ConcVecExtracellular,
560            ConcBoundaryExtracellular,
561            ConcVecIntracellular,
562        >,
563    >,
564
565    pub(crate) mvc_id: u32,
566}
567
568impl<
569    Ind,
570    Pos,
571    Vel,
572    For,
573    Inf,
574    ConcVecExtracellular,
575    ConcBoundaryExtracellular,
576    ConcVecIntracellular,
577    Vox,
578    Dom,
579    Cel,
580>
581    MultiVoxelContainer<
582        Ind,
583        Pos,
584        Vel,
585        For,
586        Inf,
587        Vox,
588        Dom,
589        Cel,
590        ConcVecExtracellular,
591        ConcBoundaryExtracellular,
592        ConcVecIntracellular,
593    >
594where
595    // TODO abstract away these trait bounds to more abstract traits
596    // these traits should be defined when specifying the individual cell components
597    // (eg. mechanics, interaction, etc...)
598    Ind: Index + Serialize + for<'a> Deserialize<'a>,
599    Vox: Voxel<Ind, Pos, Vel, For>,
600    Dom: cellular_raza_concepts::domain_old::Domain<Cel, Ind, Vox>,
601    Pos: Serialize + for<'a> Deserialize<'a>,
602    Vel: Serialize + for<'a> Deserialize<'a>,
603    For: ForceBound + Serialize + for<'a> Deserialize<'a>,
604    Inf: Clone,
605    Cel: Serialize + for<'a> Deserialize<'a> + Send + Sync,
606    ConcVecExtracellular: Serialize + for<'a> Deserialize<'a>,
607    ConcBoundaryExtracellular: Serialize + for<'a> Deserialize<'a>,
608    ConcVecIntracellular: Serialize + for<'a> Deserialize<'a> + Zero,
609{
610    #[cfg_attr(feature = "tracing", instrument(skip_all))]
611    fn update_local_functions<ConcGradientExtracellular, ConcTotalExtracellular>(
612        &mut self,
613        dt: &f64,
614    ) -> Result<(), SimulationError>
615    where
616        Pos: PositionBound,
617        Vel: VelocityBound,
618        For: ForceBound,
619        Inf: InteractionInf,
620        Cel: Agent<Pos, Vel, For, Inf>
621            + InteractionExtracellularGradient<Cel, ConcGradientExtracellular>,
622        Vox: ExtracellularMechanics<
623                Ind,
624                Pos,
625                ConcVecExtracellular,
626                ConcGradientExtracellular,
627                ConcTotalExtracellular,
628                ConcBoundaryExtracellular,
629            >,
630    {
631        self.voxels
632            .iter_mut()
633            .map(|(_, vox)| {
634                use cellular_raza_concepts::domain_old::Domain;
635                // Update all local functions inside the voxel
636                vox.update_cell_cycle(dt)?;
637
638                // TODO every voxel should apply its own boundary conditions
639                // This is now a global rule but we do not want this
640                // This should not be dependent on the domain
641                // Apply boundary conditions to the cells in the respective voxels
642                vox.cells
643                    .iter_mut()
644                    .map(|(cell, _)| self.domain.apply_boundary(cell))
645                    .collect::<Result<(), BoundaryError>>()?;
646
647                #[cfg(feature = "gradients")]
648                vox.cells
649                    .iter_mut()
650                    .map(|(cell, _)| {
651                        let gradient = vox
652                            .voxel
653                            .get_extracellular_gradient_at_point(&cell.cell.pos())?;
654                        Cel::sense_gradient(&mut cell.cell, &gradient)?;
655                        Ok(())
656                    })
657                    .collect::<Result<(), SimulationError>>()?;
658                Ok(())
659            })
660            .collect::<Result<(), SimulationError>>()
661    }
662
663    // TODO make sure that if no extracellular mechanics are in action updating is correct and the
664    // trait may be adjusted
665    #[cfg_attr(feature = "tracing", instrument(skip_all))]
666    fn update_cellular_reactions<ConcGradientExtracellular, ConcTotalExtracellular>(
667        &mut self,
668        dt: &f64,
669    ) -> Result<(), SimulationError>
670    where
671        ConcVecIntracellular: std::ops::AddAssign + Mul<f64, Output = ConcVecIntracellular>,
672        Cel: Mechanics<Pos, Vel, For>
673            + CellularReactions<ConcVecIntracellular, ConcVecExtracellular>
674            + Volume,
675        Cel: Position<Pos>,
676        Cel: Velocity<Vel>,
677        Vox: ExtracellularMechanics<
678                Ind,
679                Pos,
680                ConcVecExtracellular,
681                ConcGradientExtracellular,
682                ConcTotalExtracellular,
683                ConcBoundaryExtracellular,
684            > + Volume,
685        ConcVecExtracellular: core::ops::Mul<f64, Output = ConcVecExtracellular>,
686    {
687        self.voxels
688            .iter_mut()
689            .map(|(_, voxelbox)| {
690                voxelbox.cells.iter_mut().map(
691                    |(cellbox, _aux_storage)| -> Result<(), SimulationError> {
692                        let internal_concentration_vector = cellbox.cell.get_intracellular();
693
694                        let external_concentration_vector = voxelbox
695                            .voxel
696                            .get_extracellular_at_point(&cellbox.cell.pos())?;
697                        let (increment_intracellular, increment_extracellular) = cellbox
698                            .cell
699                            .calculate_intra_and_extracellular_reaction_increment(
700                                &internal_concentration_vector,
701                                &external_concentration_vector,
702                            )?;
703
704                        let cell_volume = cellbox.cell.get_volume();
705                        let voxel_volume = voxelbox.voxel.get_volume();
706                        let cell_to_voxel_volume = cell_volume / voxel_volume;
707
708                        voxelbox.extracellular_concentration_increments.push((
709                            cellbox.cell.pos(),
710                            increment_extracellular * cell_to_voxel_volume,
711                        ));
712                        // TODO these reactions are currently on the same timescale as the
713                        // fluid-dynamics but we should consider how this may change if we have
714                        // different time-scales here
715                        // Also the solver is currently simply an euler stepper.
716                        // This should be altered to have something like an (adaptive) Runge Kutta
717                        // or Dopri (or better)
718                        cellbox.cell.set_intracellular(
719                            internal_concentration_vector + increment_intracellular * *dt,
720                        );
721                        Ok(())
722                    },
723                )
724            })
725            .flatten()
726            .collect::<Result<(), SimulationError>>()
727    }
728
729    #[cfg_attr(feature = "tracing", instrument(skip_all))]
730    fn sort_cell_in_voxel(
731        &mut self,
732        cell: CellBox<Cel>,
733        aux_storage: AuxiliaryCellPropertyStorage<Pos, Vel, For, ConcVecIntracellular>,
734    ) -> Result<(), SimulationError> {
735        use cellular_raza_concepts::domain_old::Domain;
736        let index = self.index_to_plain_index[&self.domain.get_voxel_index(&cell)];
737
738        match self.voxels.get_mut(&index) {
739            Some(vox) => vox.cells.push((cell, aux_storage)),
740            None => {
741                let thread_index = self.plain_index_to_thread[&index];
742                match self.senders_cell.get(&thread_index) {
743                    Some(sender) => sender.send((cell, aux_storage)),
744                    None => Err(SendError((cell, aux_storage))),
745                }?;
746            }
747        }
748        Ok(())
749    }
750
751    #[cfg_attr(feature = "tracing", instrument(skip_all))]
752    fn update_fluid_mechanics_step_1<ConcGradientExtracellular, ConcTotalExtracellular>(
753        &mut self,
754    ) -> Result<(), SimulationError>
755    where
756        Vox: ExtracellularMechanics<
757                Ind,
758                Pos,
759                ConcVecExtracellular,
760                ConcGradientExtracellular,
761                ConcTotalExtracellular,
762                ConcBoundaryExtracellular,
763            >,
764    {
765        let indices_iterator = self
766            .voxels
767            .iter()
768            .map(|(ind, vox)| (ind.clone(), vox.index.clone(), vox.neighbors.clone()))
769            .collect::<Vec<_>>();
770        for (voxel_plain_index, voxel_index, neighbor_indices) in indices_iterator.into_iter() {
771            for neighbor_index in neighbor_indices {
772                match self.voxels.get(&neighbor_index) {
773                    Some(neighbor_voxel) => {
774                        let neighbor_voxel_index_raw = neighbor_voxel.index.clone();
775                        let bc = neighbor_voxel
776                            .voxel
777                            .boundary_condition_to_neighbor_voxel(&voxel_index)?;
778                        let vox = self.voxels.get_mut(&voxel_plain_index).unwrap();
779                        vox.concentration_boundaries
780                            .push((neighbor_voxel_index_raw, bc));
781                        Ok(())
782                    }
783                    None => self.senders_boundary_index
784                        [&self.plain_index_to_thread[&neighbor_index]]
785                        .send(IndexBoundaryInformation {
786                            index_original_sender: voxel_plain_index,
787                            index_original_receiver: neighbor_index,
788                            index_original_sender_raw: voxel_index.clone(),
789                        }),
790                }?;
791            }
792        }
793        Ok(())
794    }
795
796    #[cfg_attr(feature = "tracing", instrument(skip_all))]
797    fn update_fluid_mechanics_step_2<ConcGradientExtracellular, ConcTotalExtracellular>(
798        &mut self,
799    ) -> Result<(), SimulationError>
800    where
801        Vox: ExtracellularMechanics<
802                Ind,
803                Pos,
804                ConcVecExtracellular,
805                ConcGradientExtracellular,
806                ConcTotalExtracellular,
807                ConcBoundaryExtracellular,
808            >,
809    {
810        // Receive IndexBoundaryInformation and send back BoundaryConcentrationInformation
811        for index_boundary_information in self.receiver_index.try_iter() {
812            let voxel_box = self
813                .voxels
814                .get(&index_boundary_information.index_original_receiver)
815                .ok_or(IndexError(format!("")))?;
816
817            // Obtain the boundary concentrations to this voxel
818            let concentration_boundary = voxel_box.voxel.boundary_condition_to_neighbor_voxel(
819                &index_boundary_information.index_original_sender_raw,
820            )?;
821
822            // Send back the concentrations here
823            let thread_index =
824                self.plain_index_to_thread[&index_boundary_information.index_original_sender];
825            self.senders_boundary_concentrations[&thread_index].send(
826                ConcentrationBoundaryInformation {
827                    index_original_sender: index_boundary_information.index_original_sender,
828                    index_original_receiver_raw: voxel_box.voxel.get_index(),
829                    concentration_boundary,
830                },
831            )?;
832        }
833
834        Ok(())
835    }
836
837    #[cfg_attr(feature = "tracing", instrument(skip_all))]
838    fn update_fluid_mechanics_step_3<ConcGradientExtracellular, ConcTotalExtracellular>(
839        &mut self,
840        dt: &f64,
841    ) -> Result<(), SimulationError>
842    where
843        ConcVecExtracellular: cellular_raza_concepts::domain_old::Concentration,
844        ConcTotalExtracellular: cellular_raza_concepts::domain_old::Concentration,
845        Vox: ExtracellularMechanics<
846                Ind,
847                Pos,
848                ConcVecExtracellular,
849                ConcGradientExtracellular,
850                ConcTotalExtracellular,
851                ConcBoundaryExtracellular,
852            >,
853    {
854        // Update boundary conditions with new
855        for concentration_boundary_information in self.receiver_concentrations.try_iter() {
856            let vox = self
857                .voxels
858                .get_mut(&concentration_boundary_information.index_original_sender)
859                .ok_or(IndexError(format!(
860                "EngineError: Sender with plain index {} was ended up in location where index is\
861                not present anymore",
862                concentration_boundary_information.index_original_sender)))?;
863            vox.concentration_boundaries.push((
864                concentration_boundary_information.index_original_receiver_raw,
865                concentration_boundary_information.concentration_boundary,
866            ));
867        }
868
869        self.voxels
870            .iter_mut()
871            .map(|(_, voxel_box)| -> Result<(), SimulationError> {
872                let total_extracellular = voxel_box.voxel.get_total_extracellular();
873                let concentration_increment = voxel_box.voxel.calculate_increment(
874                    &total_extracellular,
875                    &voxel_box.extracellular_concentration_increments,
876                    &voxel_box.concentration_boundaries[..],
877                )?;
878                // Update the gradients before we set new extracllular because otherwise it would
879                // be inaccurate. This way the gradients are "behind" the actual concentrations by
880                // one timestep
881                #[cfg(feature = "gradients")]
882                voxel_box
883                    .voxel
884                    .update_extracellular_gradient(&voxel_box.concentration_boundaries[..])?;
885                voxel_box.voxel.set_total_extracellular(
886                    &(total_extracellular + concentration_increment * *dt),
887                )?;
888                voxel_box.extracellular_concentration_increments.drain(..);
889                voxel_box.concentration_boundaries.drain(..);
890                Ok(())
891            })
892            .collect::<Result<(), SimulationError>>()?;
893        Ok(())
894    }
895
896    // TODO the trait bounds here are too harsh. We should not be required to have Pos: Position
897    // or Vel: Velocity here at all!
898    #[cfg_attr(feature = "tracing", instrument(skip_all))]
899    fn update_cellular_mechanics_step_1(&mut self) -> Result<(), SimulationError>
900    where
901        Pos: PositionBound,
902        Vel: VelocityBound,
903        Inf: Clone,
904        For: std::fmt::Debug,
905        Cel: Interaction<Pos, Vel, For, Inf> + Mechanics<Pos, Vel, For> + Clone,
906        Cel: Position<Pos>,
907        Cel: Velocity<Vel>,
908    {
909        use cellular_raza_concepts::InteractionInformation;
910        // General Idea of this function
911        // for each cell
912        //      for each neighbor_voxel in neighbors of voxel containing cell
913        //              if neighbor_voxel is in current MultivoxelContainer
914        //                      calculate forces of current cells on cell and store
915        //                      calculate force from voxel on cell and store
916        //              else
917        //                      send PosInformation to other MultivoxelContainer
918        //
919        // for each PosInformation received from other MultivoxelContainers
920        //      calculate forces of current_cells on cell and send back
921        //
922        // for each ForceInformation received from other MultivoxelContainers
923        //      store received force
924        //
925        // for each cell in this MultiVoxelContainer
926        //      update pos and velocity with all forces obtained
927        //      Simultaneously
928
929        // Calculate forces between cells of own voxel
930        self.voxels
931            .iter_mut()
932            .map(|(_, vox)| vox.calculate_force_between_cells_internally())
933            .collect::<Result<(), CalcError>>()?;
934
935        // Calculate forces for all cells from neighbors
936        // TODO can we do this without memory allocation?
937        let key_iterator: Vec<_> = self.voxels.keys().map(|k| *k).collect();
938
939        for voxel_index in key_iterator {
940            for cell_count in 0..self.voxels[&voxel_index].cells.len() {
941                let cell_pos = self.voxels[&voxel_index].cells[cell_count].0.pos();
942                let cell_vel = self.voxels[&voxel_index].cells[cell_count].0.velocity();
943                let cell_inf = self.voxels[&voxel_index].cells[cell_count]
944                    .0
945                    .get_interaction_information();
946                let mut force = For::zero();
947                let neighbors = self.voxels[&voxel_index].neighbors.clone();
948                for neighbor_index in neighbors {
949                    match self.voxels.get_mut(&neighbor_index) {
950                        Some(vox) => Ok::<(), CalcError>(
951                            force += vox.calculate_force_between_cells_external(
952                                &cell_pos, &cell_vel, &cell_inf,
953                            )?,
954                        ),
955                        None => Ok(
956                            self.senders_pos[&self.plain_index_to_thread[&neighbor_index]].send(
957                                PosInformation {
958                                    index_sender: voxel_index,
959                                    index_receiver: neighbor_index.clone(),
960                                    pos: cell_pos.clone(),
961                                    vel: cell_vel.clone(),
962                                    info: cell_inf.clone(),
963                                    count: cell_count,
964                                },
965                            )?,
966                        ),
967                    }?;
968                }
969                self.voxels.get_mut(&voxel_index).unwrap().cells[cell_count]
970                    .1
971                    .force += force;
972            }
973        }
974
975        // Calculate custom force of voxel on cell
976        self.voxels
977            .iter_mut()
978            .map(|(_, vox)| vox.calculate_custom_force_on_cells())
979            .collect::<Result<(), CalcError>>()?;
980        Ok(())
981    }
982
983    #[cfg_attr(feature = "tracing", instrument(skip_all))]
984    fn update_cellular_mechanics_step_2(&mut self) -> Result<(), SimulationError>
985    where
986        Pos: PositionBound,
987        Vel: VelocityBound,
988        Cel: Interaction<Pos, Vel, For, Inf> + Mechanics<Pos, Vel, For> + Clone,
989        Cel: Position<Pos>,
990        Cel: Velocity<Vel>,
991    {
992        // Receive PositionInformation and send back ForceInformation
993        for pos_info in self.receiver_pos.try_iter() {
994            let vox = self
995                .voxels
996                .get_mut(&pos_info.index_receiver)
997                .ok_or(IndexError(format!(
998                    "EngineError: Voxel with index {:?} of PosInformation can not be found in this\
999                    thread.",
1000                    pos_info.index_receiver
1001                )))?;
1002            // Calculate force from cells in voxel
1003            let force = vox.calculate_force_between_cells_external(
1004                &pos_info.pos,
1005                &pos_info.vel,
1006                &pos_info.info,
1007            )?;
1008
1009            // Send back force information
1010            let thread_index = self.plain_index_to_thread[&pos_info.index_sender];
1011            self.senders_force[&thread_index].send(ForceInformation {
1012                force,
1013                count: pos_info.count,
1014                index_sender: pos_info.index_sender,
1015            })?;
1016        }
1017        Ok(())
1018    }
1019
1020    #[cfg_attr(feature = "tracing", instrument(skip_all))]
1021    fn update_cellular_mechanics_step_3(&mut self, dt: &f64) -> Result<(), SimulationError>
1022    where
1023        Pos: PositionBound,
1024        Vel: VelocityBound,
1025        Cel: Interaction<Pos, Vel, For, Inf> + Mechanics<Pos, Vel, For> + Clone,
1026        Cel: Position<Pos>,
1027        Cel: Velocity<Vel>,
1028    {
1029        // Update position and velocity of all cells with new information
1030        for obt_forces in self.receiver_force.try_iter() {
1031            let vox = self
1032                .voxels
1033                .get_mut(&obt_forces.index_sender)
1034                .ok_or(IndexError(format!(
1035                "EngineError: Sender with plain index {} was ended up in location where index is\
1036                not present anymore", 
1037                obt_forces.index_sender)))?;
1038            match vox.cells.get_mut(obt_forces.count) {
1039                Some((_, aux_storage)) => Ok(aux_storage.force += obt_forces.force),
1040                None => Err(IndexError(format!(
1041                    "EngineError: Force Information with sender index {:?} and cell at vector\
1042                    position {} could not be matched",
1043                    obt_forces.index_sender, obt_forces.count
1044                ))),
1045            }?;
1046        }
1047
1048        // Update position and velocity of cells
1049        for (_, vox) in self.voxels.iter_mut() {
1050            for (cell, aux_storage) in vox.cells.iter_mut() {
1051                // Calculate the current increment
1052                let (dx, dv) = cell.calculate_increment(aux_storage.force.clone())?;
1053                let (dx_rand, dv_rand) = cell.get_random_contribution(&mut vox.rng, *dt)?;
1054
1055                match (
1056                    aux_storage.inc_pos_back_1.clone(),
1057                    aux_storage.inc_pos_back_2.clone(),
1058                    aux_storage.inc_vel_back_1.clone(),
1059                    aux_storage.inc_vel_back_2.clone(),
1060                ) {
1061                    // If all values are present, use the Adams-Bashforth 3rd order
1062                    (
1063                        Some(inc_pos_back_1),
1064                        Some(inc_pos_back_2),
1065                        Some(inc_vel_back_1),
1066                        Some(inc_vel_back_2),
1067                    ) => {
1068                        cell.set_pos(
1069                            &(cell.pos() + dx.clone() * (23.0 / 12.0) * *dt
1070                                - inc_pos_back_1 * (16.0 / 12.0) * *dt
1071                                + inc_pos_back_2 * (5.0 / 12.0) * *dt
1072                                + dx_rand.clone()),
1073                        );
1074                        cell.set_velocity(
1075                            &(cell.velocity() + dv.clone() * (23.0 / 12.0) * *dt
1076                                - inc_vel_back_1 * (16.0 / 12.0) * *dt
1077                                + inc_vel_back_2 * (5.0 / 12.0) * *dt
1078                                + dv_rand.clone()),
1079                        );
1080                    }
1081                    // Otherwise check and use the 2nd order
1082                    (Some(inc_pos_back_1), None, Some(inc_vel_back_1), None) => {
1083                        cell.set_pos(
1084                            &(cell.pos() + dx.clone() * (3.0 / 2.0) * *dt
1085                                - inc_pos_back_1 * (1.0 / 2.0) * *dt
1086                                + dx_rand.clone()),
1087                        );
1088                        cell.set_velocity(
1089                            &(cell.velocity() + dv.clone() * (3.0 / 2.0) * *dt
1090                                - inc_vel_back_1 * (1.0 / 2.0) * *dt
1091                                + dv_rand.clone()),
1092                        );
1093                    }
1094                    // This case should only exists when the cell was first created
1095                    // Then use the Euler Method
1096                    _ => {
1097                        cell.set_pos(&(cell.pos() + dx.clone() * *dt + dx_rand.clone()));
1098                        cell.set_velocity(&(cell.velocity() + dv.clone() * *dt + dv_rand.clone()));
1099                    }
1100                }
1101
1102                // Afterwards update values in auxiliary storage
1103                aux_storage.force = For::zero();
1104                aux_storage.inc_pos_back_1 = Some(dx + dx_rand);
1105                aux_storage.inc_vel_back_1 = Some(dv + dv_rand);
1106            }
1107        }
1108        Ok(())
1109    }
1110
1111    #[cfg_attr(feature = "tracing", instrument(skip_all))]
1112    fn sort_cells_in_voxels_step_1(&mut self) -> Result<(), SimulationError>
1113    where
1114        Pos: PositionBound,
1115        Vel: VelocityBound,
1116        Cel: Mechanics<Pos, Vel, For>,
1117        Cel: Position<Pos>,
1118        Cel: Velocity<Vel>,
1119    {
1120        // Store all cells which need to find a new home in this variable
1121        let mut find_new_home_cells = Vec::<_>::new();
1122
1123        for (voxel_index, vox) in self.voxels.iter_mut() {
1124            // Drain every cell which is currently not in the correct voxel
1125            // TODO use drain_filter when stabilized
1126            let (new_voxel_cells, old_voxel_cells): (Vec<_>, Vec<_>) =
1127                vox.cells.drain(..).partition(|(c, _)| {
1128                    use cellular_raza_concepts::domain_old::Domain;
1129                    match self
1130                        .index_to_plain_index
1131                        .get(&self.domain.get_voxel_index(&c))
1132                    {
1133                        Some(ind) => ind != voxel_index,
1134                        None => panic!("Cannot find index {:?}", self.domain.get_voxel_index(&c)),
1135                    }
1136                });
1137            find_new_home_cells.extend(new_voxel_cells);
1138            vox.cells = old_voxel_cells;
1139        }
1140
1141        // Send cells to other multivoxelcontainer or keep them here
1142        for (cell, aux_storage) in find_new_home_cells {
1143            use cellular_raza_concepts::domain_old::Domain;
1144            let ind = self.domain.get_voxel_index(&cell);
1145            let new_thread_index = self.index_to_thread[&ind];
1146            let cell_index = self.index_to_plain_index[&ind];
1147            match self.voxels.get_mut(&cell_index) {
1148                // If new voxel is in current multivoxelcontainer then save them there
1149                Some(vox) => {
1150                    vox.cells.push((cell, aux_storage));
1151                    Ok(())
1152                }
1153                // Otherwise send them to the correct other multivoxelcontainer
1154                None => match self.senders_cell.get(&new_thread_index) {
1155                    Some(sender) => {
1156                        sender.send((cell, aux_storage))?;
1157                        Ok(())
1158                    }
1159                    None => Err(IndexError(format!(
1160                        "Could not correctly send cell with id {:?}",
1161                        cell.get_id()
1162                    ))),
1163                },
1164            }?;
1165        }
1166        Ok(())
1167    }
1168
1169    #[cfg_attr(feature = "tracing", instrument(skip_all))]
1170    fn sort_cells_in_voxels_step_2(&mut self) -> Result<(), SimulationError> {
1171        // Now receive new cells and insert them
1172        let mut new_cells = self.receiver_cell.try_iter().collect::<Vec<_>>();
1173        for (cell, aux_storage) in new_cells.drain(..) {
1174            self.sort_cell_in_voxel(cell, aux_storage)?;
1175        }
1176        Ok(())
1177    }
1178
1179    #[cfg_attr(feature = "tracing", instrument(skip_all))]
1180    pub(crate) fn save_cells_to_database(
1181        &mut self,
1182        iteration: &u64,
1183    ) -> Result<(), crate::storage::StorageError>
1184    where
1185        Cel: 'static,
1186        CellBox<Cel>: Clone,
1187        AuxiliaryCellPropertyStorage<Pos, Vel, For, ConcVecIntracellular>: Clone,
1188    {
1189        use crate::storage::StorageInterfaceStore;
1190        let cells = self
1191            .voxels
1192            .iter()
1193            .map(|(_, vox)| vox.cells.iter().map(|(c, _)| (c.ref_id(), c)))
1194            .flatten();
1195
1196        self.storage_cells.store_batch_elements(*iteration, cells)
1197    }
1198
1199    #[cfg_attr(feature = "tracing", instrument(skip_all))]
1200    pub(crate) fn save_voxels_to_database(
1201        &mut self,
1202        iteration: &u64,
1203    ) -> Result<(), crate::storage::StorageError>
1204    where
1205        VoxelBox<
1206            Ind,
1207            Pos,
1208            Vel,
1209            For,
1210            Vox,
1211            Cel,
1212            ConcVecExtracellular,
1213            ConcBoundaryExtracellular,
1214            ConcVecIntracellular,
1215        >: Clone + Send + Sync + 'static,
1216    {
1217        use crate::storage::StorageInterfaceStore;
1218        let voxels = self.voxels.iter().map(|(_, voxel)| (voxel.ref_id(), voxel));
1219
1220        self.storage_voxels.store_batch_elements(*iteration, voxels)
1221    }
1222
1223    // TODO find better function signature to have multiple time-scales
1224    // or split into different functions
1225    /// Runs a full update step for all simulation aspects
1226    ///
1227    /// This function is supposed to be run iteratively multiple times in order to
1228    /// completely run a full simulation.
1229    #[cfg_attr(feature = "tracing", instrument(skip_all))]
1230    pub fn run_full_update<ConcGradientExtracellular, ConcTotalExtracellular>(
1231        &mut self,
1232        dt: &f64,
1233    ) -> Result<(), SimulationError>
1234    where
1235        Inf: Send + Sync + core::fmt::Debug,
1236        Pos: PositionBound,
1237        Vel: VelocityBound,
1238        ConcVecExtracellular: cellular_raza_concepts::domain_old::Concentration,
1239        ConcTotalExtracellular: cellular_raza_concepts::domain_old::Concentration,
1240        ConcVecIntracellular: Mul<f64, Output = ConcVecIntracellular>
1241            + Add<ConcVecIntracellular, Output = ConcVecIntracellular>
1242            + AddAssign<ConcVecIntracellular>,
1243        Cel: Cycle<Cel>
1244            + Mechanics<Pos, Vel, For>
1245            + Position<Pos>
1246            + Velocity<Vel>
1247            + Interaction<Pos, Vel, For, Inf>
1248            + CellularReactions<ConcVecIntracellular, ConcVecExtracellular>
1249            + InteractionExtracellularGradient<Cel, ConcGradientExtracellular>
1250            + Volume
1251            + Clone,
1252        Vox: ExtracellularMechanics<
1253                Ind,
1254                Pos,
1255                ConcVecExtracellular,
1256                ConcGradientExtracellular,
1257                ConcTotalExtracellular,
1258                ConcBoundaryExtracellular,
1259            > + Volume,
1260    {
1261        // These methods are used for sending requests and gathering information in general
1262        // This gathers information of forces acting between cells and send between threads
1263
1264        self.update_fluid_mechanics_step_1()?;
1265
1266        // Gather boundary conditions between voxels and domain boundaries and send between threads
1267        self.update_cellular_mechanics_step_1()?;
1268
1269        // Wait for all threads to synchronize.
1270        // The goal is to have as few as possible synchronizations
1271        self.barrier.wait();
1272
1273        self.update_fluid_mechanics_step_2()?;
1274
1275        self.update_cellular_mechanics_step_2()?;
1276
1277        self.barrier.wait();
1278
1279        self.update_cellular_reactions(dt)?;
1280
1281        // These are the true update steps where cell agents are modified the order here may play a
1282        // role!
1283
1284        self.update_fluid_mechanics_step_3(dt)?;
1285
1286        self.update_cellular_mechanics_step_3(dt)?;
1287
1288        // TODO this currently also does application of domain boundaries and inclusion of new cells
1289        // which is wrong in general!
1290        self.update_local_functions(dt)?;
1291
1292        // This function needs an additional synchronization step which cannot correctly be done in
1293        // between the other ones
1294        self.sort_cells_in_voxels_step_1()?;
1295
1296        self.barrier.wait();
1297
1298        self.sort_cells_in_voxels_step_2()?;
1299        Ok(())
1300    }
1301}
cellular_raza_core/backend/cpu_os_threads/domain_decomposition.rs

cellular_raza_core/backend/cpu_os_threads/
domain_decomposition.rs
