41 Tscal operator()(
u32 i)
const {
42 using namespace shamrock::sph;
43 return rho_h(pmass, h[i], hfact);
65 const Tscal *buf_epsilon;
70 Tscal operator()(
u32 i)
const {
72 Tscal epsilon_sum = 0;
73 for (
u32 j = 0; j < nvar_dust; j++) {
74 epsilon_sum += buf_epsilon[i * nvar_dust + j];
77 using namespace shamrock::sph;
78 return (1 - epsilon_sum) * rho_h(pmass, h[i], hfact);
86 return accessed{h, epsilon, nvar_dust, pmass, hfact};
92template<
class Tvec,
template<
class>
class SPHKernel>
93template<
class RhoGetGen>
95 RhoGetGen &&rho_getter_gen) {
102 using namespace shamrock::patch;
104 using SolverConfigEOS =
typename Config::EOSConfig;
105 using SolverEOS_Isothermal =
typename SolverConfigEOS::Isothermal;
106 using SolverEOS_Adiabatic =
typename SolverConfigEOS::Adiabatic;
107 using SolverEOS_Polytropic =
typename SolverConfigEOS::Polytropic;
108 using SolverEOS_LocallyIsothermal =
typename SolverConfigEOS::LocallyIsothermal;
109 using SolverEOS_LocallyIsothermalLP07 =
typename SolverConfigEOS::LocallyIsothermalLP07;
110 using SolverEOS_LocallyIsothermalFA2014 =
typename SolverConfigEOS::LocallyIsothermalFA2014;
111 using SolverEOS_LocallyIsothermalFA2014Extended =
112 typename SolverConfigEOS::LocallyIsothermalFA2014Extended;
113 using SolverEOS_Fermi =
typename SolverConfigEOS::Fermi;
117 if (SolverEOS_Isothermal *eos_config
118 = std::get_if<SolverEOS_Isothermal>(&solver_config.eos_config.config)) {
122 storage.merged_patchdata_ghost.get().for_each([&](
u64 id,
PatchDataLayer &mpdat) {
127 auto rho_getter = rho_getter_gen(mpdat);
138 = eos_config->cs](
u32 i,
auto rho, Tscal *__restrict P, Tscal *__restrict cs) {
139 using namespace shamrock::sph;
140 Tscal rho_a = rho(i);
141 Tscal P_a = EOS::pressure(cs_cfg, rho_a);
147 SolverEOS_Adiabatic *eos_config
148 = std::get_if<SolverEOS_Adiabatic>(&solver_config.eos_config.config)) {
152 storage.merged_patchdata_ghost.get().for_each([&](
u64 id,
PatchDataLayer &mpdat) {
158 auto rho_getter = rho_getter_gen(mpdat);
168 [gamma = eos_config->gamma](
171 const Tscal *__restrict U,
173 Tscal *__restrict cs) {
174 using namespace shamrock::sph;
175 Tscal rho_a = rho(i);
176 Tscal P_a = EOS::pressure(gamma, rho_a, U[i]);
177 Tscal cs_a = EOS::cs_from_p(gamma, rho_a, P_a);
184 SolverEOS_Polytropic *eos_config
185 = std::get_if<SolverEOS_Polytropic>(&solver_config.eos_config.config)) {
189 storage.merged_patchdata_ghost.get().for_each([&](
u64 id,
PatchDataLayer &mpdat) {
194 auto rho_getter = rho_getter_gen(mpdat);
204 [K = eos_config->K, gamma = eos_config->gamma](
205 u32 i,
auto rho, Tscal *__restrict P, Tscal *__restrict cs) {
206 using namespace shamrock::sph;
207 Tscal rho_a = rho(i);
208 Tscal P_a = EOS::pressure(gamma, K, rho_a);
209 Tscal cs_a = EOS::soundspeed(gamma, K, rho_a);
216 SolverEOS_LocallyIsothermal *eos_config
217 = std::get_if<SolverEOS_LocallyIsothermal>(&solver_config.eos_config.config)) {
223 storage.merged_patchdata_ghost.get().for_each([&](
u64 id,
PatchDataLayer &mpdat) {
229 auto rho_getter = rho_getter_gen(mpdat);
242 const Tscal *__restrict U,
243 const Tscal *__restrict cs0,
245 Tscal *__restrict cs) {
246 using namespace shamrock::sph;
248 Tscal cs_out = cs0[i];
249 Tscal rho_a = rho(i);
251 Tscal P_a = EOS::pressure_from_cs(cs_out * cs_out, rho_a);
259 SolverEOS_LocallyIsothermalLP07 *eos_config
260 = std::get_if<SolverEOS_LocallyIsothermalLP07>(&solver_config.eos_config.config)) {
264 storage.merged_patchdata_ghost.get().for_each([&](
u64 id,
PatchDataLayer &mpdat) {
265 auto &mfield = storage.merged_xyzh.get().get(
id);
274 auto rho_getter = rho_getter_gen(mpdat);
276 Tscal cs0 = eos_config->cs0;
277 Tscal r0sq = eos_config->r0 * eos_config->r0;
278 Tscal mq = -eos_config->q;
291 const Tscal *__restrict U,
292 const Tvec *__restrict
xyz,
294 Tscal *__restrict cs) {
295 using namespace shamrock::sph;
298 Tscal rho_a = rho(i);
300 Tscal Rsq = sycl::dot(R, R);
301 Tscal cs_sq = EOS::soundspeed_sq(cs0 * cs0, Rsq / r0sq, mq);
302 Tscal cs_out = sycl::sqrt(cs_sq);
304 Tscal P_a = EOS::pressure_from_cs(cs_sq, rho_a);
312 SolverEOS_LocallyIsothermalFA2014 *eos_config
313 = std::get_if<SolverEOS_LocallyIsothermalFA2014>(&solver_config.eos_config.config)) {
315 Tscal _G = solver_config.get_constant_G();
319 auto &sink_parts = storage.sinks.get();
320 std::vector<Tvec> sink_pos;
321 std::vector<Tscal> sink_mass;
324 for (
auto &s : sink_parts) {
325 sink_pos.push_back(s.pos);
326 sink_mass.push_back(s.mass);
330 sycl::buffer<Tvec> sink_pos_buf{sink_pos};
331 sycl::buffer<Tscal> sink_mass_buf{sink_mass};
333 storage.merged_patchdata_ghost.get().for_each([&](
u64 id,
PatchDataLayer &mpdat) {
334 auto &mfield = storage.merged_xyzh.get().get(
id);
343 auto rho_getter = rho_getter_gen(mpdat);
351 auto rho = rho_getter.get_read_access(depends_list);
358 auto e = q.
submit(depends_list, [&](sycl::handler &cgh) {
359 sycl::accessor spos{sink_pos_buf, cgh, sycl::read_only};
360 sycl::accessor smass{sink_mass_buf, cgh, sycl::read_only};
361 u32 scount = sink_cnt;
363 Tscal h_over_r = eos_config->h_over_r;
366 cgh.parallel_for(sycl::range<1>{total_elements}, [=](sycl::item<1> item) {
367 using namespace shamrock::sph;
370 Tscal rho_a = rho(item.get_linear_id());
372 Tscal mpotential = 0;
373 for (
u32 i = 0; i < scount; i++) {
374 Tvec s_r = spos[i] - R;
375 Tscal s_m = smass[i];
376 Tscal s_r_abs = sycl::length(s_r);
377 mpotential += G * s_m / s_r_abs;
380 Tscal cs_out = h_over_r * sycl::sqrt(mpotential);
381 Tscal P_a = EOS::pressure_from_cs(cs_out * cs_out, rho_a);
390 rho_getter.complete_event_state(e);
396 SolverEOS_LocallyIsothermalFA2014Extended *eos_config
397 = std::get_if<SolverEOS_LocallyIsothermalFA2014Extended>(
398 &solver_config.eos_config.config)) {
400 Tscal _cs0 = eos_config->cs0;
401 Tscal _r0 = eos_config->r0;
402 Tscal _q = eos_config->q;
403 u32 n_sinks = eos_config->n_sinks;
407 auto &sink_parts = storage.sinks.get();
408 std::vector<Tvec> sink_pos;
409 std::vector<Tscal> sink_mass;
412 for (
auto &s : sink_parts) {
413 sink_pos.push_back(s.pos);
414 sink_mass.push_back(s.mass);
416 if (sink_pos.size() >= n_sinks) {
423 "No sinks found for the equation of state");
426 sycl::buffer<Tvec> sink_pos_buf{sink_pos};
427 sycl::buffer<Tscal> sink_mass_buf{sink_mass};
429 storage.merged_patchdata_ghost.get().for_each([&](
u64 id,
PatchDataLayer &mpdat) {
430 auto &mfield = storage.merged_xyzh.get().get(
id);
439 auto rho_getter = rho_getter_gen(mpdat);
447 auto rho = rho_getter.get_read_access(depends_list);
454 auto e = q.
submit(depends_list, [&](sycl::handler &cgh) {
455 sycl::accessor spos{sink_pos_buf, cgh, sycl::read_only};
456 sycl::accessor smass{sink_mass_buf, cgh, sycl::read_only};
457 u32 scount = sink_cnt;
463 Tscal inv_r0_q = 1. / sycl::pow(r0, q);
465 cgh.parallel_for(sycl::range<1>{total_elements}, [=](sycl::item<1> item) {
466 using namespace shamrock::sph;
469 Tscal rho_a = rho(item.get_linear_id());
471 Tscal sink_mass_sum = 0;
473 for (
u32 i = 0; i < scount; i++) {
474 Tvec s_r = spos[i] - R;
475 Tscal s_m = smass[i];
476 Tscal s_r_abs = sycl::length(s_r);
477 sink_mass_sum += s_m;
478 pot_sum += s_m / s_r_abs;
481 Tscal cs_out = cs0 * inv_r0_q * sycl::pow(pot_sum / sink_mass_sum, q);
482 Tscal P_a = EOS::pressure_from_cs(cs_out * cs_out, rho_a);
491 rho_getter.complete_event_state(e);
497 SolverEOS_Fermi *eos_config
498 = std::get_if<SolverEOS_Fermi>(&solver_config.eos_config.config)) {
502 storage.merged_patchdata_ghost.get().for_each([&](
u64 id,
PatchDataLayer &mpdat) {
507 auto rho_getter = rho_getter_gen(mpdat);
513 auto unit_sys = *solver_config.unit_sys;
515 Tscal mass = unit_sys.template to<units::kilogram>();
516 Tscal length = unit_sys.template to<units::metre>();
517 Tscal time = unit_sys.template to<units::second>();
519 Tscal pressure_unit = mass / length / (time * time);
520 Tscal density_unit = mass / (length * length * length);
521 Tscal velocity_unit = length / time;
528 [mu_e = eos_config->mu_e, density_unit, pressure_unit, velocity_unit](
529 u32 i,
auto rho, Tscal *__restrict P, Tscal *__restrict cs) {
530 Tscal rho_a = rho(i);
531 auto const res = EOS::pressure_and_soundspeed(mu_e, rho_a * density_unit);
532 P[i] = res.pressure / pressure_unit;
533 cs[i] = res.soundspeed / velocity_unit;
542template<
class Tvec,
template<
class>
class SPHKernel>
547 Tscal gpart_mass = solver_config.gpart_mass;
550 using namespace shamrock::patch;
564 if (solver_config.dust_config.has_epsilon_field()) {
567 u32 nvar_dust = solver_config.dust_config.get_dust_nvar();
571 mpdat.get_field_buf_ref<Tscal>(ihpart_interf),
572 mpdat.get_field_buf_ref<Tscal>(iepsilon_interf),
580 mpdat.get_field_buf_ref<Tscal>(ihpart_interf), gpart_mass, Kernel::hfactd};
constexpr const char * xyz
Position field (3D coordinates)
std::uint32_t u32
32 bit unsigned integer
std::uint64_t u64
64 bit unsigned integer
A buffer allocated in USM (Unified Shared Memory)
void complete_event_state(sycl::event e) const
Complete the event state of the buffer.
T * get_write_access(sham::EventList &depends_list, SourceLocation src_loc=SourceLocation{})
Get a read-write pointer to the buffer's data.
const T * get_read_access(sham::EventList &depends_list, SourceLocation src_loc=SourceLocation{}) const
Get a read-only pointer to the buffer's data.
A SYCL queue associated with a device and a context.
sycl::event submit(Fct &&fct)
Submits a kernel to the SYCL queue.
DeviceQueue & get_queue(u32 id=0)
Get a reference to a DeviceQueue.
Class to manage a list of SYCL events.
Represents a collection of objects distributed across patches identified by a u64 id.
Module for computing equation of state quantities.
void compute_eos()
Computes pressure and sound speed from equation of state.
u32 get_field_idx(const std::string &field_name) const
Get the field id if matching name & type.
PatchDataLayer container class, the layout is described in patchdata_layout.
This header file contains utility functions related to exception handling in the code.
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.
void throw_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Throw an exception and append the source location to it.
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...
void throw_unimplemented(SourceLocation loc=SourceLocation{})
Throw a std::runtime_error saying that the function is unimplemented.
namespace for math utility
namespace for the main framework
namespace containing the units library
A class that references multiple buffers or similar objects.
Adiabatic equation of state.
Isothermal equation of state.
Locally isothermal equation of state with radial dependence.
Polytropic equation of state.