Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
CartesianRender.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
20#include "shammath/AABB.hpp"
31
33
34 template<class Tvec>
35 sham::DeviceBuffer<Tvec> pixel_to_positions(
36 Tvec center, Tvec delta_x, Tvec delta_y, u32 nx, u32 ny) {
37
38 sham::DeviceBuffer<Tvec> ret{nx * ny, shamsys::instance::get_compute_scheduler_ptr()};
39
40 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
41
43 q, sham::MultiRef{}, sham::MultiRef{ret}, nx * ny, [=](u32 gid, Tvec *position) {
44 u32 ix = gid % nx;
45 u32 iy = gid / nx;
46 f64 fx = ((f64(ix) + 0.5) / nx) - 0.5;
47 f64 fy = ((f64(iy) + 0.5) / ny) - 0.5;
48 position[gid] = center + delta_x * fx + delta_y * fy;
49 });
50
51 return ret;
52 }
53
54 template<class Tvec>
55 sham::DeviceBuffer<shammath::Ray<Tvec>> pixel_to_orthographic_rays(
56 Tvec center, Tvec delta_x, Tvec delta_y, u32 nx, u32 ny) {
57
58 using Tscal = shambase::VecComponent<Tvec>;
59
60 sham::DeviceBuffer<shammath::Ray<Tvec>> ret{
61 nx * ny, shamsys::instance::get_compute_scheduler_ptr()};
62
63 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().get_queue();
64
65 Tvec e_z = sycl::cross(delta_x, delta_y);
66 Tscal len = sycl::length(e_z);
67 if (!(len > 0)) {
69 "The cross product of delta_x and delta_y is zero\n"
70 " args :"
71 " center = {}\n"
72 " delta_x = {}\n"
73 " delta_y = {}\n"
74 " nx = {}\n"
75 " ny = {}\n"
76 " -> e_z = {}\n",
77 center,
78 delta_x,
79 delta_y,
80 nx,
81 ny,
82 e_z));
83 }
84 e_z /= len;
85
87 q,
88 sham::MultiRef{},
89 sham::MultiRef{ret},
90 nx * ny,
91 [=](u32 gid, shammath::Ray<Tvec> *ray) {
92 u32 ix = gid % nx;
93 u32 iy = gid / nx;
94 f64 fx = ((f64(ix) + 0.5) / nx) - 0.5;
95 f64 fy = ((f64(iy) + 0.5) / ny) - 0.5;
96 Tvec pos_render = center + delta_x * fx + delta_y * fy;
97
98 ray[gid] = shammath::Ray<Tvec>(pos_render, e_z);
99 });
100
101 return ret;
102 }
103
104 template<class Tvec, class Tfield, template<class> class SPHKernel>
105 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_slice(
106 std::string field_name,
107 const sham::DeviceBuffer<Tvec> &positions,
108 std::optional<std::function<py::array_t<Tfield>(size_t, pybind11::dict &)>> custom_getter)
109 -> sham::DeviceBuffer<Tfield> {
110
111 if (shamcomm::world_rank() == 0) {
112 logger::info_ln(
113 "sph::CartesianRender",
114 shambase::format(
115 "compute_slice field_name: {}, positions count: {}",
116 field_name,
117 positions.get_size()));
118 }
119
120 shambase::Timer t;
121 t.start();
122
123 auto ret = RenderFieldGetter<Tvec, Tfield, SPHKernel>(context, solver_config, storage)
124 .runner_function(
125 field_name,
126 [&](auto field_getter) -> sham::DeviceBuffer<Tfield> {
127 return compute_slice(field_getter, positions);
128 },
129 custom_getter);
130
131 t.stop();
132 if (shamcomm::world_rank() == 0) {
133 logger::info_ln(
134 "sph::CartesianRender",
135 shambase::format("compute_slice took {}", t.get_time_str()));
136 }
137
138 return ret;
139 }
140
141 template<class Tvec, class Tfield, template<class> class SPHKernel>
142 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_column_integ(
143 std::string field_name,
144 const sham::DeviceBuffer<shammath::Ray<Tvec>> &rays,
145 std::optional<std::function<py::array_t<Tfield>(size_t, pybind11::dict &)>> custom_getter)
146 -> sham::DeviceBuffer<Tfield> {
147
148 if (shamcomm::world_rank() == 0) {
149 logger::info_ln(
150 "sph::CartesianRender",
151 shambase::format(
152 "compute_column_integ field_name: {}, rays count: {}",
153 field_name,
154 rays.get_size()));
155 }
156
157 shambase::Timer t;
158 t.start();
159
160 auto ret = RenderFieldGetter<Tvec, Tfield, SPHKernel>(context, solver_config, storage)
161 .runner_function(
162 field_name,
163 [&](auto field_getter) -> sham::DeviceBuffer<Tfield> {
164 return compute_column_integ(field_getter, rays);
165 },
166 custom_getter);
167
168 t.stop();
169 if (shamcomm::world_rank() == 0) {
170 logger::info_ln(
171 "sph::CartesianRender",
172 shambase::format("compute_column_integ took {}", t.get_time_str()));
173 }
174
175 return ret;
176 }
177
178 template<class Tvec, class Tfield, template<class> class SPHKernel>
179 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_azymuthal_integ(
180 std::string field_name,
181 const sham::DeviceBuffer<shammath::RingRay<Tvec>> &ring_rays,
182 std::optional<std::function<py::array_t<Tfield>(size_t, pybind11::dict &)>> custom_getter)
183 -> sham::DeviceBuffer<Tfield> {
184
185 if (shamcomm::world_rank() == 0) {
186 logger::info_ln(
187 "sph::CartesianRender",
188 shambase::format(
189 "compute_azymuthal_integ field_name: {}, ring_rays count: {}",
190 field_name,
191 ring_rays.get_size()));
192 }
193
194 shambase::Timer t;
195 t.start();
196
197 auto ret = RenderFieldGetter<Tvec, Tfield, SPHKernel>(context, solver_config, storage)
198 .runner_function(
199 field_name,
200 [&](auto field_getter) -> sham::DeviceBuffer<Tfield> {
201 return compute_azymuthal_integ(field_getter, ring_rays);
202 },
203 custom_getter);
204
205 t.stop();
206 if (shamcomm::world_rank() == 0) {
207 logger::info_ln(
208 "sph::CartesianRender",
209 shambase::format("compute_azymuthal_integ took {}", t.get_time_str()));
210 }
211
212 return ret;
213 }
214
215 template<class Tvec, class Tfield, template<class> class SPHKernel>
216 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_slice(
217 std::function<field_getter_t> field_getter, const sham::DeviceBuffer<Tvec> &positions)
218 -> sham::DeviceBuffer<Tfield> {
219
220 auto part_counts = shamrock::solvergraph::Indexes<u32>::make_shared("part_counts", "N");
221 auto positions_refs
222 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("positions", "\\mathbf{r}");
223 auto hpart_refs = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("h_part", "h");
224 auto field_data
225 = std::make_shared<shamrock::solvergraph::Field<Tfield>>(1, "field_data", "f");
226
227 shamrock::solvergraph::DDPatchDataFieldRef<Tvec> pos_dd;
228 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> h_dd;
229
230 scheduler().for_each_patchdata_nonempty(
231 [&](const shamrock::patch::Patch cur_p, shamrock::patch::PatchDataLayer &pdat) {
232 u64 id = cur_p.id_patch;
233 u32 cnt = pdat.get_obj_cnt();
234
235 part_counts->indexes.add_obj(id, std::move(cnt));
236 pos_dd.add_obj(id, std::ref(pdat.get_field<Tvec>(0)));
237 h_dd.add_obj(
238 id, std::ref(pdat.get_field<Tscal>(pdat.pdl().get_field_idx<Tscal>("hpart"))));
239 });
240
241 positions_refs->set_refs(pos_dd);
242 hpart_refs->set_refs(h_dd);
243
244 field_data->ensure_sizes(part_counts->indexes);
245
246 scheduler().for_each_patchdata_nonempty(
247 [&](const shamrock::patch::Patch cur_p, shamrock::patch::PatchDataLayer &pdat) {
248 const sham::DeviceBuffer<Tfield> &src = field_getter(cur_p, pdat);
249 field_data->get(cur_p.id_patch).overwrite(src, static_cast<u32>(src.get_size()));
250 });
251
252 auto gpart_mass = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("gpart_mass", "m");
253 gpart_mass->data = solver_config.gpart_mass;
254
255 auto tree_reduction_level
256 = shamrock::solvergraph::IDataEdge<u32>::make_shared("tree_reduction_level", "l");
257 tree_reduction_level->data = solver_config.tree_reduction_level;
258
259 auto interp_points = std::make_shared<shamrock::solvergraph::DeviceBufferEdge<Tvec>>(
260 "interp_points", "\\mathbf{q}");
261 interp_points->value.resize(positions.get_size());
262 interp_points->value.copy_from(positions);
263
264 auto interpolated_field = std::make_shared<shamrock::solvergraph::DeviceBufferEdge<Tfield>>(
265 "interpolated_field", "f_{\\rm interp}");
266
267 auto node = std::make_shared<SPHInterpolation<Tvec, Tfield, SPHKernel>>();
268 node->set_edges(
269 gpart_mass,
270 tree_reduction_level,
271 part_counts,
272 positions_refs,
273 hpart_refs,
274 field_data,
275 interp_points,
276 interpolated_field);
277 node->evaluate();
278
279 sham::DeviceBuffer<Tfield> ret{
280 interpolated_field->value.get_size(), shamsys::instance::get_compute_scheduler_ptr()};
281 ret.copy_from(interpolated_field->value);
282
283 return ret;
284 }
285
286 template<class Tvec, class Tfield, template<class> class SPHKernel>
287 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_column_integ(
288 std::function<field_getter_t> field_getter,
289 const sham::DeviceBuffer<shammath::Ray<Tvec>> &rays) -> sham::DeviceBuffer<Tfield> {
290
291 auto part_counts = shamrock::solvergraph::Indexes<u32>::make_shared("part_counts", "N");
292 auto positions_refs
293 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("positions", "\\mathbf{r}");
294 auto hpart_refs = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("h_part", "h");
295 auto field_data
296 = std::make_shared<shamrock::solvergraph::Field<Tfield>>(1, "field_data", "f");
297
298 shamrock::solvergraph::DDPatchDataFieldRef<Tvec> pos_dd;
299 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> h_dd;
300
301 scheduler().for_each_patchdata_nonempty(
302 [&](const shamrock::patch::Patch cur_p, shamrock::patch::PatchDataLayer &pdat) {
303 u64 id = cur_p.id_patch;
304 u32 cnt = pdat.get_obj_cnt();
305
306 part_counts->indexes.add_obj(id, std::move(cnt));
307 pos_dd.add_obj(id, std::ref(pdat.get_field<Tvec>(0)));
308 h_dd.add_obj(
309 id, std::ref(pdat.get_field<Tscal>(pdat.pdl().get_field_idx<Tscal>("hpart"))));
310 });
311
312 positions_refs->set_refs(pos_dd);
313 hpart_refs->set_refs(h_dd);
314
315 field_data->ensure_sizes(part_counts->indexes);
316
317 scheduler().for_each_patchdata_nonempty(
318 [&](const shamrock::patch::Patch cur_p, shamrock::patch::PatchDataLayer &pdat) {
319 const sham::DeviceBuffer<Tfield> &src = field_getter(cur_p, pdat);
320 field_data->get(cur_p.id_patch).overwrite(src, static_cast<u32>(src.get_size()));
321 });
322
323 auto gpart_mass = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("gpart_mass", "m");
324 gpart_mass->data = solver_config.gpart_mass;
325
326 auto tree_reduction_level
327 = shamrock::solvergraph::IDataEdge<u32>::make_shared("tree_reduction_level", "l");
328 tree_reduction_level->data = solver_config.tree_reduction_level;
329
330 auto rays_edge
331 = std::make_shared<shamrock::solvergraph::DeviceBufferEdge<shammath::Ray<Tvec>>>(
332 "rays", "\\mathbf{r}_{\\rm ray}");
333 rays_edge->value.resize(rays.get_size());
334 rays_edge->value.copy_from(rays);
335
336 auto interpolated_field = std::make_shared<shamrock::solvergraph::DeviceBufferEdge<Tfield>>(
337 "interpolated_field", "f_{\\rm interp}");
338
339 auto node = std::make_shared<SPHColumnInteg<Tvec, Tfield, SPHKernel>>();
340 node->set_edges(
341 gpart_mass,
342 tree_reduction_level,
343 part_counts,
344 positions_refs,
345 hpart_refs,
346 field_data,
347 rays_edge,
348 interpolated_field);
349 node->evaluate();
350
351 sham::DeviceBuffer<Tfield> ret{
352 interpolated_field->value.get_size(), shamsys::instance::get_compute_scheduler_ptr()};
353 ret.copy_from(interpolated_field->value);
354
355 return ret;
356 }
357
358 template<class Tvec, class Tfield, template<class> class SPHKernel>
359 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_azymuthal_integ(
360 std::function<field_getter_t> field_getter,
361 const sham::DeviceBuffer<shammath::RingRay<Tvec>> &ring_rays)
362 -> sham::DeviceBuffer<Tfield> {
363
364 auto part_counts = shamrock::solvergraph::Indexes<u32>::make_shared("part_counts", "N");
365 auto positions_refs
366 = std::make_shared<shamrock::solvergraph::FieldRefs<Tvec>>("positions", "\\mathbf{r}");
367 auto hpart_refs = std::make_shared<shamrock::solvergraph::FieldRefs<Tscal>>("h_part", "h");
368 auto field_data
369 = std::make_shared<shamrock::solvergraph::Field<Tfield>>(1, "field_data", "f");
370
371 shamrock::solvergraph::DDPatchDataFieldRef<Tvec> pos_dd;
372 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> h_dd;
373
374 scheduler().for_each_patchdata_nonempty(
375 [&](const shamrock::patch::Patch cur_p, shamrock::patch::PatchDataLayer &pdat) {
376 u64 id = cur_p.id_patch;
377 u32 cnt = pdat.get_obj_cnt();
378
379 part_counts->indexes.add_obj(id, std::move(cnt));
380 pos_dd.add_obj(id, std::ref(pdat.get_field<Tvec>(0)));
381 h_dd.add_obj(
382 id, std::ref(pdat.get_field<Tscal>(pdat.pdl().get_field_idx<Tscal>("hpart"))));
383 });
384
385 positions_refs->set_refs(pos_dd);
386 hpart_refs->set_refs(h_dd);
387
388 field_data->ensure_sizes(part_counts->indexes);
389
390 scheduler().for_each_patchdata_nonempty(
391 [&](const shamrock::patch::Patch cur_p, shamrock::patch::PatchDataLayer &pdat) {
392 const sham::DeviceBuffer<Tfield> &src = field_getter(cur_p, pdat);
393 field_data->get(cur_p.id_patch).overwrite(src, static_cast<u32>(src.get_size()));
394 });
395
396 auto gpart_mass = shamrock::solvergraph::IDataEdge<Tscal>::make_shared("gpart_mass", "m");
397 gpart_mass->data = solver_config.gpart_mass;
398
399 auto tree_reduction_level
400 = shamrock::solvergraph::IDataEdge<u32>::make_shared("tree_reduction_level", "l");
401 tree_reduction_level->data = solver_config.tree_reduction_level;
402
403 auto ring_rays_edge
404 = std::make_shared<shamrock::solvergraph::DeviceBufferEdge<shammath::RingRay<Tvec>>>(
405 "ring_rays", "\\mathbf{r}_{\\rm ring}");
406 ring_rays_edge->value.resize(ring_rays.get_size());
407 ring_rays_edge->value.copy_from(ring_rays);
408
409 auto interpolated_field = std::make_shared<shamrock::solvergraph::DeviceBufferEdge<Tfield>>(
410 "interpolated_field", "f_{\\rm interp}");
411
412 auto node = std::make_shared<SPHAzymuthalInteg<Tvec, Tfield, SPHKernel>>();
413 node->set_edges(
414 gpart_mass,
415 tree_reduction_level,
416 part_counts,
417 positions_refs,
418 hpart_refs,
419 field_data,
420 ring_rays_edge,
421 interpolated_field);
422 node->evaluate();
423
424 sham::DeviceBuffer<Tfield> ret{
425 interpolated_field->value.get_size(), shamsys::instance::get_compute_scheduler_ptr()};
426 ret.copy_from(interpolated_field->value);
427
428 return ret;
429 }
430
431 template<class Tvec, class Tfield, template<class> class SPHKernel>
432 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_slice(
433 std::function<field_getter_t> field_getter,
434 Tvec center,
435 Tvec delta_x,
436 Tvec delta_y,
437 u32 nx,
438 u32 ny) -> sham::DeviceBuffer<Tfield> {
439
440 auto positions = pixel_to_positions(center, delta_x, delta_y, nx, ny);
441
442 return compute_slice(field_getter, positions);
443 }
444
445 template<class Tvec, class Tfield, template<class> class SPHKernel>
446 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_column_integ(
447 std::function<field_getter_t> field_getter,
448 Tvec center,
449 Tvec delta_x,
450 Tvec delta_y,
451 u32 nx,
452 u32 ny) -> sham::DeviceBuffer<Tfield> {
453
454 auto rays = pixel_to_orthographic_rays(center, delta_x, delta_y, nx, ny);
455
456 return compute_column_integ(field_getter, rays);
457 }
458
459 template<class Tvec, class Tfield, template<class> class SPHKernel>
460 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_slice(
461 std::string field_name,
462 Tvec center,
463 Tvec delta_x,
464 Tvec delta_y,
465 u32 nx,
466 u32 ny,
467 std::optional<std::function<pybind11::array_t<Tfield>(size_t, pybind11::dict &)>>
468 custom_getter) -> sham::DeviceBuffer<Tfield> {
469 auto positions = pixel_to_positions(center, delta_x, delta_y, nx, ny);
470 return compute_slice(field_name, positions, custom_getter);
471 }
472
473 template<class Tvec, class Tfield, template<class> class SPHKernel>
474 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_column_integ(
475 std::string field_name,
476 Tvec center,
477 Tvec delta_x,
478 Tvec delta_y,
479 u32 nx,
480 u32 ny,
481 std::optional<std::function<pybind11::array_t<Tfield>(size_t, pybind11::dict &)>>
482 custom_getter) -> sham::DeviceBuffer<Tfield> {
483 auto rays = pixel_to_orthographic_rays(center, delta_x, delta_y, nx, ny);
484 return compute_column_integ(field_name, rays, custom_getter);
485 }
486
487} // namespace shammodels::sph::modules
488
489using namespace shammath;
493
497
501
Solver graph edge wrapping a global device buffer.
constexpr const char * part_counts
Particle counts per patch.
SPH azimuthal integration solver graph node.
SPH column integration solver graph node.
SPH slice interpolation solver graph node.
double f64
Alias for double.
std::uint32_t u32
32 bit unsigned integer
std::uint64_t u64
64 bit unsigned integer
void copy_from(const DeviceBuffer< T, new_target > &other, size_t copy_size)
Copies the content of another buffer to this one.
size_t get_size() const
Gets the number of elements in the buffer.
DeviceQueue & get_queue(u32 id=0)
Get a reference to a DeviceQueue.
iterator add_obj(u64 id, T &&obj)
Adds a new object to the collection.
std::string get_time_str() const
Converts the stored nanosecond time to a string representation.
Definition Timer.hpp:79
void start()
Starts the timer.
Definition Timer.hpp:51
void stop()
Stops the timer and stores the elapsed time in nanoseconds.
Definition Timer.hpp:65
u32 get_field_idx(const std::string &field_name) const
Get the field id if matching name & type.
This header file contains utility functions related to exception handling in the code.
void kernel_call(sham::DeviceQueue &q, RefIn in, RefOut in_out, u32 n, Functor &&func, SourceLocation &&callsite=SourceLocation{})
Submit a kernel to a SYCL queue.
ExcptTypes make_except_with_loc(std::string message, SourceLocation loc=SourceLocation{})
Create an exception with a message and a location.
i32 world_rank()
Gives the rank of the current process in the MPI communicator.
Definition worldInfo.cpp:40
namespace for math utility
Definition AABB.hpp:26
namespace for the sph model modules
u64 id_patch
unique key that identify the patch
Definition Patch.hpp:86