Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
ComputeEos.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
19#include "shambase/DistributedData.hpp"
20#include "shambase/exception.hpp"
21#include "shambase/memory.hpp"
22#include "shambase/narrowing.hpp"
23#include "shambackends/kernel_call.hpp"
24#include "shambackends/kernel_call_distrib.hpp"
25#include "shammath/sphkernels.hpp"
26#include "shammodels/sph/math/density.hpp"
27#include "shammodels/sph/modules/ComputeEos.hpp"
28#include "shammodels/sph/sink_edges_helper.hpp"
29#include "shamphys/eos.hpp"
30#include "shamrock/patch/PatchDataField.hpp"
31#include "shamrock/patch/PatchDataLayer.hpp"
32#include "shamrock/solvergraph/Field.hpp"
33#include "shamrock/solvergraph/FieldRefs.hpp"
34#include "shamrock/solvergraph/Indexes.hpp"
35#include "shamsolvergraph/edge/IDataEdge.hpp"
36#include "shamsolvergraph/node/INode.hpp"
37#include <memory>
38
39#define NODE_EDGES(X_RO, X_RW, X_RO_OPTIONAL, X_RW_OPTIONAL) \
40 /* ------------------- inputs ------------------- */ \
41 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, cs) \
42 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, hfactd) \
43 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, pmass) \
44 X_RO_OPTIONAL(shamrock::solvergraph::IFieldSpan<Tscal>, spans_rho) \
45 X_RO_OPTIONAL(shamrock::solvergraph::IFieldSpan<Tscal>, spans_h) \
46 X_RO(shamrock::solvergraph::Indexes<u32>, sizes) \
47 \
48 /* ------------------- outputs ------------------- */ \
49 X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, spans_pressure) \
50 X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, spans_soundspeed)
51
52namespace shammodels::common::modules {
53 template<class Tvec>
54 class ComputeEOSIsothermal : public shamrock::solvergraph::INode {
55
56 using Tscal = shambase::VecComponent<Tvec>;
57
58 public:
59 ComputeEOSIsothermal() = default;
60
61 EXPAND_NODE_EDGES_OPTIONAL(NODE_EDGES)
62
63 inline static void internal_eos(
64 const Tscal &cs, const Tscal &rho, Tscal &pressure, Tscal &soundspeed) noexcept {
65 using EOS = shamphys::EOS_Isothermal<Tscal>;
66 Tscal P_a = EOS::pressure(cs, rho);
67 pressure = P_a;
68 soundspeed = cs;
69 }
70
71 inline void _impl_evaluate_internal() {
72
73 __shamrock_stack_entry();
74
75 auto edges = get_edges();
76
77 bool has_rho = edges.spans_rho.has_value();
78 bool has_h = edges.spans_h.has_value();
79
80 // must have either rho or h
81 if ((has_rho && has_h) || (!has_rho && !has_h)) {
82 throw shambase::make_except_with_loc<std::invalid_argument>(
83 "Must have either rho or h");
84 }
85
86 edges.spans_pressure.ensure_sizes(edges.sizes.indexes);
87 edges.spans_soundspeed.ensure_sizes(edges.sizes.indexes);
88
89 Tscal cs = edges.cs.data;
90 Tscal pmass = edges.pmass.data;
91 Tscal hfactd = edges.hfactd.data;
92
93 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
94
95 auto out_refs = sham::DDMultiRef{
96 edges.spans_pressure.get_spans(), edges.spans_soundspeed.get_spans()};
97
98 if (has_rho) {
99 auto &spans_rho = edges.spans_rho.value().get();
100 spans_rho.check_sizes(edges.sizes.indexes);
101
102 sham::distributed_data_kernel_call(
103 dev_sched,
104 sham::DDMultiRef{spans_rho.get_spans()},
105 out_refs,
106 edges.sizes.indexes,
107 [cs](u32 gid, const Tscal *rho, Tscal *pressure, Tscal *soundspeed) {
108 Tscal rho_a = rho[gid];
109 internal_eos(cs, rho_a, pressure[gid], soundspeed[gid]);
110 });
111 } else if (has_h) {
112 auto &spans_h = edges.spans_h.value().get();
113 spans_h.check_sizes(edges.sizes.indexes);
114
115 sham::distributed_data_kernel_call(
116 dev_sched,
117 sham::DDMultiRef{spans_h.get_spans()},
118 out_refs,
119 edges.sizes.indexes,
120 [cs, pmass, hfactd](
121 u32 gid, const Tscal *h, Tscal *pressure, Tscal *soundspeed) {
122 using namespace shamrock::sph;
123 Tscal rho = rho_h(pmass, h[gid], hfactd);
124 internal_eos(cs, rho, pressure[gid], soundspeed[gid]);
125 });
126 }
127 }
128
129 inline virtual std::string _impl_get_label() const { return "ComputeEOSIsothermal"; };
130
131 inline virtual std::string _impl_get_tex() const { return "TODO"; };
132 };
133} // namespace shammodels::common::modules
134
135#undef NODE_EDGES
136
137#define NODE_EDGES(X_RO, X_RW, X_RO_OPTIONAL, X_RW_OPTIONAL) \
138 /* ------------------- inputs ------------------- */ \
139 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, gamma) \
140 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, hfactd) \
141 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, pmass) \
142 X_RO_OPTIONAL(shamrock::solvergraph::IFieldSpan<Tscal>, spans_rho) \
143 X_RO_OPTIONAL(shamrock::solvergraph::IFieldSpan<Tscal>, spans_h) \
144 X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, spans_uint) \
145 X_RO(shamrock::solvergraph::Indexes<u32>, sizes) \
146 \
147 /* ------------------- outputs ------------------- */ \
148 X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, spans_pressure) \
149 X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, spans_soundspeed)
150
151namespace shammodels::common::modules {
152 template<class Tvec>
153 class ComputeEOSAdiabatic : public shamrock::solvergraph::INode {
154
155 using Tscal = shambase::VecComponent<Tvec>;
156
157 public:
158 ComputeEOSAdiabatic() = default;
159
160 EXPAND_NODE_EDGES_OPTIONAL(NODE_EDGES)
161
162 inline static void internal_eos(
163 const Tscal &gamma,
164 const Tscal &rho,
165 const Tscal &uint,
166 Tscal &pressure,
167 Tscal &soundspeed) noexcept {
168 using EOS = shamphys::EOS_Adiabatic<Tscal>;
169 Tscal P_a = EOS::pressure(gamma, rho, uint);
170 Tscal cs_a = EOS::cs_from_p(gamma, rho, P_a);
171 pressure = P_a;
172 soundspeed = cs_a;
173 }
174
175 inline void _impl_evaluate_internal() {
176
177 __shamrock_stack_entry();
178
179 auto edges = get_edges();
180
181 bool has_rho = edges.spans_rho.has_value();
182 bool has_h = edges.spans_h.has_value();
183
184 // must have either rho or h
185 if ((has_rho && has_h) || (!has_rho && !has_h)) {
186 throw shambase::make_except_with_loc<std::invalid_argument>(
187 "Must have either rho or h");
188 }
189
190 edges.spans_pressure.ensure_sizes(edges.sizes.indexes);
191 edges.spans_soundspeed.ensure_sizes(edges.sizes.indexes);
192
193 Tscal gamma = edges.gamma.data;
194 Tscal pmass = edges.pmass.data;
195 Tscal hfactd = edges.hfactd.data;
196
197 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
198
199 auto out_refs = sham::DDMultiRef{
200 edges.spans_pressure.get_spans(), edges.spans_soundspeed.get_spans()};
201
202 if (has_rho) {
203 auto &spans_rho = edges.spans_rho.value().get();
204 spans_rho.check_sizes(edges.sizes.indexes);
205
206 sham::distributed_data_kernel_call(
207 dev_sched,
208 sham::DDMultiRef{spans_rho.get_spans(), edges.spans_uint.get_spans()},
209 out_refs,
210 edges.sizes.indexes,
211 [gamma](
212 u32 gid,
213 const Tscal *rho,
214 const Tscal *uint,
215 Tscal *pressure,
216 Tscal *soundspeed) {
217 Tscal rho_a = rho[gid];
218 Tscal uint_a = uint[gid];
219 internal_eos(gamma, rho_a, uint_a, pressure[gid], soundspeed[gid]);
220 });
221 } else if (has_h) {
222 auto &spans_h = edges.spans_h.value().get();
223 spans_h.check_sizes(edges.sizes.indexes);
224
225 sham::distributed_data_kernel_call(
226 dev_sched,
227 sham::DDMultiRef{spans_h.get_spans(), edges.spans_uint.get_spans()},
228 out_refs,
229 edges.sizes.indexes,
230 [gamma, pmass, hfactd](
231 u32 gid,
232 const Tscal *h,
233 const Tscal *uint,
234 Tscal *pressure,
235 Tscal *soundspeed) {
236 using namespace shamrock::sph;
237 Tscal rho = rho_h(pmass, h[gid], hfactd);
238 Tscal uint_a = uint[gid];
239 internal_eos(gamma, rho, uint_a, pressure[gid], soundspeed[gid]);
240 });
241 }
242 }
243
244 inline virtual std::string _impl_get_label() const { return "ComputeEOSAdiabatic"; };
245
246 inline virtual std::string _impl_get_tex() const { return "TODO"; };
247 };
248} // namespace shammodels::common::modules
249
250#undef NODE_EDGES
251
252#define NODE_EDGES(X_RO, X_RW, X_RO_OPTIONAL, X_RW_OPTIONAL) \
253 /* ------------------- inputs ------------------- */ \
254 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, K) \
255 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, gamma) \
256 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, hfactd) \
257 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, pmass) \
258 X_RO_OPTIONAL(shamrock::solvergraph::IFieldSpan<Tscal>, spans_rho) \
259 X_RO_OPTIONAL(shamrock::solvergraph::IFieldSpan<Tscal>, spans_h) \
260 X_RO(shamrock::solvergraph::Indexes<u32>, sizes) \
261 \
262 /* ------------------- outputs ------------------- */ \
263 X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, spans_pressure) \
264 X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, spans_soundspeed)
265
266namespace shammodels::common::modules {
267 template<class Tvec>
268 class ComputeEOSPolytropic : public shamrock::solvergraph::INode {
269
270 using Tscal = shambase::VecComponent<Tvec>;
271
272 public:
273 ComputeEOSPolytropic() = default;
274
275 EXPAND_NODE_EDGES_OPTIONAL(NODE_EDGES)
276
277 inline static void internal_eos(
278 const Tscal &K,
279 const Tscal &gamma,
280 const Tscal &rho,
281 Tscal &pressure,
282 Tscal &soundspeed) noexcept {
283 using EOS = shamphys::EOS_Polytropic<Tscal>;
284 Tscal P_a = EOS::pressure(gamma, K, rho);
285 Tscal cs_a = EOS::soundspeed(gamma, K, rho);
286 pressure = P_a;
287 soundspeed = cs_a;
288 }
289
290 inline void _impl_evaluate_internal() {
291
292 __shamrock_stack_entry();
293
294 auto edges = get_edges();
295
296 bool has_rho = edges.spans_rho.has_value();
297 bool has_h = edges.spans_h.has_value();
298
299 // must have either rho or h
300 if ((has_rho && has_h) || (!has_rho && !has_h)) {
301 throw shambase::make_except_with_loc<std::invalid_argument>(
302 "Must have either rho or h");
303 }
304
305 edges.spans_pressure.ensure_sizes(edges.sizes.indexes);
306 edges.spans_soundspeed.ensure_sizes(edges.sizes.indexes);
307
308 Tscal K = edges.K.data;
309 Tscal gamma = edges.gamma.data;
310 Tscal pmass = edges.pmass.data;
311 Tscal hfactd = edges.hfactd.data;
312
313 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
314
315 auto out_refs = sham::DDMultiRef{
316 edges.spans_pressure.get_spans(), edges.spans_soundspeed.get_spans()};
317
318 if (has_rho) {
319 auto &spans_rho = edges.spans_rho.value().get();
320 spans_rho.check_sizes(edges.sizes.indexes);
321
322 sham::distributed_data_kernel_call(
323 dev_sched,
324 sham::DDMultiRef{spans_rho.get_spans()},
325 out_refs,
326 edges.sizes.indexes,
327 [K, gamma](u32 gid, const Tscal *rho, Tscal *pressure, Tscal *soundspeed) {
328 Tscal rho_a = rho[gid];
329 internal_eos(K, gamma, rho_a, pressure[gid], soundspeed[gid]);
330 });
331 } else if (has_h) {
332 auto &spans_h = edges.spans_h.value().get();
333 spans_h.check_sizes(edges.sizes.indexes);
334
335 sham::distributed_data_kernel_call(
336 dev_sched,
337 sham::DDMultiRef{spans_h.get_spans()},
338 out_refs,
339 edges.sizes.indexes,
340 [K, gamma, pmass, hfactd](
341 u32 gid, const Tscal *h, Tscal *pressure, Tscal *soundspeed) {
342 using namespace shamrock::sph;
343 Tscal rho = rho_h(pmass, h[gid], hfactd);
344 internal_eos(K, gamma, rho, pressure[gid], soundspeed[gid]);
345 });
346 }
347 }
348
349 inline virtual std::string _impl_get_label() const { return "ComputeEOSPolytropic"; };
350
351 inline virtual std::string _impl_get_tex() const { return "TODO"; };
352 };
353} // namespace shammodels::common::modules
354
355#undef NODE_EDGES
356
357#define NODE_EDGES(X_RO, X_RW, X_RO_OPTIONAL, X_RW_OPTIONAL) \
358 /* ------------------- inputs ------------------- */ \
359 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, hfactd) \
360 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, pmass) \
361 X_RO_OPTIONAL(shamrock::solvergraph::IFieldSpan<Tscal>, spans_rho) \
362 X_RO_OPTIONAL(shamrock::solvergraph::IFieldSpan<Tscal>, spans_h) \
363 X_RO(shamrock::solvergraph::IFieldSpan<Tscal>, spans_cs0) \
364 X_RO(shamrock::solvergraph::Indexes<u32>, sizes) \
365 \
366 /* ------------------- outputs ------------------- */ \
367 X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, spans_pressure) \
368 X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, spans_soundspeed)
369
370namespace shammodels::common::modules {
371 template<class Tvec>
372 class ComputeEOSLocallyIsothermal : public shamrock::solvergraph::INode {
373
374 using Tscal = shambase::VecComponent<Tvec>;
375
376 public:
377 ComputeEOSLocallyIsothermal() = default;
378
379 EXPAND_NODE_EDGES_OPTIONAL(NODE_EDGES)
380
381 inline static void internal_eos(
382 const Tscal &cs0, const Tscal &rho, Tscal &pressure, Tscal &soundspeed) noexcept {
383 using EOS = shamphys::EOS_LocallyIsothermal<Tscal>;
384 pressure = EOS::pressure_from_cs(cs0 * cs0, rho);
385 soundspeed = cs0;
386 }
387
388 inline void _impl_evaluate_internal() {
389
390 __shamrock_stack_entry();
391
392 auto edges = get_edges();
393
394 bool has_rho = edges.spans_rho.has_value();
395 bool has_h = edges.spans_h.has_value();
396
397 // must have either rho or h
398 if ((has_rho && has_h) || (!has_rho && !has_h)) {
399 throw shambase::make_except_with_loc<std::invalid_argument>(
400 "Must have either rho or h");
401 }
402
403 edges.spans_pressure.ensure_sizes(edges.sizes.indexes);
404 edges.spans_soundspeed.ensure_sizes(edges.sizes.indexes);
405
406 Tscal pmass = edges.pmass.data;
407 Tscal hfactd = edges.hfactd.data;
408
409 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
410
411 auto out_refs = sham::DDMultiRef{
412 edges.spans_pressure.get_spans(), edges.spans_soundspeed.get_spans()};
413
414 if (has_rho) {
415 auto &spans_rho = edges.spans_rho.value().get();
416 spans_rho.check_sizes(edges.sizes.indexes);
417
418 sham::distributed_data_kernel_call(
419 dev_sched,
420 sham::DDMultiRef{spans_rho.get_spans(), edges.spans_cs0.get_spans()},
421 out_refs,
422 edges.sizes.indexes,
423 [](u32 gid,
424 const Tscal *rho,
425 const Tscal *cs0,
426 Tscal *pressure,
427 Tscal *soundspeed) {
428 Tscal rho_a = rho[gid];
429 Tscal cs0_a = cs0[gid];
430 internal_eos(cs0_a, rho_a, pressure[gid], soundspeed[gid]);
431 });
432 } else if (has_h) {
433 auto &spans_h = edges.spans_h.value().get();
434 spans_h.check_sizes(edges.sizes.indexes);
435
436 sham::distributed_data_kernel_call(
437 dev_sched,
438 sham::DDMultiRef{spans_h.get_spans(), edges.spans_cs0.get_spans()},
439 out_refs,
440 edges.sizes.indexes,
441 [pmass, hfactd](
442 u32 gid,
443 const Tscal *h,
444 const Tscal *cs0,
445 Tscal *pressure,
446 Tscal *soundspeed) {
447 using namespace shamrock::sph;
448 Tscal rho = rho_h(pmass, h[gid], hfactd);
449 Tscal cs0_a = cs0[gid];
450 internal_eos(cs0_a, rho, pressure[gid], soundspeed[gid]);
451 });
452 }
453 }
454
455 inline virtual std::string _impl_get_label() const {
456 return "ComputeEOSLocallyIsothermal";
457 };
458
459 inline virtual std::string _impl_get_tex() const { return "TODO"; };
460 };
461} // namespace shammodels::common::modules
462
463#undef NODE_EDGES
464
465#define NODE_EDGES(X_RO, X_RW, X_RO_OPTIONAL, X_RW_OPTIONAL) \
466 /* ------------------- inputs ------------------- */ \
467 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, mu_e) \
468 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, density_unit) \
469 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, pressure_unit) \
470 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, velocity_unit) \
471 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, hfactd) \
472 X_RO(shamrock::solvergraph::IDataEdge<Tscal>, pmass) \
473 X_RO_OPTIONAL(shamrock::solvergraph::IFieldSpan<Tscal>, spans_rho) \
474 X_RO_OPTIONAL(shamrock::solvergraph::IFieldSpan<Tscal>, spans_h) \
475 X_RO(shamrock::solvergraph::Indexes<u32>, sizes) \
476 \
477 /* ------------------- outputs ------------------- */ \
478 X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, spans_pressure) \
479 X_RW(shamrock::solvergraph::IFieldSpan<Tscal>, spans_soundspeed)
480
481namespace shammodels::common::modules {
482 template<class Tvec>
483 class ComputeEOSFermi : public shamrock::solvergraph::INode {
484
485 using Tscal = shambase::VecComponent<Tvec>;
486
487 public:
488 ComputeEOSFermi() = default;
489
490 EXPAND_NODE_EDGES_OPTIONAL(NODE_EDGES)
491
492 inline static void internal_eos(
493 const Tscal &mu_e,
494 const Tscal &density_unit,
495 const Tscal &pressure_unit,
496 const Tscal &velocity_unit,
497 const Tscal &rho,
498 Tscal &pressure,
499 Tscal &soundspeed) noexcept {
500 using EOS = shamphys::EOS_Fermi<Tscal>;
501 auto const res = EOS::pressure_and_soundspeed(mu_e, rho * density_unit);
502 pressure = res.pressure / pressure_unit;
503 soundspeed = res.soundspeed / velocity_unit;
504 }
505
506 inline void _impl_evaluate_internal() {
507
508 __shamrock_stack_entry();
509
510 auto edges = get_edges();
511
512 bool has_rho = edges.spans_rho.has_value();
513 bool has_h = edges.spans_h.has_value();
514
515 // must have either rho or h
516 if ((has_rho && has_h) || (!has_rho && !has_h)) {
517 throw shambase::make_except_with_loc<std::invalid_argument>(
518 "Must have either rho or h");
519 }
520
521 edges.spans_pressure.ensure_sizes(edges.sizes.indexes);
522 edges.spans_soundspeed.ensure_sizes(edges.sizes.indexes);
523
524 Tscal mu_e = edges.mu_e.data;
525 Tscal density_unit = edges.density_unit.data;
526 Tscal pressure_unit = edges.pressure_unit.data;
527 Tscal velocity_unit = edges.velocity_unit.data;
528 Tscal pmass = edges.pmass.data;
529 Tscal hfactd = edges.hfactd.data;
530
531 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
532
533 auto out_refs = sham::DDMultiRef{
534 edges.spans_pressure.get_spans(), edges.spans_soundspeed.get_spans()};
535
536 if (has_rho) {
537 auto &spans_rho = edges.spans_rho.value().get();
538 spans_rho.check_sizes(edges.sizes.indexes);
539
540 sham::distributed_data_kernel_call(
541 dev_sched,
542 sham::DDMultiRef{spans_rho.get_spans()},
543 out_refs,
544 edges.sizes.indexes,
545 [mu_e, density_unit, pressure_unit, velocity_unit](
546 u32 gid, const Tscal *rho, Tscal *pressure, Tscal *soundspeed) {
547 Tscal rho_a = rho[gid];
548 internal_eos(
549 mu_e,
550 density_unit,
551 pressure_unit,
552 velocity_unit,
553 rho_a,
554 pressure[gid],
555 soundspeed[gid]);
556 });
557 } else if (has_h) {
558 auto &spans_h = edges.spans_h.value().get();
559 spans_h.check_sizes(edges.sizes.indexes);
560
561 sham::distributed_data_kernel_call(
562 dev_sched,
563 sham::DDMultiRef{spans_h.get_spans()},
564 out_refs,
565 edges.sizes.indexes,
566 [mu_e, density_unit, pressure_unit, velocity_unit, pmass, hfactd](
567 u32 gid, const Tscal *h, Tscal *pressure, Tscal *soundspeed) {
568 using namespace shamrock::sph;
569 Tscal rho = rho_h(pmass, h[gid], hfactd);
570 internal_eos(
571 mu_e,
572 density_unit,
573 pressure_unit,
574 velocity_unit,
575 rho,
576 pressure[gid],
577 soundspeed[gid]);
578 });
579 }
580 }
581
582 inline virtual std::string _impl_get_label() const { return "ComputeEOSFermi"; };
583
584 inline virtual std::string _impl_get_tex() const { return "TODO"; };
585 };
586} // namespace shammodels::common::modules
587
588#undef NODE_EDGES
589
590template<class Tvec, template<class> class SPHKernel>
591void shammodels::sph::modules::ComputeEos<Tvec, SPHKernel>::compute_eos_internal(
592 const std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> &hfactd,
593 const std::shared_ptr<shamrock::solvergraph::IDataEdge<Tscal>> &pmass,
594 const std::optional<std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>>> &spans_rho,
595 const std::optional<std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>>> &spans_h,
596 const std::optional<std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>>> &spans_uint,
597 const std::shared_ptr<shamrock::solvergraph::Indexes<u32>> &sizes,
598 const std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> &spans_pressure,
599 const std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>> &spans_soundspeed) {
600
601 shamrock::patch::PatchDataLayerLayout &ghost_layout
602 = shambase::get_check_ref(storage.ghost_layout.get());
603
604 using namespace shamrock;
605 using namespace shamrock::patch;
606
607 using SolverConfigEOS = typename Config::EOSConfig;
608 using SolverEOS_Isothermal = typename SolverConfigEOS::Isothermal;
609 using SolverEOS_Adiabatic = typename SolverConfigEOS::Adiabatic;
610 using SolverEOS_Polytropic = typename SolverConfigEOS::Polytropic;
611 using SolverEOS_LocallyIsothermal = typename SolverConfigEOS::LocallyIsothermal;
612 using SolverEOS_LocallyIsothermalLP07 = typename SolverConfigEOS::LocallyIsothermalLP07;
613 using SolverEOS_LocallyIsothermalFA2014 = typename SolverConfigEOS::LocallyIsothermalFA2014;
614 using SolverEOS_LocallyIsothermalFA2014Extended =
615 typename SolverConfigEOS::LocallyIsothermalFA2014Extended;
616 using SolverEOS_Fermi = typename SolverConfigEOS::Fermi;
617
618 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
619 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
620
621 const auto &sizes_indexes = shambase::get_check_ref(sizes).indexes;
622
623 shambase::get_check_ref(storage.pressure).ensure_sizes(sizes_indexes);
624 shambase::get_check_ref(storage.soundspeed).ensure_sizes(sizes_indexes);
625
626 bool has_rho = spans_rho.has_value();
627 bool has_h = spans_h.has_value();
628
629 // must have either rho or h
630 if ((!has_rho || has_h) && (has_rho || !has_h)) {
631 throw shambase::make_except_with_loc<std::invalid_argument>("Must have either rho or h");
632 }
633
634 auto map_opt_span
635 = [](const std::optional<std::shared_ptr<shamrock::solvergraph::IFieldSpan<Tscal>>> &opt)
636 -> std::optional<std::reference_wrapper<shamrock::solvergraph::IFieldSpan<Tscal>>> {
637 if (!opt.has_value()) {
638 return std::nullopt;
639 }
640 return std::ref(shambase::get_check_ref(opt.value()));
641 };
642
643 struct Edges {
644 const std::optional<std::reference_wrapper<shamrock::solvergraph::IFieldSpan<Tscal>>>
645 spans_rho;
646 const std::optional<std::reference_wrapper<shamrock::solvergraph::IFieldSpan<Tscal>>>
647 spans_h;
648 const std::optional<std::reference_wrapper<shamrock::solvergraph::IFieldSpan<Tscal>>>
649 spans_uint;
650 } edges{map_opt_span(spans_rho), map_opt_span(spans_h), map_opt_span(spans_uint)};
651
652 auto out_refs = sham::DDMultiRef{
653 shambase::get_check_ref(spans_pressure).get_spans(),
654 shambase::get_check_ref(spans_soundspeed).get_spans()};
655
656 if (SolverEOS_Isothermal *eos_config
657 = std::get_if<SolverEOS_Isothermal>(&solver_config.eos_config.config)) {
658
659 auto cs = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("cs", "c_s");
660 cs->data = eos_config->cs;
661
662 shammodels::common::modules::ComputeEOSIsothermal<Tvec> node;
663 node.set_edges(
664 cs, hfactd, pmass, spans_rho, spans_h, sizes, spans_pressure, spans_soundspeed);
665 node.evaluate();
666 } else if (
667 SolverEOS_Adiabatic *eos_config
668 = std::get_if<SolverEOS_Adiabatic>(&solver_config.eos_config.config)) {
669
670 auto gamma = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("gamma", "\\gamma");
671 gamma->data = eos_config->gamma;
672
673 shammodels::common::modules::ComputeEOSAdiabatic<Tvec> node;
674 node.set_edges(
675 gamma,
676 hfactd,
677 pmass,
678 spans_rho,
679 spans_h,
680 spans_uint.value(),
681 sizes,
682 spans_pressure,
683 spans_soundspeed);
684 node.evaluate();
685 } else if (
686 SolverEOS_Polytropic *eos_config
687 = std::get_if<SolverEOS_Polytropic>(&solver_config.eos_config.config)) {
688
689 auto K = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("K", "K");
690 auto gamma = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("gamma", "\\gamma");
691 K->data = eos_config->K;
692 gamma->data = eos_config->gamma;
693
694 shammodels::common::modules::ComputeEOSPolytropic<Tvec> node;
695 node.set_edges(
696 K, gamma, hfactd, pmass, spans_rho, spans_h, sizes, spans_pressure, spans_soundspeed);
697 node.evaluate();
698 } else if (
699 [[maybe_unused]] SolverEOS_LocallyIsothermal *eos_config
700 = std::get_if<SolverEOS_LocallyIsothermal>(&solver_config.eos_config.config)) {
701
702 u32 isoundspeed_interf = ghost_layout.get_field_idx<Tscal>("soundspeed");
703
704 auto soundspeed_refs
705 = shamrock::solvergraph::FieldRefs<Tscal>::make_shared("cs0", "c_{s,0}");
706 auto refs = storage.merged_patchdata_ghost.get()
707 .template map<shamrock::solvergraph::PatchDataFieldRef<Tscal>>(
708 [&](u64 id, PatchDataLayer &mpdat)
709 -> shamrock::solvergraph::PatchDataFieldRef<Tscal> {
710 return mpdat.get_field<Tscal>(isoundspeed_interf);
711 });
712 soundspeed_refs->set_refs(refs);
713
714 shammodels::common::modules::ComputeEOSLocallyIsothermal<Tvec> node;
715 node.set_edges(
716 hfactd,
717 pmass,
718 spans_rho,
719 spans_h,
720 soundspeed_refs,
721 sizes,
722 spans_pressure,
723 spans_soundspeed);
724 node.evaluate();
725 } else if (
726 SolverEOS_LocallyIsothermalLP07 *eos_config
727 = std::get_if<SolverEOS_LocallyIsothermalLP07>(&solver_config.eos_config.config)) {
728
729 Tscal cs0 = eos_config->cs0;
730 Tscal r0sq = eos_config->r0 * eos_config->r0;
731 Tscal mq = -eos_config->q;
732
733 Tscal pmass_ = shambase::get_check_ref(pmass).data;
734 Tscal hfactd_ = shambase::get_check_ref(hfactd).data;
735
736 shamrock::solvergraph::FieldRefs<Tvec> xyz_refs{"", ""};
737 auto refs
738 = storage.merged_xyzh.get()
739 .template map<shamrock::solvergraph::PatchDataFieldRef<Tvec>>(
740 [&](u64 id,
741 PatchDataLayer &mpdat) -> shamrock::solvergraph::PatchDataFieldRef<Tvec> {
742 return mpdat.get_field<Tvec>(0);
743 });
744 xyz_refs.set_refs(refs);
745
746 using EOS = shamphys::EOS_LocallyIsothermal<Tscal>;
747
748 auto eos_internal = [](Tvec R,
749 Tscal cs0,
750 Tscal r0sq,
751 Tscal mq,
752 Tscal rho_a,
753 Tscal &pressure,
754 Tscal &soundspeed) {
755 Tscal Rsq = sycl::dot(R, R);
756 Tscal cs_sq = EOS::soundspeed_sq(cs0 * cs0, Rsq / r0sq, mq);
757 Tscal cs_out = sycl::sqrt(cs_sq);
758
759 Tscal P_a = EOS::pressure_from_cs(cs_sq, rho_a);
760
761 pressure = P_a;
762 soundspeed = cs_out;
763 };
764
765 if (has_rho) {
766 auto &spans_rho_ = edges.spans_rho.value().get();
767 spans_rho_.check_sizes(sizes_indexes);
768 sham::distributed_data_kernel_call(
769 dev_sched,
770 sham::DDMultiRef{spans_rho_.get_spans(), xyz_refs.get_spans()},
771 out_refs,
772 sizes_indexes,
773 [cs0, r0sq, mq, eos_internal](
774 u32 gid,
775 const Tscal *rho,
776 const Tvec *xyz,
777 Tscal *pressure,
778 Tscal *soundspeed) {
779 Tvec R_a = xyz[gid];
780 Tscal rho_a = rho[gid];
781 eos_internal(R_a, cs0, r0sq, mq, rho_a, pressure[gid], soundspeed[gid]);
782 });
783 } else if (has_h) {
784 auto &spans_h_ = edges.spans_h.value().get();
785 spans_h_.check_sizes(sizes_indexes);
786 sham::distributed_data_kernel_call(
787 dev_sched,
788 sham::DDMultiRef{spans_h_.get_spans(), xyz_refs.get_spans()},
789 out_refs,
790 sizes_indexes,
791 [cs0, r0sq, mq, pmass_, hfactd_, eos_internal](
792 u32 gid, const Tscal *h, const Tvec *xyz, Tscal *pressure, Tscal *soundspeed) {
793 using namespace shamrock::sph;
794 Tvec R_a = xyz[gid];
795 Tscal rho_a = rho_h(pmass_, h[gid], hfactd_);
796 eos_internal(R_a, cs0, r0sq, mq, rho_a, pressure[gid], soundspeed[gid]);
797 });
798 }
799
800 } else if (
801 SolverEOS_LocallyIsothermalFA2014 *eos_config
802 = std::get_if<SolverEOS_LocallyIsothermalFA2014>(&solver_config.eos_config.config)) {
803
804 Tscal G = solver_config.get_constant_G();
805 Tscal h_over_r = eos_config->h_over_r;
806 Tscal pmass_ = shambase::get_check_ref(pmass).data;
807 Tscal hfactd_ = shambase::get_check_ref(hfactd).data;
808
809 using EOS = shamphys::EOS_LocallyIsothermal<Tscal>;
810
811 auto &sink_pos = get_sink_pos<Tvec>(scheduler().synchronized_data);
812 auto &sink_mass = get_sink_mass<Tvec>(scheduler().synchronized_data);
813 u32 sink_cnt = shambase::narrow_or_throw<u32>(sink_pos.size());
814
815 if (sink_cnt == 0) {
816 throw shambase::make_except_with_loc<std::runtime_error>(
817 "No sinks found for the equation of state");
818 }
819
820 shamrock::solvergraph::FieldRefs<Tvec> xyz_refs{"", ""};
821 auto refs
822 = storage.merged_xyzh.get()
823 .template map<shamrock::solvergraph::PatchDataFieldRef<Tvec>>(
824 [&](u64 id,
825 PatchDataLayer &mpdat) -> shamrock::solvergraph::PatchDataFieldRef<Tvec> {
826 return mpdat.get_field<Tvec>(0);
827 });
828 xyz_refs.set_refs(refs);
829
830 sham::DeviceBuffer<Tvec> sink_pos_buf(sink_pos.size(), dev_sched);
831 sham::DeviceBuffer<Tscal> sink_mass_buf(sink_mass.size(), dev_sched);
832
833 sink_pos_buf.copy_from_stdvec(sink_pos);
834 sink_mass_buf.copy_from_stdvec(sink_mass);
835
836 auto eos_internal = [](Tvec R,
837 Tscal rho_a,
838 u32 scount,
839 auto spos,
840 auto smass,
841 Tscal G,
842 Tscal h_over_r,
843 Tscal &pressure,
844 Tscal &soundspeed) {
845 Tscal mpotential = 0;
846 for (u32 i = 0; i < scount; i++) {
847 Tvec s_r = spos[i] - R;
848 Tscal s_m = smass[i];
849 Tscal s_r_abs = sycl::length(s_r);
850 mpotential += G * s_m / s_r_abs;
851 }
852
853 Tscal cs_out = h_over_r * sycl::sqrt(mpotential);
854 Tscal P_a = EOS::pressure_from_cs(cs_out * cs_out, rho_a);
855
856 pressure = P_a;
857 soundspeed = cs_out;
858 };
859
860 if (has_rho) {
861 auto &spans_rho_ = edges.spans_rho.value().get();
862 spans_rho_.check_sizes(sizes_indexes);
863
864 sizes_indexes.for_each([&](u64 id, u32 count) {
865 sham::kernel_call(
866 q,
867 sham::MultiRef{
868 spans_rho_.get_spans().get(id),
869 xyz_refs.get_spans().get(id),
870 sink_pos_buf,
871 sink_mass_buf},
872 out_refs.get(id),
873 count,
874 [G, h_over_r, sink_cnt, eos_internal](
875 u32 gid,
876 const Tscal *rho,
877 const Tvec *xyz,
878 const Tvec *spos,
879 const Tscal *smass,
880 Tscal *pressure,
881 Tscal *soundspeed) {
882 Tvec R_a = xyz[gid];
883 Tscal rho_a = rho[gid];
884 eos_internal(
885 R_a,
886 rho_a,
887 sink_cnt,
888 spos,
889 smass,
890 G,
891 h_over_r,
892 pressure[gid],
893 soundspeed[gid]);
894 });
895 });
896
897 } else if (has_h) {
898 auto &spans_h_ = edges.spans_h.value().get();
899 spans_h_.check_sizes(sizes_indexes);
900
901 sizes_indexes.for_each([&](u64 id, u32 count) {
902 sham::kernel_call(
903 q,
904 sham::MultiRef{
905 spans_h_.get_spans().get(id),
906 xyz_refs.get_spans().get(id),
907 sink_pos_buf,
908 sink_mass_buf},
909 out_refs.get(id),
910 count,
911 [G, h_over_r, sink_cnt, pmass_, hfactd_, eos_internal](
912 u32 gid,
913 const Tscal *h,
914 const Tvec *xyz,
915 const Tvec *spos,
916 const Tscal *smass,
917 Tscal *pressure,
918 Tscal *soundspeed) {
919 using namespace shamrock::sph;
920 Tvec R_a = xyz[gid];
921 Tscal rho_a = rho_h(pmass_, h[gid], hfactd_);
922 eos_internal(
923 R_a,
924 rho_a,
925 sink_cnt,
926 spos,
927 smass,
928 G,
929 h_over_r,
930 pressure[gid],
931 soundspeed[gid]);
932 });
933 });
934 }
935
936 } else if (
937 SolverEOS_LocallyIsothermalFA2014Extended *eos_config
938 = std::get_if<SolverEOS_LocallyIsothermalFA2014Extended>(
939 &solver_config.eos_config.config)) {
940
941 Tscal cs0 = eos_config->cs0;
942 Tscal r0 = eos_config->r0;
943 Tscal q_ = eos_config->q;
944 Tscal pmass_ = shambase::get_check_ref(pmass).data;
945 Tscal hfactd_ = shambase::get_check_ref(hfactd).data;
946 u32 n_sinks = eos_config->n_sinks;
947
948 Tscal inv_r0_q = 1. / sycl::pow(r0, q_);
949
950 using EOS = shamphys::EOS_LocallyIsothermal<Tscal>;
951
952 auto &all_sink_pos = get_sink_pos<Tvec>(scheduler().synchronized_data);
953 auto &all_sink_mass = get_sink_mass<Tvec>(scheduler().synchronized_data);
954 std::vector<Tvec> sink_pos;
955 std::vector<Tscal> sink_mass;
956 u32 sink_cnt = 0;
957
958 for (size_t i = 0; i < all_sink_pos.size(); i++) {
959 sink_pos.push_back(all_sink_pos[i]);
960 sink_mass.push_back(all_sink_mass[i]);
961 sink_cnt++;
962 if (sink_pos.size() >= n_sinks) { // We only consider the first n_sinks sinks
963 break;
964 }
965 }
966
967 if (sink_cnt == 0) {
968 throw shambase::make_except_with_loc<std::runtime_error>(
969 "No sinks found for the equation of state");
970 }
971
972 shamrock::solvergraph::FieldRefs<Tvec> xyz_refs{"", ""};
973 auto refs
974 = storage.merged_xyzh.get()
975 .template map<shamrock::solvergraph::PatchDataFieldRef<Tvec>>(
976 [&](u64 id,
977 PatchDataLayer &mpdat) -> shamrock::solvergraph::PatchDataFieldRef<Tvec> {
978 return mpdat.get_field<Tvec>(0);
979 });
980 xyz_refs.set_refs(refs);
981
982 sham::DeviceBuffer<Tvec> sink_pos_buf(sink_pos.size(), dev_sched);
983 sham::DeviceBuffer<Tscal> sink_mass_buf(sink_mass.size(), dev_sched);
984
985 sink_pos_buf.copy_from_stdvec(sink_pos);
986 sink_mass_buf.copy_from_stdvec(sink_mass);
987
988 auto eos_internal = [](Tvec R,
989 Tscal rho_a,
990 u32 scount,
991 auto spos,
992 auto smass,
993 Tscal cs0,
994 Tscal inv_r0_q,
995 Tscal q,
996 Tscal &pressure,
997 Tscal &soundspeed) {
998 Tscal sink_mass_sum = 0;
999 Tscal pot_sum = 0;
1000 for (u32 i = 0; i < scount; i++) {
1001 Tvec s_r = spos[i] - R;
1002 Tscal s_m = smass[i];
1003 Tscal s_r_abs = sycl::length(s_r);
1004 sink_mass_sum += s_m;
1005 pot_sum += s_m / s_r_abs;
1006 }
1007
1008 Tscal cs_out = cs0 * inv_r0_q * sycl::pow(pot_sum / sink_mass_sum, q);
1009 Tscal P_a = EOS::pressure_from_cs(cs_out * cs_out, rho_a);
1010
1011 pressure = P_a;
1012 soundspeed = cs_out;
1013 };
1014
1015 if (has_rho) {
1016 auto &spans_rho_ = edges.spans_rho.value().get();
1017 spans_rho_.check_sizes(sizes_indexes);
1018
1019 sizes_indexes.for_each([&](u64 id, u32 count) {
1020 sham::kernel_call(
1021 q,
1022 sham::MultiRef{
1023 spans_rho_.get_spans().get(id),
1024 xyz_refs.get_spans().get(id),
1025 sink_pos_buf,
1026 sink_mass_buf},
1027 out_refs.get(id),
1028 count,
1029 [cs0, inv_r0_q, q_, sink_cnt, eos_internal](
1030 u32 gid,
1031 const Tscal *rho,
1032 const Tvec *xyz,
1033 const Tvec *spos,
1034 const Tscal *smass,
1035 Tscal *pressure,
1036 Tscal *soundspeed) {
1037 Tvec R_a = xyz[gid];
1038 Tscal rho_a = rho[gid];
1039 eos_internal(
1040 R_a,
1041 rho_a,
1042 sink_cnt,
1043 spos,
1044 smass,
1045 cs0,
1046 inv_r0_q,
1047 q_,
1048 pressure[gid],
1049 soundspeed[gid]);
1050 });
1051 });
1052 } else if (has_h) {
1053 auto &spans_h_ = edges.spans_h.value().get();
1054 spans_h_.check_sizes(sizes_indexes);
1055
1056 sizes_indexes.for_each([&](u64 id, u32 count) {
1057 sham::kernel_call(
1058 q,
1059 sham::MultiRef{
1060 spans_h_.get_spans().get(id),
1061 xyz_refs.get_spans().get(id),
1062 sink_pos_buf,
1063 sink_mass_buf},
1064 out_refs.get(id),
1065 count,
1066 [cs0, inv_r0_q, q_, sink_cnt, pmass_, hfactd_, eos_internal](
1067 u32 gid,
1068 const Tscal *h,
1069 const Tvec *xyz,
1070 const Tvec *spos,
1071 const Tscal *smass,
1072 Tscal *pressure,
1073 Tscal *soundspeed) {
1074 using namespace shamrock::sph;
1075 Tvec R_a = xyz[gid];
1076 Tscal rho_a = rho_h(pmass_, h[gid], hfactd_);
1077 eos_internal(
1078 R_a,
1079 rho_a,
1080 sink_cnt,
1081 spos,
1082 smass,
1083 cs0,
1084 inv_r0_q,
1085 q_,
1086 pressure[gid],
1087 soundspeed[gid]);
1088 });
1089 });
1090 }
1091
1092 } else if (
1093 SolverEOS_Fermi *eos_config
1094 = std::get_if<SolverEOS_Fermi>(&solver_config.eos_config.config)) {
1095
1096 using namespace shamunits;
1097 auto unit_sys = *solver_config.unit_sys;
1098
1099 Tscal mass = unit_sys.template to<units::kilogram>();
1100 Tscal length = unit_sys.template to<units::metre>();
1101 Tscal time = unit_sys.template to<units::second>();
1102
1103 auto mu_e = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("mu_e", "\\mu_e");
1104 mu_e->data = eos_config->mu_e;
1105
1106 auto density_unit
1107 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("density_unit", "\\rho_u");
1108 density_unit->data = mass / (length * length * length);
1109
1110 auto pressure_unit
1111 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("pressure_unit", "P_u");
1112 pressure_unit->data = mass / length / (time * time);
1113
1114 auto velocity_unit
1115 = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("velocity_unit", "v_u");
1116 velocity_unit->data = length / time;
1117
1118 shammodels::common::modules::ComputeEOSFermi<Tvec> node;
1119 node.set_edges(
1120 mu_e,
1121 density_unit,
1122 pressure_unit,
1123 velocity_unit,
1124 hfactd,
1125 pmass,
1126 spans_rho,
1127 spans_h,
1128 sizes,
1129 spans_pressure,
1130 spans_soundspeed);
1131 node.evaluate();
1132 } else {
1133 shambase::throw_unimplemented();
1134 }
1135}
1136
1137template<class Tvec, template<class> class SPHKernel>
1138void shammodels::sph::modules::ComputeEos<Tvec, SPHKernel>::compute_eos() {
1139
1140 NamedStackEntry stack_loc{"compute eos"};
1141
1142 Tscal gpart_mass = solver_config.gpart_mass;
1143
1144 using namespace shamrock;
1145 using namespace shamrock::patch;
1146
1147 shamrock::patch::PatchDataLayerLayout &ghost_layout
1148 = shambase::get_check_ref(storage.ghost_layout.get());
1149 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
1150 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
1151
1152 auto hfactd = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("hfactd", "hfactd");
1153 auto pmass = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("pmass", "pmass");
1154
1155 hfactd->data = Kernel::hfactd;
1156 pmass->data = gpart_mass;
1157
1158 auto sizes = storage.part_counts_with_ghost;
1159
1160 auto h_refs = shamrock::solvergraph::FieldRefs<Tscal>::make_shared("", "");
1161 {
1162 auto refs = storage.merged_patchdata_ghost.get()
1163 .template map<shamrock::solvergraph::PatchDataFieldRef<Tscal>>(
1164 [&](u64 id, PatchDataLayer &mpdat)
1165 -> shamrock::solvergraph::PatchDataFieldRef<Tscal> {
1166 return mpdat.get_field<Tscal>(ihpart_interf);
1167 });
1168 h_refs->set_refs(refs);
1169 }
1170
1171 auto uint_refs = shamrock::solvergraph::FieldRefs<Tscal>::make_shared("", "");
1172 {
1173 auto refs = storage.merged_patchdata_ghost.get()
1174 .template map<shamrock::solvergraph::PatchDataFieldRef<Tscal>>(
1175 [&](u64 id, PatchDataLayer &mpdat)
1176 -> shamrock::solvergraph::PatchDataFieldRef<Tscal> {
1177 return mpdat.get_field<Tscal>(iuint_interf);
1178 });
1179 uint_refs->set_refs(refs);
1180 }
1181
1182 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
1183
1184 if (solver_config.dust_config.has_epsilon_field()) {
1185
1186 u32 iepsilon_interf = ghost_layout.get_field_idx<Tscal>("epsilon");
1187 u32 nvar_dust = solver_config.dust_config.get_dust_nvar();
1188
1189 auto rho_g = std::make_shared<shamrock::solvergraph::Field<Tscal>>(1, "rho_g", "rho_g");
1190 auto uint_g = std::make_shared<shamrock::solvergraph::Field<Tscal>>(1, "uint_g", "uint_g");
1191
1192 rho_g->ensure_sizes(shambase::get_check_ref(sizes).indexes);
1193 uint_g->ensure_sizes(shambase::get_check_ref(sizes).indexes);
1194
1195 shamrock::solvergraph::FieldRefs<Tscal> epsilon_refs{"", ""};
1196 auto refs = storage.merged_patchdata_ghost.get()
1197 .template map<shamrock::solvergraph::PatchDataFieldRef<Tscal>>(
1198 [&](u64 id, PatchDataLayer &mpdat)
1199 -> shamrock::solvergraph::PatchDataFieldRef<Tscal> {
1200 return mpdat.get_field<Tscal>(iepsilon_interf);
1201 });
1202 epsilon_refs.set_refs(refs);
1203
1204 sham::distributed_data_kernel_call(
1205 dev_sched,
1206 sham::DDMultiRef{h_refs->get_spans(), uint_refs->get_spans(), epsilon_refs.get_spans()},
1207 sham::DDMultiRef{rho_g->get_spans(), uint_g->get_spans()},
1208 shambase::get_check_ref(sizes).indexes,
1209 [pmass = pmass->data, hfactd = hfactd->data, nvar_dust](
1210 u32 gid,
1211 const Tscal *h,
1212 const Tscal *uint,
1213 const Tscal *epsilon,
1214 Tscal *rho_g,
1215 Tscal *uint_g) {
1216 using namespace shamrock::sph;
1217 Tscal rho_a = rho_h(pmass, h[gid], hfactd);
1218 Tscal uint_a = uint[gid];
1219
1220 Tscal epsilon_sum = 0;
1221 for (u32 j = 0; j < nvar_dust; j++) {
1222 epsilon_sum += epsilon[gid * nvar_dust + j];
1223 }
1224
1225 Tscal rho_g_a = rho_a * (1 - epsilon_sum);
1226 Tscal uint_g_a = uint_a / (1 - epsilon_sum);
1227
1228 rho_g[gid] = rho_g_a;
1229 uint_g[gid] = uint_g_a;
1230 });
1231
1232 compute_eos_internal(
1233 hfactd,
1234 pmass,
1235 rho_g,
1236 std::nullopt,
1237 uint_g,
1238 sizes,
1239 storage.pressure,
1240 storage.soundspeed);
1241 } else if (solver_config.dust_config.has_s_j_field()) {
1242
1243 u32 is_j_interf = ghost_layout.get_field_idx<Tscal>("s_j");
1244 u32 nvar_dust = solver_config.dust_config.get_dust_nvar();
1245
1246 auto rho_g = std::make_shared<shamrock::solvergraph::Field<Tscal>>(1, "rho_g", "rho_g");
1247 auto uint_g = std::make_shared<shamrock::solvergraph::Field<Tscal>>(1, "uint_g", "uint_g");
1248
1249 rho_g->ensure_sizes(shambase::get_check_ref(sizes).indexes);
1250 uint_g->ensure_sizes(shambase::get_check_ref(sizes).indexes);
1251
1252 shamrock::solvergraph::FieldRefs<Tscal> s_j_refs{"", ""};
1253 auto refs = storage.merged_patchdata_ghost.get()
1254 .template map<shamrock::solvergraph::PatchDataFieldRef<Tscal>>(
1255 [&](u64 id, PatchDataLayer &mpdat)
1256 -> shamrock::solvergraph::PatchDataFieldRef<Tscal> {
1257 return mpdat.get_field<Tscal>(is_j_interf);
1258 });
1259 s_j_refs.set_refs(refs);
1260
1261 sham::distributed_data_kernel_call(
1262 dev_sched,
1263 sham::DDMultiRef{h_refs->get_spans(), uint_refs->get_spans(), s_j_refs.get_spans()},
1264 sham::DDMultiRef{rho_g->get_spans(), uint_g->get_spans()},
1265 shambase::get_check_ref(sizes).indexes,
1266 [pmass = pmass->data, hfactd = hfactd->data, nvar_dust](
1267 u32 gid,
1268 const Tscal *h,
1269 const Tscal *uint,
1270 const Tscal *s_j,
1271 Tscal *rho_g,
1272 Tscal *uint_g) {
1273 using namespace shamrock::sph;
1274 Tscal rho_a = rho_h(pmass, h[gid], hfactd);
1275 Tscal uint_a = uint[gid];
1276
1277 Tscal epsilon_sum = 0;
1278 for (u32 j = 0; j < nvar_dust; j++) {
1279 Tscal s = s_j[gid * nvar_dust + j];
1280 epsilon_sum += s * s / rho_a;
1281 }
1282
1283 Tscal rho_g_a = rho_a * (1 - epsilon_sum);
1284 Tscal uint_g_a = uint_a / (1 - epsilon_sum);
1285
1286 rho_g[gid] = rho_g_a;
1287 uint_g[gid] = uint_g_a;
1288 });
1289
1290 compute_eos_internal(
1291 hfactd,
1292 pmass,
1293 rho_g,
1294 std::nullopt,
1295 uint_g,
1296 sizes,
1297 storage.pressure,
1298 storage.soundspeed);
1299 } else {
1300
1301 compute_eos_internal(
1302 hfactd,
1303 pmass,
1304 std::nullopt,
1305 h_refs,
1306 uint_refs,
1307 sizes,
1308 storage.pressure,
1309 storage.soundspeed);
1310 }
1311}
1312
1313using namespace shammath;
1314template class shammodels::sph::modules::ComputeEos<f64_3, M4>;
1315template class shammodels::sph::modules::ComputeEos<f64_3, M6>;
1316template class shammodels::sph::modules::ComputeEos<f64_3, M8>;
1317
1318template class shammodels::sph::modules::ComputeEos<f64_3, C2>;
1319template class shammodels::sph::modules::ComputeEos<f64_3, C4>;
1320template class shammodels::sph::modules::ComputeEos<f64_3, C6>;
shammodels::gsph::names::common::xyz
constexpr const char * xyz
Position field (3D coordinates).
Definition FieldNames.hpp:33
shammodels::gsph::names::newtonian::soundspeed
constexpr const char * soundspeed
Sound speed c_s (derived from EOS).
Definition FieldNames.hpp:61
shammodels::gsph::names::newtonian::pressure
constexpr const char * pressure
Pressure P (derived from EOS).
Definition FieldNames.hpp:58
shamrock::solvergraph::PatchDataFieldRef
std::reference_wrapper< PatchDataField< T > > PatchDataFieldRef
Alias for a reference to a PatchDataField.
Definition IFieldRefs.hpp:28
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::DeviceQueue
A SYCL queue associated with a device and a context.
Definition DeviceQueue.hpp:32
shammodels::common::modules::ComputeEOSAdiabatic::_impl_get_label
virtual std::string _impl_get_label() const
get the label of the node
Definition ComputeEos.cpp:244
shammodels::common::modules::ComputeEOSAdiabatic::_impl_get_tex
virtual std::string _impl_get_tex() const
get the tex of the node
Definition ComputeEos.cpp:246
shammodels::common::modules::ComputeEOSFermi::_impl_get_label
virtual std::string _impl_get_label() const
get the label of the node
Definition ComputeEos.cpp:582
shammodels::common::modules::ComputeEOSFermi::_impl_get_tex
virtual std::string _impl_get_tex() const
get the tex of the node
Definition ComputeEos.cpp:584
shammodels::common::modules::ComputeEOSFermi::_impl_evaluate_internal
void _impl_evaluate_internal()
evaluate the node
Definition ComputeEos.cpp:506
shammodels::common::modules::ComputeEOSIsothermal::_impl_get_label
virtual std::string _impl_get_label() const
get the label of the node
Definition ComputeEos.cpp:129
shammodels::common::modules::ComputeEOSIsothermal::_impl_get_tex
virtual std::string _impl_get_tex() const
get the tex of the node
Definition ComputeEos.cpp:131
shammodels::common::modules::ComputeEOSLocallyIsothermal::_impl_get_tex
virtual std::string _impl_get_tex() const
get the tex of the node
Definition ComputeEos.cpp:459
shammodels::common::modules::ComputeEOSLocallyIsothermal::_impl_get_label
virtual std::string _impl_get_label() const
get the label of the node
Definition ComputeEos.cpp:455
shammodels::common::modules::ComputeEOSPolytropic::_impl_get_label
virtual std::string _impl_get_label() const
get the label of the node
Definition ComputeEos.cpp:349
shammodels::common::modules::ComputeEOSPolytropic::_impl_get_tex
virtual std::string _impl_get_tex() const
get the tex of the node
Definition ComputeEos.cpp:351
shammodels::sph::modules::ComputeEos
Module for computing equation of state quantities.
Definition ComputeEos.hpp:39
shammodels::sph::modules::ComputeEos::compute_eos
void compute_eos()
Computes pressure and sound speed from equation of state.
Definition ComputeEos.cpp:1138
shamrock::patch::PatchDataLayerLayout::get_field_idx
u32 get_field_idx(const std::string &field_name) const
Get the field id if matching name & type.
Definition PatchDataLayerLayout.cpp:130
shamrock::patch::PatchDataLayer
PatchDataLayer container class, the layout is described in patchdata_layout.
Definition PatchDataLayer.hpp:35
shamrock::solvergraph::FieldRefs::get_spans
virtual DDPatchDataFieldSpanPointer< T > & get_spans()
Get the DistributedData of spans attached to the underlying field.
Definition FieldRefs.hpp:49
shamrock::solvergraph::IFieldSpan
Interface for a solver graph edge representing a field as spans.
Definition IFieldSpan.hpp:37
shamrock::solvergraph::INode
Inode is node between data edges, takes multiple inputs, multiple outputs.
Definition INode.hpp:31
shamrock::solvergraph::INode::evaluate
void evaluate()
Evaluate the node.
Definition INode.hpp:156
exception.hpp
This header file contains utility functions related to exception handling in the code.
kernel_call.hpp
kernel_call_distrib.hpp
sham::distributed_data_kernel_call
void distributed_data_kernel_call(sham::DeviceScheduler_ptr dev_sched, RefIn in, RefOut in_out, const shambase::DistributedData< index_t > &thread_counts, Functor &&func)
A variant of sham::kernel_call for distributed data.
Definition kernel_call_distrib.hpp:79
sham::kernel_call
void kernel_call(sham::DeviceQueue &q, RefIn in, RefOut in_out, u32 n, Functor &&func, SourceLocation &&callsite=SourceLocation{})
Submit a kernel to a SYCL queue.
Definition kernel_call.hpp:344
shambase::get_check_ref
T & get_check_ref(const std::unique_ptr< T > &ptr, SourceLocation loc=SourceLocation())
Takes a std::unique_ptr and returns a reference to the object it holds. It throws a std::runtime_erro...
Definition memory.hpp:112
shambase::make_except_with_loc
ExcptTypes make_except_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Create an exception with a message and a location.
Definition exception.hpp:105
shambase::throw_unimplemented
void throw_unimplemented(SourceLocation loc=SourceLocation{})
Throw a std::runtime_error saying that the function is unimplemented.
Definition exception.hpp:153
shammath
namespace for math utility
Definition AABB.hpp:26
shamrock
namespace for the main framework
Definition __init__.py:1
shamunits
namespace containing the units library
Definition units_json.hpp:24
narrowing.hpp
Utilities for safe type narrowing conversions.
sink_edges_helper.hpp
Helpers to access SPH sink particles stored as SoA synchronized data edges.
sphkernels.hpp
sph kernels
__shamrock_stack_entry
#define __shamrock_stack_entry()
Macro to create a stack entry.
Definition stacktrace.hpp:257
NamedStackEntry
shambase::details::NamedBasicStackEntry NamedStackEntry
Alias for shambase::details::NamedBasicStackEntry.
Definition stacktrace.hpp:249
sham::DDMultiRef
A variant of sham::MultiRef for distributed data.
Definition kernel_call_distrib.hpp:33
sham::MultiRef
A class that references multiple buffers or similar objects.
Definition MultiRef.hpp:33
shamphys::EOS_Adiabatic
Adiabatic equation of state.
Definition eos.hpp:45
shamphys::EOS_Fermi
Fermi Gas EoS.
Definition eos.hpp:196
shamphys::EOS_Isothermal
Isothermal equation of state.
Definition eos.hpp:32
shamphys::EOS_LocallyIsothermal
Locally isothermal equation of state with radial dependence.
Definition eos.hpp:87
shamphys::EOS_Polytropic
Polytropic equation of state.
Definition eos.hpp:66
Generated on for Shamrock by doxygen 1.15.0