Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
InterpolateToFace.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
17#include "shambase/assert.hpp"
18#include "shammath/AABB.hpp"
22
23namespace {
24
25 template<class Tvec, class TgridVec, class AMRBlock>
26 class GetShift {
27
28 public:
29 using Tscal = shambase::VecComponent<Tvec>;
30 using Tgridscal = shambase::VecComponent<TgridVec>;
31
32 const Tvec *acc_aabb_block_lower;
33 const Tscal *acc_aabb_cell_size;
34
35 GetShift(const Tvec *aabb_block_lower, const Tscal *aabb_cell_size)
36 : acc_aabb_block_lower{aabb_block_lower}, acc_aabb_cell_size{aabb_cell_size} {}
37
38 shammath::AABB<Tvec> get_cell_aabb(u32 id) const {
39
40 const u32 cell_global_id = (u32) id;
41
42 const u32 block_id = cell_global_id / AMRBlock::block_size;
43 const u32 cell_loc_id = cell_global_id % AMRBlock::block_size;
44
45 // fetch current block info
46 const Tvec cblock_min = acc_aabb_block_lower[block_id];
47 const Tscal delta_cell = acc_aabb_cell_size[block_id];
48
49 std::array<u32, 3> lcoord_arr = AMRBlock::get_coord(cell_loc_id);
50 Tvec offset = Tvec{lcoord_arr[0], lcoord_arr[1], lcoord_arr[2]} * delta_cell;
51
52 Tvec aabb_min = cblock_min + offset;
53 Tvec aabb_max = aabb_min + delta_cell;
54
55 return {aabb_min, aabb_max};
56 }
57
58 std::pair<Tvec, Tvec> get_shifts(u32 id_a, u32 id_b) const {
59
60 shammath::AABB<Tvec> aabb_cell_a = get_cell_aabb(id_a);
61 shammath::AABB<Tvec> aabb_cell_b = get_cell_aabb(id_b);
62
63 shammath::AABB<Tvec> face_aabb = aabb_cell_a.get_intersect(aabb_cell_b);
64
65 Tvec face_center = face_aabb.get_center();
66
67 Tvec shift_a = face_center - aabb_cell_a.get_center();
68 Tvec shift_b = face_center - aabb_cell_b.get_center();
69
70 return {shift_a, shift_b};
71 }
72 };
73
74 template<class Tvec, class TgridVec, class AMRBlock>
75 class RhoInterpolate {
76 using Tscal = shambase::VecComponent<Tvec>;
77
78 public:
83 // For time interpolation
84 Tscal dt_interp;
89
90 class acc {
91 public:
92 GetShift<Tvec, TgridVec, AMRBlock> shift_get;
93
94 const Tscal *acc_rho_cell;
95 const Tvec *acc_grad_rho_cell;
96
97 // For time interpolation
98 const Tvec *acc_vel_cell;
99 const Tvec *acc_dx_v_cell;
100 const Tvec *acc_dy_v_cell;
101 const Tvec *acc_dz_v_cell;
102
103 Tscal dt_interp;
104
105 acc(const Tvec *aabb_block_lower,
106 const Tscal *aabb_cell_size,
107 const Tscal *rho_cell,
108 const Tvec *grad_rho_cell,
109 // For time interpolation
110 Tscal dt_interp,
111 const Tvec *vel_cell,
112 const Tvec *dx_v_cell,
113 const Tvec *dy_v_cell,
114 const Tvec *dz_v_cell)
115 : shift_get(aabb_block_lower, aabb_cell_size), acc_rho_cell{rho_cell},
116 acc_grad_rho_cell{grad_rho_cell}, dt_interp(dt_interp), acc_vel_cell{vel_cell},
117 acc_dx_v_cell{dx_v_cell}, acc_dy_v_cell{dy_v_cell}, acc_dz_v_cell{dz_v_cell} {}
118
119 Tscal get_dt_rho(
120 Tscal rho, Tvec v, Tvec grad_rho, Tvec dx_v, Tvec dy_v, Tvec dz_v) const {
121 return -(sham::dot(v, grad_rho) + rho * (dx_v[0] + dy_v[1] + dz_v[2]));
122 }
123
124 std::array<Tscal, 2> get_link_field_val(u32 id_a, u32 id_b) const {
125
126 auto [shift_a, shift_b] = shift_get.get_shifts(id_a, id_b);
127
128 Tscal rho_a = acc_rho_cell[id_a];
129 Tvec grad_rho_a = acc_grad_rho_cell[id_a];
130 Tscal rho_b = acc_rho_cell[id_b];
131 Tvec grad_rho_b = acc_grad_rho_cell[id_b];
132
133 Tvec vel_a = acc_vel_cell[id_a];
134 Tvec dx_v_a = acc_dx_v_cell[id_a];
135 Tvec dy_v_a = acc_dy_v_cell[id_a];
136 Tvec dz_v_a = acc_dz_v_cell[id_a];
137 Tvec vel_b = acc_vel_cell[id_b];
138 Tvec dx_v_b = acc_dx_v_cell[id_b];
139 Tvec dy_v_b = acc_dy_v_cell[id_b];
140 Tvec dz_v_b = acc_dz_v_cell[id_b];
141
142 // Spatial interpolate
143 Tscal rho_face_a = rho_a + sycl::dot(grad_rho_a, shift_a);
144 Tscal rho_face_b = rho_b + sycl::dot(grad_rho_b, shift_b);
145
146 // Interpolate also to half a timestep
147 rho_face_a
148 += get_dt_rho(rho_a, vel_a, grad_rho_a, dx_v_a, dy_v_a, dz_v_a) * dt_interp;
149 rho_face_b
150 += get_dt_rho(rho_b, vel_b, grad_rho_b, dx_v_b, dy_v_b, dz_v_b) * dt_interp;
151
152 return {rho_face_a, rho_face_b};
153 }
154 };
155
156 inline acc get_read_access(sham::EventList &deps) {
157 return acc(
158 aabb_block_lower.get_read_access(deps),
159 aabb_cell_size.get_read_access(deps),
160 rho_cell.get_read_access(deps),
161 grad_rho_cell.get_read_access(deps),
162 // For time interpolation
163 dt_interp,
164 vel_cell.get_read_access(deps),
165 dx_v_cell.get_read_access(deps),
166 dy_v_cell.get_read_access(deps),
167 dz_v_cell.get_read_access(deps));
168 }
169
170 inline void complete_event_state(sycl::event e) {
171 aabb_block_lower.complete_event_state(e);
172 aabb_cell_size.complete_event_state(e);
173 rho_cell.complete_event_state(e);
174 grad_rho_cell.complete_event_state(e);
175 vel_cell.complete_event_state(e);
176 dx_v_cell.complete_event_state(e);
177 dy_v_cell.complete_event_state(e);
178 dz_v_cell.complete_event_state(e);
179 }
180 };
181
182 template<class Tvec, class TgridVec, class AMRBlock>
183 class VelInterpolate {
184 using Tscal = shambase::VecComponent<Tvec>;
185
186 public:
189
194 // For time interpolation
195 Tscal dt_interp;
198
199 class acc {
200 public:
201 GetShift<Tvec, TgridVec, AMRBlock> shift_get;
202
203 const Tvec *acc_vel_cell;
204 const Tvec *acc_dx_v_cell;
205 const Tvec *acc_dy_v_cell;
206 const Tvec *acc_dz_v_cell;
207
208 // For time interpolation
209 const Tscal *acc_rho_cell;
210 const Tvec *acc_grad_P_cell;
211
212 Tscal dt_interp;
213
214 acc(const Tvec *aabb_block_lower,
215 const Tscal *aabb_cell_size,
216 const Tvec *vel_cell,
217 const Tvec *dx_v_cell,
218 const Tvec *dy_v_cell,
219 const Tvec *dz_v_cell,
220 // For time interpolation
221 Tscal dt_interp,
222 const Tscal *rho_cell,
223 const Tvec *grad_P_cell)
224 : shift_get(aabb_block_lower, aabb_cell_size), acc_vel_cell{vel_cell},
225 acc_dx_v_cell{dx_v_cell}, acc_dy_v_cell{dy_v_cell}, acc_dz_v_cell{dz_v_cell},
226 dt_interp(dt_interp), acc_rho_cell{rho_cell}, acc_grad_P_cell{grad_P_cell} {}
227
228 Tvec get_dt_v(Tvec v, Tvec dx_v, Tvec dy_v, Tvec dz_v, Tscal rho, Tvec grad_P) const {
229 return -(v[0] * dx_v + v[1] * dy_v + v[2] * dz_v + grad_P / rho);
230 }
231
232 std::array<Tvec, 2> get_link_field_val(u32 id_a, u32 id_b) const {
233
234 auto [shift_a, shift_b] = shift_get.get_shifts(id_a, id_b);
235
236 Tvec v_a = acc_vel_cell[id_a];
237 Tvec dx_vel_a = acc_dx_v_cell[id_a];
238 Tvec dy_vel_a = acc_dy_v_cell[id_a];
239 Tvec dz_vel_a = acc_dz_v_cell[id_a];
240
241 Tvec v_b = acc_vel_cell[id_b];
242 Tvec dx_vel_b = acc_dx_v_cell[id_b];
243 Tvec dy_vel_b = acc_dy_v_cell[id_b];
244 Tvec dz_vel_b = acc_dz_v_cell[id_b];
245
246 Tscal rho_a = acc_rho_cell[id_a];
247 Tvec grad_P_a = acc_grad_P_cell[id_a];
248 Tscal rho_b = acc_rho_cell[id_b];
249 Tvec grad_P_b = acc_grad_P_cell[id_b];
250
251 Tvec dx_v_a_dot_shift
252 = shift_a.x() * dx_vel_a + shift_a.y() * dy_vel_a + shift_a.z() * dz_vel_a;
253 Tvec dx_v_b_dot_shift
254 = shift_b.x() * dx_vel_b + shift_b.y() * dy_vel_b + shift_b.z() * dz_vel_b;
255
256 Tvec dt_v_a = get_dt_v(v_a, dx_vel_a, dy_vel_a, dz_vel_a, rho_a, grad_P_a);
257 Tvec dt_v_b = get_dt_v(v_b, dx_vel_b, dy_vel_b, dz_vel_b, rho_b, grad_P_b);
258
259 Tvec vel_face_a = v_a + dx_v_a_dot_shift + dt_v_a * dt_interp;
260 Tvec vel_face_b = v_b + dx_v_b_dot_shift + dt_v_b * dt_interp;
261
262 return {vel_face_a, vel_face_b};
263 }
264 };
265
266 inline acc get_read_access(sham::EventList &deps) {
267 return acc(
268 aabb_block_lower.get_read_access(deps),
269 aabb_cell_size.get_read_access(deps),
270 vel_cell.get_read_access(deps),
271 dx_v_cell.get_read_access(deps),
272 dy_v_cell.get_read_access(deps),
273 dz_v_cell.get_read_access(deps),
274 // For time interpolation
275 dt_interp,
276 rho_cell.get_read_access(deps),
277 grad_P_cell.get_read_access(deps));
278 }
279
280 inline void complete_event_state(sycl::event e) {
281 aabb_block_lower.complete_event_state(e);
282 aabb_cell_size.complete_event_state(e);
283 vel_cell.complete_event_state(e);
284 dx_v_cell.complete_event_state(e);
285 dy_v_cell.complete_event_state(e);
286 dz_v_cell.complete_event_state(e);
287 rho_cell.complete_event_state(e);
288 grad_P_cell.complete_event_state(e);
289 }
290 };
291
292 template<class Tvec, class TgridVec, class AMRBlock>
293 class PressInterpolate {
294 using Tscal = shambase::VecComponent<Tvec>;
295
296 public:
301 // For time interpolation
302 Tscal dt_interp;
303 Tscal gamma;
308
309 class acc {
310 public:
311 GetShift<Tvec, TgridVec, AMRBlock> shift_get;
312
313 const Tscal *acc_P_cell;
314 const Tvec *acc_grad_P_cell;
315
316 // For time interpolation
317 const Tvec *acc_vel_cell;
318 const Tvec *acc_dx_v_cell;
319 const Tvec *acc_dy_v_cell;
320 const Tvec *acc_dz_v_cell;
321
322 Tscal gamma;
323 Tscal dt_interp;
324
325 acc(const Tvec *aabb_block_lower,
326 const Tscal *aabb_cell_size,
327 const Tscal *P_cell,
328 const Tvec *grad_P_cell,
329 // For time interpolation
330 Tscal dt_interp,
331 Tscal gamma,
332 const Tvec *vel_cell,
333 const Tvec *dx_v_cell,
334 const Tvec *dy_v_cell,
335 const Tvec *dz_v_cell)
336 : shift_get(aabb_block_lower, aabb_cell_size), acc_P_cell{P_cell},
337 acc_grad_P_cell{grad_P_cell}, dt_interp(dt_interp), gamma(gamma),
338 acc_vel_cell{vel_cell}, acc_dx_v_cell{dx_v_cell}, acc_dy_v_cell{dy_v_cell},
339 acc_dz_v_cell{dz_v_cell} {}
340
341 Tscal get_dt_P(
342 Tscal P, Tvec grad_P, Tvec v, Tvec dx_v, Tvec dy_v, Tvec dz_v, Tscal gamma) const {
343 return -(gamma * P * (dx_v[0] + dy_v[1] + dz_v[2]) + sham::dot(v, grad_P));
344 }
345
346 std::array<Tscal, 2> get_link_field_val(u32 id_a, u32 id_b) const {
347
348 auto [shift_a, shift_b] = shift_get.get_shifts(id_a, id_b);
349
350 Tscal P_a = acc_P_cell[id_a];
351 Tvec grad_P_a = acc_grad_P_cell[id_a];
352 Tscal P_b = acc_P_cell[id_b];
353 Tvec grad_P_b = acc_grad_P_cell[id_b];
354
355 Tvec v_a = acc_vel_cell[id_a];
356 Tvec dx_v_a = acc_dx_v_cell[id_a];
357 Tvec dy_v_a = acc_dy_v_cell[id_a];
358 Tvec dz_v_a = acc_dz_v_cell[id_a];
359 Tvec v_b = acc_vel_cell[id_b];
360 Tvec dx_v_b = acc_dx_v_cell[id_b];
361 Tvec dy_v_b = acc_dy_v_cell[id_b];
362 Tvec dz_v_b = acc_dz_v_cell[id_b];
363
364 Tscal dtP_cell_a = get_dt_P(P_a, grad_P_a, v_a, dx_v_a, dy_v_a, dz_v_a, gamma);
365 Tscal dtP_cell_b = get_dt_P(P_b, grad_P_b, v_b, dx_v_b, dy_v_b, dz_v_b, gamma);
366
367 Tscal P_face_a = P_a + sycl::dot(grad_P_a, shift_a) + dtP_cell_a * dt_interp;
368 Tscal P_face_b = P_b + sycl::dot(grad_P_b, shift_b) + dtP_cell_b * dt_interp;
369
370 SHAM_ASSERT(P_face_a >= 0.0);
371 SHAM_ASSERT(P_face_b >= 0.0);
372
373 return {P_face_a, P_face_b};
374 }
375 };
376
377 inline acc get_read_access(sham::EventList &deps) {
378 return acc(
379 aabb_block_lower.get_read_access(deps),
380 aabb_cell_size.get_read_access(deps),
381 P_cell.get_read_access(deps),
382 grad_P_cell.get_read_access(deps),
383 dt_interp,
384 gamma,
385 vel_cell.get_read_access(deps),
386 dx_v_cell.get_read_access(deps),
387 dy_v_cell.get_read_access(deps),
388 dz_v_cell.get_read_access(deps));
389 }
390
391 inline void complete_event_state(sycl::event e) {
392 aabb_block_lower.complete_event_state(e);
393 aabb_cell_size.complete_event_state(e);
394 P_cell.complete_event_state(e);
395 grad_P_cell.complete_event_state(e);
396 vel_cell.complete_event_state(e);
397 dx_v_cell.complete_event_state(e);
398 dy_v_cell.complete_event_state(e);
399 dz_v_cell.complete_event_state(e);
400 }
401 };
402
403 template<class Tvec, class TgridVec, class AMRBlock>
404 class RhoDustInterpolate {
405 using Tscal = shambase::VecComponent<Tvec>;
406
407 public:
408 u32 nvar;
413 // For time interpolation
414 Tscal dt_interp;
419
420 class acc {
421 public:
422 GetShift<Tvec, TgridVec, AMRBlock> shift_get;
423 u32 nvar;
424
425 const Tscal *acc_rho_dust_cell;
426 const Tvec *acc_grad_rho_dust_cell;
427
428 // For time interpolation
429 const Tvec *acc_vel_dust_cell;
430 const Tvec *acc_dx_v_dust_cell;
431 const Tvec *acc_dy_v_dust_cell;
432 const Tvec *acc_dz_v_dust_cell;
433
434 Tscal dt_interp;
435
436 acc(u32 nvar,
437 const Tvec *aabb_block_lower,
438 const Tscal *aabb_cell_size,
439 const Tscal *rho_dust_cell,
440 const Tvec *grad_rho_dust_cell,
441 // For time interpolation
442 Tscal dt_interp,
443 const Tvec *vel_dust_cell,
444 const Tvec *dx_v_dust_cell,
445 const Tvec *dy_v_dust_cell,
446 const Tvec *dz_v_dust_cell)
447 : shift_get(aabb_block_lower, aabb_cell_size), nvar(nvar),
448 acc_rho_dust_cell{rho_dust_cell}, acc_grad_rho_dust_cell{grad_rho_dust_cell},
449 dt_interp(dt_interp), acc_vel_dust_cell{vel_dust_cell},
450 acc_dx_v_dust_cell{dx_v_dust_cell}, acc_dy_v_dust_cell{dy_v_dust_cell},
451 acc_dz_v_dust_cell{dz_v_dust_cell} {}
452
453 Tscal get_dt_rho_dust(
454 Tscal rho_dust,
455 Tvec v_dust,
456 Tvec grad_rho_dust,
457 Tvec dx_v_dust,
458 Tvec dy_v_dust,
459 Tvec dz_v_dust) const {
460 return -(
461 sham::dot(v_dust, grad_rho_dust)
462 + rho_dust * (dx_v_dust[0] + dy_v_dust[1] + dz_v_dust[2]));
463 }
464
465 std::array<Tscal, 2> get_link_field_val(u32 id_a, u32 id_b) const {
466 const u32 icell_a = id_a / nvar;
467 const u32 icell_b = id_b / nvar;
468
469 auto [shift_a, shift_b] = shift_get.get_shifts(icell_a, icell_b);
470
471 Tscal rho_dust_a = acc_rho_dust_cell[id_a];
472 Tvec grad_rho_dust_a = acc_grad_rho_dust_cell[id_a];
473 Tscal rho_dust_b = acc_rho_dust_cell[id_b];
474 Tvec grad_rho_dust_b = acc_grad_rho_dust_cell[id_b];
475
476 Tvec vel_dust_a = acc_vel_dust_cell[id_a];
477 Tvec dx_v_dust_a = acc_dx_v_dust_cell[id_a];
478 Tvec dy_v_dust_a = acc_dy_v_dust_cell[id_a];
479 Tvec dz_v_dust_a = acc_dz_v_dust_cell[id_a];
480 Tvec vel_dust_b = acc_vel_dust_cell[id_b];
481 Tvec dx_v_dust_b = acc_dx_v_dust_cell[id_b];
482 Tvec dy_v_dust_b = acc_dy_v_dust_cell[id_b];
483 Tvec dz_v_dust_b = acc_dz_v_dust_cell[id_b];
484
485 Tscal rho_dust_face_a = rho_dust_a + sycl::dot(grad_rho_dust_a, shift_a);
486 Tscal rho_dust_face_b = rho_dust_b + sycl::dot(grad_rho_dust_b, shift_b);
487
488 rho_dust_face_a += get_dt_rho_dust(
489 rho_dust_a,
490 vel_dust_a,
491 grad_rho_dust_a,
492 dx_v_dust_a,
493 dy_v_dust_a,
494 dz_v_dust_a)
495 * dt_interp;
496 rho_dust_face_b += get_dt_rho_dust(
497 rho_dust_b,
498 vel_dust_b,
499 grad_rho_dust_b,
500 dx_v_dust_b,
501 dy_v_dust_b,
502 dz_v_dust_b)
503 * dt_interp;
504
505 return {rho_dust_face_a, rho_dust_face_b};
506 }
507 };
508
509 inline acc get_read_access(sham::EventList &deps) {
510 return acc(
511 nvar,
512 aabb_block_lower.get_read_access(deps),
513 aabb_cell_size.get_read_access(deps),
514 rho_dust_cell.get_read_access(deps),
515 grad_rho_dust_cell.get_read_access(deps),
516 // For time interpolation
517 dt_interp,
518 vel_dust_cell.get_read_access(deps),
519 dx_v_dust_cell.get_read_access(deps),
520 dy_v_dust_cell.get_read_access(deps),
521 dz_v_dust_cell.get_read_access(deps));
522 }
523
524 inline void complete_event_state(sycl::event e) {
525 aabb_block_lower.complete_event_state(e);
526 aabb_cell_size.complete_event_state(e);
527 rho_dust_cell.complete_event_state(e);
528 grad_rho_dust_cell.complete_event_state(e);
529 vel_dust_cell.complete_event_state(e);
530 dx_v_dust_cell.complete_event_state(e);
531 dy_v_dust_cell.complete_event_state(e);
532 dz_v_dust_cell.complete_event_state(e);
533 }
534 };
535
536 template<class Tvec, class TgridVec, class AMRBlock>
537 class VelDustInterpolate {
538 using Tscal = shambase::VecComponent<Tvec>;
539
540 public:
541 u32 nvar;
548 // For time interpolation
549 Tscal dt_interp;
551
552 class acc {
553 public:
554 GetShift<Tvec, TgridVec, AMRBlock> shift_get;
555 u32 nvar;
556
557 const Tvec *acc_vel_dust_cell;
558 const Tvec *acc_dx_v_dust_cell;
559 const Tvec *acc_dy_v_dust_cell;
560 const Tvec *acc_dz_v_dust_cell;
561
562 // For time interpolation
563 const Tscal *acc_rho_dust_cell;
564
565 Tscal dt_interp;
566
567 acc(u32 nvar,
568 const Tvec *aabb_block_lower,
569 const Tscal *aabb_cell_size,
570 const Tvec *vel_dust_cell,
571 const Tvec *dx_v_dust_cell,
572 const Tvec *dy_v_dust_cell,
573 const Tvec *dz_v_dust_cell,
574 // For time interpolation
575 Tscal dt_interp,
576 const Tscal *rho_dust_cell)
577 : shift_get(aabb_block_lower, aabb_cell_size), nvar(nvar),
578 acc_vel_dust_cell{vel_dust_cell}, acc_dx_v_dust_cell{dx_v_dust_cell},
579 acc_dy_v_dust_cell{dy_v_dust_cell}, acc_dz_v_dust_cell{dz_v_dust_cell},
580 dt_interp(dt_interp), acc_rho_dust_cell{rho_dust_cell} {}
581
582 Tvec get_dt_v_dust(Tvec v, Tvec dx_v, Tvec dy_v, Tvec dz_v, Tscal rho) const {
583 return -(v[0] * dx_v + v[1] * dy_v + v[2] * dz_v);
584 }
585
586 std::array<Tvec, 2> get_link_field_val(u32 id_a, u32 id_b) const {
587 const u32 icell_a = id_a / nvar;
588 const u32 icell_b = id_b / nvar;
589
590 auto [shift_a, shift_b] = shift_get.get_shifts(icell_a, icell_b);
591
592 Tvec v_dust_a = acc_vel_dust_cell[id_a];
593 Tvec dx_vel_dust_a = acc_dx_v_dust_cell[id_a];
594 Tvec dy_vel_dust_a = acc_dy_v_dust_cell[id_a];
595 Tvec dz_vel_dust_a = acc_dz_v_dust_cell[id_a];
596
597 Tvec v_dust_b = acc_vel_dust_cell[id_b];
598 Tvec dx_vel_dust_b = acc_dx_v_dust_cell[id_b];
599 Tvec dy_vel_dust_b = acc_dy_v_dust_cell[id_b];
600 Tvec dz_vel_dust_b = acc_dz_v_dust_cell[id_b];
601
602 Tscal rho_dust_a = acc_rho_dust_cell[id_a];
603 Tscal rho_dust_b = acc_rho_dust_cell[id_b];
604
605 Tvec dx_v_dust_a_dot_shift = shift_a.x() * dx_vel_dust_a
606 + shift_a.y() * dy_vel_dust_a
607 + shift_a.z() * dz_vel_dust_a;
608 Tvec dx_v_dust_b_dot_shift = shift_b.x() * dx_vel_dust_b
609 + shift_b.y() * dy_vel_dust_b
610 + shift_b.z() * dz_vel_dust_b;
611
612 Tvec dt_v_dust_a = get_dt_v_dust(
613 v_dust_a, dx_vel_dust_a, dy_vel_dust_a, dz_vel_dust_a, rho_dust_a);
614 Tvec dt_v_dust_b = get_dt_v_dust(
615 v_dust_b, dx_vel_dust_b, dy_vel_dust_b, dz_vel_dust_b, rho_dust_b);
616
617 Tvec vel_dust_face_a = v_dust_a + dx_v_dust_a_dot_shift + dt_v_dust_a * dt_interp;
618 Tvec vel_dust_face_b = v_dust_b + dx_v_dust_b_dot_shift + dt_v_dust_b * dt_interp;
619
620 return {vel_dust_face_a, vel_dust_face_b};
621 }
622 };
623
624 inline acc get_read_access(sham::EventList &deps) {
625 return acc(
626 nvar,
627 aabb_block_lower.get_read_access(deps),
628 aabb_cell_size.get_read_access(deps),
629 vel_dust_cell.get_read_access(deps),
630 dx_v_dust_cell.get_read_access(deps),
631 dy_v_dust_cell.get_read_access(deps),
632 dz_v_dust_cell.get_read_access(deps),
633 // For time interpolation
634 dt_interp,
635 rho_dust_cell.get_read_access(deps));
636 }
637
638 inline void complete_event_state(sycl::event e) {
639 aabb_block_lower.complete_event_state(e);
640 aabb_cell_size.complete_event_state(e);
641 vel_dust_cell.complete_event_state(e);
642 dx_v_dust_cell.complete_event_state(e);
643 dy_v_dust_cell.complete_event_state(e);
644 dz_v_dust_cell.complete_event_state(e);
645 rho_dust_cell.complete_event_state(e);
646 }
647 };
648
649} // namespace
650
651template<class Tvec, class TgridVec>
654 StackEntry stack_loc{};
655
657
658 static constexpr u32 NsideBlockPow = 1;
660
661 SHAM_ASSERT(AMRBlock::block_size == block_size);
662
663 auto edges = get_edges();
664
665 auto dt_interp = edges.dt_interp.value;
666
667 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &rho_face_xp = edges.rho_face_xp;
668 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &rho_face_xm = edges.rho_face_xm;
669 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &rho_face_yp = edges.rho_face_yp;
670 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &rho_face_ym = edges.rho_face_ym;
671 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &rho_face_zp = edges.rho_face_zp;
672 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &rho_face_zm = edges.rho_face_zm;
673
674 rho_face_xp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::xp));
675 rho_face_xm.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::xm));
676 rho_face_yp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::yp));
677 rho_face_ym.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::ym));
678 rho_face_zp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::zp));
679 rho_face_zm.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::zm));
680
681 auto spans_block_cell_sizes = edges.spans_block_cell_sizes.get_spans();
682 auto spans_cell0block_aabb_lower = edges.spans_cell0block_aabb_lower.get_spans();
683 auto spans_rhos = edges.spans_rhos.get_spans();
684 auto spans_grad_rho = edges.spans_grad_rho.get_spans();
685 auto spans_vel = edges.spans_vel.get_spans();
686 auto spans_dx_vel = edges.spans_dx_vel.get_spans();
687 auto spans_dy_vel = edges.spans_dy_vel.get_spans();
688 auto spans_dz_vel = edges.spans_dz_vel.get_spans();
689
690 using Interp = RhoInterpolate<Tvec, TgridVec, AMRBlock>;
691 auto interpolators
692 = spans_block_cell_sizes.template map<Interp>([&](u64 id, auto &csize) -> Interp {
693 return {
694 spans_cell0block_aabb_lower.get(id),
695 spans_block_cell_sizes.get(id),
696 spans_rhos.get(id),
697 spans_grad_rho.get(id),
698 dt_interp,
699 spans_vel.get(id),
700 spans_dx_vel.get(id),
701 spans_dy_vel.get(id),
702 spans_dz_vel.get(id)};
703 });
704
705 auto graphs_xp = edges.cell_neigh_graph.get_refs_dir(Direction::xp);
706 auto graphs_xm = edges.cell_neigh_graph.get_refs_dir(Direction::xm);
707 auto graphs_yp = edges.cell_neigh_graph.get_refs_dir(Direction::yp);
708 auto graphs_ym = edges.cell_neigh_graph.get_refs_dir(Direction::ym);
709 auto graphs_zp = edges.cell_neigh_graph.get_refs_dir(Direction::zp);
710 auto graphs_zm = edges.cell_neigh_graph.get_refs_dir(Direction::zm);
711
713 = graphs_xp.template map<u32>([&](u64 id, auto &graph) {
714 return graph.get().obj_cnt;
715 });
717 = graphs_xm.template map<u32>([&](u64 id, auto &graph) {
718 return graph.get().obj_cnt;
719 });
721 = graphs_yp.template map<u32>([&](u64 id, auto &graph) {
722 return graph.get().obj_cnt;
723 });
725 = graphs_ym.template map<u32>([&](u64 id, auto &graph) {
726 return graph.get().obj_cnt;
727 });
729 = graphs_zp.template map<u32>([&](u64 id, auto &graph) {
730 return graph.get().obj_cnt;
731 });
733 = graphs_zm.template map<u32>([&](u64 id, auto &graph) {
734 return graph.get().obj_cnt;
735 });
736
738 shamsys::instance::get_compute_scheduler_ptr(),
739 sham::DDMultiRef{graphs_xp, interpolators},
740 sham::DDMultiRef{rho_face_xp.link_fields},
741 counts_xp,
742 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
743 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
744 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
745 });
746 });
748 shamsys::instance::get_compute_scheduler_ptr(),
749 sham::DDMultiRef{graphs_xm, interpolators},
750 sham::DDMultiRef{rho_face_xm.link_fields},
751 counts_xm,
752 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
753 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
754 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
755 });
756 });
758 shamsys::instance::get_compute_scheduler_ptr(),
759 sham::DDMultiRef{graphs_yp, interpolators},
760 sham::DDMultiRef{rho_face_yp.link_fields},
761 counts_yp,
762 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
763 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
764 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
765 });
766 });
768 shamsys::instance::get_compute_scheduler_ptr(),
769 sham::DDMultiRef{graphs_ym, interpolators},
770 sham::DDMultiRef{rho_face_ym.link_fields},
771 counts_ym,
772 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
773 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
774 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
775 });
776 });
778 shamsys::instance::get_compute_scheduler_ptr(),
779 sham::DDMultiRef{graphs_zp, interpolators},
780 sham::DDMultiRef{rho_face_zp.link_fields},
781 counts_zp,
782 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
783 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
784 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
785 });
786 });
788 shamsys::instance::get_compute_scheduler_ptr(),
789 sham::DDMultiRef{graphs_zm, interpolators},
790 sham::DDMultiRef{rho_face_zm.link_fields},
791 counts_zm,
792 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
793 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
794 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
795 });
796 });
797}
798
799template<class Tvec, class TgridVec>
804
806
807template<class Tvec, class TgridVec>
810 StackEntry stack_loc{};
811
813
814 static constexpr u32 NsideBlockPow = 1;
816
817 SHAM_ASSERT(AMRBlock::block_size == block_size);
818
819 auto edges = get_edges();
820
821 auto dt_interp = edges.dt_interp.value;
822
823 solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>> &vel_face_xp = edges.vel_face_xp;
824 solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>> &vel_face_xm = edges.vel_face_xm;
825 solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>> &vel_face_yp = edges.vel_face_yp;
826 solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>> &vel_face_ym = edges.vel_face_ym;
827 solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>> &vel_face_zp = edges.vel_face_zp;
828 solvergraph::NeighGraphLinkFieldEdge<std::array<Tvec, 2>> &vel_face_zm = edges.vel_face_zm;
829
830 vel_face_xp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::xp));
831 vel_face_xm.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::xm));
832 vel_face_yp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::yp));
833 vel_face_ym.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::ym));
834 vel_face_zp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::zp));
835 vel_face_zm.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::zm));
836
837 auto spans_block_cell_sizes = edges.spans_block_cell_sizes.get_spans();
838 auto spans_cell0block_aabb_lower = edges.spans_cell0block_aabb_lower.get_spans();
839 auto spans_rhos = edges.spans_rhos.get_spans();
840 auto spans_grad_P = edges.spans_grad_P.get_spans();
841 auto spans_vel = edges.spans_vel.get_spans();
842 auto spans_dx_vel = edges.spans_dx_vel.get_spans();
843 auto spans_dy_vel = edges.spans_dy_vel.get_spans();
844 auto spans_dz_vel = edges.spans_dz_vel.get_spans();
845
846 using Interp = VelInterpolate<Tvec, TgridVec, AMRBlock>;
847 auto interpolators
848 = spans_block_cell_sizes.template map<Interp>([&](u64 id, auto &csize) -> Interp {
849 return {
850 spans_cell0block_aabb_lower.get(id),
851 spans_block_cell_sizes.get(id),
852 spans_vel.get(id),
853 spans_dx_vel.get(id),
854 spans_dy_vel.get(id),
855 spans_dz_vel.get(id),
856 dt_interp,
857 spans_rhos.get(id),
858 spans_grad_P.get(id)};
859 });
860
861 auto graphs_xp = edges.cell_neigh_graph.get_refs_dir(Direction::xp);
862 auto graphs_xm = edges.cell_neigh_graph.get_refs_dir(Direction::xm);
863 auto graphs_yp = edges.cell_neigh_graph.get_refs_dir(Direction::yp);
864 auto graphs_ym = edges.cell_neigh_graph.get_refs_dir(Direction::ym);
865 auto graphs_zp = edges.cell_neigh_graph.get_refs_dir(Direction::zp);
866 auto graphs_zm = edges.cell_neigh_graph.get_refs_dir(Direction::zm);
867
869 = graphs_xp.template map<u32>([&](u64 id, auto &graph) {
870 return graph.get().obj_cnt;
871 });
873 = graphs_xm.template map<u32>([&](u64 id, auto &graph) {
874 return graph.get().obj_cnt;
875 });
877 = graphs_yp.template map<u32>([&](u64 id, auto &graph) {
878 return graph.get().obj_cnt;
879 });
881 = graphs_ym.template map<u32>([&](u64 id, auto &graph) {
882 return graph.get().obj_cnt;
883 });
885 = graphs_zp.template map<u32>([&](u64 id, auto &graph) {
886 return graph.get().obj_cnt;
887 });
889 = graphs_zm.template map<u32>([&](u64 id, auto &graph) {
890 return graph.get().obj_cnt;
891 });
892
894 shamsys::instance::get_compute_scheduler_ptr(),
895 sham::DDMultiRef{graphs_xp, interpolators},
896 sham::DDMultiRef{vel_face_xp.link_fields},
897 counts_xp,
898 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
899 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
900 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
901 });
902 });
904 shamsys::instance::get_compute_scheduler_ptr(),
905 sham::DDMultiRef{graphs_xm, interpolators},
906 sham::DDMultiRef{vel_face_xm.link_fields},
907 counts_xm,
908 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
909 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
910 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
911 });
912 });
914 shamsys::instance::get_compute_scheduler_ptr(),
915 sham::DDMultiRef{graphs_yp, interpolators},
916 sham::DDMultiRef{vel_face_yp.link_fields},
917 counts_yp,
918 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
919 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
920 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
921 });
922 });
924 shamsys::instance::get_compute_scheduler_ptr(),
925 sham::DDMultiRef{graphs_ym, interpolators},
926 sham::DDMultiRef{vel_face_ym.link_fields},
927 counts_ym,
928 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
929 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
930 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
931 });
932 });
934 shamsys::instance::get_compute_scheduler_ptr(),
935 sham::DDMultiRef{graphs_zp, interpolators},
936 sham::DDMultiRef{vel_face_zp.link_fields},
937 counts_zp,
938 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
939 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
940 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
941 });
942 });
944 shamsys::instance::get_compute_scheduler_ptr(),
945 sham::DDMultiRef{graphs_zm, interpolators},
946 sham::DDMultiRef{vel_face_zm.link_fields},
947 counts_zm,
948 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
949 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
950 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
951 });
952 });
953}
954
955template<class Tvec, class TgridVec>
960
962
963template<class Tvec, class TgridVec>
966 StackEntry stack_loc{};
967
969
970 static constexpr u32 NsideBlockPow = 1;
972
973 SHAM_ASSERT(AMRBlock::block_size == block_size);
974
975 auto edges = get_edges();
976
977 auto dt_interp = edges.dt_interp.value;
978
979 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &press_face_xp = edges.press_face_xp;
980 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &press_face_xm = edges.press_face_xm;
981 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &press_face_yp = edges.press_face_yp;
982 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &press_face_ym = edges.press_face_ym;
983 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &press_face_zp = edges.press_face_zp;
984 solvergraph::NeighGraphLinkFieldEdge<std::array<Tscal, 2>> &press_face_zm = edges.press_face_zm;
985
986 press_face_xp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::xp));
987 press_face_xm.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::xm));
988 press_face_yp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::yp));
989 press_face_ym.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::ym));
990 press_face_zp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::zp));
991 press_face_zm.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::zm));
992
993 auto spans_block_cell_sizes = edges.spans_block_cell_sizes.get_spans();
994 auto spans_cell0block_aabb_lower = edges.spans_cell0block_aabb_lower.get_spans();
995 auto spans_press = edges.spans_press.get_spans();
996 auto spans_grad_P = edges.spans_grad_P.get_spans();
997 auto spans_vel = edges.spans_vel.get_spans();
998 auto spans_dx_vel = edges.spans_dx_vel.get_spans();
999 auto spans_dy_vel = edges.spans_dy_vel.get_spans();
1000 auto spans_dz_vel = edges.spans_dz_vel.get_spans();
1001
1002 using Interp = PressInterpolate<Tvec, TgridVec, AMRBlock>;
1003 auto interpolators
1004 = spans_block_cell_sizes.template map<Interp>([&](u64 id, auto &csize) -> Interp {
1005 return {
1006 spans_cell0block_aabb_lower.get(id),
1007 spans_block_cell_sizes.get(id),
1008 spans_press.get(id),
1009 spans_grad_P.get(id),
1010 dt_interp,
1011 gamma,
1012 spans_vel.get(id),
1013 spans_dx_vel.get(id),
1014 spans_dy_vel.get(id),
1015 spans_dz_vel.get(id)};
1016 });
1017
1018 auto graphs_xp = edges.cell_neigh_graph.get_refs_dir(Direction::xp);
1019 auto graphs_xm = edges.cell_neigh_graph.get_refs_dir(Direction::xm);
1020 auto graphs_yp = edges.cell_neigh_graph.get_refs_dir(Direction::yp);
1021 auto graphs_ym = edges.cell_neigh_graph.get_refs_dir(Direction::ym);
1022 auto graphs_zp = edges.cell_neigh_graph.get_refs_dir(Direction::zp);
1023 auto graphs_zm = edges.cell_neigh_graph.get_refs_dir(Direction::zm);
1024
1026 = graphs_xp.template map<u32>([&](u64 id, auto &graph) {
1027 return graph.get().obj_cnt;
1028 });
1030 = graphs_xm.template map<u32>([&](u64 id, auto &graph) {
1031 return graph.get().obj_cnt;
1032 });
1034 = graphs_yp.template map<u32>([&](u64 id, auto &graph) {
1035 return graph.get().obj_cnt;
1036 });
1038 = graphs_ym.template map<u32>([&](u64 id, auto &graph) {
1039 return graph.get().obj_cnt;
1040 });
1042 = graphs_zp.template map<u32>([&](u64 id, auto &graph) {
1043 return graph.get().obj_cnt;
1044 });
1046 = graphs_zm.template map<u32>([&](u64 id, auto &graph) {
1047 return graph.get().obj_cnt;
1048 });
1049
1051 shamsys::instance::get_compute_scheduler_ptr(),
1052 sham::DDMultiRef{graphs_xp, interpolators},
1053 sham::DDMultiRef{press_face_xp.link_fields},
1054 counts_xp,
1055 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
1056 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
1057 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
1058 });
1059 });
1061 shamsys::instance::get_compute_scheduler_ptr(),
1062 sham::DDMultiRef{graphs_xm, interpolators},
1063 sham::DDMultiRef{press_face_xm.link_fields},
1064 counts_xm,
1065 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
1066 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
1067 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
1068 });
1069 });
1071 shamsys::instance::get_compute_scheduler_ptr(),
1072 sham::DDMultiRef{graphs_yp, interpolators},
1073 sham::DDMultiRef{press_face_yp.link_fields},
1074 counts_yp,
1075 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
1076 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
1077 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
1078 });
1079 });
1081 shamsys::instance::get_compute_scheduler_ptr(),
1082 sham::DDMultiRef{graphs_ym, interpolators},
1083 sham::DDMultiRef{press_face_ym.link_fields},
1084 counts_ym,
1085 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
1086 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
1087 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
1088 });
1089 });
1091 shamsys::instance::get_compute_scheduler_ptr(),
1092 sham::DDMultiRef{graphs_zp, interpolators},
1093 sham::DDMultiRef{press_face_zp.link_fields},
1094 counts_zp,
1095 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
1096 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
1097 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
1098 });
1099 });
1101 shamsys::instance::get_compute_scheduler_ptr(),
1102 sham::DDMultiRef{graphs_zm, interpolators},
1103 sham::DDMultiRef{press_face_zm.link_fields},
1104 counts_zm,
1105 [](u32 id_a, auto link_iter, auto compute, auto acc_link_field) {
1106 link_iter.for_each_object_link_id(id_a, [&](u32 id_b, u32 link_id) {
1107 acc_link_field[link_id] = compute.get_link_field_val(id_a, id_b);
1108 });
1109 });
1110}
1111
1112template<class Tvec, class TgridVec>
1117
1119
1120template<class Tvec, class TgridVec>
1123 StackEntry stack_loc{};
1124
1126
1127 static constexpr u32 NsideBlockPow = 1;
1129
1130 SHAM_ASSERT(AMRBlock::block_size == block_size);
1131
1132 auto edges = get_edges();
1133
1134 auto dt_interp = edges.dt_interp.value;
1135 auto ndust = this->ndust;
1136
1138 = edges.rho_dust_face_xp;
1140 = edges.rho_dust_face_xm;
1142 = edges.rho_dust_face_yp;
1144 = edges.rho_dust_face_ym;
1146 = edges.rho_dust_face_zp;
1148 = edges.rho_dust_face_zm;
1149
1150 rho_dust_face_xp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::xp));
1151 rho_dust_face_xm.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::xm));
1152 rho_dust_face_yp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::yp));
1153 rho_dust_face_ym.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::ym));
1154 rho_dust_face_zp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::zp));
1155 rho_dust_face_zm.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::zm));
1156
1157 auto spans_block_cell_sizes = edges.spans_block_cell_sizes.get_spans();
1158 auto spans_cell0block_aabb_lower = edges.spans_cell0block_aabb_lower.get_spans();
1159 auto spans_rhos_dust = edges.spans_rhos_dust.get_spans();
1160 auto spans_grad_rho_dust = edges.spans_grad_rho_dust.get_spans();
1161 auto spans_vel_dust = edges.spans_vel_dust.get_spans();
1162 auto spans_dx_vel_dust = edges.spans_dx_vel_dust.get_spans();
1163 auto spans_dy_vel_dust = edges.spans_dy_vel_dust.get_spans();
1164 auto spans_dz_vel_dust = edges.spans_dz_vel_dust.get_spans();
1165
1166 using Interp = RhoDustInterpolate<Tvec, TgridVec, AMRBlock>;
1167 auto interpolators
1168 = spans_block_cell_sizes.template map<Interp>([&](u64 id, auto &csize) -> Interp {
1169 return {
1170 ndust,
1171 spans_cell0block_aabb_lower.get(id),
1172 spans_block_cell_sizes.get(id),
1173 spans_rhos_dust.get(id),
1174 spans_grad_rho_dust.get(id),
1175 dt_interp,
1176 spans_vel_dust.get(id),
1177 spans_dx_vel_dust.get(id),
1178 spans_dy_vel_dust.get(id),
1179 spans_dz_vel_dust.get(id)};
1180 });
1181
1182 auto graphs_xp = edges.cell_neigh_graph.get_refs_dir(Direction::xp);
1183 auto graphs_xm = edges.cell_neigh_graph.get_refs_dir(Direction::xm);
1184 auto graphs_yp = edges.cell_neigh_graph.get_refs_dir(Direction::yp);
1185 auto graphs_ym = edges.cell_neigh_graph.get_refs_dir(Direction::ym);
1186 auto graphs_zp = edges.cell_neigh_graph.get_refs_dir(Direction::zp);
1187 auto graphs_zm = edges.cell_neigh_graph.get_refs_dir(Direction::zm);
1188
1190 = graphs_xp.template map<u32>([&](u64 id, auto &graph) {
1191 return graph.get().obj_cnt * ndust;
1192 });
1194 = graphs_xm.template map<u32>([&](u64 id, auto &graph) {
1195 return graph.get().obj_cnt * ndust;
1196 });
1198 = graphs_yp.template map<u32>([&](u64 id, auto &graph) {
1199 return graph.get().obj_cnt * ndust;
1200 });
1202 = graphs_ym.template map<u32>([&](u64 id, auto &graph) {
1203 return graph.get().obj_cnt * ndust;
1204 });
1206 = graphs_zp.template map<u32>([&](u64 id, auto &graph) {
1207 return graph.get().obj_cnt * ndust;
1208 });
1210 = graphs_zm.template map<u32>([&](u64 id, auto &graph) {
1211 return graph.get().obj_cnt * ndust;
1212 });
1213
1215 shamsys::instance::get_compute_scheduler_ptr(),
1216 sham::DDMultiRef{graphs_xp, interpolators},
1217 sham::DDMultiRef{rho_dust_face_xp.link_fields},
1218 counts_xp,
1219 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1220 const u32 id_cell_a = idvar_a / ndust;
1221 const u32 nvar_loc = idvar_a % ndust;
1222 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1223 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1224 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1225 });
1226 });
1228 shamsys::instance::get_compute_scheduler_ptr(),
1229 sham::DDMultiRef{graphs_xm, interpolators},
1230 sham::DDMultiRef{rho_dust_face_xm.link_fields},
1231 counts_xm,
1232 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1233 const u32 id_cell_a = idvar_a / ndust;
1234 const u32 nvar_loc = idvar_a % ndust;
1235 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1236 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1237 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1238 });
1239 });
1241 shamsys::instance::get_compute_scheduler_ptr(),
1242 sham::DDMultiRef{graphs_yp, interpolators},
1243 sham::DDMultiRef{rho_dust_face_yp.link_fields},
1244 counts_yp,
1245 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1246 const u32 id_cell_a = idvar_a / ndust;
1247 const u32 nvar_loc = idvar_a % ndust;
1248 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1249 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1250 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1251 });
1252 });
1254 shamsys::instance::get_compute_scheduler_ptr(),
1255 sham::DDMultiRef{graphs_ym, interpolators},
1256 sham::DDMultiRef{rho_dust_face_ym.link_fields},
1257 counts_ym,
1258 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1259 const u32 id_cell_a = idvar_a / ndust;
1260 const u32 nvar_loc = idvar_a % ndust;
1261 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1262 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1263 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1264 });
1265 });
1267 shamsys::instance::get_compute_scheduler_ptr(),
1268 sham::DDMultiRef{graphs_zp, interpolators},
1269 sham::DDMultiRef{rho_dust_face_zp.link_fields},
1270 counts_zp,
1271 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1272 const u32 id_cell_a = idvar_a / ndust;
1273 const u32 nvar_loc = idvar_a % ndust;
1274 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1275 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1276 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1277 });
1278 });
1280 shamsys::instance::get_compute_scheduler_ptr(),
1281 sham::DDMultiRef{graphs_zm, interpolators},
1282 sham::DDMultiRef{rho_dust_face_zm.link_fields},
1283 counts_zm,
1284 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1285 const u32 id_cell_a = idvar_a / ndust;
1286 const u32 nvar_loc = idvar_a % ndust;
1287 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1288 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1289 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1290 });
1291 });
1292}
1293
1294template<class Tvec, class TgridVec>
1299
1301
1302template<class Tvec, class TgridVec>
1305 StackEntry stack_loc{};
1306
1308
1309 static constexpr u32 NsideBlockPow = 1;
1311
1312 SHAM_ASSERT(AMRBlock::block_size == block_size);
1313
1314 auto edges = get_edges();
1315
1316 auto dt_interp = edges.dt_interp.value;
1317 auto ndust = this->ndust;
1318
1320 = edges.vel_dust_face_xp;
1322 = edges.vel_dust_face_xm;
1324 = edges.vel_dust_face_yp;
1326 = edges.vel_dust_face_ym;
1328 = edges.vel_dust_face_zp;
1330 = edges.vel_dust_face_zm;
1331
1332 vel_dust_face_xp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::xp));
1333 vel_dust_face_xm.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::xm));
1334 vel_dust_face_yp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::yp));
1335 vel_dust_face_ym.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::ym));
1336 vel_dust_face_zp.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::zp));
1337 vel_dust_face_zm.resize_according_to(edges.cell_neigh_graph.get_refs_dir(Direction::zm));
1338
1339 auto spans_block_cell_sizes = edges.spans_block_cell_sizes.get_spans();
1340 auto spans_cell0block_aabb_lower = edges.spans_cell0block_aabb_lower.get_spans();
1341 auto spans_rhos_dust = edges.spans_rhos_dust.get_spans();
1342 auto spans_vel_dust = edges.spans_vel_dust.get_spans();
1343 auto spans_dx_vel_dust = edges.spans_dx_vel_dust.get_spans();
1344 auto spans_dy_vel_dust = edges.spans_dy_vel_dust.get_spans();
1345 auto spans_dz_vel_dust = edges.spans_dz_vel_dust.get_spans();
1346
1347 using Interp = VelDustInterpolate<Tvec, TgridVec, AMRBlock>;
1348 auto interpolators
1349 = spans_block_cell_sizes.template map<Interp>([&](u64 id, auto &csize) -> Interp {
1350 return {
1351 ndust,
1352 spans_cell0block_aabb_lower.get(id),
1353 spans_block_cell_sizes.get(id),
1354 spans_vel_dust.get(id),
1355 spans_dx_vel_dust.get(id),
1356 spans_dy_vel_dust.get(id),
1357 spans_dz_vel_dust.get(id),
1358 dt_interp,
1359 spans_rhos_dust.get(id)};
1360 });
1361
1362 auto graphs_xp = edges.cell_neigh_graph.get_refs_dir(Direction::xp);
1363 auto graphs_xm = edges.cell_neigh_graph.get_refs_dir(Direction::xm);
1364 auto graphs_yp = edges.cell_neigh_graph.get_refs_dir(Direction::yp);
1365 auto graphs_ym = edges.cell_neigh_graph.get_refs_dir(Direction::ym);
1366 auto graphs_zp = edges.cell_neigh_graph.get_refs_dir(Direction::zp);
1367 auto graphs_zm = edges.cell_neigh_graph.get_refs_dir(Direction::zm);
1368
1370 = graphs_xp.template map<u32>([&](u64 id, auto &graph) {
1371 return graph.get().obj_cnt * ndust;
1372 });
1374 = graphs_xm.template map<u32>([&](u64 id, auto &graph) {
1375 return graph.get().obj_cnt * ndust;
1376 });
1378 = graphs_yp.template map<u32>([&](u64 id, auto &graph) {
1379 return graph.get().obj_cnt * ndust;
1380 });
1382 = graphs_ym.template map<u32>([&](u64 id, auto &graph) {
1383 return graph.get().obj_cnt * ndust;
1384 });
1386 = graphs_zp.template map<u32>([&](u64 id, auto &graph) {
1387 return graph.get().obj_cnt * ndust;
1388 });
1390 = graphs_zm.template map<u32>([&](u64 id, auto &graph) {
1391 return graph.get().obj_cnt * ndust;
1392 });
1393
1395 shamsys::instance::get_compute_scheduler_ptr(),
1396 sham::DDMultiRef{graphs_xp, interpolators},
1397 sham::DDMultiRef{vel_dust_face_xp.link_fields},
1398 counts_xp,
1399 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1400 const u32 id_cell_a = idvar_a / ndust;
1401 const u32 nvar_loc = idvar_a % ndust;
1402 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1403 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1404 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1405 });
1406 });
1408 shamsys::instance::get_compute_scheduler_ptr(),
1409 sham::DDMultiRef{graphs_xm, interpolators},
1410 sham::DDMultiRef{vel_dust_face_xm.link_fields},
1411 counts_xm,
1412 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1413 const u32 id_cell_a = idvar_a / ndust;
1414 const u32 nvar_loc = idvar_a % ndust;
1415 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1416 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1417 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1418 });
1419 });
1421 shamsys::instance::get_compute_scheduler_ptr(),
1422 sham::DDMultiRef{graphs_yp, interpolators},
1423 sham::DDMultiRef{vel_dust_face_yp.link_fields},
1424 counts_yp,
1425 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1426 const u32 id_cell_a = idvar_a / ndust;
1427 const u32 nvar_loc = idvar_a % ndust;
1428 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1429 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1430 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1431 });
1432 });
1434 shamsys::instance::get_compute_scheduler_ptr(),
1435 sham::DDMultiRef{graphs_ym, interpolators},
1436 sham::DDMultiRef{vel_dust_face_ym.link_fields},
1437 counts_ym,
1438 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1439 const u32 id_cell_a = idvar_a / ndust;
1440 const u32 nvar_loc = idvar_a % ndust;
1441 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1442 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1443 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1444 });
1445 });
1447 shamsys::instance::get_compute_scheduler_ptr(),
1448 sham::DDMultiRef{graphs_zp, interpolators},
1449 sham::DDMultiRef{vel_dust_face_zp.link_fields},
1450 counts_zp,
1451 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1452 const u32 id_cell_a = idvar_a / ndust;
1453 const u32 nvar_loc = idvar_a % ndust;
1454 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1455 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1456 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1457 });
1458 });
1460 shamsys::instance::get_compute_scheduler_ptr(),
1461 sham::DDMultiRef{graphs_zm, interpolators},
1462 sham::DDMultiRef{vel_dust_face_zm.link_fields},
1463 counts_zm,
1464 [ndust](u32 idvar_a, auto link_iter, auto compute, auto acc_link_field) {
1465 const u32 id_cell_a = idvar_a / ndust;
1466 const u32 nvar_loc = idvar_a % ndust;
1467 link_iter.for_each_object_link_id(id_cell_a, [&](u32 id_cell_b, u32 link_id) {
1468 acc_link_field[link_id * ndust + nvar_loc] = compute.get_link_field_val(
1469 id_cell_a * ndust + nvar_loc, id_cell_b * ndust + nvar_loc);
1470 });
1471 });
1472}
1473
1474template<class Tvec, class TgridVec>
1479
utility to manipulate AMR blocks
std::uint32_t u32
32 bit unsigned integer
std::uint64_t u64
64 bit unsigned integer
Shamrock assertion utility.
#define SHAM_ASSERT(x)
Shorthand for SHAM_ASSERT_NAMED without a message.
Definition assert.hpp:67
Class to manage a list of SYCL events.
Definition EventList.hpp:31
Represents a collection of objects distributed across patches identified by a u64 id.
virtual std::string _impl_get_tex() const
get the tex of the node
virtual std::string _impl_get_tex() const
get the tex of the node
virtual std::string _impl_get_tex() const
get the tex of the node
virtual std::string _impl_get_tex() const
get the tex of the node
virtual std::string _impl_get_tex() const
get the tex of the node
Represents a span of data within a PatchDataField.
auto get_read_access(sham::EventList &depends_list, SourceLocation src_loc=SourceLocation{}) const -> details::PatchDataFieldSpan_access_ro_dyn_nvar< T >
Returns a read-only accessor to the data in the span.
void complete_event_state(sycl::event e) const
Completes the event state of the underlying buffer.
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.
A variant of sham::MultiRef for distributed data.
Axis-Aligned bounding box.
Definition AABB.hpp:99
T get_center() const noexcept
Returns the center of the AABB.
Definition AABB.hpp:174
AABB get_intersect(AABB other) const noexcept
Compute the intersection of two AABB.
Definition AABB.hpp:234
utility class to handle AMR blocks
Definition AMRBlock.hpp:35