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