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