cellular_raza_concepts/
domain.rs

1use serde::{Deserialize, Serialize};
2use std::collections::{BTreeMap, BTreeSet};
3
4use crate::errors::{BoundaryError, DecomposeError};
5
6/// Identifier for voxels used internally to get rid of user-defined ones.
7#[cfg_attr(feature = "pyo3", pyo3::pyclass)]
8#[derive(Clone, Copy, Debug, Deserialize, Hash, PartialEq, Eq, Ord, PartialOrd, Serialize)]
9pub struct VoxelPlainIndex(pub usize);
10
11/// This is mainly used by the simulation_flow::cocmmunicator for testing purposes
12#[allow(unused)]
13#[doc(hidden)]
14impl VoxelPlainIndex {
15    pub fn new(value: usize) -> Self {
16        Self(value)
17    }
18}
19
20#[cfg(feature = "pyo3")]
21#[pyo3::pymethods]
22impl VoxelPlainIndex {
23    #[new]
24    fn __new__(index: usize) -> Self {
25        Self::new(index)
26    }
27}
28
29/// Identifier or subdomains
30#[derive(Clone, Copy, Debug, Deserialize, Hash, PartialEq, Eq, Ord, PartialOrd, Serialize)]
31#[cfg_attr(feature = "pyo3", pyo3::pyclass)]
32pub struct SubDomainPlainIndex(pub usize);
33
34/// Provides an abstraction of the physical total simulation domain.
35///
36/// [cellular_raza](https://github.com/jonaspleyer/cellular_raza) uses domain-decomposition
37/// algorithms to split up the computational workload over multiple physical regions.
38/// That's why the domain itself is mostly responsible for being deconstructed
39/// into smaller [SubDomains](SubDomain) which can then be used to numerically solve our system.
40///
41/// This trait can be automatically implemented when the [SortCells], [DomainRngSeed],
42/// and [DomainCreateSubDomains] are satisfied together with a small number of trait bounds to hash
43/// and compare indices.
44pub trait Domain<C, S, Ci = Vec<C>> {
45    /// Subdomains can be identified by their unique [SubDomainIndex](Domain::SubDomainIndex).
46    /// The backend uses this property to construct a mapping (graph) between subdomains.
47    type SubDomainIndex;
48
49    /// Similarly to the [SubDomainIndex](Domain::SubDomainIndex), voxels can be accessed by
50    /// their unique index. The backend will use this information to construct a mapping
51    /// (graph) between voxels inside their respective subdomains.
52    type VoxelIndex;
53
54    /// Deconstructs the [Domain] into its respective subdomains.
55    ///
56    /// When using the blanket implementation of this function, the following steps are carried
57    /// out:
58    /// Its functionality consists of the following steps:
59    /// 1. Decompose the Domain into [Subdomains](SubDomain)
60    /// 2. Build a neighbor map between [SubDomains](SubDomain)
61    /// 3. Sort cells to their respective [SubDomain]
62    ///
63    /// However, to increase performance or avoid trait bounds, one can also opt to implement this
64    /// trait directly.
65    fn decompose(
66        self,
67        n_subdomains: core::num::NonZeroUsize,
68        cells: Ci,
69    ) -> Result<DecomposedDomain<Self::SubDomainIndex, S, C>, DecomposeError>;
70}
71
72/// Manage the current rng seed of a [Domain]
73pub trait DomainRngSeed {
74    // fn set_rng_seed(&mut self, seed: u64);
75
76    /// Obtains the current rng seed
77    fn get_rng_seed(&self) -> u64;
78}
79
80/// Generate [SubDomains](SubDomain) from an existing [Domain]
81pub trait DomainCreateSubDomains<S> {
82    /// This should always be identical to [Domain::SubDomainIndex].
83    type SubDomainIndex;
84    /// This should always be identical to [Domain::VoxelIndex].
85    type VoxelIndex;
86
87    /// Generates at most `n_subdomains`. This function can also return a lower amount of
88    /// subdomains but never less than 1.
89    fn create_subdomains(
90        &self,
91        n_subdomains: core::num::NonZeroUsize,
92    ) -> Result<
93        impl IntoIterator<Item = (Self::SubDomainIndex, S, Vec<Self::VoxelIndex>)>,
94        DecomposeError,
95    >;
96}
97
98/// Generated by the [decompose](Domain::decompose) method. The backend will know how to
99/// deal with this type and crate a working simulation from it.
100pub struct DecomposedDomain<I, S, C> {
101    /// Number of spawned [SubDomains](SubDomain). This number is guaranteed to be
102    /// smaller or equal to the number may be different to the one given to the
103    /// [Domain::decompose] method.
104    /// Such behaviour can result from not being able to construct as many subdomains as desired.
105    /// Note that this function will attempt to construct more [SubDomains](SubDomain)
106    /// than available CPUs if given a larger number.
107    pub n_subdomains: core::num::NonZeroUsize,
108    /// Vector containing properties of individual [SubDomains](SubDomain).
109    /// Entries are [Domain::SubDomainIndex], [SubDomain], and a vector of cells.
110    // TODO can we use another iterator than Vec<(I, S, Vec<C>)>?
111    pub index_subdomain_cells: Vec<(I, S, Vec<C>)>,
112    /// Encapsulates how the subdomains are linked to each other.
113    /// Eg. two subdomains without any boundary will never appear in each others collection
114    /// of neighbors.
115    /// For the future, we might opt to change to an undirected graph rather than a [BTreeMap].
116    pub neighbor_map: BTreeMap<I, BTreeSet<I>>,
117    /// Initial seed of the simulation for random number generation.
118    pub rng_seed: u64,
119}
120
121/// Subdomains are produced by decomposing a [Domain] into multiple physical regions.
122///
123/// # Derivation
124/// ```
125/// # use cellular_raza_concepts::*;
126/// struct MySubDomain {
127///     x_min: f32,
128///     x_max: f32,
129///     n: usize,
130/// }
131///
132/// impl SubDomain for MySubDomain {
133///     type VoxelIndex = usize;
134///
135///     fn get_neighbor_voxel_indices(
136///         &self,
137///         voxel_index: &Self::VoxelIndex
138///     ) -> Vec<Self::VoxelIndex> {
139///         (voxel_index.saturating_sub(1)..voxel_index.saturating_add(1).min(self.n)+1)
140///             .filter(|k| k!=voxel_index)
141///             .collect()
142///     }
143///
144///     fn get_all_indices(&self) -> Vec<Self::VoxelIndex> {
145///         (0..self.n).collect()
146///     }
147/// }
148///
149/// #[derive(SubDomain)]
150/// struct MyNewSubDomain {
151///     #[Base]
152///     base: MySubDomain,
153/// }
154/// # let _my_sdm = MyNewSubDomain {
155/// #     base: MySubDomain {
156/// #         x_min: -20.0,
157/// #         x_max: -11.0,
158/// #         n: 20,
159/// #     }
160/// # };
161/// # assert_eq!(_my_sdm.get_all_indices(), (0..20).collect::<Vec<_>>());
162/// # assert_eq!(_my_sdm.get_neighbor_voxel_indices(&0), vec![1]);
163/// # assert_eq!(_my_sdm.get_neighbor_voxel_indices(&3), vec![2,4]);
164/// # assert_eq!(_my_sdm.get_neighbor_voxel_indices(&7), vec![6,8]);
165/// ```
166pub trait SubDomain {
167    /// Individual Voxels inside each subdomain can be accessed by this index.
168    type VoxelIndex;
169
170    /// Obtains the neighbor voxels of the specified voxel index. This function behaves similarly
171    /// to [SortCells::get_voxel_index_of] in that it also has to return
172    /// indices which are in other [SubDomains](SubDomain).
173    fn get_neighbor_voxel_indices(&self, voxel_index: &Self::VoxelIndex) -> Vec<Self::VoxelIndex>;
174
175    // fn apply_boundary(&self, cell: &mut C) -> Result<(), BoundaryError>;
176
177    /// Get all voxel indices of this [SubDomain].
178    fn get_all_indices(&self) -> Vec<Self::VoxelIndex>;
179}
180
181/// Assign an [VoxelIndex](SortCells::VoxelIndex) to a given cell.
182///
183/// This trait is used by the [Domain] and [SubDomain] trait to assign a [Domain::SubDomainIndex]
184/// and [SubDomain::VoxelIndex] respectively.
185///
186/// # [SubDomain]
187/// This trait is supposed to return the correct voxel index of the cell
188/// even if this index is inside another [SubDomain].
189/// This restriction might be lifted in the future but is still
190/// required now.
191pub trait SortCells<C> {
192    /// An index which determines to which next smaller unit the cell should be assigned.
193    type VoxelIndex;
194
195    /// If given a cell, we can sort this cell into the corresponding sub unit.
196    fn get_voxel_index_of(&self, cell: &C) -> Result<Self::VoxelIndex, BoundaryError>;
197}
198
199/// Apply boundary conditions to a cells position and velocity.
200///
201/// # Derivation
202/// ```
203/// # use cellular_raza_concepts::*;
204/// # use cellular_raza_concepts::BoundaryError;
205/// struct MyMechanics {
206///     x_min: f64,
207///     x_max: f64,
208/// }
209///
210/// impl SubDomainMechanics<f64, f64> for MyMechanics {
211///     fn apply_boundary(&self, pos: &mut f64, vel: &mut f64) -> Result<(), BoundaryError> {
212///         if *pos < self.x_min {
213///             *vel = vel.abs();
214///         }
215///         if *pos > self.x_max {
216///             *vel = -vel.abs();
217///         }
218///         *pos = pos.clamp(self.x_min, self.x_max);
219///         Ok(())
220///     }
221/// }
222///
223/// #[derive(SubDomain)]
224/// struct MySubDomain {
225///     #[Mechanics]
226///     mechanics: MyMechanics,
227/// }
228/// # let _my_sdm = MySubDomain {
229/// #     mechanics: MyMechanics {
230/// #         x_min: 1.0,
231/// #         x_max: 33.0,
232/// #     }
233/// # };
234/// # let mut pos = 0.0;
235/// # let mut vel = - 0.1;
236/// # _my_sdm.apply_boundary(&mut pos, &mut vel).unwrap();
237/// # assert_eq!(pos, 1.0);
238/// # assert_eq!(vel, 0.1);
239/// ```
240pub trait SubDomainMechanics<Pos, Vel> {
241    /// If the subdomain has boundary conditions, this function will enforce them onto the cells.
242    /// For the future, we plan to replace this function to additionally obtain information
243    /// about the previous and current location of the cell.
244    fn apply_boundary(&self, pos: &mut Pos, vel: &mut Vel) -> Result<(), BoundaryError>;
245}
246
247/// Apply a force on a cell depending on its position and velocity.
248///
249/// # Derivation
250/// ```
251/// # use cellular_raza_concepts::*;
252/// struct MyForce {
253///     damping: f64,
254/// }
255///
256/// impl SubDomainForce<f64, f64, f64, f64> for MyForce {
257///     fn calculate_custom_force(&self, _: &f64, vel: &f64, _: &f64) -> Result<f64, CalcError> {
258///         Ok(- self.damping * vel)
259///     }
260/// }
261///
262/// #[derive(SubDomain)]
263/// struct MySubDomain {
264///     #[Force]
265///     force: MyForce,
266/// }
267/// # let _my_sdm = MySubDomain {
268/// #     force: MyForce {
269/// #         damping: 0.1,
270/// #     }
271/// # };
272/// # let calculated_force = _my_sdm.calculate_custom_force(&0.0, &1.0, &0.0).unwrap();
273/// # assert_eq!(calculated_force, -0.1);
274/// ```
275pub trait SubDomainForce<Pos, Vel, For, Inf> {
276    /// Calculates the force which acts on the cell agent
277    fn calculate_custom_force(
278        &self,
279        pos: &Pos,
280        vel: &Vel,
281        inf: &Inf,
282    ) -> Result<For, crate::CalcError>;
283}
284
285/// Describes extracellular reactions and fluid dynamics
286///
287/// # Derivation
288/// ```
289/// # use cellular_raza_concepts::*;
290///
291/// #[derive(Clone, Debug)]
292/// struct MyReactions<const N: usize> {
293///     values: Vec<f32>,
294///     pos: [f32; N],
295/// }
296///
297/// impl<const N: usize> SubDomainReactions<[f32; N], Vec<f32>, f32> for MyReactions<N> {
298///     type NeighborValue = Vec<f32>;
299///     type BorderInfo = Self;
300///
301///     fn treat_increments<I, J>(
302///         &mut self,
303///         neighbors: I,
304///         sources: J,
305///     ) -> Result<(), CalcError>
306///     where
307///         I: IntoIterator<Item = Self::NeighborValue>,
308///         J: IntoIterator<Item = ([f32; N], Vec<f32>)>,
309///     {
310///         Ok(())
311///     }
312///
313///     fn update_fluid_dynamics(&mut self, dt: f32) -> Result<(), CalcError> {
314///         Ok(())
315///     }
316///
317///     fn get_extracellular_at_pos(&self, pos: &[f32; N]) -> Result<Vec<f32>, CalcError> {
318///         Ok(self.values.clone())
319///     }
320///
321///     fn get_neighbor_value(&self, border_info: Self::BorderInfo) -> Self::NeighborValue {
322///         self.values.clone()
323///     }
324///
325///     fn get_border_info(&self) -> Self::BorderInfo {
326///         self.clone()
327///     }
328/// }
329///
330/// #[derive(SubDomain)]
331/// struct DerivedSubDomain<const N: usize> {
332///     #[Reactions]
333///     reactions: MyReactions<N>,
334/// }
335/// ```
336pub trait SubDomainReactions<Pos, Re, Float> {
337    /// Extracellular value of neighbor
338    type NeighborValue;
339    /// Exchanged information to locate neighboring subdomains.
340    type BorderInfo;
341
342    /// Combines increments which have been obtained by neighbors and cell-sources
343    fn treat_increments<I, J>(&mut self, neighbors: I, sources: J) -> Result<(), crate::CalcError>
344    where
345        I: IntoIterator<Item = Self::NeighborValue>,
346        J: IntoIterator<Item = (Pos, Re)>;
347
348    /// Main update function to calculate new values of extracellular concentrations.
349    fn update_fluid_dynamics(&mut self, dt: Float) -> Result<(), crate::CalcError>;
350
351    /// Obtain extracellular concentrations at given point.
352    fn get_extracellular_at_pos(&self, pos: &Pos) -> Result<Re, crate::CalcError>;
353    /// Obtains the [SubDomainReactions::NeighborValue] which should be sent to the neighbor which
354    /// has exposed the given [SubDomainReactions::BorderInfo].
355    fn get_neighbor_value(&self, border_info: Self::BorderInfo) -> Self::NeighborValue;
356    /// Obtains the [SubDomainReactions::BorderInfo] used to retrieve the
357    /// [SubDomainReactions::NeighborValue].
358    fn get_border_info(&self) -> Self::BorderInfo;
359}
360
361/// This trait derives the different aspects of a [SubDomain].
362///
363/// It serves similarly as the [cellular_raza_concepts_derive::CellAgent] trait to quickly
364/// build new structures from already existing functionality.
365///
366/// | Attribute | Trait | Implemented |
367/// | ---  | --- |:---:|
368/// | `Base` | [SubDomain] | ✅ |
369/// | `SortCells` | [SortCells] | ✅ |
370/// | `Mechanics` | [SubDomainMechanics] | ✅ |
371/// | `Force` | [SubDomainForce] | ✅  |
372/// | `Reactions` | [SubDomainReactions] | ❌ |
373///
374/// # Example Usage
375/// ```
376/// # use cellular_raza_concepts::*;
377/// # struct MySubDomain;
378/// # impl SubDomain for MySubDomain {
379/// #     type VoxelIndex = usize;
380/// #     fn get_neighbor_voxel_indices(&self, voxel_index: &Self::VoxelIndex) -> Vec<usize> {
381/// #         Vec::new()
382/// #     }
383/// #     fn get_all_indices(&self) -> Vec<Self::VoxelIndex> {
384/// #         Vec::new()
385/// #     }
386/// # }
387/// #[derive(SubDomain)]
388/// struct MyDerivedSubDomain {
389///     #[Base]
390///     s: MySubDomain,
391/// }
392/// # let derived_subdomain = MyDerivedSubDomain {
393/// #     s: MySubDomain,
394/// # };
395/// # let all_indices = derived_subdomain.get_all_indices();
396/// # assert_eq!(all_indices.len(), 0);
397/// ```
398#[doc(inline)]
399pub use cellular_raza_concepts_derive::SubDomain;
400
401/// Derives aspects related to the simulation [Domain]
402///
403/// | Attribute | Trait | Implemented |
404/// | -- | --- |:---:|
405/// | `Base` | [Domain] | ✅ |
406/// | `DomainPartialDerive` | [Domain] | ✅ |
407/// | `DomainRngSeed` | [DomainRngSeed] | ✅ |
408/// | `DomainCreateSubDomains` | [DomainCreateSubDomains] | ✅ |
409/// | `SortCells` | [SortCells] | ✅ |
410#[doc(inline)]
411pub use cellular_raza_concepts_derive::Domain;
cellular_raza_concepts/domain.rs

cellular_raza_concepts/
domain.rs
