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
19#include "shamcomm/wrapper.hpp"
32
33template<class Tvec, class TgridVec>
34auto shammodels::zeus::Solver<Tvec, TgridVec>::evolve_once(Tscal t_current, Tscal dt_input)
35 -> Tscal {
36
37 StackEntry stack_loc{};
39 f64 mpi_timer_start = shamcomm::mpi::get_timer("total");
40
41 if (shamcomm::world_rank() == 0) {
42 logger::normal_ln("amr::Zeus", shambase::format("t = {}, dt = {}", t_current, dt_input));
43 }
44
45 shambase::Timer tstep;
46 tstep.start();
47
48 scheduler().update_local_load_value([&](shamrock::patch::Patch p) {
49 return scheduler().patch_data.owned_data.get(p.id_patch).get_obj_cnt();
50 });
51
53 _sptree.attach_buf();
54 storage.serial_patch_tree.set(std::move(_sptree));
55
56 // ghost zone exchange
57 modules::GhostZones gz(context, solver_config, storage);
58 gz.build_ghost_cache();
59
60 gz.exchange_ghost();
61
62 // compute bound received
63 // round to next pow of 2
64 // build radix trees
65 modules::AMRTree amrtree(context, solver_config, storage);
66 amrtree.build_trees();
67
68 amrtree.correct_bounding_box();
69
70 // build neigh table
71 amrtree.build_neigh_cache();
72
73 modules::ComputePressure comp_eos(context, solver_config, storage);
74 comp_eos.compute_p();
75
76 modules::FaceFlagger compute_face_flag(context, solver_config, storage);
77 compute_face_flag.flag_faces();
78 compute_face_flag.split_face_list();
79
80 // modules::DiffOperator diff_op(context,solver_config,storage);
81 // diff_op.compute_gradu();
82
83 using namespace shamrock::patch;
84 using namespace shamrock;
85 using Block = typename Config::AMRBlock;
86 AsciiSplitDump debug_dump(
87 "ghost_dump_debug" + std::to_string(t_current) + std::to_string(solver_config.use_van_leer)
88 + std::to_string(solver_config.use_consistent_transport));
89
90 bool do_debug_dump = false;
91
92 if (do_debug_dump) {
93 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
94 debug_dump.create_id(p.id_patch);
95 });
96
97 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
98 using MergedPDat = shamrock::MergedPatchData;
99 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
100 debug_dump.get_file(p.id_patch).change_table_name("Nobj_original", "u32");
101 debug_dump.get_file(p.id_patch).write_val(mpdat.original_elements);
102 debug_dump.get_file(p.id_patch).change_table_name("Nobj_total", "u32");
103 debug_dump.get_file(p.id_patch).write_val(mpdat.total_elements);
104 });
105
106 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
107 using MergedPDat = shamrock::MergedPatchData;
108 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
109
111 = shambase::get_check_ref(storage.ghost_layout.get());
112 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
113 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
114 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
115
116 sham::DeviceBuffer<TgridVec> &cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
117 sham::DeviceBuffer<TgridVec> &cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
118
119 sham::DeviceBuffer<Tscal> &rho_merged
120 = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
121 sham::DeviceBuffer<Tscal> &eint_merged
122 = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
123 sham::DeviceBuffer<Tvec> &vel_merged = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
124
125 debug_dump.get_file(p.id_patch).change_table_name("cell_min", "i64_3");
126 debug_dump.get_file(p.id_patch)
127 .write_table(cell_min.copy_to_stdvec(), mpdat.total_elements);
128 debug_dump.get_file(p.id_patch).change_table_name("cell_max", "i64_3");
129 debug_dump.get_file(p.id_patch)
130 .write_table(cell_max.copy_to_stdvec(), mpdat.total_elements);
131
132 debug_dump.get_file(p.id_patch).change_table_name("rho", "f64");
133 debug_dump.get_file(p.id_patch)
134 .write_table(
135 rho_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
136 debug_dump.get_file(p.id_patch).change_table_name("eint", "f64");
137 debug_dump.get_file(p.id_patch)
138 .write_table(
139 eint_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
140 debug_dump.get_file(p.id_patch).change_table_name("vel", "f64_3");
141 debug_dump.get_file(p.id_patch)
142 .write_table(
143 vel_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
144 });
145 }
146
147 // save velocity field
149 = shambase::get_check_ref(storage.ghost_layout.get());
150 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
151 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
152 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
153
154 shamrock::SchedulerUtility utility(scheduler());
155 storage.vel_n.set(
156 utility.save_field_custom<Tvec>("vel_n", [&](u64 id_patch) -> PatchDataField<Tvec> & {
157 using MergedPDat = shamrock::MergedPatchData;
158 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(id_patch);
159 return mpdat.pdat.get_field<Tvec>(ivel_interf);
160 }));
161
162 // prepare velocity gradients
163 modules::ValueLoader<Tvec, TgridVec, Tvec> val_load_vec(context, solver_config, storage);
164 storage.vel_n_xp.set(val_load_vec.load_value_with_gz("vel", {1, 0, 0}, "vel_n_xp"));
165 storage.vel_n_yp.set(val_load_vec.load_value_with_gz("vel", {0, 1, 0}, "vel_n_yp"));
166 storage.vel_n_zp.set(val_load_vec.load_value_with_gz("vel", {0, 0, 1}, "vel_n_zp"));
167
168 modules::ValueLoader<Tvec, TgridVec, Tscal> val_load_scal(context, solver_config, storage);
169 storage.rho_n_xm.set(val_load_scal.load_value_with_gz("rho", {-1, 0, 0}, "rho_n_xm"));
170 storage.rho_n_ym.set(val_load_scal.load_value_with_gz("rho", {0, -1, 0}, "rho_n_ym"));
171 storage.rho_n_zm.set(val_load_scal.load_value_with_gz("rho", {0, 0, -1}, "rho_n_zm"));
172
173 shamrock::ComputeField<Tscal> &pressure_field = storage.pressure.get();
174 storage.pres_n_xm.set(
175 val_load_scal.load_value_with_gz(pressure_field, {-1, 0, 0}, "pres_n_xm"));
176 storage.pres_n_ym.set(
177 val_load_scal.load_value_with_gz(pressure_field, {0, -1, 0}, "pres_n_ym"));
178 storage.pres_n_zm.set(
179 val_load_scal.load_value_with_gz(pressure_field, {0, 0, -1}, "pres_n_zm"));
180
181 modules::SourceStep src_step(context, solver_config, storage);
182 src_step.compute_forces();
183
184 if (do_debug_dump) {
185 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
186 using MergedPDat = shamrock::MergedPatchData;
187 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
188
189 sham::DeviceBuffer<Tvec> &forces_buf = storage.forces.get().get_buf_check(p.id_patch);
190
191 debug_dump.get_file(p.id_patch).change_table_name("force_press", "f64_3");
192 debug_dump.get_file(p.id_patch)
193 .write_table(
194 forces_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
195 });
196 }
197
198 src_step.apply_force(dt_input);
199
200 src_step.compute_AV();
201
202 shamrock::ComputeField<Tvec> &q_AV = storage.q_AV.get();
203 storage.q_AV_n_xm.set(val_load_vec.load_value_with_gz(q_AV, {-1, 0, 0}, "q_AV_n_xm"));
204 storage.q_AV_n_ym.set(val_load_vec.load_value_with_gz(q_AV, {0, -1, 0}, "q_AV_n_ym"));
205 storage.q_AV_n_zm.set(val_load_vec.load_value_with_gz(q_AV, {0, 0, -1}, "q_AV_n_zm"));
206
207 src_step.apply_AV(dt_input);
208 if (do_debug_dump) {
209 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
210 using MergedPDat = shamrock::MergedPatchData;
211 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
212
214 = shambase::get_check_ref(storage.ghost_layout.get());
215 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
216 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
217 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
218
219 sham::DeviceBuffer<TgridVec> &cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
220 sham::DeviceBuffer<TgridVec> &cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
221
222 sham::DeviceBuffer<Tscal> &rho_merged
223 = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
224 sham::DeviceBuffer<Tscal> &eint_merged
225 = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
226 sham::DeviceBuffer<Tvec> &vel_merged = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
227
228 debug_dump.get_file(p.id_patch).change_table_name("eint_post_source", "f64");
229 debug_dump.get_file(p.id_patch)
230 .write_table(
231 eint_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
232 debug_dump.get_file(p.id_patch).change_table_name("vel_post_source", "f64_3");
233 debug_dump.get_file(p.id_patch)
234 .write_table(
235 vel_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
236 });
237 }
238
239 src_step.compute_div_v();
240 src_step.update_eint_eos(dt_input);
241
242 if (do_debug_dump) {
243 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
244 using MergedPDat = shamrock::MergedPatchData;
245 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
246
247 sham::DeviceBuffer<Tscal> &divv = storage.div_v_n.get().get_buf_check(p.id_patch);
248
249 debug_dump.get_file(p.id_patch).change_table_name("divv_source", "f64");
250 debug_dump.get_file(p.id_patch)
251 .write_table(divv.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
252 });
253 }
254
255 storage.div_v_n.reset();
256
257 modules::WriteBack wb(context, solver_config, storage);
258 wb.write_back_merged_data();
259
260 storage.merged_patchdata_ghost.reset();
261 storage.ghost_layout.reset();
262
263 storage.vel_n.reset();
264 storage.vel_n_xp.reset();
265 storage.vel_n_yp.reset();
266 storage.vel_n_zp.reset();
267
268 storage.rho_n_xm.reset();
269 storage.rho_n_ym.reset();
270 storage.rho_n_zm.reset();
271
272 storage.pres_n_xm.reset();
273 storage.pres_n_ym.reset();
274 storage.pres_n_zm.reset();
275
276 storage.q_AV.reset();
277 storage.q_AV_n_xm.reset();
278 storage.q_AV_n_ym.reset();
279 storage.q_AV_n_zm.reset();
280
281 // transport step
282 gz.exchange_ghost();
283
284 if (do_debug_dump) {
285 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
286 using MergedPDat = shamrock::MergedPatchData;
287 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
288
290 = shambase::get_check_ref(storage.ghost_layout.get());
291 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
292 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
293 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
294
295 sham::DeviceBuffer<TgridVec> &cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
296 sham::DeviceBuffer<TgridVec> &cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
297
298 sham::DeviceBuffer<Tscal> &rho_merged
299 = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
300 sham::DeviceBuffer<Tscal> &eint_merged
301 = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
302 sham::DeviceBuffer<Tvec> &vel_merged = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
303
304 debug_dump.get_file(p.id_patch).change_table_name("eint_start_transp", "f64");
305 debug_dump.get_file(p.id_patch)
306 .write_table(
307 eint_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
308 debug_dump.get_file(p.id_patch).change_table_name("vel_start_transp", "f64_3");
309 debug_dump.get_file(p.id_patch)
310 .write_table(
311 vel_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
312 });
313 }
314
315 modules::ValueLoader<Tvec, TgridVec, Tvec> val_load_vec_v2(context, solver_config, storage);
316 storage.vel_n_xp.set(val_load_vec_v2.load_value_with_gz("vel", {1, 0, 0}, "vel_n_xp"));
317 storage.vel_n_yp.set(val_load_vec_v2.load_value_with_gz("vel", {0, 1, 0}, "vel_n_yp"));
318 storage.vel_n_zp.set(val_load_vec_v2.load_value_with_gz("vel", {0, 0, 1}, "vel_n_zp"));
319
320 if (do_debug_dump) {
321 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
322 using MergedPDat = shamrock::MergedPatchData;
323 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
324
325 sham::DeviceBuffer<Tvec> &vel_n_xp = storage.vel_n_xp.get().get_buf_check(p.id_patch);
326
327 debug_dump.get_file(p.id_patch).change_table_name("vel_n_xp", "f64_3");
328 debug_dump.get_file(p.id_patch)
329 .write_table(
330 vel_n_xp.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
331 });
332 }
333
334 /*
335 using namespace shamrock::patch;
336 using namespace shamrock;
337 using Block = typename Config::AMRBlock;
338 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchData &pdat) {
339
340 using MergedPDat = shamrock::MergedPatchData;
341 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
342
343 sycl::buffer<Tscal> &rho_merged = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
344 sycl::buffer<Tscal> &eint_merged = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
345 sycl::buffer<Tvec> &vel_merged = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
346
347 PatchData &patch_dest = scheduler().patch_data.get_pdat(p.id_patch);
348 sycl::buffer<Tscal> &rho_dest = patch_dest.get_field_buf_ref<Tscal>(irho_interf);
349 sycl::buffer<Tscal> &eint_dest = patch_dest.get_field_buf_ref<Tscal>(ieint_interf);
350 sycl::buffer<Tvec> &vel_dest = patch_dest.get_field_buf_ref<Tvec>(ivel_interf);
351
352
353 sycl::buffer<Tvec> &forces_buf = storage.vel_n_xp.get().get_buf_check(p.id_patch);
354 //sycl::buffer<Tscal> & tmp = storage.pres_n_xm.get().get_buf_check(p.id_patch);
355 //sycl::buffer<sycl::vec<Tscal, 8>> & Q_tmp = storage.Q.get().get_buf_check(p.id_patch);
356 sycl::buffer<Tscal> &buf_p = pressure_field.get_buf_check(p.id_patch);
357
358 shamsys::instance::get_compute_queue().submit([&](sycl::handler & cgh){
359
360 sycl::accessor acc_rho_src{buf_p, cgh, sycl::read_only};
361 sycl::accessor acc_eint_src{eint_merged, cgh, sycl::read_only};
362 sycl::accessor acc_vel_src{vel_merged, cgh, sycl::read_only};
363 sycl::accessor acc_vel_src_xp{forces_buf, cgh, sycl::read_only};
364 //sycl::accessor Q{Q_tmp, cgh, sycl::read_only};
365
366 sycl::accessor acc_rho_dest{rho_dest, cgh, sycl::write_only};
367 sycl::accessor acc_eint_dest{eint_dest, cgh, sycl::write_only};
368 sycl::accessor acc_vel_dest{vel_dest, cgh, sycl::write_only};
369
370 shambase::parallel_for(cgh, mpdat.original_elements*Block::block_size, "tmp copy_ack",
371 [=](u32 id){
372 //acc_rho_dest[id] = acc_rho_src[id];
373 acc_eint_dest[id] = acc_vel_src_xp[id].x();
374 acc_vel_dest[id] = acc_vel_src[id];
375 });
376 });
377
378 if (mpdat.pdat.has_nan()) {
379 logger::err_ln("[Zeus]", "nan detected in write back");
380 throw shambase::make_except_with_loc<std::runtime_error>("detected nan");
381 }
382
383 });
384
385 return 0;
386 */
387
388 modules::TransportStep transport(context, solver_config, storage);
389 transport.compute_cell_centered_momentas();
390
391 if (do_debug_dump) {
392 using Tscal8 = sycl::vec<Tscal, 8>;
393 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
394 using MergedPDat = shamrock::MergedPatchData;
395 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
396
397 sham::DeviceBuffer<Tscal8> &Q_buf = storage.Q.get().get_buf_check(p.id_patch);
398
399 debug_dump.get_file(p.id_patch).change_table_name("Q", "f64_8");
400 debug_dump.get_file(p.id_patch)
401 .write_table(Q_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
402 });
403 }
404
405 storage.vel_n_xp.reset();
406 storage.vel_n_yp.reset();
407 storage.vel_n_zp.reset();
408
409 transport.compute_limiter();
410
411 if (do_debug_dump) {
412 using Tscal8 = sycl::vec<Tscal, 8>;
413 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
414 using MergedPDat = shamrock::MergedPatchData;
415 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
416
417 sham::DeviceBuffer<Tscal8> &ax_buf = storage.a_x.get().get_buf_check(p.id_patch);
418 sham::DeviceBuffer<Tscal8> &ay_buf = storage.a_y.get().get_buf_check(p.id_patch);
419 sham::DeviceBuffer<Tscal8> &az_buf = storage.a_z.get().get_buf_check(p.id_patch);
420
421 debug_dump.get_file(p.id_patch).change_table_name("ax", "f64_8");
422 debug_dump.get_file(p.id_patch)
423 .write_table(ax_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
424 debug_dump.get_file(p.id_patch).change_table_name("ay", "f64_8");
425 debug_dump.get_file(p.id_patch)
426 .write_table(ay_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
427 debug_dump.get_file(p.id_patch).change_table_name("az", "f64_8");
428 debug_dump.get_file(p.id_patch)
429 .write_table(az_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
430 });
431 }
432
433 transport.compute_face_centered_moments(dt_input);
434
435 storage.a_x.reset();
436 storage.a_y.reset();
437 storage.a_z.reset();
438 storage.Q_xm.reset();
439 storage.Q_ym.reset();
440 storage.Q_zm.reset();
441
442 transport.exchange_face_centered_gz();
443
444 if (do_debug_dump) {
445 using Tscal8 = sycl::vec<Tscal, 8>;
446 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
447 using MergedPDat = shamrock::MergedPatchData;
448 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
449
451 = storage.Qstar_x.get().get_buf_check(p.id_patch);
453 = storage.Qstar_y.get().get_buf_check(p.id_patch);
455 = storage.Qstar_z.get().get_buf_check(p.id_patch);
456
457 debug_dump.get_file(p.id_patch).change_table_name("Qstar_x", "f64_8");
458 debug_dump.get_file(p.id_patch)
459 .write_table(
460 Qstarx_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
461 debug_dump.get_file(p.id_patch).change_table_name("Qstar_y", "f64_8");
462 debug_dump.get_file(p.id_patch)
463 .write_table(
464 Qstary_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
465 debug_dump.get_file(p.id_patch).change_table_name("Qstar_z", "f64_8");
466 debug_dump.get_file(p.id_patch)
467 .write_table(
468 Qstarz_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
469 });
470 }
471
472 transport.compute_flux();
473
474 if (do_debug_dump) {
475 using Tscal8 = sycl::vec<Tscal, 8>;
476 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
477 using MergedPDat = shamrock::MergedPatchData;
478 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
479
480 sham::DeviceBuffer<Tscal8> &Fluxx_buf = storage.Flux_x.get().get_buf_check(p.id_patch);
481 sham::DeviceBuffer<Tscal8> &Fluxy_buf = storage.Flux_y.get().get_buf_check(p.id_patch);
482 sham::DeviceBuffer<Tscal8> &Fluxz_buf = storage.Flux_z.get().get_buf_check(p.id_patch);
483
484 debug_dump.get_file(p.id_patch).change_table_name("Flux_x", "f64_8");
485 debug_dump.get_file(p.id_patch)
486 .write_table(
487 Fluxx_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
488 debug_dump.get_file(p.id_patch).change_table_name("Flux_y", "f64_8");
489 debug_dump.get_file(p.id_patch)
490 .write_table(
491 Fluxy_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
492 debug_dump.get_file(p.id_patch).change_table_name("Flux_z", "f64_8");
493 debug_dump.get_file(p.id_patch)
494 .write_table(
495 Fluxz_buf.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
496 });
497 }
498
499 transport.compute_stencil_flux();
500
501 transport.update_Q(dt_input);
502
503 transport.compute_new_qte();
504
505 if (do_debug_dump) {
506 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
507 using MergedPDat = shamrock::MergedPatchData;
508 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
509
511 = shambase::get_check_ref(storage.ghost_layout.get());
512 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
513 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
514 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
515
516 sham::DeviceBuffer<Tscal> &rho_merged
517 = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
518 sham::DeviceBuffer<Tscal> &eint_merged
519 = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
520 sham::DeviceBuffer<Tvec> &vel_merged = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
521
522 debug_dump.get_file(p.id_patch).change_table_name("rho_end_transp", "f64");
523 debug_dump.get_file(p.id_patch)
524 .write_table(
525 rho_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
526 debug_dump.get_file(p.id_patch).change_table_name("eint_end_transp", "f64");
527 debug_dump.get_file(p.id_patch)
528 .write_table(
529 eint_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
530 debug_dump.get_file(p.id_patch).change_table_name("vel_end_transp", "f64_3");
531 debug_dump.get_file(p.id_patch)
532 .write_table(
533 vel_merged.copy_to_stdvec(), mpdat.total_elements * AMRBlock::block_size);
534 });
535 }
536
537 wb.write_back_merged_data();
538
539 storage.Q.reset();
540 storage.Q_xm.reset();
541 storage.Q_ym.reset();
542 storage.Q_zm.reset();
543
544 storage.face_lists.reset();
545 storage.pressure.reset();
546 storage.trees.reset();
547 storage.merge_patch_bounds.reset();
548 storage.merged_patchdata_ghost.reset();
549 storage.ghost_layout.reset();
550 storage.ghost_zone_infos.reset();
551 storage.serial_patch_tree.reset();
552
553 if (do_debug_dump) {
554 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
555 debug_dump.get_file(p.id_patch).close();
556 });
557 }
558
559 tstep.end();
560
562
563 f64 delta_mpi_timer = shamcomm::mpi::get_timer("total") - mpi_timer_start;
564 f64 t_dev_alloc
565 = (mem_perf_infos_end.time_alloc_device - mem_perf_infos_start.time_alloc_device)
566 + (mem_perf_infos_end.time_free_device - mem_perf_infos_start.time_free_device);
567 f64 t_host_alloc = (mem_perf_infos_end.time_alloc_host - mem_perf_infos_start.time_alloc_host)
568 + (mem_perf_infos_end.time_free_host - mem_perf_infos_start.time_free_host);
569
570 u64 rank_count = scheduler().get_rank_count() * AMRBlock::block_size;
571 f64 rate = f64(rank_count) / tstep.elasped_sec();
572
573 u64 npatch = scheduler().patch_list.local.size();
574
575 std::string log_step = report_perf_timestep(
576 rate,
577 rank_count,
578 npatch,
579 tstep.elasped_sec(),
580 delta_mpi_timer,
581 t_dev_alloc,
582 t_host_alloc,
583 mem_perf_infos_end.max_allocated_byte_device,
584 mem_perf_infos_end.max_allocated_byte_host);
585
586 if (shamcomm::world_rank() == 0) {
587 logger::info_ln("amr::Zeus", log_step);
588 logger::info_ln(
589 "amr::Zeus", "estimated rate :", dt_input * (3600 / tstep.elasped_sec()), "(tsim/hr)");
590 }
591
592 storage.timings_details.reset();
593
594 return 0;
595}
596
double f64
Alias for double.
std::uint32_t u32
32 bit unsigned integer
std::uint64_t u64
64 bit unsigned integer
A class to dump a simulation state into ASCII files.
A buffer allocated in USM (Unified Shared Memory)
std::vector< T > copy_to_stdvec() const
Copy the content of the buffer to a std::vector.
Class Timer measures the time elapsed since the timer was started.
Definition time.hpp:96
void end()
Stops the timer and stores the elapsed time in nanoseconds.
Definition time.hpp:111
f64 elasped_sec() const
Converts the stored nanosecond time to a floating point representation in seconds.
Definition time.hpp:123
void start()
Starts the timer.
Definition time.hpp:106
u32 get_field_idx(const std::string &field_name) const
Get the field id if matching name & type.
PatchDataLayer container class, the layout is described in patchdata_layout.
MPI string gather / allgather helpers (declarations; implementations in shamalgs/src/collective/gathe...
MemPerfInfos get_mem_perf_info()
Retrieve the memory performance information.
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
i32 world_rank()
Gives the rank of the current process in the MPI communicator.
Definition worldInfo.cpp:40
namespace for the main framework
Definition __init__.py:1
Structure to store the performance informations about memory allocation and deallocation.
f64 time_alloc_host
Time spent allocating memory on the host.
size_t max_allocated_byte_host
max bytes allocated on the host
f64 time_free_device
Time spent deallocating memory on the device.
size_t max_allocated_byte_device
max bytes allocated on the device
f64 time_alloc_device
Time spent allocating memory on the device.
f64 time_free_host
Time spent deallocating memory on the host.
Patch object that contain generic patch information.
Definition Patch.hpp:33
f64 get_timer(std::string timername)
get a timer value
Definition wrapper.cpp:44