Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
Model.cpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
17
18#include "shambase/aliases_float.hpp"
19#include "shambase/exception.hpp"
20#include "shambase/memory.hpp"
21#include "shambase/stacktrace.hpp"
22#include "shambase/string.hpp"
23#include "shamcomm/logs.hpp"
24#include "shammath/CoordRange.hpp"
25#include "shammath/crystalLattice.hpp"
26#include "shammath/sphkernels.hpp"
27#include "shammodels/common/setup/generators.hpp"
28#include "shammodels/gsph/Model.hpp"
29#include "shammodels/gsph/modules/ComputeLoadBalanceValue.hpp"
30#include "shamrock/patch/PatchDataLayer.hpp"
31#include "shamrock/scheduler/DataInserterUtility.hpp"
32#include "shamrock/scheduler/PatchScheduler.hpp"
33#include "shamsys/NodeInstance.hpp"
34#include "shamsys/legacy/log.hpp"
35#include <functional>
36#include <stdexcept>
37#include <utility>
38#include <vector>
39
40template<class Tvec, template<class> class SPHKernel>
41void shammodels::gsph::Model<Tvec, SPHKernel>::init() {
42
43 if (solver.solver_config.scheduler_conf.split_load_value == 0) {
44 throw shambase::make_except_with_loc<std::invalid_argument>(
45 "Scheduler load value should be greater than 0");
46 }
47
48 solver.init_required_fields();
49 ctx.init_sched(
50 solver.solver_config.scheduler_conf.split_load_value,
51 solver.solver_config.scheduler_conf.merge_load_value);
52
53 using namespace shamrock::patch;
54
55 PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
56
57 sched.add_root_patch();
58
59 shamlog_debug_ln("Sys", "build local scheduler tables");
60 sched.owned_patch_id = sched.patch_list.build_local();
61 sched.patch_list.build_local_idx_map();
62 sched.update_local_load_value([&](shamrock::patch::Patch p) {
63 return sched.patch_data.owned_data.get(p.id_patch).get_obj_cnt();
64 });
65 solver.init_ghost_layout();
66
67 solver.init_solver_graph();
68
69 solver.ensure_time_state_edges();
70}
71
72template<class Tvec>
73inline void post_insert_data(PatchScheduler &sched) {
74 StackEntry stack_loc{};
75
76 sched.scheduler_step(false, false);
77
78 auto [m, M] = sched.get_box_tranform<Tvec>();
79
80 {
81 StackEntry stack_loc{};
82 SerialPatchTree<Tvec> sptree(
83 sched.patch_tree, sched.get_sim_box().get_patch_transform<Tvec>());
84 shamrock::ReattributeDataUtility reatrib(sched);
85 sptree.attach_buf();
86 reatrib.reatribute_patch_objects(sptree, "xyz");
87 sched.check_patchdata_locality_correctness();
88 }
89
90 sched.scheduler_step(true, true);
91
92 {
93 StackEntry stack_loc{};
94 SerialPatchTree<Tvec> sptree(
95 sched.patch_tree, sched.get_sim_box().get_patch_transform<Tvec>());
96
97 shamrock::ReattributeDataUtility reatrib(sched);
98 sptree.attach_buf();
99 reatrib.reatribute_patch_objects(sptree, "xyz");
100 sched.check_patchdata_locality_correctness();
101 }
102
103 std::string log = "";
104
105 using namespace shamrock::patch;
106
107 u32 smallest_count = u32_max;
108 u32 largest_count = 0;
109
110 sched.for_each_local_patchdata([&](const Patch &p, PatchDataLayer &pdat) {
111 u32 tmp = pdat.get_obj_cnt();
112 smallest_count = sham::min(tmp, smallest_count);
113 largest_count = sham::max(tmp, largest_count);
114 });
115
116 smallest_count = shamalgs::collective::allreduce_min(smallest_count);
117 largest_count = shamalgs::collective::allreduce_max(largest_count);
118
119 if (shamcomm::world_rank() == 0) {
120 logger::info_ln(
121 "Model", "current particle counts : min = ", smallest_count, "max = ", largest_count);
122 }
123}
124
125template<class Tvec, template<class> class SPHKernel>
126u64 shammodels::gsph::Model<Tvec, SPHKernel>::get_total_part_count() {
127 PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
128 return shamalgs::collective::allreduce_sum(sched.get_rank_count());
129}
130
131template<class Tvec, template<class> class SPHKernel>
132f64 shammodels::gsph::Model<Tvec, SPHKernel>::total_mass_to_part_mass(f64 totmass) {
133 return totmass / get_total_part_count();
134}
135
136template<class Tvec, template<class> class SPHKernel>
137auto shammodels::gsph::Model<Tvec, SPHKernel>::get_ideal_fcc_box(
138 Tscal dr, std::pair<Tvec, Tvec> box) -> std::pair<Tvec, Tvec> {
139 StackEntry stack_loc{};
140 auto [a, b] = generic::setup::generators::get_ideal_fcc_box<Tscal>(
141 dr, std::make_tuple(box.first, box.second));
142 return {a, b};
143}
144
145template<class Tvec, template<class> class SPHKernel>
146auto shammodels::gsph::Model<Tvec, SPHKernel>::get_ideal_hcp_box(
147 Tscal dr, std::pair<Tvec, Tvec> box) -> std::pair<Tvec, Tvec> {
148 StackEntry stack_loc{};
149 auto [a, b] = generic::setup::generators::get_ideal_fcc_box<Tscal>(
150 dr, std::make_tuple(box.first, box.second));
151 return {a, b};
152}
153
154template<class Tvec, template<class> class SPHKernel>
155void shammodels::gsph::Model<Tvec, SPHKernel>::add_cube_fcc_3d(
156 Tscal dr, std::pair<Tvec, Tvec> _box) {
157 StackEntry stack_loc{};
158
159 shammath::CoordRange<Tvec> box = _box;
160
161 using namespace shamrock::patch;
162
163 PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
164
165 std::string log = "";
166
167 auto make_sliced = [&]() {
168 std::vector<Tvec> vec_lst;
169 generic::setup::generators::add_particles_fcc(
170 dr,
171 std::make_tuple(box.lower, box.upper),
172 [&](Tvec r) {
173 return box.contain_pos(r);
174 },
175 [&](Tvec r, Tscal h) {
176 vec_lst.push_back(r);
177 });
178
179 std::vector<std::vector<Tvec>> sliced_buf;
180
181 u32 sz_buf = sched.crit_patch_split * 4;
182
183 std::vector<Tvec> cur_buf;
184 for (u32 i = 0; i < vec_lst.size(); i++) {
185 cur_buf.push_back(vec_lst[i]);
186
187 if (cur_buf.size() > sz_buf) {
188 sliced_buf.push_back(std::exchange(cur_buf, std::vector<Tvec>{}));
189 }
190 }
191
192 if (cur_buf.size() > 0) {
193 sliced_buf.push_back(std::exchange(cur_buf, std::vector<Tvec>{}));
194 }
195
196 return sliced_buf;
197 };
198
199 std::vector<std::vector<Tvec>> sliced_buf = make_sliced();
200
201 for (std::vector<Tvec> to_ins : sliced_buf) {
202
203 sched.for_each_local_patchdata([&](const Patch &p, PatchDataLayer &pdat) {
204 PatchCoordTransform<Tvec> ptransf
205 = sched.get_sim_box().template get_patch_transform<Tvec>();
206
207 shammath::CoordRange<Tvec> patch_coord = ptransf.to_obj_coord(p);
208
209 std::vector<Tvec> vec_acc;
210 for (Tvec r : to_ins) {
211 if (patch_coord.contain_pos(r)) {
212 vec_acc.push_back(r);
213 }
214 }
215
216 if (vec_acc.size() == 0) {
217 return;
218 }
219
220 log += sham::format(
221 "\n rank = {} patch id={}, add N={} particles, coords = {} {}",
222 shamcomm::world_rank(),
223 p.id_patch,
224 vec_acc.size(),
225 patch_coord.lower,
226 patch_coord.upper);
227
228 PatchDataLayer tmp(sched.get_layout_ptr_old());
229 tmp.resize(vec_acc.size());
230 tmp.fields_raz();
231
232 {
233 u32 len = vec_acc.size();
234 PatchDataField<Tvec> &f = tmp.template get_field<Tvec>(
235 sched.pdl_old().template get_field_idx<Tvec>("xyz"));
236 sycl::buffer<Tvec> buf(vec_acc.data(), len);
237 f.override(buf, len);
238 }
239
240 {
241 PatchDataField<Tscal> &f = tmp.template get_field<Tscal>(
242 sched.pdl_old().template get_field_idx<Tscal>("hpart"));
243 using Kernel = SPHKernel<Tscal>;
244 f.override(Kernel::hfactd * dr);
245 }
246
247 pdat.insert_elements(tmp);
248 });
249
250 sched.check_patchdata_locality_correctness();
251 modules::ComputeLoadBalanceValue<Tvec, SPHKernel>(ctx, solver.solver_config, solver.storage)
252 .update_load_balancing();
253 post_insert_data<Tvec>(sched);
254 }
255
256 shamlog_debug_ln("setup", log);
257}
258
259template<class Tvec, template<class> class SPHKernel>
260void shammodels::gsph::Model<Tvec, SPHKernel>::add_cube_hcp_3d(
261 Tscal dr, std::pair<Tvec, Tvec> _box) {
262 StackEntry stack_loc{};
263
264 shammath::CoordRange<Tvec> box = _box;
265
266 using namespace shamrock::patch;
267
268 PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
269
270 std::string log = "";
271
272 auto make_sliced = [&]() {
273 std::vector<Tvec> vec_lst;
274 generic::setup::generators::add_particles_fcc(
275 dr,
276 std::make_tuple(box.lower, box.upper),
277 [&](Tvec r) {
278 return box.contain_pos(r);
279 },
280 [&](Tvec r, Tscal h) {
281 vec_lst.push_back(r);
282 });
283
284 std::vector<std::vector<Tvec>> sliced_buf;
285
286 u32 sz_buf = sched.crit_patch_split * 4;
287
288 std::vector<Tvec> cur_buf;
289 for (u32 i = 0; i < vec_lst.size(); i++) {
290 cur_buf.push_back(vec_lst[i]);
291
292 if (cur_buf.size() > sz_buf) {
293 sliced_buf.push_back(std::exchange(cur_buf, std::vector<Tvec>{}));
294 }
295 }
296
297 if (cur_buf.size() > 0) {
298 sliced_buf.push_back(std::exchange(cur_buf, std::vector<Tvec>{}));
299 }
300
301 return sliced_buf;
302 };
303
304 std::vector<std::vector<Tvec>> sliced_buf = make_sliced();
305
306 for (std::vector<Tvec> to_ins : sliced_buf) {
307
308 sched.for_each_local_patchdata([&](const Patch &p, PatchDataLayer &pdat) {
309 PatchCoordTransform<Tvec> ptransf
310 = sched.get_sim_box().template get_patch_transform<Tvec>();
311
312 shammath::CoordRange<Tvec> patch_coord = ptransf.to_obj_coord(p);
313
314 std::vector<Tvec> vec_acc;
315 for (Tvec r : to_ins) {
316 if (patch_coord.contain_pos(r)) {
317 vec_acc.push_back(r);
318 }
319 }
320
321 if (vec_acc.size() == 0) {
322 return;
323 }
324
325 log += sham::format(
326 "\n rank = {} patch id={}, add N={} particles, coords = {} {}",
327 shamcomm::world_rank(),
328 p.id_patch,
329 vec_acc.size(),
330 patch_coord.lower,
331 patch_coord.upper);
332
333 PatchDataLayer tmp(sched.get_layout_ptr_old());
334 tmp.resize(vec_acc.size());
335 tmp.fields_raz();
336
337 {
338 u32 len = vec_acc.size();
339 PatchDataField<Tvec> &f = tmp.template get_field<Tvec>(
340 sched.pdl_old().template get_field_idx<Tvec>("xyz"));
341 sycl::buffer<Tvec> buf(vec_acc.data(), len);
342 f.override(buf, len);
343 }
344
345 {
346 PatchDataField<Tscal> &f = tmp.template get_field<Tscal>(
347 sched.pdl_old().template get_field_idx<Tscal>("hpart"));
348 using Kernel = SPHKernel<Tscal>;
349 f.override(Kernel::hfactd * dr);
350 }
351
352 pdat.insert_elements(tmp);
353 });
354
355 sched.check_patchdata_locality_correctness();
356 modules::ComputeLoadBalanceValue<Tvec, SPHKernel>(ctx, solver.solver_config, solver.storage)
357 .update_load_balancing();
358 post_insert_data<Tvec>(sched);
359 }
360
361 shamlog_debug_ln("setup", log);
362}
363
364// Explicit template instantiations for all supported kernel types
365template class shammodels::gsph::Model<f64_3, shammath::M4>;
366template class shammodels::gsph::Model<f64_3, shammath::M6>;
367template class shammodels::gsph::Model<f64_3, shammath::M8>;
368template class shammodels::gsph::Model<f64_3, shammath::C2>;
369template class shammodels::gsph::Model<f64_3, shammath::C4>;
370template class shammodels::gsph::Model<f64_3, shammath::C6>;
NodeInstance.hpp
Header file describing a Node Instance.
PatchScheduler.hpp
MPI scheduler.
aliases_float.hpp
f64
double f64
Alias for double.
Definition aliases_float.hpp:20
u32
std::uint32_t u32
32 bit unsigned integer
Definition aliases_int.hpp:27
u64
std::uint64_t u64
64 bit unsigned integer
Definition aliases_int.hpp:26
PatchScheduler
The MPI scheduler.
Definition PatchScheduler.hpp:86
PatchScheduler::patch_data
SchedulerPatchData patch_data
handle the data of the patches of the scheduler
Definition PatchScheduler.hpp:104
PatchScheduler::crit_patch_split
u64 crit_patch_split
splitting limit (if load value > crit_patch_split => patch split)
Definition PatchScheduler.hpp:100
PatchScheduler::patch_tree
PatchTree patch_tree
handle the tree structure of the patches
Definition PatchScheduler.hpp:105
PatchScheduler::scheduler_step
void scheduler_step(bool do_split_merge, bool do_load_balancing)
scheduler step
Definition PatchScheduler.cpp:308
PatchScheduler::patch_list
SchedulerPatchList patch_list
handle the list of the patches of the scheduler
Definition PatchScheduler.hpp:103
PatchScheduler::owned_patch_id
std::unordered_set< u64 > owned_patch_id
(owned_patch_id = patch_list.build_local())
Definition PatchScheduler.hpp:109
PatchScheduler::add_root_patch
void add_root_patch()
add patch to the scheduler
Definition PatchScheduler.cpp:193
SchedulerPatchList::build_local
std::unordered_set< u64 > build_local()
select owned patches owned by the node to rebuild local
Definition scheduler_patch_list.cpp:34
SchedulerPatchList::build_local_idx_map
void build_local_idx_map()
recompute id_patch_to_local_idx
Definition scheduler_patch_list.cpp:91
shammodels::gsph::Model
The GSPH Model class.
Definition Model.hpp:63
shammodels::gsph::Model::init
void init()
Initialise the model and all the related data structures (patch scheduler in particular).
Definition Model.cpp:41
shamrock::ReattributeDataUtility
Utility class used to move the objects between patches.
Definition ReattributeDataUtility.hpp:38
shamrock::patch::PatchDataLayer
PatchDataLayer container class, the layout is described in patchdata_layout.
Definition PatchDataLayer.hpp:35
shamrock::patch::PatchDataLayer::get_obj_cnt
u32 get_obj_cnt() const
get the number of objects (particles) stored in this layer
Definition PatchDataLayer.cpp:346
shamrock::patch::SimulationBoxInfo::get_patch_transform
PatchCoordTransform< T > get_patch_transform() const
Get a PatchCoordTransform object that describes the conversion between patch coordinates and domain c...
Definition SimBox.hpp:286
shamrock::scheduler::SchedulerPatchData::owned_data
shambase::DistributedData< PatchData > owned_data
map container for patchdata owned by the current node (layout : id_patch,data)
Definition SchedulerPatchData.hpp:58
exception.hpp
This header file contains utility functions related to exception handling in the code.
Model.hpp
GSPH Model class - high-level interface for GSPH simulations.
shambase::get_check_ref
T & get_check_ref(const std::unique_ptr< T > &ptr, SourceLocation loc=SourceLocation())
Takes a std::unique_ptr and returns a reference to the object it holds. It throws a std::runtime_erro...
Definition memory.hpp:112
shambase::make_except_with_loc
ExcptTypes make_except_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Create an exception with a message and a location.
Definition exception.hpp:105
shamcomm::world_rank
i32 world_rank()
Gives the rank of the current process in the MPI communicator.
Definition worldInfo.cpp:41
u32_max
constexpr u32 u32_max
u32 max value
Definition numeric_limits.hpp:76
shamcomm::logs::info_ln
void info_ln(std::string module_name, Types... var2)
Prints a log message with multiple arguments followed by a newline.
Definition logs.hpp:132
sphkernels.hpp
sph kernels
stacktrace.hpp
This file contains the definition for the stacktrace related functionality.
StackEntry
shambase::details::BasicStackEntry StackEntry
Alias for shambase::details::BasicStackEntry.
Definition stacktrace.hpp:242
shamrock::patch::Patch
Patch object that contain generic patch information.
Definition Patch.hpp:33
Generated on for Shamrock by doxygen 1.15.0