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
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/math/density.hpp"
26#include "shammodels/sph/math/forces.hpp"
27#include "shammodels/sph/math/mhd.hpp"
28#include "shammodels/sph/math/q_ab.hpp"
29#include "shammodels/sph/modules/NodeUpdateDerivsVaryingAlphaAV.hpp"
30#include "shammodels/sph/modules/UpdateDerivs.hpp"
31#include "shamphys/mhd.hpp"
32#include "shamrock/patch/PatchDataFieldSpan.hpp"
33#include "shamrock/solvergraph/FieldRefs.hpp"
34#include "shamrock/solvergraph/IFieldSpan.hpp"
35#include "shamrock/solvergraph/Indexes.hpp"
36#include "shamrock/solvergraph/ScalarEdge.hpp"
37
38template<class Tvec, template<class> class SPHKernel>
39void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs() {
40
41 Cfg_AV cfg_av = solver_config.artif_viscosity;
42 Cfg_MHD cfg_mhd = solver_config.mhd_config;
43
44 if (Constant *v = std::get_if<Constant>(&cfg_av.config)) {
45 update_derivs_constantAV(*v);
46 } else if (VaryingMM97 *v = std::get_if<VaryingMM97>(&cfg_av.config)) {
47 update_derivs_mm97(*v);
48 } else if (VaryingCD10 *v = std::get_if<VaryingCD10>(&cfg_av.config)) {
49 update_derivs_cd10(*v);
50 } else if (ConstantDisc *v = std::get_if<ConstantDisc>(&cfg_av.config)) {
51 update_derivs_disc_visco(*v);
52 } else if (IdealMHD *v = std::get_if<IdealMHD>(&cfg_mhd.config)) {
53 update_derivs_MHD(*v);
54 } else if (NonIdealMHD *v = std::get_if<NonIdealMHD>(&cfg_mhd.config)) {
55 shambase::throw_unimplemented();
56 } else if (NoneMHD *v = std::get_if<NoneMHD>(&cfg_mhd.config)) {
57 shambase::throw_unimplemented();
58 } else if (None *v = std::get_if<None>(&cfg_av.config)) {
59 shambase::throw_unimplemented();
60 } else {
61 shambase::throw_unimplemented();
62 }
63}
64
65template<class Tvec, template<class> class SPHKernel>
66void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_noAV(None cfg) {}
67
68template<class Tvec, template<class> class SPHKernel>
69void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_constantAV(
70 Constant cfg) {
71 StackEntry stack_loc{};
72
73 using namespace shamrock;
74 using namespace shamrock::patch;
75
76 PatchDataLayerLayout &pdl = scheduler().pdl_old();
77
78 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
79 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
80 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
81 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
82 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
83 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
84
85 shamrock::patch::PatchDataLayerLayout &ghost_layout
86 = shambase::get_check_ref(storage.ghost_layout.get());
87 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
88 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
89 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
90 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
91
92 auto &merged_xyzh = storage.merged_xyzh.get();
93 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
94 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
95
96 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
97 PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch);
98
99 sham::DeviceBuffer<Tvec> &buf_xyz
100 = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref<Tvec>(0);
101 sham::DeviceBuffer<Tvec> &buf_axyz = pdat.get_field_buf_ref<Tvec>(iaxyz);
102 sham::DeviceBuffer<Tscal> &buf_duint = pdat.get_field_buf_ref<Tscal>(iduint);
103 sham::DeviceBuffer<Tvec> &buf_vxyz = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);
104 sham::DeviceBuffer<Tscal> &buf_hpart = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);
105 sham::DeviceBuffer<Tscal> &buf_omega = mpdat.get_field_buf_ref<Tscal>(iomega_interf);
106 sham::DeviceBuffer<Tscal> &buf_uint = mpdat.get_field_buf_ref<Tscal>(iuint_interf);
107 sham::DeviceBuffer<Tscal> &buf_pressure
108 = shambase::get_check_ref(storage.pressure).get_field(cur_p.id_patch).get_buf();
109 sham::DeviceBuffer<Tscal> &buf_cs
110 = shambase::get_check_ref(storage.soundspeed).get_field(cur_p.id_patch).get_buf();
111
112 sycl::range range_npart{pdat.get_obj_cnt()};
113
114 tree::ObjectCache &pcache
115 = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);
116
118
119 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
120 sham::EventList depends_list;
121
122 auto xyz = buf_xyz.get_read_access(depends_list);
123 auto axyz = buf_axyz.get_write_access(depends_list);
124 auto du = buf_duint.get_write_access(depends_list);
125 auto vxyz = buf_vxyz.get_read_access(depends_list);
126 auto hpart = buf_hpart.get_read_access(depends_list);
127 auto omega = buf_omega.get_read_access(depends_list);
128 auto u = buf_uint.get_read_access(depends_list); // TODO rename to uint
129 auto pressure = buf_pressure.get_read_access(depends_list);
130 auto cs = buf_cs.get_read_access(depends_list);
131 auto ploop_ptrs = pcache.get_read_access(depends_list);
132
133 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
134 const Tscal pmass = solver_config.gpart_mass;
135 const Tscal alpha_u = cfg.alpha_u;
136 const Tscal alpha_AV = cfg.alpha_AV;
137 const Tscal beta_AV = cfg.beta_AV;
138
139 shamlog_debug_sycl_ln("deriv kernel", "alpha_u :", alpha_u);
140 shamlog_debug_sycl_ln("deriv kernel", "alpha_AV :", alpha_AV);
141 shamlog_debug_sycl_ln("deriv kernel", "beta_AV :", beta_AV);
142
143 // tree::ObjectIterator particle_looper(tree,cgh);
144
145 // tree::LeafCacheObjectIterator
146 // particle_looper(tree,*xyz_cell_id,leaf_cache,cgh);
147
148 tree::ObjectCacheIterator particle_looper(ploop_ptrs);
149
150 // sycl::accessor hmax_tree{tree_field_hmax, cgh, sycl::read_only};
151
152 // sycl::stream out {4096,1024,cgh};
153
154 constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;
155
156 shambase::parallel_for(cgh, pdat.get_obj_cnt(), "compute force cte AV", [=](u64 gid) {
157 u32 id_a = (u32) gid;
158
159 using namespace shamrock::sph;
160
161 Tvec sum_axyz = {0, 0, 0};
162 Tscal sum_du_a = 0;
163
164 Tscal h_a = hpart[id_a];
165 Tvec xyz_a = xyz[id_a];
166 Tvec vxyz_a = vxyz[id_a];
167 Tscal P_a = pressure[id_a];
168 Tscal omega_a = omega[id_a];
169 const Tscal u_a = u[id_a];
170
171 Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd);
172 Tscal rho_a_sq = rho_a * rho_a;
173 Tscal rho_a_inv = 1. / rho_a;
174
175 // f32 P_a = cs * cs * rho_a;
176
177 Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a);
178
179 Tscal cs_a = cs[id_a];
180
181 Tvec force_pressure{0, 0, 0};
182 Tscal tmpdU_pressure = 0;
183
184 particle_looper.for_each_object(id_a, [&](u32 id_b) {
185 // compute only omega_a
186 Tvec dr = xyz_a - xyz[id_b];
187 Tscal rab2 = sycl::dot(dr, dr);
188 Tscal h_b = hpart[id_b];
189
190 if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) {
191 return;
192 }
193
194 Tscal rab = sycl::sqrt(rab2);
195 Tvec vxyz_b = vxyz[id_b];
196 const Tscal u_b = u[id_b];
197
198 Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd);
199 Tscal P_b = pressure[id_b];
200 // f32 P_b = cs * cs * rho_b;
201 Tscal omega_b = omega[id_b];
202 Tscal cs_b = cs[id_b];
203
204 const Tscal alpha_a = alpha_AV;
205 const Tscal alpha_b = alpha_AV;
206
207 Tscal Fab_a = Kernel::dW_3d(rab, h_a);
208 Tscal Fab_b = Kernel::dW_3d(rab, h_b);
209
210 Tvec v_ab = vxyz_a - vxyz_b;
211
212 Tvec r_ab_unit = dr * sham::inv_sat_positive(rab);
213
214 // f32 P_b = cs * cs * rho_b;
215 Tscal v_ab_r_ab = sycl::dot(v_ab, r_ab_unit);
216 Tscal abs_v_ab_r_ab = sycl::fabs(v_ab_r_ab);
217
218 Tscal vsig_a = alpha_a * cs_a + beta_AV * abs_v_ab_r_ab;
219 Tscal vsig_b = alpha_b * cs_b + beta_AV * abs_v_ab_r_ab;
220
221 Tscal vsig_u = shamrock::sph::vsig_u(P_a, P_b, rho_a, rho_b);
222
223 Tscal qa_ab = shamrock::sph::q_av(rho_a, vsig_a, v_ab_r_ab);
224 Tscal qb_ab = shamrock::sph::q_av(rho_b, vsig_b, v_ab_r_ab);
225
226 add_to_derivs_sph_artif_visco_cond(
227 pmass,
228 rho_a_sq,
229 omega_a_rho_a_inv,
230 rho_a_inv,
231 rho_b,
232 omega_a,
233 omega_b,
234 Fab_a,
235 Fab_b,
236 u_a,
237 u_b,
238 P_a,
239 P_b,
240 alpha_u,
241 v_ab,
242 r_ab_unit,
243 vsig_u,
244 qa_ab,
245 qb_ab,
246 force_pressure,
247 tmpdU_pressure);
248 });
249 axyz[id_a] = force_pressure;
250 du[id_a] = tmpdU_pressure;
251 });
252 });
253
254 buf_xyz.complete_event_state(e);
255 buf_axyz.complete_event_state(e);
256 buf_duint.complete_event_state(e);
257 buf_vxyz.complete_event_state(e);
258 buf_hpart.complete_event_state(e);
259 buf_omega.complete_event_state(e);
260 buf_uint.complete_event_state(e);
261 buf_pressure.complete_event_state(e);
262 buf_cs.complete_event_state(e);
263
264 sham::EventList resulting_events;
265 resulting_events.add_event(e);
266 pcache.complete_event_state(resulting_events);
267 });
268}
269template<class Tvec, template<class> class SPHKernel>
270void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_mm97(VaryingMM97 cfg) {
271 StackEntry stack_loc{};
272
273 using namespace shamrock;
274 using namespace shamrock::patch;
275
276 PatchDataLayerLayout &pdl = scheduler().pdl_old();
277
278 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
279 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
280 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
281 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
282 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
283 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
284
285 shamrock::patch::PatchDataLayerLayout &ghost_layout
286 = shambase::get_check_ref(storage.ghost_layout.get());
287 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
288 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
289 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
290 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
291
292 auto &merged_xyzh = storage.merged_xyzh.get();
293 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
294 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
295
296 auto &part_counts = storage.part_counts;
297 auto &part_counts_with_ghost = storage.part_counts_with_ghost;
298 auto &xyz_refs = storage.positions_with_ghosts;
299 auto &pressure_field = storage.pressure;
300 auto &soundspeed_field = storage.soundspeed;
301
302 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> uint_refs
303 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("uint", "u");
304 {
305 shambase::get_check_ref(uint_refs).set_refs(
306 mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
307 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
308 return std::ref(mpdat.get_field<Tscal>(iuint_interf));
309 }));
310 }
311
312 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> vxyz_refs
313 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("vxyz", "v");
314 {
315 shambase::get_check_ref(vxyz_refs).set_refs(
316 mpdats.map<std::reference_wrapper<PatchDataField<Tvec>>>(
317 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
318 return std::ref(mpdat.get_field<Tvec>(ivxyz_interf));
319 }));
320 }
321
322 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hpart_refs
323 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("hpart", "h");
324 { // if was just reset before this call
325 shambase::get_check_ref(hpart_refs)
326 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
327 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
328 return std::ref(mpdat.get_field<Tscal>(ihpart_interf));
329 }));
330 }
331
332 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> omega_refs
333 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("omega", "omega");
334 {
335 shambase::get_check_ref(omega_refs)
336 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
337 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
338 return std::ref(mpdat.get_field<Tscal>(iomega_interf));
339 }));
340 }
341
342 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> alpha_av_refs
343 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("alpha_av", "alpha_av");
344 {
345 shambase::DistributedData<std::reference_wrapper<PatchDataField<Tscal>>> refs{};
346 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
347 refs.add_obj(
348 cur_p.id_patch, std::ref(storage.alpha_av_ghost.get().get(cur_p.id_patch)));
349 });
350 shambase::get_check_ref(alpha_av_refs).set_refs(refs);
351 }
352
353 shamrock::solvergraph::SolverGraph &solver_graph = storage.solver_graph;
354
355 auto axyz_refs = solver_graph.get_edge_ptr<shamrock::solvergraph::FieldRefs<Tvec>>("axyz");
356 auto duint_refs = solver_graph.get_edge_ptr<shamrock::solvergraph::FieldRefs<Tscal>>("duint");
357 auto gpart_mass
358 = solver_graph.get_edge_ptr<shamrock::solvergraph::ScalarEdge<Tscal>>("gpart_mass");
359
360 std::shared_ptr<shamrock::solvergraph::ScalarEdge<Tscal>> alpha_u
361 = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>("alpha_u", "alpha_u");
362 {
363 shambase::get_check_ref(alpha_u).value = cfg.alpha_u;
364 }
365 std::shared_ptr<shamrock::solvergraph::ScalarEdge<Tscal>> beta_AV
366 = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>("beta_AV", "beta_AV");
367 {
368 shambase::get_check_ref(beta_AV).value = cfg.beta_AV;
369 }
370
371 std::shared_ptr<NodeUpdateDerivsVaryingAlphaAV<Tvec, SPHKernel>> node
372 = std::make_shared<NodeUpdateDerivsVaryingAlphaAV<Tvec, SPHKernel>>();
373 {
374 node->set_edges(
375 gpart_mass,
376 alpha_u,
377 beta_AV,
378 part_counts,
379 part_counts_with_ghost,
380 xyz_refs,
381 hpart_refs,
382 vxyz_refs,
383 uint_refs,
384 omega_refs,
385 pressure_field,
386 soundspeed_field,
387 alpha_av_refs,
388 storage.neigh_cache,
389 axyz_refs,
390 duint_refs);
391 }
392 node->evaluate();
393}
394template<class Tvec, template<class> class SPHKernel>
395void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_cd10(VaryingCD10 cfg) {
396 StackEntry stack_loc{};
397
398 using namespace shamrock;
399 using namespace shamrock::patch;
400
401 PatchDataLayerLayout &pdl = scheduler().pdl_old();
402
403 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
404 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
405 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
406 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
407 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
408 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
409
410 shamrock::patch::PatchDataLayerLayout &ghost_layout
411 = shambase::get_check_ref(storage.ghost_layout.get());
412 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
413 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
414 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
415 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
416
417 auto &merged_xyzh = storage.merged_xyzh.get();
418 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
419 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
420
421 auto &part_counts = storage.part_counts;
422 auto &part_counts_with_ghost = storage.part_counts_with_ghost;
423 auto &xyz_refs = storage.positions_with_ghosts;
424 auto &pressure_field = storage.pressure;
425 auto &soundspeed_field = storage.soundspeed;
426
427 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> uint_refs
428 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("uint", "u");
429 {
430 shambase::get_check_ref(uint_refs).set_refs(
431 mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
432 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
433 return std::ref(mpdat.get_field<Tscal>(iuint_interf));
434 }));
435 }
436
437 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tvec>> vxyz_refs
438 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("vxyz", "v");
439 {
440 shambase::get_check_ref(vxyz_refs).set_refs(
441 mpdats.map<std::reference_wrapper<PatchDataField<Tvec>>>(
442 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
443 return std::ref(mpdat.get_field<Tvec>(ivxyz_interf));
444 }));
445 }
446
447 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> hpart_refs
448 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("hpart", "h");
449 { // if was just reset before this call
450 shambase::get_check_ref(hpart_refs)
451 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
452 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
453 return std::ref(mpdat.get_field<Tscal>(ihpart_interf));
454 }));
455 }
456
457 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> omega_refs
458 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("omega", "omega");
459 {
460 shambase::get_check_ref(omega_refs)
461 .set_refs(mpdats.map<std::reference_wrapper<PatchDataField<Tscal>>>(
462 [&](u64 id, shamrock::patch::PatchDataLayer &mpdat) {
463 return std::ref(mpdat.get_field<Tscal>(iomega_interf));
464 }));
465 }
466
467 std::shared_ptr<shamrock::solvergraph::FieldRefs<Tscal>> alpha_av_refs
468 = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("alpha_av", "alpha_av");
469 {
470 shambase::DistributedData<std::reference_wrapper<PatchDataField<Tscal>>> refs{};
471 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
472 refs.add_obj(
473 cur_p.id_patch, std::ref(storage.alpha_av_ghost.get().get(cur_p.id_patch)));
474 });
475 shambase::get_check_ref(alpha_av_refs).set_refs(refs);
476 }
477
478 shamrock::solvergraph::SolverGraph &solver_graph = storage.solver_graph;
479
480 auto axyz_refs = solver_graph.get_edge_ptr<shamrock::solvergraph::FieldRefs<Tvec>>("axyz");
481 auto duint_refs = solver_graph.get_edge_ptr<shamrock::solvergraph::FieldRefs<Tscal>>("duint");
482 auto gpart_mass
483 = solver_graph.get_edge_ptr<shamrock::solvergraph::ScalarEdge<Tscal>>("gpart_mass");
484
485 std::shared_ptr<shamrock::solvergraph::ScalarEdge<Tscal>> alpha_u
486 = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>("alpha_u", "alpha_u");
487 {
488 shambase::get_check_ref(alpha_u).value = cfg.alpha_u;
489 }
490 std::shared_ptr<shamrock::solvergraph::ScalarEdge<Tscal>> beta_AV
491 = std::make_shared<shamrock::solvergraph::ScalarEdge<Tscal>>("beta_AV", "beta_AV");
492 {
493 shambase::get_check_ref(beta_AV).value = cfg.beta_AV;
494 }
495
496 std::shared_ptr<NodeUpdateDerivsVaryingAlphaAV<Tvec, SPHKernel>> node
497 = std::make_shared<NodeUpdateDerivsVaryingAlphaAV<Tvec, SPHKernel>>();
498 {
499 node->set_edges(
500 gpart_mass,
501 alpha_u,
502 beta_AV,
503 part_counts,
504 part_counts_with_ghost,
505 xyz_refs,
506 hpart_refs,
507 vxyz_refs,
508 uint_refs,
509 omega_refs,
510 pressure_field,
511 soundspeed_field,
512 alpha_av_refs,
513 storage.neigh_cache,
514 axyz_refs,
515 duint_refs);
516 }
517 node->evaluate();
518}
519
520template<class Tvec, template<class> class SPHKernel>
521void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_disc_visco(
522 ConstantDisc cfg) {
523 StackEntry stack_loc{};
524
525 using namespace shamrock;
526 using namespace shamrock::patch;
527
528 PatchDataLayerLayout &pdl = scheduler().pdl_old();
529
530 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
531 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
532 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
533 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
534 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
535 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
536
537 shamrock::patch::PatchDataLayerLayout &ghost_layout
538 = shambase::get_check_ref(storage.ghost_layout.get());
539 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
540 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
541 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
542 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
543
544 auto &merged_xyzh = storage.merged_xyzh.get();
545 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
546 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
547
548 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
549 PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch);
550
551 sham::DeviceBuffer<Tvec> &buf_xyz
552 = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref<Tvec>(0);
553 sham::DeviceBuffer<Tvec> &buf_axyz = pdat.get_field_buf_ref<Tvec>(iaxyz);
554 sham::DeviceBuffer<Tscal> &buf_duint = pdat.get_field_buf_ref<Tscal>(iduint);
555 sham::DeviceBuffer<Tvec> &buf_vxyz = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);
556 sham::DeviceBuffer<Tscal> &buf_hpart = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);
557 sham::DeviceBuffer<Tscal> &buf_omega = mpdat.get_field_buf_ref<Tscal>(iomega_interf);
558 sham::DeviceBuffer<Tscal> &buf_uint = mpdat.get_field_buf_ref<Tscal>(iuint_interf);
559 sham::DeviceBuffer<Tscal> &buf_pressure
560 = shambase::get_check_ref(storage.pressure).get_field(cur_p.id_patch).get_buf();
561 sham::DeviceBuffer<Tscal> &buf_cs
562 = shambase::get_check_ref(storage.soundspeed).get_field(cur_p.id_patch).get_buf();
563
564 sycl::range range_npart{pdat.get_obj_cnt()};
565
566 tree::ObjectCache &pcache
567 = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);
568
570
571 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
572 sham::EventList depends_list;
573
574 auto xyz = buf_xyz.get_read_access(depends_list);
575 auto axyz = buf_axyz.get_write_access(depends_list);
576 auto du = buf_duint.get_write_access(depends_list);
577 auto vxyz = buf_vxyz.get_read_access(depends_list);
578 auto hpart = buf_hpart.get_read_access(depends_list);
579 auto omega = buf_omega.get_read_access(depends_list);
580 auto u = buf_uint.get_read_access(depends_list);
581 auto pressure = buf_pressure.get_read_access(depends_list);
582 auto cs = buf_cs.get_read_access(depends_list);
583 auto ploop_ptrs = pcache.get_read_access(depends_list);
584
585 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
586 const Tscal pmass = solver_config.gpart_mass;
587 const Tscal alpha_AV = cfg.alpha_AV;
588 const Tscal alpha_u = cfg.alpha_u;
589 const Tscal beta_AV = cfg.beta_AV;
590
591 shamlog_debug_sycl_ln("deriv kernel", "alpha_AV :", alpha_AV);
592 shamlog_debug_sycl_ln("deriv kernel", "alpha_u :", alpha_u);
593 shamlog_debug_sycl_ln("deriv kernel", "beta_AV :", beta_AV);
594
595 // tree::ObjectIterator particle_looper(tree,cgh);
596
597 // tree::LeafCacheObjectIterator
598 // particle_looper(tree,*xyz_cell_id,leaf_cache,cgh);
599
600 tree::ObjectCacheIterator particle_looper(ploop_ptrs);
601
602 // sycl::accessor hmax_tree{tree_field_hmax, cgh, sycl::read_only};
603
604 // sycl::stream out {4096,1024,cgh};
605
606 constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;
607
608 shambase::parallel_for(cgh, pdat.get_obj_cnt(), "compute force disc", [=](u64 gid) {
609 u32 id_a = (u32) gid;
610
611 using namespace shamrock::sph;
612
613 Tvec sum_axyz = {0, 0, 0};
614 Tscal sum_du_a = 0;
615
616 Tscal h_a = hpart[id_a];
617 Tvec xyz_a = xyz[id_a];
618 Tvec vxyz_a = vxyz[id_a];
619 Tscal P_a = pressure[id_a];
620 Tscal cs_a = cs[id_a];
621 Tscal omega_a = omega[id_a];
622 const Tscal u_a = u[id_a];
623
624 Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd);
625 Tscal rho_a_sq = rho_a * rho_a;
626 Tscal rho_a_inv = 1. / rho_a;
627
628 // f32 P_a = cs * cs * rho_a;
629
630 Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a);
631
632 Tvec force_pressure{0, 0, 0};
633 Tscal tmpdU_pressure = 0;
634
635 particle_looper.for_each_object(id_a, [&](u32 id_b) {
636 // compute only omega_a
637 Tvec dr = xyz_a - xyz[id_b];
638 Tscal rab2 = sycl::dot(dr, dr);
639 Tscal h_b = hpart[id_b];
640
641 if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) {
642 return;
643 }
644
645 Tvec vxyz_b = vxyz[id_b];
646 const Tscal u_b = u[id_b];
647 Tscal P_b = pressure[id_b];
648 Tscal omega_b = omega[id_b];
649 Tscal cs_b = cs[id_b];
650
651 Tscal rab = sycl::sqrt(rab2);
652
653 Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd);
654 const Tscal alpha_a = alpha_AV;
655 const Tscal alpha_b = alpha_AV;
656 Tscal Fab_a = Kernel::dW_3d(rab, h_a);
657 Tscal Fab_b = Kernel::dW_3d(rab, h_b);
658
659 Tvec v_ab = vxyz_a - vxyz_b;
660
661 Tvec r_ab_unit = dr * sham::inv_sat_positive(rab);
662
663 // f32 P_b = cs * cs * rho_b;
664 Tscal v_ab_r_ab = sycl::dot(v_ab, r_ab_unit);
665 Tscal abs_v_ab_r_ab = sycl::fabs(v_ab_r_ab);
666
667 Tscal vsig_a = alpha_a * cs_a + beta_AV * abs_v_ab_r_ab;
668 Tscal vsig_b = alpha_b * cs_b + beta_AV * abs_v_ab_r_ab;
669
670 Tscal vsig_u = shamrock::sph::vsig_u(P_a, P_b, rho_a, rho_b);
671
672 Tscal qa_ab = shamrock::sph::q_av_disc(
673 rho_a, h_a, rab, alpha_a, cs_a, vsig_a, v_ab_r_ab);
674 Tscal qb_ab = shamrock::sph::q_av_disc(
675 rho_b, h_b, rab, alpha_b, cs_b, vsig_b, v_ab_r_ab);
676
677 add_to_derivs_sph_artif_visco_cond(
678 pmass,
679 rho_a_sq,
680 omega_a_rho_a_inv,
681 rho_a_inv,
682 rho_b,
683 omega_a,
684 omega_b,
685 Fab_a,
686 Fab_b,
687 u_a,
688 u_b,
689 P_a,
690 P_b,
691 alpha_u,
692 v_ab,
693 r_ab_unit,
694 vsig_u,
695 qa_ab,
696 qb_ab,
697
698 force_pressure,
699 tmpdU_pressure);
700 });
701
702 axyz[id_a] = force_pressure;
703 du[id_a] = tmpdU_pressure;
704 });
705 });
706
707 buf_xyz.complete_event_state(e);
708 buf_axyz.complete_event_state(e);
709 buf_duint.complete_event_state(e);
710 buf_vxyz.complete_event_state(e);
711 buf_hpart.complete_event_state(e);
712 buf_omega.complete_event_state(e);
713 buf_uint.complete_event_state(e);
714 buf_pressure.complete_event_state(e);
715 buf_cs.complete_event_state(e);
716
717 sham::EventList resulting_events;
718 resulting_events.add_event(e);
719 pcache.complete_event_state(resulting_events);
720 });
721}
722
723template<class Tvec, template<class> class SPHKernel>
724void shammodels::sph::modules::UpdateDerivs<Tvec, SPHKernel>::update_derivs_MHD(IdealMHD cfg) {
725 StackEntry stack_loc{};
726
727 using namespace shamrock;
728 using namespace shamrock::patch;
729
730 PatchDataLayerLayout &pdl = scheduler().pdl_old();
731
732 const u32 ixyz = pdl.get_field_idx<Tvec>("xyz");
733 const u32 ivxyz = pdl.get_field_idx<Tvec>("vxyz");
734 const u32 iaxyz = pdl.get_field_idx<Tvec>("axyz");
735 const u32 iuint = pdl.get_field_idx<Tscal>("uint");
736 const u32 iduint = pdl.get_field_idx<Tscal>("duint");
737 const u32 ihpart = pdl.get_field_idx<Tscal>("hpart");
738 const u32 iB_on_rho = pdl.get_field_idx<Tvec>("B/rho");
739 const u32 idB_on_rho = pdl.get_field_idx<Tvec>("dB/rho");
740 const u32 ipsi_on_ch = pdl.get_field_idx<Tscal>("psi/ch");
741 const u32 idpsi_on_ch = pdl.get_field_idx<Tscal>("dpsi/ch");
742 const u32 idrho_dt = pdl.get_field_idx<Tscal>("drho/dt");
743
744 bool do_MHD_debug = solver_config.do_MHD_debug();
745 const u32 imag_pressure = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("mag_pressure") : -1;
746 const u32 imag_tension = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("mag_tension") : -1;
747 const u32 igas_pressure = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("gas_pressure") : -1;
748 const u32 itensile_corr = (do_MHD_debug) ? pdl.get_field_idx<Tvec>("tensile_corr") : -1;
749 const u32 ipsi_propag = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("psi_propag") : -1;
750 const u32 ipsi_diff = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("psi_diff") : -1;
751 const u32 ipsi_cons = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("psi_cons") : -1;
752 const u32 iu_mhd = (do_MHD_debug) ? pdl.get_field_idx<Tscal>("u_mhd") : -1;
753
754 // Tscal mu_0 = 1.;
755 Tscal const mu_0 = solver_config.get_constant_mu_0();
756
757 shamrock::patch::PatchDataLayerLayout &ghost_layout
758 = shambase::get_check_ref(storage.ghost_layout.get());
759 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
760 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
761 u32 ivxyz_interf = ghost_layout.get_field_idx<Tvec>("vxyz");
762 u32 iomega_interf = ghost_layout.get_field_idx<Tscal>("omega");
763 u32 iB_on_rho_interf = ghost_layout.get_field_idx<Tvec>("B/rho");
764 u32 ipsi_on_ch_interf = ghost_layout.get_field_idx<Tscal>("psi/ch");
765
766 // logger::raw_ln("charged the ghost fields.");
767
768 auto &merged_xyzh = storage.merged_xyzh.get();
769 shamrock::solvergraph::Field<Tscal> &omega = shambase::get_check_ref(storage.omega);
770 shambase::DistributedData<PatchDataLayer> &mpdats = storage.merged_patchdata_ghost.get();
771
772 scheduler().for_each_patchdata_nonempty([&](Patch cur_p, PatchDataLayer &pdat) {
773 PatchDataLayer &mpdat = mpdats.get(cur_p.id_patch);
774
775 sham::DeviceBuffer<Tvec> &buf_xyz
776 = merged_xyzh.get(cur_p.id_patch).template get_field_buf_ref<Tvec>(0);
777 sham::DeviceBuffer<Tvec> &buf_axyz = pdat.get_field_buf_ref<Tvec>(iaxyz);
778 sham::DeviceBuffer<Tscal> &buf_duint = pdat.get_field_buf_ref<Tscal>(iduint);
779 sham::DeviceBuffer<Tvec> &buf_vxyz = mpdat.get_field_buf_ref<Tvec>(ivxyz_interf);
780 sham::DeviceBuffer<Tscal> &buf_hpart = mpdat.get_field_buf_ref<Tscal>(ihpart_interf);
781 sham::DeviceBuffer<Tscal> &buf_omega = mpdat.get_field_buf_ref<Tscal>(iomega_interf);
782 sham::DeviceBuffer<Tscal> &buf_uint = mpdat.get_field_buf_ref<Tscal>(iuint_interf);
783 sham::DeviceBuffer<Tscal> &buf_pressure
784 = shambase::get_check_ref(storage.pressure).get_field(cur_p.id_patch).get_buf();
785 sham::DeviceBuffer<Tscal> &buf_cs
786 = shambase::get_check_ref(storage.soundspeed).get_field(cur_p.id_patch).get_buf();
787
788 sham::DeviceBuffer<Tvec> &buf_dB_on_rho = pdat.get_field_buf_ref<Tvec>(idB_on_rho);
789 sham::DeviceBuffer<Tscal> &buf_dpsi_on_ch = pdat.get_field_buf_ref<Tscal>(idpsi_on_ch);
790 sham::DeviceBuffer<Tscal> &buf_drho_dt = pdat.get_field_buf_ref<Tscal>(idrho_dt);
791 // logger::raw_ln("charged dB dpsi");
792
793 sham::DeviceBuffer<Tvec> &buf_B_on_rho = mpdat.get_field_buf_ref<Tvec>(iB_on_rho_interf);
794 sham::DeviceBuffer<Tscal> &buf_psi_on_ch
795 = mpdat.get_field_buf_ref<Tscal>(ipsi_on_ch_interf);
796
797 // logger::raw_ln("charged B psi");
798 // ADD curlBBBBBBBBB
799
800 sycl::range range_npart{pdat.get_obj_cnt()};
801
802 tree::ObjectCache &pcache
803 = shambase::get_check_ref(storage.neigh_cache).get_cache(cur_p.id_patch);
804
806
807 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
808 sham::EventList depends_list;
809
810 auto xyz = buf_xyz.get_read_access(depends_list);
811 auto axyz = buf_axyz.get_write_access(depends_list);
812 auto du = buf_duint.get_write_access(depends_list);
813 auto vxyz = buf_vxyz.get_read_access(depends_list);
814 auto hpart = buf_hpart.get_read_access(depends_list);
815 auto omega = buf_omega.get_read_access(depends_list);
816 auto u = buf_uint.get_read_access(depends_list);
817 auto pressure = buf_pressure.get_read_access(depends_list);
818 auto cs = buf_cs.get_read_access(depends_list);
819 auto B_on_rho = buf_B_on_rho.get_read_access(depends_list);
820 auto psi_on_ch = buf_psi_on_ch.get_read_access(depends_list);
821 auto dB_on_rho = buf_dB_on_rho.get_write_access(depends_list);
822 auto dpsi_on_ch = buf_dpsi_on_ch.get_write_access(depends_list);
823 auto drho_dt = buf_drho_dt.get_write_access(depends_list);
824
825 Tvec *mag_pressure
826 = (do_MHD_debug)
827 ? pdat.get_field_buf_ref<Tvec>(imag_pressure).get_write_access(depends_list)
828 : nullptr;
829 Tvec *mag_tension
830 = (do_MHD_debug)
831 ? pdat.get_field_buf_ref<Tvec>(imag_tension).get_write_access(depends_list)
832 : nullptr;
833 Tvec *gas_pressure
834 = (do_MHD_debug)
835 ? pdat.get_field_buf_ref<Tvec>(igas_pressure).get_write_access(depends_list)
836 : nullptr;
837 Tvec *tensile_corr
838 = (do_MHD_debug)
839 ? pdat.get_field_buf_ref<Tvec>(itensile_corr).get_write_access(depends_list)
840 : nullptr;
841
842 Tscal *psi_propag
843 = (do_MHD_debug)
844 ? pdat.get_field_buf_ref<Tscal>(ipsi_propag).get_write_access(depends_list)
845 : nullptr;
846 Tscal *psi_diff
847 = (do_MHD_debug)
848 ? pdat.get_field_buf_ref<Tscal>(ipsi_diff).get_write_access(depends_list)
849 : nullptr;
850 Tscal *psi_cons
851 = (do_MHD_debug)
852 ? pdat.get_field_buf_ref<Tscal>(ipsi_cons).get_write_access(depends_list)
853 : nullptr;
854
855 Tscal *u_mhd = (do_MHD_debug)
856 ? pdat.get_field_buf_ref<Tscal>(iu_mhd).get_write_access(depends_list)
857 : nullptr;
858
859 auto ploop_ptrs = pcache.get_read_access(depends_list);
860
861 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
862 const Tscal pmass = solver_config.gpart_mass;
863 const Tscal sigma_mhd = cfg.sigma_mhd;
864 const Tscal alpha_u = cfg.alpha_u;
865
866 shamlog_debug_ln("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@", "");
867 shamlog_debug_sycl_ln("deriv kernel", "sigma_mhd :", sigma_mhd);
868 shamlog_debug_sycl_ln("deriv kernel", "alpha_u :", alpha_u);
869 shamlog_debug_ln("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@", "");
870
871 tree::ObjectCacheIterator particle_looper(ploop_ptrs);
872
873 constexpr Tscal Rker2 = Kernel::Rkern * Kernel::Rkern;
874
875 shambase::parallel_for(cgh, pdat.get_obj_cnt(), "compute MHD", [=](u64 gid) {
876 u32 id_a = (u32) gid;
877
878 using namespace shamrock::sph;
879
880 Tvec sum_axyz = {0, 0, 0};
881 Tscal sum_du_a = 0;
882
883 Tscal h_a = hpart[id_a];
884 Tvec xyz_a = xyz[id_a];
885 Tvec vxyz_a = vxyz[id_a];
886 Tscal P_a = pressure[id_a];
887 Tscal cs_a = cs[id_a];
888 Tscal omega_a = omega[id_a];
889 const Tscal u_a = u[id_a];
890
891 Tscal rho_a = rho_h(pmass, h_a, Kernel::hfactd);
892 Tscal rho_a_sq = rho_a * rho_a;
893 Tscal rho_a_inv = 1. / rho_a;
894
895 Tvec B_a = B_on_rho[id_a] * rho_a;
896 Tscal v_alfven_a = sycl::sqrt(sycl::dot(B_a, B_a) / (mu_0 * rho_a));
897 Tscal v_shock_a = sycl::sqrt(cs_a * cs_a + v_alfven_a * v_alfven_a);
898 Tscal psi_a = psi_on_ch[id_a] * v_shock_a;
899
900 Tscal omega_a_rho_a_inv = 1 / (omega_a * rho_a);
901
902 Tvec force_pressure{0, 0, 0};
903 Tscal tmpdU_pressure = 0;
904 Tvec magnetic_eq{0, 0, 0};
905 Tscal psi_eq = 0;
906 Tscal drho_eq = 0;
907
908 Tvec mag_pressure_term{0, 0, 0};
909 Tvec mag_tension_term{0, 0, 0};
910 Tvec gas_pressure_term{0, 0, 0};
911 Tvec tensile_corr_term{0, 0, 0};
912
913 Tscal psi_propag_term = 0;
914 Tscal psi_diff_term = 0;
915 Tscal psi_cons_term = 0;
916
917 Tscal u_mhd_term = 0;
918
919 particle_looper.for_each_object(id_a, [&](u32 id_b) {
920 // compute only omega_a
921 Tvec dr = xyz_a - xyz[id_b];
922 Tscal rab2 = sycl::dot(dr, dr);
923 Tscal h_b = hpart[id_b];
924
925 if (rab2 > h_a * h_a * Rker2 && rab2 > h_b * h_b * Rker2) {
926 return;
927 }
928
929 Tvec vxyz_b = vxyz[id_b];
930 const Tscal u_b = u[id_b];
931 Tscal P_b = pressure[id_b];
932 Tscal omega_b = omega[id_b];
933 Tscal cs_b = cs[id_b];
934
935 Tscal rab = sycl::sqrt(rab2);
936
937 Tscal rho_b = rho_h(pmass, h_b, Kernel::hfactd);
938 Tvec B_b = B_on_rho[id_b] * rho_b;
939 Tscal v_alfven_b = sycl::sqrt(sycl::dot(B_b, B_b) / (mu_0 * rho_b));
940 Tscal v_shock_b = sycl::sqrt(cs_b * cs_b + v_alfven_b * v_alfven_b);
941 Tscal psi_b = psi_on_ch[id_b] * v_shock_b;
942 // const Tscal alpha_a = alpha_AV;
943 // const Tscal alpha_b = alpha_AV;
944 Tscal Fab_a = Kernel::dW_3d(rab, h_a);
945 Tscal Fab_b = Kernel::dW_3d(rab, h_b);
946
947 // Tscal sigma_mhd = 0.3;
948 shamrock::sph::mhd::add_to_derivs_spmhd<Kernel, Tvec, Tscal>(
949 pmass,
950 dr,
951 rab,
952 rho_a,
953 rho_a_sq,
954 omega_a_rho_a_inv,
955 rho_a_inv,
956 rho_b,
957 omega_a,
958 omega_b,
959 Fab_a,
960 Fab_b,
961 vxyz_a,
962 vxyz_b,
963 u_a,
964 u_b,
965 P_a,
966 P_b,
967 cs_a,
968 cs_b,
969 h_a,
970 h_b,
971
972 alpha_u,
973
974 B_a,
975 B_b,
976
977 psi_a,
978 psi_b,
979
980 mu_0,
981 sigma_mhd,
982
983 force_pressure,
984 tmpdU_pressure,
985 magnetic_eq,
986 psi_eq,
987 drho_eq,
988 mag_pressure_term,
989 mag_tension_term,
990 gas_pressure_term,
991 tensile_corr_term,
992
993 psi_propag_term,
994 psi_diff_term,
995 psi_cons_term,
996 u_mhd_term);
997 });
998
999 axyz[id_a] = force_pressure;
1000 du[id_a] = tmpdU_pressure;
1001 dB_on_rho[id_a] = magnetic_eq;
1002 dpsi_on_ch[id_a] = psi_eq - psi_a / h_a;
1003 drho_dt[id_a] = drho_eq;
1004
1005 if (do_MHD_debug) {
1006 mag_pressure[id_a] = mag_pressure_term;
1007 mag_tension[id_a] = mag_tension_term;
1008 gas_pressure[id_a] = gas_pressure_term;
1009 tensile_corr[id_a] = tensile_corr_term;
1010
1011 psi_propag[id_a] = psi_propag_term;
1012 psi_diff[id_a] = psi_diff_term;
1013 psi_cons[id_a] = -psi_a / h_a;
1014
1015 u_mhd[id_a] = u_mhd_term;
1016 }
1017 });
1018 });
1019
1020 buf_xyz.complete_event_state(e);
1021 buf_axyz.complete_event_state(e);
1022 buf_duint.complete_event_state(e);
1023 buf_vxyz.complete_event_state(e);
1024 buf_hpart.complete_event_state(e);
1025 buf_omega.complete_event_state(e);
1026 buf_uint.complete_event_state(e);
1027 buf_pressure.complete_event_state(e);
1028 buf_cs.complete_event_state(e);
1029 buf_B_on_rho.complete_event_state(e);
1030 buf_psi_on_ch.complete_event_state(e);
1031 buf_dB_on_rho.complete_event_state(e);
1032 buf_dpsi_on_ch.complete_event_state(e);
1033 buf_drho_dt.complete_event_state(e);
1034
1035 if (do_MHD_debug) {
1036 pdat.get_field_buf_ref<Tvec>(imag_pressure).complete_event_state(e);
1037 pdat.get_field_buf_ref<Tvec>(imag_tension).complete_event_state(e);
1038 pdat.get_field_buf_ref<Tvec>(igas_pressure).complete_event_state(e);
1039 pdat.get_field_buf_ref<Tvec>(itensile_corr).complete_event_state(e);
1040
1041 pdat.get_field_buf_ref<Tscal>(ipsi_propag).complete_event_state(e);
1042 pdat.get_field_buf_ref<Tscal>(ipsi_diff).complete_event_state(e);
1043 pdat.get_field_buf_ref<Tscal>(ipsi_cons).complete_event_state(e);
1044
1045 pdat.get_field_buf_ref<Tscal>(iu_mhd).complete_event_state(e);
1046 }
1047
1048 sham::EventList resulting_events;
1049 resulting_events.add_event(e);
1050 pcache.complete_event_state(resulting_events);
1051 });
1052}
1053
1054using namespace shammath;
1055template class shammodels::sph::modules::UpdateDerivs<f64_3, M4>;
1056template class shammodels::sph::modules::UpdateDerivs<f64_3, M6>;
1057template class shammodels::sph::modules::UpdateDerivs<f64_3, M8>;
1058
1059template class shammodels::sph::modules::UpdateDerivs<f64_3, C2>;
1060template class shammodels::sph::modules::UpdateDerivs<f64_3, C4>;
1061template class shammodels::sph::modules::UpdateDerivs<f64_3, C6>;
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::DeviceScheduler::get_queue
DeviceQueue & get_queue(u32 id=0)
Get a reference to a DeviceQueue.
Definition DeviceScheduler.cpp:57
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::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:56
shamrock::solvergraph::SolverGraph::get_edge_ptr
std::shared_ptr< T > get_edge_ptr(const std::string &name)
Get a typed shared pointer to an edge by name.
Definition SolverGraph.hpp:276
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)
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
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 Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8