35 sham::DeviceBuffer<Tvec> pixel_to_positions(
36 Tvec center, Tvec delta_x, Tvec delta_y,
u32 nx,
u32 ny) {
38 sham::DeviceBuffer<Tvec> ret{nx * ny, shamsys::instance::get_compute_scheduler_ptr()};
40 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().
get_queue();
43 q, sham::MultiRef{}, sham::MultiRef{ret}, nx * ny, [=](
u32 gid, Tvec *position) {
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;
55 sham::DeviceBuffer<shammath::Ray<Tvec>> pixel_to_orthographic_rays(
56 Tvec center, Tvec delta_x, Tvec delta_y,
u32 nx,
u32 ny) {
58 using Tscal = shambase::VecComponent<Tvec>;
60 sham::DeviceBuffer<shammath::Ray<Tvec>> ret{
61 nx * ny, shamsys::instance::get_compute_scheduler_ptr()};
63 sham::DeviceQueue &q = shamsys::instance::get_compute_scheduler().
get_queue();
65 Tvec e_z = sycl::cross(delta_x, delta_y);
66 Tscal len = sycl::length(e_z);
69 "The cross product of delta_x and delta_y is zero\n"
91 [=](
u32 gid, shammath::Ray<Tvec> *ray) {
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;
98 ray[gid] = shammath::Ray<Tvec>(pos_render, e_z);
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> {
113 "sph::CartesianRender",
115 "compute_slice field_name: {}, positions count: {}",
126 [&](
auto field_getter) -> sham::DeviceBuffer<Tfield> {
127 return compute_slice(field_getter, positions);
134 "sph::CartesianRender",
135 shambase::format(
"compute_slice took {}", t.
get_time_str()));
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> {
150 "sph::CartesianRender",
152 "compute_column_integ field_name: {}, rays count: {}",
163 [&](
auto field_getter) -> sham::DeviceBuffer<Tfield> {
164 return compute_column_integ(field_getter, rays);
171 "sph::CartesianRender",
172 shambase::format(
"compute_column_integ took {}", t.
get_time_str()));
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> {
187 "sph::CartesianRender",
189 "compute_azymuthal_integ field_name: {}, ring_rays count: {}",
191 ring_rays.get_size()));
200 [&](
auto field_getter) -> sham::DeviceBuffer<Tfield> {
201 return compute_azymuthal_integ(field_getter, ring_rays);
208 "sph::CartesianRender",
209 shambase::format(
"compute_azymuthal_integ took {}", t.
get_time_str()));
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> {
220 auto part_counts = shamrock::solvergraph::Indexes<u32>::make_shared(
"part_counts",
"N");
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");
225 = std::make_shared<shamrock::solvergraph::Field<Tfield>>(1,
"field_data",
"f");
227 shamrock::solvergraph::DDPatchDataFieldRef<Tvec> pos_dd;
228 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> h_dd;
230 scheduler().for_each_patchdata_nonempty(
231 [&](
const shamrock::patch::Patch cur_p, shamrock::patch::PatchDataLayer &pdat) {
233 u32 cnt = pdat.get_obj_cnt();
236 pos_dd.add_obj(
id, std::ref(pdat.get_field<Tvec>(0)));
238 id, std::ref(pdat.get_field<Tscal>(pdat.pdl().
get_field_idx<Tscal>(
"hpart"))));
241 positions_refs->set_refs(pos_dd);
242 hpart_refs->set_refs(h_dd);
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);
252 auto gpart_mass = shamrock::solvergraph::IDataEdge<Tscal>::make_shared(
"gpart_mass",
"m");
253 gpart_mass->data = solver_config.gpart_mass;
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;
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);
264 auto interpolated_field = std::make_shared<shamrock::solvergraph::DeviceBufferEdge<Tfield>>(
265 "interpolated_field",
"f_{\\rm interp}");
267 auto node = std::make_shared<SPHInterpolation<Tvec, Tfield, SPHKernel>>();
270 tree_reduction_level,
279 sham::DeviceBuffer<Tfield> ret{
280 interpolated_field->value.get_size(), shamsys::instance::get_compute_scheduler_ptr()};
281 ret.
copy_from(interpolated_field->value);
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> {
291 auto part_counts = shamrock::solvergraph::Indexes<u32>::make_shared(
"part_counts",
"N");
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");
296 = std::make_shared<shamrock::solvergraph::Field<Tfield>>(1,
"field_data",
"f");
298 shamrock::solvergraph::DDPatchDataFieldRef<Tvec> pos_dd;
299 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> h_dd;
301 scheduler().for_each_patchdata_nonempty(
302 [&](
const shamrock::patch::Patch cur_p, shamrock::patch::PatchDataLayer &pdat) {
304 u32 cnt = pdat.get_obj_cnt();
307 pos_dd.add_obj(
id, std::ref(pdat.get_field<Tvec>(0)));
309 id, std::ref(pdat.get_field<Tscal>(pdat.pdl().
get_field_idx<Tscal>(
"hpart"))));
312 positions_refs->set_refs(pos_dd);
313 hpart_refs->set_refs(h_dd);
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);
323 auto gpart_mass = shamrock::solvergraph::IDataEdge<Tscal>::make_shared(
"gpart_mass",
"m");
324 gpart_mass->data = solver_config.gpart_mass;
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;
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);
336 auto interpolated_field = std::make_shared<shamrock::solvergraph::DeviceBufferEdge<Tfield>>(
337 "interpolated_field",
"f_{\\rm interp}");
339 auto node = std::make_shared<SPHColumnInteg<Tvec, Tfield, SPHKernel>>();
342 tree_reduction_level,
351 sham::DeviceBuffer<Tfield> ret{
352 interpolated_field->value.get_size(), shamsys::instance::get_compute_scheduler_ptr()};
353 ret.
copy_from(interpolated_field->value);
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> {
364 auto part_counts = shamrock::solvergraph::Indexes<u32>::make_shared(
"part_counts",
"N");
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");
369 = std::make_shared<shamrock::solvergraph::Field<Tfield>>(1,
"field_data",
"f");
371 shamrock::solvergraph::DDPatchDataFieldRef<Tvec> pos_dd;
372 shamrock::solvergraph::DDPatchDataFieldRef<Tscal> h_dd;
374 scheduler().for_each_patchdata_nonempty(
375 [&](
const shamrock::patch::Patch cur_p, shamrock::patch::PatchDataLayer &pdat) {
377 u32 cnt = pdat.get_obj_cnt();
380 pos_dd.add_obj(
id, std::ref(pdat.get_field<Tvec>(0)));
382 id, std::ref(pdat.get_field<Tscal>(pdat.pdl().
get_field_idx<Tscal>(
"hpart"))));
385 positions_refs->set_refs(pos_dd);
386 hpart_refs->set_refs(h_dd);
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);
396 auto gpart_mass = shamrock::solvergraph::IDataEdge<Tscal>::make_shared(
"gpart_mass",
"m");
397 gpart_mass->data = solver_config.gpart_mass;
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;
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);
409 auto interpolated_field = std::make_shared<shamrock::solvergraph::DeviceBufferEdge<Tfield>>(
410 "interpolated_field",
"f_{\\rm interp}");
412 auto node = std::make_shared<SPHAzymuthalInteg<Tvec, Tfield, SPHKernel>>();
415 tree_reduction_level,
424 sham::DeviceBuffer<Tfield> ret{
425 interpolated_field->value.get_size(), shamsys::instance::get_compute_scheduler_ptr()};
426 ret.
copy_from(interpolated_field->value);
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,
438 u32 ny) -> sham::DeviceBuffer<Tfield> {
440 auto positions = pixel_to_positions(center, delta_x, delta_y, nx, ny);
442 return compute_slice(field_getter, positions);
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,
452 u32 ny) -> sham::DeviceBuffer<Tfield> {
454 auto rays = pixel_to_orthographic_rays(center, delta_x, delta_y, nx, ny);
456 return compute_column_integ(field_getter, rays);
459 template<
class Tvec,
class Tfield,
template<
class>
class SPHKernel>
460 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_slice(
461 std::string field_name,
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);
473 template<
class Tvec,
class Tfield,
template<
class>
class SPHKernel>
474 auto CartesianRender<Tvec, Tfield, SPHKernel>::compute_column_integ(
475 std::string field_name,
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);
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.
void start()
Starts the timer.
void stop()
Stops the timer and stores the elapsed time in nanoseconds.
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.
namespace for math utility
namespace for the sph model modules
u64 id_patch
unique key that identify the patch