Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
UpdateDerivs.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
17
18#include "shambase/DistributedData.hpp"
19#include "shambase/memory.hpp"
20#include "shambackends/kernel_call.hpp"
21#include "shambackends/kernel_call_distrib.hpp"
22#include "shambackends/math.hpp"
23#include "shamcomm/logs.hpp"
24#include "shammath/sphkernels.hpp"
25#include "shammodels/sph/SolverConfig.hpp"
26#include "shammodels/sph/math/density.hpp"
27#include "shammodels/sph/math/forces.hpp"
28#include "shammodels/sph/math/mhd.hpp"
29#include "shammodels/sph/math/q_ab.hpp"
30#include "shammodels/sph/modules/ComputeDustTtilde.hpp"
31#include "shammodels/sph/modules/MonoFluidTVADeltav.hpp"
32#include "shammodels/sph/modules/NodeComputePressureGrad.hpp"
33#include "shammodels/sph/modules/NodeEvolveDustCOALASourceTerm.hpp"
34#include "shammodels/sph/modules/NodeMonofluidTVAAddSourceTerm.hpp"
35#include "shammodels/sph/modules/NodeMonofluidTVASmoothSPositivityLimiter.hpp"
36#include "shammodels/sph/modules/NodeUpdateDerivsMonofluidTVA.hpp"
37#include "shammodels/sph/modules/NodeUpdateDerivsVaryingAlphaAV.hpp"
38#include "shammodels/sph/modules/NodeUpdateDerivsVaryingAlphaAVDustTVA.hpp"
39#include "shammodels/sph/modules/SetDustStoppingTimeConstant.hpp"
40#include "shammodels/sph/modules/SetDustStoppingTimeEpstein.hpp"
41#include "shammodels/sph/modules/UpdateDerivs.hpp"
42#include "shamphys/mhd.hpp"
43#include "shamrock/patch/PatchDataFieldSpan.hpp"
44#include "shamrock/solvergraph/FieldRefs.hpp"
45#include "shamrock/solvergraph/IFieldSpan.hpp"
46#include "shamrock/solvergraph/Indexes.hpp"
47#include "shamrock/solvergraph/ScalarEdge.hpp"
48#include <memory>
49#include <vector>
50
51template<class Tvec, template<class> class SPHKernel>
52void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs(Tscal dt_hydro) {
53
54 Cfg_AV cfg_av = solver_config.artif_viscosity;
55 Cfg_MHD cfg_mhd = solver_config.mhd_config;
56 DustConfig cfg_dust = solver_config.dust_config;
57
58 if (Constant *v = std::get_if<Constant>(&cfg_av.config)) {
59 update_derivs_constantAV(*v);
60 } else if (VaryingMM97 *v = std::get_if<VaryingMM97>(&cfg_av.config)) {
61 update_derivs_mm97(*v);
62 } else if (VaryingCD10 *v = std::get_if<VaryingCD10>(&cfg_av.config)) {
63 update_derivs_cd10(*v);
64 } else if (ConstantDisc *v = std::get_if<ConstantDisc>(&cfg_av.config)) {
65 update_derivs_disc_visco(*v);
66 } else if (IdealMHD *v = std::get_if<IdealMHD>(&cfg_mhd.config)) {
67 update_derivs_MHD(*v);
68 } else if (NonIdealMHD *v = std::get_if<NonIdealMHD>(&cfg_mhd.config)) {
69 shambase::throw_unimplemented();
70 } else if (NoneMHD *v = std::get_if<NoneMHD>(&cfg_mhd.config)) {
71 shambase::throw_unimplemented();
72 } else if (None *v = std::get_if<None>(&cfg_av.config)) {
73 shambase::throw_unimplemented();
74 } else {
75 shambase::throw_unimplemented();
76 }
77
78 if (cfg_dust.has_s_j_field()) {
79 // we can do it separately because the backreaction is done only through the pressure
80 update_derivs_dust_monofluid_tva_Sj(cfg_dust, dt_hydro);
81 }
82}
83
84template<class Tvec, template<class> class SPHKernel>
85void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_noAV(None cfg) {}
86
87template<class Tvec, template<class> class SPHKernel>
88void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_constantAV(
89 Constant cfg) {
90 StackEntry stack_loc{};
91
92 using namespace shamrock;
93 using namespace shamrock::patch;
94
95 PatchDataLayerLayout &pdl = scheduler().pdl_old();
96
97 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
98 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
99 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
100 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
101 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
102 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
103
104 shamrock::patch::PatchDataLayerLayout &ghost_layout
105 = shambase::get_check_ref(storage.ghost_layout.get());
106 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
107 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
108 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
109 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
110
111 auto &merged_xyzh = storage.merged_xyzh.get();
112 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
113 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
114
115 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
116 PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch);
117
118 sham::DeviceBuffer<Tvec> &buf_xyz
119 = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref<Tvec>(0);
120 sham::DeviceBuffer<Tvec> &buf_axyz = pdat.get_field_buf_ref<Tvec>(iaxyz);
121 sham::DeviceBuffer<Tscal> &buf_duint = pdat.get_field_buf_ref<Tscal>(iduint);
122 sham::DeviceBuffer<Tvec> &buf_vxyz = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);
123 sham::DeviceBuffer<Tscal> &buf_hpart = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);
124 sham::DeviceBuffer<Tscal> &buf_omega = mpdat.get_field_buf_ref<Tscal>(iomega_interf);
125 sham::DeviceBuffer<Tscal> &buf_uint = mpdat.get_field_buf_ref<Tscal>(iuint_interf);
126 sham::DeviceBuffer<Tscal> &buf_pressure
127 = shambase::get_check_ref(storage.pressure).get_field(cur_p.id_patch).get_buf();
128 sham::DeviceBuffer<Tscal> &buf_cs
129 = shambase::get_check_ref(storage.soundspeed).get_field(cur_p.id_patch).get_buf();
130
131 sycl::range range_npart{pdat.get_obj_cnt()};
132
133 tree::ObjectCache &pcache
134 = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);
135
137
138 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
139 sham::EventList depends_list;
140
141 auto xyz = buf_xyz.get_read_access(depends_list);
142 auto axyz = buf_axyz.get_write_access(depends_list);
143 auto du = buf_duint.get_write_access(depends_list);
144 auto vxyz = buf_vxyz.get_read_access(depends_list);
145 auto hpart = buf_hpart.get_read_access(depends_list);
146 auto omega = buf_omega.get_read_access(depends_list);
147 auto u = buf_uint.get_read_access(depends_list); // TODO rename to uint
148 auto pressure = buf_pressure.get_read_access(depends_list);
149 auto cs = buf_cs.get_read_access(depends_list);
150 auto ploop_ptrs = pcache.get_read_access(depends_list);
151
152 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
153 const Tscal pmass = solver_config.gpart_mass;
154 const Tscal alpha_u = cfg.alpha_u;
155 const Tscal alpha_AV = cfg.alpha_AV;
156 const Tscal beta_AV = cfg.beta_AV;
157
158 shamlog_debug_sycl_ln("deriv kernel", "alpha_u :", alpha_u);
159 shamlog_debug_sycl_ln("deriv kernel", "alpha_AV :", alpha_AV);
160 shamlog_debug_sycl_ln("deriv kernel", "beta_AV :", beta_AV);
161
162 // tree::ObjectIterator particle_looper(tree,cgh);
163
164 // tree::LeafCacheObjectIterator
165 // particle_looper(tree,*xyz_cell_id,leaf_cache,cgh);
166
167 tree::ObjectCacheIterator particle_looper(ploop_ptrs);
168
169 // sycl::accessor hmax_tree{tree_field_hmax, cgh, sycl::read_only};
170
171 // sycl::stream out {4096,1024,cgh};
172
173 constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;
174
175 shambase::parallel_for(cgh, pdat.get_obj_cnt(), "compute force cte AV", [=](u64 gid) {
176 u32 id_a = (u32) gid;
177
178 using namespace shamrock::sph;
179
180 Tvec sum_axyz = {0, 0, 0};
181 Tscal sum_du_a = 0;
182
183 Tscal h_a = hpart[id_a];
184 Tvec xyz_a = xyz[id_a];
185 Tvec vxyz_a = vxyz[id_a];
186 Tscal P_a = pressure[id_a];
187 Tscal omega_a = omega[id_a];
188 const Tscal u_a = u[id_a];
189
190 Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd);
191 Tscal rho_a_sq = rho_a * rho_a;
192 Tscal rho_a_inv = 1. / rho_a;
193
194 // f32 P_a = cs * cs * rho_a;
195
196 Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a);
197
198 Tscal cs_a = cs[id_a];
199
200 Tvec force_pressure{0, 0, 0};
201 Tscal tmpdU_pressure = 0;
202
203 particle_looper.for_each_object(id_a, [&](u32 id_b) {
204 // compute only omega_a
205 Tvec dr = xyz_a - xyz[id_b];
206 Tscal rab2 = sycl::dot(dr, dr);
207 Tscal h_b = hpart[id_b];
208
209 if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) {
210 return;
211 }
212
213 Tscal rab = sycl::sqrt(rab2);
214 Tvec vxyz_b = vxyz[id_b];
215 const Tscal u_b = u[id_b];
216
217 Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd);
218 Tscal P_b = pressure[id_b];
219 // f32 P_b = cs * cs * rho_b;
220 Tscal omega_b = omega[id_b];
221 Tscal cs_b = cs[id_b];
222
223 const Tscal alpha_a = alpha_AV;
224 const Tscal alpha_b = alpha_AV;
225
226 Tscal Fab_a = Kernel::dW_3d(rab, h_a);
227 Tscal Fab_b = Kernel::dW_3d(rab, h_b);
228
229 Tvec v_ab = vxyz_a - vxyz_b;
230
231 Tvec r_ab_unit = dr * sham::inv_sat_positive(rab);
232
233 // f32 P_b = cs * cs * rho_b;
234 Tscal v_ab_r_ab = sycl::dot(v_ab, r_ab_unit);
235 Tscal abs_v_ab_r_ab = sycl::fabs(v_ab_r_ab);
236
237 Tscal vsig_a = alpha_a * cs_a + beta_AV * abs_v_ab_r_ab;
238 Tscal vsig_b = alpha_b * cs_b + beta_AV * abs_v_ab_r_ab;
239
240 Tscal vsig_u = shamrock::sph::vsig_u(P_a, P_b, rho_a, rho_b);
241
242 Tscal qa_ab = shamrock::sph::q_av(rho_a, vsig_a, v_ab_r_ab);
243 Tscal qb_ab = shamrock::sph::q_av(rho_b, vsig_b, v_ab_r_ab);
244
245 add_to_derivs_sph_artif_visco_cond(
246 pmass,
247 rho_a_sq,
248 omega_a_rho_a_inv,
249 rho_a_inv,
250 rho_b,
251 omega_a,
252 omega_b,
253 Fab_a,
254 Fab_b,
255 u_a,
256 u_b,
257 P_a,
258 P_b,
259 alpha_u,
260 v_ab,
261 r_ab_unit,
262 vsig_u,
263 qa_ab,
264 qb_ab,
265 force_pressure,
266 tmpdU_pressure);
267 });
268 axyz[id_a] = force_pressure;
269 du[id_a] = tmpdU_pressure;
270 });
271 });
272
273 buf_xyz.complete_event_state(e);
274 buf_axyz.complete_event_state(e);
275 buf_duint.complete_event_state(e);
276 buf_vxyz.complete_event_state(e);
277 buf_hpart.complete_event_state(e);
278 buf_omega.complete_event_state(e);
279 buf_uint.complete_event_state(e);
280 buf_pressure.complete_event_state(e);
281 buf_cs.complete_event_state(e);
282
283 sham::EventList resulting_events;
284 resulting_events.add_event(e);
285 pcache.complete_event_state(resulting_events);
286 });
287}
288template<class Tvec, template<class> class SPHKernel>
289void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_mm97(VaryingMM97 cfg) {
290 StackEntry stack_loc{};
291
292 using namespace shamrock;
293 using namespace shamrock::patch;
294
295 PatchDataLayerLayout &pdl = scheduler().pdl_old();
296
297 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
298 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
299 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
300 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
301 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
302 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
303
304 shamrock::patch::PatchDataLayerLayout &ghost_layout
305 = shambase::get_check_ref(storage.ghost_layout.get());
306 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
307 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
308 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
309 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
310
311 auto &merged_xyzh = storage.merged_xyzh.get();
312 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
313 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
314
315 auto &part_counts = storage.part_counts;
316 auto &part_counts_with_ghost = storage.part_counts_with_ghost;
317 auto &xyz_refs = storage.positions_with_ghosts;
318 auto &pressure_field = storage.pressure;
319 auto &soundspeed_field = storage.soundspeed;
320
321 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> uint_refs
322 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("uint", "u");
323 {
324 shambase::get_check_ref(uint_refs).set_refs(
325 mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
326 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
327 return std::ref(mpdat.get_field<Tscal>(iuint_interf));
328 }));
329 }
330
331 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> vxyz_refs
332 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("vxyz", "v");
333 {
334 shambase::get_check_ref(vxyz_refs).set_refs(
335 mpdats.map<std::reference_wrapper<PatchDataField<Tvec>>>(
336 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
337 return std::ref(mpdat.get_field<Tvec>(ivxyz_interf));
338 }));
339 }
340
341 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hpart_refs
342 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("hpart", "h");
343 { // if was just reset before this call
344 shambase::get_check_ref(hpart_refs)
345 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
346 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
347 return std::ref(mpdat.get_field<Tscal>(ihpart_interf));
348 }));
349 }
350
351 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> omega_refs
352 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("omega", "omega");
353 {
354 shambase::get_check_ref(omega_refs)
355 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
356 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
357 return std::ref(mpdat.get_field<Tscal>(iomega_interf));
358 }));
359 }
360
361 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> alpha_av_refs
362 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("alpha_av", "alpha_av");
363 {
364 shambase::DistributedData<std::reference_wrapper<PatchDataField<Tscal>>> refs{};
365 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
366 refs.add_obj(
367 cur_p.id_patch, std::ref(storage.alpha_av_ghost.get().get(cur_p.id_patch)));
368 });
369 shambase::get_check_ref(alpha_av_refs).set_refs(refs);
370 }
371
372 shamrock::solvergraph::SolverGraph &solver_graph = storage.solver_graph;
373
374 auto axyz_refs = solver_graph.get_edge_ptr<shamrock::solvergraph::FieldRefs<Tvec>>("axyz");
375 auto duint_refs = solver_graph.get_edge_ptr<shamrock::solvergraph::FieldRefs<Tscal>>("duint");
376 auto gpart_mass
377 = solver_graph.get_edge_ptr<shamrock::solvergraph::ScalarEdge<Tscal>>("gpart_mass");
378
379 std::shared_ptr<shamrock::solvergraph::ScalarEdge<Tscal>> alpha_u
380 = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>("alpha_u", "alpha_u");
381 {
382 shambase::get_check_ref(alpha_u).value = cfg.alpha_u;
383 }
384 std::shared_ptr<shamrock::solvergraph::ScalarEdge<Tscal>> beta_AV
385 = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>("beta_AV", "beta_AV");
386 {
387 shambase::get_check_ref(beta_AV).value = cfg.beta_AV;
388 }
389
390 std::shared_ptr<NodeUpdateDerivsVaryingAlphaAV<Tvec, SPHKernel>> node
391 = std::make_shared<NodeUpdateDerivsVaryingAlphaAV<Tvec, SPHKernel>>();
392 {
393 node->set_edges(
394 gpart_mass,
395 alpha_u,
396 beta_AV,
397 part_counts,
398 part_counts_with_ghost,
399 xyz_refs,
400 hpart_refs,
401 vxyz_refs,
402 uint_refs,
403 omega_refs,
404 pressure_field,
405 soundspeed_field,
406 alpha_av_refs,
407 storage.neigh_cache,
408 axyz_refs,
409 duint_refs);
410 }
411 node->evaluate();
412}
413template<class Tvec, template<class> class SPHKernel>
414void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_cd10(VaryingCD10 cfg) {
415 StackEntry stack_loc{};
416
417 using namespace shamrock;
418 using namespace shamrock::patch;
419
420 PatchDataLayerLayout &pdl = scheduler().pdl_old();
421
422 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
423 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
424 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
425 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
426 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
427 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
428
429 shamrock::patch::PatchDataLayerLayout &ghost_layout
430 = shambase::get_check_ref(storage.ghost_layout.get());
431 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
432 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
433 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
434 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
435
436 auto &merged_xyzh = storage.merged_xyzh.get();
437 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
438 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
439
440 auto &part_counts = storage.part_counts;
441 auto &part_counts_with_ghost = storage.part_counts_with_ghost;
442 auto &xyz_refs = storage.positions_with_ghosts;
443 auto &pressure_field = storage.pressure;
444 auto &soundspeed_field = storage.soundspeed;
445
446 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> uint_refs
447 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("uint", "u");
448 {
449 shambase::get_check_ref(uint_refs).set_refs(
450 mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
451 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
452 return std::ref(mpdat.get_field<Tscal>(iuint_interf));
453 }));
454 }
455
456 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> vxyz_refs
457 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("vxyz", "v");
458 {
459 shambase::get_check_ref(vxyz_refs).set_refs(
460 mpdats.map<std::reference_wrapper<PatchDataField<Tvec>>>(
461 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
462 return std::ref(mpdat.get_field<Tvec>(ivxyz_interf));
463 }));
464 }
465
466 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hpart_refs
467 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("hpart", "h");
468 { // if was just reset before this call
469 shambase::get_check_ref(hpart_refs)
470 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
471 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
472 return std::ref(mpdat.get_field<Tscal>(ihpart_interf));
473 }));
474 }
475
476 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> omega_refs
477 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("omega", "omega");
478 {
479 shambase::get_check_ref(omega_refs)
480 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
481 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
482 return std::ref(mpdat.get_field<Tscal>(iomega_interf));
483 }));
484 }
485
486 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> alpha_av_refs
487 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("alpha_av", "alpha_av");
488 {
489 shambase::DistributedData<std::reference_wrapper<PatchDataField<Tscal>>> refs{};
490 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
491 refs.add_obj(
492 cur_p.id_patch, std::ref(storage.alpha_av_ghost.get().get(cur_p.id_patch)));
493 });
494 shambase::get_check_ref(alpha_av_refs).set_refs(refs);
495 }
496
497 shamrock::solvergraph::SolverGraph &solver_graph = storage.solver_graph;
498
499 auto axyz_refs = solver_graph.get_edge_ptr<shamrock::solvergraph::FieldRefs<Tvec>>("axyz");
500 auto duint_refs = solver_graph.get_edge_ptr<shamrock::solvergraph::FieldRefs<Tscal>>("duint");
501 auto gpart_mass
502 = solver_graph.get_edge_ptr<shamrock::solvergraph::ScalarEdge<Tscal>>("gpart_mass");
503
504 std::shared_ptr<shamrock::solvergraph::ScalarEdge<Tscal>> alpha_u
505 = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>("alpha_u", "alpha_u");
506 {
507 shambase::get_check_ref(alpha_u).value = cfg.alpha_u;
508 }
509 std::shared_ptr<shamrock::solvergraph::ScalarEdge<Tscal>> beta_AV
510 = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>("beta_AV", "beta_AV");
511 {
512 shambase::get_check_ref(beta_AV).value = cfg.beta_AV;
513 }
514
515 if (solver_config.dust_config.should_use_dust_av()) {
516 u32 ndust = solver_config.dust_config.get_dust_nvar();
517 u32 is_j_interf = ghost_layout.get_field_idx<Tscal>("s_j");
518
519 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> s_j_refs
520 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("s_j", "s_j");
521 {
522 shambase::get_check_ref(s_j_refs).set_refs(
523 mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
524 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
525 return std::ref(mpdat.get_field<Tscal>(is_j_interf));
526 }));
527 }
528
529 std::shared_ptr<NodeUpdateDerivsVaryingAlphaAVDustTVA<Tvec, SPHKernel>> node
530 = std::make_shared<NodeUpdateDerivsVaryingAlphaAVDustTVA<Tvec, SPHKernel>>(ndust);
531 {
532 node->set_edges(
533 gpart_mass,
534 alpha_u,
535 beta_AV,
536 part_counts,
537 part_counts_with_ghost,
538 xyz_refs,
539 hpart_refs,
540 vxyz_refs,
541 uint_refs,
542 omega_refs,
543 pressure_field,
544 soundspeed_field,
545 alpha_av_refs,
546 s_j_refs,
547 storage.neigh_cache,
548 axyz_refs,
549 duint_refs);
550 }
551 node->evaluate();
552 } else {
553 std::shared_ptr<NodeUpdateDerivsVaryingAlphaAV<Tvec, SPHKernel>> node
554 = std::make_shared<NodeUpdateDerivsVaryingAlphaAV<Tvec, SPHKernel>>();
555 {
556 node->set_edges(
557 gpart_mass,
558 alpha_u,
559 beta_AV,
560 part_counts,
561 part_counts_with_ghost,
562 xyz_refs,
563 hpart_refs,
564 vxyz_refs,
565 uint_refs,
566 omega_refs,
567 pressure_field,
568 soundspeed_field,
569 alpha_av_refs,
570 storage.neigh_cache,
571 axyz_refs,
572 duint_refs);
573 }
574 node->evaluate();
575 }
576}
577
578template<class Tvec, template<class> class SPHKernel>
579void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_disc_visco(
580 ConstantDisc cfg) {
581 StackEntry stack_loc{};
582
583 using namespace shamrock;
584 using namespace shamrock::patch;
585
586 PatchDataLayerLayout &pdl = scheduler().pdl_old();
587
588 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
589 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
590 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
591 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
592 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
593 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
594
595 shamrock::patch::PatchDataLayerLayout &ghost_layout
596 = shambase::get_check_ref(storage.ghost_layout.get());
597 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
598 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
599 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
600 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
601
602 auto &merged_xyzh = storage.merged_xyzh.get();
603 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
604 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
605
606 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
607 PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch);
608
609 sham::DeviceBuffer<Tvec> &buf_xyz
610 = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref<Tvec>(0);
611 sham::DeviceBuffer<Tvec> &buf_axyz = pdat.get_field_buf_ref<Tvec>(iaxyz);
612 sham::DeviceBuffer<Tscal> &buf_duint = pdat.get_field_buf_ref<Tscal>(iduint);
613 sham::DeviceBuffer<Tvec> &buf_vxyz = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);
614 sham::DeviceBuffer<Tscal> &buf_hpart = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);
615 sham::DeviceBuffer<Tscal> &buf_omega = mpdat.get_field_buf_ref<Tscal>(iomega_interf);
616 sham::DeviceBuffer<Tscal> &buf_uint = mpdat.get_field_buf_ref<Tscal>(iuint_interf);
617 sham::DeviceBuffer<Tscal> &buf_pressure
618 = shambase::get_check_ref(storage.pressure).get_field(cur_p.id_patch).get_buf();
619 sham::DeviceBuffer<Tscal> &buf_cs
620 = shambase::get_check_ref(storage.soundspeed).get_field(cur_p.id_patch).get_buf();
621
622 sycl::range range_npart{pdat.get_obj_cnt()};
623
624 tree::ObjectCache &pcache
625 = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);
626
628
629 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
630 sham::EventList depends_list;
631
632 auto xyz = buf_xyz.get_read_access(depends_list);
633 auto axyz = buf_axyz.get_write_access(depends_list);
634 auto du = buf_duint.get_write_access(depends_list);
635 auto vxyz = buf_vxyz.get_read_access(depends_list);
636 auto hpart = buf_hpart.get_read_access(depends_list);
637 auto omega = buf_omega.get_read_access(depends_list);
638 auto u = buf_uint.get_read_access(depends_list);
639 auto pressure = buf_pressure.get_read_access(depends_list);
640 auto cs = buf_cs.get_read_access(depends_list);
641 auto ploop_ptrs = pcache.get_read_access(depends_list);
642
643 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
644 const Tscal pmass = solver_config.gpart_mass;
645 const Tscal alpha_AV = cfg.alpha_AV;
646 const Tscal alpha_u = cfg.alpha_u;
647 const Tscal beta_AV = cfg.beta_AV;
648
649 shamlog_debug_sycl_ln("deriv kernel", "alpha_AV :", alpha_AV);
650 shamlog_debug_sycl_ln("deriv kernel", "alpha_u :", alpha_u);
651 shamlog_debug_sycl_ln("deriv kernel", "beta_AV :", beta_AV);
652
653 // tree::ObjectIterator particle_looper(tree,cgh);
654
655 // tree::LeafCacheObjectIterator
656 // particle_looper(tree,*xyz_cell_id,leaf_cache,cgh);
657
658 tree::ObjectCacheIterator particle_looper(ploop_ptrs);
659
660 // sycl::accessor hmax_tree{tree_field_hmax, cgh, sycl::read_only};
661
662 // sycl::stream out {4096,1024,cgh};
663
664 constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;
665
666 shambase::parallel_for(cgh, pdat.get_obj_cnt(), "compute force disc", [=](u64 gid) {
667 u32 id_a = (u32) gid;
668
669 using namespace shamrock::sph;
670
671 Tvec sum_axyz = {0, 0, 0};
672 Tscal sum_du_a = 0;
673
674 Tscal h_a = hpart[id_a];
675 Tvec xyz_a = xyz[id_a];
676 Tvec vxyz_a = vxyz[id_a];
677 Tscal P_a = pressure[id_a];
678 Tscal cs_a = cs[id_a];
679 Tscal omega_a = omega[id_a];
680 const Tscal u_a = u[id_a];
681
682 Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd);
683 Tscal rho_a_sq = rho_a * rho_a;
684 Tscal rho_a_inv = 1. / rho_a;
685
686 // f32 P_a = cs * cs * rho_a;
687
688 Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a);
689
690 Tvec force_pressure{0, 0, 0};
691 Tscal tmpdU_pressure = 0;
692
693 particle_looper.for_each_object(id_a, [&](u32 id_b) {
694 // compute only omega_a
695 Tvec dr = xyz_a - xyz[id_b];
696 Tscal rab2 = sycl::dot(dr, dr);
697 Tscal h_b = hpart[id_b];
698
699 if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) {
700 return;
701 }
702
703 Tvec vxyz_b = vxyz[id_b];
704 const Tscal u_b = u[id_b];
705 Tscal P_b = pressure[id_b];
706 Tscal omega_b = omega[id_b];
707 Tscal cs_b = cs[id_b];
708
709 Tscal rab = sycl::sqrt(rab2);
710
711 Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd);
712 const Tscal alpha_a = alpha_AV;
713 const Tscal alpha_b = alpha_AV;
714 Tscal Fab_a = Kernel::dW_3d(rab, h_a);
715 Tscal Fab_b = Kernel::dW_3d(rab, h_b);
716
717 Tvec v_ab = vxyz_a - vxyz_b;
718
719 Tvec r_ab_unit = dr * sham::inv_sat_positive(rab);
720
721 // f32 P_b = cs * cs * rho_b;
722 Tscal v_ab_r_ab = sycl::dot(v_ab, r_ab_unit);
723 Tscal abs_v_ab_r_ab = sycl::fabs(v_ab_r_ab);
724
725 Tscal vsig_a = alpha_a * cs_a + beta_AV * abs_v_ab_r_ab;
726 Tscal vsig_b = alpha_b * cs_b + beta_AV * abs_v_ab_r_ab;
727
728 Tscal vsig_u = shamrock::sph::vsig_u(P_a, P_b, rho_a, rho_b);
729
730 Tscal qa_ab = shamrock::sph::q_av_disc(
731 rho_a, h_a, rab, alpha_a, cs_a, vsig_a, v_ab_r_ab);
732 Tscal qb_ab = shamrock::sph::q_av_disc(
733 rho_b, h_b, rab, alpha_b, cs_b, vsig_b, v_ab_r_ab);
734
735 add_to_derivs_sph_artif_visco_cond(
736 pmass,
737 rho_a_sq,
738 omega_a_rho_a_inv,
739 rho_a_inv,
740 rho_b,
741 omega_a,
742 omega_b,
743 Fab_a,
744 Fab_b,
745 u_a,
746 u_b,
747 P_a,
748 P_b,
749 alpha_u,
750 v_ab,
751 r_ab_unit,
752 vsig_u,
753 qa_ab,
754 qb_ab,
755
756 force_pressure,
757 tmpdU_pressure);
758 });
759
760 axyz[id_a] = force_pressure;
761 du[id_a] = tmpdU_pressure;
762 });
763 });
764
765 buf_xyz.complete_event_state(e);
766 buf_axyz.complete_event_state(e);
767 buf_duint.complete_event_state(e);
768 buf_vxyz.complete_event_state(e);
769 buf_hpart.complete_event_state(e);
770 buf_omega.complete_event_state(e);
771 buf_uint.complete_event_state(e);
772 buf_pressure.complete_event_state(e);
773 buf_cs.complete_event_state(e);
774
775 sham::EventList resulting_events;
776 resulting_events.add_event(e);
777 pcache.complete_event_state(resulting_events);
778 });
779}
780
781template<class Tvec, template<class> class SPHKernel>
782void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(IdealMHD cfg) {
783 StackEntry stack_loc{};
784
785 using namespace shamrock;
786 using namespace shamrock::patch;
787
788 PatchDataLayerLayout &pdl = scheduler().pdl_old();
789
790 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
791 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
792 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
793 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
794 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
795 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
796 const u32 iB_on_rho = pdl.get_field_idx<Tvec>("B/rho");
797 const u32 idB_on_rho = pdl.get_field_idx<Tvec>("dB/rho");
798 const u32 ipsi_on_ch = pdl.get_field_idx<Tscal>("psi/ch");
799 const u32 idpsi_on_ch = pdl.get_field_idx<Tscal>("dpsi/ch");
800 const u32 idrho_dt = pdl.get_field_idx<Tscal>("drho/dt");
801
802 bool do_MHD_debug = solver_config.do_MHD_debug();
803 const u32 imag_pressure = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("mag_pressure") : -1;
804 const u32 imag_tension = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("mag_tension") : -1;
805 const u32 igas_pressure = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("gas_pressure") : -1;
806 const u32 itensile_corr = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("tensile_corr") : -1;
807 const u32 ipsi_propag = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("psi_propag") : -1;
808 const u32 ipsi_diff = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("psi_diff") : -1;
809 const u32 ipsi_cons = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("psi_cons") : -1;
810 const u32 iu_mhd = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("u_mhd") : -1;
811
812 // Tscal mu_0 = 1.;
813 Tscal const mu_0 = solver_config.get_constant_mu_0();
814
815 shamrock::patch::PatchDataLayerLayout &ghost_layout
816 = shambase::get_check_ref(storage.ghost_layout.get());
817 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
818 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
819 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
820 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
821 u32 iB_on_rho_interf = ghost_layout.get_field_idx<Tvec>("B/rho");
822 u32 ipsi_on_ch_interf = ghost_layout.get_field_idx<Tscal>("psi/ch");
823
824 // logger::raw_ln("charged the ghost fields.");
825
826 auto &merged_xyzh = storage.merged_xyzh.get();
827 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
828 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
829
830 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
831 PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch);
832
833 sham::DeviceBuffer<Tvec> &buf_xyz
834 = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref<Tvec>(0);
835 sham::DeviceBuffer<Tvec> &buf_axyz = pdat.get_field_buf_ref<Tvec>(iaxyz);
836 sham::DeviceBuffer<Tscal> &buf_duint = pdat.get_field_buf_ref<Tscal>(iduint);
837 sham::DeviceBuffer<Tvec> &buf_vxyz = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);
838 sham::DeviceBuffer<Tscal> &buf_hpart = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);
839 sham::DeviceBuffer<Tscal> &buf_omega = mpdat.get_field_buf_ref<Tscal>(iomega_interf);
840 sham::DeviceBuffer<Tscal> &buf_uint = mpdat.get_field_buf_ref<Tscal>(iuint_interf);
841 sham::DeviceBuffer<Tscal> &buf_pressure
842 = shambase::get_check_ref(storage.pressure).get_field(cur_p.id_patch).get_buf();
843 sham::DeviceBuffer<Tscal> &buf_cs
844 = shambase::get_check_ref(storage.soundspeed).get_field(cur_p.id_patch).get_buf();
845
846 sham::DeviceBuffer<Tvec> &buf_dB_on_rho = pdat.get_field_buf_ref<Tvec>(idB_on_rho);
847 sham::DeviceBuffer<Tscal> &buf_dpsi_on_ch = pdat.get_field_buf_ref<Tscal>(idpsi_on_ch);
848 sham::DeviceBuffer<Tscal> &buf_drho_dt = pdat.get_field_buf_ref<Tscal>(idrho_dt);
849 // logger::raw_ln("charged dB dpsi");
850
851 sham::DeviceBuffer<Tvec> &buf_B_on_rho = mpdat.get_field_buf_ref<Tvec>(iB_on_rho_interf);
852 sham::DeviceBuffer<Tscal> &buf_psi_on_ch
853 = mpdat.get_field_buf_ref<Tscal>(ipsi_on_ch_interf);
854
855 // logger::raw_ln("charged B psi");
856 // ADD curlBBBBBBBBB
857
858 sycl::range range_npart{pdat.get_obj_cnt()};
859
860 tree::ObjectCache &pcache
861 = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);
862
864
865 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
866 sham::EventList depends_list;
867
868 auto xyz = buf_xyz.get_read_access(depends_list);
869 auto axyz = buf_axyz.get_write_access(depends_list);
870 auto du = buf_duint.get_write_access(depends_list);
871 auto vxyz = buf_vxyz.get_read_access(depends_list);
872 auto hpart = buf_hpart.get_read_access(depends_list);
873 auto omega = buf_omega.get_read_access(depends_list);
874 auto u = buf_uint.get_read_access(depends_list);
875 auto pressure = buf_pressure.get_read_access(depends_list);
876 auto cs = buf_cs.get_read_access(depends_list);
877 auto B_on_rho = buf_B_on_rho.get_read_access(depends_list);
878 auto psi_on_ch = buf_psi_on_ch.get_read_access(depends_list);
879 auto dB_on_rho = buf_dB_on_rho.get_write_access(depends_list);
880 auto dpsi_on_ch = buf_dpsi_on_ch.get_write_access(depends_list);
881 auto drho_dt = buf_drho_dt.get_write_access(depends_list);
882
883 Tvec *mag_pressure
884 = (do_MHD_debug)
885 ? pdat.get_field_buf_ref<Tvec>(imag_pressure).get_write_access(depends_list)
886 : nullptr;
887 Tvec *mag_tension
888 = (do_MHD_debug)
889 ? pdat.get_field_buf_ref<Tvec>(imag_tension).get_write_access(depends_list)
890 : nullptr;
891 Tvec *gas_pressure
892 = (do_MHD_debug)
893 ? pdat.get_field_buf_ref<Tvec>(igas_pressure).get_write_access(depends_list)
894 : nullptr;
895 Tvec *tensile_corr
896 = (do_MHD_debug)
897 ? pdat.get_field_buf_ref<Tvec>(itensile_corr).get_write_access(depends_list)
898 : nullptr;
899
900 Tscal *psi_propag
901 = (do_MHD_debug)
902 ? pdat.get_field_buf_ref<Tscal>(ipsi_propag).get_write_access(depends_list)
903 : nullptr;
904 Tscal *psi_diff
905 = (do_MHD_debug)
906 ? pdat.get_field_buf_ref<Tscal>(ipsi_diff).get_write_access(depends_list)
907 : nullptr;
908 Tscal *psi_cons
909 = (do_MHD_debug)
910 ? pdat.get_field_buf_ref<Tscal>(ipsi_cons).get_write_access(depends_list)
911 : nullptr;
912
913 Tscal *u_mhd = (do_MHD_debug)
914 ? pdat.get_field_buf_ref<Tscal>(iu_mhd).get_write_access(depends_list)
915 : nullptr;
916
917 auto ploop_ptrs = pcache.get_read_access(depends_list);
918
919 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
920 const Tscal pmass = solver_config.gpart_mass;
921 const Tscal sigma_mhd = cfg.sigma_mhd;
922 const Tscal alpha_u = cfg.alpha_u;
923
924 shamlog_debug_ln("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@", "");
925 shamlog_debug_sycl_ln("deriv kernel", "sigma_mhd :", sigma_mhd);
926 shamlog_debug_sycl_ln("deriv kernel", "alpha_u :", alpha_u);
927 shamlog_debug_ln("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@", "");
928
929 tree::ObjectCacheIterator particle_looper(ploop_ptrs);
930
931 constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;
932
933 shambase::parallel_for(cgh, pdat.get_obj_cnt(), "compute MHD", [=](u64 gid) {
934 u32 id_a = (u32) gid;
935
936 using namespace shamrock::sph;
937
938 Tvec sum_axyz = {0, 0, 0};
939 Tscal sum_du_a = 0;
940
941 Tscal h_a = hpart[id_a];
942 Tvec xyz_a = xyz[id_a];
943 Tvec vxyz_a = vxyz[id_a];
944 Tscal P_a = pressure[id_a];
945 Tscal cs_a = cs[id_a];
946 Tscal omega_a = omega[id_a];
947 const Tscal u_a = u[id_a];
948
949 Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd);
950 Tscal rho_a_sq = rho_a * rho_a;
951 Tscal rho_a_inv = 1. / rho_a;
952
953 Tvec B_a = B_on_rho[id_a] * rho_a;
954 Tscal v_alfven_a = sycl::sqrt(sycl::dot(B_a, B_a) / (mu_0 * rho_a));
955 Tscal v_shock_a = sycl::sqrt(cs_a * cs_a + v_alfven_a * v_alfven_a);
956 Tscal psi_a = psi_on_ch[id_a] * v_shock_a;
957
958 Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a);
959
960 Tvec force_pressure{0, 0, 0};
961 Tscal tmpdU_pressure = 0;
962 Tvec magnetic_eq{0, 0, 0};
963 Tscal psi_eq = 0;
964 Tscal drho_eq = 0;
965
966 Tvec mag_pressure_term{0, 0, 0};
967 Tvec mag_tension_term{0, 0, 0};
968 Tvec gas_pressure_term{0, 0, 0};
969 Tvec tensile_corr_term{0, 0, 0};
970
971 Tscal psi_propag_term = 0;
972 Tscal psi_diff_term = 0;
973 Tscal psi_cons_term = 0;
974
975 Tscal u_mhd_term = 0;
976
977 particle_looper.for_each_object(id_a, [&](u32 id_b) {
978 // compute only omega_a
979 Tvec dr = xyz_a - xyz[id_b];
980 Tscal rab2 = sycl::dot(dr, dr);
981 Tscal h_b = hpart[id_b];
982
983 if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) {
984 return;
985 }
986
987 Tvec vxyz_b = vxyz[id_b];
988 const Tscal u_b = u[id_b];
989 Tscal P_b = pressure[id_b];
990 Tscal omega_b = omega[id_b];
991 Tscal cs_b = cs[id_b];
992
993 Tscal rab = sycl::sqrt(rab2);
994
995 Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd);
996 Tvec B_b = B_on_rho[id_b] * rho_b;
997 Tscal v_alfven_b = sycl::sqrt(sycl::dot(B_b, B_b) / (mu_0 * rho_b));
998 Tscal v_shock_b = sycl::sqrt(cs_b * cs_b + v_alfven_b * v_alfven_b);
999 Tscal psi_b = psi_on_ch[id_b] * v_shock_b;
1000 // const Tscal alpha_a = alpha_AV;
1001 // const Tscal alpha_b = alpha_AV;
1002 Tscal Fab_a = Kernel::dW_3d(rab, h_a);
1003 Tscal Fab_b = Kernel::dW_3d(rab, h_b);
1004
1005 // Tscal sigma_mhd = 0.3;
1006 shamrock::sph::mhd::add_to_derivs_spmhd<Kernel, Tvec, Tscal>(
1007 pmass,
1008 dr,
1009 rab,
1010 rho_a,
1011 rho_a_sq,
1012 omega_a_rho_a_inv,
1013 rho_a_inv,
1014 rho_b,
1015 omega_a,
1016 omega_b,
1017 Fab_a,
1018 Fab_b,
1019 vxyz_a,
1020 vxyz_b,
1021 u_a,
1022 u_b,
1023 P_a,
1024 P_b,
1025 cs_a,
1026 cs_b,
1027 h_a,
1028 h_b,
1029
1030 alpha_u,
1031
1032 B_a,
1033 B_b,
1034
1035 psi_a,
1036 psi_b,
1037
1038 mu_0,
1039 sigma_mhd,
1040
1041 force_pressure,
1042 tmpdU_pressure,
1043 magnetic_eq,
1044 psi_eq,
1045 drho_eq,
1046 mag_pressure_term,
1047 mag_tension_term,
1048 gas_pressure_term,
1049 tensile_corr_term,
1050
1051 psi_propag_term,
1052 psi_diff_term,
1053 psi_cons_term,
1054 u_mhd_term);
1055 });
1056
1057 axyz[id_a] = force_pressure;
1058 du[id_a] = tmpdU_pressure;
1059 dB_on_rho[id_a] = magnetic_eq;
1060 dpsi_on_ch[id_a] = psi_eq - psi_a / h_a;
1061 drho_dt[id_a] = drho_eq;
1062
1063 if (do_MHD_debug) {
1064 mag_pressure[id_a] = mag_pressure_term;
1065 mag_tension[id_a] = mag_tension_term;
1066 gas_pressure[id_a] = gas_pressure_term;
1067 tensile_corr[id_a] = tensile_corr_term;
1068
1069 psi_propag[id_a] = psi_propag_term;
1070 psi_diff[id_a] = psi_diff_term;
1071 psi_cons[id_a] = -psi_a / h_a;
1072
1073 u_mhd[id_a] = u_mhd_term;
1074 }
1075 });
1076 });
1077
1078 buf_xyz.complete_event_state(e);
1079 buf_axyz.complete_event_state(e);
1080 buf_duint.complete_event_state(e);
1081 buf_vxyz.complete_event_state(e);
1082 buf_hpart.complete_event_state(e);
1083 buf_omega.complete_event_state(e);
1084 buf_uint.complete_event_state(e);
1085 buf_pressure.complete_event_state(e);
1086 buf_cs.complete_event_state(e);
1087 buf_B_on_rho.complete_event_state(e);
1088 buf_psi_on_ch.complete_event_state(e);
1089 buf_dB_on_rho.complete_event_state(e);
1090 buf_dpsi_on_ch.complete_event_state(e);
1091 buf_drho_dt.complete_event_state(e);
1092
1093 if (do_MHD_debug) {
1094 pdat.get_field_buf_ref<Tvec>(imag_pressure).complete_event_state(e);
1095 pdat.get_field_buf_ref<Tvec>(imag_tension).complete_event_state(e);
1096 pdat.get_field_buf_ref<Tvec>(igas_pressure).complete_event_state(e);
1097 pdat.get_field_buf_ref<Tvec>(itensile_corr).complete_event_state(e);
1098
1099 pdat.get_field_buf_ref<Tscal>(ipsi_propag).complete_event_state(e);
1100 pdat.get_field_buf_ref<Tscal>(ipsi_diff).complete_event_state(e);
1101 pdat.get_field_buf_ref<Tscal>(ipsi_cons).complete_event_state(e);
1102
1103 pdat.get_field_buf_ref<Tscal>(iu_mhd).complete_event_state(e);
1104 }
1105
1106 sham::EventList resulting_events;
1107 resulting_events.add_event(e);
1108 pcache.complete_event_state(resulting_events);
1109 });
1110}
1111
1112template<class Tvec, template<class> class SPHKernel>
1113void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_dust_monofluid_tva_Sj(
1114 DustConfig cfg, Tscal dt_hydro) {
1115
1116 using MonofluidTVA = typename DustConfig::MonofluidTVA;
1117
1118 StackEntry stack_loc{};
1119
1120 using namespace shamrock;
1121 using namespace shamrock::patch;
1122
1123 PatchDataLayerLayout &pdl = scheduler().pdl_old();
1124
1125 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
1126 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
1127 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
1128 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
1129 const u32 is_j = pdl.get_field_idx<Tscal>("s_j");
1130 const u32 ids_j_dt = pdl.get_field_idx<Tscal>("ds_j_dt");
1131
1132 shamrock::patch::PatchDataLayerLayout &ghost_layout
1133 = shambase::get_check_ref(storage.ghost_layout.get());
1134 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
1135 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
1136 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
1137 u32 is_j_interf = ghost_layout.get_field_idx<Tscal>("s_j");
1138
1139 u32 ndust = cfg.get_dust_nvar();
1140
1141 auto &merged_xyzh = storage.merged_xyzh.get();
1142 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
1143 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
1144
1145 auto &part_counts = storage.part_counts;
1146 auto &part_counts_with_ghost = storage.part_counts_with_ghost;
1147 auto &xyz_refs = storage.positions_with_ghosts;
1148 auto &pressure_field = storage.pressure;
1149
1150 shamrock::solvergraph::SolverGraph &solver_graph = storage.solver_graph;
1151 auto gpart_mass
1152 = solver_graph.get_edge_ptr<shamrock::solvergraph::ScalarEdge<Tscal>>("gpart_mass");
1153
1154 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> vxyz_refs
1155 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("vxyz", "v");
1156 {
1157 shambase::get_check_ref(vxyz_refs).set_refs(
1158 mpdats.map<std::reference_wrapper<PatchDataField<Tvec>>>(
1159 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
1160 return std::ref(mpdat.get_field<Tvec>(ivxyz_interf));
1161 }));
1162 }
1163
1164 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hpart_refs
1165 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("hpart", "h");
1166 { // if was just reset before this call
1167 shambase::get_check_ref(hpart_refs)
1168 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
1169 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
1170 return std::ref(mpdat.get_field<Tscal>(ihpart_interf));
1171 }));
1172 }
1173
1174 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> omega_refs
1175 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("omega", "omega");
1176 {
1177 shambase::get_check_ref(omega_refs)
1178 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
1179 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
1180 return std::ref(mpdat.get_field<Tscal>(iomega_interf));
1181 }));
1182 }
1183
1184 // s_j_interf
1185 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> s_j_refs
1186 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("s_j", "s_j");
1187 {
1188 shambase::get_check_ref(s_j_refs).set_refs(
1189 mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
1190 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
1191 return std::ref(mpdat.get_field<Tscal>(is_j_interf));
1192 }));
1193 }
1194
1195 auto t_j_field
1196 = storage.solver_graph.template get_edge_ptr<shamrock::solvergraph::Field<Tscal>>("Ts_j");
1197
1198 std::shared_ptr<shamrock::solvergraph::Field<Tscal>> Ttilde_sj_field
1199 = std::make_shared<shamrock::solvergraph::Field<Tscal>>(ndust, "Ttilde_sj", "Ttilde_sj");
1200
1201 auto ds_j_dt_refs
1202 = solver_graph.get_edge_ptr<shamrock::solvergraph::FieldRefs<Tscal>>("ds_j_dt");
1203
1204 std::shared_ptr<ComputeDustTtilde<Tvec, SPHKernel>> node_tj
1205 = std::make_shared<ComputeDustTtilde<Tvec, SPHKernel>>(ndust);
1206 {
1207 node_tj->set_edges(
1208 gpart_mass, part_counts_with_ghost, hpart_refs, s_j_refs, t_j_field, Ttilde_sj_field);
1209 }
1210 node_tj->evaluate();
1211
1212 std::shared_ptr<NodeUpdateDerivsMonofluidTVA<Tvec, SPHKernel>> node
1213 = std::make_shared<NodeUpdateDerivsMonofluidTVA<Tvec, SPHKernel>>(ndust);
1214 {
1215 node->set_edges(
1216 gpart_mass,
1217 part_counts,
1218 part_counts_with_ghost,
1219 xyz_refs,
1220 hpart_refs,
1221 vxyz_refs,
1222 omega_refs,
1223 pressure_field,
1224 s_j_refs,
1225 Ttilde_sj_field,
1226 storage.neigh_cache,
1227 ds_j_dt_refs);
1228 }
1229 node->evaluate();
1230
1231 MonofluidTVA &cfg_monofluid_tva
1232 = shambase::get_check_ref((std::get_if<MonofluidTVA>(&cfg.current_mode)));
1233
1234 if (cfg_monofluid_tva.smooth_s_positivity_limiter) {
1235 std::shared_ptr<NodeMonofluidTVASmoothSPositivityLimiter<Tvec>> node_limiter
1236 = std::make_shared<NodeMonofluidTVASmoothSPositivityLimiter<Tvec>>(ndust);
1237 {
1238 node_limiter->set_edges(part_counts, s_j_refs, Ttilde_sj_field, ds_j_dt_refs);
1239 }
1240 node_limiter->evaluate();
1241 }
1242
1243 if (cfg_monofluid_tva.pure_diffusion_mode) {
1244 // reset accelerations & du/dt to 0
1245
1246 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
1247 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
1248
1249 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
1250 pdat.get_field_buf_ref<Tvec>(iaxyz).fill({0, 0, 0});
1251 pdat.get_field_buf_ref<Tscal>(iduint).fill(0);
1252 });
1253 }
1254}
1255
1256using namespace shammath;
1257template class shammodels::sph::modules::UpdateDerivs<f64_3, M4>;
1258template class shammodels::sph::modules::UpdateDerivs<f64_3, M6>;
1259template class shammodels::sph::modules::UpdateDerivs<f64_3, M8>;
1260
1261template class shammodels::sph::modules::UpdateDerivs<f64_3, C2>;
1262template class shammodels::sph::modules::UpdateDerivs<f64_3, C4>;
1263template class shammodels::sph::modules::UpdateDerivs<f64_3, C6>;
ComputeDustTtilde.hpp
Compute the dust combined stopping times Ttilde_sj for each dust species j see Hutchison 2018 eq 15.
shammodels::gsph::names::newtonian::axyz
constexpr const char * axyz
3-acceleration field
Definition FieldNames.hpp:46
shammodels::gsph::names::newtonian::vxyz
constexpr const char * vxyz
3-velocity field
Definition FieldNames.hpp:43
shammodels::gsph::edges::part_counts_with_ghost
constexpr const char * part_counts_with_ghost
Particle counts including ghosts.
Definition FieldNames.hpp:104
shammodels::gsph::names::common::xyz
constexpr const char * xyz
Position field (3D coordinates).
Definition FieldNames.hpp:33
shammodels::gsph::edges::part_counts
constexpr const char * part_counts
Particle counts per patch.
Definition FieldNames.hpp:101
shammodels::gsph::names::newtonian::pressure
constexpr const char * pressure
Pressure P (derived from EOS).
Definition FieldNames.hpp:58
shammodels::gsph::names::common::hpart
constexpr const char * hpart
Smoothing length field.
Definition FieldNames.hpp:34
shammodels::gsph::names::newtonian::omega
constexpr const char * omega
Grad-h correction factor \Omega.
Definition FieldNames.hpp:64
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
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_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::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::DeviceQueue
A SYCL queue associated with a device and a context.
Definition DeviceQueue.hpp:32
sham::DeviceQueue::submit
sycl::event submit(Fct &&fct)
Submits a kernel to the SYCL queue.
Definition DeviceQueue.hpp:101
sham::EventList
Class to manage a list of SYCL events.
Definition EventList.hpp:31
sham::EventList::add_event
void add_event(sycl::event e)
Add an event to the list of events.
Definition EventList.hpp:87
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:251
shambase::DistributedData::get
T & get(u64 id)
Returns a reference to an object in the collection.
Definition DistributedData.hpp:145
shamrock::patch::PatchDataLayerLayout::get_field_idx
u32 get_field_idx(const std::string &field_name) const
Get the field id if matching name & type.
Definition PatchDataLayerLayout.hpp:349
shamrock::patch::PatchDataLayer
PatchDataLayer container class, the layout is described in patchdata_layout.
Definition PatchDataLayer.hpp:35
shamrock::solvergraph::SolverGraph
A graph container for managing solver nodes and edges with type-safe access.
Definition SolverGraph.hpp:110
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
kernel_call.hpp
kernel_call_distrib.hpp
sham::inv_sat_positive
T inv_sat_positive(T v, T minvsat=T{1e-9}, T satval=T{0.}) noexcept
inverse saturated (positive numbers only)
Definition math.hpp:841
shambase::get_check_ref
T & get_check_ref(const std::unique_ptr< T > &ptr, SourceLocation loc=SourceLocation())
Takes a std::unique_ptr and returns a reference to the object it holds. It throws a std::runtime_erro...
Definition memory.hpp:110
shambase::throw_unimplemented
void throw_unimplemented(SourceLocation loc=SourceLocation{})
Throw a std::runtime_error saying that the function is unimplemented.
Definition exception.hpp:153
shammath
namespace for math utility
Definition AABB.hpp:26
shamrock
namespace for the main framework
Definition __init__.py:1
q_ab.hpp
shamrock::sph::q_av
constexpr Tscal q_av(const Tscal &rho, const Tscal &vsig, const Tscal &v_scal_rhat)
phantom_2018 eq.40
Definition q_ab.hpp:37
mhd.hpp
file containing formulas for sphmhd forces, evolution of magnetic and divergence cleaning fields.
forces.hpp
file containing formulas for sph forces
sphkernels.hpp
sph kernels
StackEntry
shambase::details::BasicStackEntry StackEntry
Alias for shambase::details::BasicStackEntry.
Definition stacktrace.hpp:242
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
Generated on for Shamrock by doxygen 1.15.0