Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
SourceStep.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 "shammodels/zeus/modules/SourceStep.hpp"
18#include "shammodels/zeus/modules/FaceFlagger.hpp"
19#include "shammodels/zeus/modules/ValueLoader.hpp"
20#include "shamrock/scheduler/InterfacesUtility.hpp"
21#include "shamrock/scheduler/SchedulerUtility.hpp"
22#include "shamsys/legacy/log.hpp"
23
24template<class Tvec, class TgridVec>
25void shammodels::zeus::modules::SourceStep<Tvec, TgridVec>::compute_forces() {
26 StackEntry stack_loc{};
27
28 using namespace shamrock::patch;
29 using namespace shamrock;
30 using namespace shammath;
31 using MergedPDat = shamrock::MergedPatchData;
32
33 using Flagger = FaceFlagger<Tvec, TgridVec>;
34
35 using Block = typename Config::AMRBlock;
36
37 ComputeField<Tscal> &rho_xm = storage.rho_n_xm.get();
38 ComputeField<Tscal> &rho_ym = storage.rho_n_ym.get();
39 ComputeField<Tscal> &rho_zm = storage.rho_n_zm.get();
40
41 ComputeField<Tscal> &pressure_field = storage.pressure.get();
42 ComputeField<Tscal> &p_xm = storage.pres_n_xm.get();
43 ComputeField<Tscal> &p_ym = storage.pres_n_ym.get();
44 ComputeField<Tscal> &p_zm = storage.pres_n_zm.get();
45
46 shamrock::SchedulerUtility utility(scheduler());
47 storage.forces.set(utility.make_compute_field<Tvec>("forces", Block::block_size, [&](u64 id) {
48 return storage.merged_patchdata_ghost.get().get(id).total_elements;
49 }));
50
51 shamrock::patch::PatchDataLayerLayout &ghost_layout
52 = shambase::get_check_ref(storage.ghost_layout.get());
53 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
54
55 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
56 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
57
58 sham::DeviceBuffer<TgridVec> &buf_cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
59 sham::DeviceBuffer<TgridVec> &buf_cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
60
61 Tscal coord_conv_fact = solver_config.grid_coord_to_pos_fact / Block::Nside;
62
63 sham::DeviceBuffer<Tscal> &buf_p = pressure_field.get_buf_check(p.id_patch);
64 sham::DeviceBuffer<Tscal> &buf_rho = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
65
66 sham::DeviceBuffer<Tscal> &buf_rho_xm = rho_xm.get_buf_check(p.id_patch);
67 sham::DeviceBuffer<Tscal> &buf_rho_ym = rho_ym.get_buf_check(p.id_patch);
68 sham::DeviceBuffer<Tscal> &buf_rho_zm = rho_zm.get_buf_check(p.id_patch);
69
70 sham::DeviceBuffer<Tscal> &buf_p_xm = p_xm.get_buf_check(p.id_patch);
71 sham::DeviceBuffer<Tscal> &buf_p_ym = p_ym.get_buf_check(p.id_patch);
72 sham::DeviceBuffer<Tscal> &buf_p_zm = p_zm.get_buf_check(p.id_patch);
73
74 sham::DeviceBuffer<Tvec> &forces_buf = storage.forces.get().get_buf_check(p.id_patch);
75
76 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
77
78 sham::EventList depends_list;
79
80 auto cell_min = buf_cell_min.get_read_access(depends_list);
81 auto cell_max = buf_cell_max.get_read_access(depends_list);
82 auto rho = buf_rho.get_read_access(depends_list);
83 auto rho_xm = buf_rho_xm.get_read_access(depends_list);
84 auto rho_ym = buf_rho_ym.get_read_access(depends_list);
85 auto rho_zm = buf_rho_zm.get_read_access(depends_list);
86 auto press = buf_p.get_read_access(depends_list);
87 auto p_xm = buf_p_xm.get_read_access(depends_list);
88 auto p_ym = buf_p_ym.get_read_access(depends_list);
89 auto p_zm = buf_p_zm.get_read_access(depends_list);
90
91 auto grad_p = forces_buf.get_write_access(depends_list);
92
93 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
94 shambase::parallel_for(
95 cgh, mpdat.total_elements * Block::block_size, "compute grad p", [=](u64 id_a) {
96 u32 block_id = id_a / Block::block_size;
97 Tvec d_cell
98 = (cell_max[block_id] - cell_min[block_id]).template convert<Tscal>()
99 * coord_conv_fact;
100
101 // clang-format off
102 Tscal rho_i_j_k = rho[id_a];
103 Tscal rho_im1_j_k = rho_xm[id_a];
104 Tscal rho_i_jm1_k = rho_ym[id_a];
105 Tscal rho_i_j_km1 = rho_zm[id_a];
106
107 Tscal p_i_j_k = press[id_a];
108 Tscal p_im1_j_k = p_xm[id_a];
109 Tscal p_i_jm1_k = p_ym[id_a];
110 Tscal p_i_j_km1 = p_zm[id_a];
111
112 Tvec dp = {
113 p_i_j_k - p_im1_j_k,
114 p_i_j_k - p_i_jm1_k,
115 p_i_j_k - p_i_j_km1
116 };
117
118 Tvec avg_rho =
119 Tvec{
120 rho_i_j_k + rho_im1_j_k,
121 rho_i_j_k + rho_i_jm1_k,
122 rho_i_j_k + rho_i_j_km1
123 } * Tscal{0.5};
124
125 Tvec grad_p_source_term = dp / (avg_rho * d_cell);
126
127 //sycl::ext::oneapi::experimental::printf("(%f %f %f) (%f %f %f) (%f %f %f)\n"
128 // , dp.x(),dp.y(),dp.z()
129 // , avg_rho.x(),avg_rho.y(),avg_rho.z()
130 // , grad_p_source_term.x(),grad_p_source_term.y(),grad_p_source_term.z()
131 // );
132 //grad_p[id_a] = grad_p_source_term;
133 grad_p[id_a] = -grad_p_source_term;
134 // clang-format on
135 });
136 });
137
138 buf_cell_min.complete_event_state(e);
139 buf_cell_max.complete_event_state(e);
140 buf_rho.complete_event_state(e);
141 buf_rho_xm.complete_event_state(e);
142 buf_rho_ym.complete_event_state(e);
143 buf_rho_zm.complete_event_state(e);
144 buf_p.complete_event_state(e);
145 buf_p_xm.complete_event_state(e);
146 buf_p_ym.complete_event_state(e);
147 buf_p_zm.complete_event_state(e);
148
149 forces_buf.complete_event_state(e);
150 });
151
152 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
153 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
154
155 sham::DeviceBuffer<TgridVec> &buf_cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
156 sham::DeviceBuffer<TgridVec> &buf_cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
157
158 sham::DeviceBuffer<Tvec> &forces_buf = storage.forces.get().get_buf_check(p.id_patch);
159
160 sham::DeviceBuffer<Tscal> &buf_rho = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
161
162 Tscal coord_conv_fact = solver_config.grid_coord_to_pos_fact;
163
164 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
165
166 sham::EventList depends_list;
167 auto cell_min = buf_cell_min.get_read_access(depends_list);
168 auto cell_max = buf_cell_max.get_read_access(depends_list);
169 auto forces = forces_buf.get_write_access(depends_list);
170 auto rho = buf_rho.get_read_access(depends_list);
171
172 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
173 shambase::parallel_for(cgh, pdat.get_obj_cnt(), "add ext force", [=](u64 id_a) {
174 Tvec block_min = cell_min[id_a].template convert<Tscal>();
175 Tvec block_max = cell_max[id_a].template convert<Tscal>();
176 Tvec delta_cell = (block_max - block_min) / Block::side_size;
177 Tvec delta_cell_h = delta_cell * Tscal(0.5);
178
179 Block::for_each_cell_in_block(delta_cell, [=](u32 lid, Tvec delta) {
180 auto get_ext_force = [](Tvec r) {
181 Tscal d = sycl::length(r);
182 return r / (d * d * d + 1e-5);
183 };
184
185 // forces[id_a * Block::block_size + lid] +=
186 // get_ext_force(block_min + delta + delta_cell_h);
187 });
188 });
189 });
190
191 buf_cell_min.complete_event_state(e);
192 buf_cell_max.complete_event_state(e);
193 forces_buf.complete_event_state(e);
194 buf_rho.complete_event_state(e);
195
196 if (storage.forces.get().get_field(p.id_patch).has_nan()) {
197 logger::err_ln("[Zeus]", "nan detected in forces");
198 throw shambase::make_except_with_loc<std::runtime_error>("detected nan");
199 }
200 // logger::raw_ln(storage.forces.get().get_field(p.id_patch).compute_max());
201 });
202}
203
204template<class Tvec, class TgridVec>
205void shammodels::zeus::modules::SourceStep<Tvec, TgridVec>::apply_force(Tscal dt) {
206 StackEntry stack_loc{};
207
208 using namespace shamrock::patch;
209 using namespace shamrock;
210 using namespace shammath;
211 using MergedPDat = shamrock::MergedPatchData;
212
213 using Flagger = FaceFlagger<Tvec, TgridVec>;
214
215 using Block = typename Config::AMRBlock;
216
217 shamrock::patch::PatchDataLayerLayout &ghost_layout
218 = shambase::get_check_ref(storage.ghost_layout.get());
219 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
220
221 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
222 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
223
224 sham::DeviceBuffer<Tvec> &forces_buf = storage.forces.get().get_buf_check(p.id_patch);
225 sham::DeviceBuffer<Tvec> &vel_buf = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
226
227 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
228
229 sham::EventList depends_list;
230
231 auto forces = forces_buf.get_read_access(depends_list);
232 auto vel = vel_buf.get_write_access(depends_list);
233
234 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
235 shambase::parallel_for(
236 cgh, mpdat.total_elements * Block::block_size, "add ext force", [=](u64 id_a) {
237 vel[id_a] += dt * forces[id_a];
238 });
239 });
240
241 forces_buf.complete_event_state(e);
242 vel_buf.complete_event_state(e);
243
244 // logger::raw_ln(storage.forces.get().get_field(p.id_patch).compute_max());
245 });
246
247 storage.forces.reset();
248}
249
250template<class Tvec, class TgridVec>
251void shammodels::zeus::modules::SourceStep<Tvec, TgridVec>::compute_AV() {
252 StackEntry stack_loc{};
253
254 using namespace shamrock::patch;
255 using namespace shamrock;
256 using namespace shammath;
257 using MergedPDat = shamrock::MergedPatchData;
258
259 using Block = typename Config::AMRBlock;
260
261 ComputeField<Tvec> &vel_n = storage.vel_n.get();
262 ComputeField<Tvec> &vel_n_xp = storage.vel_n_xp.get();
263 ComputeField<Tvec> &vel_n_yp = storage.vel_n_yp.get();
264 ComputeField<Tvec> &vel_n_zp = storage.vel_n_zp.get();
265
266 shamrock::SchedulerUtility utility(scheduler());
267 storage.q_AV.set(utility.make_compute_field<Tvec>("q_AV", Block::block_size, [&](u64 id) {
268 return storage.merged_patchdata_ghost.get().get(id).total_elements;
269 }));
270
271 shamrock::patch::PatchDataLayerLayout &ghost_layout
272 = shambase::get_check_ref(storage.ghost_layout.get());
273 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
274 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
275
276 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
277 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
278
279 sham::DeviceBuffer<TgridVec> &buf_cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
280 sham::DeviceBuffer<TgridVec> &buf_cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
281
282 Tscal coord_conv_fact = solver_config.grid_coord_to_pos_fact / Block::Nside;
283
284 sham::DeviceBuffer<Tscal> &buf_rho = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
285
286 sham::DeviceBuffer<Tvec> &buf_vel = vel_n.get_buf_check(p.id_patch);
287 sham::DeviceBuffer<Tvec> &buf_vel_xp = vel_n_xp.get_buf_check(p.id_patch);
288 sham::DeviceBuffer<Tvec> &buf_vel_yp = vel_n_yp.get_buf_check(p.id_patch);
289 sham::DeviceBuffer<Tvec> &buf_vel_zp = vel_n_zp.get_buf_check(p.id_patch);
290
291 sham::DeviceBuffer<Tvec> &q_AV_buf = storage.q_AV.get().get_buf_check(p.id_patch);
292
293 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
294
295 sham::EventList depends_list;
296
297 auto cell_min = buf_cell_min.get_read_access(depends_list);
298 auto cell_max = buf_cell_max.get_read_access(depends_list);
299 auto rho = buf_rho.get_read_access(depends_list);
300 auto vel = buf_vel.get_read_access(depends_list);
301 auto vel_xp = buf_vel_xp.get_read_access(depends_list);
302 auto vel_yp = buf_vel_yp.get_read_access(depends_list);
303 auto vel_zp = buf_vel_zp.get_read_access(depends_list);
304 auto q_AV = q_AV_buf.get_write_access(depends_list);
305
306 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
307 shambase::parallel_for(
308 cgh, mpdat.total_elements * Block::block_size, "compute AV", [=](u64 id_a) {
309 u32 block_id = id_a / Block::block_size;
310 Tvec d_cell
311 = (cell_max[block_id] - cell_min[block_id]).template convert<Tscal>()
312 * coord_conv_fact;
313
314 // clang-format off
315 Tscal rho_i_j_k = rho[id_a];
316
317 Tvec vel_i_j_k = vel[id_a];
318 Tvec vel_ip1_j_k = vel_xp[id_a];
319 Tvec vel_i_jp1_k = vel_yp[id_a];
320 Tvec vel_i_j_kp1 = vel_zp[id_a];
321
322 Tvec dv = {
323 vel_ip1_j_k.x() - vel_i_j_k.x(),
324 vel_i_jp1_k.y() - vel_i_j_k.y(),
325 vel_i_j_kp1.z() - vel_i_j_k.z()
326 };
327
328 dv = sham::negative_part(dv);
329
330 constexpr Tscal C2 = 3;
331
332 q_AV[id_a] = C2*rho_i_j_k*(dv*dv);
333 // clang-format on
334 });
335 });
336
337 buf_cell_min.complete_event_state(e);
338 buf_cell_max.complete_event_state(e);
339 buf_rho.complete_event_state(e);
340 buf_vel.complete_event_state(e);
341 buf_vel_xp.complete_event_state(e);
342 buf_vel_yp.complete_event_state(e);
343 buf_vel_zp.complete_event_state(e);
344 q_AV_buf.complete_event_state(e);
345 });
346}
347
348template<class Tvec, class TgridVec>
349void shammodels::zeus::modules::SourceStep<Tvec, TgridVec>::apply_AV(Tscal dt) {
350 StackEntry stack_loc{};
351
352 using namespace shamrock::patch;
353 using namespace shamrock;
354 using namespace shammath;
355 using MergedPDat = shamrock::MergedPatchData;
356
357 using Block = typename Config::AMRBlock;
358
359 ComputeField<Tvec> &q_AV_n = storage.q_AV.get();
360 ComputeField<Tvec> &q_AV_n_xm = storage.q_AV_n_xm.get();
361 ComputeField<Tvec> &q_AV_n_ym = storage.q_AV_n_ym.get();
362 ComputeField<Tvec> &q_AV_n_zm = storage.q_AV_n_zm.get();
363
364 ComputeField<Tscal> &rho_xm = storage.rho_n_xm.get();
365 ComputeField<Tscal> &rho_ym = storage.rho_n_ym.get();
366 ComputeField<Tscal> &rho_zm = storage.rho_n_zm.get();
367
368 ComputeField<Tvec> &vel_n = storage.vel_n.get();
369 ComputeField<Tvec> &vel_n_xp = storage.vel_n_xp.get();
370 ComputeField<Tvec> &vel_n_yp = storage.vel_n_yp.get();
371 ComputeField<Tvec> &vel_n_zp = storage.vel_n_zp.get();
372
373 shamrock::patch::PatchDataLayerLayout &ghost_layout
374 = shambase::get_check_ref(storage.ghost_layout.get());
375 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
376 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
377 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
378
379 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
380 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
381
382 sham::DeviceBuffer<TgridVec> &buf_cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
383 sham::DeviceBuffer<TgridVec> &buf_cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
384
385 Tscal coord_conv_fact = solver_config.grid_coord_to_pos_fact / Block::Nside;
386
387 sham::DeviceBuffer<Tscal> &buf_rho = mpdat.pdat.get_field_buf_ref<Tscal>(irho_interf);
388 sham::DeviceBuffer<Tvec> &buf_vel = mpdat.pdat.get_field_buf_ref<Tvec>(ivel_interf);
389
390 sham::DeviceBuffer<Tvec> &buf_vel_xp = vel_n_xp.get_buf_check(p.id_patch);
391 sham::DeviceBuffer<Tvec> &buf_vel_yp = vel_n_yp.get_buf_check(p.id_patch);
392 sham::DeviceBuffer<Tvec> &buf_vel_zp = vel_n_zp.get_buf_check(p.id_patch);
393
394 sham::DeviceBuffer<Tscal> &buf_rho_xm = rho_xm.get_buf_check(p.id_patch);
395 sham::DeviceBuffer<Tscal> &buf_rho_ym = rho_ym.get_buf_check(p.id_patch);
396 sham::DeviceBuffer<Tscal> &buf_rho_zm = rho_zm.get_buf_check(p.id_patch);
397
398 sham::DeviceBuffer<Tvec> &q_AV_buf = storage.q_AV.get().get_buf_check(p.id_patch);
399 sham::DeviceBuffer<Tvec> &buf_q_AV_n_xm = q_AV_n_xm.get_buf_check(p.id_patch);
400 sham::DeviceBuffer<Tvec> &buf_q_AV_n_ym = q_AV_n_ym.get_buf_check(p.id_patch);
401 sham::DeviceBuffer<Tvec> &buf_q_AV_n_zm = q_AV_n_zm.get_buf_check(p.id_patch);
402
403 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
404
405 sham::EventList depends_list;
406
407 auto cell_min = buf_cell_min.get_read_access(depends_list);
408 auto cell_max = buf_cell_max.get_read_access(depends_list);
409 auto rho = buf_rho.get_read_access(depends_list);
410 auto rho_xm = buf_rho_xm.get_read_access(depends_list);
411 auto rho_ym = buf_rho_ym.get_read_access(depends_list);
412 auto rho_zm = buf_rho_zm.get_read_access(depends_list);
413 auto q_AV = q_AV_buf.get_read_access(depends_list);
414 auto q_AV_xm = buf_q_AV_n_xm.get_read_access(depends_list);
415 auto q_AV_ym = buf_q_AV_n_ym.get_read_access(depends_list);
416 auto q_AV_zm = buf_q_AV_n_zm.get_read_access(depends_list);
417 auto vel = buf_vel.get_write_access(depends_list);
418
419 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
420 shambase::parallel_for(
421 cgh, mpdat.total_elements * Block::block_size, "add vel AV", [=](u64 id_a) {
422 u32 block_id = id_a / Block::block_size;
423 Tvec d_cell
424 = (cell_max[block_id] - cell_min[block_id]).template convert<Tscal>()
425 * coord_conv_fact;
426
427 // clang-format off
428 Tscal rho_i_j_k = rho[id_a];
429 Tscal rho_im1_j_k = rho_xm[id_a];
430 Tscal rho_i_jm1_k = rho_ym[id_a];
431 Tscal rho_i_j_km1 = rho_zm[id_a];
432
433 Tvec q_i_j_k = q_AV[id_a];
434 Tvec q_im1_j_k = q_AV_xm[id_a];
435 Tvec q_i_jm1_k = q_AV_ym[id_a];
436 Tvec q_i_j_km1 = q_AV_zm[id_a];
437
438 Tvec avg_rho =
439 Tvec{
440 rho_i_j_k + rho_im1_j_k,
441 rho_i_j_k + rho_i_jm1_k,
442 rho_i_j_k + rho_i_j_km1
443 } * Tscal{0.5};
444
445 Tvec dq = {
446 q_i_j_k.x()-q_im1_j_k.x(),
447 q_i_j_k.y()-q_i_jm1_k.y(),
448 q_i_j_k.z()-q_i_j_km1.z()
449 };
450
451 vel[id_a] += - dt*(dq)/ (avg_rho * d_cell);
452 // clang-format on
453 });
454 });
455
456 buf_cell_min.complete_event_state(e);
457 buf_cell_max.complete_event_state(e);
458 buf_rho.complete_event_state(e);
459 buf_rho_xm.complete_event_state(e);
460 buf_rho_ym.complete_event_state(e);
461 buf_rho_zm.complete_event_state(e);
462 q_AV_buf.complete_event_state(e);
463 buf_q_AV_n_xm.complete_event_state(e);
464 buf_q_AV_n_ym.complete_event_state(e);
465 buf_q_AV_n_zm.complete_event_state(e);
466 buf_vel.complete_event_state(e);
467 });
468
469 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
470 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
471
472 sham::DeviceBuffer<TgridVec> &buf_cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
473 sham::DeviceBuffer<TgridVec> &buf_cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
474
475 Tscal coord_conv_fact = solver_config.grid_coord_to_pos_fact / Block::Nside;
476
477 sham::DeviceBuffer<Tscal> &buf_eint = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
478
479 sham::DeviceBuffer<Tvec> &buf_vel = vel_n.get_buf_check(p.id_patch);
480 sham::DeviceBuffer<Tvec> &buf_vel_xp = vel_n_xp.get_buf_check(p.id_patch);
481 sham::DeviceBuffer<Tvec> &buf_vel_yp = vel_n_yp.get_buf_check(p.id_patch);
482 sham::DeviceBuffer<Tvec> &buf_vel_zp = vel_n_zp.get_buf_check(p.id_patch);
483
484 sham::DeviceBuffer<Tvec> &q_AV_buf = storage.q_AV.get().get_buf_check(p.id_patch);
485 sham::DeviceBuffer<Tvec> &buf_q_AV_n_xm = q_AV_n_xm.get_buf_check(p.id_patch);
486 sham::DeviceBuffer<Tvec> &buf_q_AV_n_ym = q_AV_n_ym.get_buf_check(p.id_patch);
487 sham::DeviceBuffer<Tvec> &buf_q_AV_n_zm = q_AV_n_zm.get_buf_check(p.id_patch);
488
489 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
490
491 sham::EventList depends_list;
492 auto cell_min = buf_cell_min.get_read_access(depends_list);
493 auto cell_max = buf_cell_max.get_read_access(depends_list);
494 auto vel = buf_vel.get_read_access(depends_list);
495 auto vel_xp = buf_vel_xp.get_read_access(depends_list);
496 auto vel_yp = buf_vel_yp.get_read_access(depends_list);
497 auto vel_zp = buf_vel_zp.get_read_access(depends_list);
498 auto q_AV = q_AV_buf.get_read_access(depends_list);
499 auto q_AV_xm = buf_q_AV_n_xm.get_read_access(depends_list);
500 auto q_AV_ym = buf_q_AV_n_ym.get_read_access(depends_list);
501 auto q_AV_zm = buf_q_AV_n_zm.get_read_access(depends_list);
502 auto eint = buf_eint.get_write_access(depends_list);
503
504 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
505 shambase::parallel_for(
506 cgh, pdat.get_obj_cnt() * Block::block_size, "add eint AV", [=](u64 id_a) {
507 u32 block_id = id_a / Block::block_size;
508 Tvec d_cell
509 = (cell_max[block_id] - cell_min[block_id]).template convert<Tscal>()
510 * coord_conv_fact;
511
512 // clang-format off
513 Tvec vel_i_j_k = vel[id_a];
514 Tvec vel_ip1_j_k = vel_xp[id_a];
515 Tvec vel_i_jp1_k = vel_yp[id_a];
516 Tvec vel_i_j_kp1 = vel_zp[id_a];
517
518 Tvec q_i_j_k = q_AV[id_a];
519
520 Tvec dv = {
521 vel_ip1_j_k.x() - vel_i_j_k.x(),
522 vel_i_jp1_k.y() - vel_i_j_k.y(),
523 vel_i_j_kp1.z() - vel_i_j_k.z()
524 };
525
526 eint[id_a] += -dt*sycl::dot(q_i_j_k,dv/ d_cell);
527 // clang-format on
528 });
529 });
530
531 buf_cell_min.complete_event_state(e);
532 buf_cell_max.complete_event_state(e);
533 buf_vel.complete_event_state(e);
534 buf_vel_xp.complete_event_state(e);
535 buf_vel_yp.complete_event_state(e);
536 buf_vel_zp.complete_event_state(e);
537 q_AV_buf.complete_event_state(e);
538 buf_q_AV_n_xm.complete_event_state(e);
539 buf_q_AV_n_ym.complete_event_state(e);
540 buf_q_AV_n_zm.complete_event_state(e);
541 buf_eint.complete_event_state(e);
542 });
543}
544
545template<class Tvec, class TgridVec>
546void shammodels::zeus::modules::SourceStep<Tvec, TgridVec>::compute_div_v() {
547 StackEntry stack_loc{};
548 using namespace shamrock::patch;
549 using namespace shamrock;
550 using namespace shammath;
551 using MergedPDat = shamrock::MergedPatchData;
552
553 using Flagger = FaceFlagger<Tvec, TgridVec>;
554
555 using Block = typename Config::AMRBlock;
556
557 shamrock::SchedulerUtility utility(scheduler());
558 storage.div_v_n.set(
559 utility.make_compute_field<Tscal>("div_v_n", Block::block_size, [&](u64 id) {
560 return storage.merged_patchdata_ghost.get().get(id).total_elements;
561 }));
562
563 ComputeField<Tvec> &vel_n = storage.vel_n.get();
564 ComputeField<Tvec> &vel_n_xp = storage.vel_n_xp.get();
565 ComputeField<Tvec> &vel_n_yp = storage.vel_n_yp.get();
566 ComputeField<Tvec> &vel_n_zp = storage.vel_n_zp.get();
567
568 ComputeField<Tscal> &div_v = storage.div_v_n.get();
569
570 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
571 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
572
573 sham::DeviceBuffer<TgridVec> &buf_cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
574 sham::DeviceBuffer<TgridVec> &buf_cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
575
576 Tscal coord_conv_fact = solver_config.grid_coord_to_pos_fact / Block::Nside;
577
578 sham::DeviceBuffer<Tvec> &buf_vel = vel_n.get_buf_check(p.id_patch);
579 sham::DeviceBuffer<Tvec> &buf_vel_xp = vel_n_xp.get_buf_check(p.id_patch);
580 sham::DeviceBuffer<Tvec> &buf_vel_yp = vel_n_yp.get_buf_check(p.id_patch);
581 sham::DeviceBuffer<Tvec> &buf_vel_zp = vel_n_zp.get_buf_check(p.id_patch);
582
583 sham::DeviceBuffer<Tscal> &buf_div_v = div_v.get_buf_check(p.id_patch);
584
585 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
586
587 sham::EventList depends_list;
588 auto cell_min = buf_cell_min.get_read_access(depends_list);
589 auto cell_max = buf_cell_max.get_read_access(depends_list);
590
591 auto vel = buf_vel.get_read_access(depends_list);
592 auto vel_xp = buf_vel_xp.get_read_access(depends_list);
593 auto vel_yp = buf_vel_yp.get_read_access(depends_list);
594 auto vel_zp = buf_vel_zp.get_read_access(depends_list);
595
596 auto divv = buf_div_v.get_write_access(depends_list);
597
598 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
599 shambase::parallel_for(
600 cgh, pdat.get_obj_cnt() * Block::block_size, "compute divv", [=](u64 id_a) {
601 u32 block_id = id_a / Block::block_size;
602 Tvec d_cell
603 = (cell_max[block_id] - cell_min[block_id]).template convert<Tscal>()
604 * coord_conv_fact;
605
606 // clang-format off
607 Tvec vel_i_j_k = vel[id_a];
608 Tvec vel_ip1_j_k = vel_xp[id_a];
609 Tvec vel_i_jp1_k = vel_yp[id_a];
610 Tvec vel_i_j_kp1 = vel_zp[id_a];
611
612 Tvec dv = {
613 vel_ip1_j_k.x() - vel_i_j_k.x(),
614 vel_i_jp1_k.y() - vel_i_j_k.y(),
615 vel_i_j_kp1.z() - vel_i_j_k.z()
616 };
617
618 divv[id_a] += sycl::dot(dv,Tvec{1,1,1}/ d_cell);
619 // clang-format on
620 });
621 });
622
623 buf_cell_min.complete_event_state(e);
624 buf_cell_max.complete_event_state(e);
625 buf_vel.complete_event_state(e);
626 buf_vel_xp.complete_event_state(e);
627 buf_vel_yp.complete_event_state(e);
628 buf_vel_zp.complete_event_state(e);
629 buf_div_v.complete_event_state(e);
630 });
631}
632
633template<class Tvec, class TgridVec>
634void shammodels::zeus::modules::SourceStep<Tvec, TgridVec>::update_eint_eos(Tscal dt) {
635 StackEntry stack_loc{};
636
637 using namespace shamrock::patch;
638 using namespace shamrock;
639 using namespace shammath;
640 using MergedPDat = shamrock::MergedPatchData;
641
642 using Block = typename Config::AMRBlock;
643
644 ComputeField<Tscal> &div_v = storage.div_v_n.get();
645
646 shamrock::patch::PatchDataLayerLayout &ghost_layout
647 = shambase::get_check_ref(storage.ghost_layout.get());
648 u32 irho_interf = ghost_layout.get_field_idx<Tscal>("rho");
649 u32 ieint_interf = ghost_layout.get_field_idx<Tscal>("eint");
650 u32 ivel_interf = ghost_layout.get_field_idx<Tvec>("vel");
651
652 scheduler().for_each_patchdata_nonempty([&](Patch p, PatchDataLayer &pdat) {
653 MergedPDat &mpdat = storage.merged_patchdata_ghost.get().get(p.id_patch);
654
655 sham::DeviceBuffer<TgridVec> &buf_cell_min = mpdat.pdat.get_field_buf_ref<TgridVec>(0);
656 sham::DeviceBuffer<TgridVec> &buf_cell_max = mpdat.pdat.get_field_buf_ref<TgridVec>(1);
657
658 Tscal coord_conv_fact = solver_config.grid_coord_to_pos_fact / Block::Nside;
659
660 sham::DeviceBuffer<Tscal> &buf_eint = mpdat.pdat.get_field_buf_ref<Tscal>(ieint_interf);
661
662 sham::DeviceBuffer<Tscal> &buf_divv = div_v.get_buf_check(p.id_patch);
663
664 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
665
666 sham::EventList depends_list;
667 auto cell_min = buf_cell_min.get_read_access(depends_list);
668 auto cell_max = buf_cell_max.get_read_access(depends_list);
669
670 auto divv = buf_divv.get_read_access(depends_list);
671 auto eint = buf_eint.get_write_access(depends_list);
672
673 auto e = q.submit(depends_list, [&](sycl::handler &cgh) {
674 Tscal fact = (dt / 2.) * (solver_config.eos_gamma - 1);
675
676 shambase::parallel_for(
677 cgh, pdat.get_obj_cnt() * Block::block_size, "evolve eint", [=](u64 id_a) {
678 u32 block_id = id_a / Block::block_size;
679 Tvec d_cell
680 = (cell_max[block_id] - cell_min[block_id]).template convert<Tscal>()
681 * coord_conv_fact;
682
683 Tscal factdivv = divv[id_a] * fact;
684
685 eint[id_a] *= (1 - factdivv) / (1 + factdivv);
686 });
687 });
688
689 buf_cell_min.complete_event_state(e);
690 buf_cell_max.complete_event_state(e);
691 buf_divv.complete_event_state(e);
692 buf_eint.complete_event_state(e);
693 });
694}
695
696template class shammodels::zeus::modules::SourceStep<f64_3, i64_3>;
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::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
shammodels::zeus::modules::FaceFlagger
flag faces with a lookup index for the orientation
Definition FaceFlagger.hpp:39
shammodels::zeus::modules::SourceStep::compute_forces
void compute_forces()
compute general forces (pressure + external and store them into SolverStorage::forces)
Definition SourceStep.cpp:25
shammodels::zeus::modules::SourceStep::apply_force
void apply_force(Tscal dt)
apply the generalized forces
Definition SourceStep.cpp:205
shammodels::zeus::modules::SourceStep::compute_AV
void compute_AV()
Compute the values of the artificial viscosity terms ( equations 33,34)
Definition SourceStep.cpp:251
shamrock::SchedulerUtility::make_compute_field
ComputeField< T > make_compute_field(std::string new_name, u32 nvar)
create a compute field and init it to zeros
Definition SchedulerUtility.hpp:227
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
shambase::throw_with_loc
void throw_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Throw an exception and append the source location to it.
Definition exception.hpp:132
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
shammath
namespace for math utility
Definition AABB.hpp:26
shamrock
namespace for the main framework
Definition __init__.py:1
shammodels::amr::AMRBlock
utility class to handle AMR blocks
Definition AMRBlock.hpp:35
shamrock::patch::Patch
Patch object that contain generic patch information.
Definition Patch.hpp:33
Generated on Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8