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
18
20#include "shambase/memory.hpp"
22#include "shamcomm/logs.hpp"
69#include <memory>
70
71template<class Tvec, class TgridVec>
72class PatchDataLayerToVtk : public shamrock::solvergraph::INode {
73 bool write_id_patch;
74 bool write_world_rank;
75 using Tscal = shambase::VecComponent<Tvec>;
76 u32 block_size;
77
78 public:
79 PatchDataLayerToVtk(bool write_id_patch, bool write_world_rank, u32 block_size)
80 : write_id_patch(write_id_patch), write_world_rank(write_world_rank),
81 block_size(block_size) {}
82
83 struct Edges {
84 // inputs
86 const shamrock::solvergraph::IPatchDataLayerRefs &patch_data_layers;
87 };
88
89 inline void set_edges(
90 std::shared_ptr<shamrock::solvergraph::IDataEdge<std::string>> filename,
91 std::shared_ptr<shamrock::solvergraph::IPatchDataLayerRefs> patch_data_layers) {
92 __internal_set_ro_edges({filename, patch_data_layers});
94 }
95
96 inline Edges get_edges() {
97 return Edges{
100 };
101 }
102
105
106 auto edges = get_edges();
107
108 auto &filename = edges.filename;
109 auto &patch_data_layers = edges.patch_data_layers;
110
111 // Compute the number of fields to generate
112 auto get_field_count = [&]() {
113 u32 field_count = 0;
114
115 {
116 u64 id_patch = patch_data_layers.get_const_refs().get_ids().front();
117 auto &pdat = patch_data_layers.get(id_patch);
118
119 pdat.for_each_field_any([&](auto &field) {
120 field_count++;
121 });
122 }
123
124 if (write_id_patch) {
125 field_count++;
126 }
127 if (write_world_rank) {
128 field_count++;
129 }
130
131 return field_count - 2; // to remove the block infos
132 };
133
134 auto get_layout = [&]() -> const shamrock::patch::PatchDataLayerLayout & {
135 u64 id_patch = patch_data_layers.get_const_refs().get_ids().front();
136 const shamrock::patch::PatchDataLayer &pdat = patch_data_layers.get(id_patch);
137 return pdat.pdl();
138 };
139
140 shamrock::LegacyVtkWriter writer(filename.data, true, shamrock::UnstructuredGrid);
141
142 u32 field_count = get_field_count();
143
144 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
145 auto &q = shambase::get_check_ref(dev_sched).get_queue();
146
147 sham::DeviceBuffer<TgridVec> pos_min_block(0, dev_sched);
148 sham::DeviceBuffer<TgridVec> pos_max_block(0, dev_sched);
149
150 patch_data_layers.get_const_refs().for_each(
151 [&](u64 id_patch, const std::reference_wrapper<shamrock::patch::PatchDataLayer> &pdat) {
152 auto &pdat_ref = pdat.get();
153 auto &buf_pos_min = pdat_ref.get_field_buf_ref<TgridVec>(0);
154 auto &buf_pos_max = pdat_ref.get_field_buf_ref<TgridVec>(1);
155 pos_min_block.append(buf_pos_min);
156 pos_max_block.append(buf_pos_max);
157 });
158
159 u64 num_obj = pos_min_block.get_size();
160
161 sham::DeviceBuffer<Tvec> pos_max_cell(num_obj * block_size, dev_sched);
162 sham::DeviceBuffer<Tvec> pos_min_cell(num_obj * block_size, dev_sched);
163
164 if (num_obj > 0) {
165
167
168 if (Block::block_size != block_size) {
170 "block_size mismatch, got {} expected {}", Block::block_size, block_size));
171 }
172
174 q,
175 sham::MultiRef{pos_min_block, pos_max_block},
176 sham::MultiRef{pos_min_cell, pos_max_cell},
177 num_obj,
178 [](u32 id_a,
179 const TgridVec *__restrict ptr_block_min,
180 const TgridVec *__restrict ptr_block_max,
181 Tvec *cell_min,
182 Tvec *cell_max) {
183 Tvec block_min = ptr_block_min[id_a].template convert<Tscal>();
184 Tvec block_max = ptr_block_max[id_a].template convert<Tscal>();
185
186 Tvec delta_cell = (block_max - block_min) / Block::side_size;
187 for (u32 ix = 0; ix < Block::side_size; ix++) {
188 for (u32 iy = 0; iy < Block::side_size; iy++) {
189 for (u32 iz = 0; iz < Block::side_size; iz++) {
190 u32 i = Block::get_index({ix, iy, iz});
191 Tvec delta_val = delta_cell * Tvec{ix, iy, iz};
192 cell_min[id_a * Block::block_size + i] = block_min + delta_val;
193 cell_max[id_a * Block::block_size + i]
194 = block_min + (delta_cell) + delta_val;
195 }
196 }
197 }
198 });
199 }
200
201 auto pos_min_cell_buf = pos_min_cell.copy_to_sycl_buffer();
202 auto pos_max_cell_buf = pos_max_cell.copy_to_sycl_buffer();
203 writer.write_voxel_cells(pos_min_cell_buf, pos_max_cell_buf, num_obj * block_size);
204
205 writer.add_cell_data_section();
206 writer.add_field_data_section(field_count);
207
208 const shamrock::patch::PatchDataLayerLayout &layout = get_layout();
209
210 layout.for_each_field_any([&](auto &field_desc) {
211 using f_t = typename std::remove_reference<decltype(field_desc)>::type::field_T;
212 u32 nvar = field_desc.nvar;
213 std::string field_name = field_desc.name;
214
215 u32 idx = layout.get_field_idx<f_t>(field_name);
216
217 if (nvar == 1) {
218 // this the block info and i'll skip it for now
219 } else if (nvar != block_size) {
221 } else {
222 sham::DeviceBuffer<f_t> data(0, dev_sched);
223
224 patch_data_layers.get_const_refs().for_each(
225 [&](u64 id_patch,
226 const std::reference_wrapper<shamrock::patch::PatchDataLayer> &pdat) {
227 auto &pdat_ref = pdat.get();
228 auto &buf_field = pdat_ref.get_field_buf_ref<f_t>(idx);
229 data.append(buf_field);
230 });
231
232 auto tmp_buf = data.copy_to_sycl_buffer();
233 writer.write_field(field_name, tmp_buf, num_obj * block_size);
234 }
235 });
236
237 if (write_id_patch) {
238 using f_t = u32;
239 sham::DeviceBuffer<f_t> data(0, dev_sched);
240
241 patch_data_layers.get_const_refs().for_each(
242 [&](u64 id_patch,
243 const std::reference_wrapper<shamrock::patch::PatchDataLayer> &pdat) {
244 auto buf_field
245 = sham::DeviceBuffer<f_t>(pdat.get().get_obj_cnt() * block_size, dev_sched);
246 buf_field.fill(id_patch);
247 data.append(buf_field);
248 });
249
250 auto tmp_buf = data.copy_to_sycl_buffer();
251 writer.write_field("id_patch", tmp_buf, num_obj * block_size);
252 }
253 if (write_world_rank) {
254 using f_t = u32;
255 sham::DeviceBuffer<f_t> data(0, dev_sched);
256
257 patch_data_layers.get_const_refs().for_each(
258 [&](u64 id_patch,
259 const std::reference_wrapper<shamrock::patch::PatchDataLayer> &pdat) {
260 auto buf_field
261 = sham::DeviceBuffer<f_t>(pdat.get().get_obj_cnt() * block_size, dev_sched);
262 buf_field.fill(shamcomm::world_rank());
263 data.append(buf_field);
264 });
265 auto tmp_buf = data.copy_to_sycl_buffer();
266 writer.write_field("world_rank", tmp_buf, num_obj * block_size);
267 }
268 }
269
270 std::string _impl_get_label() { return "PatchDataLayerToVtk"; }
271
272 std::string _impl_get_tex() { return "TODO"; }
273};
274
275template<class Tvec, class TgridVec>
277
278 bool enable_mem_free = false;
279
280 auto get_optional_free_mem = [&](auto &bind_to, auto &add_to) {
281 if (enable_mem_free) {
283 node.set_edges(bind_to);
284 add_to.push_back(std::make_shared<decltype(node)>(std::move(node)));
285 }
286 };
287
288 {
289 storage.ghost_layout = std::make_shared<shamrock::patch::PatchDataLayerLayout>();
291 = shambase::get_check_ref(storage.ghost_layout);
292
293 ghost_layout.add_field<TgridVec>("cell_min", 1);
294 ghost_layout.add_field<TgridVec>("cell_max", 1);
295 ghost_layout.add_field<Tscal>("rho", AMRBlock::block_size);
296 ghost_layout.add_field<Tscal>("rhoetot", AMRBlock::block_size);
297 ghost_layout.add_field<Tvec>("rhovel", AMRBlock::block_size);
298
299 if (solver_config.is_dust_on()) {
300 auto ndust = solver_config.dust_config.ndust;
301 ghost_layout.add_field<Tscal>("rho_dust", ndust * AMRBlock::block_size);
302 ghost_layout.add_field<Tvec>("rhovel_dust", ndust * AMRBlock::block_size);
303 }
304
305 if (solver_config.is_gravity_on()) {
306 ghost_layout.add_field<Tscal>("phi", AMRBlock::block_size);
307 }
308
309 if (solver_config.is_gas_passive_scalar_on()) {
310 u32 npscal_gas = solver_config.npscal_gas_config.npscal_gas;
311 ghost_layout.add_field<Tscal>("rho_gas_pscal", npscal_gas * AMRBlock::block_size);
312 }
313 }
314
318
319 using namespace shamrock::solvergraph;
320
321 SolverGraph &graph = storage.solver_graph;
322
323 graph.register_edge("sptree", SerialPatchTreeRefEdge<TgridVec>("sptree", "sptree"));
324
325 graph.register_edge(
326 "global_patch_boxes",
327 ScalarsEdge<shammath::AABB<TgridVec>>("global_patch_boxes", "global_patch_boxes"));
328
329 graph.register_node(
330 "set_sptree",
332 edge.patch_tree = std::ref(storage.serial_patch_tree.get());
333 }));
334
336 .set_edges(graph.get_edge_ptr<SerialPatchTreeRefEdge<TgridVec>>("sptree"));
337
338 storage.local_patch_ids
339 = std::make_shared<shamrock::solvergraph::IDataEdge<std::vector<u64>>>("", "");
340
341 storage.sim_box_edge
342 = std::make_shared<shamrock::solvergraph::ScalarEdge<shammath::AABB<TgridVec>>>(
343 "sim_box", "sim_box");
344
345 storage.exchange_gz_edge = std::make_shared<shamrock::solvergraph::PatchDataLayerDDShared>(
346 "exchange_gz_edge", "exchange_gz_edge");
347
348 storage.idx_in_ghost = std::make_shared<shamrock::solvergraph::DDSharedBuffers<u32>>(
349 "idx_in_ghost", "idx_in_ghost");
350
351 storage.ghost_layers_candidates_edge = std::make_shared<
353 "ghost_layers_candidates", "ghost_layers_candidates");
354
355 storage.patch_rank_owner = std::make_shared<shamrock::solvergraph::RankGetter>(
356 [&](u64 patch_id) -> u32 {
357 return scheduler().get_patch_rank_owner(patch_id);
358 },
359 "patch_rank_owner",
360 "rank");
361
362 storage.source_patches = std::make_shared<shamrock::solvergraph::PatchDataLayerRefs>(
363 "source_patches", "P_{\\rm source}");
364
365 storage.merged_patchdata_ghost = std::make_shared<shamrock::solvergraph::PatchDataLayerEdge>(
366 "merged_patchdata_ghost", "patchdata_{\\rm ghost}", storage.ghost_layout);
367
368 storage.block_counts
369 = std::make_shared<shamrock::solvergraph::Indexes<u32>>("block_count", "N_{\\rm block}");
370
371 storage.block_counts_with_ghost = std::make_shared<shamrock::solvergraph::Indexes<u32>>(
372 "block_count_with_ghost", "N_{\\rm block, with ghost}");
373
374 // merged ghost spans
375 storage.refs_block_min = std::make_shared<shamrock::solvergraph::FieldRefs<TgridVec>>(
376 "block_min", "\\mathbf{r}_{\\rm block, min}");
377 storage.refs_block_max = std::make_shared<shamrock::solvergraph::FieldRefs<TgridVec>>(
378 "block_max", "\\mathbf{r}_{\\rm block, max}");
379
380 storage.refs_rho = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("rho", "\\rho");
381 storage.refs_rhov
382 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("rhovel", "(\\rho \\mathbf{v})");
383 storage.refs_rhoe
384 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("rhoetot", "(\\rho E)");
385
386 if (solver_config.is_dust_on()) {
387 storage.refs_rho_dust = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>(
388 "rho_dust", "\\rho_{\\rm dust}");
389 storage.refs_rhov_dust = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>(
390 "rhovel_dust", "(\\rho_{\\rm dust} \\mathbf{v}_{\\rm dust})");
391 }
392
393 // will be filled by NodeConsToPrimGas
394 storage.vel = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
395 AMRBlock::block_size, "vel", "\\mathbf{v}");
396 storage.press
397 = std::make_shared<shamrock::solvergraph::Field<Tscal>>(AMRBlock::block_size, "P", "P");
398
399 if (!solver_config.amr_mode.old_amr) { // TODO disable also if amr is none
400 storage.rho_primitive = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
401 AMRBlock::block_size, "rho-prim", "rho-prim");
402 }
403
404 // will be filled only if valid refinement criterion is provided
405 using AMRmode_None = typename AMRMode<Tvec, TgridVec>::None;
406 if (std::get_if<AMRmode_None>(&solver_config.amr_mode.config) == nullptr) {
407 storage.rho_snap = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
408 AMRBlock::block_size, "rho_snap", "rho_snap");
409 storage.rhoe_snap = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
410 AMRBlock::block_size, "rhoe_snap", "rhoe_snap");
411 storage.rho_vel_snap = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
412 AMRBlock::block_size, "rhov_snap", "rhov_snap");
413 }
414
415 if (solver_config.is_dust_on()) {
416 u32 ndust = solver_config.dust_config.ndust;
417
418 // will be filled by NodeConsToPrimDust
419 storage.vel_dust = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
420 AMRBlock::block_size * ndust, "vel_dust", "{\\mathbf{v}_{\\rm dust}}");
421 }
422
423 storage.trees
424 = std::make_shared<solvergraph::TreeEdge<u_morton, TgridVec>>("trees", "\\text{trees}");
425
426 storage.block_graph_edge = std::make_shared<
428 "block_graph_edge", "\\text{block graph edge}");
429
430 storage.cell_graph_edge = std::make_shared<
432 "cell_graph_edge", "\\text{cell graph edge}");
433
434 // will be filled by NodeComputeCellAABB
435 storage.block_cell_sizes = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
436 1, "block_cell_sizes", "s_{\\rm cell}");
437 storage.cell0block_aabb_lower = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
438 1, "cell0block_aabb_lower", "\\mathbf{s}_{\\rm inf,block}");
439
440 if (solver_config.is_coordinate_field_required()) {
441 // will be filled by NodeComputeCoordinates
442 storage.coordinates = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
443 AMRBlock::block_size, "coordinates", "\\mathbf{xyz}");
444 }
445
446 if (solver_config.amr_mode.need_level_zero_compute()) {
447 // get blocks at level0 sizes for all patches
448 storage.level0_size = std::make_shared<shamrock::solvergraph::ScalarsEdge<TgridVec>>(
449 "level0_amr", "level0_amr");
450 }
451
452 if (solver_config.amr_mode.need_amr_level_compute()) {
453 using TgridUint = typename std::make_unsigned<shambase::VecComponent<TgridVec>>::type;
454 storage.amr_block_levels
455 = std::make_shared<shamrock::solvergraph::Field<TgridUint>>(1, "", "");
456 }
457
458 storage.grad_rho = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
459 AMRBlock::block_size, "grad_rho", "\\nabla \\rho");
460 storage.dx_v = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
461 AMRBlock::block_size, "dx_v", "\\nabla_x \\mathbf{v}");
462 storage.dy_v = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
463 AMRBlock::block_size, "dy_v", "\\nabla_y \\mathbf{v}");
464 storage.dz_v = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
465 AMRBlock::block_size, "dz_v", "\\nabla_z \\mathbf{v}");
466 storage.grad_P = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
467 AMRBlock::block_size, "grad_P", "\\nabla P");
468
469 // will be filled by NodeEulerTimeDerivativeGas
470 storage.euler_dt_rho = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
471 AMRBlock::block_size, "euler_dt_rho", "\\partial_t \\rho");
472 storage.euler_dt_vel = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
473 AMRBlock::block_size, "euler_dt_vel", "\\partial_t \\mathbf{v}");
474 storage.euler_dt_press = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
475 AMRBlock::block_size, "euler_dt_press", "\\partial_t P");
476
477 if (solver_config.is_dust_on()) {
478 u32 ndust = solver_config.dust_config.ndust;
479 storage.grad_rho_dust = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
480 AMRBlock::block_size * ndust, "grad_rho_dust", "\\nabla \\rho_{\\rm dust}");
481 storage.dx_v_dust = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
482 AMRBlock::block_size * ndust, "dx_v_dust", "\\nabla_x \\mathbf{v}_{\\rm dust}");
483 storage.dy_v_dust = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
484 AMRBlock::block_size * ndust, "dy_v_dust", "\\nabla_y \\mathbf{v}_{\\rm dust}");
485 storage.dz_v_dust = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
486 AMRBlock::block_size * ndust, "dz_v_dust", "\\nabla_z \\mathbf{v}_{\\rm dust}");
487
488 // will be filled by NodeEulerTimeDerivativeDust
489 storage.euler_dt_rho_dust = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
490 AMRBlock::block_size * ndust, "euler_dt_rho_dust", "\\partial_t \\rho_{\\rm dust}");
491 storage.euler_dt_vel_dust = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
492 AMRBlock::block_size * ndust,
493 "euler_dt_vel_dust",
494 "\\partial_t \\mathbf{v}_{\\rm dust}");
495 }
496
497 {
498
499 storage.rho_face_xp
500 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
501 "rho_face_xp", "rho_face_xp", 1);
502 storage.rho_face_xm
503 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
504 "rho_face_xm", "rho_face_xm", 1);
505 storage.rho_face_yp
506 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
507 "rho_face_yp", "rho_face_yp", 1);
508 storage.rho_face_ym
509 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
510 "rho_face_ym", "rho_face_ym", 1);
511 storage.rho_face_zp
512 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
513 "rho_face_zp", "rho_face_zp", 1);
514 storage.rho_face_zm
515 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
516 "rho_face_zm", "rho_face_zm", 1);
517
518 storage.vel_face_xp
519 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
520 "vel_face_xp", "vel_face_xp", 1);
521 storage.vel_face_xm
522 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
523 "vel_face_xm", "vel_face_xm", 1);
524 storage.vel_face_yp
525 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
526 "vel_face_yp", "vel_face_yp", 1);
527 storage.vel_face_ym
528 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
529 "vel_face_ym", "vel_face_ym", 1);
530 storage.vel_face_zp
531 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
532 "vel_face_zp", "vel_face_zp", 1);
533 storage.vel_face_zm
534 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
535 "vel_face_zm", "vel_face_zm", 1);
536
537 storage.press_face_xp
538 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
539 "press_face_xp", "press_face_xp", 1);
540 storage.press_face_xm
541 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
542 "press_face_xm", "press_face_xm", 1);
543 storage.press_face_yp
544 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
545 "press_face_yp", "press_face_yp", 1);
546 storage.press_face_ym
547 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
548 "press_face_ym", "press_face_ym", 1);
549 storage.press_face_zp
550 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
551 "press_face_zp", "press_face_zp", 1);
552 storage.press_face_zm
553 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
554 "press_face_zm", "press_face_zm", 1);
555 }
556
557 if (solver_config.is_dust_on()) {
558 u32 ndust = solver_config.dust_config.ndust;
559
560 storage.rho_dust_face_xp
561 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
562 "rho_dust_face_xp", "rho_dust_face_xp", ndust);
563 storage.rho_dust_face_xm
564 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
565 "rho_dust_face_xm", "rho_dust_face_xm", ndust);
566 storage.rho_dust_face_yp
567 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
568 "rho_dust_face_yp", "rho_dust_face_yp", ndust);
569 storage.rho_dust_face_ym
570 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
571 "rho_dust_face_ym", "rho_dust_face_ym", ndust);
572 storage.rho_dust_face_zp
573 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
574 "rho_dust_face_zp", "rho_dust_face_zp", ndust);
575 storage.rho_dust_face_zm
576 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>>>(
577 "rho_dust_face_zm", "rho_dust_face_zm", ndust);
578
579 storage.vel_dust_face_xp
580 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
581 "vel_dust_face_xp", "vel_dust_face_xp", ndust);
582 storage.vel_dust_face_xm
583 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
584 "vel_dust_face_xm", "vel_dust_face_xm", ndust);
585 storage.vel_dust_face_yp
586 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
587 "vel_dust_face_yp", "vel_dust_face_yp", ndust);
588 storage.vel_dust_face_ym
589 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
590 "vel_dust_face_ym", "vel_dust_face_ym", ndust);
591 storage.vel_dust_face_zp
592 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
593 "vel_dust_face_zp", "vel_dust_face_zp", ndust);
594 storage.vel_dust_face_zm
595 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>>>(
596 "vel_dust_face_zm", "vel_dust_face_zm", ndust);
597 }
598
599 if (solver_config.should_compute_rho_mean()) {
600 storage.cell_mass = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
601 AMRBlock::block_size, "cell_mass", "m");
602 storage.rho_mean
603 = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>("rho_mean", "< \\rho >");
604 storage.simulation_volume = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>(
605 "simulation_volume", "V_{\\rm sim}");
606 }
607
608 storage.dt_over2
609 = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>("dt_half", "dt_{half}");
610
611 {
612
613 storage.flux_rho_face_xm = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
614 "flux_rho_face_xm", "flux_rho_face_xm", 1);
615
616 storage.flux_rho_face_xp = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
617 "flux_rho_face_xp", "flux_rho_face_xp", 1);
618
619 storage.flux_rho_face_ym = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
620 "flux_rho_face_ym", "flux_rho_face_ym", 1);
621
622 storage.flux_rho_face_yp = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
623 "flux_rho_face_yp", "flux_rho_face_yp", 1);
624
625 storage.flux_rho_face_zm = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
626 "flux_rho_face_zm", "flux_rho_face_zm", 1);
627
628 storage.flux_rho_face_zp = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
629 "flux_rho_face_zp", "flux_rho_face_zp", 1);
630
631 storage.flux_rhov_face_xm = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
632 "flux_rhov_face_xm", "flux_rhov_face_xm", 1);
633
634 storage.flux_rhov_face_xp = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
635 "flux_rhov_face_xp", "flux_rhov_face_xp", 1);
636
637 storage.flux_rhov_face_ym = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
638 "flux_rhov_face_ym", "flux_rhov_face_ym", 1);
639
640 storage.flux_rhov_face_yp = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
641 "flux_rhov_face_yp", "flux_rhov_face_yp", 1);
642
643 storage.flux_rhov_face_zm = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
644 "flux_rhov_face_zm", "flux_rhov_face_zm", 1);
645
646 storage.flux_rhov_face_zp = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
647 "flux_rhov_face_zp", "flux_rhov_face_zp", 1);
648
649 storage.flux_rhoe_face_xm = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
650 "flux_rhoe_face_xm", "flux_rhoe_face_xm", 1);
651
652 storage.flux_rhoe_face_xp = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
653 "flux_rhoe_face_xp", "flux_rhoe_face_xp", 1);
654
655 storage.flux_rhoe_face_ym = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
656 "flux_rhoe_face_ym", "flux_rhoe_face_ym", 1);
657
658 storage.flux_rhoe_face_yp = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
659 "flux_rhoe_face_yp", "flux_rhoe_face_yp", 1);
660
661 storage.flux_rhoe_face_zm = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
662 "flux_rhoe_face_zm", "flux_rhoe_face_zm", 1);
663
664 storage.flux_rhoe_face_zp = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
665 "flux_rhoe_face_zp", "flux_rhoe_face_zp", 1);
666 }
667
668 if (solver_config.is_dust_on()) {
669 u32 ndust = solver_config.dust_config.ndust;
670
671 storage.flux_rho_dust_face_xm
672 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
673 "flux_rho_dust_face_xm", "flux_rho_dust_face_xm", ndust);
674
675 storage.flux_rho_dust_face_xp
676 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
677 "flux_rho_dust_face_xp", "flux_rho_dust_face_xp", ndust);
678
679 storage.flux_rho_dust_face_ym
680 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
681 "flux_rho_dust_face_ym", "flux_rho_dust_face_ym", ndust);
682
683 storage.flux_rho_dust_face_yp
684 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
685 "flux_rho_dust_face_yp", "flux_rho_dust_face_yp", ndust);
686
687 storage.flux_rho_dust_face_zm
688 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
689 "flux_rho_dust_face_zm", "flux_rho_dust_face_zm", ndust);
690
691 storage.flux_rho_dust_face_zp
692 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tscal>>(
693 "flux_rho_dust_face_zp", "flux_rho_dust_face_zp", ndust);
694
695 storage.flux_rhov_dust_face_xm
696 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
697 "flux_rhov_dust_face_xm", "flux_rhov_dust_face_xm", ndust);
698
699 storage.flux_rhov_dust_face_xp
700 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
701 "flux_rhov_dust_face_xp", "flux_rhov_dust_face_xp", ndust);
702
703 storage.flux_rhov_dust_face_ym
704 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
705 "flux_rhov_dust_face_ym", "flux_rhov_dust_face_ym", ndust);
706
707 storage.flux_rhov_dust_face_yp
708 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
709 "flux_rhov_dust_face_yp", "flux_rhov_dust_face_yp", ndust);
710
711 storage.flux_rhov_dust_face_zm
712 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
713 "flux_rhov_dust_face_zm", "flux_rhov_dust_face_zm", ndust);
714
715 storage.flux_rhov_dust_face_zp
716 = std::make_shared<solvergraph::NeighGraphLinkFieldEdge<Tvec>>(
717 "flux_rhov_dust_face_zp", "flux_rhov_dust_face_zp", ndust);
718 }
719
720 storage.dtrho = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
721 AMRBlock::block_size, "dtrho", "dtrho");
722 storage.dtrhov = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
723 AMRBlock::block_size, "dtrhov", "dtrhov");
724 storage.dtrhoe = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
725 AMRBlock::block_size, "dtrhoe", "dtrhoe");
726
727 if (solver_config.is_dust_on()) {
728 u32 ndust = solver_config.dust_config.ndust;
729 storage.dtrho_dust = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
730 AMRBlock::block_size * ndust, "dtrho_dust", "dtrho_dust");
731 storage.dtrhov_dust = std::make_shared<shamrock::solvergraph::Field<Tvec>>(
732 AMRBlock::block_size * ndust, "dtrhov_dust", "dtrhov_dust");
733 }
734
738 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> solver_sequence;
739
740 solver_sequence.push_back(graph.get_node_ptr_base("set_sptree"));
741
742 auto cfg_bc_to_geom = [](BCConfig::GhostType ghost_type) {
743 switch (ghost_type) {
744 case BCConfig::GhostType::Periodic : return modules::GhostType::Periodic;
745 case BCConfig::GhostType::Reflective: return modules::GhostType::Reflective;
746 case BCConfig::GhostType::Outflow : return modules::GhostType::Reflective;
747 default:
749 "Unsupported ghost type: " + std::to_string(static_cast<int>(ghost_type)));
750 }
751 };
752
753 modules::GhostLayerGenMode ghost_layer_gen_mode{
754 cfg_bc_to_geom(solver_config.bc_config.get_x()),
755 cfg_bc_to_geom(solver_config.bc_config.get_y()),
756 cfg_bc_to_geom(solver_config.bc_config.get_z())};
757
758 // if outflow we want zero gradient so we skip the vector transformation in TransformGhostLayer
759 bool transform_vec_x = solver_config.bc_config.get_x() != BCConfig::GhostType::Outflow;
760 bool transform_vec_y = solver_config.bc_config.get_y() != BCConfig::GhostType::Outflow;
761 bool transform_vec_z = solver_config.bc_config.get_z() != BCConfig::GhostType::Outflow;
762
763 { // Ghost zone finder
764
765 modules::FindGhostLayerCandidates<TgridVec> find_ghost_layer_candidates(
766 ghost_layer_gen_mode);
767 find_ghost_layer_candidates.set_edges(
768 storage.local_patch_ids,
769 storage.sim_box_edge,
771 graph.get_edge_ptr<ScalarsEdge<shammath::AABB<TgridVec>>>("global_patch_boxes"),
772 storage.ghost_layers_candidates_edge);
773 solver_sequence.push_back(
774 std::make_shared<decltype(find_ghost_layer_candidates)>(
775 std::move(find_ghost_layer_candidates)));
776
777 modules::FindGhostLayerIndices<TgridVec> find_ghost_layer_indices(ghost_layer_gen_mode);
778 find_ghost_layer_indices.set_edges(
779 storage.sim_box_edge,
780 storage.source_patches,
781 storage.ghost_layers_candidates_edge,
782 graph.get_edge_ptr<ScalarsEdge<shammath::AABB<TgridVec>>>("global_patch_boxes"),
783 storage.idx_in_ghost);
784 solver_sequence.push_back(
785 std::make_shared<decltype(find_ghost_layer_indices)>(
786 std::move(find_ghost_layer_indices)));
787 }
788
789 { // Ghost zone exchange
790 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> gz_xchg_sequence;
791
792 auto &ghost_layout_ptr = storage.ghost_layout;
793 {
794 auto copy_fields = std::make_shared<shamrock::solvergraph::CopyPatchDataLayerFields>(
795 scheduler().get_layout_ptr_old(), ghost_layout_ptr);
796
797 copy_fields->set_edges(storage.source_patches, storage.merged_patchdata_ghost);
798 gz_xchg_sequence.push_back(std::move(copy_fields));
799 }
800
801 {
802 auto extract_gz_node
803 = std::make_shared<shammodels::basegodunov::modules::ExtractGhostLayer>(
804 ghost_layout_ptr);
805
806 extract_gz_node->set_edges(
807 storage.merged_patchdata_ghost, storage.idx_in_ghost, storage.exchange_gz_edge);
808 gz_xchg_sequence.push_back(std::move(extract_gz_node));
809 }
810
811 {
812 auto transform_gz_node = std::make_shared<
814 ghost_layer_gen_mode,
815 transform_vec_x,
816 transform_vec_y,
817 transform_vec_z,
818 ghost_layout_ptr);
819
820 transform_gz_node->set_edges(
821 storage.sim_box_edge,
822 storage.ghost_layers_candidates_edge,
823 storage.exchange_gz_edge);
824 gz_xchg_sequence.push_back(std::move(transform_gz_node));
825 }
826
827 {
828 auto exchange_gz_node
829 = std::make_shared<shamrock::solvergraph::ExchangeGhostLayer>(ghost_layout_ptr);
830 exchange_gz_node->set_edges(storage.patch_rank_owner, storage.exchange_gz_edge);
831 gz_xchg_sequence.push_back(std::move(exchange_gz_node));
832 }
833
834 {
835 auto fuse_gz_node
836 = std::make_shared<shammodels::basegodunov::modules::FuseGhostLayer>();
837 fuse_gz_node->set_edges(storage.exchange_gz_edge, storage.merged_patchdata_ghost);
838 gz_xchg_sequence.push_back(std::move(fuse_gz_node));
839 }
840
841 // enable this to debug GZ
842 //{
843 // auto filename_edge = std::make_shared<shamrock::solvergraph::IDataEdge<std::string>>(
844 // "debug_fuse.vtk", "debug_fuse.vtk");
845 // filename_edge->data = "debug_fuse.vtk";
846 //
847 // auto patch_data_layer_to_vtk_node
848 // = std::make_shared<PatchDataLayerToVtk<Tvec, TgridVec>>(true, true, 8);
849 // patch_data_layer_to_vtk_node->set_edges(filename_edge,
850 // storage.merged_patchdata_ghost);
851 // gz_xchg_sequence.push_back(std::move(patch_data_layer_to_vtk_node));
852 //}
853
855 "Ghost zone exchange", std::move(gz_xchg_sequence));
856 solver_sequence.push_back(std::make_shared<decltype(seq)>(std::move(seq)));
857 }
858
859 { // attach fields with ghosts
860
861 { // set element counts
862 shamrock::solvergraph::ExtractCounts extract_counts_node
864 extract_counts_node.set_edges(storage.source_patches, storage.block_counts);
865 solver_sequence.push_back(
866 std::make_shared<decltype(extract_counts_node)>(std::move(extract_counts_node)));
867 }
868
869 { // set element counts
870 shamrock::solvergraph::ExtractCounts extract_counts_node
872 extract_counts_node.set_edges(
873 storage.merged_patchdata_ghost, storage.block_counts_with_ghost);
874 solver_sequence.push_back(
875 std::make_shared<decltype(extract_counts_node)>(std::move(extract_counts_node)));
876 }
877
878 { // Attach spans to block coords
881 attach_block_min.set_edges(storage.merged_patchdata_ghost, storage.refs_block_min);
882 solver_sequence.push_back(
883 std::make_shared<decltype(attach_block_min)>(std::move(attach_block_min)));
884
887 attach_block_max.set_edges(storage.merged_patchdata_ghost, storage.refs_block_max);
888 solver_sequence.push_back(
889 std::make_shared<decltype(attach_block_max)>(std::move(attach_block_max)));
890 }
891
892 { // attach spans to gas field with ghosts
894 = shamrock::solvergraph::GetFieldRefFromLayer<Tscal>(storage.ghost_layout, "rho");
895 attach_rho.set_edges(storage.merged_patchdata_ghost, storage.refs_rho);
896 solver_sequence.push_back(
897 std::make_shared<decltype(attach_rho)>(std::move(attach_rho)));
898
900 = shamrock::solvergraph::GetFieldRefFromLayer<Tvec>(storage.ghost_layout, "rhovel");
901 attach_rhov.set_edges(storage.merged_patchdata_ghost, storage.refs_rhov);
902 solver_sequence.push_back(
903 std::make_shared<decltype(attach_rhov)>(std::move(attach_rhov)));
904
907 storage.ghost_layout, "rhoetot");
908 attach_rhoe.set_edges(storage.merged_patchdata_ghost, storage.refs_rhoe);
909 solver_sequence.push_back(
910 std::make_shared<decltype(attach_rhoe)>(std::move(attach_rhoe)));
911 }
912
913 if (solver_config.is_dust_on()) { // attach spans to dust field with ghosts
916 storage.ghost_layout, "rho_dust");
917 attach_rho_dust.set_edges(storage.merged_patchdata_ghost, storage.refs_rho_dust);
918 solver_sequence.push_back(
919 std::make_shared<decltype(attach_rho_dust)>(std::move(attach_rho_dust)));
920
923 storage.ghost_layout, "rhovel_dust");
924 attach_rhov_dust.set_edges(storage.merged_patchdata_ghost, storage.refs_rhov_dust);
925 solver_sequence.push_back(
926 std::make_shared<decltype(attach_rhov_dust)>(std::move(attach_rhov_dust)));
927 }
928 }
929
930 { // build trees
931
933 node.set_edges(
934 storage.block_counts_with_ghost,
935 storage.refs_block_min,
936 storage.refs_block_max,
937 storage.trees);
938
939 solver_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
940 }
941
942 { // build neigh tables
943 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> neigh_table_sequence;
944
946 node1.set_edges(
947 storage.block_counts_with_ghost,
948 storage.refs_block_min,
949 storage.refs_block_max,
950 storage.trees,
951 storage.block_graph_edge);
952 node1.evaluate();
953
955 node2.set_edges(
956 storage.block_counts_with_ghost,
957 storage.refs_block_min,
958 storage.refs_block_max,
959 storage.block_graph_edge,
960 storage.cell_graph_edge);
961 node2.evaluate();
962
963 neigh_table_sequence.push_back(std::make_shared<decltype(node1)>(std::move(node1)));
964 get_optional_free_mem(storage.trees, neigh_table_sequence);
965 neigh_table_sequence.push_back(std::make_shared<decltype(node2)>(std::move(node2)));
966 get_optional_free_mem(storage.block_graph_edge, neigh_table_sequence);
967
969 "Compute neigh table", std::move(neigh_table_sequence));
970 solver_sequence.push_back(std::make_shared<decltype(seq)>(std::move(seq)));
971 }
972
973 { // Compute cell infos
974
976 AMRBlock::Nside, solver_config.grid_coord_to_pos_fact};
977
978 node.set_edges(
979 storage.block_counts_with_ghost,
980 storage.refs_block_min,
981 storage.refs_block_max,
982 storage.block_cell_sizes,
983 storage.cell0block_aabb_lower);
984 solver_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
985 }
986
987 if (solver_config.is_coordinate_field_required()) { // Compute coordinates
988
990 AMRBlock::block_size,
991 AMRBlock::Nside,
992 solver_config.grid_coord_to_pos_fact,
993 };
994
995 node_coordinates.set_edges(
996 storage.block_counts,
997 storage.refs_block_min,
998 storage.refs_block_max,
999 storage.coordinates);
1000
1001 solver_sequence.push_back(
1002 std::make_shared<decltype(node_coordinates)>(std::move(node_coordinates)));
1003 }
1004
1005 if (solver_config.amr_mode.need_level_zero_compute()) { // compute level0 sizes in patch (to be
1006 // enabled later when needed)
1007 modules::ComputeLevel0CellSize<TgridVec> node_level0_sizes{};
1008 node_level0_sizes.set_edges(
1009 graph.get_edge_ptr<ScalarsEdge<shammath::AABB<TgridVec>>>("global_patch_boxes"),
1010 storage.source_patches,
1011 storage.level0_size);
1012 solver_sequence.push_back(
1013 std::make_shared<decltype(node_level0_sizes)>(std::move(node_level0_sizes)));
1014 }
1015
1016 if (solver_config.amr_mode.need_amr_level_compute()) { // compute block amr level in patch
1017 modules::ComputeAMRLevel<TgridVec> node_amr_level{};
1018 node_amr_level.set_edges(
1019 storage.block_counts,
1020 storage.level0_size,
1021 storage.refs_block_min,
1022 storage.refs_block_max,
1023 storage.amr_block_levels);
1024 solver_sequence.push_back(
1025 std::make_shared<decltype(node_amr_level)>(std::move(node_amr_level)));
1026 }
1027
1028 if (solver_config.should_compute_rho_mean()) {
1029 modules::NodeComputeMass<Tvec, TgridVec> node{AMRBlock::block_size};
1030 node.set_edges(
1031 storage.block_counts, storage.block_cell_sizes, storage.refs_rho, storage.cell_mass);
1032 solver_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1033
1034 modules::NodeComputeSumOverV<Tscal> node2{AMRBlock::block_size};
1035 node2.set_edges(
1036 storage.block_counts, storage.cell_mass, storage.simulation_volume, storage.rho_mean);
1037 solver_sequence.push_back(std::make_shared<decltype(node2)>(std::move(node2)));
1038 }
1039
1040 { // Build ConsToPrim node
1041 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> const_to_prim_sequence;
1042
1043 {
1044 modules::NodeConsToPrimGas<Tvec> node{AMRBlock::block_size, solver_config.eos_gamma};
1045 node.set_edges(
1046 storage.block_counts_with_ghost,
1047 storage.refs_rho,
1048 storage.refs_rhov,
1049 storage.refs_rhoe,
1050 storage.vel,
1051 storage.press);
1052
1053 const_to_prim_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1054 }
1055
1056 if (solver_config.is_dust_on()) {
1057 u32 ndust = solver_config.dust_config.ndust;
1058 modules::NodeConsToPrimDust<Tvec> node{AMRBlock::block_size, ndust};
1059 node.set_edges(
1060 storage.block_counts_with_ghost,
1061 storage.refs_rho_dust,
1062 storage.refs_rhov_dust,
1063 storage.vel_dust);
1064
1065 const_to_prim_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1066 }
1067
1069 "Cons to Prim", std::move(const_to_prim_sequence));
1070 solver_sequence.push_back(std::make_shared<decltype(seq)>(std::move(seq)));
1071 }
1072
1073 { // Build slope limited gradients
1074
1075 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> grad_sequence;
1076
1077 {
1079 AMRBlock::block_size, 1, solver_config.slope_config};
1080 node.set_edges(
1081 storage.block_counts_with_ghost,
1082 storage.cell_graph_edge,
1083 storage.block_cell_sizes,
1084 storage.refs_rho,
1085 storage.grad_rho);
1086 grad_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1087 }
1088
1089 {
1091 AMRBlock::block_size, 1, solver_config.slope_config};
1092 node.set_edges(
1093 storage.block_counts_with_ghost,
1094 storage.cell_graph_edge,
1095 storage.block_cell_sizes,
1096 storage.vel,
1097 storage.dx_v,
1098 storage.dy_v,
1099 storage.dz_v);
1100 grad_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1101 }
1102 {
1104 AMRBlock::block_size, 1, solver_config.slope_config};
1105 node.set_edges(
1106 storage.block_counts_with_ghost,
1107 storage.cell_graph_edge,
1108 storage.block_cell_sizes,
1109 storage.press,
1110 storage.grad_P);
1111 grad_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1112 }
1113
1114 if (solver_config.is_dust_on()) {
1115 u32 ndust = solver_config.dust_config.ndust;
1117 AMRBlock::block_size, ndust, solver_config.slope_config};
1118 node.set_edges(
1119 storage.block_counts_with_ghost,
1120 storage.cell_graph_edge,
1121 storage.block_cell_sizes,
1122 storage.refs_rho_dust,
1123 storage.grad_rho_dust);
1124 grad_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1125
1127 AMRBlock::block_size, ndust, solver_config.slope_config};
1128 node2.set_edges(
1129 storage.block_counts_with_ghost,
1130 storage.cell_graph_edge,
1131 storage.block_cell_sizes,
1132 storage.vel_dust,
1133 storage.dx_v_dust,
1134 storage.dy_v_dust,
1135 storage.dz_v_dust);
1136 grad_sequence.push_back(std::make_shared<decltype(node2)>(std::move(node2)));
1137 }
1138
1140 "Slope limited gradients", std::move(grad_sequence));
1141 solver_sequence.push_back(std::make_shared<decltype(seq)>(std::move(seq)));
1142 }
1143
1144 { // Euler time derivatives of the primitive state
1145 // Hoisted out of the face interpolation nodes: they only depend on cell
1146 // local quantities, so computing them once per cell here avoids
1147 // re-fetching the velocity gradients for every face link.
1148 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> dt_prim_sequence;
1149
1150 {
1152 AMRBlock::block_size, solver_config.eos_gamma};
1153 node.set_edges(
1154 storage.block_counts_with_ghost,
1155 storage.refs_rho,
1156 storage.vel,
1157 storage.press,
1158 storage.grad_rho,
1159 storage.dx_v,
1160 storage.dy_v,
1161 storage.dz_v,
1162 storage.grad_P,
1163 storage.euler_dt_rho,
1164 storage.euler_dt_vel,
1165 storage.euler_dt_press);
1166 dt_prim_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1167 }
1168
1169 if (solver_config.is_dust_on()) {
1170 u32 ndust = solver_config.dust_config.ndust;
1171 modules::NodeEulerTimeDerivativeDust<Tvec> node{AMRBlock::block_size, ndust};
1172 node.set_edges(
1173 storage.block_counts_with_ghost,
1174 storage.refs_rho_dust,
1175 storage.vel_dust,
1176 storage.grad_rho_dust,
1177 storage.dx_v_dust,
1178 storage.dy_v_dust,
1179 storage.dz_v_dust,
1180 storage.euler_dt_rho_dust,
1181 storage.euler_dt_vel_dust);
1182 dt_prim_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1183 }
1184
1186 "Euler time derivatives", std::move(dt_prim_sequence));
1187 solver_sequence.push_back(std::make_shared<decltype(seq)>(std::move(seq)));
1188 }
1189
1190 { // interpolate to face
1191 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> interp_sequence;
1192 {
1193 modules::InterpolateToFaceRho<Tvec, TgridVec> node{AMRBlock::block_size};
1194 node.set_edges(
1195 storage.dt_over2,
1196 storage.cell_graph_edge,
1197 storage.block_cell_sizes,
1198 storage.cell0block_aabb_lower,
1199 storage.refs_rho,
1200 storage.grad_rho,
1201 storage.euler_dt_rho,
1202 storage.rho_face_xp,
1203 storage.rho_face_xm,
1204 storage.rho_face_yp,
1205 storage.rho_face_ym,
1206 storage.rho_face_zp,
1207 storage.rho_face_zm);
1208 interp_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1209 }
1210
1211 {
1212 modules::InterpolateToFaceVel<Tvec, TgridVec> node{AMRBlock::block_size};
1213 node.set_edges(
1214 storage.dt_over2,
1215 storage.cell_graph_edge,
1216 storage.block_cell_sizes,
1217 storage.cell0block_aabb_lower,
1218 storage.vel,
1219 storage.dx_v,
1220 storage.dy_v,
1221 storage.dz_v,
1222 storage.euler_dt_vel,
1223 storage.vel_face_xp,
1224 storage.vel_face_xm,
1225 storage.vel_face_yp,
1226 storage.vel_face_ym,
1227 storage.vel_face_zp,
1228 storage.vel_face_zm);
1229 interp_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1230 }
1231
1232 {
1233 modules::InterpolateToFacePress<Tvec, TgridVec> node{AMRBlock::block_size};
1234 node.set_edges(
1235 storage.dt_over2,
1236 storage.cell_graph_edge,
1237 storage.block_cell_sizes,
1238 storage.cell0block_aabb_lower,
1239 storage.press,
1240 storage.grad_P,
1241 storage.euler_dt_press,
1242 storage.press_face_xp,
1243 storage.press_face_xm,
1244 storage.press_face_yp,
1245 storage.press_face_ym,
1246 storage.press_face_zp,
1247 storage.press_face_zm);
1248 interp_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1249 }
1250
1251 if (solver_config.is_dust_on()) {
1252 u32 ndust = solver_config.dust_config.ndust;
1253 modules::InterpolateToFaceRhoDust<Tvec, TgridVec> node{AMRBlock::block_size, ndust};
1254 node.set_edges(
1255 storage.dt_over2,
1256 storage.cell_graph_edge,
1257 storage.block_cell_sizes,
1258 storage.cell0block_aabb_lower,
1259 storage.refs_rho_dust,
1260 storage.grad_rho_dust,
1261 storage.euler_dt_rho_dust,
1262 storage.rho_dust_face_xp,
1263 storage.rho_dust_face_xm,
1264 storage.rho_dust_face_yp,
1265 storage.rho_dust_face_ym,
1266 storage.rho_dust_face_zp,
1267 storage.rho_dust_face_zm);
1268 interp_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1269 }
1270
1271 if (solver_config.is_dust_on()) {
1272 u32 ndust = solver_config.dust_config.ndust;
1273 modules::InterpolateToFaceVelDust<Tvec, TgridVec> node{AMRBlock::block_size, ndust};
1274 node.set_edges(
1275 storage.dt_over2,
1276 storage.cell_graph_edge,
1277 storage.block_cell_sizes,
1278 storage.cell0block_aabb_lower,
1279 storage.vel_dust,
1280 storage.dx_v_dust,
1281 storage.dy_v_dust,
1282 storage.dz_v_dust,
1283 storage.euler_dt_vel_dust,
1284 storage.vel_dust_face_xp,
1285 storage.vel_dust_face_xm,
1286 storage.vel_dust_face_yp,
1287 storage.vel_dust_face_ym,
1288 storage.vel_dust_face_zp,
1289 storage.vel_dust_face_zm);
1290 interp_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1291 }
1292
1294 "Interpolate to face", std::move(interp_sequence));
1295 solver_sequence.push_back(std::make_shared<decltype(seq)>(std::move(seq)));
1296 }
1297
1298 { // flux
1299
1300 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> flux_sequence;
1301
1302 if (solver_config.riemann_config == Rusanov) {
1304 node(
1305 "Gas flux compute",
1306 solver_config.eos_gamma,
1307 storage.cell_graph_edge,
1308 storage.rho_face_xp,
1309 storage.rho_face_xm,
1310 storage.rho_face_yp,
1311 storage.rho_face_ym,
1312 storage.rho_face_zp,
1313 storage.rho_face_zm,
1314 storage.vel_face_xp,
1315 storage.vel_face_xm,
1316 storage.vel_face_yp,
1317 storage.vel_face_ym,
1318 storage.vel_face_zp,
1319 storage.vel_face_zm,
1320 storage.press_face_xp,
1321 storage.press_face_xm,
1322 storage.press_face_yp,
1323 storage.press_face_ym,
1324 storage.press_face_zp,
1325 storage.press_face_zm,
1326 storage.flux_rho_face_xp,
1327 storage.flux_rho_face_xm,
1328 storage.flux_rho_face_yp,
1329 storage.flux_rho_face_ym,
1330 storage.flux_rho_face_zp,
1331 storage.flux_rho_face_zm,
1332 storage.flux_rhov_face_xp,
1333 storage.flux_rhov_face_xm,
1334 storage.flux_rhov_face_yp,
1335 storage.flux_rhov_face_ym,
1336 storage.flux_rhov_face_zp,
1337 storage.flux_rhov_face_zm,
1338 storage.flux_rhoe_face_xp,
1339 storage.flux_rhoe_face_xm,
1340 storage.flux_rhoe_face_yp,
1341 storage.flux_rhoe_face_ym,
1342 storage.flux_rhoe_face_zp,
1343 storage.flux_rhoe_face_zm);
1344 flux_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1345 } else if (solver_config.riemann_config == HLL) {
1347 "Gas flux compute",
1348 solver_config.eos_gamma,
1349 storage.cell_graph_edge,
1350 storage.rho_face_xp,
1351 storage.rho_face_xm,
1352 storage.rho_face_yp,
1353 storage.rho_face_ym,
1354 storage.rho_face_zp,
1355 storage.rho_face_zm,
1356 storage.vel_face_xp,
1357 storage.vel_face_xm,
1358 storage.vel_face_yp,
1359 storage.vel_face_ym,
1360 storage.vel_face_zp,
1361 storage.vel_face_zm,
1362 storage.press_face_xp,
1363 storage.press_face_xm,
1364 storage.press_face_yp,
1365 storage.press_face_ym,
1366 storage.press_face_zp,
1367 storage.press_face_zm,
1368 storage.flux_rho_face_xp,
1369 storage.flux_rho_face_xm,
1370 storage.flux_rho_face_yp,
1371 storage.flux_rho_face_ym,
1372 storage.flux_rho_face_zp,
1373 storage.flux_rho_face_zm,
1374 storage.flux_rhov_face_xp,
1375 storage.flux_rhov_face_xm,
1376 storage.flux_rhov_face_yp,
1377 storage.flux_rhov_face_ym,
1378 storage.flux_rhov_face_zp,
1379 storage.flux_rhov_face_zm,
1380 storage.flux_rhoe_face_xp,
1381 storage.flux_rhoe_face_xm,
1382 storage.flux_rhoe_face_yp,
1383 storage.flux_rhoe_face_ym,
1384 storage.flux_rhoe_face_zp,
1385 storage.flux_rhoe_face_zm);
1386 flux_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1387 } else if (solver_config.riemann_config == HLLC) {
1389 "Gas flux compute",
1390 solver_config.eos_gamma,
1391 storage.cell_graph_edge,
1392 storage.rho_face_xp,
1393 storage.rho_face_xm,
1394 storage.rho_face_yp,
1395 storage.rho_face_ym,
1396 storage.rho_face_zp,
1397 storage.rho_face_zm,
1398 storage.vel_face_xp,
1399 storage.vel_face_xm,
1400 storage.vel_face_yp,
1401 storage.vel_face_ym,
1402 storage.vel_face_zp,
1403 storage.vel_face_zm,
1404 storage.press_face_xp,
1405 storage.press_face_xm,
1406 storage.press_face_yp,
1407 storage.press_face_ym,
1408 storage.press_face_zp,
1409 storage.press_face_zm,
1410 storage.flux_rho_face_xp,
1411 storage.flux_rho_face_xm,
1412 storage.flux_rho_face_yp,
1413 storage.flux_rho_face_ym,
1414 storage.flux_rho_face_zp,
1415 storage.flux_rho_face_zm,
1416 storage.flux_rhov_face_xp,
1417 storage.flux_rhov_face_xm,
1418 storage.flux_rhov_face_yp,
1419 storage.flux_rhov_face_ym,
1420 storage.flux_rhov_face_zp,
1421 storage.flux_rhov_face_zm,
1422 storage.flux_rhoe_face_xp,
1423 storage.flux_rhoe_face_xm,
1424 storage.flux_rhoe_face_yp,
1425 storage.flux_rhoe_face_ym,
1426 storage.flux_rhoe_face_zp,
1427 storage.flux_rhoe_face_zm);
1428 flux_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1429 } else {
1430 shambase::throw_unimplemented("unknown flux mode");
1431 }
1432
1433 if (solver_config.is_dust_on()) {
1434 u32 ndust = solver_config.dust_config.ndust;
1435 if (solver_config.dust_config.dust_riemann_config == DHLL) {
1436 modules::
1437 NodeComputeFluxDustMode<Tvec, TgridVec, modules::DustRiemannSolverMode::DHLL>
1438 node(
1439 "Dust flux compute",
1440 ndust,
1441 storage.cell_graph_edge,
1442 storage.rho_dust_face_xp,
1443 storage.rho_dust_face_xm,
1444 storage.rho_dust_face_yp,
1445 storage.rho_dust_face_ym,
1446 storage.rho_dust_face_zp,
1447 storage.rho_dust_face_zm,
1448 storage.vel_dust_face_xp,
1449 storage.vel_dust_face_xm,
1450 storage.vel_dust_face_yp,
1451 storage.vel_dust_face_ym,
1452 storage.vel_dust_face_zp,
1453 storage.vel_dust_face_zm,
1454 storage.flux_rho_dust_face_xp,
1455 storage.flux_rho_dust_face_xm,
1456 storage.flux_rho_dust_face_yp,
1457 storage.flux_rho_dust_face_ym,
1458 storage.flux_rho_dust_face_zp,
1459 storage.flux_rho_dust_face_zm,
1460 storage.flux_rhov_dust_face_xp,
1461 storage.flux_rhov_dust_face_xm,
1462 storage.flux_rhov_dust_face_yp,
1463 storage.flux_rhov_dust_face_ym,
1464 storage.flux_rhov_dust_face_zp,
1465 storage.flux_rhov_dust_face_zm);
1466 flux_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1467 } else if (solver_config.dust_config.dust_riemann_config == HB) {
1469 node(
1470 "Dust flux compute",
1471 ndust,
1472 storage.cell_graph_edge,
1473 storage.rho_dust_face_xp,
1474 storage.rho_dust_face_xm,
1475 storage.rho_dust_face_yp,
1476 storage.rho_dust_face_ym,
1477 storage.rho_dust_face_zp,
1478 storage.rho_dust_face_zm,
1479 storage.vel_dust_face_xp,
1480 storage.vel_dust_face_xm,
1481 storage.vel_dust_face_yp,
1482 storage.vel_dust_face_ym,
1483 storage.vel_dust_face_zp,
1484 storage.vel_dust_face_zm,
1485 storage.flux_rho_dust_face_xp,
1486 storage.flux_rho_dust_face_xm,
1487 storage.flux_rho_dust_face_yp,
1488 storage.flux_rho_dust_face_ym,
1489 storage.flux_rho_dust_face_zp,
1490 storage.flux_rho_dust_face_zm,
1491 storage.flux_rhov_dust_face_xp,
1492 storage.flux_rhov_dust_face_xm,
1493 storage.flux_rhov_dust_face_yp,
1494 storage.flux_rhov_dust_face_ym,
1495 storage.flux_rhov_dust_face_zp,
1496 storage.flux_rhov_dust_face_zm);
1497 flux_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1498 } else {
1499 shambase::throw_unimplemented("unknown flux mode");
1500 }
1501 }
1502
1503 shamrock::solvergraph::OperationSequence seq("Compute fluxes", std::move(flux_sequence));
1504 solver_sequence.push_back(std::make_shared<decltype(seq)>(std::move(seq)));
1505 }
1506
1507 {
1509 AMRBlock::block_size, solver_config.grid_coord_to_pos_fact};
1510 node.set_edges(
1511 storage.block_counts,
1512 storage.cell_graph_edge,
1513 storage.block_cell_sizes,
1514 storage.cell0block_aabb_lower,
1515 storage.flux_rho_face_xp,
1516 storage.flux_rho_face_xm,
1517 storage.flux_rho_face_yp,
1518 storage.flux_rho_face_ym,
1519 storage.flux_rho_face_zp,
1520 storage.flux_rho_face_zm,
1521 storage.flux_rhov_face_xp,
1522 storage.flux_rhov_face_xm,
1523 storage.flux_rhov_face_yp,
1524 storage.flux_rhov_face_ym,
1525 storage.flux_rhov_face_zp,
1526 storage.flux_rhov_face_zm,
1527 storage.flux_rhoe_face_xp,
1528 storage.flux_rhoe_face_xm,
1529 storage.flux_rhoe_face_yp,
1530 storage.flux_rhoe_face_ym,
1531 storage.flux_rhoe_face_zp,
1532 storage.flux_rhoe_face_zm,
1533 storage.dtrho,
1534 storage.dtrhov,
1535 storage.dtrhoe);
1536 solver_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1537 }
1538
1539 if (solver_config.is_dust_on()) {
1541 AMRBlock::block_size,
1542 solver_config.grid_coord_to_pos_fact,
1543 solver_config.dust_config.ndust};
1544 node.set_edges(
1545 storage.block_counts,
1546 storage.cell_graph_edge,
1547 storage.block_cell_sizes,
1548 storage.cell0block_aabb_lower,
1549 storage.flux_rho_dust_face_xp,
1550 storage.flux_rho_dust_face_xm,
1551 storage.flux_rho_dust_face_yp,
1552 storage.flux_rho_dust_face_ym,
1553 storage.flux_rho_dust_face_zp,
1554 storage.flux_rho_dust_face_zm,
1555 storage.flux_rhov_dust_face_xp,
1556 storage.flux_rhov_dust_face_xm,
1557 storage.flux_rhov_dust_face_yp,
1558 storage.flux_rhov_dust_face_ym,
1559 storage.flux_rhov_dust_face_zp,
1560 storage.flux_rhov_dust_face_zm,
1561 storage.dtrho_dust,
1562 storage.dtrhov_dust);
1563 solver_sequence.push_back(std::make_shared<decltype(node)>(std::move(node)));
1564 }
1565
1566 shamrock::solvergraph::OperationSequence seq("Solver", std::move(solver_sequence));
1567 storage.solver_sequence = std::make_shared<decltype(seq)>(std::move(seq));
1568
1569 if (false) {
1570 logger::raw_ln(" -- tex:\n" + shambase::get_check_ref(storage.solver_sequence).get_tex());
1572 " -- dot:\n" + shambase::get_check_ref(storage.solver_sequence).get_dot_graph());
1573 }
1574}
1575
1576template<class Tvec, class TgridVec>
1578
1579 StackEntry stack_loc{};
1580
1581 // has to be first since there is a barrier that may mess the other timers
1582 shamsys::SystemMetrics system_metrics_start = shamsys::get_system_metrics();
1583
1585 f64 mpi_timer_start = shamcomm::mpi::get_timer("total");
1586
1587 for (auto &callbacks : timestep_callbacks) {
1588 if (callbacks.step_begin_callback) {
1589 shambase::get_check_ref(callbacks.step_begin_callback)();
1590 }
1591 }
1592
1593 Tscal t_current = get_time();
1594 Tscal dt_input = get_dt();
1595
1596 if (shamcomm::world_rank() == 0) {
1597 logger::normal_ln("amr::Godunov", sham::format("t = {}, dt = {}", t_current, dt_input));
1598 }
1599
1600 if (solver_config.face_half_time_interpolation) {
1601 shambase::get_check_ref(storage.dt_over2).value = dt_input / 2.0;
1602 }
1603
1604 shambase::Timer tstep;
1605 tstep.start();
1606
1607 // Scheduler step
1608 auto update_load_val = [&]() {
1609 shamlog_debug_ln("ComputeLoadBalanceValue", "update load balancing");
1610 scheduler().update_local_load_value([&](shamrock::patch::Patch p) {
1611 return scheduler().patch_data.owned_data.get(p.id_patch).get_obj_cnt();
1612 });
1613 };
1614 update_load_val();
1615 scheduler().scheduler_step(true, true);
1616 update_load_val();
1617 scheduler().scheduler_step(false, false);
1618
1619 if (solver_config.should_compute_rho_mean()) {
1620 auto [bmin, bmax] = scheduler().template get_box_volume<TgridVec>();
1621 Tscal dxfact = solver_config.grid_coord_to_pos_fact;
1622 Tvec dV = (bmax - bmin).template convert<Tscal>() * dxfact;
1623 Tscal Vsim = dV.x() * dV.y() * dV.z();
1624 shambase::get_check_ref(storage.simulation_volume).value = Vsim;
1625 }
1626
1628
1629 scheduler().for_each_patchdata_nonempty(
1631 storage.source_patches->patchdatas.add_obj(p.id_patch, std::ref(pdat));
1632 });
1633
1634 {
1635 shamrock::patch::SimulationBoxInfo &sim_box = scheduler().get_sim_box();
1636 auto [bmin, bmax] = sim_box.get_bounding_box<TgridVec>();
1637
1638 shambase::get_check_ref(storage.sim_box_edge).value = shammath::AABB<TgridVec>(bmin, bmax);
1639 }
1640
1641 SerialPatchTree<TgridVec> _sptree = SerialPatchTree<TgridVec>::build(scheduler());
1642 _sptree.attach_buf();
1643 storage.serial_patch_tree.set(std::move(_sptree));
1644
1648
1649 {
1650 auto &sim_box = scheduler().get_sim_box();
1651 auto transf = sim_box.template get_patch_transform<TgridVec>();
1652
1653 auto &global_patch_boxes_edge
1654 = shambase::get_check_ref(storage.solver_graph.template get_edge_ptr<
1656 "global_patch_boxes"));
1657
1658 global_patch_boxes_edge.values = {};
1659
1660 scheduler().for_each_global_patch([&](const shamrock::patch::Patch &p) {
1661 auto pbounds = transf.to_obj_coord(p);
1662 global_patch_boxes_edge.values.add_obj(
1663 p.id_patch, shammath::AABB<TgridVec>{pbounds.lower, pbounds.upper});
1664 });
1665 }
1666
1667 {
1668 auto &sim_box = scheduler().get_sim_box();
1669 auto transf = sim_box.template get_patch_transform<TgridVec>();
1670
1671 auto &local_patch_ids = shambase::get_check_ref(storage.local_patch_ids);
1672
1673 local_patch_ids.data = {};
1674
1675 scheduler().for_each_local_patch([&](const shamrock::patch::Patch &p) {
1676 local_patch_ids.data.push_back(p.id_patch);
1677 });
1678 }
1679
1683
1684 // Solvergraph evaluation
1685 {
1686 shambase::get_check_ref(storage.solver_sequence).evaluate();
1687 }
1688
1689 // RK2 + flux lim
1690 if (solver_config.drag_config.drag_solver_config == DragSolverMode::NoDrag) {
1691 modules::TimeIntegrator dt_integ(context, solver_config, storage);
1692 dt_integ.forward_euler(dt_input);
1693 } else if (solver_config.drag_config.drag_solver_config == DragSolverMode::IRK1) {
1694 modules::DragIntegrator drag_integ(context, solver_config, storage);
1695 drag_integ.involve_with_no_src(dt_input);
1696 drag_integ.enable_irk1_drag_integrator(dt_input);
1697 } else if (solver_config.drag_config.drag_solver_config == DragSolverMode::EXPO) {
1698 modules::DragIntegrator drag_integ(context, solver_config, storage);
1699 drag_integ.involve_with_no_src(dt_input);
1700 drag_integ.enable_expo_drag_integrator(dt_input);
1701 } else {
1703 }
1704
1705 if (!solver_config.amr_mode.old_amr) {
1707 node_copy_rho.set_edges(storage.refs_rho, storage.rho_primitive);
1708 node_copy_rho.evaluate();
1709 }
1710
1711 // TODO: check if we can drop that
1712 if (dt_input > 0) {
1713 modules::AMRGridRefinementHandler refinement(context, solver_config, storage);
1714 if (solver_config.amr_mode.old_amr) {
1715 refinement.update_refinement_old();
1716 } else {
1717 using AMRmode_None = typename AMRMode<Tvec, TgridVec>::None;
1718 if (std::get_if<AMRmode_None>(&solver_config.amr_mode.config) == nullptr) {
1720 node_cpy_rho_snap.set_edges(storage.refs_rho, storage.rho_snap);
1721 node_cpy_rho_snap.evaluate();
1723 node_cpy_rhoe_snap.set_edges(storage.refs_rhoe, storage.rhoe_snap);
1724 node_cpy_rhoe_snap.evaluate();
1726 node_cpy_rho_vel_snap.set_edges(storage.refs_rhov, storage.rho_vel_snap);
1727 node_cpy_rho_vel_snap.evaluate();
1728 }
1729
1730 // Then refinement pass
1731 refinement.update_refinement_new();
1732 }
1733 }
1734
1735 modules::ComputeCFL cfl_compute(context, solver_config, storage);
1736 f64 new_dt = cfl_compute.compute_cfl();
1737
1738 // if new physics like dust is added then use the smallest dt
1739 if (solver_config.is_dust_on())
1740 new_dt = std::min(new_dt, cfl_compute.compute_dust_cfl());
1741
1742 set_next_dt(new_dt);
1743 set_time(t_current + dt_input);
1744
1745 if (solver_config.drag_config.drag_solver_config != DragSolverMode::NoDrag) {
1746 storage.rho_next_no_drag.reset();
1747 storage.rhov_next_no_drag.reset();
1748 storage.rhoe_next_no_drag.reset();
1749 storage.rho_d_next_no_drag.reset();
1750 storage.rhov_d_next_no_drag.reset();
1751 }
1752
1753 storage.merge_patch_bounds.reset();
1754
1755 storage.ghost_zone_infos.reset();
1756
1757 storage.serial_patch_tree.reset();
1758
1759 shambase::get_check_ref(storage.source_patches).free_alloc();
1760
1761 shambase::get_check_ref(storage.exchange_gz_edge).free_alloc();
1762 shambase::get_check_ref(storage.idx_in_ghost).free_alloc();
1763
1764 shambase::get_check_ref(storage.ghost_layers_candidates_edge).free_alloc();
1765
1766 tstep.stop();
1767
1768 for (auto it = timestep_callbacks.rbegin(); it != timestep_callbacks.rend(); ++it) {
1769 if (it->step_end_callback) {
1770 shambase::get_check_ref(it->step_end_callback)();
1771 }
1772 }
1773
1775
1776 shamsys::SystemMetrics system_metrics_end = shamsys::get_system_metrics();
1777 shamsys::SystemMetrics system_metrics_delta = system_metrics_end - system_metrics_start;
1778
1779 f64 delta_mpi_timer = shamcomm::mpi::get_timer("total") - mpi_timer_start;
1780 f64 t_dev_alloc
1781 = (mem_perf_infos_end.time_alloc_device - mem_perf_infos_start.time_alloc_device)
1782 + (mem_perf_infos_end.time_free_device - mem_perf_infos_start.time_free_device);
1783 f64 t_host_alloc = (mem_perf_infos_end.time_alloc_host - mem_perf_infos_start.time_alloc_host)
1784 + (mem_perf_infos_end.time_free_host - mem_perf_infos_start.time_free_host);
1785
1786 u64 rank_count = scheduler().get_rank_count() * AMRBlock::block_size;
1787 f64 rate = f64(rank_count) / tstep.elapsed_sec();
1788
1789 u64 npatch = scheduler().patch_list.local.size();
1790
1791 std::string log_step = report_perf_timestep(
1792 rate,
1793 rank_count,
1794 npatch,
1795 tstep.elapsed_sec(),
1796 delta_mpi_timer,
1797 t_dev_alloc,
1798 t_host_alloc,
1799 mem_perf_infos_end.max_allocated_byte_device,
1800 mem_perf_infos_end.max_allocated_byte_host,
1801 system_metrics_delta,
1802 shamsys::has_reporter());
1803
1804 if (shamcomm::world_rank() == 0) {
1805 logger::info_ln("amr::RAMSES", log_step);
1807 "amr::RAMSES",
1808 "estimated rate :",
1809 dt_input * (3600 / tstep.elapsed_sec()),
1810 "(tsim/hr)");
1811 }
1812
1813 solve_logs.register_log(
1814 {t_current, // f64 solver_t;
1815 dt_input, // f64 solver_dt;
1816 shamcomm::world_rank(), // i32 world_rank;
1817 rank_count, // u64 rank_count;
1818 rate, // f64 rate;
1819 tstep.elapsed_sec(), // f64 elapsed_sec;
1821 system_metrics_delta});
1822
1823 storage.timings_details.reset();
1824}
1825
1826template<class Tvec, class TgridVec>
1827void shammodels::basegodunov::Solver<Tvec, TgridVec>::do_debug_vtk_dump(std::string filename) {
1828
1829 StackEntry stack_loc{};
1830 shamrock::LegacyVtkWriter writer(filename, true, shamrock::UnstructuredGrid);
1831
1832 PatchScheduler &sched = shambase::get_check_ref(context.sched);
1833
1834 u32 block_size = Solver::AMRBlock::block_size;
1835
1836 u64 num_obj = sched.get_rank_count();
1837
1838 std::unique_ptr<sycl::buffer<TgridVec>> pos1 = sched.rankgather_field<TgridVec>(0);
1839 std::unique_ptr<sycl::buffer<TgridVec>> pos2 = sched.rankgather_field<TgridVec>(1);
1840
1841 sycl::buffer<Tvec> pos_min_cell(num_obj * block_size);
1842 sycl::buffer<Tvec> pos_max_cell(num_obj * block_size);
1843
1844 shamsys::instance::get_compute_queue().submit([&, block_size](sycl::handler &cgh) {
1845 sycl::accessor acc_p1{shambase::get_check_ref(pos1), cgh, sycl::read_only};
1846 sycl::accessor acc_p2{shambase::get_check_ref(pos2), cgh, sycl::read_only};
1847 sycl::accessor cell_min{pos_min_cell, cgh, sycl::write_only, sycl::no_init};
1848 sycl::accessor cell_max{pos_max_cell, cgh, sycl::write_only, sycl::no_init};
1849
1850 using Block = typename Solver::AMRBlock;
1851
1852 shambase::parallel_for(cgh, num_obj, "rescale cells", [=](u64 id_a) {
1853 Tvec block_min = acc_p1[id_a].template convert<Tscal>();
1854 Tvec block_max = acc_p2[id_a].template convert<Tscal>();
1855
1856 Tvec delta_cell = (block_max - block_min) / Block::side_size;
1857#pragma unroll
1858 for (u32 ix = 0; ix < Block::side_size; ix++) {
1859#pragma unroll
1860 for (u32 iy = 0; iy < Block::side_size; iy++) {
1861#pragma unroll
1862 for (u32 iz = 0; iz < Block::side_size; iz++) {
1863 u32 i = Block::get_index({ix, iy, iz});
1864 Tvec delta_val = delta_cell * Tvec{ix, iy, iz};
1865 cell_min[id_a * block_size + i] = block_min + delta_val;
1866 cell_max[id_a * block_size + i] = block_min + (delta_cell) + delta_val;
1867 }
1868 }
1869 }
1870 });
1871 });
1872
1873 writer.write_voxel_cells(pos_min_cell, pos_max_cell, num_obj * block_size);
1874
1875 writer.add_cell_data_section();
1876 writer.add_field_data_section(3);
1877
1878 std::unique_ptr<sycl::buffer<Tscal>> fields_rho = sched.rankgather_field<Tscal>(2);
1879 writer.write_field("rho", fields_rho, num_obj * block_size);
1880
1881 std::unique_ptr<sycl::buffer<Tvec>> fields_vel = sched.rankgather_field<Tvec>(3);
1882 writer.write_field("rhovel", fields_vel, num_obj * block_size);
1883
1884 std::unique_ptr<sycl::buffer<Tscal>> fields_eint = sched.rankgather_field<Tscal>(4);
1885 writer.write_field("rhoetot", fields_eint, num_obj * block_size);
1886 /*
1887 std::unique_ptr<sycl::buffer<Tvec>> grad_rho
1888 = storage.grad_rho.get().rankgather_computefield(sched);
1889 writer.write_field("grad_rho", grad_rho, num_obj * block_size);
1890
1891 std::unique_ptr<sycl::buffer<Tvec>> dx_v =
1892 storage.dx_v.get().rankgather_computefield(sched); writer.write_field("dx_v", dx_v,
1893 num_obj * block_size);
1894
1895 std::unique_ptr<sycl::buffer<Tvec>> dy_v =
1896 storage.dy_v.get().rankgather_computefield(sched); writer.write_field("dy_v", dy_v,
1897 num_obj * block_size);
1898
1899 std::unique_ptr<sycl::buffer<Tvec>> dz_v =
1900 storage.dz_v.get().rankgather_computefield(sched); writer.write_field("dz_v", dz_v,
1901 num_obj * block_size);
1902
1903 std::unique_ptr<sycl::buffer<Tvec>> grad_P
1904 = storage.grad_P.get().rankgather_computefield(sched);
1905 writer.write_field("grad_P", grad_P, num_obj * block_size);
1906 */
1907
1908 /*
1909 std::unique_ptr<sycl::buffer<Tscal>> dtrho =
1910 storage.dtrho.get().rankgather_computefield(sched); writer.write_field("dtrho", dtrho,
1911 num_obj * block_size);
1912
1913 std::unique_ptr<sycl::buffer<Tvec>> dtrhov
1914 = storage.dtrhov.get().rankgather_computefield(sched);
1915 writer.write_field("dtrhov", dtrhov, num_obj * block_size);
1916
1917 std::unique_ptr<sycl::buffer<Tscal>> dtrhoe
1918 = storage.dtrhoe.get().rankgather_computefield(sched);
1919 writer.write_field("dtrhoe", dtrhoe, num_obj * block_size);
1920 */
1921}
1922
Field variant object to instanciate a variant on the patch types.
Computes the coordinates of each cell.
Field variant object to instanciate a variant on the patch types.
Field variant object to instanciate a variant on the patch types.
Defines the CopyPatchDataField class for copying fields between patch data field references.
Defines the CopyPatchDataLayerFields class for copying fields between patch data layers.
Per cell Euler time derivatives of the dust primitive state.
Per cell Euler time derivatives of the gas primitive state.
Solver graph node for exchanging ghost layer data between distributed processes.
Defines the ExtractCounts class for extracting object counts from patch data layer references.
Extract the ghost layer from the patch data layers.
Field variant object to instanciate a variant on the patch types.
A solver graph node to fuse a ghost layer into a set of patch data layers.
Defines the GetFieldRefFromLayer class for extracting field references from patch data layers.
Field variant object to instanciate a variant on the patch types.
Field variant object to instanciate a variant on the patch types.
sycl::queue & get_compute_queue(u32 id=0)
Node that applies a custom function to modify connected edges.
Shared distributed data layer for patch data management in solver graphs.
Defines the PatchDataLayerEdge class for managing patch data layer edges.
Declare a class to register and retrieve nodes and edges from a unique container.
Sum the fluxes into the time derivative fields for Dust.
Sum the fluxes into the time derivative fields for Hydro.
Field variant object to instanciate a variant on the patch types.
double f64
Alias for double.
std::uint32_t u32
32 bit unsigned integer
std::uint64_t u64
64 bit unsigned integer
void _impl_evaluate_internal()
evaluate the node
Definition Solver.cpp:103
The MPI scheduler.
A buffer allocated in USM (Unified Shared Memory).
void append(const DeviceBuffer &other)
Append the content of another buffer to this one.
size_t get_size() const
Gets the number of elements in the buffer.
sycl::buffer< T > copy_to_sycl_buffer() const
Copy the content of the buffer to a new SYCL buffer.
Class Timer measures the time elapsed since the timer was started.
Definition Timer.hpp:36
f64 elapsed_sec() const
Converts the stored nanosecond time to a floating point representation in seconds.
Definition Timer.hpp:88
void start()
Starts the timer.
Definition Timer.hpp:51
void stop()
Stops the timer and stores the elapsed time in nanoseconds.
Definition Timer.hpp:65
Compute the Euler time derivatives of the dust primitive state per cell.
Compute the Euler time derivatives of the gas primitive state per cell.
u32 get_field_idx(const std::string &field_name) const
Get the field id if matching name & type.
void add_field(const std::string &field_name, u32 nvar, SourceLocation loc=SourceLocation{})
add a field of type T to the layout
void for_each_field_any(Functor &&func) const
for each visit of each field
PatchDataLayer container class, the layout is described in patchdata_layout.
Store the information related to the size of the simulation box to convert patch integer coordinates ...
Definition SimBox.hpp:36
std::tuple< T, T > get_bounding_box() const
Get the stored bounding box of the domain.
Definition SimBox.hpp:248
A solver graph node that copies field data from source field references to target fields.
Inode is node between data edges, takes multiple inputs, multiple outputs.
Definition INode.hpp:31
void evaluate()
Evaluate the node.
Definition INode.hpp:156
void __internal_set_rw_edges(std::vector< std::shared_ptr< IEdge > > new_rw_edges)
Set the read write edges.
Definition INode.hpp:249
void __internal_set_ro_edges(std::vector< std::shared_ptr< IEdge > > new_ro_edges)
Set the read only edges.
Definition INode.hpp:238
virtual std::string _impl_get_label() const =0
get the label of the node
const T & get_ro_edge(int slot)
Get a read only edge and cast it to the type T.
Definition INode.hpp:91
virtual std::string _impl_get_tex() const =0
get the tex of the node
A node that simply frees the allocation of the connected node.
A node that applies a custom function to modify connected edges.
std::optional< std::reference_wrapper< SerialPatchTree< Tvec > > > patch_tree
The patch tree.
A graph container for managing solver nodes and edges with type-safe access.
std::shared_ptr< INode > & get_node_ptr_base(const std::string &name)
Retrieve a node by name as a shared pointer to the base interface.
std::shared_ptr< T > get_edge_ptr(const std::string &name)
Get a typed shared pointer to an edge by name.
std::shared_ptr< T > register_edge(const std::string &name, T &&edge)
Register an edge with automatic type deduction and shared pointer creation.
std::shared_ptr< T > register_node(const std::string &name, T &&node)
Register a node with automatic type deduction and shared pointer creation.
T & get_node_ref(const std::string &name)
Get a typed reference to a node by name.
This header file contains utility functions related to exception handling in the code.
MPI string gather / allgather helpers (declarations; implementations in shamalgs/src/collective/gathe...
MemPerfInfos get_mem_perf_info()
Retrieve the memory performance information.
void kernel_call(sham::DeviceQueue &q, RefIn in, RefOut in_out, u32 n, Functor &&func, SourceLocation &&callsite=SourceLocation{})
Submit a kernel to a SYCL queue.
void throw_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Throw an exception and append the source location to it.
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
void throw_unimplemented(SourceLocation loc=SourceLocation{})
Throw a std::runtime_error saying that the function is unimplemented.
i32 world_rank()
Gives the rank of the current process in the MPI communicator.
Definition worldInfo.cpp:41
@ HB
Huang and Bai. Pressureless Riemann solver by Huang and Bai (2022) in Athena++.
@ DHLL
Dust HLL. This is merely the HLL solver for dust. It's then a Rusanov like.
void normal_ln(std::string module_name, Types... var2)
Prints a log message with multiple arguments followed by a newline.
Definition logs.hpp:132
void raw_ln(Types... var2)
Prints a log message with multiple arguments followed by a newline.
Definition logs.hpp:89
void info_ln(std::string module_name, Types... var2)
Prints a log message with multiple arguments followed by a newline.
Definition logs.hpp:132
#define __shamrock_stack_entry()
Macro to create a stack entry.
shambase::details::BasicStackEntry StackEntry
Alias for shambase::details::BasicStackEntry.
f64 get_wtime()
Returns the current wall clock time in seconds.
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.
A class that references multiple buffers or similar objects.
Definition MultiRef.hpp:33
Axis-Aligned bounding box.
Definition AABB.hpp:99
utility class to handle AMR blocks
Definition AMRBlock.hpp:35
Patch object that contain generic patch information.
Definition Patch.hpp:33
f64 get_timer(std::string timername)
get a timer value
Definition wrapper.cpp:46