|
14 | 14 | #include <memory> |
15 | 15 | #include <utility> |
16 | 16 |
|
17 | | -namespace |
18 | | -{ |
19 | | - template <typename ListType> struct PostInitFunctor |
20 | | - { |
21 | | - |
22 | | - PostInitFunctor(ListType _list, double _new_weigth) |
23 | | - : list(std::move(_list)), new_weigth(_new_weigth) |
24 | | - { |
25 | | - } |
26 | | - |
27 | | - KOKKOS_INLINE_FUNCTION void operator()(std::size_t i_particle) const |
28 | | - { |
29 | | - list._owned_data(i_particle).properties.weight = new_weigth; |
30 | | - } |
31 | | - ListType list; |
32 | | - double new_weigth; |
33 | | - }; |
34 | | -} // namespace |
35 | | - |
36 | 17 | namespace MC |
37 | 18 | { |
38 | 19 |
|
39 | | - constexpr bool uniform_init = true; |
40 | | - |
41 | | - namespace |
42 | | - { |
43 | | - |
44 | | - /** |
45 | | - * @brief helper function to get initial particle weight |
46 | | - */ |
47 | | - inline double get_initial_weight(double scale_factor, |
48 | | - double particle_concentration, |
49 | | - double total_volume, |
50 | | - double mass_cell, |
51 | | - size_t n_particles) |
52 | | - { |
53 | | - // Scale factor is a fine tunning adjustement in case of misprediction of |
54 | | - // particle weight particle_concentration is expected to be the real cell |
55 | | - // concentration in g/L (kg/m3) Total volume is expected to be in m3 As a |
56 | | - // result we can calculate the mass carried by each MC particle |
57 | | - return scale_factor * particle_concentration * total_volume / |
58 | | - static_cast<double>(n_particles) / mass_cell; |
59 | | - } |
| 20 | + // constexpr bool uniform_init = true; |
| 21 | + |
| 22 | + // namespace |
| 23 | + // { |
| 24 | + // /** |
| 25 | + // * @brief helper function to get initial particle weight |
| 26 | + // */ |
| 27 | + // inline double get_initial_weight(double scale_factor, |
| 28 | + // double particle_concentration, |
| 29 | + // double total_volume, |
| 30 | + // double mass_cell, |
| 31 | + // size_t n_particles) |
| 32 | + // { |
| 33 | + // // Scale factor is a fine tunning adjustement in case of misprediction of |
| 34 | + // // particle weight particle_concentration is expected to be the real cell |
| 35 | + // // concentration in g/L (kg/m3) Total volume is expected to be in m3 As a |
| 36 | + // // result we can calculate the mass carried by each MC particle |
| 37 | + // return scale_factor * particle_concentration * total_volume / |
| 38 | + // static_cast<double>(n_particles) / mass_cell; |
| 39 | + // } |
60 | 40 |
|
| 41 | + // /** |
| 42 | + // */ |
| 43 | + // template <ParticleModel Model> |
| 44 | + // void impl_init(std::unique_ptr<MonteCarloUnit>& unit, |
| 45 | + // uint64_t particle_per_process, |
| 46 | + // double& total_mass) |
| 47 | + // { |
| 48 | + |
| 49 | + // auto container = ParticlesContainer<Model>(particle_per_process); |
| 50 | + // auto& list = container.get_compute(); |
| 51 | + |
| 52 | + // auto rng = unit->rng; |
| 53 | + // auto particle_rng = list.rng_instance; |
| 54 | + |
| 55 | + // list.set_allocation_factor(AutoGenerated::default_particle_container_allocation_factor); |
| 56 | + // const auto n_compartments = unit->domain.getNumberCompartments(); |
| 57 | + |
| 58 | + // uint64_t min_c = 0; |
| 59 | + // uint64_t max_c = n_compartments; |
| 60 | + |
| 61 | + // if (!uniform_init) |
| 62 | + // { |
| 63 | + // min_c = 0; |
| 64 | + // max_c = 1; |
| 65 | + // } |
| 66 | + |
| 67 | + // /* |
| 68 | + // * The mass of each cell in the reactor can be calculated after the model initialization. |
| 69 | + // * To ensure that each cell has a unique weight based on the total mass, the following |
| 70 | + // formula |
| 71 | + // * is used: weight = XV / m_tot Where: XV - represents a certain property or value |
| 72 | + // related to |
| 73 | + // * the cell (e.g., volume, particle count, etc.). m_tot - the total mass of the cell or |
| 74 | + // * reactor, which needs to be determined first. |
| 75 | + // * |
| 76 | + // * In order to compute the weight correctly, the initialization process needs to be split |
| 77 | + // into |
| 78 | + // * two phases: |
| 79 | + // * 1. The first phase is to calculate the total mass of the cell (m_tot). |
| 80 | + // * 2. The second phase is to apply the newly calculated weight to the cell using the |
| 81 | + // formula |
| 82 | + // * above. |
| 83 | + // * |
| 84 | + // * This split ensures that the mass is determined before the weight, as the weight is |
| 85 | + // * dependent on the total mass. |
| 86 | + // */ |
| 87 | + |
| 88 | + // Kokkos::parallel_reduce( |
| 89 | + // "mc_init_first", |
| 90 | + // Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, particle_per_process), |
| 91 | + // KOKKOS_LAMBDA(const int i, double& local_mass) { |
| 92 | + // auto p = Particle<Model>(1.); |
| 93 | + // p.properties.id = i; |
| 94 | + // const uint64_t location = rng.uniform_u(min_c, max_c); |
| 95 | + // p.properties.current_container = location; |
| 96 | + // p.init(particle_rng); |
| 97 | + // const double mass_i = p.data.mass(); |
| 98 | + // local_mass += mass_i; |
| 99 | + // list.set(i, std::move(p)); |
| 100 | + // }, |
| 101 | + // total_mass); |
| 102 | + // Kokkos::fence(); |
| 103 | + |
| 104 | + // unit->container = container; |
| 105 | + // } |
| 106 | + |
| 107 | + // } // namespace |
| 108 | + void impl_init(double& total_mass,uint64_t n_particles,MonteCarloUnit& unit,AutoGenerated::ContainerVariant&& container); |
61 | 109 | /** |
| 110 | + * @brief Helper function to initialize a MonteCarloUnit. |
| 111 | + * |
| 112 | + * Since MonteCarloUnit is not a generic type and the model type is resolved |
| 113 | + * at runtime, no constructors are defined to avoid carrying template |
| 114 | + * functions when using the unit. This function wraps the constructor |
| 115 | + * externally, providing a convenient way to initialize a MonteCarloUnit with |
| 116 | + * the appropriate model type and parameters. |
| 117 | + * |
| 118 | + * @tparam Model The particle model type, specified at compile time. |
| 119 | + * |
| 120 | + * |
| 121 | + * @return A unique pointer to the initialized MonteCarloUnit. |
62 | 122 | */ |
63 | 123 | template <ParticleModel Model> |
64 | | - void impl_init(std::unique_ptr<MonteCarloUnit>& unit, |
65 | | - uint64_t particle_per_process, |
66 | | - double& total_mass) |
| 124 | + std::unique_ptr<MonteCarloUnit> init(const ExecInfo& info, |
| 125 | + uint64_t n_particles, |
| 126 | + std::span<double> volumes, |
| 127 | + const NeighborsView<HostSpace>& neighbors, |
| 128 | + double& total_mass) |
67 | 129 | { |
68 | | - |
69 | | - auto container = ParticlesContainer<Model>(particle_per_process); |
70 | | - |
71 | | - // auto& compartments = unit->domain.data(); |
72 | | - // auto n_cells = unit->domain.n_cells(); |
73 | | - auto& list = container.get_compute(); |
74 | | - |
75 | | - auto rng = unit->rng; |
76 | | - auto particle_rng = list.rng_instance; |
77 | | - |
78 | | - list.set_allocation_factor(AutoGenerated::default_particle_container_allocation_factor); |
79 | | - const auto n_compartments = unit->domain.getNumberCompartments(); |
80 | | - |
81 | | - uint64_t min_c = 0; |
82 | | - uint64_t max_c = n_compartments; |
83 | | - |
84 | | - if (!uniform_init) |
85 | | - { |
86 | | - min_c = 0; |
87 | | - max_c = 1; |
88 | | - } |
89 | | - |
90 | | - // double reactor_total_mass = x0 * unit->domain.getTotalVolume(); |
91 | | - |
92 | | - // Kokkos::View<double,ComputeSpace> view_reactor_total_mass("view_reactor_total_mass"); |
93 | | - // Kokkos::deep_copy(view_reactor_total_mass,reactor_total_mass); |
94 | | - |
95 | | - /* |
96 | | - * The mass of each cell in the reactor can be calculated after the model initialization. |
97 | | - * To ensure that each cell has a unique weight based on the total mass, the following formula |
98 | | - * is used: weight = XV / m_tot Where: XV - represents a certain property or value related to |
99 | | - * the cell (e.g., volume, particle count, etc.). m_tot - the total mass of the cell or |
100 | | - * reactor, which needs to be determined first. |
101 | | - * |
102 | | - * In order to compute the weight correctly, the initialization process needs to be split into |
103 | | - * two phases: |
104 | | - * 1. The first phase is to calculate the total mass of the cell (m_tot). |
105 | | - * 2. The second phase is to apply the newly calculated weight to the cell using the formula |
106 | | - * above. |
107 | | - * |
108 | | - * This split ensures that the mass is determined before the weight, as the weight is |
109 | | - * dependent on the total mass. |
110 | | - */ |
111 | | - |
112 | | - Kokkos::parallel_reduce( |
113 | | - "mc_init_first", |
114 | | - Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, particle_per_process), |
115 | | - KOKKOS_LAMBDA(const int i, double& local_mass) { |
116 | | - auto p = Particle<Model>(1.); |
117 | | - p.properties.id = i; |
118 | | - const uint64_t location = rng.uniform_u(min_c, max_c); |
119 | | - p.properties.current_container = location; |
120 | | - |
121 | | - // Kokkos::atomic_increment(&n_cells(location)); |
122 | | - // Kokkos::atomic_increment(&compartments(location).n_cells); |
123 | | - p.init(particle_rng); |
124 | | - const double mass_i = p.data.mass(); |
125 | | - // Add the mass of the particle to the local mass variable |
126 | | - local_mass += mass_i; |
127 | | - // p.properties.weight = view_reactor_total_mass()/mass_i; |
128 | | - // Store the particle in the list |
129 | | - list.set(i, std::move(p)); |
130 | | - }, |
131 | | - total_mass // This is the result that will accumulate after the reduction |
132 | | - ); |
133 | | - Kokkos::fence(); |
134 | | - |
135 | | - unit->container = container; |
| 130 | + auto unit = std::make_unique<MonteCarloUnit>(); |
| 131 | + unit->domain = ReactorDomain(volumes, neighbors); |
| 132 | + auto container = ParticlesContainer<Model>(n_particles); |
| 133 | + impl_init(total_mass, n_particles, *unit, std::move(container)); |
| 134 | + return unit; |
| 135 | + |
| 136 | + // auto& list = container.get_compute(); |
| 137 | + |
| 138 | + // auto rng = unit->rng; |
| 139 | + // auto particle_rng = list.rng_instance; |
| 140 | + |
| 141 | + // list.set_allocation_factor(AutoGenerated::default_particle_container_allocation_factor); |
| 142 | + // const auto n_compartments = unit->domain.getNumberCompartments(); |
| 143 | + |
| 144 | + // uint64_t min_c = 0; |
| 145 | + // uint64_t max_c = n_compartments; |
| 146 | + |
| 147 | + // if (!uniform_init) |
| 148 | + // { |
| 149 | + // min_c = 0; |
| 150 | + // max_c = 1; |
| 151 | + // } |
| 152 | + |
| 153 | + // /* |
| 154 | + // * The mass of each cell in the reactor can be calculated after the model initialization. |
| 155 | + // * To ensure that each cell has a unique weight based on the total mass, the following formula |
| 156 | + // * is used: weight = XV / m_tot Where: XV - represents a certain property or value related to |
| 157 | + // * the cell (e.g., volume, particle count, etc.). m_tot - the total mass of the cell or |
| 158 | + // * reactor, which needs to be determined first. |
| 159 | + // * |
| 160 | + // * In order to compute the weight correctly, the initialization process needs to be split into |
| 161 | + // * two phases: |
| 162 | + // * 1. The first phase is to calculate the total mass of the cell (m_tot). |
| 163 | + // * 2. The second phase is to apply the newly calculated weight to the cell using the formula |
| 164 | + // * above. |
| 165 | + // * |
| 166 | + // * This split ensures that the mass is determined before the weight, as the weight is |
| 167 | + // * dependent on the total mass. |
| 168 | + // */ |
| 169 | + |
| 170 | + // Kokkos::parallel_reduce( |
| 171 | + // "mc_init_first", |
| 172 | + // Kokkos::RangePolicy<Kokkos::DefaultExecutionSpace>(0, n_particles), |
| 173 | + // KOKKOS_LAMBDA(const int i, double& local_mass) { |
| 174 | + // auto p = Particle<Model>(1.); |
| 175 | + // p.properties.id = i; |
| 176 | + // const uint64_t location = rng.uniform_u(min_c, max_c); |
| 177 | + // p.properties.current_container = location; |
| 178 | + // p.init(particle_rng); |
| 179 | + // const double mass_i = p.data.mass(); |
| 180 | + // local_mass += mass_i; |
| 181 | + // list.set(i, std::move(p)); |
| 182 | + // }, |
| 183 | + // total_mass); |
| 184 | + // Kokkos::fence(); |
| 185 | + |
| 186 | + // unit->container = container; |
| 187 | + // return unit; |
136 | 188 | } |
137 | 189 |
|
138 | | - } // namespace |
139 | | - |
140 | | - /** |
141 | | - * @brief Helper function to initialize a MonteCarloUnit. |
142 | | - * |
143 | | - * Since MonteCarloUnit is not a generic type and the model type is resolved |
144 | | - * at runtime, no constructors are defined to avoid carrying template |
145 | | - * functions when using the unit. This function wraps the constructor |
146 | | - * externally, providing a convenient way to initialize a MonteCarloUnit with |
147 | | - * the appropriate model type and parameters. |
148 | | - * |
149 | | - * @tparam Model The particle model type, specified at compile time. |
150 | | - * |
151 | | - * |
152 | | - * @return A unique pointer to the initialized MonteCarloUnit. |
153 | | - */ |
154 | | - template <ParticleModel Model> |
155 | | - std::unique_ptr<MonteCarloUnit> init(const ExecInfo& info, |
156 | | - uint64_t n_particles, |
157 | | - std::span<double> volumes, |
158 | | - const NeighborsView<HostSpace>& neighbors, |
159 | | - double& total_mass) |
160 | | - { |
161 | | - auto unit = std::make_unique<MonteCarloUnit>(); |
162 | | - |
163 | | - unit->domain = ReactorDomain(volumes, neighbors); |
164 | | - |
165 | | - // // Note: use size_t because number of particle represent actually an array size. |
166 | | - // std::size_t particle_per_process = n_particles / info.n_rank; |
167 | | - |
168 | | - // const std::size_t remainder = n_particles % info.n_rank; |
169 | | - // if (remainder != 0 && info.current_rank == info.n_rank - 1) |
170 | | - // { |
171 | | - // particle_per_process += remainder; |
172 | | - // } |
173 | | - // constexpr double scale_factor = 1.; |
174 | | - // const double weight = get_initial_weight( |
175 | | - // scale_factor, x0, unit->domain.getTotalVolume(), initial_mass_cell, n_particles); |
176 | | - |
177 | | - impl_init<Model>(unit, n_particles, total_mass); |
178 | | - |
179 | | - return unit; |
180 | | - } |
| 190 | + void post_init_weight(std::unique_ptr<MonteCarloUnit>& unit, double x0, double total_mass); |
181 | 191 |
|
182 | | - inline void post_init_weight(std::unique_ptr<MonteCarloUnit>& unit, double x0, double total_mass) |
183 | | - { |
184 | | - |
185 | | - // const double new_weight = (x0 * unit->domain.getTotalVolume()) / (total_mass); |
186 | | - auto functor = [total_mass, x0, &unit](auto& container) |
187 | | - { |
188 | | - auto& list = container.get_compute(); |
189 | | - const std::size_t n_p = list.size(); |
190 | | - // const double new_weight = static_cast<double>(n_p)*(x0 * unit->domain.getTotalVolume()) / |
191 | | - // (total_mass); |
192 | | - const double new_weight = (x0 * unit->domain.getTotalVolume()) / (total_mass); |
193 | | - KOKKOS_ASSERT(new_weight > 0); |
194 | | - |
195 | | - Kokkos::View<double, ComputeSpace> view_new_weight("view_new_weight"); |
196 | | - Kokkos::deep_copy(view_new_weight, new_weight); |
197 | | - Kokkos::parallel_for("mc_init_apply", |
198 | | - Kokkos::RangePolicy<ComputeSpace>(0, n_p), |
199 | | - PostInitFunctor(list, new_weight)); |
200 | | - unit->init_weight = new_weight; |
201 | | - }; |
202 | | - |
203 | | - std::visit(functor, unit->container); |
204 | | - } |
205 | | - |
206 | | -} // namespace MC |
| 192 | + } // namespace MC |
207 | 193 |
|
208 | 194 | #endif //__MC_INIT_HPP__ |
0 commit comments