Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
SPHUtilities.cpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
18#include "shammodels/sph/SPHUtilities.hpp"
19#include "shammath/sphkernels.hpp"
20#include "shamtree/TreeTraversal.hpp"
21
22using namespace shamrock::sph;
23
24namespace shammodels::sph {
25
26 template<class vec, class SPHKernel>
27 void SPHUtilities<vec, SPHKernel>::iterate_smoothing_length_cache(
28
29 sham::DeviceBuffer<vec> &merged_r,
30 sham::DeviceBuffer<flt> &hnew,
31 sham::DeviceBuffer<flt> &hold,
32 sham::DeviceBuffer<flt> &eps_h,
33 sycl::range<1> update_range,
34 shamrock::tree::ObjectCache &neigh_cache,
35
36 flt gpart_mass,
37 flt h_evol_max,
38 flt h_evol_iter_max
39
40 ) {
41
42 sham::DeviceQueue &queue = shamsys::instance::get_compute_scheduler().get_queue();
43
44 sham::EventList depends_list;
45
46 auto r = merged_r.get_read_access(depends_list);
47 auto h_old = hold.get_read_access(depends_list);
48 auto ploop_ptrs = neigh_cache.get_read_access(depends_list);
49 auto eps = eps_h.get_write_access(depends_list);
50 auto h_new = hnew.get_write_access(depends_list);
51
52 auto e = queue.submit(depends_list, [&](sycl::handler &cgh) {
53 // tree::ObjectIterator particle_looper(tree,cgh);
54
55 shamrock::tree::ObjectCacheIterator particle_looper(ploop_ptrs);
56
57 // sycl::accessor omega {omega_h, cgh, sycl::write_only, sycl::no_init};
58
59 const flt part_mass = gpart_mass;
60 const flt h_max_tot_max_evol = h_evol_max;
61 const flt h_max_evol_p = h_evol_iter_max;
62 const flt h_max_evol_m = 1 / h_evol_iter_max;
63
64 shambase::parallel_for(cgh, update_range.size(), "iter h", [=](u32 id_a) {
65 if (eps[id_a] > 1e-6) {
66
67 vec xyz_a = r[id_a]; // could be recovered from lambda
68
69 flt h_a = h_new[id_a];
70 flt dint = h_a * h_a * Rkern * Rkern;
71
72 vec inter_box_a_min = xyz_a - h_a * Rkern;
73 vec inter_box_a_max = xyz_a + h_a * Rkern;
74
75 flt rho_sum = 0;
76 flt sumdWdh = 0;
77
78 // particle_looper.rtree_for([&](u32, vec bmin,vec bmax) -> bool {
79 // return
80 // shammath::domain_are_connected(bmin,bmax,inter_box_a_min,inter_box_a_max);
81 // },[&](u32 id_b){
82 particle_looper.for_each_object(id_a, [&](u32 id_b) {
83 vec dr = xyz_a - r[id_b];
84 flt rab2 = sycl::dot(dr, dr);
85
86 if (rab2 > dint) {
87 return;
88 }
89
90 flt rab = sycl::sqrt(rab2);
91
92 rho_sum += part_mass * SPHKernel::W_3d(rab, h_a);
93 sumdWdh += part_mass * SPHKernel::dhW_3d(rab, h_a);
94 });
95
96 using namespace shamrock::sph;
97
98 flt rho_ha = rho_h(part_mass, h_a, SPHKernel::hfactd);
99 flt new_h = newton_iterate_new_h(rho_ha, rho_sum, sumdWdh, h_a);
100
101 if (new_h < h_a * h_max_evol_m)
102 new_h = h_max_evol_m * h_a;
103 if (new_h > h_a * h_max_evol_p)
104 new_h = h_max_evol_p * h_a;
105
106 flt ha_0 = h_old[id_a];
107
108 if (new_h < ha_0 * h_max_tot_max_evol) {
109 h_new[id_a] = new_h;
110 eps[id_a] = sycl::fabs(new_h - h_a) / ha_0;
111 } else {
112 h_new[id_a] = ha_0 * h_max_tot_max_evol;
113 eps[id_a] = -1;
114 }
115 }
116 });
117 });
118
119 merged_r.complete_event_state(e);
120 eps_h.complete_event_state(e);
121 hnew.complete_event_state(e);
122 hold.complete_event_state(e);
123
124 sham::EventList resulting_events;
125 resulting_events.add_event(e);
126 neigh_cache.complete_event_state(resulting_events);
127 }
128
129 template<class vec, class SPHKernel, class u_morton>
130 void SPHTreeUtilities<vec, SPHKernel, u_morton>::iterate_smoothing_length_tree(
131 sycl::buffer<vec> &merged_r,
132 sycl::buffer<flt> &hnew,
133 sycl::buffer<flt> &hold,
134 sycl::buffer<flt> &eps_h,
135 sycl::range<1> update_range,
136 RadixTree<u_morton, vec> &tree,
137
138 flt gpart_mass,
139 flt h_evol_max,
140 flt h_evol_iter_max) {
141
142 shamsys::instance::get_compute_queue()
143 .submit([&](sycl::handler &cgh) {
144 shamrock::tree::ObjectIterator particle_looper(tree, cgh);
145 // shamrock::tree::ObjectCacheIterator particle_looper(neigh_cache, cgh);
146
147 sycl::accessor eps{eps_h, cgh, sycl::read_write};
148 sycl::accessor r{merged_r, cgh, sycl::read_only};
149 sycl::accessor h_new{hnew, cgh, sycl::read_write};
150 sycl::accessor h_old{hold, cgh, sycl::read_only};
151 // sycl::accessor omega {omega_h, cgh, sycl::write_only, sycl::no_init};
152
153 const flt part_mass = gpart_mass;
154 const flt h_max_tot_max_evol = h_evol_max;
155 const flt h_max_evol_p = h_evol_iter_max;
156 const flt h_max_evol_m = 1 / h_evol_iter_max;
157
158 shambase::parallel_for(cgh, update_range.size(), "iter h", [=](u32 id_a) {
159 if (eps[id_a] > 1e-6) {
160
161 vec xyz_a = r[id_a]; // could be recovered from lambda
162
163 flt h_a = h_new[id_a];
164 flt dint = h_a * h_a * Rkern * Rkern;
165
166 vec inter_box_a_min = xyz_a - h_a * Rkern;
167 vec inter_box_a_max = xyz_a + h_a * Rkern;
168
169 flt rho_sum = 0;
170 flt sumdWdh = 0;
171
172 particle_looper.rtree_for(
173 [&](u32, vec bmin, vec bmax) -> bool {
174 return shammath::domain_are_connected(
175 bmin, bmax, inter_box_a_min, inter_box_a_max);
176 },
177 [&](u32 id_b) {
178 // particle_looper.for_each_object(id_a, [&](u32 id_b) {
179 vec dr = xyz_a - r[id_b];
180 flt rab2 = sycl::dot(dr, dr);
181
182 if (rab2 > dint) {
183 return;
184 }
185
186 flt rab = sycl::sqrt(rab2);
187
188 rho_sum += part_mass * SPHKernel::W_3d(rab, h_a);
189 sumdWdh += part_mass * SPHKernel::dhW_3d(rab, h_a);
190 });
191
192 using namespace shamrock::sph;
193
194 flt rho_ha = rho_h(part_mass, h_a, SPHKernel::hfactd);
195 flt new_h = newton_iterate_new_h(rho_ha, rho_sum, sumdWdh, h_a);
196
197 if (new_h < h_a * h_max_evol_m)
198 new_h = h_max_evol_m * h_a;
199 if (new_h > h_a * h_max_evol_p)
200 new_h = h_max_evol_p * h_a;
201
202 flt ha_0 = h_old[id_a];
203
204 if (new_h < ha_0 * h_max_tot_max_evol) {
205 h_new[id_a] = new_h;
206 eps[id_a] = sycl::fabs(new_h - h_a) / ha_0;
207 } else {
208 h_new[id_a] = ha_0 * h_max_tot_max_evol;
209 eps[id_a] = -1;
210 }
211 }
212 });
213 })
214 .wait();
215 }
216
217 template<class vec, class SPHKernel>
218 void SPHUtilities<vec, SPHKernel>::compute_omega(
219 sham::DeviceBuffer<vec> &merged_r,
220 sham::DeviceBuffer<flt> &h_part,
221 sham::DeviceBuffer<flt> &omega_h,
222 sycl::range<1> part_range,
223 shamrock::tree::ObjectCache &neigh_cache,
224 flt gpart_mass) {
225 using namespace shamrock::tree;
226
227 sham::DeviceQueue &queue = shamsys::instance::get_compute_scheduler().get_queue();
228
229 sham::EventList depends_list;
230
231 auto r = merged_r.get_read_access(depends_list);
232 auto hpart = h_part.get_read_access(depends_list);
233 auto omega = omega_h.get_write_access(depends_list);
234 auto ploop_ptrs = neigh_cache.get_read_access(depends_list);
235
236 auto e = queue.submit(depends_list, [&](sycl::handler &cgh) {
237 // tree::ObjectIterator particle_looper(tree,cgh);
238
239 ObjectCacheIterator particle_looper(ploop_ptrs);
240
241 const flt part_mass = gpart_mass;
242
243 shambase::parallel_for(cgh, part_range.size(), "compute omega", [=](u32 id_a) {
244 vec xyz_a = r[id_a]; // could be recovered from lambda
245
246 flt h_a = hpart[id_a];
247 flt dint = h_a * h_a * Rkern * Rkern;
248
249 // vec inter_box_a_min = xyz_a - h_a * Rkern;
250 // vec inter_box_a_max = xyz_a + h_a * Rkern;
251
252 flt rho_sum = 0;
253 flt part_omega_sum = 0;
254
255 // particle_looper.rtree_for([&](u32, vec bmin,vec bmax) -> bool {
256 // return
257 // shammath::domain_are_connected(bmin,bmax,inter_box_a_min,inter_box_a_max);
258 // },[&](u32 id_b){
259 particle_looper.for_each_object(id_a, [&](u32 id_b) {
260 vec dr = xyz_a - r[id_b];
261 flt rab2 = sycl::dot(dr, dr);
262
263 if (rab2 > dint) {
264 return;
265 }
266
267 flt rab = sycl::sqrt(rab2);
268
269 rho_sum += part_mass * SPHKernel::W_3d(rab, h_a);
270 part_omega_sum += part_mass * SPHKernel::dhW_3d(rab, h_a);
271 });
272
273 using namespace shamrock::sph;
274
275 flt rho_ha = rho_h(part_mass, h_a, SPHKernel::hfactd);
276 flt omega_a = 1 + (h_a / (3 * rho_ha)) * part_omega_sum;
277 omega[id_a] = omega_a;
278
279 // logger::raw(shambase::format("pmass {}, rho_a {}, omega_a {}\n",
280 // part_mass,rho_ha, omega_a));
281 });
282 });
283
284 merged_r.complete_event_state(e);
285 h_part.complete_event_state(e);
286 omega_h.complete_event_state(e);
287
288 sham::EventList resulting_events;
289 resulting_events.add_event(e);
290 neigh_cache.complete_event_state(resulting_events);
291 }
292
293 template class SPHUtilities<f64_3, shammath::M4<f64>>;
294 template class SPHUtilities<f64_3, shammath::M6<f64>>;
295 template class SPHUtilities<f64_3, shammath::M8<f64>>;
296
297 template class SPHUtilities<f64_3, shammath::C2<f64>>;
298 template class SPHUtilities<f64_3, shammath::C4<f64>>;
299 template class SPHUtilities<f64_3, shammath::C6<f64>>;
300
301 template class SPHTreeUtilities<f64_3, shammath::M4<f64>, u32>;
302 template class SPHTreeUtilities<f64_3, shammath::M6<f64>, u64>;
303 template class SPHTreeUtilities<f64_3, shammath::M8<f64>, u64>;
304
305 template class SPHTreeUtilities<f64_3, shammath::C2<f64>, u32>;
306 template class SPHTreeUtilities<f64_3, shammath::C4<f64>, u64>;
307 template class SPHTreeUtilities<f64_3, shammath::C6<f64>, u64>;
308
309} // namespace shammodels::sph
shammodels::gsph::edges::eps_h
constexpr const char * eps_h
Epsilon h references (for h-iteration convergence)
Definition FieldNames.hpp:131
shammodels::gsph::edges::h_new
constexpr const char * h_new
New smoothing length references (for h-iteration)
Definition FieldNames.hpp:128
shammodels::gsph::edges::h_old
constexpr const char * h_old
Old smoothing length references (for h-iteration)
Definition FieldNames.hpp:125
shammodels::gsph::names::common::hpart
constexpr const char * hpart
Smoothing length field.
Definition FieldNames.hpp:34
shammodels::gsph::edges::neigh_cache
constexpr const char * neigh_cache
Neighbor cache.
Definition FieldNames.hpp:110
shammodels::gsph::names::newtonian::omega
constexpr const char * omega
Grad-h correction factor \Omega.
Definition FieldNames.hpp:64
shamsys::instance::get_compute_queue
sycl::queue & get_compute_queue(u32 id=0)
Definition NodeInstance.cpp:370
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
RadixTree
The radix tree.
Definition RadixTree.hpp:50
sham::DeviceBuffer
A buffer allocated in USM (Unified Shared Memory)
Definition DeviceBuffer.hpp:106
sham::DeviceBuffer::complete_event_state
void complete_event_state(sycl::event e) const
Complete the event state of the buffer.
Definition DeviceBuffer.hpp:368
sham::DeviceBuffer::get_write_access
T * get_write_access(sham::EventList &depends_list, SourceLocation src_loc=SourceLocation{})
Get a read-write pointer to the buffer's data.
Definition DeviceBuffer.hpp:349
sham::DeviceBuffer::get_read_access
const T * get_read_access(sham::EventList &depends_list, SourceLocation src_loc=SourceLocation{}) const
Get a read-only pointer to the buffer's data.
Definition DeviceBuffer.hpp:331
sham::DeviceQueue
A SYCL queue associated with a device and a context.
Definition DeviceQueue.hpp:32
sham::DeviceQueue::submit
sycl::event submit(Fct &&fct)
Submits a kernel to the SYCL queue.
Definition DeviceQueue.hpp:101
sham::DeviceScheduler::get_queue
DeviceQueue & get_queue(u32 id=0)
Get a reference to a DeviceQueue.
Definition DeviceScheduler.cpp:57
sham::EventList
Class to manage a list of SYCL events.
Definition EventList.hpp:31
sham::EventList::add_event
void add_event(sycl::event e)
Add an event to the list of events.
Definition EventList.hpp:87
shammodels::sph
namespace for the sph model
Definition BasicSPHGhosts.hpp:37
sphkernels.hpp
sph kernels
Generated on Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8