Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
integrators.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
16
17#include "shambase/exception.hpp"
18#include "shambackends/math.hpp"
19#include "shambackends/sycl_utils.hpp"
20#include "shamrock/math/integrators.hpp"
21#include <algorithm>
22
23namespace integ = shamrock::integrators;
24namespace util = shamrock::utilities;
25
29
30template<class flt, class T>
31void integ::forward_euler(
32 sham::DeviceQueue &queue,
33 sham::DeviceBuffer<T> &buf_val,
34 sham::DeviceBuffer<T> &buf_der,
35 sycl::range<1> elem_range,
36 flt dt) {
37
38 sham::EventList depends_list;
39
40 auto acc_u = buf_val.get_write_access(depends_list);
41 auto acc_du = buf_der.get_read_access(depends_list);
42
43 auto e = queue.submit(depends_list, [&](sycl::handler &cgh) {
44 cgh.parallel_for(elem_range, [=](sycl::item<1> item) {
45 u32 gid = (u32) item.get_id();
46 acc_u[item] = acc_u[item] + (dt) *acc_du[item];
47 });
48 });
49
50 buf_val.complete_event_state(e);
51 buf_der.complete_event_state(e);
52}
53
54#ifndef DOXYGEN
55template void integ::forward_euler(
56 sham::DeviceQueue &queue,
57 sham::DeviceBuffer<f32_3> &buf_val,
58 sham::DeviceBuffer<f32_3> &buf_der,
59 sycl::range<1> elem_range,
60 f32 dt);
61
62template void integ::forward_euler(
63 sham::DeviceQueue &queue,
64 sham::DeviceBuffer<f32> &buf_val,
65 sham::DeviceBuffer<f32> &buf_der,
66 sycl::range<1> elem_range,
67 f32 dt);
68
69template void integ::forward_euler(
70 sham::DeviceQueue &queue,
71 sham::DeviceBuffer<f64_3> &buf_val,
72 sham::DeviceBuffer<f64_3> &buf_der,
73 sycl::range<1> elem_range,
74 f64 dt);
75
76template void integ::forward_euler(
77 sham::DeviceQueue &queue,
78 sham::DeviceBuffer<f64> &buf_val,
79 sham::DeviceBuffer<f64> &buf_der,
80 sycl::range<1> elem_range,
81 f64 dt);
82#endif
86
87template<class flt, class T>
88void integ::leapfrog_corrector(
89 sham::DeviceQueue &queue,
90 sham::DeviceBuffer<T> &buf_val,
91 sham::DeviceBuffer<T> &buf_der,
92 sham::DeviceBuffer<T> &buf_der_old,
93 sham::DeviceBuffer<flt> &buf_eps_sq,
94 sycl::range<1> elem_range,
95 flt hdt) {
96
97 sham::EventList depends_list;
98
99 auto acc_u = buf_val.get_write_access(depends_list);
100 auto acc_du = buf_der.get_read_access(depends_list);
101 auto acc_du_old = buf_der_old.get_read_access(depends_list);
102 auto acc_epsilon_sq = buf_eps_sq.get_write_access(depends_list);
103
104 auto e = queue.submit(depends_list, [&](sycl::handler &cgh) {
105 cgh.parallel_for(elem_range, [=](sycl::item<1> item) {
106 u32 gid = (u32) item.get_id();
107
108 T incr = (hdt) * (acc_du[item] - acc_du_old[item]);
109
110 acc_u[item] = acc_u[item] + incr;
111 acc_epsilon_sq[item] = sycl::dot(incr, incr);
112 });
113 });
114
115 buf_val.complete_event_state(e);
116 buf_der.complete_event_state(e);
117 buf_der_old.complete_event_state(e);
118 buf_eps_sq.complete_event_state(e);
119}
120
121template<class flt, class T>
122void integ::leapfrog_corrector_positive_only(
123 sham::DeviceQueue &queue,
124 sham::DeviceBuffer<T> &buf_val,
125 sham::DeviceBuffer<T> &buf_der,
126 sham::DeviceBuffer<T> &buf_der_old,
127 sham::DeviceBuffer<flt> &buf_eps_sq,
128 sycl::range<1> elem_range,
129 flt hdt) {
130
131 sham::EventList depends_list;
132
133 auto acc_u = buf_val.get_write_access(depends_list);
134 auto acc_du = buf_der.get_read_access(depends_list);
135 auto acc_du_old = buf_der_old.get_read_access(depends_list);
136 auto acc_epsilon_sq = buf_eps_sq.get_write_access(depends_list);
137
138 auto e = queue.submit(depends_list, [&](sycl::handler &cgh) {
139 cgh.parallel_for(elem_range, [=](sycl::item<1> item) {
140 u32 gid = (u32) item.get_id();
141
142 T incr = (hdt) * (acc_du[item] - acc_du_old[item]);
143
144 acc_u[item] = sham::max(acc_u[item] + incr, 0.0);
145 acc_epsilon_sq[item] = sycl::dot(incr, incr);
146 });
147 });
148
149 buf_val.complete_event_state(e);
150 buf_der.complete_event_state(e);
151 buf_der_old.complete_event_state(e);
152 buf_eps_sq.complete_event_state(e);
153}
154
155#ifndef DOXYGEN
156template void integ::leapfrog_corrector(
157 sham::DeviceQueue &queue,
158 sham::DeviceBuffer<f32_3> &buf_val,
159 sham::DeviceBuffer<f32_3> &buf_der,
160 sham::DeviceBuffer<f32_3> &buf_der_old,
161 sham::DeviceBuffer<f32> &buf_eps_sq,
162 sycl::range<1> elem_range,
163 f32 hdt);
164
165template void integ::leapfrog_corrector(
166 sham::DeviceQueue &queue,
167 sham::DeviceBuffer<f32> &buf_val,
168 sham::DeviceBuffer<f32> &buf_der,
169 sham::DeviceBuffer<f32> &buf_der_old,
170 sham::DeviceBuffer<f32> &buf_eps_sq,
171 sycl::range<1> elem_range,
172 f32 hdt);
173
174template void integ::leapfrog_corrector(
175 sham::DeviceQueue &queue,
176 sham::DeviceBuffer<f64_3> &buf_val,
177 sham::DeviceBuffer<f64_3> &buf_der,
178 sham::DeviceBuffer<f64_3> &buf_der_old,
179 sham::DeviceBuffer<f64> &buf_eps_sq,
180 sycl::range<1> elem_range,
181 f64 hdt);
182
183template void integ::leapfrog_corrector(
184 sham::DeviceQueue &queue,
185 sham::DeviceBuffer<f64> &buf_val,
186 sham::DeviceBuffer<f64> &buf_der,
187 sham::DeviceBuffer<f64> &buf_der_old,
188 sham::DeviceBuffer<f64> &buf_eps_sq,
189 sycl::range<1> elem_range,
190 f64 hdt);
191
192template void integ::leapfrog_corrector_positive_only(
193 sham::DeviceQueue &queue,
194 sham::DeviceBuffer<f64> &buf_val,
195 sham::DeviceBuffer<f64> &buf_der,
196 sham::DeviceBuffer<f64> &buf_der_old,
197 sham::DeviceBuffer<f64> &buf_eps_sq,
198 sycl::range<1> elem_range,
199 f64 hdt);
200
201#endif
205
206template<class T>
207void util::sycl_position_modulo(
208 sham::DeviceQueue &queue,
209 sham::DeviceBuffer<T> &buf_xyz,
210 sycl::range<1> elem_range,
211 std::pair<T, T> box) {
212
213 sham::EventList depends_list;
214 auto xyz = buf_xyz.get_write_access(depends_list);
215
216 auto e = queue.submit(depends_list, [&](sycl::handler &cgh) {
217 T box_min = std::get<0>(box);
218 T box_max = std::get<1>(box);
219 T delt = box_max - box_min;
220
221 cgh.parallel_for(elem_range, [=](sycl::item<1> item) {
222 u32 gid = (u32) item.get_id();
223
224 T r = xyz[gid] - box_min;
225
226 r = sycl::fmod(r, delt);
227 r += delt;
228 r = sycl::fmod(r, delt);
229 r += box_min;
230
231 xyz[gid] = r;
232 });
233 });
234
235 buf_xyz.complete_event_state(e);
236}
237
238#ifndef DOXYGEN
239template void util::sycl_position_modulo(
240 sham::DeviceQueue &queue,
241 sham::DeviceBuffer<f32_3> &buf_xyz,
242 sycl::range<1> elem_range,
243 std::pair<f32_3, f32_3> box);
244
245template void util::sycl_position_modulo(
246 sham::DeviceQueue &queue,
247 sham::DeviceBuffer<f64_3> &buf_xyz,
248 sycl::range<1> elem_range,
249 std::pair<f64_3, f64_3> box);
250#endif
254
255template<class T>
256void util::sycl_position_sheared_modulo(
257 sham::DeviceQueue &queue,
258 sham::DeviceBuffer<T> &buf_xyz,
259 sham::DeviceBuffer<T> &buf_vxyz,
260 sycl::range<1> elem_range,
261 std::pair<T, T> box,
262 i32_3 shear_base,
263 i32_3 shear_dir,
264 shambase::VecComponent<T> shear_value,
265 shambase::VecComponent<T> shear_speed) {
266
267 sham::EventList depends_list;
268 auto xyz = buf_xyz.get_write_access(depends_list);
269 auto vxyz = buf_vxyz.get_write_access(depends_list);
270
271 auto e
272 = queue.submit(depends_list, [&, shear_base, shear_value, shear_speed](sycl::handler &cgh) {
273 T box_min = std::get<0>(box);
274 T box_max = std::get<1>(box);
275 T delt = box_max - box_min;
276
277 cgh.parallel_for(elem_range, [=](sycl::item<1> item) {
278 u32 gid = (u32) item.get_id();
279
280 T r = xyz[gid] - box_min;
281
282 T roff = r / delt;
283
284 // T dn = sycl::trunc(roff);
285
286 // auto d = sycl::dot(dn,shear_base.convert<shambase::VecComponent<T>>());
287
288 //*
289 auto cnt_per = [](shambase::VecComponent<T> v) -> int {
290 return (v >= 0) ? int(v) : (int(v) - 1);
291 };
292
293 i32 xoff = cnt_per(roff.x());
294 i32 yoff = cnt_per(roff.y());
295 i32 zoff = cnt_per(roff.z());
296
297 i32 dx = xoff * shear_base.x();
298 i32 dy = yoff * shear_base.y();
299 i32 dz = zoff * shear_base.z();
300
301 i32 d = dx + dy + dz;
302 //*/
303
304 T shift
305 = {(d * shear_dir.x()) * shear_value,
306 (d * shear_dir.y()) * shear_value,
307 (d * shear_dir.z()) * shear_value};
308
309 T shift_speed
310 = {(d * shear_dir.x()) * shear_speed,
311 (d * shear_dir.y()) * shear_speed,
312 (d * shear_dir.z()) * shear_speed};
313
314 vxyz[gid] -= shift_speed;
315 r -= shift;
316
317 r = sycl::fmod(r, delt);
318 r += delt;
319 r = sycl::fmod(r, delt);
320 r += box_min;
321
322 xyz[gid] = r;
323 });
324 });
325
326 buf_xyz.complete_event_state(e);
327 buf_vxyz.complete_event_state(e);
328}
329
330#ifndef DOXYGEN
331template void util::sycl_position_sheared_modulo(
332 sham::DeviceQueue &queue,
333 sham::DeviceBuffer<f32_3> &buf_xyz,
334 sham::DeviceBuffer<f32_3> &buf_vxyz,
335 sycl::range<1> elem_range,
336 std::pair<f32_3, f32_3> box,
337 i32_3 shear_base,
338 i32_3 shear_dir,
339 f32 shear_value,
340 f32 shear_speed);
341
342template void util::sycl_position_sheared_modulo(
343 sham::DeviceQueue &queue,
344 sham::DeviceBuffer<f64_3> &buf_xyz,
345 sham::DeviceBuffer<f64_3> &buf_vxyz,
346 sycl::range<1> elem_range,
347 std::pair<f64_3, f64_3> box,
348 i32_3 shear_base,
349 i32_3 shear_dir,
350 f64 shear_value,
351 f64 shear_speed);
352#endif
356
357template<class T>
358void util::swap_fields(
359 sham::DeviceQueue &queue, sham::DeviceBuffer<T> &b1, sham::DeviceBuffer<T> &b2, u32 cnt) {
360
361 sham::EventList depends_list;
362 auto acc1 = b1.get_write_access(depends_list);
363 auto acc2 = b2.get_write_access(depends_list);
364
365 auto e = queue.submit(depends_list, [&](sycl::handler &cgh) {
366 cgh.parallel_for(sycl::range<1>{cnt}, [=](sycl::item<1> item) {
367 T v1 = acc1[item];
368 T v2 = acc2[item];
369
370 acc1[item] = v2;
371 acc2[item] = v1;
372 });
373 });
374
375 b1.complete_event_state(e);
376 b2.complete_event_state(e);
377}
378
379#ifndef DOXYGEN
380template void util::swap_fields(
381 sham::DeviceQueue &queue,
382 sham::DeviceBuffer<f32_3> &b1,
383 sham::DeviceBuffer<f32_3> &b2,
384 u32 cnt);
385
386template void util::swap_fields(
387 sham::DeviceQueue &queue, sham::DeviceBuffer<f32> &b1, sham::DeviceBuffer<f32> &b2, u32 cnt);
388
389template void util::swap_fields(
390 sham::DeviceQueue &queue,
391 sham::DeviceBuffer<f64_3> &b1,
392 sham::DeviceBuffer<f64_3> &b2,
393 u32 cnt);
394
395template void util::swap_fields(
396 sham::DeviceQueue &queue, sham::DeviceBuffer<f64> &b1, sham::DeviceBuffer<f64> &b2, u32 cnt);
397#endif
f64
double f64
Alias for double.
Definition aliases_float.hpp:20
f32
float f32
Alias for float.
Definition aliases_float.hpp:19
u32
std::uint32_t u32
32 bit unsigned integer
Definition aliases_int.hpp:27
i32
std::int32_t i32
32 bit integer
Definition aliases_int.hpp:23
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::EventList
Class to manage a list of SYCL events.
Definition EventList.hpp:31
exception.hpp
This header file contains utility functions related to exception handling in the code.
shamrock::integrators::forward_euler
void forward_euler(sham::DeviceQueue &queue, sham::DeviceBuffer< T > &buf_val, sham::DeviceBuffer< T > &buf_der, sycl::range< 1 > elem_range, flt dt)
Perform forward Euler integration step.
Definition integrators.cpp:31
shamrock::utilities::sycl_position_sheared_modulo
void sycl_position_sheared_modulo(sham::DeviceQueue &queue, sham::DeviceBuffer< T > &buf_xyz, sham::DeviceBuffer< T > &buf_vxyz, sycl::range< 1 > elem_range, std::pair< T, T > box, i32_3 shear_base, i32_3 shear_dir, shambase::VecComponent< T > shear_value, shambase::VecComponent< T > shear_speed)
Apply periodic boundary conditions with shearing.
Definition integrators.cpp:256
shamrock::utilities::sycl_position_modulo
void sycl_position_modulo(sham::DeviceQueue &queue, sham::DeviceBuffer< T > &buf_xyz, sycl::range< 1 > elem_range, std::pair< T, T > box)
Apply periodic boundary conditions to positions.
Definition integrators.cpp:207
shamrock::integrators::leapfrog_corrector
void leapfrog_corrector(sham::DeviceQueue &queue, sham::DeviceBuffer< T > &buf_val, sham::DeviceBuffer< T > &buf_der, sham::DeviceBuffer< T > &buf_der_old, sham::DeviceBuffer< flt > &buf_eps_sq, sycl::range< 1 > elem_range, flt hdt)
Perform leapfrog corrector step with adaptive softening.
Definition integrators.cpp:88
shamrock::utilities::swap_fields
void swap_fields(sham::DeviceQueue &queue, sham::DeviceBuffer< T > &b1, sham::DeviceBuffer< T > &b2, u32 cnt)
Swap contents of two buffers.
Definition integrators.cpp:358
sycl_utils.hpp
Generated on for Shamrock by doxygen 1.15.0