Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
Solver.cpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
16
17#include "shambase/assert.hpp"
18#include "shambase/exception.hpp"
19#include "shambase/logs/loglevels.hpp"
20#include "shambase/memory.hpp"
21#include "shambase/numeric_limits.hpp"
22#include "shambase/string.hpp"
23#include "shambase/tabulate.hpp"
24#include "shambase/time.hpp"
25#include "shamalgs/collective/exchanges.hpp"
26#include "shamalgs/collective/gather_str.hpp"
27#include "shamalgs/collective/reduction.hpp"
28#include "shamalgs/reduction.hpp"
29#include "shambackends/MemPerfInfos.hpp"
30#include "shambackends/details/memoryHandle.hpp"
31#include "shambackends/kernel_call.hpp"
32#include "shambackends/math.hpp"
33#include "shamcomm/logs.hpp"
34#include "shamcomm/worldInfo.hpp"
35#include "shamcomm/wrapper.hpp"
36#include "shammath/sphkernels.hpp"
37#include "shammodels/common/modules/ComputeGravWave.hpp"
38#include "shammodels/common/modules/ForwardEuler.hpp"
39#include "shammodels/common/modules/ForwardEulerPositive.hpp"
40#include "shammodels/common/timestep_report.hpp"
41#include "shammodels/sph/BasicSPHGhosts.hpp"
42#include "shammodels/sph/SPHUtilities.hpp"
43#include "shammodels/sph/Solver.hpp"
44#include "shammodels/sph/SolverConfig.hpp"
45#include "shammodels/sph/io/PhantomDump.hpp"
46#include "shammodels/sph/math/density.hpp"
47#include "shammodels/sph/math/forces.hpp"
48#include "shammodels/sph/math/q_ab.hpp"
49#include "shammodels/sph/modules/BallabioTsLimiter.hpp"
50#include "shammodels/sph/modules/BuildTrees.hpp"
51#include "shammodels/sph/modules/ComputeCFLCourant.hpp"
52#include "shammodels/sph/modules/ComputeCFLDivBCleaning.hpp"
53#include "shammodels/sph/modules/ComputeCFLDust1Fluid.hpp"
54#include "shammodels/sph/modules/ComputeCFLDustDrift.hpp"
55#include "shammodels/sph/modules/ComputeCFLForce.hpp"
56#include "shammodels/sph/modules/ComputeCFLSinkSink.hpp"
57#include "shammodels/sph/modules/ComputeEos.hpp"
58#include "shammodels/sph/modules/ComputeLoadBalanceValue.hpp"
59#include "shammodels/sph/modules/ComputeLuminosity.hpp"
60#include "shammodels/sph/modules/ComputeNeighStats.hpp"
61#include "shammodels/sph/modules/ComputeOmega.hpp"
62#include "shammodels/sph/modules/ConservativeCheck.hpp"
63#include "shammodels/sph/modules/DiffOperator.hpp"
64#include "shammodels/sph/modules/DiffOperatorDtDivv.hpp"
65#include "shammodels/sph/modules/ExternalForces.hpp"
66#include "shammodels/sph/modules/GetParticlesOutsideSphere.hpp"
67#include "shammodels/sph/modules/IterateSmoothingLengthDensity.hpp"
68#include "shammodels/sph/modules/IterateSmoothingLengthDensityNeighLim.hpp"
69#include "shammodels/sph/modules/KillParticles.hpp"
70#include "shammodels/sph/modules/LoopSmoothingLengthIter.hpp"
71#include "shammodels/sph/modules/MonoFluidTVADeltav.hpp"
72#include "shammodels/sph/modules/NeighbourCache.hpp"
73#include "shammodels/sph/modules/NodeComputePressureGrad.hpp"
74#include "shammodels/sph/modules/NodeMonofluidTVADustDensityClamp.hpp"
75#include "shammodels/sph/modules/ParticleReordering.hpp"
76#include "shammodels/sph/modules/SetDustStoppingTimeConstant.hpp"
77#include "shammodels/sph/modules/SetDustStoppingTimeEpstein.hpp"
78#include "shammodels/sph/modules/SinkParticlesAccreteQuantities.hpp"
79#include "shammodels/sph/modules/SinkParticlesEvictAccretedParticles.hpp"
80#include "shammodels/sph/modules/SinkParticlesFlagAccreteHard.hpp"
81#include "shammodels/sph/modules/SinkSelfGravityHost.hpp"
82#include "shammodels/sph/modules/UpdateDerivs.hpp"
83#include "shammodels/sph/modules/UpdateViscosity.hpp"
84#include "shammodels/sph/modules/io/VTKDump.hpp"
85#include "shammodels/sph/modules/self_gravity/SGDirectPlummer.hpp"
86#include "shammodels/sph/modules/self_gravity/SGFMMPlummer.hpp"
87#include "shammodels/sph/modules/self_gravity/SGMMPlummer.hpp"
88#include "shammodels/sph/modules/self_gravity/SGSFMMPlummer.hpp"
89#include "shammodels/sph/sink_edges_helper.hpp"
90#include "shammodels/sph/solvergraph/NeighCache.hpp"
91#include "shamphys/mhd.hpp"
92#include "shamrock/patch/Patch.hpp"
93#include "shamrock/patch/PatchDataLayer.hpp"
94#include "shamrock/patch/PatchDataLayerLayout.hpp"
95#include "shamrock/scheduler/ComputeField.hpp"
96#include "shamrock/scheduler/InterfacesUtility.hpp"
97#include "shamrock/scheduler/PatchScheduler.hpp"
98#include "shamrock/scheduler/ReattributeDataUtility.hpp"
99#include "shamrock/scheduler/SchedulerUtility.hpp"
100#include "shamrock/scheduler/SerialPatchTree.hpp"
101#include "shamrock/solvergraph/CopyPatchDataFieldFromLayer.hpp"
102#include "shamrock/solvergraph/DistributedBuffers.hpp"
103#include "shamrock/solvergraph/Field.hpp"
104#include "shamrock/solvergraph/FieldRefs.hpp"
105#include "shamrock/solvergraph/GetFieldRefFromLayer.hpp"
106#include "shamrock/solvergraph/GetObjCntFromLayer.hpp"
107#include "shamrock/solvergraph/IFieldRefs.hpp"
108#include "shamrock/solvergraph/Indexes.hpp"
109#include "shamrock/solvergraph/PatchDataLayerRefs.hpp"
110#include "shamrock/solvergraph/RankGetter.hpp"
111#include "shamrock/solvergraph/ScalarsEdge.hpp"
112#include "shamsolvergraph/SolverGraph.hpp"
113#include "shamsolvergraph/edge/IDataEdge.hpp"
114#include "shamsolvergraph/edge/IDataEdgeSerializable.hpp"
115#include "shamsolvergraph/node/ForwardEulerHost.hpp"
116#include "shamsolvergraph/node/ForwardEulerHost2Deriv.hpp"
117#include "shamsolvergraph/node/NodeFreeAlloc.hpp"
118#include "shamsolvergraph/node/NodeMapEdge.hpp"
119#include "shamsolvergraph/node/NodeSetEdge.hpp"
120#include "shamsolvergraph/node/OperationIf.hpp"
121#include "shamsolvergraph/node/OperationSequence.hpp"
122#include "shamsys/NodeInstance.hpp"
123#include "shamsys/legacy/log.hpp"
124#include "shamsys/system_metrics.hpp"
125#include "shamtree/KarrasRadixTreeField.hpp"
126#include "shamtree/TreeTraversalCache.hpp"
127#include <memory>
128#include <stdexcept>
129#include <vector>
130
131namespace shambase {
132
133 template<class T>
134 std::shared_ptr<T> to_shared(T &&t) {
135 return std::make_shared<T>(std::forward<T>(t));
136 }
137} // namespace shambase
138
139namespace shammodels::sph {
140
141 namespace {
143 template<class Tvec>
144 std::shared_ptr<shamrock::solvergraph::INode> build_self_gravity_node(
145 SelfGravConfig &self_grav_config,
146 std::shared_ptr<shamrock::solvergraph::Indexes<u32>> sizes,
147 std::shared_ptr<shamrock::solvergraph::IDataEdge<shambase::VecComponent<Tvec>>>
148 gpart_mass,
149 std::shared_ptr<shamrock::solvergraph::IDataEdge<shambase::VecComponent<Tvec>>>
150 constant_G,
151 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> field_xyz,
152 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> field_axyz_ext) {
153
154 using Tscal = shambase::VecComponent<Tvec>;
155
156 Tscal eps_grav = shambase::get_check_ref(
157 std::get_if<SelfGravConfig::SofteningPlummer>(
158 &self_grav_config.softening_mode))
159 .epsilon;
160
161 std::shared_ptr<shamrock::solvergraph::INode> sg_inode;
162
163 if (self_grav_config.is_none()) {
164 throw shambase::make_except_with_loc<std::runtime_error>(
165 "How did you get there ?\?\?!!!");
166 } else if (self_grav_config.is_direct()) {
167
168 SelfGravConfig::Direct &direct_config = shambase::get_check_ref(
169 std::get_if<SelfGravConfig::Direct>(&self_grav_config.config));
170
171 modules::SGDirectPlummer<Tvec> self_gravity_direct_node(
172 eps_grav, direct_config.reference_mode);
173 self_gravity_direct_node.set_edges(
174 sizes, gpart_mass, constant_G, field_xyz, field_axyz_ext);
175
176 sg_inode = shambase::to_shared(std::move(self_gravity_direct_node));
177
178 } else if (self_grav_config.is_mm()) {
179
180 SelfGravConfig::MM &mm_config = shambase::get_check_ref(
181 std::get_if<SelfGravConfig::MM>(&self_grav_config.config));
182
183 auto run_sg_mm = [&](auto mm_order_tag) {
184 constexpr u32 order = decltype(mm_order_tag)::value;
185 modules::SGMMPlummer<Tvec, order> self_gravity_mm_node(
186 eps_grav, mm_config.opening_angle, mm_config.reduction_level);
187 self_gravity_mm_node.set_edges(
188 sizes, gpart_mass, constant_G, field_xyz, field_axyz_ext);
189 sg_inode = shambase::to_shared(std::move(self_gravity_mm_node));
190 };
191
192 switch (mm_config.order) {
193 case 1 : run_sg_mm(std::integral_constant<u32, 1>{}); break;
194 case 2 : run_sg_mm(std::integral_constant<u32, 2>{}); break;
195 case 3 : run_sg_mm(std::integral_constant<u32, 3>{}); break;
196 case 4 : run_sg_mm(std::integral_constant<u32, 4>{}); break;
197 case 5 : run_sg_mm(std::integral_constant<u32, 5>{}); break;
198 default: shambase::throw_unimplemented();
199 }
200
201 } else if (self_grav_config.is_fmm()) {
202
203 SelfGravConfig::FMM &fmm_config = shambase::get_check_ref(
204 std::get_if<SelfGravConfig::FMM>(&self_grav_config.config));
205
206 auto run_sg_fmm = [&](auto fmm_order_tag) {
207 constexpr u32 order = decltype(fmm_order_tag)::value;
208 modules::SGFMMPlummer<Tvec, order> self_gravity_fmm_node(
209 eps_grav, fmm_config.opening_angle, fmm_config.reduction_level);
210 self_gravity_fmm_node.set_edges(
211 sizes, gpart_mass, constant_G, field_xyz, field_axyz_ext);
212 sg_inode = shambase::to_shared(std::move(self_gravity_fmm_node));
213 };
214
215 switch (fmm_config.order) {
216 case 1 : run_sg_fmm(std::integral_constant<u32, 1>{}); break;
217 case 2 : run_sg_fmm(std::integral_constant<u32, 2>{}); break;
218 case 3 : run_sg_fmm(std::integral_constant<u32, 3>{}); break;
219 case 4 : run_sg_fmm(std::integral_constant<u32, 4>{}); break;
220 case 5 : run_sg_fmm(std::integral_constant<u32, 5>{}); break;
221 default: shambase::throw_unimplemented();
222 }
223
224 } else if (self_grav_config.is_sfmm()) {
225
226 SelfGravConfig::SFMM &sfmm_config = shambase::get_check_ref(
227 std::get_if<SelfGravConfig::SFMM>(&self_grav_config.config));
228
229 auto run_sg_sfmm = [&](auto sfmm_order_tag) {
230 constexpr u32 order = decltype(sfmm_order_tag)::value;
231 modules::SGSFMMPlummer<Tvec, order> self_gravity_sfmm_node(
232 eps_grav,
233 sfmm_config.opening_angle,
234 sfmm_config.leaf_lowering,
235 sfmm_config.reduction_level);
236 self_gravity_sfmm_node.set_edges(
237 sizes, gpart_mass, constant_G, field_xyz, field_axyz_ext);
238 sg_inode = shambase::to_shared(std::move(self_gravity_sfmm_node));
239 };
240
241 switch (sfmm_config.order) {
242 case 1 : run_sg_sfmm(std::integral_constant<u32, 1>{}); break;
243 case 2 : run_sg_sfmm(std::integral_constant<u32, 2>{}); break;
244 case 3 : run_sg_sfmm(std::integral_constant<u32, 3>{}); break;
245 case 4 : run_sg_sfmm(std::integral_constant<u32, 4>{}); break;
246 case 5 : run_sg_sfmm(std::integral_constant<u32, 5>{}); break;
247 default: shambase::throw_unimplemented();
248 }
249
250 } else {
251 throw shambase::make_except_with_loc<std::runtime_error>(
252 "Self gravity config not supported, current state is : \n"
253 + nlohmann::json{self_grav_config}.dump(4));
254 }
255
256 return sg_inode;
257 }
258 } // namespace
259
260} // namespace shammodels::sph
261
262template<class Tvec, template<class> class Kern>
263void shammodels::sph::Solver<Tvec, Kern>::init_solver_graph() {
264
265 PatchScheduler &sched = scheduler();
266
267 auto &sync_data = sched.synchronized_data;
268
269 shamrock::patch::PatchDataLayerLayout &pdl = scheduler().pdl_old();
270 bool has_B_field = solver_config.has_field_B_on_rho();
271 bool has_psi_field = solver_config.has_field_psi_on_ch();
272 bool has_epsilon_field = solver_config.dust_config.has_epsilon_field();
273 bool has_deltav_field = solver_config.dust_config.has_deltav_field();
274 bool has_s_j_field = solver_config.dust_config.has_s_j_field();
275
276 using namespace shamrock::solvergraph;
277
278 SolverGraph &solver_graph = storage.solver_graph;
279
280 solver_graph.register_edge(
281 "scheduler_patchdata", PatchDataLayerRefs("patchdatas", "\\mathbb{U}_{\\rm patch}"));
282 solver_graph.register_edge("part_counts", Indexes<u32>("Npart_patch", "N_{\\rm part}_p"));
283
284 solver_graph.register_edge("dt_half", IDataEdge<Tscal>("dt_half", "\\frac{dt}{2}"));
285 solver_graph.register_edge("gpart_mass", IDataEdge<Tscal>("m", "m"));
286
287 solver_graph.register_edge("xyz", FieldRefs<Tvec>("xyz", "\\mathbf{r}"));
288 solver_graph.register_edge("vxyz", FieldRefs<Tvec>("vxyz", "\\mathbf{v}"));
289 solver_graph.register_edge("axyz", FieldRefs<Tvec>("axyz", "\\mathbf{a}"));
290 solver_graph.register_edge("uint", FieldRefs<Tscal>("uint", "u_{\\rm int}"));
291 solver_graph.register_edge("duint", FieldRefs<Tscal>("duint", "du_{\\rm int}"));
292 solver_graph.register_edge("hpart", FieldRefs<Tscal>("hpart", "h_{\\rm part}"));
293
294 if (has_B_field) {
295 solver_graph.register_edge("B/rho", FieldRefs<Tvec>("B/rho", "B_{\\rho}"));
296 solver_graph.register_edge("dB/rho", FieldRefs<Tvec>("dB/rho", "dB_{\\rho}"));
297 }
298 if (has_psi_field) {
299 solver_graph.register_edge("psi/ch", FieldRefs<Tscal>("psi/ch", "\\psi_{\\rm ch}"));
300 solver_graph.register_edge("dpsi/ch", FieldRefs<Tscal>("dpsi/ch", "d\\psi_{\\rm ch}"));
301 }
302 if (has_epsilon_field) {
303 solver_graph.register_edge("epsilon", FieldRefs<Tscal>("epsilon", "\\epsilon"));
304 solver_graph.register_edge("dtepsilon", FieldRefs<Tscal>("dtepsilon", "d\\epsilon"));
305 }
306 if (has_deltav_field) {
307 solver_graph.register_edge("deltav", FieldRefs<Tvec>("deltav", "\\Delta v"));
308 solver_graph.register_edge("dtdeltav", FieldRefs<Tvec>("dtdeltav", "d\\Delta v"));
309 }
310 if (has_s_j_field) {
311 solver_graph.register_edge("s_j", FieldRefs<Tscal>("s_j", "S_j"));
312 solver_graph.register_edge("ds_j_dt", FieldRefs<Tscal>("ds_j_dt", "dS_j/dt"));
313
314 u32 ndust = solver_config.dust_config.get_dust_nvar();
315 solver_graph.register_edge("Ts_j", Field<Tscal>(ndust, "Ts_j", "Ts_j"));
316 }
317
318 {
319 auto set_gpart_mass = solver_graph.register_node(
320 "set_gpart_mass", NodeSetEdge<IDataEdge<Tscal>>([&](IDataEdge<Tscal> &gpart_mass) {
321 gpart_mass.data = solver_config.gpart_mass;
322 }));
323 shambase::get_check_ref(set_gpart_mass)
324 .set_edges(solver_graph.get_edge_ptr<IDataEdge<Tscal>>("gpart_mass"));
325 }
326
328 // attach fields to scheduler
330 {
331 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> attach_field_sequence;
332
333 {
334 auto set_scheduler_patchdata = solver_graph.register_node(
335 "set_scheduler_patchdata",
336 NodeSetEdge<PatchDataLayerRefs>([&](PatchDataLayerRefs &scheduler_patchdata) {
337 scheduler_patchdata.free_alloc();
338 scheduler().for_each_patchdata_nonempty(
339 [&](const shamrock::patch::Patch &p,
340 shamrock::patch::PatchDataLayer &pdat) {
341 scheduler_patchdata.patchdatas.add_obj(p.id_patch, std::ref(pdat));
342 });
343 }));
344 shambase::get_check_ref(set_scheduler_patchdata)
345 .set_edges(solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"));
346 attach_field_sequence.push_back(set_scheduler_patchdata);
347 }
348
349 {
350 auto attach_part_counts
351 = solver_graph.register_node("attach_part_counts", GetObjCntFromLayer{});
352 shambase::get_check_ref(attach_part_counts)
353 .set_edges(
354 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
355 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"));
356 attach_field_sequence.push_back(attach_part_counts);
357 }
358
359 {
360 auto attach_xyz
361 = solver_graph.register_node("attach_xyz", GetFieldRefFromLayer<Tvec>(pdl, "xyz"));
362 shambase::get_check_ref(attach_xyz)
363 .set_edges(
364 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
365 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("xyz"));
366 attach_field_sequence.push_back(attach_xyz);
367 }
368
369 {
370 auto attach_vxyz = solver_graph.register_node(
371 "attach_vxyz", GetFieldRefFromLayer<Tvec>(pdl, "vxyz"));
372 shambase::get_check_ref(attach_vxyz)
373 .set_edges(
374 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
375 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("vxyz"));
376 attach_field_sequence.push_back(attach_vxyz);
377 }
378
379 {
380 auto attach_axyz = solver_graph.register_node(
381 "attach_axyz", GetFieldRefFromLayer<Tvec>(pdl, "axyz"));
382 shambase::get_check_ref(attach_axyz)
383 .set_edges(
384 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
385 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("axyz"));
386 attach_field_sequence.push_back(attach_axyz);
387 }
388
389 {
390 auto attach_uint = solver_graph.register_node(
391 "attach_uint", GetFieldRefFromLayer<Tscal>(pdl, "uint"));
392 shambase::get_check_ref(attach_uint)
393 .set_edges(
394 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
395 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("uint"));
396 attach_field_sequence.push_back(attach_uint);
397 }
398
399 {
400 auto attach_duint = solver_graph.register_node(
401 "attach_duint", GetFieldRefFromLayer<Tscal>(pdl, "duint"));
402 shambase::get_check_ref(attach_duint)
403 .set_edges(
404 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
405 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("duint"));
406 attach_field_sequence.push_back(attach_duint);
407 }
408
409 {
410 auto attach_hpart = solver_graph.register_node(
411 "attach_hpart", GetFieldRefFromLayer<Tscal>(pdl, "hpart"));
412 shambase::get_check_ref(attach_hpart)
413 .set_edges(
414 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
415 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("hpart"));
416 attach_field_sequence.push_back(attach_hpart);
417 }
418
419 if (has_B_field) {
420 auto attach_B_on_rho = solver_graph.register_node(
421 "attach_B_on_rho", GetFieldRefFromLayer<Tvec>(pdl, "B/rho"));
422 shambase::get_check_ref(attach_B_on_rho)
423 .set_edges(
424 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
425 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("B/rho"));
426 attach_field_sequence.push_back(attach_B_on_rho);
427 }
428
429 if (has_B_field) {
430 auto attach_dB_on_rho = solver_graph.register_node(
431 "attach_dB_on_rho", GetFieldRefFromLayer<Tvec>(pdl, "dB/rho"));
432 shambase::get_check_ref(attach_dB_on_rho)
433 .set_edges(
434 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
435 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("dB/rho"));
436 attach_field_sequence.push_back(attach_dB_on_rho);
437 }
438
439 if (has_psi_field) {
440 auto attach_psi_on_ch = solver_graph.register_node(
441 "attach_psi_on_ch", GetFieldRefFromLayer<Tscal>(pdl, "psi/ch"));
442 shambase::get_check_ref(attach_psi_on_ch)
443 .set_edges(
444 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
445 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("psi/ch"));
446 attach_field_sequence.push_back(attach_psi_on_ch);
447 }
448
449 if (has_psi_field) {
450 auto attach_dpsi_on_ch = solver_graph.register_node(
451 "attach_dpsi_on_ch", GetFieldRefFromLayer<Tscal>(pdl, "dpsi/ch"));
452 shambase::get_check_ref(attach_dpsi_on_ch)
453 .set_edges(
454 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
455 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("dpsi/ch"));
456 attach_field_sequence.push_back(attach_dpsi_on_ch);
457 }
458
459 if (has_epsilon_field) {
460 auto attach_epsilon = solver_graph.register_node(
461 "attach_epsilon", GetFieldRefFromLayer<Tscal>(pdl, "epsilon"));
462 shambase::get_check_ref(attach_epsilon)
463 .set_edges(
464 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
465 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("epsilon"));
466 attach_field_sequence.push_back(attach_epsilon);
467 }
468
469 if (has_epsilon_field) {
470 auto attach_dtepsilon = solver_graph.register_node(
471 "attach_dtepsilon", GetFieldRefFromLayer<Tscal>(pdl, "dtepsilon"));
472 shambase::get_check_ref(attach_dtepsilon)
473 .set_edges(
474 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
475 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("dtepsilon"));
476 attach_field_sequence.push_back(attach_dtepsilon);
477 }
478
479 if (has_deltav_field) {
480 auto attach_deltav = solver_graph.register_node(
481 "attach_deltav", GetFieldRefFromLayer<Tvec>(pdl, "deltav"));
482 shambase::get_check_ref(attach_deltav)
483 .set_edges(
484 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
485 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("deltav"));
486 attach_field_sequence.push_back(attach_deltav);
487 }
488
489 if (has_deltav_field) {
490 auto attach_dtdeltav = solver_graph.register_node(
491 "attach_dtdeltav", GetFieldRefFromLayer<Tvec>(pdl, "dtdeltav"));
492 shambase::get_check_ref(attach_dtdeltav)
493 .set_edges(
494 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
495 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("dtdeltav"));
496 attach_field_sequence.push_back(attach_dtdeltav);
497 }
498
499 if (has_s_j_field) {
500 auto attach_s_j
501 = solver_graph.register_node("attach_s_j", GetFieldRefFromLayer<Tscal>(pdl, "s_j"));
502 shambase::get_check_ref(attach_s_j)
503 .set_edges(
504 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
505 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("s_j"));
506 attach_field_sequence.push_back(attach_s_j);
507 }
508
509 if (has_s_j_field) {
510 auto attach_ds_j_dt = solver_graph.register_node(
511 "attach_ds_j_dt", GetFieldRefFromLayer<Tscal>(pdl, "ds_j_dt"));
512 shambase::get_check_ref(attach_ds_j_dt)
513 .set_edges(
514 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"),
515 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("ds_j_dt"));
516 attach_field_sequence.push_back(attach_ds_j_dt);
517 }
518 solver_graph.register_node(
519 "attach fields to scheduler",
520 OperationSequence("attach fields", std::move(attach_field_sequence)));
521 }
522
524 // leapfrog predictor
526
527 {
528
529 auto make_half_step_sequence = [&](std::string prefix) {
530 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> half_step_sequence;
531
532 {
533 auto half_step_vxyz = solver_graph.register_node(
534 prefix + "_vxyz", shammodels::common::modules::ForwardEuler<Tvec>{});
535 shambase::get_check_ref(half_step_vxyz)
536 .set_edges(
537 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"),
538 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("axyz"),
539 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
540 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("vxyz"));
541 half_step_sequence.push_back(half_step_vxyz);
542 }
543
544 {
545 auto half_step_uint = solver_graph.register_node(
546 prefix + "_uint", shammodels::common::modules::ForwardEuler<Tscal>{});
547 shambase::get_check_ref(half_step_uint)
548 .set_edges(
549 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"),
550 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("duint"),
551 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
552 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("uint"));
553 half_step_sequence.push_back(half_step_uint);
554 }
555
556 if (has_B_field) {
557 auto half_step_B_on_rho = solver_graph.register_node(
558 prefix + "_B_on_rho", shammodels::common::modules::ForwardEuler<Tvec>{});
559 shambase::get_check_ref(half_step_B_on_rho)
560 .set_edges(
561 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"),
562 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("dB/rho"),
563 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
564 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("B/rho"));
565 half_step_sequence.push_back(half_step_B_on_rho);
566 }
567
568 if (has_psi_field) {
569 auto half_step_psi_on_ch = solver_graph.register_node(
570 prefix + "_psi_on_ch", shammodels::common::modules::ForwardEuler<Tscal>{});
571 shambase::get_check_ref(half_step_psi_on_ch)
572 .set_edges(
573 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"),
574 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("dpsi/ch"),
575 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
576 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("psi/ch"));
577 half_step_sequence.push_back(half_step_psi_on_ch);
578 }
579
580 if (has_epsilon_field) {
581 auto half_step_epsilon = solver_graph.register_node(
582 prefix + "_epsilon", shammodels::common::modules::ForwardEuler<Tscal>{});
583 shambase::get_check_ref(half_step_epsilon)
584 .set_edges(
585 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"),
586 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("dtepsilon"),
587 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
588 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("epsilon"));
589 half_step_sequence.push_back(half_step_epsilon);
590 }
591
592 if (has_deltav_field) {
593 auto half_step_deltav = solver_graph.register_node(
594 prefix + "_deltav", shammodels::common::modules::ForwardEuler<Tvec>{});
595 shambase::get_check_ref(half_step_deltav)
596 .set_edges(
597 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"),
598 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("dtdeltav"),
599 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
600 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("deltav"));
601 half_step_sequence.push_back(half_step_deltav);
602 }
603
604 if (has_s_j_field) {
605 u32 ndust = solver_config.dust_config.get_dust_nvar();
606
607 auto &cfg = solver_config.dust_config.get_monofluid_tva();
608
609 if (cfg.ensure_s_j_positivity) {
610 auto half_step_s_j = solver_graph.register_node(
611 prefix + "_s_j",
612 shammodels::common::modules::ForwardEulerPositive<Tscal>(ndust));
613 shambase::get_check_ref(half_step_s_j)
614 .set_edges(
615 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"),
616 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("ds_j_dt"),
617 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
618 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("s_j"));
619 half_step_sequence.push_back(half_step_s_j);
620 } else {
621 auto half_step_s_j = solver_graph.register_node(
622 prefix + "_s_j", shammodels::common::modules::ForwardEuler<Tscal>(ndust));
623 shambase::get_check_ref(half_step_s_j)
624 .set_edges(
625 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"),
626 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("ds_j_dt"),
627 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
628 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("s_j"));
629 half_step_sequence.push_back(half_step_s_j);
630 }
631
632 if (cfg.should_clamp_dust_density()) {
633 auto hfactd_edge = IDataEdge<Tscal>::make_shared("hfactd", "hfactd");
634 hfactd_edge->data = Kernel::hfactd;
635
636 auto clamp_frac_edge
637 = IDataEdge<Tscal>::make_shared("clamp_frac", "clamp_frac");
638 clamp_frac_edge->data = cfg.get_clamp_dust_frac();
639
640 auto half_step_s_j_density_clamp = solver_graph.register_node(
641 prefix + "_s_j_density_clamp",
642 shammodels::sph::modules::NodeMonofluidTVADustDensityClamp<Tvec>(ndust));
643 shambase::get_check_ref(half_step_s_j_density_clamp)
644 .set_edges(
645 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
646 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("gpart_mass"),
647 hfactd_edge,
648 clamp_frac_edge,
649 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("hpart"),
650 solver_graph.get_edge_ptr<FieldRefs<Tscal>>("s_j"));
651 half_step_sequence.push_back(half_step_s_j_density_clamp);
652 }
653 }
654
655 return OperationSequence("half step", std::move(half_step_sequence));
656 };
657
658 solver_graph.register_node("half_step1", make_half_step_sequence("half_step1"));
659 solver_graph.register_node("half_step2", make_half_step_sequence("half_step2"));
660
661 {
662 auto full_step_xyz = solver_graph.register_node(
663 "full_step_xyz", shammodels::common::modules::ForwardEuler<Tvec>{});
664 shambase::get_check_ref(full_step_xyz)
665 .set_edges(
666 sync_data.get_edge_ptr<IDataEdge<Tscal>>("dt"),
667 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("vxyz"),
668 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
669 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("xyz"));
670 }
671
672 {
673 auto leapfrog_predictor = solver_graph.register_node(
674 "leapfrog predictor",
675 OperationSequence(
676 "leapfrog predictor",
677 {
678 solver_graph.get_node_ptr_base("half_step1"),
679 solver_graph.get_node_ptr_base("full_step_xyz"),
680 solver_graph.get_node_ptr_base("half_step2"),
681 }));
682 }
683 }
684
686 // Part killing step
688 bool do_part_killing_step = solver_config.particle_killing.kill_list.size() > 0;
689
690 if (do_part_killing_step) {
691
692 auto patchdatas = solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata");
693 auto xyz_edge = solver_graph.get_edge_ptr<FieldRefs<Tvec>>("xyz");
694
695 auto part_to_remove = solver_graph.register_edge(
696 "part_to_remove", DistributedBuffers<u32>("part_to_remove", "part_to_remove"));
697
698 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> part_kill_sequence{};
699
700 {
701
702 auto empty_part_to_remove
703 = solver_graph.register_node("empty_part_to_remove", NodeFreeAlloc{});
704 shambase::get_check_ref(empty_part_to_remove).set_edges(part_to_remove);
705 part_kill_sequence.push_back(empty_part_to_remove);
706 }
707
708 using kill_t = typename ParticleKillingConfig<Tvec>::kill_t;
709 using kill_sphere = typename ParticleKillingConfig<Tvec>::Sphere;
710
711 // selectors
712 for (kill_t &kill_obj : solver_config.particle_killing.kill_list) {
713 if (kill_sphere *kill_info = std::get_if<kill_sphere>(&kill_obj)) {
714
715 modules::GetParticlesOutsideSphere<Tvec> node_selector(
716 kill_info->center, kill_info->radius);
717 node_selector.set_edges(xyz_edge, part_to_remove);
718
719 part_kill_sequence.push_back(
720 std::make_shared<decltype(node_selector)>(std::move(node_selector)));
721 }
722 }
723
724 { // killing
725 modules::KillParticles node_killer{};
726 node_killer.set_edges(part_to_remove, patchdatas);
727
728 part_kill_sequence.push_back(
729 std::make_shared<decltype(node_killer)>(std::move(node_killer)));
730 }
731
732 // update part counts and spans since particles have been killed and thus
733 // patches can become empty
734 part_kill_sequence.push_back(solver_graph.get_node_ptr_base("attach fields to scheduler"));
735
736 solver_graph.register_node(
737 "part killing step",
738 OperationSequence("part killing step", std::move(part_kill_sequence)));
739 }
740
741 {
742 auto dt_to_half_dt = solver_graph.register_node(
743 "dt_to_half_dt",
744 NodeMapEdge<IDataEdge<Tscal>, IDataEdge<Tscal>>{
745 [](const IDataEdge<Tscal> &dt, IDataEdge<Tscal> &half_dt) {
746 half_dt.data = dt.data / 2;
747 }});
748 shambase::get_check_ref(dt_to_half_dt)
749 .set_edges(
750 sync_data.get_edge_ptr<IDataEdge<Tscal>>("dt"),
751 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"));
752 }
753
754 storage.part_counts
755 = std::make_shared<shamrock::solvergraph::Indexes<u32>>("part_counts", "N_{\\rm part}");
756
757 storage.part_counts_with_ghost = std::make_shared<shamrock::solvergraph::Indexes<u32>>(
758 "part_counts_with_ghost", "N_{\\rm part, with ghost}");
759
760 storage.patch_rank_owner = std::make_shared<shamrock::solvergraph::RankGetter>(
761 [&](u64 patch_id) -> u32 {
762 return scheduler().get_patch_rank_owner(patch_id);
763 },
764 "patch_rank_owner",
765 "rank");
766
767 // merged ghost spans
768 storage.positions_with_ghosts
769 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("part_pos", "\\mathbf{r}");
770 storage.hpart_with_ghosts
771 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("h_part", "h");
772
773 storage.neigh_cache
774 = std::make_shared<shammodels::sph::solvergraph::NeighCache>("neigh_cache", "neigh");
775
776 storage.omega = std::make_shared<shamrock::solvergraph::Field<Tscal>>(1, "omega", "\\Omega");
777
778 if (solver_config.has_field_alphaAV()) {
779 storage.alpha_av_updated = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
780 1, "alpha_av_updated", "\\alpha_{\\rm AV}");
781 }
782
783 storage.pressure = std::make_shared<shamrock::solvergraph::Field<Tscal>>(1, "pressure", "P");
784 storage.soundspeed
785 = std::make_shared<shamrock::solvergraph::Field<Tscal>>(1, "soundspeed", "c_s");
786
787 storage.exchange_gz_alpha
788 = std::make_shared<shamrock::solvergraph::ExchangeGhostField<Tscal>>();
789 storage.exchange_gz_node
790 = std::make_shared<shamrock::solvergraph::ExchangeGhostLayer>(storage.ghost_layout);
791 storage.exchange_gz_positions
792 = std::make_shared<shamrock::solvergraph::ExchangeGhostLayer>(storage.xyzh_ghost_layout);
793
795 // sink accretion
797 {
798 solver_graph.register_edge("has_sinks", IDataEdge<bool>("has_sinks", "has_sinks"));
799
800 auto set_has_sinks = solver_graph.register_node(
801 "set_has_sinks", NodeSetEdge<IDataEdge<bool>>([&](IDataEdge<bool> &has_sinks_edge) {
802 has_sinks_edge.data = has_sinks<Tvec>(sync_data);
803 }));
804 shambase::get_check_ref(set_has_sinks)
805 .set_edges(solver_graph.get_edge_ptr<IDataEdge<bool>>("has_sinks"));
806
807 auto free_xyz = solver_graph.register_node("free_xyz_refs", NodeFreeAlloc{});
808 shambase::get_check_ref(free_xyz).set_edges(
809 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("xyz"));
810
811 auto free_vxyz = solver_graph.register_node("free_vxyz_refs", NodeFreeAlloc{});
812 shambase::get_check_ref(free_vxyz).set_edges(
813 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("vxyz"));
814
815 auto free_axyz = solver_graph.register_node("free_axyz_refs", NodeFreeAlloc{});
816 shambase::get_check_ref(free_axyz).set_edges(
817 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("axyz"));
818
819 // sink synchronized edges, kept around as they are used by several nodes below
820 auto sink_positions
821 = sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_pos");
822 auto sink_velocities
823 = sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_vel");
824 auto sink_accelerations
825 = sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_acc_sph");
826 auto sink_angmom = sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>(
827 "sink_angular_momentum");
828 auto sink_mass
829 = sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tscal>>>("sink_mass");
830 auto sink_accr_radii = sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tscal>>>(
831 "sink_accretion_radius");
832
833 solver_graph.register_edge(
834 "sink_accretion_table", Field<u32>(1, "sink_accretion_table", "\\mathrm{acc}"));
835
836 auto flag_node = solver_graph.register_node(
837 "flag_accrete_hard", modules::SinkParticlesFlagAccreteHard<Tvec>{});
838 shambase::get_check_ref(flag_node).set_edges(
839 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
840 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("xyz"),
841 sink_positions,
842 sink_accr_radii,
843 solver_graph.get_edge_ptr<Field<u32>>("sink_accretion_table"));
844
845 auto qty_node = solver_graph.register_node(
846 "accrete_quantities", modules::SinkParticlesAccreteQuantities<Tvec>{});
847 shambase::get_check_ref(qty_node).set_edges(
848 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("gpart_mass"),
849 sync_data.get_edge_ptr<IDataEdge<Tscal>>("dt"),
850 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
851 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("xyz"),
852 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("vxyz"),
853 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("axyz"),
854 solver_graph.get_edge_ptr<Field<u32>>("sink_accretion_table"),
855 sink_positions,
856 sink_velocities,
857 sink_accelerations,
858 sink_angmom,
859 sink_mass);
860
861 auto evict_node = solver_graph.register_node(
862 "evict_accreted_particles", modules::SinkParticlesEvictAccretedParticles<Tvec>{});
863 shambase::get_check_ref(evict_node)
864 .set_edges(
865 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
866 solver_graph.get_edge_ptr<Field<u32>>("sink_accretion_table"),
867 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"));
868
869 auto if_has_sinks = solver_graph.register_node(
870 "if_has_sinks",
871 OperationSequence(
872 "if_has_sinks",
873 {
874 // the "time_step" sequence (set_gpart_mass, attach fields to scheduler, ...)
875 // has not run yet at this stage of the timestep
876 solver_graph.get_node_ptr_base("set_gpart_mass"),
877 solver_graph.get_node_ptr_base("set_scheduler_patchdata"),
878 solver_graph.get_node_ptr_base("attach_part_counts"),
879 solver_graph.get_node_ptr_base("attach_xyz"),
880 solver_graph.get_node_ptr_base("attach_vxyz"),
881 solver_graph.get_node_ptr_base("attach_axyz"),
882 // Actually perform the accretion
883 flag_node,
884 qty_node,
885 evict_node,
886 // free the refs since the particle counts may have changed
887 free_xyz,
888 free_vxyz,
889 free_axyz,
890 }));
891
892 // register the actual node that will be used
893 solver_graph.register_node(
894 "sink accretion",
895 OperationSequence(
896 "sink accretion",
897 {
898 set_has_sinks,
899 if_has_sinks,
900 }));
901 }
902
904 // sink ext force (pairwise self-gravity between sink particles)
906 {
907 solver_graph.register_edge("sink_ext_force_G", IDataEdge<Tscal>("G", "G"));
908 solver_graph.register_edge(
909 "sink_ext_force_epsilon", IDataEdge<Tscal>("epsilon_grav_sink", "\\epsilon"));
910
911 auto set_G = solver_graph.register_node(
912 "set_sink_ext_force_G", NodeSetEdge<IDataEdge<Tscal>>([&](IDataEdge<Tscal> &g_edge) {
913 g_edge.data = solver_config.get_constant_G();
914 }));
915 shambase::get_check_ref(set_G).set_edges(
916 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("sink_ext_force_G"));
917
918 auto set_epsilon = solver_graph.register_node(
919 "set_sink_ext_force_epsilon",
920 NodeSetEdge<IDataEdge<Tscal>>([&](IDataEdge<Tscal> &epsilon_edge) {
921 epsilon_edge.data = 1e-9;
922 }));
923 shambase::get_check_ref(set_epsilon)
924 .set_edges(solver_graph.get_edge_ptr<IDataEdge<Tscal>>("sink_ext_force_epsilon"));
925
926 auto reset_acc_ext = solver_graph.register_node(
927 "reset_sink_acc_ext",
928 NodeSetEdge<IDataEdgeSerializable<std::vector<Tvec>>>(
929 [](IDataEdgeSerializable<std::vector<Tvec>> &acc_ext) {
930 for (Tvec &a : acc_ext.data) {
931 a = Tvec{};
932 }
933 }));
934 shambase::get_check_ref(reset_acc_ext)
935 .set_edges(
936 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_acc_ext"));
937
938 auto self_gravity
939 = solver_graph.register_node("sink_self_gravity", modules::SinkSelfGravityHost<Tvec>{});
940 shambase::get_check_ref(self_gravity)
941 .set_edges(
942 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("sink_ext_force_G"),
943 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("sink_ext_force_epsilon"),
944 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_pos"),
945 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tscal>>>("sink_mass"),
946 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_acc_ext"));
947
948 auto ext_force_body = solver_graph.register_node(
949 "sink_ext_force_body",
950 OperationSequence(
951 "sink ext force body",
952 {
953 set_G,
954 set_epsilon,
955 reset_acc_ext,
956 self_gravity,
957 }));
958
959 // register the actual node that will be used, gated on the same "has_sinks" edge
960 // maintained by the "sink accretion" section above
961 auto sink_ext_force = solver_graph.register_node(
962 "sink ext force", OperationIf("sink ext force", ext_force_body));
963 shambase::get_check_ref(sink_ext_force)
964 .set_edges(solver_graph.get_edge_ptr<IDataEdge<bool>>("has_sinks"));
965 }
966
968 // sink predictor step (leapfrog kick-drift of the sink particles themselves)
970 {
971 solver_graph.register_edge(
972 "sink_predictor_dt_half", IDataEdge<Tscal>("dt_half", "\\frac{dt}{2}"));
973
974 auto sink_predictor_dt_to_half_dt = solver_graph.register_node(
975 "sink_predictor_dt_to_half_dt",
976 NodeMapEdge<IDataEdge<Tscal>, IDataEdge<Tscal>>{
977 [](const IDataEdge<Tscal> &dt, IDataEdge<Tscal> &half_dt) {
978 half_dt.data = dt.data / 2;
979 }});
980 shambase::get_check_ref(sink_predictor_dt_to_half_dt)
981 .set_edges(
982 sync_data.get_edge_ptr<IDataEdge<Tscal>>("dt"),
983 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("sink_predictor_dt_half"));
984
985 auto sink_predictor_vel_update = solver_graph.register_node(
986 "sink_predictor_vel_update", ForwardEulerHost2Deriv<Tvec, Tscal>{});
987 shambase::get_check_ref(sink_predictor_vel_update)
988 .set_edges(
989 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("sink_predictor_dt_half"),
990 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_acc_sph"),
991 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_acc_ext"),
992 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_vel"));
993
994 auto sink_predictor_pos_update = solver_graph.register_node(
995 "sink_predictor_pos_update", ForwardEulerHost<Tvec, Tscal>{});
996 shambase::get_check_ref(sink_predictor_pos_update)
997 .set_edges(
998 sync_data.get_edge_ptr<IDataEdge<Tscal>>("dt"),
999 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_vel"),
1000 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_pos"));
1001
1002 auto sink_predictor_body = solver_graph.register_node(
1003 "sink_predictor_body",
1004 OperationSequence(
1005 "sink predictor body",
1006 {
1007 // recompute the sink self-gravity at the current (pre-predictor) sink
1008 // positions before using it to kick the sink velocities
1009 solver_graph.get_node_ptr_base("sink ext force"),
1010 sink_predictor_dt_to_half_dt,
1011 sink_predictor_vel_update,
1012 sink_predictor_pos_update,
1013 }));
1014
1015 // register the actual node that will be used, gated on the same "has_sinks" edge
1016 // maintained by the "sink accretion" section above
1017 auto sink_predictor = solver_graph.register_node(
1018 "sink predictor", OperationIf("sink predictor", sink_predictor_body));
1019 shambase::get_check_ref(sink_predictor)
1020 .set_edges(solver_graph.get_edge_ptr<IDataEdge<bool>>("has_sinks"));
1021 }
1022
1024 // sink corrector step (leapfrog kick of the sink particles themselves)
1026 {
1027 auto sink_corrector_vel_update = solver_graph.register_node(
1028 "sink_corrector_vel_update", ForwardEulerHost2Deriv<Tvec, Tscal>{});
1029 shambase::get_check_ref(sink_corrector_vel_update)
1030 .set_edges(
1031 solver_graph.get_edge_ptr<IDataEdge<Tscal>>("dt_half"),
1032 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_acc_sph"),
1033 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_acc_ext"),
1034 sync_data.get_edge_ptr<IDataEdgeSerializable<std::vector<Tvec>>>("sink_vel"));
1035
1036 // register the actual node that will be used, gated on the same "has_sinks" edge
1037 // maintained by the "sink accretion" section above
1038 auto sink_corrector = solver_graph.register_node(
1039 "sink corrector", OperationIf("sink corrector", sink_corrector_vel_update));
1040 shambase::get_check_ref(sink_corrector)
1041 .set_edges(solver_graph.get_edge_ptr<IDataEdge<bool>>("has_sinks"));
1042 }
1043
1045 // external force (point mass) accretion
1047 {
1048 using EF_PointMass = typename Config::ExtForceConfig::PointMass;
1049 using EF_PN_PW = typename Config::ExtForceConfig::PN_PW;
1050 using EF_LenseThirring = typename Config::ExtForceConfig::LenseThirring;
1051
1052 // Collect the accretion centers and radii of every external force able to accrete.
1053 // This is re-done at every timestep, hence a lambda rather than a one shot list.
1054 auto collect_accretors = [&]() {
1055 std::vector<Tvec> positions{};
1056 std::vector<Tscal> radii{};
1057 for (auto &var_force : solver_config.ext_force_config.ext_forces) {
1058 if (EF_PointMass *ext_force = std::get_if<EF_PointMass>(&var_force.val)) {
1059 positions.push_back(ext_force->central_pos);
1060 radii.push_back(ext_force->Racc);
1061 } else if (EF_PN_PW *ext_force = std::get_if<EF_PN_PW>(&var_force.val)) {
1062 positions.push_back(ext_force->central_pos);
1063 radii.push_back(ext_force->Racc);
1064 } else if (
1065 EF_LenseThirring *ext_force = std::get_if<EF_LenseThirring>(&var_force.val)) {
1066 positions.push_back(ext_force->central_pos);
1067 radii.push_back(ext_force->Racc);
1068 }
1069 }
1070 return std::pair{std::move(positions), std::move(radii)};
1071 };
1072
1073 solver_graph.register_edge(
1074 "ext_force_accretion_pos",
1075 IDataEdge<std::vector<Tvec>>("ext_force_accretion_pos", "\\mathbf{r}_{\\rm acc, ext}"));
1076 solver_graph.register_edge(
1077 "ext_force_accretion_racc",
1078 IDataEdge<std::vector<Tscal>>("ext_force_accretion_racc", "R_{\\rm acc, ext}"));
1079 solver_graph.register_edge(
1080 "ext_force_accretion_table",
1081 Field<u32>(1, "ext_force_accretion_table", "\\mathrm{acc}"));
1082 solver_graph.register_edge(
1083 "has_ext_force_accretion",
1084 IDataEdge<bool>("has_ext_force_accretion", "\\rm has\\_ext\\_force\\_accretion"));
1085
1086 auto set_accretion_pos = solver_graph.register_node(
1087 "set_ext_force_accretion_pos",
1088 NodeSetEdge<IDataEdge<std::vector<Tvec>>>(
1089 [collect_accretors](IDataEdge<std::vector<Tvec>> &accretion_pos) {
1090 accretion_pos.data = std::get<0>(collect_accretors());
1091 }));
1092 shambase::get_check_ref(set_accretion_pos)
1093 .set_edges(
1094 solver_graph.get_edge_ptr<IDataEdge<std::vector<Tvec>>>("ext_force_accretion_pos"));
1095
1096 auto set_accretion_racc = solver_graph.register_node(
1097 "set_ext_force_accretion_racc",
1098 NodeSetEdge<IDataEdge<std::vector<Tscal>>>(
1099 [collect_accretors](IDataEdge<std::vector<Tscal>> &accretion_racc) {
1100 accretion_racc.data = std::get<1>(collect_accretors());
1101 }));
1102 shambase::get_check_ref(set_accretion_racc)
1103 .set_edges(solver_graph.get_edge_ptr<IDataEdge<std::vector<Tscal>>>(
1104 "ext_force_accretion_racc"));
1105
1106 auto set_has_accretion = solver_graph.register_node(
1107 "set_has_ext_force_accretion",
1108 NodeMapEdge<IDataEdge<std::vector<Tvec>>, IDataEdge<bool>>{
1109 [](const IDataEdge<std::vector<Tvec>> &accretion_pos,
1110 IDataEdge<bool> &has_ext_force_accretion) {
1111 has_ext_force_accretion.data = !accretion_pos.data.empty();
1112 }});
1113 shambase::get_check_ref(set_has_accretion)
1114 .set_edges(
1115 solver_graph.get_edge_ptr<IDataEdge<std::vector<Tvec>>>("ext_force_accretion_pos"),
1116 solver_graph.get_edge_ptr<IDataEdge<bool>>("has_ext_force_accretion"));
1117
1118 // reuse the sink accretion nodes, without the quantity accretion step since external
1119 // forces have no dynamical state to conserve onto
1120 auto flag_node = solver_graph.register_node(
1121 "ext_force_flag_accrete_hard", modules::SinkParticlesFlagAccreteHard<Tvec>{});
1122 shambase::get_check_ref(flag_node).set_edges(
1123 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
1124 solver_graph.get_edge_ptr<FieldRefs<Tvec>>("xyz"),
1125 solver_graph.get_edge_ptr<IDataEdge<std::vector<Tvec>>>("ext_force_accretion_pos"),
1126 solver_graph.get_edge_ptr<IDataEdge<std::vector<Tscal>>>("ext_force_accretion_racc"),
1127 solver_graph.get_edge_ptr<Field<u32>>("ext_force_accretion_table"));
1128
1129 auto evict_node = solver_graph.register_node(
1130 "ext_force_evict_accreted", modules::SinkParticlesEvictAccretedParticles<Tvec>{});
1131 shambase::get_check_ref(evict_node)
1132 .set_edges(
1133 solver_graph.get_edge_ptr<Indexes<u32>>("part_counts"),
1134 solver_graph.get_edge_ptr<Field<u32>>("ext_force_accretion_table"),
1135 solver_graph.get_edge_ptr<PatchDataLayerRefs>("scheduler_patchdata"));
1136
1137 auto accretion_body = solver_graph.register_node(
1138 "ext_force_accretion_body",
1139 OperationSequence(
1140 "ext force accretion",
1141 {
1142 // the "time_step" sequence (attach fields to scheduler, ...) has not run
1143 // yet at this stage of the timestep
1144 solver_graph.get_node_ptr_base("set_scheduler_patchdata"),
1145 solver_graph.get_node_ptr_base("attach_part_counts"),
1146 solver_graph.get_node_ptr_base("attach_xyz"),
1147 // Actually perform the accretion
1148 flag_node,
1149 evict_node,
1150 // free the refs since the particle counts may have changed
1151 solver_graph.get_node_ptr_base("free_xyz_refs"),
1152 }));
1153
1154 auto if_has_accretion = solver_graph.register_node(
1155 "if_has_ext_force_accretion",
1156 OperationIf("if_has_ext_force_accretion", accretion_body));
1157 shambase::get_check_ref(if_has_accretion)
1158 .set_edges(solver_graph.get_edge_ptr<IDataEdge<bool>>("has_ext_force_accretion"));
1159
1160 // register the actual node that will be used
1161 solver_graph.register_node(
1162 "point mass accretion",
1163 OperationSequence(
1164 "point mass accretion",
1165 {
1166 set_accretion_pos,
1167 set_accretion_racc,
1168 set_has_accretion,
1169 if_has_accretion,
1170 }));
1171 }
1172
1173 {
1174 std::vector<std::shared_ptr<shamrock::solvergraph::INode>> seq{};
1175
1176 seq.push_back(solver_graph.get_node_ptr_base("sink accretion"));
1177 seq.push_back(solver_graph.get_node_ptr_base("point mass accretion"));
1178 seq.push_back(solver_graph.get_node_ptr_base("sink predictor"));
1179 seq.push_back(solver_graph.get_node_ptr_base("dt_to_half_dt"));
1180 seq.push_back(solver_graph.get_node_ptr_base("set_gpart_mass"));
1181 seq.push_back(solver_graph.get_node_ptr_base("attach fields to scheduler"));
1182 seq.push_back(solver_graph.get_node_ptr_base("leapfrog predictor"));
1183 if (do_part_killing_step) {
1184 seq.push_back(solver_graph.get_node_ptr_base("part killing step"));
1185 }
1186 seq.push_back(solver_graph.get_node_ptr_base("sink ext force"));
1187
1188 storage.solver_sequence = solver_graph.register_node(
1189 "time_step", OperationSequence("time step", std::move(seq)));
1190 }
1191
1193 // self gravity sequence
1195 if (solver_config.self_grav_config.is_sg_on()) {
1196
1197 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
1198
1199 auto constant_G = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("", "");
1200
1201 shamrock::solvergraph::NodeSetEdge<shamrock::solvergraph::IDataEdge<Tscal>> set_constant_G(
1202 [&](shamrock::solvergraph::IDataEdge<Tscal> &constant_G) {
1203 constant_G.data = solver_config.get_constant_G();
1204 });
1205
1206 set_constant_G.set_edges(constant_G);
1207
1208 auto field_xyz = shamrock::solvergraph::FieldRefs<Tvec>::make_shared("", "");
1209
1210 shamrock::solvergraph::NodeSetEdge<shamrock::solvergraph::FieldRefs<Tvec>> set_field_xyz(
1211 [&, ixyz](shamrock::solvergraph::FieldRefs<Tvec> &field_xyz_edge) {
1212 shamrock::solvergraph::DDPatchDataFieldRef<Tvec> field_xyz_refs = {};
1213 scheduler().for_each_patchdata_nonempty(
1214 [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) {
1215 auto &field = pdat.get_field<Tvec>(ixyz);
1216 field_xyz_refs.add_obj(p.id_patch, std::ref(field));
1217 });
1218 field_xyz_edge.set_refs(field_xyz_refs);
1219 });
1220 set_field_xyz.set_edges(field_xyz);
1221
1222 const u32 iaxyz_ext = pdl.get_field_idx<Tvec>("axyz_ext");
1223
1224 auto field_axyz_ext = shamrock::solvergraph::FieldRefs<Tvec>::make_shared("", "");
1225
1226 shamrock::solvergraph::NodeSetEdge<shamrock::solvergraph::FieldRefs<Tvec>>
1227 set_field_axyz_ext(
1228 [&, iaxyz_ext](shamrock::solvergraph::FieldRefs<Tvec> &field_axyz_ext_edge) {
1229 shamrock::solvergraph::DDPatchDataFieldRef<Tvec> field_axyz_ext_refs = {};
1230 scheduler().for_each_patchdata_nonempty(
1231 [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) {
1232 auto &field = pdat.get_field<Tvec>(iaxyz_ext);
1233 field_axyz_ext_refs.add_obj(p.id_patch, std::ref(field));
1234 });
1235 field_axyz_ext_edge.set_refs(field_axyz_ext_refs);
1236 });
1237 set_field_axyz_ext.set_edges(field_axyz_ext);
1238
1239 auto sizes = shamrock::solvergraph::Indexes<u32>::make_shared("", "");
1240
1241 shamrock::solvergraph::NodeSetEdge<shamrock::solvergraph::Indexes<u32>> set_sizes(
1242 [&](shamrock::solvergraph::Indexes<u32> &sizes) {
1243 sizes.indexes = {};
1244 scheduler().for_each_patchdata_nonempty(
1245 [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) {
1246 sizes.indexes.add_obj(p.id_patch, pdat.get_obj_cnt());
1247 });
1248 });
1249 set_sizes.set_edges(sizes);
1250
1251 auto gpart_mass = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("", "");
1252
1253 shamrock::solvergraph::NodeSetEdge<shamrock::solvergraph::IDataEdge<Tscal>> set_gpart_mass(
1254 [&](shamrock::solvergraph::IDataEdge<Tscal> &gpart_mass) {
1255 gpart_mass.data = solver_config.gpart_mass;
1256 });
1257
1258 set_gpart_mass.set_edges(gpart_mass);
1259
1260 std::shared_ptr<shamrock::solvergraph::INode> sg_inode = build_self_gravity_node<Tvec>(
1261 solver_config.self_grav_config,
1262 sizes,
1263 gpart_mass,
1264 constant_G,
1265 field_xyz,
1266 field_axyz_ext);
1267
1268 solver_graph.register_node(
1269 "self gravity sequence",
1270 OperationSequence(
1271 "self gravity",
1272 {
1273 shambase::to_shared(std::move(set_gpart_mass)),
1274 shambase::to_shared(std::move(set_constant_G)),
1275 shambase::to_shared(std::move(set_field_xyz)),
1276 shambase::to_shared(std::move(set_field_axyz_ext)),
1277 shambase::to_shared(std::move(set_sizes)),
1278 sg_inode,
1279 }));
1280 }
1281}
1282
1283template<class Tvec, template<class> class Kern>
1284void shammodels::sph::Solver<Tvec, Kern>::vtk_do_dump(
1285 std::string filename, bool add_patch_world_id) {
1286
1287 modules::VTKDump(context, solver_config).do_dump(filename, add_patch_world_id);
1288}
1289
1291// Debug interface dump
1293
1294namespace shammodels::sph {
1295
1296 template<class Tvec>
1297 struct Debug_ph_dump {
1298 using Tscal = shambase::VecComponent<Tvec>;
1299
1300 u64 nobj;
1301 f64 gpart_mass;
1302
1303 sycl::buffer<Tvec> &buf_xyz;
1304 sycl::buffer<Tscal> &buf_hpart;
1305 sycl::buffer<Tvec> &buf_vxyz;
1306 };
1307
1308 template<class Tvec>
1309 void fill_blocks(PhantomDumpBlock &block, Debug_ph_dump<Tvec> &info) {
1310
1311 using Tscal = shambase::VecComponent<Tvec>;
1312 std::vector<Tvec> xyz = shamalgs::memory::buf_to_vec(info.buf_xyz, info.nobj);
1313
1314 u64 xid = block.get_ref_fort_real("x");
1315 u64 yid = block.get_ref_fort_real("y");
1316 u64 zid = block.get_ref_fort_real("z");
1317
1318 for (auto vec : xyz) {
1319 block.blocks_fort_real[xid].vals.push_back(vec.x());
1320 block.blocks_fort_real[yid].vals.push_back(vec.y());
1321 block.blocks_fort_real[zid].vals.push_back(vec.z());
1322 }
1323
1324 std::vector<Tscal> h = shamalgs::memory::buf_to_vec(info.buf_hpart, info.nobj);
1325 u64 hid = block.get_ref_f32("h");
1326 for (auto h_ : h) {
1327 block.blocks_f32[hid].vals.push_back(h_);
1328 }
1329
1330 std::vector<Tvec> vxyz = shamalgs::memory::buf_to_vec(info.buf_vxyz, info.nobj);
1331
1332 u64 vxid = block.get_ref_fort_real("vx");
1333 u64 vyid = block.get_ref_fort_real("vy");
1334 u64 vzid = block.get_ref_fort_real("vz");
1335
1336 for (auto vec : vxyz) {
1337 block.blocks_fort_real[vxid].vals.push_back(vec.x());
1338 block.blocks_fort_real[vyid].vals.push_back(vec.y());
1339 block.blocks_fort_real[vzid].vals.push_back(vec.z());
1340 }
1341
1342 block.tot_count = block.blocks_fort_real[xid].vals.size();
1343 }
1344
1345 template<class Tvec>
1346 shammodels::sph::PhantomDump make_interface_debug_phantom_dump(Debug_ph_dump<Tvec> info) {
1347
1348 using Tscal = shambase::VecComponent<Tvec>;
1349 PhantomDump dump;
1350
1351 dump.override_magic_number();
1352 dump.iversion = 1;
1353 dump.fileid = sham::format("{:100s}", "FT:Phantom Shamrock writer");
1354
1355 u32 Ntot = info.nobj;
1356 dump.table_header_fort_int.add("nparttot", Ntot);
1357 dump.table_header_fort_int.add("ntypes", 8);
1358 dump.table_header_fort_int.add("npartoftype", Ntot);
1359 dump.table_header_fort_int.add("npartoftype", 0);
1360 dump.table_header_fort_int.add("npartoftype", 0);
1361 dump.table_header_fort_int.add("npartoftype", 0);
1362 dump.table_header_fort_int.add("npartoftype", 0);
1363 dump.table_header_fort_int.add("npartoftype", 0);
1364 dump.table_header_fort_int.add("npartoftype", 0);
1365 dump.table_header_fort_int.add("npartoftype", 0);
1366
1367 dump.table_header_i64.add("nparttot", Ntot);
1368 dump.table_header_i64.add("ntypes", 8);
1369 dump.table_header_i64.add("npartoftype", Ntot);
1370 dump.table_header_i64.add("npartoftype", 0);
1371 dump.table_header_i64.add("npartoftype", 0);
1372 dump.table_header_i64.add("npartoftype", 0);
1373 dump.table_header_i64.add("npartoftype", 0);
1374 dump.table_header_i64.add("npartoftype", 0);
1375 dump.table_header_i64.add("npartoftype", 0);
1376 dump.table_header_i64.add("npartoftype", 0);
1377
1378 dump.table_header_fort_int.add("nblocks", 1);
1379 dump.table_header_fort_int.add("nptmass", 0);
1380 dump.table_header_fort_int.add("ndustlarge", 0);
1381 dump.table_header_fort_int.add("ndustsmall", 0);
1382 dump.table_header_fort_int.add("idust", 7);
1383 dump.table_header_fort_int.add("idtmax_n", 1);
1384 dump.table_header_fort_int.add("idtmax_frac", 0);
1385 dump.table_header_fort_int.add("idumpfile", 0);
1386 dump.table_header_fort_int.add("majorv", 2023);
1387 dump.table_header_fort_int.add("minorv", 0);
1388 dump.table_header_fort_int.add("microv", 0);
1389 dump.table_header_fort_int.add("isink", 0);
1390
1391 dump.table_header_i32.add("iexternalforce", 0);
1392 dump.table_header_i32.add("ieos", 2);
1393 dump.table_header_fort_real.add("gamma", 1.66667);
1394 dump.table_header_fort_real.add("RK2", 0);
1395 dump.table_header_fort_real.add("polyk2", 0);
1396 dump.table_header_fort_real.add("qfacdisc", 0.75);
1397 dump.table_header_fort_real.add("qfacdisc2", 0.75);
1398
1399 dump.table_header_fort_real.add("time", 0);
1400 dump.table_header_fort_real.add("dtmax", 0.1);
1401
1402 dump.table_header_fort_real.add("rhozero", 0);
1403 dump.table_header_fort_real.add("hfact", 1.2);
1404 dump.table_header_fort_real.add("tolh", 0.0001);
1405 dump.table_header_fort_real.add("C_cour", 0);
1406 dump.table_header_fort_real.add("C_force", 0);
1407 dump.table_header_fort_real.add("alpha", 0);
1408 dump.table_header_fort_real.add("alphau", 1);
1409 dump.table_header_fort_real.add("alphaB", 1);
1410
1411 dump.table_header_fort_real.add("massoftype", info.gpart_mass);
1412 dump.table_header_fort_real.add("massoftype", 0);
1413 dump.table_header_fort_real.add("massoftype", 0);
1414 dump.table_header_fort_real.add("massoftype", 0);
1415 dump.table_header_fort_real.add("massoftype", 0);
1416 dump.table_header_fort_real.add("massoftype", 0);
1417 dump.table_header_fort_real.add("massoftype", 0);
1418 dump.table_header_fort_real.add("massoftype", 0);
1419
1420 dump.table_header_fort_real.add("Bextx", 0);
1421 dump.table_header_fort_real.add("Bexty", 0);
1422 dump.table_header_fort_real.add("Bextz", 0);
1423 dump.table_header_fort_real.add("dum", 0);
1424
1425 dump.table_header_fort_real.add("get_conserv", -1);
1426 dump.table_header_fort_real.add("etot_in", 0.59762);
1427 dump.table_header_fort_real.add("angtot_in", 0.0189694);
1428 dump.table_header_fort_real.add("totmom_in", 0.0306284);
1429
1430 dump.table_header_f64.add("udist", 1);
1431 dump.table_header_f64.add("umass", 1);
1432 dump.table_header_f64.add("utime", 1);
1433 dump.table_header_f64.add("umagfd", 3.54491);
1434
1435 PhantomDumpBlock block_part;
1436
1437 fill_blocks(block_part, info);
1438
1439 dump.blocks.push_back(std::move(block_part));
1440
1441 return dump;
1442 }
1443
1444} // namespace shammodels::sph
1445
1446template<class Tvec, template<class> class Kern>
1447void shammodels::sph::Solver<Tvec, Kern>::gen_serial_patch_tree() {
1448 StackEntry stack_loc{};
1449
1450 SerialPatchTree<Tvec> _sptree = SerialPatchTree<Tvec>::build(scheduler());
1451 _sptree.attach_buf();
1452 storage.serial_patch_tree.set(std::move(_sptree));
1453}
1454
1460template<class Tvec, template<class> class Kern>
1461void shammodels::sph::Solver<Tvec, Kern>::apply_position_boundary(Tscal time_val) {
1462
1463 StackEntry stack_loc{};
1464
1465 shamlog_debug_ln("SphSolver", "apply position boundary");
1466
1467 PatchScheduler &sched = scheduler();
1468
1469 shamrock::SchedulerUtility integrators(sched);
1470 shamrock::ReattributeDataUtility reatrib(sched);
1471
1472 auto &pdl = sched.pdl_old();
1473
1474 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
1475 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
1476 auto [bmin, bmax] = sched.get_box_volume<Tvec>();
1477
1478 using SolverConfigBC = typename Config::BCConfig;
1479 using SolverBCFree = typename SolverConfigBC::Free;
1480 using SolverBCPeriodic = typename SolverConfigBC::Periodic;
1481 using SolverBCShearingPeriodic = typename SolverConfigBC::ShearingPeriodic;
1482 if (SolverBCFree *c = std::get_if<SolverBCFree>(&solver_config.boundary_config.config)) {
1483 if (shamcomm::world_rank() == 0) {
1484 logger::info_ln("PositionUpdated", "free boundaries skipping geometry update");
1485 }
1486 } else if (
1487 SolverBCPeriodic *c
1488 = std::get_if<SolverBCPeriodic>(&solver_config.boundary_config.config)) {
1489 integrators.fields_apply_periodicity(ixyz, std::pair{bmin, bmax});
1490 } else if (
1491 SolverBCShearingPeriodic *c
1492 = std::get_if<SolverBCShearingPeriodic>(&solver_config.boundary_config.config)) {
1493 integrators.fields_apply_shearing_periodicity(
1494 ixyz,
1495 ivxyz,
1496 std::pair{bmin, bmax},
1497 c->shear_base,
1498 c->shear_dir,
1499 c->shear_speed * time_val,
1500 c->shear_speed);
1501 }
1502
1503 reatrib.reatribute_patch_objects(storage.serial_patch_tree.get(), "xyz");
1504}
1505
1506template<class Tvec, template<class> class Kern>
1507void shammodels::sph::Solver<Tvec, Kern>::build_ghost_cache() {
1508
1509 StackEntry stack_loc{};
1510
1511 using SPHUtils = sph::SPHUtilities<Tvec, Kernel>;
1512 SPHUtils sph_utils(scheduler());
1513
1514 storage.ghost_patch_cache.set(sph_utils.build_interf_cache(
1515 storage.ghost_handler.get(),
1516 storage.serial_patch_tree.get(),
1517 solver_config.htol_up_coarse_cycle));
1518
1519 // storage.ghost_handler.get().gen_debug_patch_ghost(storage.ghost_patch_cache.get());
1520}
1521
1522template<class Tvec, template<class> class Kern>
1523void shammodels::sph::Solver<Tvec, Kern>::clear_ghost_cache() {
1524 StackEntry stack_loc{};
1525 storage.ghost_patch_cache.reset();
1526}
1527
1528template<class Tvec, template<class> class Kern>
1529void shammodels::sph::Solver<Tvec, Kern>::merge_position_ghost() {
1530
1531 StackEntry stack_loc{};
1532
1533 storage.merged_xyzh.set(storage.ghost_handler.get().build_comm_merge_positions(
1534 storage.ghost_patch_cache.get(),
1535 storage.exchange_gz_positions,
1536 solver_config.show_ghost_zone_graph));
1537
1538 { // set element counts
1539 shambase::get_check_ref(storage.part_counts).indexes
1540 = storage.merged_xyzh.get().template map<u32>(
1541 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
1542 return scheduler().patch_data.get_pdat(id).get_obj_cnt();
1543 });
1544 }
1545
1546 { // set element counts
1547 shambase::get_check_ref(storage.part_counts_with_ghost).indexes
1548 = storage.merged_xyzh.get().template map<u32>(
1549 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
1550 return mpdat.get_obj_cnt();
1551 });
1552 }
1553
1554 { // Attach spans to block coords
1555 shambase::get_check_ref(storage.positions_with_ghosts)
1556 .set_refs(storage.merged_xyzh.get()
1557 .template map<std::reference_wrapper<PatchDataField<Tvec>>>(
1558 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
1559 return std::ref(mpdat.get_field<Tvec>(0));
1560 }));
1561
1562 shambase::get_check_ref(storage.hpart_with_ghosts)
1563 .set_refs(storage.merged_xyzh.get()
1564 .template map<std::reference_wrapper<PatchDataField<Tscal>>>(
1565 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
1566 return std::ref(mpdat.get_field<Tscal>(1));
1567 }));
1568 }
1569}
1570
1571template<class Tvec, template<class> class Kern>
1575
1576template<class Tvec, template<class> class Kern>
1577void shammodels::sph::Solver<Tvec, Kern>::clear_merged_pos_trees() {
1578 StackEntry stack_loc{};
1579 storage.merged_pos_trees.reset();
1580}
1581
1582template<class Tvec, template<class> class Kern>
1583void shammodels::sph::Solver<Tvec, Kern>::sph_prestep(Tscal time_val, Tscal dt) {
1584 StackEntry stack_loc{};
1585
1586 using namespace shamrock;
1587 using namespace shamrock::patch;
1588
1589 using RTree = shamtree::CompressedLeafBVH<u_morton, Tvec, 3>;
1590 using SPHUtils = sph::SPHUtilities<Tvec, Kernel>;
1591
1592 SPHUtils sph_utils(scheduler());
1593 shamrock::SchedulerUtility utility(scheduler());
1594
1595 PatchDataLayerLayout &pdl = scheduler().pdl_old();
1596 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
1597
1598 ComputeField<Tscal> _epsilon_h, _h_old;
1599
1600 auto should_set_omega_mask = std::make_shared<shamrock::solvergraph::Field<u32>>(
1601 1, "should_set_omega_mask", "should_set_omega_mask");
1602
1603 u32 hstep_cnt = 0;
1604 u32 hstep_max = solver_config.h_max_subcycles_count;
1605 for (; hstep_cnt < hstep_max; hstep_cnt++) {
1606
1607 gen_ghost_handler(time_val + dt);
1608 build_ghost_cache();
1609 merge_position_ghost();
1610 build_merged_pos_trees();
1611 compute_presteps_rint();
1612 start_neighbors_cache();
1613
1614 _epsilon_h = utility.make_compute_field<Tscal>("epsilon_h", 1, Tscal(100));
1615 _h_old = utility.save_field<Tscal>(ihpart, "h_old");
1616
1617 Tscal max_eps_h;
1618
1619 if (solver_config.gpart_mass == 0) {
1620 shambase::throw_with_loc<std::runtime_error>(sham::format(
1621 "invalid gpart_mass {}, this configuration can not converge.\n"
1622 "Please set it using either model.set_particle_mass(pmass) or "
1623 "cfg.set_particle_mass(pmass)",
1624 solver_config.gpart_mass));
1625 }
1626
1627 // sizes
1628 std::shared_ptr<shamrock::solvergraph::Indexes<u32>> sizes
1629 = std::make_shared<shamrock::solvergraph::Indexes<u32>>("", "");
1630 scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
1631 sizes->indexes.add_obj(p.id_patch, pdat.get_obj_cnt());
1632 });
1633
1634 // neigh cache
1635 auto &neigh_cache = storage.neigh_cache;
1636
1637 // positions
1638 auto &pos_merged = storage.positions_with_ghosts;
1639
1640 // old smoothing length field
1641 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hold
1642 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("", "");
1643 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> hold_refs = {};
1644 scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
1645 auto &field = _h_old.get_field(p.id_patch);
1646 hold_refs.add_obj(p.id_patch, std::ref(field));
1647 });
1648 hold->set_refs(hold_refs);
1649
1650 // new smoothing length field
1651 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hnew
1652 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("", "");
1653 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> hnew_refs = {};
1654 scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
1655 auto &field = pdat.get_field<Tscal>(ihpart);
1656 hnew_refs.add_obj(p.id_patch, std::ref(field));
1657 });
1658 hnew->set_refs(hnew_refs);
1659
1660 // epsilon field
1661 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> eps_h
1662 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("", "");
1663 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> eps_h_refs = {};
1664 scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
1665 auto &field = _epsilon_h.get_field(p.id_patch);
1666 eps_h_refs.add_obj(p.id_patch, std::ref(field));
1667 });
1668 eps_h->set_refs(eps_h_refs);
1669
1670 std::shared_ptr<shamrock::solvergraph::INode> smth_h_iter_ptr;
1671
1672 using h_conf_density_based = typename SmoothingLengthConfig::DensityBased;
1673 using h_conf_neigh_lim = typename SmoothingLengthConfig::DensityBasedNeighLim;
1674
1675 if (h_conf_density_based *conf
1676 = std::get_if<h_conf_density_based>(&solver_config.smoothing_length_config.config)) {
1677 std::shared_ptr<shammodels::sph::modules::IterateSmoothingLengthDensity<Tvec, Kernel>>
1678 smth_h_iter = std::make_shared<
1679 shammodels::sph::modules::IterateSmoothingLengthDensity<Tvec, Kernel>>(
1680 solver_config.gpart_mass,
1681 solver_config.htol_up_coarse_cycle,
1682 solver_config.htol_up_fine_cycle,
1683 solver_config.epsilon_h);
1684 smth_h_iter->set_edges(sizes, neigh_cache, pos_merged, hold, hnew, eps_h);
1685 smth_h_iter_ptr = smth_h_iter;
1686 } else if (
1687 h_conf_neigh_lim *conf
1688 = std::get_if<h_conf_neigh_lim>(&solver_config.smoothing_length_config.config)) {
1689 std::shared_ptr<
1690 shammodels::sph::modules::IterateSmoothingLengthDensityNeighLim<Tvec, Kernel>>
1691 smth_h_iter_neigh_lim = std::make_shared<
1692 shammodels::sph::modules::IterateSmoothingLengthDensityNeighLim<Tvec, Kernel>>(
1693 solver_config.gpart_mass,
1694 solver_config.htol_up_coarse_cycle,
1695 solver_config.htol_up_fine_cycle,
1696 conf->max_neigh_count,
1697 solver_config.epsilon_h);
1698 smth_h_iter_neigh_lim->set_edges(
1699 sizes, neigh_cache, pos_merged, hold, hnew, eps_h, should_set_omega_mask);
1700 smth_h_iter_ptr = smth_h_iter_neigh_lim;
1701 } else {
1702 shambase::throw_with_loc<std::runtime_error>("Invalid smoothing length configuration");
1703 }
1704 // iterate smoothing length
1705
1706 std::shared_ptr<shamrock::solvergraph::IDataEdge<bool>> is_converged
1707 = shamrock::solvergraph::IDataEdge<bool>::make_shared("", "");
1708
1709 shammodels::sph::modules::LoopSmoothingLengthIter<Tvec> loop_smth_h_iter(
1710 smth_h_iter_ptr, solver_config.epsilon_h, solver_config.h_iter_per_subcycles, false);
1711 loop_smth_h_iter.set_edges(eps_h, is_converged);
1712
1713 loop_smth_h_iter.evaluate();
1714
1715 if (!is_converged->data) {
1716
1717 Tscal largest_h = 0;
1718
1719 scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
1720 largest_h = sham::max(largest_h, pdat.get_field<Tscal>(ihpart).compute_max());
1721 });
1722 Tscal global_largest_h = shamalgs::collective::allreduce_max(largest_h);
1723
1724 std::string add_info = "";
1725 u64 cnt_unconverged = 0;
1726 scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
1727 auto res
1728 = _epsilon_h.get_field(p.id_patch).get_ids_buf_where([](auto access, u32 id) {
1729 return access[id] == -1;
1730 });
1731
1732 if (hstep_cnt == hstep_max - 1) {
1733 if (std::get<0>(res)) {
1734 add_info += "\n patch " + std::to_string(p.id_patch) + " ";
1735 add_info += "errored parts : \n";
1736 sycl::buffer<u32> &idx_err = *std::get<0>(res);
1737
1738 sham::DeviceBuffer<Tvec> &xyz = pdat.get_field_buf_ref<Tvec>(0);
1739 sham::DeviceBuffer<Tscal> &hpart = pdat.get_field_buf_ref<Tscal>(ihpart);
1740
1741 auto pos = xyz.copy_to_stdvec();
1742 auto h = hpart.copy_to_stdvec();
1743
1744 {
1745 sycl::host_accessor acc{idx_err};
1746 for (u32 i = 0; i < idx_err.size(); i++) {
1747 add_info += sham::format(
1748 "{} - pos : {}, hpart : {}\n", acc[i], pos[acc[i]], h[acc[i]]);
1749 }
1750 }
1751 }
1752 }
1753
1754 cnt_unconverged += std::get<1>(res);
1755 });
1756
1757 u64 global_cnt_unconverged = shamalgs::collective::allreduce_sum(cnt_unconverged);
1758
1759 if (shamcomm::world_rank() == 0) {
1760 logger::warn_ln(
1761 "Smoothinglength",
1762 "smoothing length is not converged, rerunning the iterator ...\n largest h "
1763 "=",
1764 global_largest_h,
1765 "unconverged cnt =",
1766 global_cnt_unconverged,
1767 add_info);
1768 }
1769
1770 reset_ghost_handler();
1771 clear_ghost_cache();
1772
1773 shambase::get_check_ref(storage.part_counts).free_alloc();
1774 shambase::get_check_ref(storage.part_counts_with_ghost).free_alloc();
1775 shambase::get_check_ref(storage.positions_with_ghosts).free_alloc();
1776 shambase::get_check_ref(storage.hpart_with_ghosts).free_alloc();
1777
1778 storage.merged_xyzh.reset();
1779
1780 clear_merged_pos_trees();
1781 reset_presteps_rint();
1782 reset_neighbors_cache();
1783
1784 // scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchData &pdat) {
1785 // pdat.synchronize_buf();
1786 // });
1787
1788 continue;
1789 }
1790
1791 // The hpart is not valid anymore in ghost zones since we iterated it's value
1792 shambase::get_check_ref(storage.hpart_with_ghosts).free_alloc();
1793
1794 _epsilon_h.reset();
1795 _h_old.reset();
1796 break;
1797 }
1798
1799 if (hstep_cnt == hstep_max) {
1800 logger::err_ln("SPH", "the h iterator is not converged after", hstep_cnt, "iterations");
1801 }
1802
1803 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hnew_edge
1804 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("", "");
1805 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> hnew_refs = {};
1806 scheduler().for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
1807 auto &field = pdat.get_field<Tscal>(ihpart);
1808 hnew_refs.add_obj(p.id_patch, std::ref(field));
1809 });
1810 hnew_edge->set_refs(hnew_refs);
1811
1812 modules::NodeComputeOmega<Tvec, Kern> compute_omega{solver_config.gpart_mass};
1813 compute_omega.set_edges(
1814 storage.part_counts,
1815 storage.neigh_cache,
1816 storage.positions_with_ghosts,
1817 hnew_edge,
1818 storage.omega);
1819 compute_omega.evaluate();
1820
1821 if (solver_config.smoothing_length_config.is_density_based_neigh_lim()) {
1822 // if the h limiter is triggered, omega does not hold it's sense of dh/dr anymore
1823 // so we set it to 1, this effectively is equivalent of disabling the energy correction
1824 // term corresponding to dh/dr
1825 modules::SetWhenMask<Tscal> set_omega_mask{1};
1826 set_omega_mask.set_edges(storage.part_counts, should_set_omega_mask, storage.omega);
1827 set_omega_mask.evaluate();
1828 }
1829}
1830
1831template<class Tvec, template<class> class Kern>
1832void shammodels::sph::Solver<Tvec, Kern>::init_ghost_layout() {
1833
1834 storage.ghost_layout = std::make_shared<shamrock::patch::PatchDataLayerLayout>();
1835
1836 shamrock::patch::PatchDataLayerLayout &ghost_layout
1837 = shambase::get_check_ref(storage.ghost_layout);
1838
1839 solver_config.set_ghost_layout(ghost_layout);
1840
1841 storage.xyzh_ghost_layout = std::make_shared<shamrock::patch::PatchDataLayerLayout>();
1842 storage.xyzh_ghost_layout->template add_field<Tvec>("xyz", 1);
1843 storage.xyzh_ghost_layout->template add_field<Tscal>("hpart", 1);
1844}
1845
1846template<class Tvec, template<class> class Kern>
1847void shammodels::sph::Solver<Tvec, Kern>::compute_presteps_rint() {
1848
1849 StackEntry stack_loc{};
1850
1851 auto &xyzh_merged = storage.merged_xyzh.get();
1852 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
1853
1854 storage.rtree_rint_field.set(
1855 storage.merged_pos_trees.get().template map<shamtree::KarrasRadixTreeField<Tscal>>(
1856 [&](u64 id, RTree &rtree) -> shamtree::KarrasRadixTreeField<Tscal> {
1857 shamrock::patch::PatchDataLayer &tmp = xyzh_merged.get(id);
1858 auto &buf = tmp.get_field_buf_ref<Tscal>(1);
1859 auto buf_int = shamtree::new_empty_karras_radix_tree_field<Tscal>();
1860
1861 auto ret = shamtree::compute_tree_field_max_field<Tscal>(
1862 rtree.structure,
1863 rtree.reduced_morton_set.get_leaf_cell_iterator(),
1864 std::move(buf_int),
1865 buf);
1866
1867 // the old tree used to increase the size of the hmax of the tree nodes by the
1868 // tolerance so we do it also with the new tree, maybe we should move that somewhere
1869 // else.
1870 sham::kernel_call(
1871 dev_sched->get_queue(),
1872 sham::MultiRef{},
1873 sham::MultiRef{ret.buf_field},
1874 ret.buf_field.get_size(),
1875 [htol = solver_config.htol_up_coarse_cycle](u32 i, Tscal *h_tree) {
1876 h_tree[i] *= htol;
1877 });
1878
1879 return std::move(ret);
1880 }));
1881}
1882
1883template<class Tvec, template<class> class Kern>
1884void shammodels::sph::Solver<Tvec, Kern>::reset_presteps_rint() {
1885 storage.rtree_rint_field.reset();
1886}
1887
1888template<class Tvec, template<class> class Kern>
1889void shammodels::sph::Solver<Tvec, Kern>::start_neighbors_cache() {
1890 shammodels::sph::modules::NeighbourCache<Tvec, u_morton, Kern> neigh_cache_builder(
1891 context, solver_config, storage);
1892
1893 switch (solver_config.neigh_cache_strategy) {
1894 case NeighCacheStrategy::SingleStage: neigh_cache_builder.start_neighbors_cache(); break;
1895 case NeighCacheStrategy::TwoStage: neigh_cache_builder.start_neighbors_cache_2stages(); break;
1896 default: shambase::throw_unimplemented("unknown neighbours cache strategy");
1897 }
1898
1899 if (solver_config.show_neigh_stats) {
1900 auto &pos_merged = storage.positions_with_ghosts;
1901 auto &neigh_cache = storage.neigh_cache;
1902 auto &hpart_with_ghosts = storage.hpart_with_ghosts;
1903 auto &part_counts = storage.part_counts;
1904
1905 modules::ComputeNeighStats<Tvec> compute_neigh_stats(Kernel::Rkern);
1906
1907 compute_neigh_stats.set_edges(part_counts, neigh_cache, pos_merged, hpart_with_ghosts);
1908 compute_neigh_stats.evaluate();
1909 }
1910}
1911
1912template<class Tvec, template<class> class Kern>
1913void shammodels::sph::Solver<Tvec, Kern>::reset_neighbors_cache() {
1914 // storage.neighbors_cache.reset();
1915}
1916
1917template<class Tvec, template<class> class Kern>
1918void shammodels::sph::Solver<Tvec, Kern>::communicate_merge_ghosts_fields() {
1919
1920 StackEntry stack_loc{};
1921
1922 shambase::Timer timer_interf;
1923 timer_interf.start();
1924
1925 using namespace shamrock;
1926 using namespace shamrock::patch;
1927
1928 bool has_alphaAV_field = solver_config.has_field_alphaAV();
1929 bool has_soundspeed_field = solver_config.ghost_has_soundspeed();
1930
1931 bool has_B_field = solver_config.has_field_B_on_rho();
1932 bool has_psi_field = solver_config.has_field_psi_on_ch();
1933 bool has_curlB_field = solver_config.has_field_curlB();
1934 bool has_epsilon_field = solver_config.dust_config.has_epsilon_field();
1935 bool has_deltav_field = solver_config.dust_config.has_deltav_field();
1936 bool has_s_j_field = solver_config.dust_config.has_s_j_field();
1937
1938 PatchDataLayerLayout &pdl = scheduler().pdl_old();
1939 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
1940 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
1941 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
1942 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
1943 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
1944 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
1945
1946 const u32 ialpha_AV = (has_alphaAV_field) ? pdl.get_field_idx<Tscal>("alpha_AV") : 0;
1947 const u32 isoundspeed = (has_soundspeed_field) ? pdl.get_field_idx<Tscal>("soundspeed") : 0;
1948
1949 const u32 iB_on_rho = (has_B_field) ? pdl.get_field_idx<Tvec>("B/rho") : 0;
1950 const u32 idB_on_rho = (has_B_field) ? pdl.get_field_idx<Tvec>("dB/rho") : 0;
1951 const u32 ipsi_on_ch = (has_psi_field) ? pdl.get_field_idx<Tscal>("psi/ch") : 0;
1952 const u32 idpsi_on_ch = (has_psi_field) ? pdl.get_field_idx<Tscal>("dpsi/ch") : 0;
1953 const u32 icurlB = (has_curlB_field) ? pdl.get_field_idx<Tvec>("curlB") : 0;
1954
1955 bool do_MHD_debug = solver_config.do_MHD_debug();
1956 const u32 imag_pressure = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("mag_pressure") : -1;
1957 const u32 imag_tension = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("mag_tension") : -1;
1958 const u32 igas_pressure = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("gas_pressure") : -1;
1959 const u32 itensile_corr = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("tensile_corr") : -1;
1960 const u32 ipsi_propag = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("psi_propag") : -1;
1961 const u32 ipsi_diff = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("psi_diff") : -1;
1962 const u32 ipsi_cons = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("psi_cons") : -1;
1963 const u32 iu_mhd = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("u_mhd") : -1;
1964
1965 const u32 iepsilon = (has_epsilon_field) ? pdl.get_field_idx<Tscal>("epsilon") : 0;
1966 const u32 ideltav = (has_deltav_field) ? pdl.get_field_idx<Tvec>("deltav") : 0;
1967 const u32 is_j = (has_s_j_field) ? pdl.get_field_idx<Tscal>("s_j") : 0;
1968
1969 auto &ghost_layout_ptr = storage.ghost_layout;
1970 shamrock::patch::PatchDataLayerLayout &ghost_layout = shambase::get_check_ref(ghost_layout_ptr);
1971 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
1972 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
1973 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
1974 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
1975
1976 const u32 iaxyz_interf
1977 = (solver_config.has_axyz_in_ghost()) ? ghost_layout.get_field_idx<Tvec>("axyz") : 0;
1978
1979 const u32 isoundspeed_interf
1980 = (has_soundspeed_field) ? ghost_layout.get_field_idx<Tscal>("soundspeed") : 0;
1981
1982 const u32 iB_interf = (has_B_field) ? ghost_layout.get_field_idx<Tvec>("B/rho") : 0;
1983 const u32 ipsi_interf = (has_psi_field) ? ghost_layout.get_field_idx<Tscal>("psi/ch") : 0;
1984 const u32 icurlB_interf = (has_curlB_field) ? ghost_layout.get_field_idx<Tvec>("curlB") : 0;
1985
1986 const u32 iepsilon_interf
1987 = (has_epsilon_field) ? ghost_layout.get_field_idx<Tscal>("epsilon") : 0;
1988 const u32 ideltav_interf = (has_deltav_field) ? ghost_layout.get_field_idx<Tvec>("deltav") : 0;
1989 const u32 is_j_interf = (has_s_j_field) ? ghost_layout.get_field_idx<Tscal>("s_j") : 0;
1990
1991 using InterfaceBuildInfos = typename sph::BasicSPHGhostHandler<Tvec>::InterfaceBuildInfos;
1992
1993 sph::BasicSPHGhostHandler<Tvec> &ghost_handle = storage.ghost_handler.get();
1994 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
1995
1996 auto pdat_interf = ghost_handle.template build_interface_native<PatchDataLayer>(
1997 storage.ghost_patch_cache.get(),
1998 [&](u64 sender, u64, InterfaceBuildInfos binfo, sham::DeviceBuffer<u32> &buf_idx, u32 cnt) {
1999 PatchDataLayer pdat(ghost_layout_ptr);
2000
2001 pdat.reserve(cnt);
2002
2003 return pdat;
2004 });
2005
2006 ghost_handle.template modify_interface_native<PatchDataLayer>(
2007 storage.ghost_patch_cache.get(),
2008 pdat_interf,
2009 [&](u64 sender,
2010 u64,
2011 InterfaceBuildInfos binfo,
2012 sham::DeviceBuffer<u32> &buf_idx,
2013 u32 cnt,
2014 PatchDataLayer &pdat) {
2015 PatchDataLayer &sender_patch = scheduler().patch_data.get_pdat(sender);
2016 PatchDataField<Tscal> &sender_omega = omega.get(sender);
2017
2018 sender_patch.get_field<Tscal>(ihpart).append_subset_to(
2019 buf_idx, cnt, pdat.get_field<Tscal>(ihpart_interf));
2020 sender_patch.get_field<Tscal>(iuint).append_subset_to(
2021 buf_idx, cnt, pdat.get_field<Tscal>(iuint_interf));
2022
2023 if (solver_config.has_axyz_in_ghost()) {
2024 sender_patch.get_field<Tvec>(iaxyz).append_subset_to(
2025 buf_idx, cnt, pdat.get_field<Tvec>(iaxyz_interf));
2026 }
2027
2028 sender_patch.get_field<Tvec>(ivxyz).append_subset_to(
2029 buf_idx, cnt, pdat.get_field<Tvec>(ivxyz_interf));
2030
2031 sender_omega.append_subset_to(buf_idx, cnt, pdat.get_field<Tscal>(iomega_interf));
2032
2033 if (has_soundspeed_field) {
2034 sender_patch.get_field<Tscal>(isoundspeed)
2035 .append_subset_to(buf_idx, cnt, pdat.get_field<Tscal>(isoundspeed_interf));
2036 }
2037
2038 if (has_B_field) {
2039 sender_patch.get_field<Tvec>(iB_on_rho).append_subset_to(
2040 buf_idx, cnt, pdat.get_field<Tvec>(iB_interf));
2041 }
2042
2043 if (has_psi_field) {
2044 sender_patch.get_field<Tscal>(ipsi_on_ch)
2045 .append_subset_to(buf_idx, cnt, pdat.get_field<Tscal>(ipsi_interf));
2046 }
2047
2048 if (has_curlB_field) {
2049 sender_patch.get_field<Tvec>(icurlB).append_subset_to(
2050 buf_idx, cnt, pdat.get_field<Tvec>(icurlB_interf));
2051 }
2052
2053 if (has_epsilon_field) {
2054 sender_patch.get_field<Tscal>(iepsilon).append_subset_to(
2055 buf_idx, cnt, pdat.get_field<Tscal>(iepsilon_interf));
2056 }
2057
2058 if (has_deltav_field) {
2059 sender_patch.get_field<Tvec>(ideltav).append_subset_to(
2060 buf_idx, cnt, pdat.get_field<Tvec>(ideltav_interf));
2061 }
2062
2063 if (has_s_j_field) {
2064 sender_patch.get_field<Tscal>(is_j).append_subset_to(
2065 buf_idx, cnt, pdat.get_field<Tscal>(is_j_interf));
2066 }
2067 });
2068
2069 ghost_handle.template modify_interface_native<PatchDataLayer>(
2070 storage.ghost_patch_cache.get(),
2071 pdat_interf,
2072 [&](u64 sender,
2073 u64,
2074 InterfaceBuildInfos binfo,
2075 sham::DeviceBuffer<u32> &buf_idx,
2076 u32 cnt,
2077 PatchDataLayer &pdat) {
2078 if (sycl::length(binfo.offset_speed) > 0) {
2079 pdat.get_field<Tvec>(ivxyz_interf).apply_offset(binfo.offset_speed);
2080 }
2081 });
2082
2083 shambase::DistributedDataShared<PatchDataLayer> interf_pdat = ghost_handle.communicate_pdat(
2084 ghost_layout_ptr,
2085 std::move(pdat_interf),
2086 storage.exchange_gz_node,
2087 solver_config.show_ghost_zone_graph);
2088
2089 std::map<u64, u64> sz_interf_map;
2090 interf_pdat.for_each([&](u64 s, u64 r, PatchDataLayer &pdat_interf) {
2091 sz_interf_map[r] += pdat_interf.get_obj_cnt();
2092 });
2093
2094 storage.merged_patchdata_ghost.set(
2095 ghost_handle.template merge_native<PatchDataLayer, PatchDataLayer>(
2096 std::move(interf_pdat),
2097 [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) {
2098 PatchDataLayer pdat_new(ghost_layout_ptr);
2099
2100 u32 or_elem = pdat.get_obj_cnt();
2101 pdat_new.reserve(or_elem + sz_interf_map[p.id_patch]);
2102 u32 total_elements = or_elem;
2103
2104 PatchDataField<Tscal> &cur_omega = omega.get(p.id_patch);
2105
2106 pdat_new.get_field<Tscal>(ihpart_interf).insert(pdat.get_field<Tscal>(ihpart));
2107 pdat_new.get_field<Tscal>(iuint_interf).insert(pdat.get_field<Tscal>(iuint));
2108 pdat_new.get_field<Tvec>(ivxyz_interf).insert(pdat.get_field<Tvec>(ivxyz));
2109
2110 if (solver_config.has_axyz_in_ghost()) {
2111 pdat_new.get_field<Tvec>(iaxyz_interf).insert(pdat.get_field<Tvec>(iaxyz));
2112 }
2113
2114 pdat_new.get_field<Tscal>(iomega_interf).insert(cur_omega);
2115
2116 if (has_soundspeed_field) {
2117 pdat_new.get_field<Tscal>(isoundspeed_interf)
2118 .insert(pdat.get_field<Tscal>(isoundspeed));
2119 }
2120
2121 if (has_B_field) {
2122 pdat_new.get_field<Tvec>(iB_interf).insert(pdat.get_field<Tvec>(iB_on_rho));
2123 }
2124
2125 if (has_psi_field) {
2126 pdat_new.get_field<Tscal>(ipsi_interf)
2127 .insert(pdat.get_field<Tscal>(ipsi_on_ch));
2128 }
2129
2130 if (has_curlB_field) {
2131 pdat_new.get_field<Tvec>(icurlB_interf).insert(pdat.get_field<Tvec>(icurlB));
2132 }
2133
2134 if (has_epsilon_field) {
2135 pdat_new.get_field<Tscal>(iepsilon_interf)
2136 .insert(pdat.get_field<Tscal>(iepsilon));
2137 }
2138
2139 if (has_deltav_field) {
2140 pdat_new.get_field<Tvec>(ideltav_interf).insert(pdat.get_field<Tvec>(ideltav));
2141 }
2142
2143 if (has_s_j_field) {
2144 pdat_new.get_field<Tscal>(is_j_interf).insert(pdat.get_field<Tscal>(is_j));
2145 }
2146
2147 pdat_new.check_field_obj_cnt_match();
2148
2149 return pdat_new;
2150 },
2151 [](PatchDataLayer &pdat, PatchDataLayer &pdat_interf) {
2152 pdat.insert_elements(pdat_interf);
2153 }));
2154
2155 timer_interf.stop();
2156 storage.timings_details.interface += timer_interf.elapsed_sec();
2157}
2158
2159template<class Tvec, template<class> class Kern>
2160void shammodels::sph::Solver<Tvec, Kern>::reset_merge_ghosts_fields() {
2161 storage.merged_patchdata_ghost.reset();
2162}
2163
2165// start artificial viscosity section //////////////////////////////////////////////////////////////
2167
2168template<class Tvec, template<class> class Kern>
2169void shammodels::sph::Solver<Tvec, Kern>::update_artificial_viscosity(Tscal dt) {
2170
2171 sph::modules::UpdateViscosity<Tvec, Kern>(context, solver_config, storage)
2172 .update_artificial_viscosity(dt);
2173}
2174
2176// end artificial viscosity section ////////////////////////////////////////////////////////////////
2178
2179template<class Tvec, template<class> class Kern>
2184
2185template<class Tvec, template<class> class Kern>
2186void shammodels::sph::Solver<Tvec, Kern>::reset_eos_fields() {
2187 shambase::get_check_ref(storage.pressure).free_alloc();
2188 shambase::get_check_ref(storage.soundspeed).free_alloc();
2189}
2190
2191template<class Tvec, template<class> class Kern>
2192void shammodels::sph::Solver<Tvec, Kern>::prepare_corrector() {
2193
2194 StackEntry stack_loc{};
2195
2196 using namespace shamrock;
2197 using namespace shamrock::patch;
2198 shamrock::SchedulerUtility utility(scheduler());
2199 PatchDataLayerLayout &pdl = scheduler().pdl_old();
2200
2201 bool has_B_field = solver_config.has_field_B_on_rho();
2202 bool has_psi_field = solver_config.has_field_psi_on_ch();
2203 bool has_epsilon_field = solver_config.dust_config.has_epsilon_field();
2204 bool has_deltav_field = solver_config.dust_config.has_deltav_field();
2205 bool has_s_j_field = solver_config.dust_config.has_s_j_field();
2206
2207 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
2208 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
2209 const u32 idB_on_rho = (has_B_field) ? pdl.get_field_idx<Tvec>("dB/rho") : 0;
2210 const u32 idpsi_on_ch = (has_psi_field) ? pdl.get_field_idx<Tscal>("dpsi/ch") : 0;
2211
2212 shamlog_debug_ln("sph::BasicGas", "save old fields");
2213 storage.old_axyz.set(utility.save_field<Tvec>(iaxyz, "axyz_old"));
2214 storage.old_duint.set(utility.save_field<Tscal>(iduint, "duint_old"));
2215
2216 if (has_B_field) {
2217 storage.old_dB_on_rho.set(utility.save_field<Tvec>(idB_on_rho, "dB/rho_old"));
2218 }
2219 if (has_psi_field) {
2220 storage.old_dpsi_on_ch.set(utility.save_field<Tscal>(idpsi_on_ch, "dpsi/ch_old"));
2221 }
2222 if (has_epsilon_field) {
2223 storage.old_dtepsilon.set(
2224 utility.save_field<Tscal>(pdl.get_field_idx<Tscal>("dtepsilon"), "dtepsilon_old"));
2225 }
2226 if (has_deltav_field) {
2227 storage.old_dtdeltav.set(
2228 utility.save_field<Tvec>(pdl.get_field_idx<Tvec>("dtdeltav"), "dtdeltav_old"));
2229 }
2230 if (has_s_j_field) {
2231 storage.old_ds_j_dt.set(
2232 utility.save_field<Tscal>(pdl.get_field_idx<Tscal>("ds_j_dt"), "ds_j_dt_old"));
2233 }
2234}
2235
2236template<class T>
2237void map_field_refs(
2238 PatchScheduler &sched, u32 field_idx, shamrock::solvergraph::FieldRefs<T> &refs) {
2239
2240 using namespace shamrock::solvergraph;
2241 using namespace shamrock::patch;
2242
2243 shamrock::solvergraph::DDPatchDataFieldRef<T> field_refs = {};
2244 sched.for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
2245 auto &field = pdat.get_field<T>(field_idx);
2246 field_refs.add_obj(p.id_patch, std::ref(field));
2247 });
2248 refs.set_refs(field_refs);
2249}
2250
2251template<class T>
2252void map_field_refs_ext(
2253 PatchScheduler &sched,
2254 shambase::DistributedData<shamrock::patch::PatchDataLayer> &mpdats,
2255 u32 field_idx,
2256 shamrock::solvergraph::FieldRefs<T> &refs) {
2257
2258 using namespace shamrock::solvergraph;
2259 using namespace shamrock::patch;
2260
2261 shamrock::solvergraph::DDPatchDataFieldRef<T> field_refs = {};
2262 sched.for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
2263 PatchDataLayer &mpdat = mpdats.get(p.id_patch);
2264 auto &field = mpdat.get_field<T>(field_idx);
2265 field_refs.add_obj(p.id_patch, std::ref(field));
2266 });
2267 refs.set_refs(field_refs);
2268}
2269
2270template<class T>
2271void map_field_refs_ext(
2272 PatchScheduler &sched,
2273 shamrock::ComputeField<T> &field_data,
2274 shamrock::solvergraph::FieldRefs<T> &refs) {
2275
2276 using namespace shamrock::solvergraph;
2277 using namespace shamrock::patch;
2278
2279 shamrock::solvergraph::DDPatchDataFieldRef<T> field_refs = {};
2280 sched.for_each_patchdata_nonempty([&](const Patch p, PatchDataLayer &pdat) {
2281 auto &field = field_data.get_field(p.id_patch);
2282 field_refs.add_obj(p.id_patch, std::ref(field));
2283 });
2284 refs.set_refs(field_refs);
2285}
2286
2287template<class Tvec, template<class> class Kern>
2288void shammodels::sph::Solver<Tvec, Kern>::update_derivs(Tscal dt_hydro) {
2289
2290 // if one fluid is enabled time to compute the stopping times
2291 if (solver_config.dust_config.has_s_j_field()) {
2292
2293 auto &cfg = solver_config.dust_config;
2294 u32 ndust = cfg.get_dust_nvar();
2295
2296 using DustConfig = typename Config::DustConfig;
2297
2298 using None = typename DustConfig::None;
2299 using ConstantStoppingTimes = typename DustConfig::ConstantStoppingTimes;
2300 using EpsteinDrag = typename DustConfig::EpsteinDrag;
2301
2302 shamrock::patch::PatchDataLayerLayout &ghost_layout
2303 = shambase::get_check_ref(storage.ghost_layout.get());
2304 shamrock::patch::PatchDataLayerLayout &pdl = scheduler().pdl_old();
2305
2306 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
2307
2308 auto &part_counts_with_ghost = storage.part_counts_with_ghost;
2309 auto &part_counts = storage.part_counts;
2310
2311 shambase::DistributedData<shamrock::patch::PatchDataLayer> &mpdats
2312 = storage.merged_patchdata_ghost.get();
2313
2314 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hpart_refs
2315 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("hpart", "h");
2316 { // if was just reset before this call
2317 shambase::get_check_ref(hpart_refs)
2318 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
2319 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
2320 return std::ref(mpdat.get_field<Tscal>(ihpart_interf));
2321 }));
2322 }
2323
2324 auto gpart_mass
2325 = storage.solver_graph.template get_edge_ptr<shamrock::solvergraph::IDataEdge<Tscal>>(
2326 "gpart_mass");
2327 auto t_j_field
2328 = storage.solver_graph.template get_edge_ptr<shamrock::solvergraph::Field<Tscal>>(
2329 "Ts_j");
2330
2331 if (std::holds_alternative<None>(cfg.dust_drag_mode)) {
2332
2333 throw "bro WTF";
2334
2335 } else if (
2336 ConstantStoppingTimes *cfg_drag
2337 = std::get_if<ConstantStoppingTimes>(&cfg.dust_drag_mode)) {
2338
2339 std::shared_ptr<shamrock::solvergraph::IDataEdge<std::vector<Tscal>>> input_t_j
2340 = shamrock::solvergraph::IDataEdge<std::vector<Tscal>>::make_shared("", "");
2341 input_t_j->data = cfg_drag->stopping_times;
2342
2343 std::shared_ptr<modules::SetDustStoppingTimeConstant<Tvec>> node_set_tj
2344 = std::make_shared<modules::SetDustStoppingTimeConstant<Tvec>>(ndust);
2345 {
2346 node_set_tj->set_edges(input_t_j, part_counts_with_ghost, t_j_field);
2347 }
2348 node_set_tj->evaluate();
2349
2350 } else if (EpsteinDrag *cfg_drag = std::get_if<EpsteinDrag>(&cfg.dust_drag_mode)) {
2351
2352 std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> input_gamma
2353 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("", "");
2354 input_gamma->data = cfg_drag->gamma;
2355
2356 std::shared_ptr<shamrock::solvergraph::IDataEdge<std::vector<Tscal>>> input_sgrain_j
2357 = shamrock::solvergraph::IDataEdge<std::vector<Tscal>>::make_shared("", "");
2358 input_sgrain_j->data = cfg_drag->grains_sizes;
2359
2360 std::shared_ptr<shamrock::solvergraph::IDataEdge<std::vector<Tscal>>> input_rho_grain_j
2361 = shamrock::solvergraph::IDataEdge<std::vector<Tscal>>::make_shared("", "");
2362 input_rho_grain_j->data = cfg_drag->grains_densities;
2363
2364 std::shared_ptr<modules::SetDustStoppingTimeEpstein<Tvec, Kern>> node_set_tj
2365 = std::make_shared<modules::SetDustStoppingTimeEpstein<Tvec, Kern>>(ndust);
2366 {
2367 node_set_tj->set_edges(
2368 gpart_mass,
2369 input_gamma,
2370 input_sgrain_j,
2371 input_rho_grain_j,
2372 part_counts_with_ghost,
2373 hpart_refs,
2374 storage.soundspeed,
2375 t_j_field);
2376 }
2377 node_set_tj->evaluate();
2378 }
2379
2380 if (cfg.ballabio_ts_limiter) {
2381 std::shared_ptr<modules::BallabioTsLimiter<Tvec>> node_ballabio_ts_limiter
2382 = std::make_shared<modules::BallabioTsLimiter<Tvec>>(ndust);
2383 {
2384 node_ballabio_ts_limiter->set_edges(
2385 part_counts_with_ghost, hpart_refs, storage.soundspeed, t_j_field);
2386 }
2387 node_ballabio_ts_limiter->evaluate();
2388 }
2389
2390 // delta v computation (for CFL or other uses e.g. COALA)
2391 auto &pressure_field = storage.pressure;
2392 auto &xyz_refs = storage.positions_with_ghosts;
2393
2394 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> omega_refs
2395 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("omega", "omega");
2396 {
2397 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
2398 shambase::get_check_ref(omega_refs)
2399 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
2400 [iomega_interf](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
2401 return std::ref(mpdat.get_field<Tscal>(iomega_interf));
2402 }));
2403 }
2404
2405 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> s_j_refs
2406 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("s_j", "s_j");
2407 {
2408 u32 is_j_interf = ghost_layout.get_field_idx<Tscal>("s_j");
2409 shambase::get_check_ref(s_j_refs).set_refs(
2410 mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
2411 [is_j_interf](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
2412 return std::ref(mpdat.get_field<Tscal>(is_j_interf));
2413 }));
2414 }
2415
2416 std::shared_ptr<shamrock::solvergraph::Field<Tvec>> grad_P_on_rho
2417 = std::make_shared<shamrock::solvergraph::Field<Tvec>>(1, "grad P/rho", "grad P/rho");
2418
2419 u32 idelta_v = pdl.get_field_idx<Tvec>("delta_v");
2420
2421 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> delta_v
2422 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("Delta v", "Delta v");
2423 map_field_refs(scheduler(), idelta_v, *delta_v);
2424
2425 auto press_grad_node = std::make_shared<modules::NodeComputePressureGrad<Tvec, Kern>>();
2426 auto delta_v_node = std::make_shared<modules::MonoFluidTVADeltav<Tvec, Kern>>(ndust);
2427
2428 press_grad_node->set_edges(
2429 gpart_mass,
2430 part_counts,
2431 part_counts_with_ghost,
2432 xyz_refs,
2433 hpart_refs,
2434 omega_refs,
2435 pressure_field,
2436 storage.neigh_cache,
2437 grad_P_on_rho);
2438
2439 delta_v_node->set_edges(
2440 gpart_mass, part_counts, hpart_refs, grad_P_on_rho, s_j_refs, t_j_field, delta_v);
2441
2442 press_grad_node->evaluate();
2443 delta_v_node->evaluate();
2444 }
2445
2446 modules::UpdateDerivs<Tvec, Kern> derivs(context, solver_config, storage);
2447 derivs.update_derivs(dt_hydro);
2448
2449 modules::ExternalForces<Tvec, Kern> ext_forces(context, solver_config, storage);
2450 ext_forces.add_ext_forces();
2451}
2452
2453template<class Tvec, template<class> class Kern>
2454bool shammodels::sph::Solver<Tvec, Kern>::apply_corrector(Tscal dt, u64 Npart_all) {
2455 return false;
2456}
2457
2458template<class Tvec, template<class> class Kern>
2459void shammodels::sph::Solver<Tvec, Kern>::update_sync_load_values() {
2460 modules::ComputeLoadBalanceValue<Tvec, Kern>(context, solver_config, storage)
2461 .update_load_balancing();
2462 scheduler().scheduler_step(false, false);
2463}
2464
2465template<class Tvec, template<class> class Kern>
2466shammodels::sph::TimestepLog shammodels::sph::Solver<Tvec, Kern>::evolve_once() {
2467
2468 // has to be first since there is a barrier that may mess the other timers
2469 shamsys::SystemMetrics system_metrics_start = shamsys::get_system_metrics();
2470
2471 sham::MemPerfInfos mem_perf_infos_start = sham::details::get_mem_perf_info();
2472 f64 mpi_timer_start = shamcomm::mpi::get_timer("total");
2473
2474 for (auto &callbacks : timestep_callbacks) {
2475 if (callbacks.step_begin_callback) {
2476 shambase::get_check_ref(callbacks.step_begin_callback)();
2477 }
2478 }
2479
2480 Tscal t_current = get_time();
2481 Tscal dt = get_dt_sph();
2482
2483 StackEntry stack_loc{};
2484
2485 if (shamcomm::world_rank() == 0) {
2486 shamcomm::logs::raw_ln(
2487 sham::format("---------------- t = {}, dt = {} ----------------", t_current, dt));
2488 }
2489
2490 shambase::Timer tstep;
2491 tstep.start();
2492
2493 // if(shamcomm::world_rank() == 0) std::cout << scheduler().dump_status() << std::endl;
2494 modules::ComputeLoadBalanceValue<Tvec, Kern>(context, solver_config, storage)
2495 .update_load_balancing();
2496 scheduler().scheduler_step(true, true);
2497 modules::ComputeLoadBalanceValue<Tvec, Kern>(context, solver_config, storage)
2498 .update_load_balancing();
2499 // if(shamcomm::world_rank() == 0) std::cout << scheduler().dump_status() << std::endl;
2500 scheduler().scheduler_step(false, false);
2501 // if(shamcomm::world_rank() == 0) std::cout << scheduler().dump_status() << std::endl;
2502
2504
2505 using namespace shamrock;
2506 using namespace shamrock::patch;
2507
2508 bool has_B_field = solver_config.has_field_B_on_rho();
2509 bool has_psi_field = solver_config.has_field_psi_on_ch();
2510 bool has_epsilon_field = solver_config.dust_config.has_epsilon_field();
2511 bool has_deltav_field = solver_config.dust_config.has_deltav_field();
2512 bool has_s_j_field = solver_config.dust_config.has_s_j_field();
2513
2514 PatchDataLayerLayout &pdl = scheduler().pdl_old();
2515
2516 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
2517 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
2518 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
2519 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
2520 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
2521 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
2522 const u32 iB_on_rho = (has_B_field) ? pdl.get_field_idx<Tvec>("B/rho") : 0;
2523 const u32 idB_on_rho = (has_B_field) ? pdl.get_field_idx<Tvec>("dB/rho") : 0;
2524 const u32 ipsi_on_ch = (has_psi_field) ? pdl.get_field_idx<Tscal>("psi/ch") : 0;
2525 const u32 idpsi_on_ch = (has_psi_field) ? pdl.get_field_idx<Tscal>("dpsi/ch") : 0;
2526 const u32 iepsilon = (has_epsilon_field) ? pdl.get_field_idx<Tscal>("epsilon") : 0;
2527 const u32 idtepsilon = (has_epsilon_field) ? pdl.get_field_idx<Tscal>("dtepsilon") : 0;
2528 const u32 is_j = (has_s_j_field) ? pdl.get_field_idx<Tscal>("s_j") : 0;
2529 const u32 ids_j_dt = (has_s_j_field) ? pdl.get_field_idx<Tscal>("ds_j_dt") : 0;
2530 const u32 ideltav = (has_deltav_field) ? pdl.get_field_idx<Tvec>("deltav") : 0;
2531 const u32 idtdeltav = (has_deltav_field) ? pdl.get_field_idx<Tvec>("dtdeltav") : 0;
2532
2533 shamrock::SchedulerUtility utility(scheduler());
2534
2535 {
2536 // beginning of SolverGraph migration
2537
2538 using namespace shamrock::solvergraph;
2539
2540 SolverGraph &solver_graph = storage.solver_graph;
2541
2543 // Solver evaluation
2545
2546 shambase::get_check_ref(storage.solver_sequence).evaluate();
2547 }
2548
2549 modules::ExternalForces<Tvec, Kern> ext_forces(context, solver_config, storage);
2550 ext_forces.compute_ext_forces_indep_v();
2551
2552 gen_serial_patch_tree();
2553
2554 apply_position_boundary(t_current + dt);
2555
2556 u64 Npart_all = scheduler().get_total_obj_count();
2557
2558 if (solver_config.enable_particle_reordering
2559 && solve_logs.step_count % solver_config.particle_reordering_step_freq == 0) {
2560 logger::info_ln("SPH", "Reordering particles at step ", solve_logs.step_count);
2561 modules::ParticleReordering<Tvec, u_morton, Kern>(context, solver_config, storage)
2562 .reorder_particles();
2563 }
2564
2565 {
2566 // update part counts and spans since particles have been moved and thus
2567 // new patch can be non-empty/empty
2568 using namespace shamrock::solvergraph;
2569 SolverGraph &solver_graph = storage.solver_graph;
2570 solver_graph.get_node_ref_base("attach fields to scheduler").evaluate();
2571 }
2572
2573 sph_prestep(t_current, dt);
2574
2575 using RTree = shamtree::CompressedLeafBVH<u_morton, Tvec, 3>;
2576
2577 // Here we will add self grav to the external forces indep of vel (this will be moved into a
2578 // sperate module later)
2579 if (solver_config.self_grav_config.is_sg_on()) {
2580 using namespace shamrock::solvergraph;
2581 SolverGraph &solver_graph = storage.solver_graph;
2582 solver_graph.get_node_ref_base("self gravity sequence").evaluate();
2583 }
2584
2585 sph::BasicSPHGhostHandler<Tvec> &ghost_handle = storage.ghost_handler.get();
2586 auto &merged_xyzh = storage.merged_xyzh.get();
2587 shambase::DistributedData<RTree> &trees = storage.merged_pos_trees.get();
2588 // ComputeField<Tscal> &omega = storage.omega.get();
2589
2590 shamrock::patch::PatchDataLayerLayout &ghost_layout
2591 = shambase::get_check_ref(storage.ghost_layout.get());
2592 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
2593 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
2594 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
2595 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
2596 u32 iB_on_rho_interf = (has_B_field) ? ghost_layout.get_field_idx<Tvec>("B/rho") : 0;
2597 u32 ipsi_on_rho_interf = (has_psi_field) ? ghost_layout.get_field_idx<Tscal>("psi/ch") : 0;
2598
2599 using RTreeField = RadixTreeField<Tscal>;
2600 shambase::DistributedData<RTreeField> rtree_field_h;
2601
2602 Tscal next_cfl = 0;
2603
2604 u32 corrector_iter_cnt = 0;
2605 bool need_rerun_corrector = false;
2606 do {
2607
2608 reset_merge_ghosts_fields();
2609 reset_eos_fields();
2610
2611 if (corrector_iter_cnt == 50) {
2612 throw shambase::make_except_with_loc<std::runtime_error>(
2613 "the corrector has made over 50 loops, either their is a bug, either you are using "
2614 "a dt that is too large");
2615 }
2616
2617 // communicate fields
2618 communicate_merge_ghosts_fields();
2619
2620 if (solver_config.has_field_alphaAV()) {
2621
2622 std::shared_ptr<shamrock::solvergraph::PatchDataLayerRefs> patchdatas
2623 = std::make_shared<shamrock::solvergraph::PatchDataLayerRefs>(
2624 "patchdata_layer_ref", "patchdata_layer_ref");
2625
2626 auto node_set_edge = scheduler().get_node_set_edge_patchdata_layer_refs();
2627 node_set_edge->set_edges(patchdatas);
2628 node_set_edge->evaluate();
2629
2630 shamrock::solvergraph::CopyPatchDataFieldFromLayer<Tscal> node_copy(
2631 scheduler().get_layout_ptr_old(), "alpha_AV");
2632 node_copy.set_edges(patchdatas, storage.alpha_av_updated);
2633 node_copy.evaluate();
2634 }
2635
2636 if (solver_config.has_field_dtdivv()) {
2637
2638 if (solver_config.combined_dtdiv_divcurlv_compute) {
2639 if (solver_config.has_field_dtdivv()) {
2640 sph::modules::DiffOperatorDtDivv<Tvec, Kern>(context, solver_config, storage)
2641 .update_dtdivv(true);
2642 }
2643 } else {
2644
2645 if (solver_config.has_field_divv()) {
2646 sph::modules::DiffOperators<Tvec, Kern>(context, solver_config, storage)
2647 .update_divv();
2648 }
2649
2650 if (solver_config.has_field_curlv()) {
2651 sph::modules::DiffOperators<Tvec, Kern>(context, solver_config, storage)
2652 .update_curlv();
2653 }
2654
2655 if (solver_config.has_field_dtdivv()) {
2656 sph::modules::DiffOperatorDtDivv<Tvec, Kern>(context, solver_config, storage)
2657 .update_dtdivv(false);
2658 }
2659 }
2660
2661 } else {
2662 if (solver_config.has_field_divv()) {
2663 sph::modules::DiffOperators<Tvec, Kern>(context, solver_config, storage)
2664 .update_divv();
2665 }
2666
2667 if (solver_config.has_field_curlv()) {
2668 sph::modules::DiffOperators<Tvec, Kern>(context, solver_config, storage)
2669 .update_curlv();
2670 }
2671 }
2672
2673 // if (solver_config.has_field_divB()) {
2674 // sph::modules::DiffOperatorsB<Tvec, Kern>(context, solver_config, storage)
2675 // .update_divB();
2676 // }
2677
2678 // if (solver_config.has_field_curlB()) {
2679 // sph::modules::DiffOperatorsB<Tvec, Kern>(context, solver_config, storage)
2680 // .update_curlB();
2681 // }
2682 update_artificial_viscosity(dt);
2683
2684 if (solver_config.has_field_alphaAV()) {
2685
2686 shamrock::solvergraph::Field<Tscal> &comp_field_send
2687 = shambase::get_check_ref(storage.alpha_av_updated);
2688
2689 using InterfaceBuildInfos =
2690 typename sph::BasicSPHGhostHandler<Tvec>::InterfaceBuildInfos;
2691
2692 shambase::Timer time_interf;
2693 time_interf.start();
2694
2695 auto field_interf = ghost_handle.template build_interface_native<PatchDataField<Tscal>>(
2696 storage.ghost_patch_cache.get(),
2697 [&](u64 sender,
2698 u64 /*receiver*/,
2699 InterfaceBuildInfos binfo,
2700 sham::DeviceBuffer<u32> &buf_idx,
2701 u32 cnt) -> PatchDataField<Tscal> {
2702 PatchDataField<Tscal> &sender_field = comp_field_send.get_field(sender);
2703
2704 return sender_field.make_new_from_subset(buf_idx, cnt);
2705 });
2706
2707 shambase::DistributedDataShared<PatchDataField<Tscal>> interf_pdat
2708 = ghost_handle.communicate_pdatfield(
2709 std::move(field_interf), 1, storage.exchange_gz_alpha);
2710
2711 shambase::DistributedData<PatchDataField<Tscal>> merged_field
2712 = ghost_handle.template merge_native<PatchDataField<Tscal>, PatchDataField<Tscal>>(
2713 std::move(interf_pdat),
2714 [&](const shamrock::patch::Patch p, shamrock::patch::PatchDataLayer &pdat) {
2715 PatchDataField<Tscal> &receiver_field
2716 = comp_field_send.get_field(p.id_patch);
2717 return receiver_field.duplicate();
2718 },
2719 [](PatchDataField<Tscal> &mpdat, PatchDataField<Tscal> &pdat_interf) {
2720 mpdat.insert(pdat_interf);
2721 });
2722
2723 time_interf.stop();
2724 storage.timings_details.interface += time_interf.elapsed_sec();
2725
2726 storage.alpha_av_ghost.set(std::move(merged_field));
2727 }
2728
2729 // compute pressure
2730 compute_eos_fields();
2731
2732 constexpr bool debug_interfaces = false;
2733 if constexpr (debug_interfaces) {
2734
2735 if (solver_config.do_debug_dump) {
2736
2737 shambase::DistributedData<MergedPatchData> &mpdat
2738 = storage.merged_patchdata_ghost.get();
2739
2740 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
2741 MergedPatchData &merged_patch = mpdat.get(cur_p.id_patch);
2742 PatchDataLayer &mpdat = merged_patch.pdat;
2743
2744 sycl::buffer<Tvec> &buf_xyz = shambase::get_check_ref(
2745 merged_xyzh.get(cur_p.id_patch).field_pos.get_buf());
2746 sycl::buffer<Tvec> &buf_vxyz = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);
2747 sycl::buffer<Tscal> &buf_hpart = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);
2748
2749 u32 total_elements = shambase::get_check_ref(storage.part_counts_with_ghost)
2750 .indexes.get(cur_p.id_patch);
2751 SHAM_ASSERT(merged_patch.total_elements == total_elements);
2752
2753 Debug_ph_dump<Tvec> info{
2754 total_elements,
2755 solver_config.gpart_mass,
2756
2757 buf_xyz,
2758 buf_hpart,
2759 buf_vxyz};
2760
2761 make_interface_debug_phantom_dump(info).gen_file().write_to_file(
2762 solver_config.debug_dump_filename);
2763 logger::raw_ln("writing : ", solver_config.debug_dump_filename);
2764 });
2765 }
2766 }
2767
2768 // compute force
2769 shamlog_debug_ln("sph::BasicGas", "compute force");
2770
2771 // save old acceleration
2772 prepare_corrector();
2773
2774 update_derivs(dt);
2775
2777 // Gravitational Wave emission
2779 bool compute_GW = solver_config.compute_gw;
2780
2781 if (compute_GW) {
2782 using namespace shamrock::solvergraph;
2783 using GW = shammodels::common::modules::ComputeGravWave<Tvec>;
2784
2785 auto central_pos = IDataEdge<Tvec>::make_shared("x_0", "\\mathbf{x}_0");
2786 central_pos->data = Tvec{0, 0, 0};
2787
2788 auto central_vel = IDataEdge<Tvec>::make_shared("v_0", "\\mathbf{v}_0");
2789 central_vel->data = Tvec{0, 0, 0};
2790
2791 auto central_acc = IDataEdge<Tvec>::make_shared("a_0", "\\mathbf{a}_0");
2792 central_acc->data = Tvec{0, 0, 0};
2793
2794 auto gw_prefactor = IDataEdge<Tscal>::make_shared("gw_prefactor", "gw_prefactor");
2795 gw_prefactor->data = Tscal(1); // should be G/c^2D
2796
2797 auto theta_gw = IDataEdge<Tscal>::make_shared("theta_gw", "\\theta_{\\rm gw}");
2798 theta_gw->data = Tscal(0);
2799
2800 auto phi_gw = IDataEdge<Tscal>::make_shared("phi_gw", "\\phi_{\\rm gw}");
2801 phi_gw->data = Tscal(0);
2802
2803 ComputeField<Tscal> gw_mass_field
2804 = utility.make_compute_field<Tscal>("gw_mass", 1, solver_config.gpart_mass);
2805
2806 auto spans_masses = std::make_shared<FieldRefs<Tscal>>("m", "m");
2807 map_field_refs_ext(scheduler(), gw_mass_field, *spans_masses);
2808
2809 const u32 iaxyz_ext = pdl.get_field_idx<Tvec>("axyz_ext");
2810 auto spans_accel_ext
2811 = std::make_shared<FieldRefs<Tvec>>("axyz_ext", "\\mathbf{a}_{\\rm ext}");
2812 map_field_refs(scheduler(), iaxyz_ext, *spans_accel_ext);
2813
2814 auto ddq = IDataEdge<typename GW::Tddq>::make_shared("ddq", "\\ddot{Q}");
2815 auto ddq_xy = IDataEdge<typename GW::Tddqxy>::make_shared("ddq_xy", "\\ddot{Q}_{xy}");
2816 auto hx = IDataEdge<typename GW::Th>::make_shared("hx", "h_x");
2817 auto hp = IDataEdge<typename GW::Th>::make_shared("hp", "h_+");
2818
2819 GW node_computeGW{};
2820 node_computeGW.set_edges(
2821 storage.solver_graph.template get_edge_ptr<FieldRefs<Tvec>>("xyz"),
2822 storage.solver_graph.template get_edge_ptr<FieldRefs<Tvec>>("vxyz"),
2823 storage.solver_graph.template get_edge_ptr<FieldRefs<Tvec>>("axyz"),
2824 spans_masses,
2825 spans_accel_ext,
2826 central_pos,
2827 central_vel,
2828 central_acc,
2829 gw_prefactor,
2830 theta_gw,
2831 phi_gw,
2832 storage.part_counts,
2833 ddq,
2834 ddq_xy,
2835 hx,
2836 hp);
2837
2838 node_computeGW.evaluate();
2839
2840 // TODO: send that somewhere rather than doing a print
2841 logger::raw_ln("################## hx = ", hx->data);
2842 logger::raw_ln("################## hp = ", hp->data);
2843 logger::raw_ln("################## ddq = ", ddq->data);
2844 logger::raw_ln("################## ddq_xy = ", ddq_xy->data);
2845 }
2846
2847 bool has_luminosity = solver_config.compute_luminosity;
2848
2849 if (has_luminosity) {
2850 const u32 iluminosity = pdl.get_field_idx<Tscal>("luminosity");
2851
2852 shambase::get_check_ref(storage.hpart_with_ghosts)
2853 .set_refs(storage.merged_xyzh.get()
2854 .template map<std::reference_wrapper<PatchDataField<Tscal>>>(
2855 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
2856 return std::ref(mpdat.get_field<Tscal>(
2857 1)); // hpart is at index 1 in merged_xyzh
2858 }));
2859
2860 shambase::get_check_ref(storage.hpart_with_ghosts)
2861 .set_refs(storage.merged_xyzh.get()
2862 .template map<std::reference_wrapper<PatchDataField<Tscal>>>(
2863 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
2864 return std::ref(mpdat.get_field<Tscal>(1));
2865 }));
2866
2867 auto uint_with_ghost = shamrock::solvergraph::FieldRefs<Tscal>::make_shared("", "");
2868 shamrock::solvergraph::NodeSetEdge<shamrock::solvergraph::FieldRefs<Tscal>>
2869 set_uint_with_ghost_refs(
2870 [&](shamrock::solvergraph::FieldRefs<Tscal> &field_uint_with_ghost_edge) {
2871 shambase::DistributedData<PatchDataLayer> &mpdats
2872 = storage.merged_patchdata_ghost.get();
2873
2874 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> field_uint_with_ghost_refs
2875 = {};
2876
2877 scheduler().for_each_patchdata_nonempty(
2878 [&](const Patch p, PatchDataLayer &pdat) {
2879 PatchDataLayer &mpdat = mpdats.get(p.id_patch);
2880
2881 auto &field = mpdat.get_field<Tscal>(iuint_interf);
2882 field_uint_with_ghost_refs.add_obj(p.id_patch, std::ref(field));
2883 });
2884
2885 field_uint_with_ghost_edge.set_refs(field_uint_with_ghost_refs);
2886 });
2887
2888 set_uint_with_ghost_refs.set_edges(uint_with_ghost);
2889
2890 auto omega_with_ghost = shamrock::solvergraph::FieldRefs<Tscal>::make_shared("", "");
2891 shamrock::solvergraph::NodeSetEdge<shamrock::solvergraph::FieldRefs<Tscal>>
2892 set_omega_with_ghost_refs([&](shamrock::solvergraph::FieldRefs<Tscal>
2893 &field_omega_with_ghost_edge) {
2894 shambase::DistributedData<PatchDataLayer> &mpdats
2895 = storage.merged_patchdata_ghost.get();
2896
2897 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> field_omega_with_ghost_refs
2898 = {};
2899
2900 scheduler().for_each_patchdata_nonempty(
2901 [&](const Patch p, PatchDataLayer &pdat) {
2902 PatchDataLayer &mpdat = mpdats.get(p.id_patch);
2903
2904 auto &field = mpdat.get_field<Tscal>(iomega_interf);
2905 field_omega_with_ghost_refs.add_obj(p.id_patch, std::ref(field));
2906 });
2907
2908 field_omega_with_ghost_edge.set_refs(field_omega_with_ghost_refs);
2909 });
2910
2911 set_omega_with_ghost_refs.set_edges(omega_with_ghost);
2912
2913 auto luminosity = shamrock::solvergraph::FieldRefs<Tscal>::make_shared("", "");
2914
2915 shamrock::solvergraph::NodeSetEdge<shamrock::solvergraph::FieldRefs<Tscal>>
2916 set_luminosity_refs(
2917 [&](shamrock::solvergraph::FieldRefs<Tscal> &field_luminosity_edge) {
2918 shambase::DistributedData<PatchDataLayer> &mpdats
2919 = storage.merged_patchdata_ghost.get();
2920
2921 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> field_luminosity_refs
2922 = {};
2923
2924 scheduler().for_each_patchdata_nonempty(
2925 [&](const Patch p, PatchDataLayer &pdat) {
2926 auto &field = pdat.get_field<Tscal>(iluminosity);
2927 field_luminosity_refs.add_obj(p.id_patch, std::ref(field));
2928 });
2929 field_luminosity_edge.set_refs(field_luminosity_refs);
2930 });
2931
2932 set_luminosity_refs.set_edges(luminosity);
2933
2934 set_uint_with_ghost_refs.evaluate();
2935 set_omega_with_ghost_refs.evaluate();
2936 set_luminosity_refs.evaluate();
2937
2938 Tscal alpha_u = solver_config.artif_viscosity.get_alpha_u().value();
2939
2940 modules::NodeComputeLuminosity<Tvec, Kern> compute_luminosity{
2941 solver_config.gpart_mass, alpha_u};
2942
2943 compute_luminosity.set_edges(
2944 storage.part_counts,
2945 storage.part_counts_with_ghost,
2946 storage.neigh_cache,
2947 storage.positions_with_ghosts,
2948 storage.hpart_with_ghosts,
2949 omega_with_ghost,
2950 uint_with_ghost,
2951 storage.pressure,
2952 luminosity);
2953
2954 compute_luminosity.evaluate();
2955 }
2956
2957 modules::ConservativeCheck<Tvec, Kern> cv_check(context, solver_config, storage);
2958 cv_check.check_conservation();
2959
2960 ComputeField<Tscal> vepsilon_v_sq
2961 = utility.make_compute_field<Tscal>("vmean epsilon_v^2", 1);
2962 ComputeField<Tscal> uepsilon_u_sq
2963 = utility.make_compute_field<Tscal>("umean epsilon_u^2", 1);
2964
2965 // corrector
2966 shamlog_debug_ln("sph::BasicGas", "leapfrog corrector");
2967 utility.fields_leapfrog_corrector<Tvec>(
2968 ivxyz, iaxyz, storage.old_axyz.get(), vepsilon_v_sq, dt / 2);
2969 utility.fields_leapfrog_corrector<Tscal>(
2970 iuint, iduint, storage.old_duint.get(), uepsilon_u_sq, dt / 2);
2971
2972 if (solver_config.has_field_B_on_rho()) {
2973 ComputeField<Tscal> BOR_epsilon_BOR_sq
2974 = utility.make_compute_field<Tscal>("B/rho epsilon_B/rho^2", 1);
2975 utility.fields_leapfrog_corrector<Tvec>(
2976 iB_on_rho, idB_on_rho, storage.old_dB_on_rho.get(), BOR_epsilon_BOR_sq, dt / 2);
2977 }
2978 if (solver_config.has_field_B_on_rho()) {
2979 ComputeField<Tscal> POC_epsilon_POC_sq
2980 = utility.make_compute_field<Tscal>("psi/ch epsilon_psi/ch^2", 1);
2981 utility.fields_leapfrog_corrector<Tscal>(
2982 ipsi_on_ch, idpsi_on_ch, storage.old_dpsi_on_ch.get(), POC_epsilon_POC_sq, dt / 2);
2983 }
2984
2985 if (solver_config.dust_config.has_epsilon_field()) {
2986 ComputeField<Tscal> epsilon_epsilon_sq
2987 = utility.make_compute_field<Tscal>("epsilon epsilon^2", 1);
2988 utility.fields_leapfrog_corrector<Tscal>(
2989 iepsilon, idtepsilon, storage.old_dtepsilon.get(), epsilon_epsilon_sq, dt / 2);
2990 }
2991
2992 if (solver_config.dust_config.has_deltav_field()) {
2993 ComputeField<Tscal> epsilon_deltav_sq
2994 = utility.make_compute_field<Tscal>("deltav deltav^2", 1);
2995 utility.fields_leapfrog_corrector<Tvec>(
2996 ideltav, idtdeltav, storage.old_dtdeltav.get(), epsilon_deltav_sq, dt / 2);
2997 }
2998
2999 if (solver_config.dust_config.has_s_j_field()) {
3000 ComputeField<Tscal> s_j_s_j_sq = utility.make_compute_field<Tscal>(
3001 "s_j s_j^2", solver_config.dust_config.get_dust_nvar());
3002 bool ensure_positivity
3003 = solver_config.dust_config.get_monofluid_tva().ensure_s_j_positivity;
3004 if (ensure_positivity) {
3005 utility.fields_leapfrog_corrector_positive_only<Tscal>(
3006 is_j, ids_j_dt, storage.old_ds_j_dt.get(), s_j_s_j_sq, dt / 2);
3007 } else {
3008 utility.fields_leapfrog_corrector<Tscal>(
3009 is_j, ids_j_dt, storage.old_ds_j_dt.get(), s_j_s_j_sq, dt / 2);
3010 }
3011
3012 auto &monofluid_tva_cfg = solver_config.dust_config.get_monofluid_tva();
3013 if (monofluid_tva_cfg.should_clamp_dust_density()) {
3014 auto hfactd_edge
3015 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("hfactd", "hfactd");
3016 hfactd_edge->data = Kernel::hfactd;
3017
3018 auto clamp_frac_edge = shamrock::solvergraph::IDataEdge<Tscal>::make_shared(
3019 "clamp_frac", "clamp_frac");
3020 clamp_frac_edge->data = monofluid_tva_cfg.get_clamp_dust_frac();
3021
3022 shammodels::sph::modules::NodeMonofluidTVADustDensityClamp<Tvec> density_clamp(
3023 solver_config.dust_config.get_dust_nvar());
3024 density_clamp.set_edges(
3025 storage.solver_graph.template get_edge_ptr<shamrock::solvergraph::Indexes<u32>>(
3026 "part_counts"),
3027 storage.solver_graph
3028 .template get_edge_ptr<shamrock::solvergraph::IDataEdge<Tscal>>(
3029 "gpart_mass"),
3030 hfactd_edge,
3031 clamp_frac_edge,
3032 storage.solver_graph
3033 .template get_edge_ptr<shamrock::solvergraph::FieldRefs<Tscal>>("hpart"),
3034 storage.solver_graph
3035 .template get_edge_ptr<shamrock::solvergraph::FieldRefs<Tscal>>("s_j"));
3036 density_clamp.evaluate();
3037 }
3038 }
3039
3040 storage.old_axyz.reset();
3041 storage.old_duint.reset();
3042 if (solver_config.has_field_B_on_rho()) {
3043 storage.old_dB_on_rho.reset();
3044 }
3045 if (solver_config.has_field_B_on_rho()) {
3046 storage.old_dpsi_on_ch.reset();
3047 }
3048
3049 if (solver_config.dust_config.has_epsilon_field()) {
3050 storage.old_dtepsilon.reset();
3051 }
3052
3053 if (solver_config.dust_config.has_deltav_field()) {
3054 storage.old_dtdeltav.reset();
3055 }
3056
3057 if (solver_config.dust_config.has_s_j_field()) {
3058 storage.old_ds_j_dt.reset();
3059 }
3060
3061 Tscal rank_veps_v = sycl::sqrt(vepsilon_v_sq.compute_rank_max());
3063 // compute means //////////////////////////
3065
3066 Tscal sum_vsq = utility.compute_rank_dot_sum<Tvec>(ivxyz);
3067
3068 Tscal vmean_sq = shamalgs::collective::allreduce_sum(sum_vsq) / Tscal(Npart_all);
3069
3070 Tscal vmean = sycl::sqrt(vmean_sq);
3071
3072 Tscal rank_eps_v = rank_veps_v / vmean;
3073
3074 if (vmean <= 0) {
3075 rank_eps_v = 0;
3076 }
3077
3078 Tscal eps_v = shamalgs::collective::allreduce_max(rank_eps_v);
3079
3080 shamlog_debug_ln("BasicGas", "epsilon v :", eps_v);
3081
3082 if (eps_v > 1e-2) {
3083 if (shamcomm::world_rank() == 0) {
3084 logger::warn_ln(
3085 "BasicGasSPH",
3086 sham::format(
3087 "the corrector tolerance are broken the step will "
3088 "be re rerunned\n eps_v = {}",
3089 eps_v));
3090 }
3091 need_rerun_corrector = true;
3092 set_cfl_multipler(get_cfl_multipler() / 2);
3093
3094 // logger::info_ln("rerun corrector ...");
3095 } else {
3096 need_rerun_corrector = false;
3097 }
3098
3099 if (!need_rerun_corrector) {
3100
3101 storage.solver_graph.get_node_ref_base("sink corrector").evaluate();
3102
3103 // write back alpha av field
3104 if (solver_config.has_field_alphaAV()) {
3105
3106 const u32 ialpha_AV = pdl.get_field_idx<Tscal>("alpha_AV");
3107 shamrock::solvergraph::Field<Tscal> &alpha_av_updated
3108 = shambase::get_check_ref(storage.alpha_av_updated);
3109
3110 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
3111 sham::DeviceBuffer<Tscal> &buf_alpha_av
3112 = pdat.get_field<Tscal>(ialpha_AV).get_buf();
3113 sham::DeviceBuffer<Tscal> &buf_alpha_av_updated
3114 = alpha_av_updated.get_field(cur_p.id_patch).get_buf();
3115
3116 auto &q = shamsys::instance::get_compute_scheduler().get_queue();
3117 sham::EventList depends_list;
3118
3119 auto alpha_av = buf_alpha_av.get_write_access(depends_list);
3120 auto alpha_av_updated = buf_alpha_av_updated.get_read_access(depends_list);
3121
3122 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
3123 shambase::parallel_for(
3124 cgh, pdat.get_obj_cnt(), "write back alpha_av", [=](i32 id_a) {
3125 alpha_av[id_a] = alpha_av_updated[id_a];
3126 });
3127 });
3128
3129 buf_alpha_av.complete_event_state(e);
3130 buf_alpha_av_updated.complete_event_state(e);
3131 });
3132 }
3133
3134 shamlog_debug_ln("BasicGas", "computing next CFL");
3135
3136 // Update element counts
3137 shambase::get_check_ref(storage.part_counts).indexes
3138 = storage.merged_xyzh.get().template map<u32>(
3139 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
3140 return scheduler().patch_data.get_pdat(id).get_obj_cnt();
3141 });
3142
3143 std::shared_ptr<shamrock::solvergraph::Field<Tscal>> vsig_max_dt
3144 = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
3145 1, "vsig_a", "v_{\\rm sig}");
3146 vsig_max_dt->ensure_sizes(shambase::get_check_ref(storage.part_counts).indexes);
3147
3148 std::shared_ptr<shamrock::solvergraph::Field<Tscal>> vclean_dt;
3149 if (has_psi_field) {
3150 vclean_dt = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
3151 1, "vclean_a", "v_{\\rm clean}");
3152 vclean_dt->ensure_sizes(shambase::get_check_ref(storage.part_counts).indexes);
3153 }
3154
3155 shambase::DistributedData<PatchDataLayer> &mpdats
3156 = storage.merged_patchdata_ghost.get();
3157
3158 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
3159 PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch);
3160
3161 sham::DeviceBuffer<Tvec> &buf_xyz
3162 = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref<Tvec>(0);
3163 sham::DeviceBuffer<Tvec> &buf_vxyz = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);
3164 sham::DeviceBuffer<Tscal> &buf_hpart
3165 = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);
3166 sham::DeviceBuffer<Tscal> &buf_uint = mpdat.get_field_buf_ref<Tscal>(iuint_interf);
3167 sham::DeviceBuffer<Tscal> &buf_pressure
3168 = shambase::get_check_ref(storage.pressure).get_field(cur_p.id_patch).get_buf();
3169 sham::DeviceBuffer<Tscal> &cs_buf = shambase::get_check_ref(storage.soundspeed)
3170 .get_field(cur_p.id_patch)
3171 .get_buf();
3172
3173 sham::DeviceBuffer<Tscal> &vsig_buf = vsig_max_dt->get_buf(cur_p.id_patch);
3174
3175 sycl::range range_npart{pdat.get_obj_cnt()};
3176
3177 tree::ObjectCache &pcache
3178 = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);
3179
3181
3182 {
3183
3184 auto &q = shamsys::instance::get_compute_scheduler().get_queue();
3185 sham::EventList depends_list;
3186
3187 auto xyz = buf_xyz.get_read_access(depends_list);
3188 auto vxyz = buf_vxyz.get_read_access(depends_list);
3189 auto hpart = buf_hpart.get_read_access(depends_list);
3190 auto u = buf_uint.get_read_access(depends_list);
3191 auto pressure = buf_pressure.get_read_access(depends_list);
3192 auto cs = cs_buf.get_read_access(depends_list);
3193 auto vsig = vsig_buf.get_write_access(depends_list);
3194 auto particle_looper_ptrs = pcache.get_read_access(depends_list);
3195
3196 NamedStackEntry tmppp{"compute vsig"};
3197 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
3198 const Tscal pmass = solver_config.gpart_mass;
3199 const Tscal alpha_u = 1.0;
3200 const Tscal alpha_AV = 1.0;
3201 const Tscal beta_AV = 2.0;
3202
3203 tree::ObjectCacheIterator particle_looper(particle_looper_ptrs);
3204
3205 constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;
3206
3207 shambase::parallel_for(
3208 cgh, pdat.get_obj_cnt(), "compute vsig", [=](i32 id_a) {
3209 using namespace shamrock::sph;
3210
3211 Tvec sum_axyz = {0, 0, 0};
3212 Tscal sum_du_a = 0;
3213 Tscal h_a = hpart[id_a];
3214
3215 Tvec xyz_a = xyz[id_a];
3216 Tvec vxyz_a = vxyz[id_a];
3217
3218 Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd);
3219 Tscal rho_a_sq = rho_a * rho_a;
3220 Tscal rho_a_inv = 1. / rho_a;
3221
3222 Tscal P_a = pressure[id_a];
3223
3224 const Tscal u_a = u[id_a];
3225
3226 Tscal cs_a = cs[id_a];
3227
3228 Tscal vsig_max = 0;
3229
3230 particle_looper.for_each_object(id_a, [&](u32 id_b) {
3231 // compute only omega_a
3232 Tvec dr = xyz_a - xyz[id_b];
3233 Tscal rab2 = sycl::dot(dr, dr);
3234 Tscal h_b = hpart[id_b];
3235
3236 if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) {
3237 return;
3238 }
3239
3240 Tscal rab = sycl::sqrt(rab2);
3241 Tvec vxyz_b = vxyz[id_b];
3242 Tvec v_ab = vxyz_a - vxyz_b;
3243 const Tscal u_b = u[id_b];
3244
3245 Tvec r_ab_unit = dr / rab;
3246
3247 if (rab < 1e-9) {
3248 r_ab_unit = {0, 0, 0};
3249 }
3250
3251 Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd);
3252 Tscal P_b = pressure[id_b];
3253 Tscal cs_b = cs[id_b];
3254 Tscal v_ab_r_ab = sycl::dot(v_ab, r_ab_unit);
3255 Tscal abs_v_ab_r_ab = sycl::fabs(v_ab_r_ab);
3256
3258 // internal energy update
3259 // scalar : f32 | vector : f32_3
3260 const Tscal alpha_a = alpha_AV;
3261 const Tscal alpha_b = alpha_AV;
3262
3263 Tscal vsig_a = alpha_a * cs_a + beta_AV * abs_v_ab_r_ab;
3264
3265 vsig_max = sycl::fmax(vsig_max, vsig_a);
3266 });
3267
3268 vsig[id_a] = vsig_max;
3269 });
3270 });
3271
3272 if (has_psi_field) {
3273 NamedStackEntry tmppp{"compute vclean"};
3274 Tscal const mu_0 = solver_config.get_constant_mu_0();
3275 sham::DeviceBuffer<Tscal> &vclean_buf = vclean_dt->get_buf(cur_p.id_patch);
3276
3277 Tvec *B_on_rho = mpdat.get_field_buf_ref<Tvec>(iB_on_rho_interf)
3278 .get_write_access(depends_list);
3279
3280 auto vclean = vclean_buf.get_write_access(depends_list);
3281
3282 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
3283 const Tscal pmass = solver_config.gpart_mass;
3284
3285 tree::ObjectCacheIterator particle_looper(particle_looper_ptrs);
3286
3287 constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;
3288
3289 shambase::parallel_for(
3290 cgh, pdat.get_obj_cnt(), "compute vclean", [=](i32 id_a) {
3291 using namespace shamrock::sph;
3292
3293 Tscal h_a = hpart[id_a];
3294 Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd);
3295 const Tscal u_a = u[id_a];
3296 Tscal cs_a = cs[id_a];
3297 Tvec B_a = B_on_rho[id_a] * rho_a;
3298
3299 Tscal vclean_a = shamphys::MHD_physics<Tvec, Tscal>::v_shock(
3300 cs_a, B_a, rho_a, mu_0);
3301
3302 vclean[id_a] = vclean_a;
3303 });
3304 });
3305 mpdat.get_field_buf_ref<Tvec>(iB_on_rho_interf).complete_event_state(e);
3306 vclean_buf.complete_event_state(e);
3307 };
3308
3309 buf_xyz.complete_event_state(e);
3310 buf_vxyz.complete_event_state(e);
3311 buf_hpart.complete_event_state(e);
3312 buf_uint.complete_event_state(e);
3313 buf_pressure.complete_event_state(e);
3314 cs_buf.complete_event_state(e);
3315 vsig_buf.complete_event_state(e);
3316
3317 sham::EventList resulting_events;
3318 resulting_events.add_event(e);
3319 pcache.complete_event_state(resulting_events);
3320 }
3321 });
3322
3323 std::shared_ptr<shamrock::solvergraph::Field<Tscal>> cfl_dt
3324 = std::make_shared<shamrock::solvergraph::Field<Tscal>>(
3325 1, "cfl_dt", "\\Delta t_{cfl}");
3326 cfl_dt->ensure_sizes(shambase::get_check_ref(storage.part_counts).indexes);
3327
3328 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> axyz_refs
3329 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("axyz", "\\mathbf{a}");
3330 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hpart_refs
3331 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("hpart", "h");
3332
3333 map_field_refs(scheduler(), iaxyz, *axyz_refs);
3334 map_field_refs_ext(scheduler(), mpdats, ihpart_interf, *hpart_refs);
3335
3336 auto &q = shamsys::instance::get_compute_scheduler().get_queue();
3337
3338 auto reset_dt_part_field = [&]() {
3339 if (solver_config.should_save_dt_to_fields()) {
3340 const u32 idt_part = pdl.get_field_idx<Tscal>("dt_part");
3341 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
3342 sham::DeviceBuffer<Tscal> &buf_dt_part
3343 = pdat.get_field_buf_ref<Tscal>(idt_part);
3344 buf_dt_part.fill(shambase::get_infty<Tscal>());
3345 });
3346 }
3347 };
3348
3349 auto save_dt_min_to_dt_part = [&]() {
3350 if (solver_config.should_save_dt_to_fields()) {
3351 const u32 idt_part = pdl.get_field_idx<Tscal>("dt_part");
3352 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
3353 sham::DeviceBuffer<Tscal> &buf_dt_part
3354 = pdat.get_field_buf_ref<Tscal>(idt_part);
3355 sham::DeviceBuffer<Tscal> &buf_dt = cfl_dt->get_buf(cur_p.id_patch);
3356
3357 sham::kernel_call(
3358 q,
3359 sham::MultiRef{buf_dt},
3360 sham::MultiRef{buf_dt_part},
3361 pdat.get_obj_cnt(),
3362 [](u32 id_a, const Tscal *dt, Tscal *dt_part) {
3363 dt_part[id_a] = sycl::min(dt_part[id_a], dt[id_a]);
3364 });
3365 });
3366 }
3367 };
3368
3369 // reset the cfl_dt field
3370 auto reset_cfl_dt = [&]() {
3371 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
3372 cfl_dt->get_buf(cur_p.id_patch).fill(shambase::get_infty<Tscal>());
3373 });
3374 };
3375
3376 Tscal C_cour = solver_config.cfl_config.cfl_cour * get_cfl_multipler();
3377 Tscal C_force = solver_config.cfl_config.cfl_force * get_cfl_multipler();
3378 Tscal eta_phi = solver_config.cfl_config.eta_sink;
3379
3380 std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> C_cour_edge
3381 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("C_cour", "C_{cour}");
3382 C_cour_edge->data = C_cour;
3383 std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> C_force_edge
3384 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("C_force", "C_{force}");
3385 C_force_edge->data = C_force;
3386 std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> eta_phi_edge
3387 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("eta_phi", "\\eta_{\\phi}");
3388 eta_phi_edge->data = eta_phi;
3389
3390 std::shared_ptr<ComputeCFLCourant<Tscal>> compute_cfl_courant
3391 = std::make_shared<ComputeCFLCourant<Tscal>>();
3392 compute_cfl_courant->set_edges(
3393 storage.part_counts, C_cour_edge, hpart_refs, vsig_max_dt, cfl_dt);
3394
3395 std::shared_ptr<ComputeCFLForce<Tvec>> compute_cfl_force
3396 = std::make_shared<ComputeCFLForce<Tvec>>();
3397 compute_cfl_force->set_edges(
3398 storage.part_counts, C_force_edge, hpart_refs, axyz_refs, cfl_dt);
3399
3400 std::shared_ptr<ComputeCFLDivBCleaning<Tscal>> compute_cfl_divB_cleaning;
3401 if (has_psi_field) {
3402 compute_cfl_divB_cleaning = std::make_shared<ComputeCFLDivBCleaning<Tscal>>();
3403 compute_cfl_divB_cleaning->set_edges(
3404 storage.part_counts, C_cour_edge, hpart_refs, vclean_dt, cfl_dt);
3405 }
3406
3407 std::shared_ptr<ComputeCFLDust1Fluid<Tvec>> compute_cfl_dust1_fluid;
3408 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> s_j_refs;
3409 std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> hfactd_edge;
3410
3411 if (solver_config.dust_config.has_s_j_field()) {
3412 u32 ndust = solver_config.dust_config.get_dust_nvar();
3413
3414 compute_cfl_dust1_fluid = std::make_shared<ComputeCFLDust1Fluid<Tvec>>(ndust);
3415
3416 auto t_j_field
3417 = storage.solver_graph
3418 .template get_edge_ptr<shamrock::solvergraph::Field<Tscal>>("Ts_j");
3419
3420 auto pmass_edge
3421 = storage.solver_graph
3422 .template get_edge_ptr<shamrock::solvergraph::IDataEdge<Tscal>>(
3423 "gpart_mass");
3424
3425 s_j_refs = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("s_j", "s_j");
3426
3427 hfactd_edge
3428 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("hfactd", "hfactd");
3429 hfactd_edge->data = Kernel::hfactd;
3430
3431 map_field_refs(scheduler(), is_j, *s_j_refs);
3432
3433 std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> C_1fluid_edge
3434 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared(
3435 "C_1fluid", "C_{1fluid}");
3436 C_1fluid_edge->data
3437 = solver_config.dust_config.get_monofluid_tva().C_1_fluid * get_cfl_multipler();
3438
3439 compute_cfl_dust1_fluid->set_edges(
3440 storage.part_counts,
3441 C_1fluid_edge,
3442 pmass_edge,
3443 hfactd_edge,
3444 hpart_refs,
3445 storage.soundspeed,
3446 s_j_refs,
3447 t_j_field,
3448 cfl_dt);
3449 }
3450
3451 std::shared_ptr<ComputeCFLDustDrift<Tvec>> compute_cfl_dust_drift;
3452 std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> C_drift_edge;
3453 std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> cfl_density_threshold_edge;
3454 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> delta_v_refs;
3455
3456 if (solver_config.dust_config.has_s_j_field()) {
3457 u32 ndust = solver_config.dust_config.get_dust_nvar();
3458
3459 compute_cfl_dust_drift = std::make_shared<ComputeCFLDustDrift<Tvec>>(ndust);
3460
3461 delta_v_refs = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>(
3462 "delta_v", "delta_v");
3463 const u32 idelta_v = pdl.get_field_idx<Tvec>("delta_v");
3464 map_field_refs(scheduler(), idelta_v, *delta_v_refs);
3465
3466 auto &cfg_monofluid_tva = solver_config.dust_config.get_monofluid_tva();
3467
3468 C_drift_edge
3469 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("C_drift", "C_{drift}");
3470 C_drift_edge->data = cfg_monofluid_tva.C_drift * get_cfl_multipler();
3471
3472 cfl_density_threshold_edge = shamrock::solvergraph::IDataEdge<Tscal>::make_shared(
3473 "cfl_density_threshold", "cfl_density_threshold");
3474 cfl_density_threshold_edge->data = cfg_monofluid_tva.cfl_density_threshold;
3475
3476 auto pmass_edge
3477 = storage.solver_graph
3478 .template get_edge_ptr<shamrock::solvergraph::IDataEdge<Tscal>>(
3479 "gpart_mass");
3480
3481 compute_cfl_dust_drift->set_edges(
3482 storage.part_counts,
3483 C_drift_edge,
3484 cfl_density_threshold_edge,
3485 pmass_edge,
3486 hfactd_edge,
3487 hpart_refs,
3488 s_j_refs,
3489 delta_v_refs,
3490 cfl_dt);
3491 }
3492
3493 bool show_cfl_detail = solver_config.show_cfl_detail;
3494 std::vector<std::pair<std::string, Tscal>> cfl_detail;
3495
3496 auto save_cfl_detail = [&](const char *key) {
3497 if (show_cfl_detail) {
3498 save_dt_min_to_dt_part();
3499 cfl_detail.push_back(
3500 {std::string(key), cfl_dt->get_native().compute_rank_min()});
3501 reset_cfl_dt();
3502 }
3503 };
3504
3505 reset_dt_part_field();
3506 reset_cfl_dt();
3507
3508 compute_cfl_courant->evaluate();
3509 save_cfl_detail("courant");
3510
3511 compute_cfl_force->evaluate();
3512 save_cfl_detail("force");
3513
3514 if (has_psi_field) {
3515 compute_cfl_divB_cleaning->evaluate();
3516 save_cfl_detail("divB_cleaning");
3517 }
3518
3519 if (solver_config.dust_config.has_s_j_field()) {
3520 compute_cfl_dust1_fluid->evaluate();
3521 save_cfl_detail("dust1_fluid");
3522
3523 compute_cfl_dust_drift->evaluate();
3524 save_cfl_detail("dust_drift");
3525 }
3526
3527 if (!show_cfl_detail) {
3528 save_dt_min_to_dt_part();
3529 cfl_detail.push_back({"all SPH", cfl_dt->get_native().compute_rank_min()});
3530 }
3531
3532 auto &sync = scheduler().synchronized_data;
3533 auto &pos = get_sink_pos<Tvec>(sync);
3534 if (!pos.empty()) {
3535 // sink sink CFL
3536
3537 std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> G_edge
3538 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("G", "G");
3539 G_edge->data = solver_config.get_constant_G();
3540
3541 std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> sink_sink_cfl
3542 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared(
3543 "sink_sink_cfl", "\\Delta t_{\\rm sink-sink}");
3544
3545 using SinkVecEdge = shamrock::solvergraph::IDataEdgeSerializable<std::vector<Tvec>>;
3546 using SinkScalEdge
3547 = shamrock::solvergraph::IDataEdgeSerializable<std::vector<Tscal>>;
3548
3549 ComputeCFLSinkSink<Tvec> compute_cfl_sink_sink{};
3550 compute_cfl_sink_sink.set_edges(
3551 G_edge,
3552 C_force_edge,
3553 eta_phi_edge,
3554 sync.template get_edge_ptr<SinkVecEdge>("sink_pos"),
3555 sync.template get_edge_ptr<SinkScalEdge>("sink_mass"),
3556 sync.template get_edge_ptr<SinkVecEdge>("sink_acc_ext"),
3557 sink_sink_cfl);
3558 compute_cfl_sink_sink.evaluate();
3559
3560 cfl_detail.push_back({"sink_sink", sink_sink_cfl->data});
3561 }
3562
3563 Tscal rank_dt = shambase::get_infty<Tscal>();
3564 for (auto &[key, value] : cfl_detail) {
3565 rank_dt = sham::min(rank_dt, value);
3566 }
3567
3568 if (show_cfl_detail) {
3569 for (auto &[key, value] : cfl_detail) {
3570 value = shamalgs::collective::allreduce_min(value);
3571 }
3572
3573 if (shamcomm::world_rank() == 0) {
3574 shambase::table table(2);
3575 table.add_double_rule();
3576 table.add_data({"key", "value"}, shambase::table::center);
3577 table.add_double_rule();
3578 for (auto &[key, value] : cfl_detail) {
3579 table.add_data(
3580 {key, sham::format("{:.2e}", value)}, shambase::table::right);
3581 }
3582 table.add_rule();
3583 logger::info_ln("sph::Model", "CFL detail :", table.render());
3584 }
3585 }
3586
3587 next_cfl = shamalgs::collective::allreduce_min(rank_dt);
3588
3589 if (shamcomm::world_rank() == 0) {
3590 logger::info_ln(
3591 "sph::Model", "cfl dt =", next_cfl, "cfl multiplier :", get_cfl_multipler());
3592 }
3593
3594 // this should not be needed idealy, but we need the pressure on the ghosts and
3595 // we don't want to communicate it as it can be recomputed from the other fields
3596 // hence we copy the soundspeed at the end of the step to a field in the patchdata
3597 if (solver_config.has_field_soundspeed()) {
3598
3599 const u32 isoundspeed = pdl.get_field_idx<Tscal>("soundspeed");
3600
3601 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
3602 sham::DeviceBuffer<Tscal> &buf_cs = pdat.get_field_buf_ref<Tscal>(isoundspeed);
3603 sham::DeviceBuffer<Tscal> &buf_cs_in
3604 = shambase::get_check_ref(storage.soundspeed)
3605 .get_field(cur_p.id_patch)
3606 .get_buf();
3607
3608 sycl::range range_npart{pdat.get_obj_cnt()};
3609
3611
3612 auto &q = shamsys::instance::get_compute_scheduler().get_queue();
3613 sham::EventList depends_list;
3614
3615 auto cs_in = buf_cs_in.get_read_access(depends_list);
3616 auto cs = buf_cs.get_write_access(depends_list);
3617
3618 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
3619 const Tscal pmass = solver_config.gpart_mass;
3620
3621 cgh.parallel_for(
3622 sycl::range<1>{pdat.get_obj_cnt()}, [=](sycl::item<1> item) {
3623 cs[item] = cs_in[item];
3624 });
3625 });
3626
3627 buf_cs_in.complete_event_state(e);
3628 buf_cs.complete_event_state(e);
3629 });
3630 }
3631
3632 } // if (!need_rerun_corrector) {
3633
3634 corrector_iter_cnt++;
3635
3636 if (solver_config.has_field_alphaAV()) {
3637 storage.alpha_av_ghost.reset();
3638 }
3639 } while (need_rerun_corrector);
3640
3641 reset_merge_ghosts_fields();
3642 reset_eos_fields();
3643
3644 // if delta too big jump to compute force
3645
3646 tstep.stop();
3647
3648 for (auto it = timestep_callbacks.rbegin(); it != timestep_callbacks.rend(); ++it) {
3649 if (it->step_end_callback) {
3650 shambase::get_check_ref(it->step_end_callback)();
3651 }
3652 }
3653
3654 f64 delta_mpi_timer = shamcomm::mpi::get_timer("total") - mpi_timer_start;
3655 sham::MemPerfInfos mem_perf_infos_end = sham::details::get_mem_perf_info();
3656
3658 shamsys::SystemMetrics system_metrics_end = shamsys::get_system_metrics();
3659 shamsys::SystemMetrics system_metrics_delta = system_metrics_end - system_metrics_start;
3660
3661 f64 t_dev_alloc
3662 = (mem_perf_infos_end.time_alloc_device - mem_perf_infos_start.time_alloc_device)
3663 + (mem_perf_infos_end.time_free_device - mem_perf_infos_start.time_free_device);
3664 f64 t_host_alloc = (mem_perf_infos_end.time_alloc_host - mem_perf_infos_start.time_alloc_host)
3665 + (mem_perf_infos_end.time_free_host - mem_perf_infos_start.time_free_host);
3666
3667 u64 rank_count = scheduler().get_rank_count();
3668 f64 rate = f64(rank_count) / tstep.elapsed_sec();
3669
3670 u64 npatch = scheduler().patch_list.local.size();
3671
3672 // logger::info_ln("SPHSolver", "process rate : ", rate, "particle.s-1");
3673
3674 std::string log_step = report_perf_timestep(
3675 rate,
3676 rank_count,
3677 npatch,
3678 tstep.elapsed_sec(),
3679 delta_mpi_timer,
3680 t_dev_alloc,
3681 t_host_alloc,
3682 mem_perf_infos_end.max_allocated_byte_device,
3683 mem_perf_infos_end.max_allocated_byte_host,
3684 system_metrics_delta,
3685 shamsys::has_reporter());
3686
3687 if (shamcomm::world_rank() == 0) {
3688 logger::info_ln("sph::Model", log_step);
3689 logger::info_ln(
3690 "sph::Model", "estimated rate :", dt * (3600 / tstep.elapsed_sec()), "(tsim/hr)");
3691 }
3692
3693 solve_logs.register_log(
3694 {t_current, // f64 solver_t;
3695 dt, // f64 solver_dt;
3696 shamcomm::world_rank(), // i32 world_rank;
3697 rank_count, // u64 rank_count;
3698 rate, // f64 rate;
3699 tstep.elapsed_sec(), // f64 elapsed_sec;
3700 shambase::details::get_wtime(),
3701 system_metrics_delta});
3702
3703 storage.timings_details.reset();
3704
3705 reset_serial_patch_tree();
3706 reset_ghost_handler();
3707
3708 shambase::get_check_ref(storage.part_counts).free_alloc();
3709 shambase::get_check_ref(storage.part_counts_with_ghost).free_alloc();
3710 shambase::get_check_ref(storage.positions_with_ghosts).free_alloc();
3711 shambase::get_check_ref(storage.hpart_with_ghosts).free_alloc();
3712 storage.merged_xyzh.reset();
3713 shambase::get_check_ref(storage.omega).free_alloc();
3714 clear_merged_pos_trees();
3715 clear_ghost_cache();
3716 reset_presteps_rint();
3717 reset_neighbors_cache();
3718
3719 shambase::get_check_ref(storage.neigh_cache).free_alloc();
3720
3721 set_next_dt(next_cfl);
3722 set_time(t_current + dt);
3723
3724 auto get_next_cfl_mult = [&]() {
3725 Tscal cfl_m = get_cfl_multipler();
3726 Tscal stiff = solver_config.cfl_config.cfl_multiplier_stiffness;
3727
3728 return (cfl_m * stiff + 1.) / (stiff + 1.);
3729 };
3730
3731 set_cfl_multipler(get_next_cfl_mult());
3732
3733 TimestepLog log;
3734 log.rank = shamcomm::world_rank();
3735 log.rate = rate;
3736 log.npart = rank_count;
3737 log.tcompute = tstep.elapsed_sec();
3738
3739 return log;
3740}
3741
3742using namespace shammath;
3743
3744template class shammodels::sph::Solver<f64_3, M4>;
3745template class shammodels::sph::Solver<f64_3, M6>;
3746template class shammodels::sph::Solver<f64_3, M8>;
3747
3748template class shammodels::sph::Solver<f64_3, C2>;
3749template class shammodels::sph::Solver<f64_3, C4>;
3750template class shammodels::sph::Solver<f64_3, C6>;
ComputeCFLSinkSink.hpp
Host-side CFL condition from the pairwise potential between sink particles.
ComputeGravWave.hpp
Compute the gravitational wave quadrupole. Based on Toscani et. al. 2021.
ComputeNeighStats.hpp
A module to compute and display statistics on neighbor counts for SPH particles.
CopyPatchDataFieldFromLayer.hpp
Defines the CopyPatchDataFieldFromLayer class for copying fields between patch data layers.
DistributedBuffers.hpp
Defines the DistributedBuffers class for managing distributed device buffers in a solver graph.
ForwardEulerHost2Deriv.hpp
Host-side (std::vector) forward Euler integration node with two derivative contributions.
ForwardEulerHost.hpp
Host-side (std::vector) forward Euler integration node.
ForwardEulerPositive.hpp
Implements a forward Euler integration step as a solver graph node.
ForwardEuler.hpp
Implements a forward Euler integration step as a solver graph node.
GetFieldRefFromLayer.hpp
Defines the GetFieldRefFromLayer class for extracting field references from patch data layers.
GetObjCntFromLayer.hpp
Defines the GetObjCntFromLayer class for extracting object counts from patch data layers.
GetParticlesOutsideSphere.hpp
Declares the GetParticlesOutsideSphere module for removing particles.
shamrock::solvergraph::DDPatchDataFieldRef
shambase::DistributedData< PatchDataFieldRef< T > > DDPatchDataFieldRef
Alias for a DistributedData of PatchDataFieldRefs.
Definition IFieldRefs.hpp:32
IterateSmoothingLengthDensityNeighLim.hpp
Declares the IterateSmoothingLengthDensityNeighLim module for iterating smoothing length based on the...
IterateSmoothingLengthDensity.hpp
Declares the IterateSmoothingLengthDensity module for iterating smoothing length based on the SPH den...
KillParticles.hpp
Declares the KillParticles module for removing particles.
LoopSmoothingLengthIter.hpp
Declares the LoopSmoothingLengthIter module for looping over the smoothing length iteration until con...
NodeFreeAlloc.hpp
Field variant object to instanciate a variant on the patch types.
NodeInstance.hpp
Header file describing a Node Instance.
NodeMapEdge.hpp
Node that maps a read-only input edge into a read-write output edge.
NodeSetEdge.hpp
Node that applies a custom function to modify connected edges.
OperationIf.hpp
Meta node that evaluates an optional then-node when a bool edge is true, and an optional else-node wh...
PatchDataLayerRefs.hpp
Defines the PatchDataLayerRefs class for managing distributed references to patch data layers.
PatchScheduler.hpp
MPI scheduler.
Patch.hpp
Header file for the patch struct and related function.
SinkParticlesAccreteQuantities.hpp
Accrete flagged SPH particles onto sinks (mass, CoM, spin, etc.).
SinkParticlesEvictAccretedParticles.hpp
Remove SPH particles flagged for sink accretion from patch data.
SinkParticlesFlagAccreteHard.hpp
Flag SPH particles inside sink accretion radii into an accretion table.
SinkSelfGravityHost.hpp
Host-side pairwise gravitational self-interaction between sink particles.
SolverGraph.hpp
Declare a class to register and retrieve nodes and edges from a unique container.
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
i32
std::int32_t i32
32 bit integer
Definition aliases_int.hpp:23
assert.hpp
Shamrock assertion utility.
SHAM_ASSERT
#define SHAM_ASSERT(x)
Shorthand for SHAM_ASSERT_NAMED without a message.
Definition assert.hpp:67
ComputeCFLSinkSink
Host-side (N^2) sink-sink CFL condition.
Definition ComputeCFLSinkSink.hpp:48
PatchScheduler
The MPI scheduler.
Definition PatchScheduler.hpp:86
PatchScheduler::synchronized_data
SynchronizedData synchronized_data
data that is synchroneous across all ranks
Definition PatchScheduler.hpp:106
sham::DeviceBuffer
A buffer allocated in USM (Unified Shared Memory).
Definition DeviceBuffer.hpp:106
sham::DeviceBuffer::complete_event_state
void complete_event_state(sycl::event e) const
Complete the event state of the buffer.
Definition DeviceBuffer.hpp:368
sham::DeviceBuffer::get_queue
DeviceQueue & get_queue() const
Gets the DeviceQueue associated with the held allocation.
Definition DeviceBuffer.hpp:457
sham::DeviceBuffer::get_write_access
T * get_write_access(sham::EventList &depends_list, SourceLocation src_loc=SourceLocation{})
Get a read-write pointer to the buffer's data.
Definition DeviceBuffer.hpp:349
sham::DeviceBuffer::fill
void fill(T value, std::array< size_t, 2 > idx_range)
Fill a subpart of the buffer with a given value.
Definition DeviceBuffer.hpp:951
sham::DeviceBuffer::get_read_access
const T * get_read_access(sham::EventList &depends_list, SourceLocation src_loc=SourceLocation{}) const
Get a read-only pointer to the buffer's data.
Definition DeviceBuffer.hpp:331
sham::EventList
Class to manage a list of SYCL events.
Definition EventList.hpp:32
sham::EventList::add_event
void add_event(sycl::event e)
Add an event to the list of events.
Definition EventList.hpp:88
shambase::DistributedDataShared
Container for objects shared between two distributed data elements.
Definition DistributedDataShared.hpp:52
shambase::DistributedDataShared::for_each
void for_each(std::function< void(u64, u64, T &)> &&f)
Apply a function to all stored objects.
Definition DistributedDataShared.hpp:110
shambase::DistributedData
Represents a collection of objects distributed across patches identified by a u64 id.
Definition DistributedData.hpp:43
shambase::DistributedData::add_obj
iterator add_obj(u64 id, T &&obj)
Adds a new object to the collection.
Definition DistributedData.hpp:67
shambase::DistributedData::map
DistributedData< Tmap > map(std::function< Tmap(u64, T &)> map_func)
Apply a function to all objects in the collection and return a new collection containing the results.
Definition DistributedData.hpp:250
shambase::DistributedData::get
T & get(u64 id)
Returns a reference to an object in the collection.
Definition DistributedData.hpp:145
shambase::Timer
Class Timer measures the time elapsed since the timer was started.
Definition Timer.hpp:36
shambase::Timer::elapsed_sec
f64 elapsed_sec() const
Converts the stored nanosecond time to a floating point representation in seconds.
Definition Timer.hpp:88
shambase::Timer::start
void start()
Starts the timer.
Definition Timer.hpp:51
shambase::Timer::stop
void stop()
Stops the timer and stores the elapsed time in nanoseconds.
Definition Timer.hpp:65
shammath::vec
Vector class based on std::array storage and mdspan.
Definition matrix.hpp:98
shammodels::sph::SPHUtilities
handle basic utilities dealing with SPH
Definition SPHUtilities.hpp:65
shammodels::sph::Solver
The shamrock SPH model.
Definition Solver.hpp:151
shammodels::sph::Solver::reset_presteps_rint
void reset_presteps_rint()
Resets tree radius interval field.
Definition Solver.cpp:1884
shammodels::sph::Solver::reset_merge_ghosts_fields
void reset_merge_ghosts_fields()
Resets merged ghost field data.
Definition Solver.cpp:2160
shammodels::sph::Solver::update_sync_load_values
void update_sync_load_values()
Updates load balancing values and synchronizes patch ownership.
Definition Solver.cpp:2459
shammodels::sph::Solver::apply_corrector
bool apply_corrector(Tscal dt, u64 Npart_all)
Definition Solver.cpp:2454
shammodels::sph::Solver::merge_position_ghost
void merge_position_ghost()
Merges ghost particle positions from neighboring patches.
Definition Solver.cpp:1529
shammodels::sph::Solver::reset_eos_fields
void reset_eos_fields()
Frees memory allocated for EOS fields.
Definition Solver.cpp:2186
shammodels::sph::Solver::prepare_corrector
void prepare_corrector()
Saves old derivative fields for predictor-corrector integration.
Definition Solver.cpp:2192
shammodels::sph::Solver::build_ghost_cache
void build_ghost_cache()
Builds ghost particle interface cache for inter-patch communication.
Definition Solver.cpp:1507
shammodels::sph::Solver::update_artificial_viscosity
void update_artificial_viscosity(Tscal dt)
Updates artificial viscosity coefficients for shock capturing.
Definition Solver.cpp:2169
shammodels::sph::Solver::evolve_once
TimestepLog evolve_once()
Performs one complete SPH timestep evolution.
Definition Solver.cpp:2466
shammodels::sph::Solver::vtk_do_dump
void vtk_do_dump(std::string filename, bool add_patch_world_id)
Writes VTK dump file for visualization.
Definition Solver.cpp:1284
shammodels::sph::Solver::update_derivs
void update_derivs(Tscal dt_hydro)
Updates time derivatives and applies external forces.
Definition Solver.cpp:2288
shammodels::sph::Solver::build_merged_pos_trees
void build_merged_pos_trees()
Builds spatial BVH trees for merged positions including ghosts.
Definition Solver.cpp:1572
shammodels::sph::Solver::clear_merged_pos_trees
void clear_merged_pos_trees()
Clears merged position trees to free memory.
Definition Solver.cpp:1577
shammodels::sph::Solver::init_solver_graph
void init_solver_graph()
Initializes the solver graph for computation pipeline.
Definition Solver.cpp:263
shammodels::sph::Solver::sph_prestep
void sph_prestep(Tscal time_val, Tscal dt)
Performs pre-step operations for SPH timestep.
Definition Solver.cpp:1583
shammodels::sph::Solver::compute_presteps_rint
void compute_presteps_rint()
Computes maximum smoothing length in tree nodes for neighbor search.
Definition Solver.cpp:1847
shammodels::sph::Solver::compute_eos_fields
void compute_eos_fields()
Computes equation of state fields (pressure, sound speed).
Definition Solver.cpp:2180
shammodels::sph::Solver::apply_position_boundary
void apply_position_boundary(Tscal time_val)
Applies position-based boundary conditions.
Definition Solver.cpp:1461
shammodels::sph::Solver::reset_neighbors_cache
void reset_neighbors_cache()
Resets neighbor cache.
Definition Solver.cpp:1913
shammodels::sph::Solver::communicate_merge_ghosts_fields
void communicate_merge_ghosts_fields()
Communicates and merges ghost particle fields across processes.
Definition Solver.cpp:1918
shammodels::sph::Solver::clear_ghost_cache
void clear_ghost_cache()
Clears ghost particle cache to free memory.
Definition Solver.cpp:1523
shammodels::sph::Solver::init_ghost_layout
void init_ghost_layout()
Initializes data layout for ghost particle fields.
Definition Solver.cpp:1832
shammodels::sph::Solver::start_neighbors_cache
void start_neighbors_cache()
Builds neighbor particle cache for SPH calculations.
Definition Solver.cpp:1889
shammodels::sph::modules::BuildTrees
Module for constructing spatial tree structures for SPH neighbor searches.
Definition BuildTrees.hpp:33
shammodels::sph::modules::BuildTrees::build_merged_pos_trees
void build_merged_pos_trees()
Builds compressed leaf BVH trees for merged particle positions including ghosts.
Definition BuildTrees.cpp:27
shammodels::sph::modules::ComputeEos
Module for computing equation of state quantities.
Definition ComputeEos.hpp:39
shammodels::sph::modules::ComputeEos::compute_eos
void compute_eos()
Computes pressure and sound speed from equation of state.
Definition ComputeEos.cpp:1138
shammodels::sph::modules::ConservativeCheck
Module for checking conservation of physical quantities.
Definition ConservativeCheck.hpp:33
shammodels::sph::modules::ConservativeCheck::check_conservation
void check_conservation()
Verifies conservation of mass, momentum, and energy.
Definition ConservativeCheck.cpp:26
shammodels::sph::modules::ExternalForces::add_ext_forces
void add_ext_forces()
add external forces to the particle acceleration, note that forces dependant on velocity shlould be a...
Definition ExternalForces.cpp:333
shammodels::sph::modules::ExternalForces::compute_ext_forces_indep_v
void compute_ext_forces_indep_v()
is ran once per timestep, it computes the forces that are independant of velocity
Definition ExternalForces.cpp:48
shammodels::sph::modules::ParticleReordering
Module for reordering particles to improve cache locality.
Definition ParticleReordering.hpp:34
shammodels::sph::modules::ParticleReordering::reorder_particles
void reorder_particles()
Reorders particles by Morton code for improved memory access patterns.
Definition ParticleReordering.cpp:22
shammodels::sph::modules::SinkSelfGravityHost
Host-side pairwise (N^2) gravitational self-interaction between sink particles.
Definition SinkSelfGravityHost.hpp:42
shammodels::sph::modules::VTKDump
Module for writing VTK format output files.
Definition VTKDump.hpp:33
shammodels::sph::modules::VTKDump::do_dump
void do_dump(std::string filename, bool add_patch_world_id)
Writes particle data to VTK file for visualization.
Definition VTKDump.cpp:37
shamrock::ReattributeDataUtility
Utility class used to move the objects between patches.
Definition ReattributeDataUtility.hpp:38
shamrock::ReattributeDataUtility::reatribute_patch_objects
void reatribute_patch_objects(SerialPatchTree< T > &sptree, std::string position_field)
Reattribute objects based on a given position field.
Definition ReattributeDataUtility.hpp:216
shamrock::SchedulerUtility::make_compute_field
ComputeField< T > make_compute_field(std::string new_name, u32 nvar)
create a compute field and init it to zeros
Definition SchedulerUtility.hpp:264
shamrock::SchedulerUtility::save_field
ComputeField< T > save_field(u32 field_idx, std::string new_name)
save a field in patchdata to a compute field
Definition SchedulerUtility.hpp:231
shamrock::patch::PatchDataLayerLayout::get_field_idx
u32 get_field_idx(const std::string &field_name) const
Get the field id if matching name & type.
Definition PatchDataLayerLayout.cpp:130
shamrock::patch::PatchDataLayer
PatchDataLayer container class, the layout is described in patchdata_layout.
Definition PatchDataLayer.hpp:35
shamrock::patch::PatchDataLayer::get_obj_cnt
u32 get_obj_cnt() const
get the number of objects (particles) stored in this layer
Definition PatchDataLayer.cpp:346
shamrock::solvergraph::DistributedBuffers
Interface for a solver graph edge representing a field as spans.
Definition DistributedBuffers.hpp:39
shamrock::solvergraph::ForwardEulerHost2Deriv
Forward Euler integration of a host-side (std::vector) field with two derivative contributions.
Definition ForwardEulerHost2Deriv.hpp:42
shamrock::solvergraph::ForwardEulerHost
Forward Euler integration of a host-side (std::vector) field.
Definition ForwardEulerHost.hpp:39
shamrock::solvergraph::IFieldRefs::get_field
PatchDataField< T > & get_field(u64 id) const
Get the underlying PatchDataField at the given id.
Definition IFieldRefs.hpp:54
shamrock::solvergraph::INode::evaluate
void evaluate()
Evaluate the node.
Definition INode.hpp:156
shamrock::solvergraph::NodeFreeAlloc
A node that simply frees the allocation of the connected node.
Definition NodeFreeAlloc.hpp:32
shamrock::solvergraph::NodeMapEdge
A node that maps an input edge into an output edge.
Definition NodeMapEdge.hpp:46
shamrock::solvergraph::NodeSetEdge
A node that applies a custom function to modify connected edges.
Definition NodeSetEdge.hpp:46
shamrock::solvergraph::NodeSetEdge::set_edges
void set_edges(std::shared_ptr< IEdge > to_set)
Set the edges of the node.
Definition NodeSetEdge.hpp:67
shamrock::solvergraph::OperationIf
Conditional meta node: if condition is true, evaluate then_node when one was provided; otherwise eval...
Definition OperationIf.hpp:34
shamrock::solvergraph::PatchDataLayerRefs::free_alloc
virtual void free_alloc() override
Free allocated memory.
Definition PatchDataLayerRefs.hpp:50
shamrock::solvergraph::SolverGraph
A graph container for managing solver nodes and edges with type-safe access.
Definition SolverGraph.hpp:110
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:202
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:391
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:341
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:321
shamrock::solvergraph::SolverGraph::get_node_ref_base
INode & get_node_ref_base(const std::string &name)
Get a reference to a node by name through the base interface.
Definition SolverGraph.hpp:278
shamtree::CompressedLeafBVH
A Compressed Leaf Bounding Volume Hierarchy (CLBVH) for neighborhood queries.
Definition CompressedLeafBVH.hpp:46
shamtree::KarrasRadixTreeField
A data structure representing a Karras Radix Tree Field.
Definition KarrasRadixTreeField.hpp:57
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::details::get_mem_perf_info
MemPerfInfos get_mem_perf_info()
Retrieve the memory performance information.
Definition internal_alloc.cpp:123
kernel_call.hpp
memoryHandle.hpp
This file contains the declaration of the memory handling and its methods.
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:344
shamalgs::memory::buf_to_vec
std::vector< T > buf_to_vec(sycl::buffer< T > &buf, u32 len)
Convert a sycl::buffer to a std::vector.
Definition memory.cpp:34
shambase
namespace for basic c++ utilities
Definition sycl_builtins.hpp:25
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:112
shambase::make_except_with_loc
ExcptTypes make_except_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Create an exception with a message and a location.
Definition exception.hpp:105
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:41
shammath
namespace for math utility
Definition AABB.hpp:26
shammodels::sph
namespace for the sph model
Definition BasicSPHGhosts.hpp:37
shammodels::sph::has_sinks
bool has_sinks(shamrock::solvergraph::SolverGraphSerializable &sync)
Check whether any sinks are present by inspecting sink_pos only.
Definition sink_edges_helper.hpp:152
shammodels::sph::get_sink_pos
std::vector< Tvec > & get_sink_pos(shamrock::solvergraph::SolverGraphSerializable &sync)
Named SoA getters (edges must already exist; call ensure_sink_edges first). Prefer these when a funct...
Definition sink_edges_helper.hpp:105
shammodels::SingleStage
@ SingleStage
Single tree traversal per particle.
Definition enum_NeighCacheStrategy.hpp:32
shammodels::TwoStage
@ TwoStage
Two stage neighbours search (see shamrock paper).
Definition enum_NeighCacheStrategy.hpp:33
shamrock
namespace for the main framework
Definition __init__.py:1
numeric_limits.hpp
q_ab.hpp
shamcomm::logs::info
void info(std::string module_name, Types... var2)
Prints a log message with multiple arguments.
Definition logs.hpp:132
shamcomm::logs::raw_ln
void raw_ln(Types... var2)
Prints a log message with multiple arguments followed by a newline.
Definition logs.hpp:89
shamcomm::logs::info_ln
void info_ln(std::string module_name, Types... var2)
Prints a log message with multiple arguments followed by a newline.
Definition logs.hpp:132
shamcomm::logs::warn_ln
void warn_ln(std::string module_name, Types... var2)
Prints a log message with multiple arguments followed by a newline.
Definition logs.hpp:132
shamcomm::logs::err_ln
void err_ln(std::string module_name, Types... var2)
Prints a log message with multiple arguments followed by a newline.
Definition logs.hpp:132
sink_edges_helper.hpp
Helpers to access SPH sink particles stored as SoA synchronized data edges.
forces.hpp
file containing formulas for sph forces
sphkernels.hpp
sph kernels
NamedStackEntry
shambase::details::NamedBasicStackEntry NamedStackEntry
Alias for shambase::details::NamedBasicStackEntry.
Definition stacktrace.hpp:249
StackEntry
shambase::details::BasicStackEntry StackEntry
Alias for shambase::details::BasicStackEntry.
Definition stacktrace.hpp:242
shambase::details::get_wtime
f64 get_wtime()
Returns the current wall clock time in seconds.
Definition stacktrace.cpp:182
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 MultiRef.hpp:33
shammodels::sph::Debug_ph_dump
Definition Solver.cpp:1297
shammodels::sph::PhantomDumpBlock
A class to represent a single block of data in a Phantom dump.
Definition PhantomDump.hpp:185
shammodels::sph::PhantomDumpBlock::get_ref_f32
u64 get_ref_f32(std::string s)
Gets the index of a block of type f32 with the given name.
Definition PhantomDump.cpp:243
shammodels::sph::PhantomDumpBlock::get_ref_fort_real
u64 get_ref_fort_real(std::string s)
Gets the index of a block of type fort_real with the given name.
Definition PhantomDump.cpp:219
shammodels::sph::PhantomDumpBlock::tot_count
i64 tot_count
The total number of values in the block.
Definition PhantomDump.hpp:187
shammodels::sph::PhantomDumpBlock::blocks_fort_real
std::vector< PhantomDumpBlockArray< fort_real > > blocks_fort_real
The blocks of values of type fort_real.
Definition PhantomDump.hpp:205
shammodels::sph::PhantomDumpBlock::blocks_f32
std::vector< PhantomDumpBlockArray< f32 > > blocks_f32
The blocks of values of type f32.
Definition PhantomDump.hpp:207
shammodels::sph::PhantomDump
Class representing a Phantom dump file.
Definition PhantomDump.hpp:347
shammodels::sph::PhantomDump::override_magic_number
void override_magic_number()
Overrides the magic numbers used in the PhantomDump struct.
Definition PhantomDump.hpp:368
shammodels::sph::SolverConfig::BCConfig
BCConfig< Tvec > BCConfig
Configuration of the boundary conditions.
Definition SolverConfig.hpp:956
shamrock::patch::Patch
Patch object that contain generic patch information.
Definition Patch.hpp:33
shamrock::patch::Patch::id_patch
u64 id_patch
unique key that identify the patch
Definition Patch.hpp:86
system_metrics.hpp
timestep_report.hpp
worldInfo.hpp
Functions related to the MPI communicator.
shamcomm::mpi::get_timer
f64 get_timer(std::string timername)
get a timer value
Definition wrapper.cpp:46
Generated on for Shamrock by doxygen 1.15.0