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/narrowing.hpp"
22#include "shambackends/kernel_call.hpp"
23#include "shambackends/kernel_call_distrib.hpp"
24#include "shammath/sphkernels.hpp"
25#include "shammodels/sph/math/density.hpp"
26#include "shammodels/sph/modules/ComputeEos.hpp"
27#include "shammodels/sph/sink_edges_helper.hpp"
28#include "shamphys/eos.hpp"
29#include "shamrock/patch/PatchDataField.hpp"
30#include "shamrock/patch/PatchDataLayer.hpp"
31#include "shamsolvergraph/edge/IDataEdge.hpp"
32
33template<class Tvec, template<class> class SPHKernel>
34void shammodels::sph::modules::ComputeEos<Tvec, SPHKernel>::compute_eos_internal(
35 const shamrock::solvergraph::IDataEdge<Tscal> &hfactd,
36 const shamrock::solvergraph::IDataEdge<Tscal> &pmass,
37 const std::optional<std::reference_wrapper<const shamrock::solvergraph::IFieldSpan<Tscal>>>
38 spans_rho,
39 const std::optional<std::reference_wrapper<const shamrock::solvergraph::IFieldSpan<Tscal>>>
40 spans_h,
41 const std::optional<std::reference_wrapper<const shamrock::solvergraph::IFieldSpan<Tscal>>>
42 spans_uint,
43 const shamrock::solvergraph::Indexes<u32> &sizes,
44 shamrock::solvergraph::IFieldSpan<Tscal> &spans_pressure,
45 shamrock::solvergraph::IFieldSpan<Tscal> &spans_soundspeed) {
46
47 shamrock::patch::PatchDataLayerLayout &ghost_layout
48 = shambase::get_check_ref(storage.ghost_layout.get());
49
50 using namespace shamrock;
51 using namespace shamrock::patch;
52
53 using SolverConfigEOS = typename Config::EOSConfig;
54 using SolverEOS_Isothermal = typename SolverConfigEOS::Isothermal;
55 using SolverEOS_Adiabatic = typename SolverConfigEOS::Adiabatic;
56 using SolverEOS_Polytropic = typename SolverConfigEOS::Polytropic;
57 using SolverEOS_LocallyIsothermal = typename SolverConfigEOS::LocallyIsothermal;
58 using SolverEOS_LocallyIsothermalLP07 = typename SolverConfigEOS::LocallyIsothermalLP07;
59 using SolverEOS_LocallyIsothermalFA2014 = typename SolverConfigEOS::LocallyIsothermalFA2014;
60 using SolverEOS_LocallyIsothermalFA2014Extended =
61 typename SolverConfigEOS::LocallyIsothermalFA2014Extended;
62 using SolverEOS_Fermi = typename SolverConfigEOS::Fermi;
63
64 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
65 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
66
67 shambase::get_check_ref(storage.pressure).ensure_sizes(sizes.indexes);
68 shambase::get_check_ref(storage.soundspeed).ensure_sizes(sizes.indexes);
69
70 bool has_rho = spans_rho.has_value();
71 bool has_h = spans_h.has_value();
72
73 // must have either rho or h
74 if ((!has_rho || has_h) && (has_rho || !has_h)) {
75 throw shambase::make_except_with_loc<std::invalid_argument>("Must have either rho or h");
76 }
77
78 auto out_refs = sham::DDMultiRef{spans_pressure.get_spans(), spans_soundspeed.get_spans()};
79
80 if (SolverEOS_Isothermal *eos_config
81 = std::get_if<SolverEOS_Isothermal>(&solver_config.eos_config.config)) {
82
83 Tscal cs = eos_config->cs;
84 Tscal pmass_ = pmass.data;
85 Tscal hfactd_ = hfactd.data;
86
87 using EOS = shamphys::EOS_Isothermal<Tscal>;
88
89 auto eos_internal = [](Tscal cs, Tscal rho, Tscal &pressure, Tscal &soundspeed) {
90 using EOS = shamphys::EOS_Isothermal<Tscal>;
91 pressure = EOS::pressure(cs, rho);
92 soundspeed = cs;
93 };
94
95 if (has_rho) {
96 auto &spans_rho_ = spans_rho.value().get();
97 spans_rho_.check_sizes(sizes.indexes);
98 sham::distributed_data_kernel_call(
99 dev_sched,
100 sham::DDMultiRef{spans_rho_.get_spans()},
101 out_refs,
102 sizes.indexes,
103 [cs, eos_internal](u32 gid, const Tscal *rho, Tscal *pressure, Tscal *soundspeed) {
104 Tscal rho_a = rho[gid];
105 eos_internal(cs, rho_a, pressure[gid], soundspeed[gid]);
106 });
107 } else if (has_h) {
108 auto &spans_h_ = spans_h.value().get();
109 spans_h_.check_sizes(sizes.indexes);
110 sham::distributed_data_kernel_call(
111 dev_sched,
112 sham::DDMultiRef{spans_h_.get_spans()},
113 out_refs,
114 sizes.indexes,
115 [cs, pmass_, hfactd_, eos_internal](
116 u32 gid, const Tscal *h, Tscal *pressure, Tscal *soundspeed) {
117 using namespace shamrock::sph;
118 Tscal rho = rho_h(pmass_, h[gid], hfactd_);
119 eos_internal(cs, rho, pressure[gid], soundspeed[gid]);
120 });
121 }
122 } else if (
123 SolverEOS_Adiabatic *eos_config
124 = std::get_if<SolverEOS_Adiabatic>(&solver_config.eos_config.config)) {
125
126 Tscal gamma = eos_config->gamma;
127 Tscal pmass_ = pmass.data;
128 Tscal hfactd_ = hfactd.data;
129
130 using EOS = shamphys::EOS_Adiabatic<Tscal>;
131
132 auto eos_internal
133 = [](Tscal gamma, Tscal rho, Tscal uint, Tscal &pressure, Tscal &soundspeed) {
134 using EOS = shamphys::EOS_Adiabatic<Tscal>;
135 Tscal P_a = EOS::pressure(gamma, rho, uint);
136 Tscal cs_a = EOS::cs_from_p(gamma, rho, P_a);
137 pressure = P_a;
138 soundspeed = cs_a;
139 };
140
141 if (has_rho) {
142 auto &spans_rho_ = spans_rho.value().get();
143 spans_rho_.check_sizes(sizes.indexes);
144 sham::distributed_data_kernel_call(
145 dev_sched,
146 sham::DDMultiRef{spans_rho_.get_spans(), spans_uint.value().get().get_spans()},
147 out_refs,
148 sizes.indexes,
149 [gamma, pmass_, hfactd_, eos_internal](
150 u32 gid,
151 const Tscal *rho,
152 const Tscal *uint,
153 Tscal *pressure,
154 Tscal *soundspeed) {
155 Tscal rho_a = rho[gid];
156 Tscal uint_a = uint[gid];
157 eos_internal(gamma, rho_a, uint_a, pressure[gid], soundspeed[gid]);
158 });
159 } else if (has_h) {
160 auto &spans_h_ = spans_h.value().get();
161 spans_h_.check_sizes(sizes.indexes);
162 sham::distributed_data_kernel_call(
163 dev_sched,
164 sham::DDMultiRef{spans_h_.get_spans(), spans_uint.value().get().get_spans()},
165 out_refs,
166 sizes.indexes,
167 [gamma, pmass_, hfactd_, eos_internal](
168 u32 gid,
169 const Tscal *h,
170 const Tscal *uint,
171 Tscal *pressure,
172 Tscal *soundspeed) {
173 using namespace shamrock::sph;
174 Tscal rho = rho_h(pmass_, h[gid], hfactd_);
175 Tscal uint_a = uint[gid];
176 eos_internal(gamma, rho, uint_a, pressure[gid], soundspeed[gid]);
177 });
178 }
179
180 } else if (
181 SolverEOS_Polytropic *eos_config
182 = std::get_if<SolverEOS_Polytropic>(&solver_config.eos_config.config)) {
183
184 using EOS = shamphys::EOS_Polytropic<Tscal>;
185
186 Tscal K = eos_config->K;
187 Tscal gamma = eos_config->gamma;
188 Tscal pmass_ = pmass.data;
189 Tscal hfactd_ = hfactd.data;
190
191 auto eos_internal
192 = [](Tscal K, Tscal gamma, Tscal rho_a, Tscal &pressure, Tscal &soundspeed) {
193 using EOS = shamphys::EOS_Polytropic<Tscal>;
194 Tscal P_a = EOS::pressure(gamma, K, rho_a);
195 Tscal cs_a = EOS::soundspeed(gamma, K, rho_a);
196 pressure = P_a;
197 soundspeed = cs_a;
198 };
199
200 if (has_rho) {
201 auto &spans_rho_ = spans_rho.value().get();
202 spans_rho_.check_sizes(sizes.indexes);
203 sham::distributed_data_kernel_call(
204 dev_sched,
205 sham::DDMultiRef{spans_rho_.get_spans()},
206 out_refs,
207 sizes.indexes,
208 [K, gamma, pmass_, hfactd_, eos_internal](
209 u32 gid, const Tscal *rho, Tscal *pressure, Tscal *soundspeed) {
210 Tscal rho_a = rho[gid];
211 eos_internal(K, gamma, rho_a, pressure[gid], soundspeed[gid]);
212 });
213 } else if (has_h) {
214 auto &spans_h_ = spans_h.value().get();
215 spans_h_.check_sizes(sizes.indexes);
216 sham::distributed_data_kernel_call(
217 dev_sched,
218 sham::DDMultiRef{spans_h_.get_spans()},
219 out_refs,
220 sizes.indexes,
221 [K, gamma, pmass_, hfactd_, eos_internal](
222 u32 gid, const Tscal *h, Tscal *pressure, Tscal *soundspeed) {
223 using namespace shamrock::sph;
224 Tscal rho = rho_h(pmass_, h[gid], hfactd_);
225 eos_internal(K, gamma, rho, pressure[gid], soundspeed[gid]);
226 });
227 }
228
229 } else if (
230 SolverEOS_LocallyIsothermal *eos_config
231 = std::get_if<SolverEOS_LocallyIsothermal>(&solver_config.eos_config.config)) {
232
233 using EOS = shamphys::EOS_LocallyIsothermal<Tscal>;
234
235 u32 isoundspeed_interf = ghost_layout.get_field_idx<Tscal>("soundspeed");
236
237 shamrock::solvergraph::FieldRefs<Tscal> soundspeed_refs{"", ""};
238 auto refs = storage.merged_patchdata_ghost.get()
239 .template map<shamrock::solvergraph::PatchDataFieldRef<Tscal>>(
240 [&](u64 id, PatchDataLayer &mpdat)
241 -> shamrock::solvergraph::PatchDataFieldRef<Tscal> {
242 return mpdat.get_field<Tscal>(isoundspeed_interf);
243 });
244 soundspeed_refs.set_refs(refs);
245
246 Tscal pmass_ = pmass.data;
247 Tscal hfactd_ = hfactd.data;
248
249 auto eos_internal = [](Tscal cs0, Tscal rho_a, Tscal &pressure, Tscal &soundspeed) {
250 using EOS = shamphys::EOS_LocallyIsothermal<Tscal>;
251 pressure = EOS::pressure_from_cs(cs0 * cs0, rho_a);
252 soundspeed = cs0;
253 };
254
255 if (has_rho) {
256 auto &spans_rho_ = spans_rho.value().get();
257 spans_rho_.check_sizes(sizes.indexes);
258 sham::distributed_data_kernel_call(
259 dev_sched,
260 sham::DDMultiRef{spans_rho_.get_spans(), soundspeed_refs.get_spans()},
261 out_refs,
262 sizes.indexes,
263 [eos_internal](
264 u32 gid,
265 const Tscal *rho,
266 const Tscal *cs0,
267 Tscal *pressure,
268 Tscal *soundspeed) {
269 Tscal rho_a = rho[gid];
270 Tscal cs0_a = cs0[gid];
271 eos_internal(cs0_a, rho_a, pressure[gid], soundspeed[gid]);
272 });
273 } else if (has_h) {
274 auto &spans_h_ = spans_h.value().get();
275 spans_h_.check_sizes(sizes.indexes);
276 sham::distributed_data_kernel_call(
277 dev_sched,
278 sham::DDMultiRef{spans_h_.get_spans(), soundspeed_refs.get_spans()},
279 out_refs,
280 sizes.indexes,
281 [pmass_, hfactd_, eos_internal](
282 u32 gid, const Tscal *h, const Tscal *cs0, Tscal *pressure, Tscal *soundspeed) {
283 using namespace shamrock::sph;
284 Tscal rho = rho_h(pmass_, h[gid], hfactd_);
285 Tscal cs0_a = cs0[gid];
286 eos_internal(cs0_a, rho, pressure[gid], soundspeed[gid]);
287 });
288 }
289
290 } else if (
291 SolverEOS_LocallyIsothermalLP07 *eos_config
292 = std::get_if<SolverEOS_LocallyIsothermalLP07>(&solver_config.eos_config.config)) {
293
294 Tscal cs0 = eos_config->cs0;
295 Tscal r0sq = eos_config->r0 * eos_config->r0;
296 Tscal mq = -eos_config->q;
297
298 Tscal pmass_ = pmass.data;
299 Tscal hfactd_ = hfactd.data;
300
301 shamrock::solvergraph::FieldRefs<Tvec> xyz_refs{"", ""};
302 auto refs
303 = storage.merged_xyzh.get()
304 .template map<shamrock::solvergraph::PatchDataFieldRef<Tvec>>(
305 [&](u64 id,
306 PatchDataLayer &mpdat) -> shamrock::solvergraph::PatchDataFieldRef<Tvec> {
307 return mpdat.get_field<Tvec>(0);
308 });
309 xyz_refs.set_refs(refs);
310
311 using EOS = shamphys::EOS_LocallyIsothermal<Tscal>;
312
313 auto eos_internal = [](Tvec R,
314 Tscal cs0,
315 Tscal r0sq,
316 Tscal mq,
317 Tscal rho_a,
318 Tscal &pressure,
319 Tscal &soundspeed) {
320 Tscal Rsq = sycl::dot(R, R);
321 Tscal cs_sq = EOS::soundspeed_sq(cs0 * cs0, Rsq / r0sq, mq);
322 Tscal cs_out = sycl::sqrt(cs_sq);
323
324 Tscal P_a = EOS::pressure_from_cs(cs_sq, rho_a);
325
326 pressure = P_a;
327 soundspeed = cs_out;
328 };
329
330 if (has_rho) {
331 auto &spans_rho_ = spans_rho.value().get();
332 spans_rho_.check_sizes(sizes.indexes);
333 sham::distributed_data_kernel_call(
334 dev_sched,
335 sham::DDMultiRef{spans_rho_.get_spans(), xyz_refs.get_spans()},
336 out_refs,
337 sizes.indexes,
338 [cs0, r0sq, mq, eos_internal](
339 u32 gid,
340 const Tscal *rho,
341 const Tvec *xyz,
342 Tscal *pressure,
343 Tscal *soundspeed) {
344 Tvec R_a = xyz[gid];
345 Tscal rho_a = rho[gid];
346 eos_internal(R_a, cs0, r0sq, mq, rho_a, pressure[gid], soundspeed[gid]);
347 });
348 } else if (has_h) {
349 auto &spans_h_ = spans_h.value().get();
350 spans_h_.check_sizes(sizes.indexes);
351 sham::distributed_data_kernel_call(
352 dev_sched,
353 sham::DDMultiRef{spans_h_.get_spans(), xyz_refs.get_spans()},
354 out_refs,
355 sizes.indexes,
356 [cs0, r0sq, mq, pmass_, hfactd_, eos_internal](
357 u32 gid, const Tscal *h, const Tvec *xyz, Tscal *pressure, Tscal *soundspeed) {
358 using namespace shamrock::sph;
359 Tvec R_a = xyz[gid];
360 Tscal rho_a = rho_h(pmass_, h[gid], hfactd_);
361 eos_internal(R_a, cs0, r0sq, mq, rho_a, pressure[gid], soundspeed[gid]);
362 });
363 }
364
365 } else if (
366 SolverEOS_LocallyIsothermalFA2014 *eos_config
367 = std::get_if<SolverEOS_LocallyIsothermalFA2014>(&solver_config.eos_config.config)) {
368
369 Tscal G = solver_config.get_constant_G();
370 Tscal h_over_r = eos_config->h_over_r;
371 Tscal pmass_ = pmass.data;
372 Tscal hfactd_ = hfactd.data;
373
374 using EOS = shamphys::EOS_LocallyIsothermal<Tscal>;
375
376 auto &sink_pos = get_sink_pos<Tvec>(scheduler().synchronized_data);
377 auto &sink_mass = get_sink_mass<Tvec>(scheduler().synchronized_data);
378 u32 sink_cnt = shambase::narrow_or_throw<u32>(sink_pos.size());
379
380 if (sink_cnt == 0) {
381 throw shambase::make_except_with_loc<std::runtime_error>(
382 "No sinks found for the equation of state");
383 }
384
385 shamrock::solvergraph::FieldRefs<Tvec> xyz_refs{"", ""};
386 auto refs
387 = storage.merged_xyzh.get()
388 .template map<shamrock::solvergraph::PatchDataFieldRef<Tvec>>(
389 [&](u64 id,
390 PatchDataLayer &mpdat) -> shamrock::solvergraph::PatchDataFieldRef<Tvec> {
391 return mpdat.get_field<Tvec>(0);
392 });
393 xyz_refs.set_refs(refs);
394
395 sham::DeviceBuffer<Tvec> sink_pos_buf(sink_pos.size(), dev_sched);
396 sham::DeviceBuffer<Tscal> sink_mass_buf(sink_mass.size(), dev_sched);
397
398 sink_pos_buf.copy_from_stdvec(sink_pos);
399 sink_mass_buf.copy_from_stdvec(sink_mass);
400
401 auto eos_internal = [](Tvec R,
402 Tscal rho_a,
403 u32 scount,
404 auto spos,
405 auto smass,
406 Tscal G,
407 Tscal h_over_r,
408 Tscal &pressure,
409 Tscal &soundspeed) {
410 Tscal mpotential = 0;
411 for (u32 i = 0; i < scount; i++) {
412 Tvec s_r = spos[i] - R;
413 Tscal s_m = smass[i];
414 Tscal s_r_abs = sycl::length(s_r);
415 mpotential += G * s_m / s_r_abs;
416 }
417
418 Tscal cs_out = h_over_r * sycl::sqrt(mpotential);
419 Tscal P_a = EOS::pressure_from_cs(cs_out * cs_out, rho_a);
420
421 pressure = P_a;
422 soundspeed = cs_out;
423 };
424
425 if (has_rho) {
426 auto &spans_rho_ = spans_rho.value().get();
427 spans_rho_.check_sizes(sizes.indexes);
428
429 sizes.indexes.for_each([&](u64 id, u32 count) {
430 sham::kernel_call(
431 q,
432 sham::MultiRef{
433 spans_rho_.get_spans().get(id),
434 xyz_refs.get_spans().get(id),
435 sink_pos_buf,
436 sink_mass_buf},
437 out_refs.get(id),
438 count,
439 [G, h_over_r, sink_cnt, eos_internal](
440 u32 gid,
441 const Tscal *rho,
442 const Tvec *xyz,
443 const Tvec *spos,
444 const Tscal *smass,
445 Tscal *pressure,
446 Tscal *soundspeed) {
447 Tvec R_a = xyz[gid];
448 Tscal rho_a = rho[gid];
449 eos_internal(
450 R_a,
451 rho_a,
452 sink_cnt,
453 spos,
454 smass,
455 G,
456 h_over_r,
457 pressure[gid],
458 soundspeed[gid]);
459 });
460 });
461
462 } else if (has_h) {
463 auto &spans_h_ = spans_h.value().get();
464 spans_h_.check_sizes(sizes.indexes);
465
466 sizes.indexes.for_each([&](u64 id, u32 count) {
467 sham::kernel_call(
468 q,
469 sham::MultiRef{
470 spans_h_.get_spans().get(id),
471 xyz_refs.get_spans().get(id),
472 sink_pos_buf,
473 sink_mass_buf},
474 out_refs.get(id),
475 count,
476 [G, h_over_r, sink_cnt, pmass_, hfactd_, eos_internal](
477 u32 gid,
478 const Tscal *h,
479 const Tvec *xyz,
480 const Tvec *spos,
481 const Tscal *smass,
482 Tscal *pressure,
483 Tscal *soundspeed) {
484 using namespace shamrock::sph;
485 Tvec R_a = xyz[gid];
486 Tscal rho_a = rho_h(pmass_, h[gid], hfactd_);
487 eos_internal(
488 R_a,
489 rho_a,
490 sink_cnt,
491 spos,
492 smass,
493 G,
494 h_over_r,
495 pressure[gid],
496 soundspeed[gid]);
497 });
498 });
499 }
500
501 } else if (
502 SolverEOS_LocallyIsothermalFA2014Extended *eos_config
503 = std::get_if<SolverEOS_LocallyIsothermalFA2014Extended>(
504 &solver_config.eos_config.config)) {
505
506 Tscal cs0 = eos_config->cs0;
507 Tscal r0 = eos_config->r0;
508 Tscal q_ = eos_config->q;
509 Tscal pmass_ = pmass.data;
510 Tscal hfactd_ = hfactd.data;
511 u32 n_sinks = eos_config->n_sinks;
512
513 Tscal inv_r0_q = 1. / sycl::pow(r0, q_);
514
515 using EOS = shamphys::EOS_LocallyIsothermal<Tscal>;
516
517 auto &all_sink_pos = get_sink_pos<Tvec>(scheduler().synchronized_data);
518 auto &all_sink_mass = get_sink_mass<Tvec>(scheduler().synchronized_data);
519 std::vector<Tvec> sink_pos;
520 std::vector<Tscal> sink_mass;
521 u32 sink_cnt = 0;
522
523 for (size_t i = 0; i < all_sink_pos.size(); i++) {
524 sink_pos.push_back(all_sink_pos[i]);
525 sink_mass.push_back(all_sink_mass[i]);
526 sink_cnt++;
527 if (sink_pos.size() >= n_sinks) { // We only consider the first n_sinks sinks
528 break;
529 }
530 }
531
532 if (sink_cnt == 0) {
533 throw shambase::make_except_with_loc<std::runtime_error>(
534 "No sinks found for the equation of state");
535 }
536
537 shamrock::solvergraph::FieldRefs<Tvec> xyz_refs{"", ""};
538 auto refs
539 = storage.merged_xyzh.get()
540 .template map<shamrock::solvergraph::PatchDataFieldRef<Tvec>>(
541 [&](u64 id,
542 PatchDataLayer &mpdat) -> shamrock::solvergraph::PatchDataFieldRef<Tvec> {
543 return mpdat.get_field<Tvec>(0);
544 });
545 xyz_refs.set_refs(refs);
546
547 sham::DeviceBuffer<Tvec> sink_pos_buf(sink_pos.size(), dev_sched);
548 sham::DeviceBuffer<Tscal> sink_mass_buf(sink_mass.size(), dev_sched);
549
550 sink_pos_buf.copy_from_stdvec(sink_pos);
551 sink_mass_buf.copy_from_stdvec(sink_mass);
552
553 auto eos_internal = [](Tvec R,
554 Tscal rho_a,
555 u32 scount,
556 auto spos,
557 auto smass,
558 Tscal cs0,
559 Tscal inv_r0_q,
560 Tscal q,
561 Tscal &pressure,
562 Tscal &soundspeed) {
563 Tscal sink_mass_sum = 0;
564 Tscal pot_sum = 0;
565 for (u32 i = 0; i < scount; i++) {
566 Tvec s_r = spos[i] - R;
567 Tscal s_m = smass[i];
568 Tscal s_r_abs = sycl::length(s_r);
569 sink_mass_sum += s_m;
570 pot_sum += s_m / s_r_abs;
571 }
572
573 Tscal cs_out = cs0 * inv_r0_q * sycl::pow(pot_sum / sink_mass_sum, q);
574 Tscal P_a = EOS::pressure_from_cs(cs_out * cs_out, rho_a);
575
576 pressure = P_a;
577 soundspeed = cs_out;
578 };
579
580 if (has_rho) {
581 auto &spans_rho_ = spans_rho.value().get();
582 spans_rho_.check_sizes(sizes.indexes);
583
584 sizes.indexes.for_each([&](u64 id, u32 count) {
585 sham::kernel_call(
586 q,
587 sham::MultiRef{
588 spans_rho_.get_spans().get(id),
589 xyz_refs.get_spans().get(id),
590 sink_pos_buf,
591 sink_mass_buf},
592 out_refs.get(id),
593 count,
594 [cs0, inv_r0_q, q_, sink_cnt, eos_internal](
595 u32 gid,
596 const Tscal *rho,
597 const Tvec *xyz,
598 const Tvec *spos,
599 const Tscal *smass,
600 Tscal *pressure,
601 Tscal *soundspeed) {
602 Tvec R_a = xyz[gid];
603 Tscal rho_a = rho[gid];
604 eos_internal(
605 R_a,
606 rho_a,
607 sink_cnt,
608 spos,
609 smass,
610 cs0,
611 inv_r0_q,
612 q_,
613 pressure[gid],
614 soundspeed[gid]);
615 });
616 });
617 } else if (has_h) {
618 auto &spans_h_ = spans_h.value().get();
619 spans_h_.check_sizes(sizes.indexes);
620
621 sizes.indexes.for_each([&](u64 id, u32 count) {
622 sham::kernel_call(
623 q,
624 sham::MultiRef{
625 spans_h_.get_spans().get(id),
626 xyz_refs.get_spans().get(id),
627 sink_pos_buf,
628 sink_mass_buf},
629 out_refs.get(id),
630 count,
631 [cs0, inv_r0_q, q_, sink_cnt, pmass_, hfactd_, eos_internal](
632 u32 gid,
633 const Tscal *h,
634 const Tvec *xyz,
635 const Tvec *spos,
636 const Tscal *smass,
637 Tscal *pressure,
638 Tscal *soundspeed) {
639 using namespace shamrock::sph;
640 Tvec R_a = xyz[gid];
641 Tscal rho_a = rho_h(pmass_, h[gid], hfactd_);
642 eos_internal(
643 R_a,
644 rho_a,
645 sink_cnt,
646 spos,
647 smass,
648 cs0,
649 inv_r0_q,
650 q_,
651 pressure[gid],
652 soundspeed[gid]);
653 });
654 });
655 }
656
657 } else if (
658 SolverEOS_Fermi *eos_config
659 = std::get_if<SolverEOS_Fermi>(&solver_config.eos_config.config)) {
660
661 using namespace shamunits;
662 auto unit_sys = *solver_config.unit_sys;
663
664 Tscal mass = unit_sys.template to<units::kilogram>();
665 Tscal length = unit_sys.template to<units::metre>();
666 Tscal time = unit_sys.template to<units::second>();
667
668 Tscal pressure_unit = mass / length / (time * time);
669 Tscal density_unit = mass / (length * length * length);
670 Tscal velocity_unit = length / time;
671
672 Tscal mu_e = eos_config->mu_e;
673
674 Tscal pmass_ = pmass.data;
675 Tscal hfactd_ = hfactd.data;
676
677 auto eos_internal = [density_unit, pressure_unit, velocity_unit](
678 Tscal mu_e, Tscal rho_a, Tscal &pressure, Tscal &soundspeed) {
679 using EOS = shamphys::EOS_Fermi<Tscal>;
680 auto const res = EOS::pressure_and_soundspeed(mu_e, rho_a * density_unit);
681 pressure = res.pressure / pressure_unit;
682 soundspeed = res.soundspeed / velocity_unit;
683 };
684
685 if (has_rho) {
686 auto &spans_rho_ = spans_rho.value().get();
687 spans_rho_.check_sizes(sizes.indexes);
688 sham::distributed_data_kernel_call(
689 dev_sched,
690 sham::DDMultiRef{spans_rho_.get_spans()},
691 out_refs,
692 sizes.indexes,
693 [mu_e,
694 eos_internal](u32 gid, const Tscal *rho, Tscal *pressure, Tscal *soundspeed) {
695 Tscal rho_a = rho[gid];
696 eos_internal(mu_e, rho_a, pressure[gid], soundspeed[gid]);
697 });
698 } else if (has_h) {
699 auto &spans_h_ = spans_h.value().get();
700 spans_h_.check_sizes(sizes.indexes);
701 sham::distributed_data_kernel_call(
702 dev_sched,
703 sham::DDMultiRef{spans_h_.get_spans()},
704 out_refs,
705 sizes.indexes,
706 [mu_e, pmass_, hfactd_, eos_internal](
707 u32 gid, const Tscal *h, Tscal *pressure, Tscal *soundspeed) {
708 using namespace shamrock::sph;
709 Tscal rho_a = rho_h(pmass_, h[gid], hfactd_);
710 eos_internal(mu_e, rho_a, pressure[gid], soundspeed[gid]);
711 });
712 }
713
714 } else {
715 shambase::throw_unimplemented();
716 }
717}
718
719template<class Tvec, template<class> class SPHKernel>
720void shammodels::sph::modules::ComputeEos<Tvec, SPHKernel>::compute_eos() {
721
722 NamedStackEntry stack_loc{"compute eos"};
723
724 Tscal gpart_mass = solver_config.gpart_mass;
725
726 using namespace shamrock;
727 using namespace shamrock::patch;
728
729 shamrock::patch::PatchDataLayerLayout &ghost_layout
730 = shambase::get_check_ref(storage.ghost_layout.get());
731 u32 ihpart_interf = ghost_layout.get_field_idx<Tscal>("hpart");
732 u32 iuint_interf = ghost_layout.get_field_idx<Tscal>("uint");
733
734 shamrock::solvergraph::IDataEdge<Tscal> hfactd("hfactd", "hfactd");
735 shamrock::solvergraph::IDataEdge<Tscal> pmass("pmass", "pmass");
736
737 hfactd.data = Kernel::hfactd;
738 pmass.data = gpart_mass;
739
740 auto &sizes = shambase::get_check_ref(storage.part_counts_with_ghost);
741
742 shamrock::solvergraph::FieldRefs<Tscal> h_refs{"", ""};
743 {
744 auto refs = storage.merged_patchdata_ghost.get()
745 .template map<shamrock::solvergraph::PatchDataFieldRef<Tscal>>(
746 [&](u64 id, PatchDataLayer &mpdat)
747 -> shamrock::solvergraph::PatchDataFieldRef<Tscal> {
748 return mpdat.get_field<Tscal>(ihpart_interf);
749 });
750 h_refs.set_refs(refs);
751 }
752
753 shamrock::solvergraph::FieldRefs<Tscal> uint_refs{"", ""};
754 {
755 auto refs = storage.merged_patchdata_ghost.get()
756 .template map<shamrock::solvergraph::PatchDataFieldRef<Tscal>>(
757 [&](u64 id, PatchDataLayer &mpdat)
758 -> shamrock::solvergraph::PatchDataFieldRef<Tscal> {
759 return mpdat.get_field<Tscal>(iuint_interf);
760 });
761 uint_refs.set_refs(refs);
762 }
763
764 auto dev_sched = shamsys::instance::get_compute_scheduler_ptr();
765
766 if (solver_config.dust_config.has_epsilon_field()) {
767
768 u32 iepsilon_interf = ghost_layout.get_field_idx<Tscal>("epsilon");
769 u32 nvar_dust = solver_config.dust_config.get_dust_nvar();
770
771 shamrock::solvergraph::Field<Tscal> rho_g
772 = shamrock::solvergraph::Field<Tscal>(1, "rho_g", "rho_g");
773 shamrock::solvergraph::Field<Tscal> uint_g
774 = shamrock::solvergraph::Field<Tscal>(1, "uint_g", "uint_g");
775
776 rho_g.ensure_sizes(sizes.indexes);
777 uint_g.ensure_sizes(sizes.indexes);
778
779 shamrock::solvergraph::FieldRefs<Tscal> epsilon_refs{"", ""};
780 auto refs = storage.merged_patchdata_ghost.get()
781 .template map<shamrock::solvergraph::PatchDataFieldRef<Tscal>>(
782 [&](u64 id, PatchDataLayer &mpdat)
783 -> shamrock::solvergraph::PatchDataFieldRef<Tscal> {
784 return mpdat.get_field<Tscal>(iepsilon_interf);
785 });
786 epsilon_refs.set_refs(refs);
787
788 sham::distributed_data_kernel_call(
789 dev_sched,
790 sham::DDMultiRef{h_refs.get_spans(), uint_refs.get_spans(), epsilon_refs.get_spans()},
791 sham::DDMultiRef{rho_g.get_spans(), uint_g.get_spans()},
792 sizes.indexes,
793 [pmass = pmass.data, hfactd = hfactd.data, nvar_dust](
794 u32 gid,
795 const Tscal *h,
796 const Tscal *uint,
797 const Tscal *epsilon,
798 Tscal *rho_g,
799 Tscal *uint_g) {
800 using namespace shamrock::sph;
801 Tscal rho_a = rho_h(pmass, h[gid], hfactd);
802 Tscal uint_a = uint[gid];
803
804 Tscal epsilon_sum = 0;
805 for (u32 j = 0; j < nvar_dust; j++) {
806 epsilon_sum += epsilon[gid * nvar_dust + j];
807 }
808
809 Tscal rho_g_a = rho_a * (1 - epsilon_sum);
810 Tscal uint_g_a = uint_a / (1 - epsilon_sum);
811
812 rho_g[gid] = rho_g_a;
813 uint_g[gid] = uint_g_a;
814 });
815
816 compute_eos_internal(
817 hfactd,
818 pmass,
819 rho_g,
820 std::nullopt,
821 uint_g,
822 sizes,
823 shambase::get_check_ref(storage.pressure),
824 shambase::get_check_ref(storage.soundspeed));
825 } else if (solver_config.dust_config.has_s_j_field()) {
826
827 u32 is_j_interf = ghost_layout.get_field_idx<Tscal>("s_j");
828 u32 nvar_dust = solver_config.dust_config.get_dust_nvar();
829
830 shamrock::solvergraph::Field<Tscal> rho_g
831 = shamrock::solvergraph::Field<Tscal>(1, "rho_g", "rho_g");
832 shamrock::solvergraph::Field<Tscal> uint_g
833 = shamrock::solvergraph::Field<Tscal>(1, "uint_g", "uint_g");
834
835 rho_g.ensure_sizes(sizes.indexes);
836 uint_g.ensure_sizes(sizes.indexes);
837
838 shamrock::solvergraph::FieldRefs<Tscal> s_j_refs{"", ""};
839 auto refs = storage.merged_patchdata_ghost.get()
840 .template map<shamrock::solvergraph::PatchDataFieldRef<Tscal>>(
841 [&](u64 id, PatchDataLayer &mpdat)
842 -> shamrock::solvergraph::PatchDataFieldRef<Tscal> {
843 return mpdat.get_field<Tscal>(is_j_interf);
844 });
845 s_j_refs.set_refs(refs);
846
847 sham::distributed_data_kernel_call(
848 dev_sched,
849 sham::DDMultiRef{h_refs.get_spans(), uint_refs.get_spans(), s_j_refs.get_spans()},
850 sham::DDMultiRef{rho_g.get_spans(), uint_g.get_spans()},
851 sizes.indexes,
852 [pmass = pmass.data, hfactd = hfactd.data, nvar_dust](
853 u32 gid,
854 const Tscal *h,
855 const Tscal *uint,
856 const Tscal *s_j,
857 Tscal *rho_g,
858 Tscal *uint_g) {
859 using namespace shamrock::sph;
860 Tscal rho_a = rho_h(pmass, h[gid], hfactd);
861 Tscal uint_a = uint[gid];
862
863 Tscal epsilon_sum = 0;
864 for (u32 j = 0; j < nvar_dust; j++) {
865 Tscal s = s_j[gid * nvar_dust + j];
866 epsilon_sum += s * s / rho_a;
867 }
868
869 Tscal rho_g_a = rho_a * (1 - epsilon_sum);
870 Tscal uint_g_a = uint_a / (1 - epsilon_sum);
871
872 rho_g[gid] = rho_g_a;
873 uint_g[gid] = uint_g_a;
874 });
875
876 compute_eos_internal(
877 hfactd,
878 pmass,
879 rho_g,
880 std::nullopt,
881 uint_g,
882 sizes,
883 shambase::get_check_ref(storage.pressure),
884 shambase::get_check_ref(storage.soundspeed));
885 } else {
886
887 compute_eos_internal(
888 hfactd,
889 pmass,
890 std::nullopt,
891 h_refs,
892 uint_refs,
893 sizes,
894 shambase::get_check_ref(storage.pressure),
895 shambase::get_check_ref(storage.soundspeed));
896 }
897}
898
899using namespace shammath;
900template class shammodels::sph::modules::ComputeEos<f64_3, M4>;
901template class shammodels::sph::modules::ComputeEos<f64_3, M6>;
902template class shammodels::sph::modules::ComputeEos<f64_3, M8>;
903
904template class shammodels::sph::modules::ComputeEos<f64_3, C2>;
905template class shammodels::sph::modules::ComputeEos<f64_3, C4>;
906template class shammodels::sph::modules::ComputeEos<f64_3, C6>;
shammodels::gsph::names::newtonian::uint
constexpr const char * uint
Specific internal energy u.
Definition FieldNames.hpp:49
shammodels::gsph::names::common::xyz
constexpr const char * xyz
Position field (3D coordinates).
Definition FieldNames.hpp:33
shammodels::gsph::edges::sizes
constexpr const char * sizes
Temporary sizes for h-iteration.
Definition FieldNames.hpp:113
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::sph::modules::ComputeEos
Module for computing equation of state quantities.
Definition ComputeEos.hpp:35
shammodels::sph::modules::ComputeEos::compute_eos
void compute_eos()
Computes pressure and sound speed from equation of state.
Definition ComputeEos.cpp:720
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::FieldRefs::get_spans
virtual DDPatchDataFieldSpanPointer< T > & get_spans()
Get the DistributedData of spans attached to the underlying field.
Definition FieldRefs.hpp:49
shamrock::solvergraph::Field::get_spans
virtual DDPatchDataFieldSpanPointer< T > & get_spans()
Get the DistributedData of spans attached to the underlying field.
Definition Field.hpp:58
shamrock::solvergraph::Field::ensure_sizes
virtual void ensure_sizes(const shambase::DistributedData< u32 > &sizes)
Ensure that the sizes of the patches in the field match the given sizes (Can resize the underlying fi...
Definition Field.hpp:92
shamrock::solvergraph::IFieldSpan
Interface for a solver graph edge representing a field as spans.
Definition IFieldSpan.hpp:37
shamrock::solvergraph::IFieldSpan::get_spans
virtual DDPatchDataFieldSpanPointer< T > & get_spans()=0
Get the DistributedData of spans attached to the underlying field.
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:110
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
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