Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
Solver.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
16
17#include "shammodels/zeus/Solver.hpp"
18#include "shamalgs/collective/gather_str.hpp"
19#include "shamalgs/collective/reduction.hpp"
20#include "shamcomm/wrapper.hpp"
21#include "shammodels/common/timestep_report.hpp"
22#include "shammodels/zeus/modules/AMRTree.hpp"
23#include "shammodels/zeus/modules/ComputePressure.hpp"
24#include "shammodels/zeus/modules/DiffOperator.hpp"
25#include "shammodels/zeus/modules/FaceFlagger.hpp"
26#include "shammodels/zeus/modules/GhostZones.hpp"
27#include "shammodels/zeus/modules/SourceStep.hpp"
28#include "shammodels/zeus/modules/TransportStep.hpp"
29#include "shammodels/zeus/modules/ValueLoader.hpp"
30#include "shammodels/zeus/modules/WriteBack.hpp"
31#include "shamrock/io/AsciiSplitDump.hpp"
32#include "shamrock/scheduler/SchedulerUtility.hpp"
33
34template<class Tvec, class TgridVec>
35auto shammodels::zeus::Solver<Tvec, TgridVec>::evolve_once(Tscal t_current, Tscal dt_input)
36 -> Tscal {
37
38 StackEntry stack_loc{};
39 sham::MemPerfInfos mem_perf_infos_start = sham::details::get_mem_perf_info();
40 f64 mpi_timer_start = shamcomm::mpi::get_timer("total");
41
42 if (shamcomm::world_rank() == 0) {
43 logger::normal_ln("amr::Zeus", sham::format("t = {}, dt = {}", t_current, dt_input));
44 }
45
46 shambase::Timer tstep;
47 tstep.start();
48
49 scheduler().update_local_load_value([&](shamrock::patch::Patch p) {
50 return scheduler().patch_data.owned_data.get(p.id_patch).get_obj_cnt();
51 });
52
53 SerialPatchTree<TgridVec> _sptree = SerialPatchTree<TgridVec>::build(scheduler());
54 _sptree.attach_buf();
55 storage.serial_patch_tree.set(std::move(_sptree));
56
57 // ghost zone exchange
58 modules::GhostZones gz(context, solver_config, storage);
59 gz.build_ghost_cache();
60
61 gz.exchange_ghost();
62
63 // compute bound received
64 // round to next pow of 2
65 // build radix trees
66 modules::AMRTree amrtree(context, solver_config, storage);
67 amrtree.build_trees();
68
69 amrtree.correct_bounding_box();
70
71 // build neigh table
72 amrtree.build_neigh_cache();
73
74 modules::ComputePressure comp_eos(context, solver_config, storage);
75 comp_eos.compute_p();
76
77 modules::FaceFlagger compute_face_flag(context, solver_config, storage);
78 compute_face_flag.flag_faces();
79 compute_face_flag.split_face_list();
80
81 // modules::DiffOperator diff_op(context,solver_config,storage);
82 // diff_op.compute_gradu();
83
84 using namespace shamrock::patch;
85 using namespace shamrock;
86 using Block = typename Config::AMRBlock;
87 AsciiSplitDump debug_dump(
88 "ghost_dump_debug" + std::to_string(t_current) + std::to_string(solver_config.use_van_leer)
89 + std::to_string(solver_config.use_consistent_transport));
90
91 bool do_debug_dump = false;
92
93 if (do_debug_dump) {
94 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
95 debug_dump.create_id(p.id_patch);
96 });
97
98 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
99 using MergedPDat = shamrock::MergedPatchData;
100 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
101 debug_dump.get_file(p.id_patch).change_table_name("Nobj_original", "u32");
102 debug_dump.get_file(p.id_patch).write_val(mpdat.original_elements);
103 debug_dump.get_file(p.id_patch).change_table_name("Nobj_total", "u32");
104 debug_dump.get_file(p.id_patch).write_val(mpdat.total_elements);
105 });
106
107 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
108 using MergedPDat = shamrock::MergedPatchData;
109 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
110
111 shamrock::patch::PatchDataLayerLayout &ghost_layout
112 = shambase::get_check_ref(storage.ghost_layout.get());
113 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
114 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
115 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
116
117 sham::DeviceBuffer<TgridVec> &cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
118 sham::DeviceBuffer<TgridVec> &cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
119
120 sham::DeviceBuffer<Tscal> &rho_merged
121 = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
122 sham::DeviceBuffer<Tscal> &eint_merged
123 = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
124 sham::DeviceBuffer<Tvec> &vel_merged = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
125
126 debug_dump.get_file(p.id_patch).change_table_name("cell_min", "i64_3");
127 debug_dump.get_file(p.id_patch)
128 .write_table(cell_min.copy_to_stdvec(), mpdat.total_elements);
129 debug_dump.get_file(p.id_patch).change_table_name("cell_max", "i64_3");
130 debug_dump.get_file(p.id_patch)
131 .write_table(cell_max.copy_to_stdvec(), mpdat.total_elements);
132
133 debug_dump.get_file(p.id_patch).change_table_name("rho", "f64");
134 debug_dump.get_file(p.id_patch)
135 .write_table(
136 rho_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
137 debug_dump.get_file(p.id_patch).change_table_name("eint", "f64");
138 debug_dump.get_file(p.id_patch)
139 .write_table(
140 eint_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
141 debug_dump.get_file(p.id_patch).change_table_name("vel", "f64_3");
142 debug_dump.get_file(p.id_patch)
143 .write_table(
144 vel_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
145 });
146 }
147
148 // save velocity field
149 shamrock::patch::PatchDataLayerLayout &ghost_layout
150 = shambase::get_check_ref(storage.ghost_layout.get());
151 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
152 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
153 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
154
155 shamrock::SchedulerUtility utility(scheduler());
156 storage.vel_n.set(
157 utility.save_field_custom<Tvec>("vel_n", [&](u64 id_patch) -> PatchDataField<Tvec> & {
158 using MergedPDat = shamrock::MergedPatchData;
159 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(id_patch);
160 return mpdat.pdat.get_field<Tvec>(ivel_interf);
161 }));
162
163 // prepare velocity gradients
164 modules::ValueLoader<Tvec, TgridVec, Tvec> val_load_vec(context, solver_config, storage);
165 storage.vel_n_xp.set(val_load_vec.load_value_with_gz("vel", {1, 0, 0}, "vel_n_xp"));
166 storage.vel_n_yp.set(val_load_vec.load_value_with_gz("vel", {0, 1, 0}, "vel_n_yp"));
167 storage.vel_n_zp.set(val_load_vec.load_value_with_gz("vel", {0, 0, 1}, "vel_n_zp"));
168
169 modules::ValueLoader<Tvec, TgridVec, Tscal> val_load_scal(context, solver_config, storage);
170 storage.rho_n_xm.set(val_load_scal.load_value_with_gz("rho", {-1, 0, 0}, "rho_n_xm"));
171 storage.rho_n_ym.set(val_load_scal.load_value_with_gz("rho", {0, -1, 0}, "rho_n_ym"));
172 storage.rho_n_zm.set(val_load_scal.load_value_with_gz("rho", {0, 0, -1}, "rho_n_zm"));
173
174 shamrock::ComputeField<Tscal> &pressure_field = storage.pressure.get();
175 storage.pres_n_xm.set(
176 val_load_scal.load_value_with_gz(pressure_field, {-1, 0, 0}, "pres_n_xm"));
177 storage.pres_n_ym.set(
178 val_load_scal.load_value_with_gz(pressure_field, {0, -1, 0}, "pres_n_ym"));
179 storage.pres_n_zm.set(
180 val_load_scal.load_value_with_gz(pressure_field, {0, 0, -1}, "pres_n_zm"));
181
182 modules::SourceStep src_step(context, solver_config, storage);
183 src_step.compute_forces();
184
185 if (do_debug_dump) {
186 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
187 using MergedPDat = shamrock::MergedPatchData;
188 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
189
190 sham::DeviceBuffer<Tvec> &forces_buf = storage.forces.get().get_buf_check(p.id_patch);
191
192 debug_dump.get_file(p.id_patch).change_table_name("force_press", "f64_3");
193 debug_dump.get_file(p.id_patch)
194 .write_table(
195 forces_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
196 });
197 }
198
199 src_step.apply_force(dt_input);
200
201 src_step.compute_AV();
202
203 shamrock::ComputeField<Tvec> &q_AV = storage.q_AV.get();
204 storage.q_AV_n_xm.set(val_load_vec.load_value_with_gz(q_AV, {-1, 0, 0}, "q_AV_n_xm"));
205 storage.q_AV_n_ym.set(val_load_vec.load_value_with_gz(q_AV, {0, -1, 0}, "q_AV_n_ym"));
206 storage.q_AV_n_zm.set(val_load_vec.load_value_with_gz(q_AV, {0, 0, -1}, "q_AV_n_zm"));
207
208 src_step.apply_AV(dt_input);
209 if (do_debug_dump) {
210 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
211 using MergedPDat = shamrock::MergedPatchData;
212 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
213
214 shamrock::patch::PatchDataLayerLayout &ghost_layout
215 = shambase::get_check_ref(storage.ghost_layout.get());
216 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
217 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
218 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
219
220 sham::DeviceBuffer<TgridVec> &cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
221 sham::DeviceBuffer<TgridVec> &cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
222
223 sham::DeviceBuffer<Tscal> &rho_merged
224 = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
225 sham::DeviceBuffer<Tscal> &eint_merged
226 = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
227 sham::DeviceBuffer<Tvec> &vel_merged = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
228
229 debug_dump.get_file(p.id_patch).change_table_name("eint_post_source", "f64");
230 debug_dump.get_file(p.id_patch)
231 .write_table(
232 eint_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
233 debug_dump.get_file(p.id_patch).change_table_name("vel_post_source", "f64_3");
234 debug_dump.get_file(p.id_patch)
235 .write_table(
236 vel_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
237 });
238 }
239
240 src_step.compute_div_v();
241 src_step.update_eint_eos(dt_input);
242
243 if (do_debug_dump) {
244 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
245 using MergedPDat = shamrock::MergedPatchData;
246 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
247
248 sham::DeviceBuffer<Tscal> &divv = storage.div_v_n.get().get_buf_check(p.id_patch);
249
250 debug_dump.get_file(p.id_patch).change_table_name("divv_source", "f64");
251 debug_dump.get_file(p.id_patch)
252 .write_table(divv.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
253 });
254 }
255
256 storage.div_v_n.reset();
257
258 modules::WriteBack wb(context, solver_config, storage);
259 wb.write_back_merged_data();
260
261 storage.merged_patchdata_ghost.reset();
262 storage.ghost_layout.reset();
263
264 storage.vel_n.reset();
265 storage.vel_n_xp.reset();
266 storage.vel_n_yp.reset();
267 storage.vel_n_zp.reset();
268
269 storage.rho_n_xm.reset();
270 storage.rho_n_ym.reset();
271 storage.rho_n_zm.reset();
272
273 storage.pres_n_xm.reset();
274 storage.pres_n_ym.reset();
275 storage.pres_n_zm.reset();
276
277 storage.q_AV.reset();
278 storage.q_AV_n_xm.reset();
279 storage.q_AV_n_ym.reset();
280 storage.q_AV_n_zm.reset();
281
282 // transport step
283 gz.exchange_ghost();
284
285 if (do_debug_dump) {
286 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
287 using MergedPDat = shamrock::MergedPatchData;
288 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
289
290 shamrock::patch::PatchDataLayerLayout &ghost_layout
291 = shambase::get_check_ref(storage.ghost_layout.get());
292 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
293 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
294 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
295
296 sham::DeviceBuffer<TgridVec> &cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
297 sham::DeviceBuffer<TgridVec> &cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
298
299 sham::DeviceBuffer<Tscal> &rho_merged
300 = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
301 sham::DeviceBuffer<Tscal> &eint_merged
302 = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
303 sham::DeviceBuffer<Tvec> &vel_merged = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
304
305 debug_dump.get_file(p.id_patch).change_table_name("eint_start_transp", "f64");
306 debug_dump.get_file(p.id_patch)
307 .write_table(
308 eint_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
309 debug_dump.get_file(p.id_patch).change_table_name("vel_start_transp", "f64_3");
310 debug_dump.get_file(p.id_patch)
311 .write_table(
312 vel_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
313 });
314 }
315
316 modules::ValueLoader<Tvec, TgridVec, Tvec> val_load_vec_v2(context, solver_config, storage);
317 storage.vel_n_xp.set(val_load_vec_v2.load_value_with_gz("vel", {1, 0, 0}, "vel_n_xp"));
318 storage.vel_n_yp.set(val_load_vec_v2.load_value_with_gz("vel", {0, 1, 0}, "vel_n_yp"));
319 storage.vel_n_zp.set(val_load_vec_v2.load_value_with_gz("vel", {0, 0, 1}, "vel_n_zp"));
320
321 if (do_debug_dump) {
322 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
323 using MergedPDat = shamrock::MergedPatchData;
324 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
325
326 sham::DeviceBuffer<Tvec> &vel_n_xp = storage.vel_n_xp.get().get_buf_check(p.id_patch);
327
328 debug_dump.get_file(p.id_patch).change_table_name("vel_n_xp", "f64_3");
329 debug_dump.get_file(p.id_patch)
330 .write_table(
331 vel_n_xp.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
332 });
333 }
334
335 /*
336 using namespace shamrock::patch;
337 using namespace shamrock;
338 using Block = typename Config::AMRBlock;
339 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchData &pdat) {
340
341 using MergedPDat = shamrock::MergedPatchData;
342 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
343
344 sycl::buffer<Tscal> &rho_merged = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
345 sycl::buffer<Tscal> &eint_merged = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
346 sycl::buffer<Tvec> &vel_merged = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
347
348 PatchData &patch_dest = scheduler().patch_data.get_pdat(p.id_patch);
349 sycl::buffer<Tscal> &rho_dest = patch_dest.get_field_buf_ref<Tscal>(irho_interf);
350 sycl::buffer<Tscal> &eint_dest = patch_dest.get_field_buf_ref<Tscal>(ieint_interf);
351 sycl::buffer<Tvec> &vel_dest = patch_dest.get_field_buf_ref<Tvec>(ivel_interf);
352
353
354 sycl::buffer<Tvec> &forces_buf = storage.vel_n_xp.get().get_buf_check(p.id_patch);
355 //sycl::buffer<Tscal> & tmp = storage.pres_n_xm.get().get_buf_check(p.id_patch);
356 //sycl::buffer<sycl::vec<Tscal, 8>> & Q_tmp = storage.Q.get().get_buf_check(p.id_patch);
357 sycl::buffer<Tscal> &buf_p = pressure_field.get_buf_check(p.id_patch);
358
359 shamsys::instance::get_compute_queue().submit([&](sycl::handler & cgh){
360
361 sycl::accessor acc_rho_src{buf_p, cgh, sycl::read_only};
362 sycl::accessor acc_eint_src{eint_merged, cgh, sycl::read_only};
363 sycl::accessor acc_vel_src{vel_merged, cgh, sycl::read_only};
364 sycl::accessor acc_vel_src_xp{forces_buf, cgh, sycl::read_only};
365 //sycl::accessor Q{Q_tmp, cgh, sycl::read_only};
366
367 sycl::accessor acc_rho_dest{rho_dest, cgh, sycl::write_only};
368 sycl::accessor acc_eint_dest{eint_dest, cgh, sycl::write_only};
369 sycl::accessor acc_vel_dest{vel_dest, cgh, sycl::write_only};
370
371 shambase::parallel_for(cgh, mpdat.original_elements*Block::block_size, "tmp copy_ack",
372 [=](u32 id){
373 //acc_rho_dest[id] = acc_rho_src[id];
374 acc_eint_dest[id] = acc_vel_src_xp[id].x();
375 acc_vel_dest[id] = acc_vel_src[id];
376 });
377 });
378
379 if (mpdat.pdat.has_nan()) {
380 logger::err_ln("[Zeus]", "nan detected in write back");
381 throw shambase::make_except_with_loc<std::runtime_error>("detected nan");
382 }
383
384 });
385
386 return 0;
387 */
388
389 modules::TransportStep transport(context, solver_config, storage);
390 transport.compute_cell_centered_momentas();
391
392 if (do_debug_dump) {
393 using Tscal8 = sycl::vec<Tscal, 8>;
394 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
395 using MergedPDat = shamrock::MergedPatchData;
396 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
397
398 sham::DeviceBuffer<Tscal8> &Q_buf = storage.Q.get().get_buf_check(p.id_patch);
399
400 debug_dump.get_file(p.id_patch).change_table_name("Q", "f64_8");
401 debug_dump.get_file(p.id_patch)
402 .write_table(Q_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
403 });
404 }
405
406 storage.vel_n_xp.reset();
407 storage.vel_n_yp.reset();
408 storage.vel_n_zp.reset();
409
410 transport.compute_limiter();
411
412 if (do_debug_dump) {
413 using Tscal8 = sycl::vec<Tscal, 8>;
414 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
415 using MergedPDat = shamrock::MergedPatchData;
416 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
417
418 sham::DeviceBuffer<Tscal8> &ax_buf = storage.a_x.get().get_buf_check(p.id_patch);
419 sham::DeviceBuffer<Tscal8> &ay_buf = storage.a_y.get().get_buf_check(p.id_patch);
420 sham::DeviceBuffer<Tscal8> &az_buf = storage.a_z.get().get_buf_check(p.id_patch);
421
422 debug_dump.get_file(p.id_patch).change_table_name("ax", "f64_8");
423 debug_dump.get_file(p.id_patch)
424 .write_table(ax_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
425 debug_dump.get_file(p.id_patch).change_table_name("ay", "f64_8");
426 debug_dump.get_file(p.id_patch)
427 .write_table(ay_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
428 debug_dump.get_file(p.id_patch).change_table_name("az", "f64_8");
429 debug_dump.get_file(p.id_patch)
430 .write_table(az_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
431 });
432 }
433
434 transport.compute_face_centered_moments(dt_input);
435
436 storage.a_x.reset();
437 storage.a_y.reset();
438 storage.a_z.reset();
439 storage.Q_xm.reset();
440 storage.Q_ym.reset();
441 storage.Q_zm.reset();
442
443 transport.exchange_face_centered_gz();
444
445 if (do_debug_dump) {
446 using Tscal8 = sycl::vec<Tscal, 8>;
447 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
448 using MergedPDat = shamrock::MergedPatchData;
449 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
450
451 sham::DeviceBuffer<Tscal8> &Qstarx_buf
452 = storage.Qstar_x.get().get_buf_check(p.id_patch);
453 sham::DeviceBuffer<Tscal8> &Qstary_buf
454 = storage.Qstar_y.get().get_buf_check(p.id_patch);
455 sham::DeviceBuffer<Tscal8> &Qstarz_buf
456 = storage.Qstar_z.get().get_buf_check(p.id_patch);
457
458 debug_dump.get_file(p.id_patch).change_table_name("Qstar_x", "f64_8");
459 debug_dump.get_file(p.id_patch)
460 .write_table(
461 Qstarx_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
462 debug_dump.get_file(p.id_patch).change_table_name("Qstar_y", "f64_8");
463 debug_dump.get_file(p.id_patch)
464 .write_table(
465 Qstary_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
466 debug_dump.get_file(p.id_patch).change_table_name("Qstar_z", "f64_8");
467 debug_dump.get_file(p.id_patch)
468 .write_table(
469 Qstarz_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
470 });
471 }
472
473 transport.compute_flux();
474
475 if (do_debug_dump) {
476 using Tscal8 = sycl::vec<Tscal, 8>;
477 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
478 using MergedPDat = shamrock::MergedPatchData;
479 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
480
481 sham::DeviceBuffer<Tscal8> &Fluxx_buf = storage.Flux_x.get().get_buf_check(p.id_patch);
482 sham::DeviceBuffer<Tscal8> &Fluxy_buf = storage.Flux_y.get().get_buf_check(p.id_patch);
483 sham::DeviceBuffer<Tscal8> &Fluxz_buf = storage.Flux_z.get().get_buf_check(p.id_patch);
484
485 debug_dump.get_file(p.id_patch).change_table_name("Flux_x", "f64_8");
486 debug_dump.get_file(p.id_patch)
487 .write_table(
488 Fluxx_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
489 debug_dump.get_file(p.id_patch).change_table_name("Flux_y", "f64_8");
490 debug_dump.get_file(p.id_patch)
491 .write_table(
492 Fluxy_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
493 debug_dump.get_file(p.id_patch).change_table_name("Flux_z", "f64_8");
494 debug_dump.get_file(p.id_patch)
495 .write_table(
496 Fluxz_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
497 });
498 }
499
500 transport.compute_stencil_flux();
501
502 transport.update_Q(dt_input);
503
504 transport.compute_new_qte();
505
506 if (do_debug_dump) {
507 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
508 using MergedPDat = shamrock::MergedPatchData;
509 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
510
511 shamrock::patch::PatchDataLayerLayout &ghost_layout
512 = shambase::get_check_ref(storage.ghost_layout.get());
513 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
514 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
515 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
516
517 sham::DeviceBuffer<Tscal> &rho_merged
518 = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
519 sham::DeviceBuffer<Tscal> &eint_merged
520 = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
521 sham::DeviceBuffer<Tvec> &vel_merged = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
522
523 debug_dump.get_file(p.id_patch).change_table_name("rho_end_transp", "f64");
524 debug_dump.get_file(p.id_patch)
525 .write_table(
526 rho_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
527 debug_dump.get_file(p.id_patch).change_table_name("eint_end_transp", "f64");
528 debug_dump.get_file(p.id_patch)
529 .write_table(
530 eint_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
531 debug_dump.get_file(p.id_patch).change_table_name("vel_end_transp", "f64_3");
532 debug_dump.get_file(p.id_patch)
533 .write_table(
534 vel_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
535 });
536 }
537
538 wb.write_back_merged_data();
539
540 storage.Q.reset();
541 storage.Q_xm.reset();
542 storage.Q_ym.reset();
543 storage.Q_zm.reset();
544
545 storage.face_lists.reset();
546 storage.pressure.reset();
547 storage.trees.reset();
548 storage.merge_patch_bounds.reset();
549 storage.merged_patchdata_ghost.reset();
550 storage.ghost_layout.reset();
551 storage.ghost_zone_infos.reset();
552 storage.serial_patch_tree.reset();
553
554 if (do_debug_dump) {
555 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
556 debug_dump.get_file(p.id_patch).close();
557 });
558 }
559
560 tstep.stop();
561
562 sham::MemPerfInfos mem_perf_infos_end = sham::details::get_mem_perf_info();
563
564 f64 delta_mpi_timer = shamcomm::mpi::get_timer("total") - mpi_timer_start;
565 f64 t_dev_alloc
566 = (mem_perf_infos_end.time_alloc_device - mem_perf_infos_start.time_alloc_device)
567 + (mem_perf_infos_end.time_free_device - mem_perf_infos_start.time_free_device);
568 f64 t_host_alloc = (mem_perf_infos_end.time_alloc_host - mem_perf_infos_start.time_alloc_host)
569 + (mem_perf_infos_end.time_free_host - mem_perf_infos_start.time_free_host);
570
571 u64 rank_count = scheduler().get_rank_count() * AMRBlock::block_size;
572 f64 rate = f64(rank_count) / tstep.elapsed_sec();
573
574 u64 npatch = scheduler().patch_list.local.size();
575
576 std::string log_step = report_perf_timestep(
577 rate,
578 rank_count,
579 npatch,
580 tstep.elapsed_sec(),
581 delta_mpi_timer,
582 t_dev_alloc,
583 t_host_alloc,
584 mem_perf_infos_end.max_allocated_byte_device,
585 mem_perf_infos_end.max_allocated_byte_host);
586
587 if (shamcomm::world_rank() == 0) {
588 logger::info_ln("amr::Zeus", log_step);
589 logger::info_ln(
590 "amr::Zeus", "estimated rate :", dt_input * (3600 / tstep.elapsed_sec()), "(tsim/hr)");
591 }
592
593 storage.timings_details.reset();
594
595 // CFL timestep computation for the next iteration
596 PatchDataLayerLayout &pdl = scheduler().pdl_old();
597 const u32 irho = pdl.get_field_idx<Tscal>("rho");
598 const u32 ieint = pdl.get_field_idx<Tscal>("eint");
599 const u32 ivel = pdl.get_field_idx<Tvec>("vel");
600
601 ComputeField<Tscal> cfl_dt = utility.make_compute_field<Tscal>("cfl_dt", AMRBlock::block_size);
602
603 Tscal gamma = solver_config.eos_gamma;
604 Tscal Csafe = solver_config.Csafe;
605 Tscal dxfact = solver_config.grid_coord_to_pos_fact;
606 Tscal one_over_Nside = 1. / AMRBlock::Nside;
607
608 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
609 u32 cell_count = pdat.get_obj_cnt() * AMRBlock::block_size;
610
611 sham::DeviceBuffer<TgridVec> &buf_block_min = pdat.get_field_buf_ref<TgridVec>(0);
612 sham::DeviceBuffer<TgridVec> &buf_block_max = pdat.get_field_buf_ref<TgridVec>(1);
613 sham::DeviceBuffer<Tscal> &buf_rho = pdat.get_field_buf_ref<Tscal>(irho);
614 sham::DeviceBuffer<Tscal> &buf_eint = pdat.get_field_buf_ref<Tscal>(ieint);
615 sham::DeviceBuffer<Tvec> &buf_vel = pdat.get_field_buf_ref<Tvec>(ivel);
616 sham::DeviceBuffer<Tscal> &cfl_dt_buf = cfl_dt.get_buf_check(cur_p.id_patch);
617
618 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
619
620 sham::EventList depends_list;
621 auto acc_cfl_dt = cfl_dt_buf.get_write_access(depends_list);
622 auto acc_block_min = buf_block_min.get_read_access(depends_list);
623 auto acc_block_max = buf_block_max.get_read_access(depends_list);
624 auto acc_rho = buf_rho.get_read_access(depends_list);
625 auto acc_eint = buf_eint.get_read_access(depends_list);
626 auto acc_vel = buf_vel.get_read_access(depends_list);
627
628 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
629 shambase::parallel_for(cgh, cell_count, "compute_cfl_zeus", [=](u64 gid) {
630 const u32 block_id = (u32) gid / AMRBlock::block_size;
631
632 TgridVec lower = acc_block_min[block_id];
633 TgridVec upper = acc_block_max[block_id];
634 Tvec lower_flt = lower.template convert<Tscal>() * dxfact;
635 Tvec upper_flt = upper.template convert<Tscal>() * dxfact;
636 Tvec block_cell_size = (upper_flt - lower_flt) * one_over_Nside;
637 Tscal dx = block_cell_size.x();
638
639 Tscal rho = acc_rho[gid];
640 Tscal eint = acc_eint[gid];
641 Tvec vel = acc_vel[gid];
642
643 Tscal press = (gamma - 1) * eint;
644 Tscal cs = sycl::sqrt(gamma * press / rho);
645
646 constexpr Tscal div = 1.; // this can be lowered if it is too unstable later
647
648 Tscal dt_vel = dx / sycl::length(vel);
649 Tscal dt_cs = dx / cs;
650 Tscal dt
651 = Csafe * div * sycl::rsqrt(sycl::pown(dt_vel, -2) + sycl::pown(dt_cs, -2));
652
653 acc_cfl_dt[gid] = dt;
654 });
655 });
656
657 cfl_dt_buf.complete_event_state(e);
658 buf_block_min.complete_event_state(e);
659 buf_block_max.complete_event_state(e);
660 buf_rho.complete_event_state(e);
661 buf_eint.complete_event_state(e);
662 buf_vel.complete_event_state(e);
663 });
664
665 Tscal rank_dt = cfl_dt.compute_rank_min();
666 Tscal next_dt = shamalgs::collective::allreduce_min(rank_dt);
667
668 if (shamcomm::world_rank() == 0) {
669 logger::info_ln("amr::Zeus", "cfl dt =", next_dt);
670 }
671
672 return next_dt;
673}
674
675template class shammodels::zeus::Solver<f64_3, i64_3>;
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
AsciiSplitDump
A class to dump a simulation state into ASCII files.
Definition AsciiSplitDump.hpp:32
sham::DeviceBuffer
A buffer allocated in USM (Unified Shared Memory).
Definition DeviceBuffer.hpp:106
sham::DeviceBuffer::complete_event_state
void complete_event_state(sycl::event e) const
Complete the event state of the buffer.
Definition DeviceBuffer.hpp:368
sham::DeviceBuffer::copy_to_stdvec
std::vector< T > copy_to_stdvec() const
Copy the content of the buffer to a std::vector.
Definition DeviceBuffer.hpp:510
sham::DeviceBuffer::get_read_access
const T * get_read_access(sham::EventList &depends_list, SourceLocation src_loc=SourceLocation{}) const
Get a read-only pointer to the buffer's data.
Definition DeviceBuffer.hpp:331
sham::DeviceQueue
A SYCL queue associated with a device and a context.
Definition DeviceQueue.hpp:32
sham::DeviceQueue::submit
sycl::event submit(Fct &&fct)
Submits a kernel to the SYCL queue.
Definition DeviceQueue.hpp:101
sham::EventList
Class to manage a list of SYCL events.
Definition EventList.hpp:32
shambase::Timer
Class Timer measures the time elapsed since the timer was started.
Definition Timer.hpp:36
shambase::Timer::elapsed_sec
f64 elapsed_sec() const
Converts the stored nanosecond time to a floating point representation in seconds.
Definition Timer.hpp:88
shambase::Timer::start
void start()
Starts the timer.
Definition Timer.hpp:51
shambase::Timer::stop
void stop()
Stops the timer and stores the elapsed time in nanoseconds.
Definition Timer.hpp:65
shamrock::patch::PatchDataLayerLayout::get_field_idx
u32 get_field_idx(const std::string &field_name) const
Get the field id if matching name & type.
Definition PatchDataLayerLayout.cpp:130
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
gather_str.hpp
MPI string gather / allgather helpers (declarations; implementations in shamalgs/src/collective/gathe...
sham::details::get_mem_perf_info
MemPerfInfos get_mem_perf_info()
Retrieve the memory performance information.
Definition internal_alloc.cpp:123
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
shamcomm::world_rank
i32 world_rank()
Gives the rank of the current process in the MPI communicator.
Definition worldInfo.cpp:41
shamrock
namespace for the main framework
Definition __init__.py:1
shamcomm::logs::normal_ln
void normal_ln(std::string module_name, Types... var2)
Prints a log message with multiple arguments followed by a newline.
Definition logs.hpp:132
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
StackEntry
shambase::details::BasicStackEntry StackEntry
Alias for shambase::details::BasicStackEntry.
Definition stacktrace.hpp:242
sham::MemPerfInfos
Structure to store the performance informations about memory allocation and deallocation.
Definition MemPerfInfos.hpp:30
sham::MemPerfInfos::time_alloc_host
f64 time_alloc_host
Time spent allocating memory on the host.
Definition MemPerfInfos.hpp:32
sham::MemPerfInfos::max_allocated_byte_host
size_t max_allocated_byte_host
max bytes allocated on the host
Definition MemPerfInfos.hpp:53
sham::MemPerfInfos::time_free_device
f64 time_free_device
Time spent deallocating memory on the device.
Definition MemPerfInfos.hpp:41
sham::MemPerfInfos::max_allocated_byte_device
size_t max_allocated_byte_device
max bytes allocated on the device
Definition MemPerfInfos.hpp:55
sham::MemPerfInfos::time_alloc_device
f64 time_alloc_device
Time spent allocating memory on the device.
Definition MemPerfInfos.hpp:34
sham::MemPerfInfos::time_free_host
f64 time_free_host
Time spent deallocating memory on the host.
Definition MemPerfInfos.hpp:39
shamrock::patch::Patch
Patch object that contain generic patch information.
Definition Patch.hpp:33
shamrock::patch::Patch::id_patch
u64 id_patch
unique key that identify the patch
Definition Patch.hpp:86
timestep_report.hpp
shamcomm::mpi::get_timer
f64 get_timer(std::string timername)
get a timer value
Definition wrapper.cpp:46
Generated on for Shamrock by doxygen 1.15.0