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 "shamrock/patch/PatchDataLayer.hpp"
30#include "shamrock/scheduler/DataInserterUtility.hpp"
31#include "shamrock/scheduler/PatchScheduler.hpp"
32#include "shamsys/NodeInstance.hpp"
33#include "shamsys/legacy/log.hpp"
34#include <functional>
35#include <stdexcept>
36#include <utility>
37#include <vector>
38
39template<class Tvec, template<class> class SPHKernel>
40void shammodels::gsph::Model<Tvec, SPHKernel>::init() {
41
42 if (solver.solver_config.scheduler_conf.split_load_value == 0) {
43 throw shambase::make_except_with_loc<std::invalid_argument>(
44 "Scheduler load value should be greater than 0");
45 }
46
47 solver.init_required_fields();
48 ctx.init_sched(
49 solver.solver_config.scheduler_conf.split_load_value,
50 solver.solver_config.scheduler_conf.merge_load_value);
51
52 using namespace shamrock::patch;
53
54 PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
55
56 sched.add_root_patch();
57
58 shamlog_debug_ln("Sys", "build local scheduler tables");
59 sched.owned_patch_id = sched.patch_list.build_local();
60 sched.patch_list.build_local_idx_map();
61 sched.update_local_load_value([&](shamrock::patch::Patch p) {
62 return sched.patch_data.owned_data.get(p.id_patch).get_obj_cnt();
63 });
64 solver.init_ghost_layout();
65
66 solver.init_solver_graph();
67
68 solver.ensure_time_state_edges();
69}
70
71template<class Tvec, template<class> class SPHKernel>
72u64 shammodels::gsph::Model<Tvec, SPHKernel>::get_total_part_count() {
73 PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
74 return shamalgs::collective::allreduce_sum(sched.get_rank_count());
75}
76
77template<class Tvec, template<class> class SPHKernel>
78f64 shammodels::gsph::Model<Tvec, SPHKernel>::total_mass_to_part_mass(f64 totmass) {
79 return totmass / get_total_part_count();
80}
81
82template<class Tvec, template<class> class SPHKernel>
83auto shammodels::gsph::Model<Tvec, SPHKernel>::get_ideal_fcc_box(
84 Tscal dr, std::pair<Tvec, Tvec> box) -> std::pair<Tvec, Tvec> {
85 StackEntry stack_loc{};
86 auto [a, b] = generic::setup::generators::get_ideal_fcc_box<Tscal>(
87 dr, std::make_tuple(box.first, box.second));
88 return {a, b};
89}
90
91template<class Tvec, template<class> class SPHKernel>
92auto shammodels::gsph::Model<Tvec, SPHKernel>::get_ideal_hcp_box(
93 Tscal dr, std::pair<Tvec, Tvec> box) -> std::pair<Tvec, Tvec> {
94 StackEntry stack_loc{};
95 auto [a, b] = generic::setup::generators::get_ideal_fcc_box<Tscal>(
96 dr, std::make_tuple(box.first, box.second));
97 return {a, b};
98}
99
100template<class Tvec, template<class> class SPHKernel>
101void shammodels::gsph::Model<Tvec, SPHKernel>::add_cube_fcc_3d(
102 Tscal dr, std::pair<Tvec, Tvec> _box) {
103 StackEntry stack_loc{};
104
105 shammath::CoordRange<Tvec> box = _box;
106
107 using namespace shamrock::patch;
108
109 PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
110
111 std::string log = "";
112
113 auto make_sliced = [&]() {
114 std::vector<Tvec> vec_lst;
115 generic::setup::generators::add_particles_fcc(
116 dr,
117 std::make_tuple(box.lower, box.upper),
118 [&](Tvec r) {
119 return box.contain_pos(r);
120 },
121 [&](Tvec r, Tscal h) {
122 vec_lst.push_back(r);
123 });
124
125 std::vector<std::vector<Tvec>> sliced_buf;
126
127 u32 sz_buf = sched.crit_patch_split * 4;
128
129 std::vector<Tvec> cur_buf;
130 for (u32 i = 0; i < vec_lst.size(); i++) {
131 cur_buf.push_back(vec_lst[i]);
132
133 if (cur_buf.size() > sz_buf) {
134 sliced_buf.push_back(std::exchange(cur_buf, std::vector<Tvec>{}));
135 }
136 }
137
138 if (cur_buf.size() > 0) {
139 sliced_buf.push_back(std::exchange(cur_buf, std::vector<Tvec>{}));
140 }
141
142 return sliced_buf;
143 };
144
145 std::vector<std::vector<Tvec>> sliced_buf = make_sliced();
146
147 for (std::vector<Tvec> to_ins : sliced_buf) {
148
149 sched.for_each_local_patchdata([&](const Patch &p, PatchDataLayer &pdat) {
150 PatchCoordTransform<Tvec> ptransf
151 = sched.get_sim_box().template get_patch_transform<Tvec>();
152
153 shammath::CoordRange<Tvec> patch_coord = ptransf.to_obj_coord(p);
154
155 std::vector<Tvec> vec_acc;
156 for (Tvec r : to_ins) {
157 if (patch_coord.contain_pos(r)) {
158 vec_acc.push_back(r);
159 }
160 }
161
162 if (vec_acc.size() == 0) {
163 return;
164 }
165
166 log += shambase::format(
167 "\n rank = {} patch id={}, add N={} particles, coords = {} {}",
168 shamcomm::world_rank(),
169 p.id_patch,
170 vec_acc.size(),
171 patch_coord.lower,
172 patch_coord.upper);
173
174 PatchDataLayer tmp(sched.get_layout_ptr_old());
175 tmp.resize(vec_acc.size());
176 tmp.fields_raz();
177
178 {
179 u32 len = vec_acc.size();
180 PatchDataField<Tvec> &f = tmp.template get_field<Tvec>(
181 sched.pdl_old().template get_field_idx<Tvec>("xyz"));
182 sycl::buffer<Tvec> buf(vec_acc.data(), len);
183 f.override(buf, len);
184 }
185
186 {
187 PatchDataField<Tscal> &f = tmp.template get_field<Tscal>(
188 sched.pdl_old().template get_field_idx<Tscal>("hpart"));
189 using Kernel = SPHKernel<Tscal>;
190 f.override(Kernel::hfactd * dr);
191 }
192
193 pdat.insert_elements(tmp);
194 });
195
196 sched.check_patchdata_locality_correctness();
197 sched.scheduler_step(true, true);
198 }
199
200 sched.owned_patch_id = sched.patch_list.build_local();
201 sched.patch_list.build_local_idx_map();
202 sched.update_local_load_value([&](Patch p) {
203 return sched.patch_data.owned_data.get(p.id_patch).get_obj_cnt();
204 });
205
206 shamlog_debug_ln("setup", log);
207}
208
209template<class Tvec, template<class> class SPHKernel>
210void shammodels::gsph::Model<Tvec, SPHKernel>::add_cube_hcp_3d(
211 Tscal dr, std::pair<Tvec, Tvec> _box) {
212 StackEntry stack_loc{};
213
214 shammath::CoordRange<Tvec> box = _box;
215
216 using namespace shamrock::patch;
217
218 PatchScheduler &sched = shambase::get_check_ref(ctx.sched);
219
220 std::string log = "";
221
222 auto make_sliced = [&]() {
223 std::vector<Tvec> vec_lst;
224 generic::setup::generators::add_particles_fcc(
225 dr,
226 std::make_tuple(box.lower, box.upper),
227 [&](Tvec r) {
228 return box.contain_pos(r);
229 },
230 [&](Tvec r, Tscal h) {
231 vec_lst.push_back(r);
232 });
233
234 std::vector<std::vector<Tvec>> sliced_buf;
235
236 u32 sz_buf = sched.crit_patch_split * 4;
237
238 std::vector<Tvec> cur_buf;
239 for (u32 i = 0; i < vec_lst.size(); i++) {
240 cur_buf.push_back(vec_lst[i]);
241
242 if (cur_buf.size() > sz_buf) {
243 sliced_buf.push_back(std::exchange(cur_buf, std::vector<Tvec>{}));
244 }
245 }
246
247 if (cur_buf.size() > 0) {
248 sliced_buf.push_back(std::exchange(cur_buf, std::vector<Tvec>{}));
249 }
250
251 return sliced_buf;
252 };
253
254 std::vector<std::vector<Tvec>> sliced_buf = make_sliced();
255
256 for (std::vector<Tvec> to_ins : sliced_buf) {
257
258 sched.for_each_local_patchdata([&](const Patch &p, PatchDataLayer &pdat) {
259 PatchCoordTransform<Tvec> ptransf
260 = sched.get_sim_box().template get_patch_transform<Tvec>();
261
262 shammath::CoordRange<Tvec> patch_coord = ptransf.to_obj_coord(p);
263
264 std::vector<Tvec> vec_acc;
265 for (Tvec r : to_ins) {
266 if (patch_coord.contain_pos(r)) {
267 vec_acc.push_back(r);
268 }
269 }
270
271 if (vec_acc.size() == 0) {
272 return;
273 }
274
275 log += shambase::format(
276 "\n rank = {} patch id={}, add N={} particles, coords = {} {}",
277 shamcomm::world_rank(),
278 p.id_patch,
279 vec_acc.size(),
280 patch_coord.lower,
281 patch_coord.upper);
282
283 PatchDataLayer tmp(sched.get_layout_ptr_old());
284 tmp.resize(vec_acc.size());
285 tmp.fields_raz();
286
287 {
288 u32 len = vec_acc.size();
289 PatchDataField<Tvec> &f = tmp.template get_field<Tvec>(
290 sched.pdl_old().template get_field_idx<Tvec>("xyz"));
291 sycl::buffer<Tvec> buf(vec_acc.data(), len);
292 f.override(buf, len);
293 }
294
295 {
296 PatchDataField<Tscal> &f = tmp.template get_field<Tscal>(
297 sched.pdl_old().template get_field_idx<Tscal>("hpart"));
298 using Kernel = SPHKernel<Tscal>;
299 f.override(Kernel::hfactd * dr);
300 }
301
302 pdat.insert_elements(tmp);
303 });
304
305 sched.check_patchdata_locality_correctness();
306 sched.scheduler_step(true, true);
307 }
308
309 sched.owned_patch_id = sched.patch_list.build_local();
310 sched.patch_list.build_local_idx_map();
311 sched.update_local_load_value([&](Patch p) {
312 return sched.patch_data.owned_data.get(p.id_patch).get_obj_cnt();
313 });
314
315 shamlog_debug_ln("setup", log);
316}
317
318// Explicit template instantiations for all supported kernel types
319template class shammodels::gsph::Model<f64_3, shammath::M4>;
320template class shammodels::gsph::Model<f64_3, shammath::M6>;
321template class shammodels::gsph::Model<f64_3, shammath::M8>;
322template class shammodels::gsph::Model<f64_3, shammath::C2>;
323template class shammodels::gsph::Model<f64_3, shammath::C4>;
324template 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::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:40
shamrock::patch::PatchDataLayer
PatchDataLayer container class, the layout is described in patchdata_layout.
Definition PatchDataLayer.hpp:35
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:110
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:40
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