cellular_raza_core/backend/chili/
setup.rs

1use cellular_raza_concepts::*;
2
3use serde::{Deserialize, Serialize};
4
5/// Struct containing all necessary information to construct a fully working simulation and run it.
6#[derive(Clone, Deserialize, Serialize)]
7pub struct SimulationSetup<C, D> {
8    /// The initial cells contained in the simulation. In the future we hope to provide an abstract type
9    /// in order to allow for pure iterators to be stored here
10    pub cells: Vec<C>,
11    /// The physical simulation domain which is specified by the [Domain] trait.
12    pub domain: D,
13}
14
15/// Specify settings surrounding execution and storage
16#[derive(Clone, Deserialize, Serialize)]
17pub struct Settings<T, const INIT: bool> {
18    /// Number of threads used for executing simulation in parallel
19    pub n_threads: core::num::NonZeroUsize,
20    // TODO replace this with timestepper in the future
21    /// Specify how time is advanced during the simulation
22    pub time: T,
23    /// Define storage properties
24    pub storage: crate::storage::StorageBuilder<INIT>,
25    /// Determines if progress bar should be shown during execution
26    pub progressbar: Option<String>,
27}
28
29impl<C, D> SimulationSetup<C, D> {
30    /// Insert more cells into the setup after having already initialized the setup.
31    pub fn insert_cells<I>(&mut self, cells: I)
32    where
33        I: IntoIterator<Item = C>,
34    {
35        self.cells.extend(cells);
36    }
37
38    /// Decomposes the struct into a [DecomposedDomain] which can be taken by the backend and turned into multiple subdomains.
39    pub fn decompose<S>(
40        self,
41        n_subdomains: core::num::NonZeroUsize,
42    ) -> Result<DecomposedDomain<D::SubDomainIndex, S, C>, DecomposeError>
43    where
44        D: Domain<C, S>,
45        S: SubDomain,
46    {
47        self.domain.decompose(n_subdomains, self.cells)
48    }
49
50    // TODO add a function which will automatically generate the correct number of subdomains
51    // and simply checks beforehand how many cpu threads are available
52    /// Similar to [decompose](SimulationSetup::decompose) method but does not require to specify
53    /// how many subdomains should be chosen. It will attempt to retrieve resources available to the system
54    /// and spawn threads which are either pre-calculated, read from an existing file or acquired by auto-tuning.
55    pub fn decompose_auto_tune<S>(
56        self,
57    ) -> Result<DecomposedDomain<D::SubDomainIndex, S, C>, DecomposeError>
58    where
59        D: Domain<C, S>,
60        S: SubDomain,
61    {
62        todo!();
63        // let max_n_threads = std::thread::available_parallelism()?;
64        // self.decompose(max_n_threads.into())
65    }
66}
67
68#[cfg(test)]
69mod test {
70    use super::*;
71
72    #[test]
73    fn test_config_insert_cells() {
74        let mut config = SimulationSetup {
75            cells: vec![1, 2, 3, 4],
76            domain: "hello".to_owned(),
77        };
78
79        let new_cells = vec![5, 6, 7];
80        config.insert_cells(new_cells);
81        assert_eq!(config.cells, vec![1, 2, 3, 4, 5, 6, 7]);
82    }
83
84    struct TestDomain {
85        min: f64,
86        max: f64,
87        n_voxels: usize,
88        rng_seed: u64,
89    }
90
91    #[allow(unused)]
92    struct VoxelIndex(usize);
93
94    struct TestSubDomain {
95        min: f64,
96        max: f64,
97        voxels: std::collections::BTreeMap<usize, [f64; 2]>,
98        reflect_at_boundary: (bool, bool),
99        total_voxels: usize,
100    }
101
102    impl Domain<f64, TestSubDomain> for TestDomain {
103        type SubDomainIndex = usize;
104        type VoxelIndex = VoxelIndex;
105
106        fn decompose(
107            self,
108            n_subdomains: core::num::NonZeroUsize,
109            cells: Vec<f64>,
110        ) -> Result<DecomposedDomain<Self::SubDomainIndex, TestSubDomain, f64>, DecomposeError>
111        {
112            let mut cells = cells;
113            let mut index_subdomain_cells = Vec::new();
114            let n_subdomains = n_subdomains.get();
115            let n_return_subdomains = n_subdomains.min(self.n_voxels);
116            let dx = (self.max - self.min) / (self.n_voxels as f64);
117
118            let voxel_distributions = (0..self.n_voxels)
119                .map(|i| (i, (i * n_return_subdomains).div_euclid(self.n_voxels)))
120                .fold(std::collections::BTreeMap::new(), |mut acc, x| {
121                    let entry = acc.entry(x.1).or_insert(Vec::new());
122                    entry.push(x.0);
123                    acc
124                });
125
126            let mut n_total_voxels = 0;
127            for (subdomain_index, voxel_indices) in voxel_distributions {
128                let lower = if subdomain_index == 0 {
129                    self.min
130                } else {
131                    self.min + n_total_voxels as f64 * dx
132                };
133                let upper = if subdomain_index == n_subdomains - 1 {
134                    self.max
135                } else {
136                    self.min + (n_total_voxels + voxel_indices.len()) as f64 * dx
137                };
138                n_total_voxels += voxel_indices.len();
139                let (cells_in_subdomain, other_cells): (Vec<_>, Vec<_>) =
140                    cells.into_iter().partition(|&x| lower <= x && x < upper);
141                cells = other_cells;
142
143                index_subdomain_cells.push((
144                    subdomain_index,
145                    TestSubDomain {
146                        min: lower,
147                        max: upper,
148                        voxels: voxel_indices
149                            .into_iter()
150                            .map(|voxel_index| {
151                                (
152                                    voxel_index,
153                                    [
154                                        self.min + voxel_index as f64 * dx,
155                                        self.min + (voxel_index + 1) as f64 * dx,
156                                    ],
157                                )
158                            })
159                            .collect(),
160                        reflect_at_boundary: (
161                            subdomain_index == 0,
162                            subdomain_index == n_subdomains - 1,
163                        ),
164                        total_voxels: self.n_voxels,
165                    },
166                    cells_in_subdomain,
167                ));
168            }
169
170            let n_subdomains = index_subdomain_cells.len();
171            let decomposed_domain = DecomposedDomain {
172                n_subdomains: n_subdomains.try_into().unwrap(),
173                index_subdomain_cells,
174                neighbor_map: (0..n_subdomains)
175                    .map(|i| {
176                        (
177                            i,
178                            std::collections::BTreeSet::from([
179                                if i == 0 { n_subdomains } else { i - 1 },
180                                i + 1,
181                            ]),
182                        )
183                    })
184                    .collect(),
185                rng_seed: self.rng_seed,
186            };
187            Ok(decomposed_domain)
188        }
189    }
190
191    impl SubDomain for TestSubDomain {
192        type VoxelIndex = usize;
193
194        fn get_all_indices(&self) -> Vec<Self::VoxelIndex> {
195            self.voxels.iter().map(|(&i, _)| i).collect()
196        }
197
198        fn get_neighbor_voxel_indices(
199            &self,
200            voxel_index: &Self::VoxelIndex,
201        ) -> Vec<Self::VoxelIndex> {
202            let mut neighbors = Vec::new();
203            if voxel_index > &0 && voxel_index <= &(self.total_voxels - 1) {
204                neighbors.push(voxel_index - 1);
205            }
206            if voxel_index >= &0 && voxel_index < &(self.total_voxels - 1) {
207                neighbors.push(voxel_index + 1);
208            }
209            neighbors
210        }
211    }
212
213    impl SortCells<f64> for TestSubDomain {
214        type VoxelIndex = usize;
215
216        fn get_voxel_index_of(
217            &self,
218            cell: &f64,
219        ) -> Result<Self::VoxelIndex, cellular_raza_concepts::BoundaryError> {
220            for (index, voxel) in self.voxels.iter() {
221                if cell >= &voxel[0] && cell <= &voxel[1] {
222                    return Ok(*index);
223                }
224            }
225            Err(cellular_raza_concepts::BoundaryError(
226                "Could not find voxel which contains cell".into(),
227            ))
228        }
229    }
230
231    impl SubDomainMechanics<f64, f64> for TestSubDomain {
232        fn apply_boundary(
233            &self,
234            pos: &mut f64,
235            _vel: &mut f64,
236        ) -> Result<(), cellular_raza_concepts::BoundaryError> {
237            if self.reflect_at_boundary.0 && *pos < self.min {
238                *pos = self.min;
239            } else if self.reflect_at_boundary.1 && *pos > self.max {
240                *pos = self.max;
241            }
242            Ok(())
243        }
244    }
245
246    #[test]
247    fn test_config_to_subdomains() {
248        // Define the testdomain with cells
249        let min = 0.0;
250        let max = 100.0;
251        let config = SimulationSetup {
252            cells: vec![1.0, 20.0, 26.0, 41.0, 56.0, 84.0, 95.0],
253            domain: TestDomain {
254                min,
255                max,
256                n_voxels: 8,
257                rng_seed: 0,
258            },
259        };
260
261        let n_subdomains = core::num::NonZeroUsize::new(4).unwrap();
262        let decomposed_domain = config.decompose(n_subdomains).unwrap();
263        for (_, subdomain, cells) in decomposed_domain.index_subdomain_cells.into_iter() {
264            assert!(!cells.is_empty());
265            // Test if each cell is really contained in the subdomain it is supposed to be.
266            for cell in cells {
267                assert!(cell >= subdomain.min);
268                assert!(cell <= subdomain.max);
269            }
270        }
271    }
272
273    #[test]
274    fn test_apply_boundary() {
275        // Define the testdomain with cells
276        let min = 0.0;
277        let max = 100.0;
278        let config = SimulationSetup {
279            cells: vec![1.0, 20.0, 25.0, 50.0, 88.0],
280            domain: TestDomain {
281                min,
282                max,
283                n_voxels: 8,
284                rng_seed: 0,
285            },
286        };
287
288        let n_subdomains = core::num::NonZeroUsize::new(4).unwrap();
289        let decomposed_domain = config.decompose(n_subdomains).unwrap();
290        let mut cell_outside = -10.0;
291        for (_, subdomain, cells) in decomposed_domain.index_subdomain_cells.into_iter() {
292            for cell in cells.into_iter() {
293                let mut cell_prev = cell;
294                let mut _nothing = cell;
295                subdomain
296                    .apply_boundary(&mut cell_prev, &mut _nothing)
297                    .unwrap();
298                assert_eq!(cell_prev, cell);
299            }
300            let mut _nothing = cell_outside;
301            subdomain
302                .apply_boundary(&mut cell_outside, &mut _nothing)
303                .unwrap();
304            assert!(cell_outside >= min);
305            assert!(cell_outside <= max);
306        }
307    }
308
309    #[test]
310    fn test_neighbor_indices() {
311        // Define the testdomain with cells
312        let min = 0.0;
313        let max = 100.0;
314        let config = SimulationSetup {
315            cells: vec![1.0, 20.0, 26.0, 41.0, 56.0, 84.0, 95.0],
316            domain: TestDomain {
317                min,
318                max,
319                n_voxels: 8,
320                rng_seed: 0,
321            },
322        };
323
324        let n_subdomains = core::num::NonZeroUsize::new(4).unwrap();
325        let decomposed_domain = config.decompose(n_subdomains).unwrap();
326        for (_, subdomain, _) in decomposed_domain.index_subdomain_cells.into_iter() {
327            assert_eq!(vec![1], subdomain.get_neighbor_voxel_indices(&0));
328            assert_eq!(vec![0, 2], subdomain.get_neighbor_voxel_indices(&1));
329            assert_eq!(vec![1, 3], subdomain.get_neighbor_voxel_indices(&2));
330            assert_eq!(vec![2, 4], subdomain.get_neighbor_voxel_indices(&3));
331            assert_eq!(vec![3, 5], subdomain.get_neighbor_voxel_indices(&4));
332            assert_eq!(vec![4, 6], subdomain.get_neighbor_voxel_indices(&5));
333            assert_eq!(vec![5, 7], subdomain.get_neighbor_voxel_indices(&6));
334            assert_eq!(vec![6], subdomain.get_neighbor_voxel_indices(&7));
335            assert_eq!(vec![0_usize; 0], subdomain.get_neighbor_voxel_indices(&8));
336            assert_eq!(vec![0_usize; 0], subdomain.get_neighbor_voxel_indices(&9));
337            assert_eq!(vec![0_usize; 0], subdomain.get_neighbor_voxel_indices(&10));
338        }
339    }
340}