From ddf2cf29a41a07d6a29c1a923fa1af518b30fb52 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 14:44:10 +0100 Subject: [PATCH 01/23] Add cuBQL submodule and CMake setup --- .gitmodules | 3 +++ CMakeLists.txt | 57 +++++++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 57 insertions(+), 3 deletions(-) diff --git a/.gitmodules b/.gitmodules index 77e2a79a..dd0e33bd 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,3 +20,6 @@ [submodule "vendor/GPRT"] path = vendor/GPRT url = https://github.com/gprt-org/GPRT.git +[submodule "vendor/cuBQL"] + path = vendor/cuBQL + url = https://github.com/NVIDIA/cuBQL diff --git a/CMakeLists.txt b/CMakeLists.txt index df2d6cb3..7196b384 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,6 +7,7 @@ option(XDG_ENABLE_LIBMESH "Enable support for the libMesh mesh library" OFF option(XDG_LINK_MPI "Link with MPI (for dependency compatibility)" OFF) option(XDG_ENABLE_EMBREE "Enable support for the Embree ray tracing library" ON) option(XDG_ENABLE_GPRT "Enable support for the GPRT ray tracing library" OFF) +option(XDG_ENABLE_CUBQL "Enable support for the cuBQL ray tracing library" OFF) option(XDG_BUILD_TESTS "Enable C++ unit testing" ON) option(XDG_BUILD_TOOLS "Enable tools and miniapps" ON) @@ -20,6 +21,10 @@ if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose build type" FORCE) endif() +if(DEFINED XDG_CMAKE_PRESET) + message(STATUS "XDG CMake preset: ${XDG_CMAKE_PRESET}") +endif() + # Compiler options (things in this section may not be platform-portable) set(CMAKE_POSITION_INDEPENDENT_CODE ON) @@ -111,6 +116,12 @@ if(XDG_ENABLE_GPRT) ) endif() +if(XDG_ENABLE_CUBQL) + list(APPEND VENDOR_PATHS + vendor/cuBQL + ) +endif() + if(GIT_FOUND AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/.git") option(XDG_GIT_SUBMODULE "Check submodules during build" ON) if(XDG_GIT_SUBMODULE) @@ -175,11 +186,12 @@ if (NOT XDG_ENABLE_MOAB AND NOT XDG_ENABLE_LIBMESH) endif() # Ensure at least one ray tracing backend is enabled -if (NOT XDG_ENABLE_EMBREE AND NOT XDG_ENABLE_GPRT) +if (NOT XDG_ENABLE_EMBREE AND NOT XDG_ENABLE_GPRT AND NOT XDG_ENABLE_CUBQL) message(FATAL_ERROR "No ray tracing backend enabled. Enable at least one of:\n" " -DXDG_ENABLE_EMBREE=ON\n" - " -DXDG_ENABLE_GPRT=ON") + " -DXDG_ENABLE_GPRT=ON\n" + " -DXDG_ENABLE_CUBQL=ON") endif() # GPRT @@ -188,6 +200,11 @@ if (XDG_ENABLE_GPRT) add_subdirectory(vendor/GPRT) endif() +if (XDG_ENABLE_CUBQL) + set(CUBQL_OMP ON CACHE BOOL "Build cuBQL with OpenMP target offload" FORCE) + add_subdirectory(vendor/cuBQL) +endif() + list(APPEND xdg_sources src/geometry/measure.cpp src/geometry/plucker.cpp @@ -221,6 +238,23 @@ dbl_deviceCode endif() +if (XDG_ENABLE_CUBQL) +list(APPEND xdg_sources +src/cuBQL/triangles.cpp +src/cuBQL/intersection.cpp +src/cuBQL/ray_tracer.cpp +) + +# We need a precompile definition to switch to using the cuBQL math types in the shared +# plucker intersection code. The compile definition is used in dp__math.h +set_source_files_properties( + src/cuBQL/ray_tracer.cpp + src/cuBQL/intersection.cpp + PROPERTIES COMPILE_DEFINITIONS XDG_DP_MATH_CUBQL +) + +endif() + if (XDG_ENABLE_LIBMESH) list(APPEND xdg_sources src/libmesh/mesh_manager.cpp @@ -281,7 +315,12 @@ if (${CMAKE_BUILD_TYPE} MATCHES "Debug") endif() # attempt to find OpenMP and include it if found -find_package(OpenMP) +if (XDG_ENABLE_CUBQL) + find_package(OpenMP REQUIRED) +else() + find_package(OpenMP) +endif() + if (OpenMP_CXX_FOUND) target_link_libraries(xdg PRIVATE OpenMP::OpenMP_CXX) target_compile_definitions(xdg PRIVATE XDG_HAVE_OPENMP) @@ -319,6 +358,18 @@ if (XDG_ENABLE_GPRT) target_link_options(xdg PRIVATE -Wl,--unresolved-symbols=ignore-in-shared-libs) endif() +if (XDG_ENABLE_CUBQL) + target_compile_definitions(xdg PUBLIC XDG_ENABLE_CUBQL) + target_link_libraries(xdg PRIVATE $) + # TODO: Stop relying on LD_LIBRARY_PATH for LLVM OpenMP offload runtimes. + # As a temporary measure whilst figuring out our way around cuBQL this is + # okay but in the long run we should aim for a more robust solutions here. + # Clang injects libomptarget when offload flags are supplied by presets, but + # CMake does not currently add that compiler runtime directory to RPATH. + # Add targeted BUILD_RPATH handling for libomptarget/libomp, and decide + # whether install RPATH should remain environment-module based or be opt-in. +endif() + target_link_libraries(xdg PRIVATE fmt::fmt) # ========================== From f17419ccf9fcfc9e714c7f5d09239db9dd890e35 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 14:52:51 +0100 Subject: [PATCH 02/23] Added CuBQLRayTracer backend support - New CuBQLRayTracer class with working create_surface_tree and basic ray queries - Updated plucker_ray_tri intersect to play nice with openmp target offload regions - Implemented new structs for cuBQL BVH objects with XDG data - Implemented mixed precision BVH traversal algorithm for cuBQL - Updated appropriate constants with new cuBQL specifics - Ensured cuBQL properly wired through to xdg public API --- include/xdg/constants.h | 8 +- include/xdg/cuBQL/cuBQL_backend.h | 30 ++ include/xdg/cuBQL/intersection.h | 89 ++++++ include/xdg/cuBQL/ray_tracer.h | 94 +++++++ include/xdg/cuBQL/triangles.h | 126 +++++++++ include/xdg/geometry/dp_math.h | 22 +- include/xdg/geometry/plucker.h | 16 +- include/xdg/ray_tracers.h | 6 +- src/cuBQL/intersection.cpp | 210 ++++++++++++++ src/cuBQL/ray_tracer.cpp | 437 ++++++++++++++++++++++++++++++ src/cuBQL/triangles.cpp | 52 ++++ src/xdg.cpp | 14 + 12 files changed, 1090 insertions(+), 14 deletions(-) create mode 100644 include/xdg/cuBQL/cuBQL_backend.h create mode 100644 include/xdg/cuBQL/intersection.h create mode 100644 include/xdg/cuBQL/ray_tracer.h create mode 100644 include/xdg/cuBQL/triangles.h create mode 100644 src/cuBQL/intersection.cpp create mode 100644 src/cuBQL/ray_tracer.cpp create mode 100644 src/cuBQL/triangles.cpp diff --git a/include/xdg/constants.h b/include/xdg/constants.h index 38df9db1..e143d537 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -54,7 +54,8 @@ enum class MeshLibrary { // Ray Tracing library identifier enum class RTLibrary { EMBREE, - GPRT + GPRT, + CUBQL }; static const std::map MESH_LIB_TO_STR = @@ -67,7 +68,8 @@ static const std::map MESH_LIB_TO_STR = static const std::map RT_LIB_TO_STR = { {RTLibrary::EMBREE, "EMBREE"}, - {RTLibrary::GPRT, "GPRT"} + {RTLibrary::GPRT, "GPRT"}, + {RTLibrary::CUBQL, "CUBQL"} }; // Mesh identifer type @@ -148,4 +150,4 @@ struct formatter : fmt::formatter { } -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/cuBQL/cuBQL_backend.h b/include/xdg/cuBQL/cuBQL_backend.h new file mode 100644 index 00000000..b348802f --- /dev/null +++ b/include/xdg/cuBQL/cuBQL_backend.h @@ -0,0 +1,30 @@ +#ifndef _XDG_CUBQL_BACKEND_H +#define _XDG_CUBQL_BACKEND_H + +#include +#include +#include +#include + +#include + +// Guards to prevent CUDA headers from being included in host code, which causes +// failed compilation with LLVM-clang. +#if defined(__CUDA_ARCH__) && !defined(__CUDACC__) +#undef __CUDA_ARCH__ +#endif + +#include "xdg/error.h" + +namespace xdg::cubql { + +struct Context { + int gpuID {0}; + int hostID {omp_get_initial_device()}; +}; + + + +} // namespace xdg::cubql + +#endif // include guard diff --git a/include/xdg/cuBQL/intersection.h b/include/xdg/cuBQL/intersection.h new file mode 100644 index 00000000..50d95141 --- /dev/null +++ b/include/xdg/cuBQL/intersection.h @@ -0,0 +1,89 @@ +#ifndef _XDG_CUBQL_INTERSECTION_H +#define _XDG_CUBQL_INTERSECTION_H + +// Guards to prevent CUDA headers from being included in host code, which causes +// failed compilation with LLVM-clang. +#if defined(__CUDA_ARCH__) && !defined(__CUDACC__) +#undef __CUDA_ARCH__ +#endif + +#include +#include + +#include "xdg/constants.h" +#include "xdg/cuBQL/triangles.h" +#include "cuBQL/math/vec.h" + +namespace xdg { + +struct CuBQLRay { + cuBQL::vec3d origin; + cuBQL::vec3d direction; + double tMin {0.0}; + double tMax {INFTY}; + MeshID volume {ID_NONE}; // volume we are tracing ray against +}; + +/* POD SurfaceRay struct for external population*/ +// struct CuBQLSurfaceRay { +// double origin[3]; +// double direction[3]; +// uint32_t volume_slot; +// uint32_t enabled; +// const MeshID* exclude_primitives; +// int32_t exclude_count; +// }; + +// TODO - Consider whether this is useful/necessary as its own struct +// struct CuBQLExcludeList { +// const MeshID* primitives {nullptr}; +// int count {0}; +// }; + +struct CuBQLSurfaceHit { + double distance {INFTY}; + MeshID surface {ID_NONE}; + MeshID primitive {ID_NONE}; + PointInVolume piv {OUTSIDE}; + + bool hit_found() const { return primitive != ID_NONE; } +}; + +inline bool orientation_cull(double normal_dot_direction, + HitOrientation orientation) +{ + if (orientation == HitOrientation::ANY) return false; + + if (orientation == HitOrientation::EXITING && normal_dot_direction < 0.0) { + return true; + } else if (orientation == HitOrientation::ENTERING && normal_dot_direction >= 0.0) { + return true; + } + + return false; +} + +/* +Wrapper for launching a single ray intersection query against the surface tree, with Host<->Device staging of ray and hit data +Performs host side staging and transfer hit data back to host after device side traversal +*/ +void +intersect_surface_tree_scalar(const cubql::Context& context, + const CuBQLVolumeTLAS& volume_tlas, + const CuBQLRay& ray, + CuBQLSurfaceHit& hit, + HitOrientation hit_orientation, + const std::vector* exclude_primitives); + + +void +intersect_surface_tree_batch(const cubql::Context& context, + const CuBQLVolumeTLAS::DD* d_volume_to_tlas, + const CuBQLRay* d_rays, + CuBQLSurfaceHit* d_hits, + std::size_t num_rays, + HitOrientation hit_orientation); + +} // namespace xdg + +#endif // include guard diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h new file mode 100644 index 00000000..c422ecb3 --- /dev/null +++ b/include/xdg/cuBQL/ray_tracer.h @@ -0,0 +1,94 @@ +#ifndef _XDG_CUBQL_RAY_TRACING_INTERFACE_H +#define _XDG_CUBQL_RAY_TRACING_INTERFACE_H + +#include +#include +#include +#include +#include + +#include "xdg/constants.h" +#include "xdg/geometry_data.h" +#include "xdg/cuBQL/triangles.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/ray.h" +#include "xdg/ray_tracing_interface.h" + +namespace xdg { + +struct CuBQLRay; +struct CuBQLSurfaceHit; + +class CuBQLRayTracer : public RayTracer { +public: + CuBQLRayTracer(); + ~CuBQLRayTracer() override; + + RTLibrary library() const override { return RTLibrary::CUBQL; } + + void init() override; + + std::pair + register_volume(const std::shared_ptr& mesh_manager, + MeshID volume) override; + + TreeID create_surface_tree(const std::shared_ptr& mesh_manager, + MeshID volume) override; + + TreeID create_element_tree(const std::shared_ptr& mesh_manager, + MeshID volume) override; + + void create_global_surface_tree() override; + + void create_global_element_tree() override; + + MeshID find_element(const Position& point) const override; + + MeshID find_element(TreeID tree, const Position& point) const override; + + bool point_in_volume(TreeID tree, + const Position& point, + const Direction* direction = nullptr, + const std::vector* exclude_primitives = nullptr) const override; + + std::pair ray_fire(TreeID tree, + const Position& origin, + const Direction& direction, + const double dist_limit = INFTY, + HitOrientation orientation = HitOrientation::EXITING, + std::vector* const exclude_primitives = nullptr) override; + + void ray_fire_batch(const CuBQLRay* d_rays, + CuBQLSurfaceHit* d_hits, + std::size_t num_rays, + HitOrientation orientation = HitOrientation::EXITING); + + std::pair closest(TreeID tree, + const Position& origin) override; + + bool occluded(TreeID tree, + const Position& origin, + const Direction& direction, + double& dist) const override; + +private: + CuBQLSurfaceBLAS + register_surface(const std::shared_ptr& mesh_manager, + MeshID surface_id, + double bounding_box_bump); + + void upload_volume_to_tlas_table_(); + + cubql::Context context_; + + std::unordered_map tree_to_volume_tlas_; + std::unordered_map surface_to_blas_map_; + + std::vector volume_to_tlas_; + CuBQLVolumeTLAS::DD* d_volume_to_tlas_ {nullptr}; + bool initialized_ {false}; +}; + +} // namespace xdg + +#endif // include guard diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h new file mode 100644 index 00000000..20e53c91 --- /dev/null +++ b/include/xdg/cuBQL/triangles.h @@ -0,0 +1,126 @@ +#ifndef _XDG_CUBQL_TRIANGLES_H +#define _XDG_CUBQL_TRIANGLES_H + +#include +#include + +// Guards to prevent CUDA headers from being included in host code, which causes +// failed compilation with LLVM-clang. +#if defined(__CUDA_ARCH__) && !defined(__CUDACC__) +#undef __CUDA_ARCH__ +#endif + +#include "cuBQL/bvh.h" +#include "cuBQL/math/vec.h" +#include "xdg/constants.h" +#include "xdg/cuBQL/cuBQL_backend.h" + +namespace xdg { + +/* + Owns the triangle buffers for one topological surface. The nested DD type is + the compact device-data view copied into OpenMP target regions instead of the + full host-side owner. +*/ +struct CuBQLSurfaceMesh { + struct DD { + // Topological metadata + MeshID surface_id {ID_NONE}; + + // Geometric data + const cuBQL::vec3d* vertices {nullptr}; + const cuBQL::vec3i* indices {nullptr}; + const MeshID* primitive_refs {nullptr}; + }; + + // Topological metadata + MeshID surface_id {ID_NONE}; + + // Device buffers for triangle data + cuBQL::vec3d* d_vertices {nullptr}; + cuBQL::vec3i* d_indices {nullptr}; + MeshID* d_primitive_refs {nullptr}; + + uint32_t num_vertices {0}; + uint32_t num_triangles {0}; + int gpu_id {0}; + + // Accessor for Device Data struct, which is passed to cuBQL BVH traversal/intersection functions + DD get_device_data() const + { + return { + surface_id, + d_vertices, + d_indices, + d_primitive_refs + }; + } + + void release(); +}; + +/* + Owns a cuBQL BVH used as a Bottom Level Acceleration Structure over surface triangles. + The nested DD type is the compact device-data view used during traversal. +*/ +struct CuBQLSurfaceBLAS { + struct DD { + CuBQLSurfaceMesh::DD mesh; // Mesh data device handle + cuBQL::bvh3f bvh; // BLAS device handle + }; + + cuBQL::bvh3f bvh; // BLAS host handle + CuBQLSurfaceMesh mesh; // Surface mesh host owner + + uint32_t num_prims {0}; + int gpu_id {0}; + + DD get_device_data() const + { + return {mesh.get_device_data(), bvh}; + } + + void release(); +}; + +/* + Owns a cuBQL BVH used as a Top-Level Acceleration Structure for one topological volume. + The TLAS groups the surface BLASes that bound that volume and stores + per-volume relationship metadata for each surface instance. +*/ +struct CuBQLVolumeTLAS { + /* + TLAS-local instance payload. The same surface BLAS can participate in + different volume TLASes with different sense, so reverse_sense belongs on + the volume-surface relationship rather than on the reusable surface mesh + or BLAS geometry. + */ + struct SurfaceInstanceDD { + CuBQLSurfaceBLAS::DD surface_blas; + bool reverse_sense {false}; // value set in create_surface_tree based on parent vols + }; + + struct DD { + MeshID volume_id {ID_NONE}; + const SurfaceInstanceDD* surface_instances {nullptr}; + cuBQL::bvh3f bvh; // TLAS device handle + }; + + MeshID volume_id {ID_NONE}; + cuBQL::bvh3f bvh; // TLAS host handle + SurfaceInstanceDD* d_surface_instances {nullptr}; + + uint32_t num_surface_instances {0}; + int gpu_id {0}; + + DD get_device_data() const + { + return {volume_id, d_surface_instances, bvh}; + } + + void release(); +}; + +} // namespace xdg + +#endif // include guard diff --git a/include/xdg/geometry/dp_math.h b/include/xdg/geometry/dp_math.h index 5b162d3f..5dfedffd 100644 --- a/include/xdg/geometry/dp_math.h +++ b/include/xdg/geometry/dp_math.h @@ -22,7 +22,27 @@ namespace dp { static const double INFTY = 1.7976931348623157e+308; // std::numeric_limits::max() is not available in slang } -#else +// TODO - Is this the right way to handle this for cubql openmp target offload compilation? +// In theory, we can compile the C++ pathway that Embree uses but will the vec3da types be +// omptarget friendly? +// For now I have defined this separate compilation pathway which is enabled with a new +// precompile definition only set when when compiling the CuBQLRayTracer that maps to cuBQL's math types +#elif defined(XDG_DP_MATH_CUBQL) +#include "cuBQL/math/common.h" + +// C++ compilation with cuBQL math types, map dp::vec3 -> cuBQL::vec3d +namespace dp { + typedef cuBQL::vec3d vec3; + + inline double dot(vec3 a, vec3 b) { return cuBQL::dot(a, b); } + inline vec3 cross(vec3 a, vec3 b) { return cuBQL::cross(a, b); } + inline double abs(double a) { return cuBQL::abst(a); } + + static constexpr double DBL_ZERO_TOL = 20.0 * 2.2204460492503131e-16; // same as 20 * std::numeric_limits::epsilon() + constexpr double INFTY {std::numeric_limits::max()}; +} + +#else #include "xdg/vec3da.h" // C++ compilation map dp::vec3 -> xdg::Vec3da diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index acf0ebfa..8b7581c9 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -25,8 +25,6 @@ struct PluckerIntersectionResult { double t = 0.0; // Distance along the ray to the intersection point }; -static constexpr PluckerIntersectionResult EXIT_EARLY = {false, 0.0}; - /* Function to return the vertex with the lowest coordinates. To force the same ray-edge computation, the Plücker test needs to use consistent edge representation. This would be more simple with MOAB handles instead of @@ -81,7 +79,7 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], // If orientation is set, confirm that sign of plucker_coordinate indicate // correct orientation of intersection if (useOrientation && orientation * plucker_coord0 > 0) { - return EXIT_EARLY; + return {false, 0.0}; } // Determine the value of the second Plucker coordinate from edge 1 @@ -92,13 +90,13 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], // correct orientation of intersection if (useOrientation) { if (orientation * plucker_coord1 > 0) { - return EXIT_EARLY; + return {false, 0.0}; } // If the orientation is not specified, all plucker_coords must be the same // sign or zero. } else if ((0.0 < plucker_coord0 && 0.0 > plucker_coord1) || (0.0 > plucker_coord0 && 0.0 < plucker_coord1)) { - return EXIT_EARLY; + return {false, 0.0}; } // Determine the value of the third Plucker coordinate from edge 2 @@ -109,7 +107,7 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], // correct orientation of intersection if (useOrientation) { if (orientation * plucker_coord2 > 0) { - return EXIT_EARLY; + return {false, 0.0}; } // If the orientation is not specified, all plucker_coords must be the same // sign or zero. @@ -117,12 +115,12 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], (0.0 > plucker_coord1 && 0.0 < plucker_coord2) || (0.0 < plucker_coord0 && 0.0 > plucker_coord2) || (0.0 > plucker_coord0 && 0.0 < plucker_coord2)) { - return EXIT_EARLY; + return {false, 0.0}; } // check for coplanar case to avoid dividing by zero if (0.0 == plucker_coord0 && 0.0 == plucker_coord1 && 0.0 == plucker_coord2) { - return EXIT_EARLY; + return {false, 0.0}; } // get the distance to intersection @@ -155,7 +153,7 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], } // is the intersection within distance limits? - if (dist_out < tMin || dist_out > tMax) return EXIT_EARLY; + if (dist_out < tMin || dist_out > tMax) return {false, 0.0}; return {true, dist_out}; } diff --git a/include/xdg/ray_tracers.h b/include/xdg/ray_tracers.h index b4816571..2031baf0 100644 --- a/include/xdg/ray_tracers.h +++ b/include/xdg/ray_tracers.h @@ -5,4 +5,8 @@ #ifdef XDG_ENABLE_GPRT #include "xdg/gprt/ray_tracer.h" -#endif \ No newline at end of file +#endif + +#ifdef XDG_ENABLE_CUBQL +#include "xdg/cuBQL/ray_tracer.h" +#endif diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp new file mode 100644 index 00000000..cc5b1bdb --- /dev/null +++ b/src/cuBQL/intersection.cpp @@ -0,0 +1,210 @@ +#include + +#include "xdg/cuBQL/intersection.h" +#include "xdg/geometry/plucker.h" +#include "xdg/error.h" + +#include "cuBQL/math/Ray.h" +#include "cuBQL/traversal/rayQueries.h" + +namespace xdg { + +// Core traversal and intersection routine for a single ray against a given volume tlas +#pragma omp declare target +static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, + CuBQLRay intersection_ray, + CuBQLSurfaceHit* hit, + int orientation, + const MeshID* exclude_primitives, + int exclude_count) +{ + // cuBQL traverses the FP32 BVH with an FP32 ray; the original CuBQLRay + // remains the FP64 source of truth for the final triangle intersection. + cuBQL::ray3f traversal_ray; + traversal_ray.origin = cuBQL::vec3f(intersection_ray.origin); + traversal_ray.direction = cuBQL::vec3f(intersection_ray.direction); + // TODO: Is this truncation safe enough for tmin and tmax? Pretty sure embree/gprt does a similar truncation for ray bounds + traversal_ray.tMin = static_cast(intersection_ray.tMin); + traversal_ray.tMax = static_cast(hit->distance); + + CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; + + auto enter_blas = [=, &surface_instance, &traversal_ray] + (cuBQL::ray3f& out_ray, cuBQL::bvh3f& out_bvh, int instance_id) + { + surface_instance = volume_tlas.surface_instances[instance_id]; + out_ray = traversal_ray; + out_bvh = surface_instance.surface_blas.bvh; + }; + + auto intersect_prim = [=, &traversal_ray, &surface_instance] + (uint32_t prim_id) -> float + { + const CuBQLSurfaceMesh::DD mesh = surface_instance.surface_blas.mesh; + const MeshID primitive_ref = mesh.primitive_refs[prim_id]; + + for (int i = 0; i < exclude_count; ++i) { + if (exclude_primitives[i] == primitive_ref) { + return traversal_ray.tMax; + } + } + + const cuBQL::vec3i index = mesh.indices[prim_id]; + + cuBQL::vec3d vertices[3] = { + mesh.vertices[index.x], + mesh.vertices[index.y], + mesh.vertices[index.z] + }; + + cuBQL::vec3d normal = cuBQL::cross(vertices[1] - vertices[0], + vertices[2] - vertices[0]); + + if (surface_instance.reverse_sense) { + normal = -normal; + } + + const double normal_dot_direction = dot(normal, intersection_ray.direction); + + if (orientation_cull(normal_dot_direction, + static_cast(orientation))) { + return traversal_ray.tMax; + } + + auto intersection = plucker_ray_tri_intersect(vertices, + intersection_ray.origin, + intersection_ray.direction, + hit->distance, + intersection_ray.tMin, + false, + 0); + + if (intersection.hit) { + hit->distance = intersection.t; + hit->surface = mesh.surface_id; + hit->primitive = primitive_ref; + hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; + traversal_ray.tMax = static_cast(intersection.t); + } + + // Return value is only the FP32 traversal shrink distance. The accepted hit + // distance stored above remains the FP64 Plucker result. + return traversal_ray.tMax; + }; + + auto leave_blas = []() -> void {}; + + cuBQL::shrinkingRayQuery::twoLevel::forEachPrim(enter_blas, + leave_blas, + intersect_prim, + volume_tlas.bvh, + traversal_ray); +} +#pragma omp end declare target + +void +intersect_surface_tree_scalar(const cubql::Context& context, + const CuBQLVolumeTLAS& volume_tlas, + const CuBQLRay& ray, + CuBQLSurfaceHit& surface_hit, + HitOrientation hit_orientation, + const std::vector* exclude_primitives) +{ + const int gpu_id = context.gpuID; + + MeshID* d_exclude_primitives = nullptr; + int exclude_count = 0; + if (exclude_primitives && !exclude_primitives->empty()) { + exclude_count = static_cast(exclude_primitives->size()); + d_exclude_primitives = static_cast + (omp_target_alloc(exclude_count * sizeof(MeshID), gpu_id)); + omp_target_memcpy(d_exclude_primitives, + exclude_primitives->data(), + exclude_count * sizeof(MeshID), + 0, + 0, + gpu_id, + context.hostID); + } + + auto* d_surface_hit = static_cast + (omp_target_alloc(sizeof(CuBQLSurfaceHit), gpu_id)); + + surface_hit.distance = ray.tMax; + omp_target_memcpy(d_surface_hit, + &surface_hit, + sizeof(CuBQLSurfaceHit), + 0, + 0, + gpu_id, + context.hostID); + + const auto volume_tlas_dd = volume_tlas.get_device_data(); + const int orientation = static_cast(hit_orientation); + + #pragma omp target device(gpu_id) \ + is_device_ptr(d_exclude_primitives, d_surface_hit) + { + intersect_surface_tree(volume_tlas_dd, + ray, + d_surface_hit, + orientation, + d_exclude_primitives, + exclude_count); + } + + omp_target_memcpy(&surface_hit, + d_surface_hit, + sizeof(CuBQLSurfaceHit), + 0, + 0, + context.hostID, + gpu_id); + + omp_target_free(d_surface_hit, gpu_id); + + if (d_exclude_primitives) { + omp_target_free(d_exclude_primitives, gpu_id); + } + + return; +} + +void +intersect_surface_tree_batch(const cubql::Context& context, + const CuBQLVolumeTLAS::DD* d_volume_to_tlas, + const CuBQLRay* d_rays, + CuBQLSurfaceHit* d_hits, + std::size_t num_rays, + HitOrientation hit_orientation) +{ + + if (num_rays == 0) return; + + if (!d_volume_to_tlas || !d_rays || !d_hits) { + fatal_error("Invalid cuBQL batch intersection buffers"); + } + + const int gpu_id = context.gpuID; + + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(d_volume_to_tlas, d_rays, d_hits) + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + const CuBQLRay ray = d_rays[ray_id]; + const CuBQLVolumeTLAS::DD volume_tlas = d_volume_to_tlas[ray.volume]; + + CuBQLSurfaceHit hit; + hit.distance = ray.tMax; + + intersect_surface_tree(volume_tlas, + ray, + &hit, + static_cast(hit_orientation), + nullptr, + 0); + + d_hits[ray_id] = hit; + } +} + +} // namespace xdg diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp new file mode 100644 index 00000000..c41bc4f2 --- /dev/null +++ b/src/cuBQL/ray_tracer.cpp @@ -0,0 +1,437 @@ +#include "xdg/cuBQL/ray_tracer.h" +#include "xdg/cuBQL/intersection.h" +#include "xdg/error.h" +#include "xdg/geometry/plucker.h" +#include "xdg/available_device_probe.h" + +#include +#include "cuBQL/builder/omp.h" +#include "cuBQL/math/Ray.h" +#include "cuBQL/queries/triangleData/Triangle.h" +#include "cuBQL/queries/triangleData/math/rayTriangleIntersections.h" +#include "cuBQL/traversal/rayQueries.h" + +#include + +namespace xdg { + +CuBQLRayTracer::CuBQLRayTracer() +{ + if (!system_has_omp_target_device()) { + fatal_error("No OpenMP target capable device found; cannot initialize cuBQL ray tracer."); + } + + context_.gpuID = 0; // TODO - support selecting among multiple OpenMP target devices. + context_.hostID = omp_get_initial_device(); +} + +CuBQLRayTracer::~CuBQLRayTracer() +{ + if (d_volume_to_tlas_) { + omp_target_free(d_volume_to_tlas_, context_.gpuID); + d_volume_to_tlas_ = nullptr; + } + + for (auto& [tree, tlas] : tree_to_volume_tlas_) { + tlas.release(); + } + + for (auto& [surface, blas] : surface_to_blas_map_) { + blas.release(); + } +} + +void CuBQLRayTracer::init() +{ + upload_volume_to_tlas_table_(); + initialized_ = true; +} + +void CuBQLRayTracer::upload_volume_to_tlas_table_() +{ + if (d_volume_to_tlas_) { + omp_target_free(d_volume_to_tlas_, context_.gpuID); + d_volume_to_tlas_ = nullptr; + } + + if (volume_to_tlas_.empty()) { + return; + } + + d_volume_to_tlas_ = static_cast + (omp_target_alloc(volume_to_tlas_.size() * sizeof(CuBQLVolumeTLAS::DD), context_.gpuID)); + omp_target_memcpy(d_volume_to_tlas_, + volume_to_tlas_.data(), + volume_to_tlas_.size() * sizeof(CuBQLVolumeTLAS::DD), + 0, + 0, + context_.gpuID, + context_.hostID); +} + +std::pair +CuBQLRayTracer::register_volume(const std::shared_ptr& mesh_manager, + MeshID volume) +{ + TreeID surface_tree = create_surface_tree(mesh_manager, volume); + TreeID element_tree = create_element_tree(mesh_manager, volume); + return {surface_tree, element_tree}; +} + +CuBQLSurfaceBLAS +CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manager, + MeshID surface_id, + double bounding_box_bump) +{ + auto num_faces = mesh_manager->num_surface_faces(surface_id); + auto vertices = mesh_manager->get_surface_vertices(surface_id); + auto indices = mesh_manager->get_surface_connectivity(surface_id); + + std::vector h_vertices; + h_vertices.reserve(vertices.size()); + for (const auto& vertex : vertices) { + h_vertices.emplace_back(vertex.x, vertex.y, vertex.z); + } + + std::vector h_indices; + h_indices.reserve(indices.size() / 3); + for (size_t i = 0; i < indices.size(); i += 3) { + h_indices.emplace_back(indices[i], indices[i + 1], indices[i + 2]); + } + + std::vector h_primitive_refs = mesh_manager->get_surface_faces(surface_id); + + // TODO- think about how to better handle omp transfer calls. AutoUploadArrays is one option + auto* d_vertices = static_cast + (omp_target_alloc(h_vertices.size() * sizeof(cuBQL::vec3d), context_.gpuID)); + omp_target_memcpy(d_vertices, + h_vertices.data(), + h_vertices.size() * sizeof(cuBQL::vec3d), + 0, + 0, + context_.gpuID, + context_.hostID); + + auto* d_indices = static_cast + (omp_target_alloc(h_indices.size() * sizeof(cuBQL::vec3i), context_.gpuID)); + omp_target_memcpy(d_indices, + h_indices.data(), + h_indices.size() * sizeof(cuBQL::vec3i), + 0, + 0, + context_.gpuID, + context_.hostID); + + auto* d_primitive_refs = static_cast + (omp_target_alloc(h_primitive_refs.size() * sizeof(MeshID), context_.gpuID)); + omp_target_memcpy(d_primitive_refs, + h_primitive_refs.data(), + h_primitive_refs.size() * sizeof(MeshID), + 0, + 0, + context_.gpuID, + context_.hostID); + + auto* d_aabbs = static_cast + (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3f), context_.gpuID)); + const auto num_primitives = static_cast(h_indices.size()); + + // TODO - Abstract this out into its own bounding_box creation function + #pragma omp target device(context_.gpuID) is_device_ptr(d_vertices, d_indices, d_aabbs) \ + firstprivate(bounding_box_bump) + #pragma omp teams distribute parallel for + for (uint32_t primID = 0; primID < num_primitives; ++primID) { + cuBQL::vec3i indices = d_indices[primID]; + + cuBQL::vec3d A = d_vertices[indices.x]; + cuBQL::vec3d B = d_vertices[indices.y]; + cuBQL::vec3d C = d_vertices[indices.z]; + + cuBQL::box3d aabb; + aabb.extend(A); + aabb.extend(B); + aabb.extend(C); + + const cuBQL::vec3d bump(bounding_box_bump); + aabb.lower = aabb.lower - bump; + aabb.upper = aabb.upper + bump; + + d_aabbs[primID] = cuBQL::box3f(aabb); + } + + cuBQL::BuildConfig blasBuildParams; + // TODO - Try setting leaf params to 1 to see what it does + // Check what default is for CUDA + cuBQL::bvh3f bvh; + cuBQL::build_omp_target(bvh, d_aabbs, num_faces, blasBuildParams, context_.gpuID); + + omp_target_free(d_aabbs, context_.gpuID); + + CuBQLSurfaceMesh surface_mesh; + surface_mesh.surface_id = surface_id; + surface_mesh.d_vertices = d_vertices; + surface_mesh.d_indices = d_indices; + surface_mesh.d_primitive_refs = d_primitive_refs; + surface_mesh.num_vertices = h_vertices.size(); + surface_mesh.num_triangles = num_faces; + surface_mesh.gpu_id = context_.gpuID; + + CuBQLSurfaceBLAS surface_blas; + surface_blas.bvh = bvh; + surface_blas.mesh = surface_mesh; + surface_blas.num_prims = num_faces; + surface_blas.gpu_id = context_.gpuID; + + return surface_blas; +} + +TreeID +CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_manager, + MeshID volume_id) +{ + // TODO - Right now each CuBQLRayTracer instance has a single "Context" which holds a single GPU_ID + // so this will need to be reworked in the future to handle multi-gpus + + SurfaceTreeID tree = next_surface_tree_id(); + surface_trees_.push_back(tree); + auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); + std::vector h_tlas_boxes; + std::vector h_surface_instances; + h_tlas_boxes.reserve(volume_surfaces.size()); + h_surface_instances.reserve(volume_surfaces.size()); + + for (const auto &surf : volume_surfaces) { + auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); + const double max_parent_bbox_bump = std::max(bounding_box_bump(mesh_manager, forward_parent), + bounding_box_bump(mesh_manager, reverse_parent)); + + if (!surface_to_blas_map_.count(surf)) { + surface_to_blas_map_[surf] = register_surface(mesh_manager, surf, max_parent_bbox_bump); + } + + CuBQLSurfaceBLAS& surface_blas = surface_to_blas_map_.at(surf); + + // Store BLAS bounding boxes to build TLAS + const auto surface_bounding_box = mesh_manager->surface_bounding_box(surf); + cuBQL::box3d surface_bounds_dp; + surface_bounds_dp.lower = cuBQL::vec3d(surface_bounding_box.min_x, + surface_bounding_box.min_y, + surface_bounding_box.min_z); + surface_bounds_dp.upper = cuBQL::vec3d(surface_bounding_box.max_x, + surface_bounding_box.max_y, + surface_bounding_box.max_z); + + const cuBQL::vec3d bump(max_parent_bbox_bump); + surface_bounds_dp.lower = surface_bounds_dp.lower - bump; + surface_bounds_dp.upper = surface_bounds_dp.upper + bump; + cuBQL::box3f surface_bounds(surface_bounds_dp); + + CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; + surface_instance.surface_blas = surface_blas.get_device_data(); + + // Sense setting for each surface instance in the TLAS + if (volume_id == forward_parent) { + surface_instance.reverse_sense = false; + } else if (volume_id == reverse_parent) { + surface_instance.reverse_sense = true; + } else { + fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); + } + + h_tlas_boxes.push_back(surface_bounds); + h_surface_instances.push_back(surface_instance); + } + + if (h_surface_instances.empty()) { + fatal_error("Volume {} has no surfaces; cannot build cuBQL surface tree", volume_id); + } + + auto* d_tlas_boxes = static_cast + (omp_target_alloc(h_tlas_boxes.size() * sizeof(cuBQL::box3f), context_.gpuID)); + omp_target_memcpy(d_tlas_boxes, + h_tlas_boxes.data(), + h_tlas_boxes.size() * sizeof(cuBQL::box3f), + 0, + 0, + context_.gpuID, + context_.hostID); + + auto* d_surface_instances = static_cast + (omp_target_alloc(h_surface_instances.size() * sizeof(CuBQLVolumeTLAS::SurfaceInstanceDD), context_.gpuID)); + omp_target_memcpy(d_surface_instances, + h_surface_instances.data(), + h_surface_instances.size() * sizeof(CuBQLVolumeTLAS::SurfaceInstanceDD), + 0, + 0, + context_.gpuID, + context_.hostID); + + cuBQL::BuildConfig tlasBuildParams; + tlasBuildParams.makeLeafThreshold = 1; + tlasBuildParams.maxAllowedLeafSize = 1; + + CuBQLVolumeTLAS volume_tlas; + volume_tlas.volume_id = volume_id; // store meshid in the TLAS object for easier mapping between the two + volume_tlas.num_surface_instances = static_cast(h_surface_instances.size()); + volume_tlas.gpu_id = context_.gpuID; + volume_tlas.d_surface_instances = d_surface_instances; + cuBQL::build_omp_target(volume_tlas.bvh, + d_tlas_boxes, + volume_tlas.num_surface_instances, + tlasBuildParams, + context_.gpuID); + + omp_target_free(d_tlas_boxes, context_.gpuID); + + // Still required for lifetime and scalar calls which need to resolve TreeID->volume_tlas on CPU side. + auto result = tree_to_volume_tlas_.emplace(tree, std::move(volume_tlas)); + auto it = result.first; + + // Keep a dense host-side MeshID -> TLAS device-data table for prepared queries. + // The TLAS object in tree_to_volume_tlas_ owns the device allocations; this table + // only stores lightweight DD views indexed by volume ID. Upload to device once in + // init(), unless a volume is registered after initialization. + + const auto volume_index = static_cast(volume_id); + if (volume_index >= volume_to_tlas_.size()) { + volume_to_tlas_.resize(volume_index + 1); + } + volume_to_tlas_[volume_index] = it->second.get_device_data(); + + if (initialized_) { + upload_volume_to_tlas_table_(); + } + + return tree; +} + +TreeID +CuBQLRayTracer::create_element_tree(const std::shared_ptr&, + MeshID) +{ + warning("Element trees not currently supported with cuBQL ray tracer"); + return TREE_NONE; +} + +void CuBQLRayTracer::create_global_surface_tree() +{ + warning("Global surface trees not currently supported with cuBQL ray tracer"); +} + +void CuBQLRayTracer::create_global_element_tree() +{ + warning("Global element trees not currently supported with cuBQL ray tracer"); +} + +MeshID CuBQLRayTracer::find_element(const Position&) const +{ + fatal_error("Element queries not currently supported with cuBQL ray tracer"); + return ID_NONE; +} + +MeshID CuBQLRayTracer::find_element(TreeID, const Position&) const +{ + fatal_error("Element queries not currently supported with cuBQL ray tracer"); + return ID_NONE; +} + +bool CuBQLRayTracer::point_in_volume(TreeID tree, + const Position& point, + const Direction* direction, + const std::vector* exclude_primitives) const +{ + const auto& context = context_; + const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); + + // Use provided direction or if Direction == nulptr use default direction + Direction directionUsed = (direction != nullptr) ? Direction{direction->x, direction->y, direction->z} + : Direction{1. / std::sqrt(2.0), 1. / std::sqrt(2.0), 0.0}; + + + CuBQLRay ray; + ray.origin = cuBQL::vec3d(point.x, point.y, point.z); + ray.direction = cuBQL::vec3d(directionUsed.x, directionUsed.y, directionUsed.z); + ray.tMin = 0.0; + ray.tMax = INFTY; + + CuBQLSurfaceHit surface_hit; + + // TODO - Maybe we can come up with a better name for this + intersect_surface_tree_scalar(context, volume_tlas, ray, surface_hit, HitOrientation::ANY, exclude_primitives); + + // if the ray hit nothing the point must be outside the volume + if (surface_hit.primitive == ID_NONE) return false; + + return surface_hit.piv == INSIDE; +} + +std::pair +CuBQLRayTracer::ray_fire(TreeID tree, + const Position& origin, + const Direction& direction, + const double tmax, + HitOrientation hitOrientation, + std::vector* const exclude_primitives) +{ + const auto& context = context_; + const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); + + CuBQLRay ray; + ray.origin = cuBQL::vec3d(origin.x, origin.y, origin.z); + ray.direction = cuBQL::vec3d(direction.x, direction.y, direction.z); + ray.tMin = 0.0; + ray.tMax = tmax; + + CuBQLSurfaceHit surface_hit; + + // TODO - Maybe we can come up with a better name for this + intersect_surface_tree_scalar(context, volume_tlas, ray, surface_hit, hitOrientation, exclude_primitives); + + if (surface_hit.primitive == ID_NONE) { + return {INFTY, ID_NONE}; + } + + if (exclude_primitives) { + exclude_primitives->push_back(surface_hit.primitive); + } + + return {surface_hit.distance, surface_hit.surface}; +} + +void +CuBQLRayTracer::ray_fire_batch(const CuBQLRay* d_rays, + CuBQLSurfaceHit* d_hits, + std::size_t num_rays, + HitOrientation orientation) +{ + if (num_rays == 0) return; + + if (!d_volume_to_tlas_) { + fatal_error("cuBQL volume TLAS lookup table has not been uploaded"); + } + + intersect_surface_tree_batch(context_, + d_volume_to_tlas_, + d_rays, + d_hits, + num_rays, + orientation); +} + +std::pair +CuBQLRayTracer::closest(TreeID, const Position&) +{ + fatal_error("Closest queries not currently supported with cuBQL ray tracer"); + return {INFTY, ID_NONE}; +} + +bool CuBQLRayTracer::occluded(TreeID, + const Position&, + const Direction&, + double&) const +{ + fatal_error("Occlusion queries not currently supported with cuBQL ray tracer"); + return false; +} + +} // namespace xdg diff --git a/src/cuBQL/triangles.cpp b/src/cuBQL/triangles.cpp new file mode 100644 index 00000000..73cf1fe3 --- /dev/null +++ b/src/cuBQL/triangles.cpp @@ -0,0 +1,52 @@ +#include "xdg/cuBQL/triangles.h" + +#include + +namespace xdg { + +void CuBQLSurfaceMesh::release() +{ + if (d_vertices) { + omp_target_free(d_vertices, gpu_id); + d_vertices = nullptr; + } + if (d_indices) { + omp_target_free(d_indices, gpu_id); + d_indices = nullptr; + } + if (d_primitive_refs) { + omp_target_free(d_primitive_refs, gpu_id); + d_primitive_refs = nullptr; + } +} + +void CuBQLSurfaceBLAS::release() +{ + if (bvh.primIDs) { + omp_target_free(bvh.primIDs, gpu_id); + bvh.primIDs = nullptr; + } + if (bvh.nodes) { + omp_target_free(bvh.nodes, gpu_id); + bvh.nodes = nullptr; + } + mesh.release(); +} + +void CuBQLVolumeTLAS::release() +{ + if (bvh.primIDs) { + omp_target_free(bvh.primIDs, gpu_id); + bvh.primIDs = nullptr; + } + if (bvh.nodes) { + omp_target_free(bvh.nodes, gpu_id); + bvh.nodes = nullptr; + } + if (d_surface_instances) { + omp_target_free(d_surface_instances, gpu_id); + d_surface_instances = nullptr; + } +} + +} // namespace xdg diff --git a/src/xdg.cpp b/src/xdg.cpp index cc738008..18737bf2 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -31,6 +31,14 @@ XDG::XDG(std::shared_ptr mesh_manager, RTLibrary ray_tracing_lib) #else fatal_error("This build was not compiled with GPRT support (XDG_ENABLE_GPRT=OFF)."); #endif + + case RTLibrary::CUBQL: + #ifdef XDG_ENABLE_CUBQL + set_ray_tracing_interface(std::make_shared()); + break; + #else + fatal_error("This build was not compiled with cuBQL support (XDG_ENABLE_CUBQL=OFF)."); + #endif } } @@ -84,6 +92,9 @@ std::shared_ptr XDG::create(MeshLibrary mesh_lib, RTLibrary ray_tracing_lib #ifdef XDG_ENABLE_GPRT if (ray_tracing_lib == RTLibrary::GPRT) return std::make_shared(); #endif + #ifdef XDG_ENABLE_CUBQL + if (ray_tracing_lib == RTLibrary::CUBQL) return std::make_shared(); + #endif // If no supported ray tracing library throw an error std::string msg = fmt::format("Invalid ray tracing library '{}'. Supported:", RT_LIB_TO_STR.at(ray_tracing_lib)); @@ -93,6 +104,9 @@ std::shared_ptr XDG::create(MeshLibrary mesh_lib, RTLibrary ray_tracing_lib #ifdef XDG_ENABLE_GPRT msg += " GPRT"; #endif + #ifdef XDG_ENABLE_CUBQL + msg += " CUBQL"; + #endif fatal_error(msg); }; From 77cb6e4442ea2ab717d13c4f1178dac2d0f8d098 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 14:57:42 +0100 Subject: [PATCH 03/23] Updated vk_device_probe to also probe for available openmp device --- ...ulkan_probe.h => available_device_probe.h} | 48 ++++++++++++++++++- src/gprt/ray_tracer.cpp | 2 +- 2 files changed, 48 insertions(+), 2 deletions(-) rename include/xdg/{gprt/vulkan_probe.h => available_device_probe.h} (64%) diff --git a/include/xdg/gprt/vulkan_probe.h b/include/xdg/available_device_probe.h similarity index 64% rename from include/xdg/gprt/vulkan_probe.h rename to include/xdg/available_device_probe.h index 8e0ff9cf..8aa6c5a1 100644 --- a/include/xdg/gprt/vulkan_probe.h +++ b/include/xdg/available_device_probe.h @@ -1,8 +1,12 @@ #pragma once -#ifdef XDG_ENABLE_GPRT #include "xdg/error.h" +// -------------------------------------------------------------------------------------- +// Vulkan probe functions to check for ray tracing capable devices at runtime +// -------------------------------------------------------------------------------------- +#ifdef XDG_ENABLE_GPRT + #include #include #include @@ -90,5 +94,47 @@ inline bool system_has_vk_device() missing); return false; } +#endif + +// -------------------------------------------------------------------------------------- +// OpenMP target probe functions to check for devices capable of running cuBQL at runtime +// -------------------------------------------------------------------------------------- + +#ifdef XDG_ENABLE_CUBQL + +#include + +inline bool system_has_omp_target_device() +{ + const int device_count = omp_get_num_devices(); + if (device_count <= 0) { + warning("No OpenMP target devices found; cuBQL ray tracer unavailable."); + return false; + } + + const int host_id = omp_get_initial_device(); + for (int device_id = 0; device_id < device_count; ++device_id) { + int value = 1; + void* d_value = omp_target_alloc(sizeof(int), device_id); + if (!d_value) continue; + + const int copy_result = omp_target_memcpy(d_value, + &value, + sizeof(int), + 0, + 0, + device_id, + host_id); + omp_target_free(d_value, device_id); + + if (copy_result == 0) { + write_message("Found OpenMP target device {}.", device_id); + return true; + } + } + + warning("OpenMP target devices were found, but none accepted target allocation; cuBQL ray tracer unavailable."); + return false; +} #endif diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index 15d21a7b..9db035ec 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -1,6 +1,6 @@ #include "xdg/gprt/ray_tracer.h" #include "gprt/gprt.h" -#include "xdg/gprt/vulkan_probe.h" +#include "xdg/available_device_probe.h" namespace xdg { From 346caf12d567573331e3346b03e2010525607751 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:07:32 +0100 Subject: [PATCH 04/23] Add CMakePresets.json for openmp target offload flags --- CMakePresets.json | 68 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 CMakePresets.json diff --git a/CMakePresets.json b/CMakePresets.json new file mode 100644 index 00000000..09b86db1 --- /dev/null +++ b/CMakePresets.json @@ -0,0 +1,68 @@ +{ + "version": 3, + "configurePresets": [ + { + "name": "base", + "generator": "Unix Makefiles", + "binaryDir": "${sourceDir}/build/${presetName}", + "cacheVariables": { + "CMAKE_CXX_FLAGS": "$env{PRESET_CXX_FLAGS} $env{LLVM_CXX_FLAGS} $env{COMMON_CXX_FLAGS}", + "XDG_CMAKE_PRESET": "${presetName}" + }, + "environment": { + "COMMON_CXX_FLAGS": "" + } + }, + { + "name": "llvm", + "inherits": ["base"], + "cacheVariables": { + "CMAKE_C_COMPILER": "clang", + "CMAKE_CXX_COMPILER": "clang++" + }, + "environment": { + "LLVM_CXX_FLAGS": "-fopenmp -fopenmp-cuda-mode" + } + }, + { + "name": "cubql_llvm_ada", + "inherits": ["llvm"], + "displayName": "cuBQL LLVM OpenMP offload RTX 2000 Ada", + "cacheVariables": { + "XDG_ENABLE_CUBQL": "ON" + }, + "environment": { + "PRESET_CXX_FLAGS": "-fopenmp-targets=nvptx64 -Xopenmp-target -march=sm_80" + } + }, + { + "name": "nvhpc", + "inherits": ["base"], + "cacheVariables": { + "CMAKE_C_COMPILER": "nvc", + "CMAKE_CXX_COMPILER": "nvc++" + } + }, + { + "name": "cubql_nvhpc_ada", + "inherits": ["nvhpc"], + "displayName": "cuBQL NVHPC OpenMP offload RTX 2000 Ada", + "cacheVariables": { + "XDG_ENABLE_CUBQL": "ON" + }, + "environment": { + "PRESET_CXX_FLAGS": "-mp=gpu -Minfo=mp -gpu=cc89" + } + } + ], + "buildPresets": [ + { + "name": "cubql_llvm_ada", + "configurePreset": "cubql_llvm_ada" + }, + { + "name": "cubql_nvhpc_ada", + "configurePreset": "cubql_nvhpc_ada" + } + ] +} From b9ebba228f00dccae3f1774b11517f0f40db4e74 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:08:50 +0100 Subject: [PATCH 05/23] Added tests for cuBQL backend --- tests/CMakeLists.txt | 3 ++ tests/test_point_in_volume.cpp | 3 +- tests/test_ray_fire.cpp | 3 +- tests/test_ray_tracer_cross_check.cpp | 57 ++++++++++++++++++++++++++- tests/util.h | 16 +++++++- 5 files changed, 77 insertions(+), 5 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 21691467..bdb1f7e8 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -54,6 +54,9 @@ foreach(test ${TEST_NAMES}) if (XDG_ENABLE_MOAB) target_link_libraries(${test} PRIVATE MOAB) endif() + if (XDG_ENABLE_CUBQL) + target_link_libraries(${test} PRIVATE $) + endif() set_target_properties(${test} PROPERTIES BUILD_RPATH "$") catch_discover_tests(${test} diff --git a/tests/test_point_in_volume.cpp b/tests/test_point_in_volume.cpp index ae34e823..65874f9f 100644 --- a/tests/test_point_in_volume.cpp +++ b/tests/test_point_in_volume.cpp @@ -16,7 +16,8 @@ using namespace xdg::test; TEMPLATE_TEST_CASE("Point-in-volume on MeshMock", "[piv][mock]", Embree_Raytracer, - GPRT_Raytracer) + GPRT_Raytracer, + CuBQL_Raytracer) { constexpr auto rt_backend = TestType::value; diff --git a/tests/test_ray_fire.cpp b/tests/test_ray_fire.cpp index be816c36..70a5353f 100644 --- a/tests/test_ray_fire.cpp +++ b/tests/test_ray_fire.cpp @@ -18,7 +18,8 @@ using namespace xdg::test; TEMPLATE_TEST_CASE("Ray Fire on MeshMock (per-backend sections)", "[rayfire][mock]", Embree_Raytracer, - GPRT_Raytracer) + GPRT_Raytracer, + CuBQL_Raytracer) { // Generate one test run per enabled backend constexpr auto rt_backend = TestType::value; diff --git a/tests/test_ray_tracer_cross_check.cpp b/tests/test_ray_tracer_cross_check.cpp index c75b0273..d2407edf 100644 --- a/tests/test_ray_tracer_cross_check.cpp +++ b/tests/test_ray_tracer_cross_check.cpp @@ -57,7 +57,7 @@ TEST_CASE("Test Pincell RT libraries Cross-Check ray_fire queries", "[moab][rayf } srand48(12345); // set fixed seed for rng - std::vector directions(1000); + std::vector directions(100); for (auto &dir : directions) { dir = rand_dir(); } @@ -86,4 +86,59 @@ TEST_CASE("Test Pincell RT libraries Cross-Check ray_fire queries", "[moab][rayf } } +TEST_CASE("Test Pincell RT libraries Cross-Check point_in_volume queries", "[moab][piv][cross-check]") +{ + const auto rt_cases = make_rt_cases("pincell.h5m"); + if (rt_cases.size() < 2) { + SKIP("Fewer than two ray tracing backends are available; skipping cross-check."); + } + + const std::array directions { + Direction {1.0, 0.0, 0.0}, // axis-aligned x ray + Direction {0.0, 1.0, 0.0}, // axis-aligned y ray + Direction {0.0, 0.0, 1.0}, // axis-aligned z ray + Direction {0.371390676, 0.557086014, 0.742781353} // non-axis ray + }; + + // Points generated by codex with a h5dump to test various edge cases around the pincell geometry + const std::array points { + Position {0.0, 0.0, 0.0}, // pincell center + Position {8.5, 0.25, 0.0}, // inside inner cylinder, r < 9 + Position {9.5, 0.25, 0.0}, // between cylinder radii, 9 < r < 10 + Position {10.5, 0.25, 0.0}, // outside outer cylinder, r > 10 + Position {20.0, 0.25, 0.0}, // inside square cell away from cylinder + Position {24.5, 0.25, 0.0}, // just inside x = 25 square boundary + Position {25.5, 0.25, 0.0}, // between x = 25 and outer x = 27.5 + Position {27.0, 0.25, 0.0}, // just inside outer x = 27.5 boundary + Position {28.0, 0.25, 0.0}, // outside outer x boundary + Position {0.25, 0.0, 19.5}, // just inside positive z = 20 cap + Position {0.25, 0.0, 20.5}, // just outside positive z = 20 cap + Position {0.25, 0.0, -19.5}, // just inside negative z = -20 cap + Position {0.25, 0.0, -20.5}, // just outside negative z = -20 cap + Position {0.25, 24.5, 0.0}, // just inside y = 25 square boundary + Position {0.25, 28.0, 0.0} // outside outer y boundary + }; + + const auto& reference_case = rt_cases.front(); + + for (const auto& volume : reference_case.xdg->mesh_manager()->volumes()) { + for (const auto& point : points) { + for (const auto& direction : directions) { + const auto reference_result = + reference_case.xdg->point_in_volume(volume, point, &direction); + + for (size_t i = 1; i < rt_cases.size(); ++i) { + const auto& candidate = rt_cases[i]; + const auto candidate_result = + candidate.xdg->point_in_volume(volume, point, &direction); + + CAPTURE(volume, point, direction, reference_case.name, candidate.name, + reference_result, candidate_result); + REQUIRE(candidate_result == reference_result); + } + } + } + } +} + // TODO - Add all of the other queries diff --git a/tests/util.h b/tests/util.h index a9841086..96946aa2 100644 --- a/tests/util.h +++ b/tests/util.h @@ -8,7 +8,7 @@ #include "xdg/constants.h" #include "xdg/ray_tracers.h" #include "xdg/mesh_managers.h" -#include "xdg/gprt/vulkan_probe.h" +#include "xdg/available_device_probe.h" namespace xdg::test { @@ -17,7 +17,7 @@ using LibMesh_Interface = std::integral_constant; using GPRT_Raytracer = std::integral_constant; - +using CuBQL_Raytracer = std::integral_constant; } // namespace xdg::test namespace Catch { @@ -49,6 +49,13 @@ inline bool ray_tracer_available(xdg::RTLibrary rt) { #else return false; #endif + + case xdg::RTLibrary::CUBQL: + #ifdef XDG_ENABLE_CUBQL + return system_has_omp_target_device(); + #else + return false; + #endif } return false; @@ -118,5 +125,10 @@ create_raytracer(xdg::RTLibrary rt) { return std::make_shared(); #endif + #ifdef XDG_ENABLE_CUBQL + if (rt == xdg::RTLibrary::CUBQL) + return std::make_shared(); + #endif + return nullptr; } From 428626d333149cfa32075e2790b0c1b92eec39d3 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:09:38 +0100 Subject: [PATCH 06/23] Added tools setup for cuBQL + updated ray_fire tool --- tools/CMakeLists.txt | 4 ++++ tools/ray_fire.cpp | 6 ++++-- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 50a0c650..883cbe84 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -38,3 +38,7 @@ if (TARGET overlap-check) ${CMAKE_CURRENT_SOURCE_DIR}/overlap.cpp ) endif() + +if (XDG_ENABLE_CUBQL AND TARGET ray-benchmark) + target_link_libraries(ray-benchmark PRIVATE $) +endif() diff --git a/tools/ray_fire.cpp b/tools/ray_fire.cpp index c75bbe82..d827b335 100644 --- a/tools/ray_fire.cpp +++ b/tools/ray_fire.cpp @@ -40,7 +40,7 @@ int main(int argc, char** argv) { .default_value("MOAB"); args.add_argument("-r", "--rt-library") - .help("Ray tracing library to use. One of (EMBREE, GPRT)") + .help("Ray tracing library to use. One of (EMBREE, GPRT, CUBQL)") .default_value("EMBREE"); try { @@ -60,6 +60,8 @@ if (rt_str == "EMBREE") rt_lib = RTLibrary::EMBREE; else if (rt_str == "GPRT") rt_lib = RTLibrary::GPRT; +else if (rt_str == "CUBQL") + rt_lib = RTLibrary::CUBQL; else fatal_error("Invalid ray tracing library '{}' specified", rt_str); @@ -68,7 +70,7 @@ if (mesh_str == "MOAB") mesh_lib = MeshLibrary::MOAB; else if (mesh_str == "LIBMESH") { mesh_lib = MeshLibrary::LIBMESH; - if (rt_lib == RTLibrary::GPRT) + if (rt_lib == RTLibrary::GPRT || rt_lib == RTLibrary::CUBQL) fatal_error("LibMesh is not currently supported with GPRT"); } else From e2eff1e117dc0ebd85d4eec26ae3a63abeca99e7 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:10:13 +0100 Subject: [PATCH 07/23] Updated particle sim tool to work with cuBQL backend --- tools/particle_sim.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp index e4460d3f..af3dcf9c 100644 --- a/tools/particle_sim.cpp +++ b/tools/particle_sim.cpp @@ -48,7 +48,7 @@ args.add_argument("-m", "--mesh-library") .default_value("MOAB"); args.add_argument("-r", "--rt-library") - .help("Ray tracing library to use. One of (EMBREE, GPRT)") + .help("Ray tracing library to use. One of (EMBREE, GPRT, CUBQL)") .default_value("EMBREE"); try { args.parse_args(argc, argv); @@ -73,6 +73,8 @@ if (rt_str == "EMBREE") rt_lib = RTLibrary::EMBREE; else if (rt_str == "GPRT") rt_lib = RTLibrary::GPRT; +else if (rt_str == "CUBQL") + rt_lib = RTLibrary::CUBQL; else fatal_error("Invalid ray tracing library '{}' specified", rt_str); From 2d588d0fa1903da270cb58f49b61a1534ee13a2e Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:10:33 +0100 Subject: [PATCH 08/23] Updated ray_benchmark tool to work with cuBQL backend --- tools/ray_benchmark.cpp | 151 ++++++++++++++++++++++++++++++---------- 1 file changed, 116 insertions(+), 35 deletions(-) diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index dea28351..847389d7 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -19,6 +19,11 @@ #include "ray_benchmark.h" +#ifdef XDG_ENABLE_CUBQL +#include +#include "xdg/cuBQL/intersection.h" +#include "xdg/cuBQL/ray_tracer.h" +#endif using namespace xdg; @@ -60,7 +65,7 @@ int main(int argc, char** argv) .default_value("MOAB"); args.add_argument("-rt", "--rt-library") - .help("Ray tracing library to use. Currently implemented: EMBREE") + .help("Ray tracing library to use. Currently implemented: EMBREE, CUBQL") .default_value("EMBREE"); args.add_argument("-l", "--list") @@ -98,6 +103,8 @@ int main(int argc, char** argv) RTLibrary rt_lib; if (rt_str == "EMBREE") { rt_lib = RTLibrary::EMBREE; + } else if (rt_str == "CUBQL") { + rt_lib = RTLibrary::CUBQL; } else { fatal_error("Ray tracing library '{}' is not implemented in this benchmark tool yet", rt_str); } @@ -171,53 +178,127 @@ int main(int argc, char** argv) setup_timer.stop(); + const auto num_faces = mesh_manager->num_volume_faces(volume); + std::size_t num_hits = 0; + if (num_rays < 1) fatal_error("Number of rays must be greater than 0"); + if (rt_lib == RTLibrary::EMBREE) { rt_label += " (" + std::to_string(XDGConfig::config().n_threads()) + " CPU threads)"; - } - const auto num_faces = mesh_manager->num_volume_faces(volume); - - - // Generate random rays from source - generation_timer.start(); - std::vector origins(num_rays); - std::vector directions(num_rays); - - #pragma omp parallel for schedule(runtime) - for (std::size_t i = 0; i < num_rays; ++i) { - std::uint32_t state = seed ^ static_cast(i); - auto sample = tools::benchmark::random_spherical_source(origin.x, - origin.y, - origin.z, - state, - source_radius); - origins[i] = Position(sample.position[0], - sample.position[1], - sample.position[2]); - directions[i] = Direction(sample.direction[0], - sample.direction[1], - sample.direction[2]); - } - generation_timer.stop(); + // Generate random rays from source + generation_timer.start(); + std::vector origins(num_rays); + std::vector directions(num_rays); + + #pragma omp parallel for schedule(runtime) + for (std::size_t i = 0; i < num_rays; ++i) { + std::uint32_t state = seed ^ static_cast(i); + auto sample = tools::benchmark::random_spherical_source(origin.x, + origin.y, + origin.z, + state, + source_radius); + origins[i] = Position(sample.position[0], + sample.position[1], + sample.position[2]); + directions[i] = Direction(sample.direction[0], + sample.direction[1], + sample.direction[2]); + } + generation_timer.stop(); - // Trace rays - trace_timer.start(); + std::vector hit_surfaces(num_rays, ID_NONE); - std::size_t num_hits = 0; + trace_timer.start(); + #pragma omp parallel for schedule(runtime) + for (std::size_t i = 0; i < num_rays; ++i) { + hit_surfaces[i] = xdg->ray_fire(volume, origins[i], directions[i]).second; // just return surface id of hit + } + trace_timer.stop(); - #pragma omp parallel for schedule(runtime) reduction(+:num_hits) - for (std::size_t i = 0; i < num_rays; ++i) { - const auto hit = xdg->ray_fire(volume, origins[i], directions[i]); - if (hit.second != ID_NONE) num_hits++; + // Count hits outside of timing region + for (std::size_t i = 0; i < num_rays; ++i) { + if (hit_surfaces[i] != ID_NONE) num_hits++; + } + } + else if (rt_lib == RTLibrary::CUBQL) { + #ifndef XDG_ENABLE_CUBQL + fatal_error("This build was not compiled with cuBQL support (XDG_ENABLE_CUBQL=OFF)."); + #else + auto rti = std::dynamic_pointer_cast(xdg->ray_tracing_interface()); + + // Generate random rays directly on the target device. + generation_timer.start(); + + const int gpu_id = omp_get_default_device(); + + CuBQLRay* d_rays = static_cast( + omp_target_alloc(num_rays * sizeof(CuBQLRay), gpu_id)); + + if (!d_rays) { + fatal_error("Failed to allocate cuBQL ray buffer"); + } + + const double origin_x = origin.x; + const double origin_y = origin.y; + const double origin_z = origin.z; + + #pragma omp target teams distribute parallel for device(gpu_id) is_device_ptr(d_rays) + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + std::uint32_t state = seed ^ static_cast(ray_id); + + auto sample = tools::benchmark::random_spherical_source(origin_x, + origin_y, + origin_z, + state, + source_radius); + + CuBQLRay ray; + ray.origin = cuBQL::vec3d(sample.position[0], + sample.position[1], + sample.position[2]); + ray.direction = cuBQL::vec3d(sample.direction[0], + sample.direction[1], + sample.direction[2]); + ray.tMin = 0.0; + ray.tMax = INFTY; + ray.volume = volume; + + d_rays[ray_id] = ray; + } + + CuBQLSurfaceHit* d_hits = static_cast( + omp_target_alloc(num_rays * sizeof(CuBQLSurfaceHit), gpu_id)); + + if (!d_hits) { + omp_target_free(d_rays, gpu_id); + fatal_error("Failed to allocate cuBQL hit buffer"); + } + + generation_timer.stop(); + + // Trace rays and count hits on the target device. + trace_timer.start(); + rti->ray_fire_batch(d_rays, d_hits, num_rays); + trace_timer.stop(); + + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(d_hits) reduction(+:num_hits) + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + if (d_hits[ray_id].primitive != ID_NONE) num_hits++; + } + + omp_target_free(d_hits, gpu_id); + omp_target_free(d_rays, gpu_id); + #endif } - - trace_timer.stop(); const std::size_t num_misses = num_rays - num_hits; const double hit_fraction = num_rays > 0 ? static_cast(num_hits) / static_cast(num_rays) : 0.0; + const double generation_time = generation_timer.elapsed(); const double trace_time = trace_timer.elapsed(); const double end_to_end_time = generation_time + trace_time; From 8dea065a9d4771e5a7a0af6f854912cc177d0823 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:15:25 +0100 Subject: [PATCH 09/23] Update to latest cuBQL version in submodule --- vendor/cuBQL | 1 + 1 file changed, 1 insertion(+) create mode 160000 vendor/cuBQL diff --git a/vendor/cuBQL b/vendor/cuBQL new file mode 160000 index 00000000..d1bfc3c2 --- /dev/null +++ b/vendor/cuBQL @@ -0,0 +1 @@ +Subproject commit d1bfc3c2e3533c14c647e967a7994c4ef379e52e From 005f0f24503df70e64a1d457c0d5866cf373ba91 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 22 Jun 2026 16:28:38 +0100 Subject: [PATCH 10/23] Reverting back to working cuBQL commit as latest seems to break ray queries at runtime in xdg --- vendor/cuBQL | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vendor/cuBQL b/vendor/cuBQL index d1bfc3c2..5a651b37 160000 --- a/vendor/cuBQL +++ b/vendor/cuBQL @@ -1 +1 @@ -Subproject commit d1bfc3c2e3533c14c647e967a7994c4ef379e52e +Subproject commit 5a651b3787e0cfae123eb427088fa73921e5a5fe From c5dc7b90fba92b84a07691dce5d9aa5afaca107f Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 25 Jun 2026 13:15:33 +0100 Subject: [PATCH 11/23] Added a public facing ray hit buffer object for downstream applications to populate --- include/xdg/cuBQL/intersection.h | 14 +----- include/xdg/cuBQL/ray_tracer.h | 10 ++-- include/xdg/device_ray.h | 31 ++++++++++++ include/xdg/ray_tracing_interface.h | 11 ++++- include/xdg/xdg.h | 18 +++++++ src/cuBQL/intersection.cpp | 51 ++++++++++++++------ src/cuBQL/ray_tracer.cpp | 47 ++++++++++++++---- src/ray_tracing_interface.cpp | 20 +++++++- src/xdg.cpp | 16 ++++++ tools/ray_benchmark.cpp | 75 +++++++++++------------------ 10 files changed, 202 insertions(+), 91 deletions(-) create mode 100644 include/xdg/device_ray.h diff --git a/include/xdg/cuBQL/intersection.h b/include/xdg/cuBQL/intersection.h index 50d95141..42301f1d 100644 --- a/include/xdg/cuBQL/intersection.h +++ b/include/xdg/cuBQL/intersection.h @@ -11,6 +11,7 @@ #include #include "xdg/constants.h" +#include "xdg/device_ray.h" #include "xdg/cuBQL/triangles.h" #include "cuBQL/math/vec.h" @@ -24,16 +25,6 @@ struct CuBQLRay { MeshID volume {ID_NONE}; // volume we are tracing ray against }; -/* POD SurfaceRay struct for external population*/ -// struct CuBQLSurfaceRay { -// double origin[3]; -// double direction[3]; -// uint32_t volume_slot; -// uint32_t enabled; -// const MeshID* exclude_primitives; -// int32_t exclude_count; -// }; - // TODO - Consider whether this is useful/necessary as its own struct // struct CuBQLExcludeList { // const MeshID* primitives {nullptr}; @@ -79,8 +70,7 @@ intersect_surface_tree_scalar(const cubql::Context& context, void intersect_surface_tree_batch(const cubql::Context& context, const CuBQLVolumeTLAS::DD* d_volume_to_tlas, - const CuBQLRay* d_rays, - CuBQLSurfaceHit* d_hits, + XDGRayHit* d_ray_hits, std::size_t num_rays, HitOrientation hit_orientation); diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index c422ecb3..4d1ee84e 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -58,10 +58,12 @@ class CuBQLRayTracer : public RayTracer { HitOrientation orientation = HitOrientation::EXITING, std::vector* const exclude_primitives = nullptr) override; - void ray_fire_batch(const CuBQLRay* d_rays, - CuBQLSurfaceHit* d_hits, - std::size_t num_rays, - HitOrientation orientation = HitOrientation::EXITING); + void ray_fire_batch(const XDGRayHitBuffer& ray_hits, + HitOrientation hit_orientation = HitOrientation::EXITING) const override; + + XDGRayHitBuffer allocate_ray_hits(std::size_t count) const override; + + void free_ray_hits(XDGRayHitBuffer& ray_hits) const override; std::pair closest(TreeID tree, const Position& origin) override; diff --git a/include/xdg/device_ray.h b/include/xdg/device_ray.h new file mode 100644 index 00000000..e3fa5963 --- /dev/null +++ b/include/xdg/device_ray.h @@ -0,0 +1,31 @@ +#ifndef _XDG_DEVICE_RAY_H +#define _XDG_DEVICE_RAY_H + +#include +#include + +namespace xdg { + +struct XDGRayHit { + double origin[3]; + double direction[3]; + double t_min; + double t_max; + std::int32_t volume; + + double distance; + std::int32_t surface; + std::int32_t primitive; + std::int32_t point_in_volume; +}; + +// Light wrapper for count and device id associated with pointer +struct XDGRayHitBuffer { + XDGRayHit* data {nullptr}; + std::size_t count {0}; + int device_id {-1}; +}; + +} // namespace xdg + +#endif // include guard diff --git a/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h index e0cf3cb7..d593cbb9 100644 --- a/include/xdg/ray_tracing_interface.h +++ b/include/xdg/ray_tracing_interface.h @@ -1,11 +1,13 @@ #ifndef _XDG_RAY_TRACING_INTERFACE_H #define _XDG_RAY_TRACING_INTERFACE_H +#include #include #include #include #include "xdg/constants.h" +#include "xdg/device_ray.h" #include "xdg/mesh_manager_interface.h" #include "xdg/primitive_ref.h" #include "xdg/geometry_data.h" @@ -85,6 +87,13 @@ class RayTracer { HitOrientation orientation = HitOrientation::EXITING, std::vector* const exclude_primitives = nullptr) = 0; + virtual XDGRayHitBuffer allocate_ray_hits(std::size_t count) const; + + virtual void free_ray_hits(XDGRayHitBuffer& ray_hits) const; + + virtual void ray_fire_batch(const XDGRayHitBuffer& ray_hits, + HitOrientation hit_orientation = HitOrientation::EXITING) const; + /** * @brief Finds the element containing a given point using the global element tree. * @@ -149,4 +158,4 @@ class RayTracer { } // namespace xdg -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/xdg.h b/include/xdg/xdg.h index 82c1f8f2..bcae8a30 100644 --- a/include/xdg/xdg.h +++ b/include/xdg/xdg.h @@ -1,6 +1,7 @@ #ifndef _XDG_INTERFACE_H #define _XDG_INTERFACE_H +#include #include #include @@ -75,6 +76,23 @@ std::pair ray_fire(MeshID volume, HitOrientation orientation = HitOrientation::EXITING, std::vector* const exclude_primitives = nullptr) const; +//! Allocates a backend-owned device buffer for batch ray-fire records. +//! @param count Number of XDGRayHit records to allocate. +//! @return Device buffer handle containing the pointer, record count, and device id. +//! Release the buffer with free_ray_hits(). +XDGRayHitBuffer allocate_ray_hits(std::size_t count) const; + +//! Releases a device ray-hit buffer allocated by allocate_ray_hits(). +//! @param ray_hits Buffer handle to release. The handle is cleared after release. +void free_ray_hits(XDGRayHitBuffer& ray_hits) const; + +//! Fires all rays stored in a device ray-hit buffer using the selected backend. +//! @param ray_hits Device buffer whose XDGRayHit records have been populated by the caller. +//! Hit result fields are written back into the same records. +//! @param hit_orientation Orientation filter applied to every ray in the batch. +void ray_fire_batch(const XDGRayHitBuffer& ray_hits, + HitOrientation hit_orientation = HitOrientation::EXITING) const; + std::pair closest(MeshID volume, const Position& origin) const; diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index cc5b1bdb..2503eb47 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -173,37 +173,56 @@ intersect_surface_tree_scalar(const cubql::Context& context, void intersect_surface_tree_batch(const cubql::Context& context, const CuBQLVolumeTLAS::DD* d_volume_to_tlas, - const CuBQLRay* d_rays, - CuBQLSurfaceHit* d_hits, + XDGRayHit* d_ray_hits, std::size_t num_rays, HitOrientation hit_orientation) { - if (num_rays == 0) return; - if (!d_volume_to_tlas || !d_rays || !d_hits) { + if (!d_volume_to_tlas || !d_ray_hits) { fatal_error("Invalid cuBQL batch intersection buffers"); } const int gpu_id = context.gpuID; #pragma omp target teams distribute parallel for device(gpu_id) \ - is_device_ptr(d_volume_to_tlas, d_rays, d_hits) + is_device_ptr(d_volume_to_tlas, d_ray_hits) for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { - const CuBQLRay ray = d_rays[ray_id]; - const CuBQLVolumeTLAS::DD volume_tlas = d_volume_to_tlas[ray.volume]; + XDGRayHit ray_hit = d_ray_hits[ray_id]; CuBQLSurfaceHit hit; - hit.distance = ray.tMax; - - intersect_surface_tree(volume_tlas, - ray, - &hit, - static_cast(hit_orientation), - nullptr, - 0); + hit.distance = ray_hit.t_max; + hit.surface = ID_NONE; + hit.primitive = ID_NONE; + hit.piv = OUTSIDE; + + if (ray_hit.volume != ID_NONE) { + CuBQLRay ray; + ray.origin = cuBQL::vec3d(ray_hit.origin[0], + ray_hit.origin[1], + ray_hit.origin[2]); + ray.direction = cuBQL::vec3d(ray_hit.direction[0], + ray_hit.direction[1], + ray_hit.direction[2]); + ray.tMin = ray_hit.t_min; + ray.tMax = ray_hit.t_max; + ray.volume = ray_hit.volume; + + const CuBQLVolumeTLAS::DD volume_tlas = d_volume_to_tlas[ray.volume]; + + intersect_surface_tree(volume_tlas, + ray, + &hit, + static_cast(hit_orientation), + nullptr, + 0); + } - d_hits[ray_id] = hit; + ray_hit.distance = hit.distance; + ray_hit.surface = hit.surface; + ray_hit.primitive = hit.primitive; + ray_hit.point_in_volume = static_cast(hit.piv); + d_ray_hits[ray_id] = ray_hit; } } diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index c41bc4f2..21e0c5e0 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -398,13 +398,43 @@ CuBQLRayTracer::ray_fire(TreeID tree, return {surface_hit.distance, surface_hit.surface}; } +XDGRayHitBuffer CuBQLRayTracer::allocate_ray_hits(std::size_t count) const +{ + if (count == 0) { + warning("Request to allocate 0 cuBQL XDG ray-hit buffer; returning empty buffer"); + return {}; + } + + auto* d_ray_hits = static_cast + (omp_target_alloc(count * sizeof(XDGRayHit), context_.gpuID)); + + if (!d_ray_hits) { + fatal_error("Failed to allocate cuBQL XDG ray-hit buffer"); + } + + return {d_ray_hits, count, context_.gpuID}; +} + +void CuBQLRayTracer::free_ray_hits(XDGRayHitBuffer& ray_hits) const +{ + if (!ray_hits.data) { + warning("Request to free empty cuBQL XDG ray-hit buffer; ignoring"); + return; + } + + omp_target_free(ray_hits.data, ray_hits.device_id); + ray_hits = {}; +} + void -CuBQLRayTracer::ray_fire_batch(const CuBQLRay* d_rays, - CuBQLSurfaceHit* d_hits, - std::size_t num_rays, - HitOrientation orientation) +CuBQLRayTracer::ray_fire_batch(const XDGRayHitBuffer& ray_hits, + HitOrientation hit_orientation) const { - if (num_rays == 0) return; + if (ray_hits.count == 0) return; + + if (!ray_hits.data) { + fatal_error("Invalid cuBQL XDG ray-hit buffer"); + } if (!d_volume_to_tlas_) { fatal_error("cuBQL volume TLAS lookup table has not been uploaded"); @@ -412,10 +442,9 @@ CuBQLRayTracer::ray_fire_batch(const CuBQLRay* d_rays, intersect_surface_tree_batch(context_, d_volume_to_tlas_, - d_rays, - d_hits, - num_rays, - orientation); + ray_hits.data, + ray_hits.count, + hit_orientation); } std::pair diff --git a/src/ray_tracing_interface.cpp b/src/ray_tracing_interface.cpp index 66e0adb8..fe9fb6eb 100644 --- a/src/ray_tracing_interface.cpp +++ b/src/ray_tracing_interface.cpp @@ -1,5 +1,6 @@ #include #include "xdg/ray_tracing_interface.h" +#include "xdg/error.h" // Any methods which are identical for all RT backends should be defined here @@ -7,6 +8,23 @@ namespace xdg { RayTracer::~RayTracer() {} +XDGRayHitBuffer RayTracer::allocate_ray_hits(std::size_t) const +{ + fatal_error("Selected ray tracer does not support device batch ray fire"); + return {}; +} + +void RayTracer::free_ray_hits(XDGRayHitBuffer&) const +{ + fatal_error("Selected ray tracer does not support device batch ray fire"); +} + +void RayTracer::ray_fire_batch(const XDGRayHitBuffer&, + HitOrientation) const +{ + fatal_error("Selected ray tracer does not support device batch ray fire"); +} + SurfaceTreeID RayTracer::next_surface_tree_id() { return ++next_surface_tree_id_; @@ -23,4 +41,4 @@ const double RayTracer::bounding_box_bump(const std::shared_ptr mes return std::max(volume_bounding_box.dilation(), numerical_precision_); } -} // namespace xdg \ No newline at end of file +} // namespace xdg diff --git a/src/xdg.cpp b/src/xdg.cpp index 18737bf2..50ae3c5b 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -260,6 +260,22 @@ XDG::ray_fire(MeshID volume, return ray_tracing_interface()->ray_fire(scene, origin, direction, dist_limit, orientation, exclude_primitives); } +XDGRayHitBuffer XDG::allocate_ray_hits(std::size_t count) const +{ + return ray_tracing_interface()->allocate_ray_hits(count); +} + +void XDG::free_ray_hits(XDGRayHitBuffer& ray_hits) const +{ + ray_tracing_interface()->free_ray_hits(ray_hits); +} + +void XDG::ray_fire_batch(const XDGRayHitBuffer& ray_hits, + HitOrientation hit_orientation) const +{ + ray_tracing_interface()->ray_fire_batch(ray_hits, hit_orientation); +} + std::pair XDG::closest(MeshID volume, const Position& origin) const { diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index 847389d7..002cdbf7 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -19,12 +19,6 @@ #include "ray_benchmark.h" -#ifdef XDG_ENABLE_CUBQL -#include -#include "xdg/cuBQL/intersection.h" -#include "xdg/cuBQL/ray_tracer.h" -#endif - using namespace xdg; int main(int argc, char** argv) @@ -223,29 +217,20 @@ int main(int argc, char** argv) } } else if (rt_lib == RTLibrary::CUBQL) { - #ifndef XDG_ENABLE_CUBQL - fatal_error("This build was not compiled with cuBQL support (XDG_ENABLE_CUBQL=OFF)."); - #else - auto rti = std::dynamic_pointer_cast(xdg->ray_tracing_interface()); - // Generate random rays directly on the target device. generation_timer.start(); - const int gpu_id = omp_get_default_device(); - - CuBQLRay* d_rays = static_cast( - omp_target_alloc(num_rays * sizeof(CuBQLRay), gpu_id)); - - if (!d_rays) { - fatal_error("Failed to allocate cuBQL ray buffer"); - } + XDGRayHitBuffer ray_hits = xdg->allocate_ray_hits(num_rays); + XDGRayHit* d_ray_hits = ray_hits.data; + const std::size_t ray_count = ray_hits.count; + const int gpu_id = ray_hits.device_id; const double origin_x = origin.x; const double origin_y = origin.y; const double origin_z = origin.z; - #pragma omp target teams distribute parallel for device(gpu_id) is_device_ptr(d_rays) - for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + #pragma omp target teams distribute parallel for device(gpu_id) is_device_ptr(d_ray_hits) + for (std::size_t ray_id = 0; ray_id < ray_count; ++ray_id) { std::uint32_t state = seed ^ static_cast(ray_id); auto sample = tools::benchmark::random_spherical_source(origin_x, @@ -254,44 +239,38 @@ int main(int argc, char** argv) state, source_radius); - CuBQLRay ray; - ray.origin = cuBQL::vec3d(sample.position[0], - sample.position[1], - sample.position[2]); - ray.direction = cuBQL::vec3d(sample.direction[0], - sample.direction[1], - sample.direction[2]); - ray.tMin = 0.0; - ray.tMax = INFTY; - ray.volume = volume; - - d_rays[ray_id] = ray; - } - - CuBQLSurfaceHit* d_hits = static_cast( - omp_target_alloc(num_rays * sizeof(CuBQLSurfaceHit), gpu_id)); - - if (!d_hits) { - omp_target_free(d_rays, gpu_id); - fatal_error("Failed to allocate cuBQL hit buffer"); + XDGRayHit ray_hit; + ray_hit.origin[0] = sample.position[0]; + ray_hit.origin[1] = sample.position[1]; + ray_hit.origin[2] = sample.position[2]; + ray_hit.direction[0] = sample.direction[0]; + ray_hit.direction[1] = sample.direction[1]; + ray_hit.direction[2] = sample.direction[2]; + ray_hit.t_min = 0.0; + ray_hit.t_max = INFTY; + ray_hit.volume = volume; + ray_hit.distance = INFTY; + ray_hit.surface = ID_NONE; + ray_hit.primitive = ID_NONE; + ray_hit.point_in_volume = OUTSIDE; + + d_ray_hits[ray_id] = ray_hit; } generation_timer.stop(); // Trace rays and count hits on the target device. trace_timer.start(); - rti->ray_fire_batch(d_rays, d_hits, num_rays); + xdg->ray_fire_batch(ray_hits); trace_timer.stop(); #pragma omp target teams distribute parallel for device(gpu_id) \ - is_device_ptr(d_hits) reduction(+:num_hits) - for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { - if (d_hits[ray_id].primitive != ID_NONE) num_hits++; + is_device_ptr(d_ray_hits) reduction(+:num_hits) + for (std::size_t ray_id = 0; ray_id < ray_count; ++ray_id) { + if (d_ray_hits[ray_id].surface != ID_NONE) num_hits++; } - omp_target_free(d_hits, gpu_id); - omp_target_free(d_rays, gpu_id); - #endif + xdg->free_ray_hits(ray_hits); } const std::size_t num_misses = num_rays - num_hits; From 80ac8e09c4e41e173b735c5a25286478f1ec56a7 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 25 Jun 2026 17:19:47 +0100 Subject: [PATCH 12/23] Minor optimizations and code clarity changes to cuBQL traversal kernel --- src/cuBQL/intersection.cpp | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index 2503eb47..eb394d79 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -11,6 +11,11 @@ namespace xdg { // Core traversal and intersection routine for a single ray against a given volume tlas #pragma omp declare target +static inline float reject_candidate(const cuBQL::ray3f& traversal_ray) +{ + return traversal_ray.tMax; +} + static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, CuBQLRay intersection_ray, CuBQLSurfaceHit* hit, @@ -45,7 +50,7 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, for (int i = 0; i < exclude_count; ++i) { if (exclude_primitives[i] == primitive_ref) { - return traversal_ray.tMax; + return reject_candidate(traversal_ray); } } @@ -60,15 +65,15 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, cuBQL::vec3d normal = cuBQL::cross(vertices[1] - vertices[0], vertices[2] - vertices[0]); + double normal_dot_direction = dot(normal, intersection_ray.direction); + if (surface_instance.reverse_sense) { - normal = -normal; + normal_dot_direction = -normal_dot_direction; } - const double normal_dot_direction = dot(normal, intersection_ray.direction); - if (orientation_cull(normal_dot_direction, static_cast(orientation))) { - return traversal_ray.tMax; + return reject_candidate(traversal_ray); } auto intersection = plucker_ray_tri_intersect(vertices, @@ -89,7 +94,7 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, // Return value is only the FP32 traversal shrink distance. The accepted hit // distance stored above remains the FP64 Plucker result. - return traversal_ray.tMax; + return reject_candidate(traversal_ray); }; auto leave_blas = []() -> void {}; @@ -188,7 +193,7 @@ intersect_surface_tree_batch(const cubql::Context& context, #pragma omp target teams distribute parallel for device(gpu_id) \ is_device_ptr(d_volume_to_tlas, d_ray_hits) for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { - XDGRayHit ray_hit = d_ray_hits[ray_id]; + const XDGRayHit ray_hit = d_ray_hits[ray_id]; CuBQLSurfaceHit hit; hit.distance = ray_hit.t_max; @@ -218,11 +223,10 @@ intersect_surface_tree_batch(const cubql::Context& context, 0); } - ray_hit.distance = hit.distance; - ray_hit.surface = hit.surface; - ray_hit.primitive = hit.primitive; - ray_hit.point_in_volume = static_cast(hit.piv); - d_ray_hits[ray_id] = ray_hit; + d_ray_hits[ray_id].distance = hit.distance; + d_ray_hits[ray_id].surface = hit.surface; + d_ray_hits[ray_id].primitive = hit.primitive; + d_ray_hits[ray_id].point_in_volume = static_cast(hit.piv); } } From 5c9569f40944c57cda5a710ffcac3bb30d6f697f Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 14 Jul 2026 11:51:35 +0100 Subject: [PATCH 13/23] Updated ray payload for ray_fire_batch to include surface crossing metadata - Boundary condition metadata now stored on surface instances - Next volume topology also stored on surface instances - Primitive normal stored on hits to be used in reflective BCs - Updated rayhit structs throughout to carry this new data --- include/xdg/cuBQL/intersection.h | 3 +++ include/xdg/cuBQL/triangles.h | 2 ++ include/xdg/device_ray.h | 3 +++ include/xdg/shared_enums.h | 9 ++++++++- src/cuBQL/intersection.cpp | 12 +++++++++++- src/cuBQL/ray_tracer.cpp | 20 +++++++++++++++++++- 6 files changed, 46 insertions(+), 3 deletions(-) diff --git a/include/xdg/cuBQL/intersection.h b/include/xdg/cuBQL/intersection.h index 42301f1d..70d596fc 100644 --- a/include/xdg/cuBQL/intersection.h +++ b/include/xdg/cuBQL/intersection.h @@ -36,6 +36,9 @@ struct CuBQLSurfaceHit { MeshID surface {ID_NONE}; MeshID primitive {ID_NONE}; PointInVolume piv {OUTSIDE}; + MeshID next_volume {ID_NONE}; + SurfaceBoundaryCondition boundary_condition {SurfaceBoundaryCondition::UNSET}; + cuBQL::vec3d normal {0.0}; bool hit_found() const { return primitive != ID_NONE; } }; diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h index 20e53c91..87ca1e0d 100644 --- a/include/xdg/cuBQL/triangles.h +++ b/include/xdg/cuBQL/triangles.h @@ -98,6 +98,8 @@ struct CuBQLVolumeTLAS { struct SurfaceInstanceDD { CuBQLSurfaceBLAS::DD surface_blas; bool reverse_sense {false}; // value set in create_surface_tree based on parent vols + MeshID next_volume {ID_NONE}; + SurfaceBoundaryCondition boundary_condition {UNSET}; }; struct DD { diff --git a/include/xdg/device_ray.h b/include/xdg/device_ray.h index e3fa5963..ee57d582 100644 --- a/include/xdg/device_ray.h +++ b/include/xdg/device_ray.h @@ -17,6 +17,9 @@ struct XDGRayHit { std::int32_t surface; std::int32_t primitive; std::int32_t point_in_volume; + std::int32_t next_volume; + std::int32_t boundary_condition; + double normal[3]; }; // Light wrapper for count and device id associated with pointer diff --git a/include/xdg/shared_enums.h b/include/xdg/shared_enums.h index f6198f60..aa0f6fc1 100644 --- a/include/xdg/shared_enums.h +++ b/include/xdg/shared_enums.h @@ -14,6 +14,13 @@ namespace xdg { ENTERING = 1, }; + enum SurfaceBoundaryCondition : int { + UNSET = -1, + TRANSMISSION = 0, // Cross into the next volume + VACUUM = 1, // Kill particle on crossing + REFLECTIVE = 2 // Reflect particle using the surface normal + }; + } -#endif // XDG_SHARED_ENUMS_H \ No newline at end of file +#endif // XDG_SHARED_ENUMS_H diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index eb394d79..edd0d0b8 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -83,12 +83,16 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, intersection_ray.tMin, false, 0); - + + // store ray payload if hit found if (intersection.hit) { hit->distance = intersection.t; hit->surface = mesh.surface_id; hit->primitive = primitive_ref; hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; + hit->next_volume = surface_instance.next_volume; + hit->boundary_condition = surface_instance.boundary_condition; + hit->normal = normal; traversal_ray.tMax = static_cast(intersection.t); } @@ -200,6 +204,7 @@ intersect_surface_tree_batch(const cubql::Context& context, hit.surface = ID_NONE; hit.primitive = ID_NONE; hit.piv = OUTSIDE; + hit.next_volume = ID_NONE; if (ray_hit.volume != ID_NONE) { CuBQLRay ray; @@ -227,6 +232,11 @@ intersect_surface_tree_batch(const cubql::Context& context, d_ray_hits[ray_id].surface = hit.surface; d_ray_hits[ray_id].primitive = hit.primitive; d_ray_hits[ray_id].point_in_volume = static_cast(hit.piv); + d_ray_hits[ray_id].next_volume = hit.next_volume; + d_ray_hits[ray_id].boundary_condition = static_cast(hit.boundary_condition); + d_ray_hits[ray_id].normal[0] = hit.normal.x; + d_ray_hits[ray_id].normal[1] = hit.normal.y; + d_ray_hits[ray_id].normal[2] = hit.normal.z; } } diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 21e0c5e0..f8f64280 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -229,15 +229,33 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; surface_instance.surface_blas = surface_blas.get_device_data(); - // Sense setting for each surface instance in the TLAS + // Store per-instance topology for this volume and surface. if (volume_id == forward_parent) { surface_instance.reverse_sense = false; + surface_instance.next_volume = reverse_parent; } else if (volume_id == reverse_parent) { surface_instance.reverse_sense = true; + surface_instance.next_volume = forward_parent; } else { fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); } + // Store boundary-condition metadata for this surface instance. + const auto property = mesh_manager->get_surface_property( + surf, PropertyType::BOUNDARY_CONDITION); + + if (property.value == "vacuum") { + surface_instance.boundary_condition = VACUUM; + } else if (property.value == "reflecting" || + property.value == "reflective") { + surface_instance.boundary_condition = REFLECTIVE; + } else if (property.value == "transmission") { + surface_instance.boundary_condition = TRANSMISSION; + } else { + fatal_error("Unsupported boundary condition '{}' on surface {}", + property.value, surf); + } + h_tlas_boxes.push_back(surface_bounds); h_surface_instances.push_back(surface_instance); } From 012451e8db266195b7d06b7ad7e0375355a3d21c Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 17 Jul 2026 12:32:55 +0100 Subject: [PATCH 14/23] Switch cuBQL to flattened per-volume BVHs - Replace the two-level BLAS/TLAS scheme with single-level traversal - Introduce CuBQLSurfaceMesh and CuBQLVolumeGroup following DPRT terminology - Build flattened primitive references and AABBs in create_surface_tree - Update scalar and batch intersection paths for flattened traversal --- CMakePresets.json | 2 +- include/xdg/cuBQL/intersection.h | 4 +- include/xdg/cuBQL/ray_tracer.h | 15 +- include/xdg/cuBQL/triangles.h | 97 ++++------ include/xdg/device_ray.h | 1 + src/cuBQL/intersection.cpp | 84 ++++---- src/cuBQL/ray_tracer.cpp | 323 +++++++++++++++---------------- src/cuBQL/triangles.cpp | 37 ++-- tools/ray_benchmark.cpp | 1 + 9 files changed, 273 insertions(+), 291 deletions(-) diff --git a/CMakePresets.json b/CMakePresets.json index 09b86db1..ad89fa82 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -32,7 +32,7 @@ "XDG_ENABLE_CUBQL": "ON" }, "environment": { - "PRESET_CXX_FLAGS": "-fopenmp-targets=nvptx64 -Xopenmp-target -march=sm_80" + "PRESET_CXX_FLAGS": "-fopenmp-targets=nvptx64 -Xopenmp-target -march=sm_89" } }, { diff --git a/include/xdg/cuBQL/intersection.h b/include/xdg/cuBQL/intersection.h index 70d596fc..9a941a4f 100644 --- a/include/xdg/cuBQL/intersection.h +++ b/include/xdg/cuBQL/intersection.h @@ -63,7 +63,7 @@ Performs host side staging and transfer hit data back to host after device side */ void intersect_surface_tree_scalar(const cubql::Context& context, - const CuBQLVolumeTLAS& volume_tlas, + const CuBQLVolumeGroup& volume_group, const CuBQLRay& ray, CuBQLSurfaceHit& hit, HitOrientation hit_orientation, @@ -72,7 +72,7 @@ intersect_surface_tree_scalar(const cubql::Context& context, void intersect_surface_tree_batch(const cubql::Context& context, - const CuBQLVolumeTLAS::DD* d_volume_to_tlas, + const CuBQLVolumeGroup::DD* d_volume_to_group, XDGRayHit* d_ray_hits, std::size_t num_rays, HitOrientation hit_orientation); diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index 4d1ee84e..ebd4bb88 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -74,20 +74,19 @@ class CuBQLRayTracer : public RayTracer { double& dist) const override; private: - CuBQLSurfaceBLAS + CuBQLSurfaceMesh register_surface(const std::shared_ptr& mesh_manager, - MeshID surface_id, - double bounding_box_bump); + MeshID surface_id); - void upload_volume_to_tlas_table_(); + void upload_volume_to_group_table_(); cubql::Context context_; - std::unordered_map tree_to_volume_tlas_; - std::unordered_map surface_to_blas_map_; + std::unordered_map tree_to_volume_group_; + std::unordered_map surface_to_mesh_; - std::vector volume_to_tlas_; - CuBQLVolumeTLAS::DD* d_volume_to_tlas_ {nullptr}; + std::vector volume_to_group_; + CuBQLVolumeGroup::DD* d_volume_to_group_ {nullptr}; bool initialized_ {false}; }; diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h index 87ca1e0d..0f190d4d 100644 --- a/include/xdg/cuBQL/triangles.h +++ b/include/xdg/cuBQL/triangles.h @@ -2,7 +2,6 @@ #define _XDG_CUBQL_TRIANGLES_H #include -#include // Guards to prevent CUDA headers from being included in host code, which causes // failed compilation with LLVM-clang. @@ -17,10 +16,9 @@ namespace xdg { -/* - Owns the triangle buffers for one topological surface. The nested DD type is - the compact device-data view copied into OpenMP target regions instead of the - full host-side owner. +/** + Owns the primitive buffers for one topological surface. The nested DD type is + the compact, non-owning device-data view used in OpenMP target regions. */ struct CuBQLSurfaceMesh { struct DD { @@ -30,99 +28,84 @@ struct CuBQLSurfaceMesh { // Geometric data const cuBQL::vec3d* vertices {nullptr}; const cuBQL::vec3i* indices {nullptr}; - const MeshID* primitive_refs {nullptr}; + const MeshID* primitive_ids {nullptr}; }; // Topological metadata MeshID surface_id {ID_NONE}; - // Device buffers for triangle data + // Device buffers for primitive data cuBQL::vec3d* d_vertices {nullptr}; cuBQL::vec3i* d_indices {nullptr}; - MeshID* d_primitive_refs {nullptr}; + MeshID* d_primitive_ids {nullptr}; - uint32_t num_vertices {0}; - uint32_t num_triangles {0}; + std::uint32_t num_vertices {0}; + std::uint32_t num_primitives {0}; int gpu_id {0}; - // Accessor for Device Data struct, which is passed to cuBQL BVH traversal/intersection functions + // Return the non-owning view used by device traversal and intersection code. DD get_device_data() const { return { surface_id, d_vertices, d_indices, - d_primitive_refs + d_primitive_ids }; } void release(); }; -/* - Owns a cuBQL BVH used as a Bottom Level Acceleration Structure over surface triangles. - The nested DD type is the compact device-data view used during traversal. +/** + Owns the flattened cuBQL BVH for one topological volume. The BVH is built over + all primitives belonging to the volume's surfaces, while PrimRef maps each + BVH primitive back to its surface and surface-local primitive. The nested DD + type is the compact, non-owning device-data view used during traversal. */ -struct CuBQLSurfaceBLAS { - struct DD { - CuBQLSurfaceMesh::DD mesh; // Mesh data device handle - cuBQL::bvh3f bvh; // BLAS device handle - }; - - cuBQL::bvh3f bvh; // BLAS host handle - CuBQLSurfaceMesh mesh; // Surface mesh host owner - - uint32_t num_prims {0}; - int gpu_id {0}; - - DD get_device_data() const - { - return {mesh.get_device_data(), bvh}; - } - - void release(); -}; +struct CuBQLVolumeGroup { + struct SurfaceDD { + CuBQLSurfaceMesh::DD mesh; -/* - Owns a cuBQL BVH used as a Top-Level Acceleration Structure for one topological volume. - The TLAS groups the surface BLASes that bound that volume and stores - per-volume relationship metadata for each surface instance. -*/ -struct CuBQLVolumeTLAS { - /* - TLAS-local instance payload. The same surface BLAS can participate in - different volume TLASes with different sense, so reverse_sense belongs on - the volume-surface relationship rather than on the reusable surface mesh - or BLAS geometry. - */ - struct SurfaceInstanceDD { - CuBQLSurfaceBLAS::DD surface_blas; - bool reverse_sense {false}; // value set in create_surface_tree based on parent vols + bool reverse_sense {false}; MeshID next_volume {ID_NONE}; SurfaceBoundaryCondition boundary_condition {UNSET}; }; + // Identifies a primitive within the volume group's local surface array. + // TODO - Should this exist outside of volume group as an indpendent struct? + // TODO - can we think of a better name to distinguish between CuBQLSurfaceMesh::primitive_ids and this struct? + struct PrimRef { + std::uint32_t surface_index {0}; + std::uint32_t primitive_index {0}; + }; + struct DD { - MeshID volume_id {ID_NONE}; - const SurfaceInstanceDD* surface_instances {nullptr}; - cuBQL::bvh3f bvh; // TLAS device handle + const SurfaceDD* surfaces {nullptr}; + const PrimRef* prim_refs {nullptr}; + cuBQL::bvh3f bvh; }; - MeshID volume_id {ID_NONE}; - cuBQL::bvh3f bvh; // TLAS host handle - SurfaceInstanceDD* d_surface_instances {nullptr}; + cuBQL::bvh3f bvh; + SurfaceDD* d_surfaces {nullptr}; + PrimRef* d_prim_refs {nullptr}; - uint32_t num_surface_instances {0}; + std::uint32_t num_surfaces {0}; + std::uint32_t num_primitives {0}; int gpu_id {0}; + // Return the non-owning view used by device traversal and intersection code. DD get_device_data() const { - return {volume_id, d_surface_instances, bvh}; + return {d_surfaces, d_prim_refs, bvh}; } void release(); }; +// Future acceleration structure over instances of volume groups. +struct CuBQLInstanceGroup; + } // namespace xdg #endif // include guard diff --git a/include/xdg/device_ray.h b/include/xdg/device_ray.h index ee57d582..8b113aa9 100644 --- a/include/xdg/device_ray.h +++ b/include/xdg/device_ray.h @@ -12,6 +12,7 @@ struct XDGRayHit { double t_min; double t_max; std::int32_t volume; + std::int32_t last_hit_primitive {-1}; double distance; std::int32_t surface; diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index edd0d0b8..49a3a151 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -9,17 +9,19 @@ namespace xdg { -// Core traversal and intersection routine for a single ray against a given volume tlas +// Core traversal and intersection routine for a single ray against a flattened +// volume group. #pragma omp declare target static inline float reject_candidate(const cuBQL::ray3f& traversal_ray) { return traversal_ray.tMax; } -static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, +static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, CuBQLRay intersection_ray, CuBQLSurfaceHit* hit, int orientation, + MeshID last_hit_primitive, const MeshID* exclude_primitives, int exclude_count) { @@ -32,34 +34,34 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, traversal_ray.tMin = static_cast(intersection_ray.tMin); traversal_ray.tMax = static_cast(hit->distance); - CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; - - auto enter_blas = [=, &surface_instance, &traversal_ray] - (cuBQL::ray3f& out_ray, cuBQL::bvh3f& out_bvh, int instance_id) - { - surface_instance = volume_tlas.surface_instances[instance_id]; - out_ray = traversal_ray; - out_bvh = surface_instance.surface_blas.bvh; - }; - - auto intersect_prim = [=, &traversal_ray, &surface_instance] - (uint32_t prim_id) -> float + auto intersect_prim = [=, &traversal_ray] + (std::uint32_t bvh_primitive_index) -> float { - const CuBQLSurfaceMesh::DD mesh = surface_instance.surface_blas.mesh; - const MeshID primitive_ref = mesh.primitive_refs[prim_id]; + const auto ref = volume_group.prim_refs[bvh_primitive_index]; + const auto surface = volume_group.surfaces[ref.surface_index]; + const auto mesh = surface.mesh; + const auto local_index = ref.primitive_index; + const MeshID primitive_id = mesh.primitive_ids[local_index]; + + // Reject the previously hit primitive to avoid immediate self-intersection. + if (primitive_id == last_hit_primitive) { + return reject_candidate(traversal_ray); + } + // Scalar queries may provide an arbitrary primitive exclusion history. + // TODO - Think about how to provide arbitrary history checks for batch queries. for (int i = 0; i < exclude_count; ++i) { - if (exclude_primitives[i] == primitive_ref) { + if (exclude_primitives[i] == primitive_id) { return reject_candidate(traversal_ray); } } - const cuBQL::vec3i index = mesh.indices[prim_id]; + const cuBQL::vec3i vertex_indices = mesh.indices[local_index]; cuBQL::vec3d vertices[3] = { - mesh.vertices[index.x], - mesh.vertices[index.y], - mesh.vertices[index.z] + mesh.vertices[vertex_indices.x], + mesh.vertices[vertex_indices.y], + mesh.vertices[vertex_indices.z] }; cuBQL::vec3d normal = cuBQL::cross(vertices[1] - vertices[0], @@ -67,7 +69,7 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, double normal_dot_direction = dot(normal, intersection_ray.direction); - if (surface_instance.reverse_sense) { + if (surface.reverse_sense) { normal_dot_direction = -normal_dot_direction; } @@ -88,10 +90,10 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, if (intersection.hit) { hit->distance = intersection.t; hit->surface = mesh.surface_id; - hit->primitive = primitive_ref; + hit->primitive = primitive_id; hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; - hit->next_volume = surface_instance.next_volume; - hit->boundary_condition = surface_instance.boundary_condition; + hit->next_volume = surface.next_volume; + hit->boundary_condition = surface.boundary_condition; hit->normal = normal; traversal_ray.tMax = static_cast(intersection.t); } @@ -101,19 +103,16 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, return reject_candidate(traversal_ray); }; - auto leave_blas = []() -> void {}; - - cuBQL::shrinkingRayQuery::twoLevel::forEachPrim(enter_blas, - leave_blas, - intersect_prim, - volume_tlas.bvh, - traversal_ray); + // Single level traversal call for a shrinking ray query against the flattened BVH of the volume group. + cuBQL::shrinkingRayQuery::forEachPrim(intersect_prim, + volume_group.bvh, + traversal_ray); } #pragma omp end declare target void intersect_surface_tree_scalar(const cubql::Context& context, - const CuBQLVolumeTLAS& volume_tlas, + const CuBQLVolumeGroup& volume_group, const CuBQLRay& ray, CuBQLSurfaceHit& surface_hit, HitOrientation hit_orientation, @@ -148,16 +147,17 @@ intersect_surface_tree_scalar(const cubql::Context& context, gpu_id, context.hostID); - const auto volume_tlas_dd = volume_tlas.get_device_data(); + const auto volume_group_dd = volume_group.get_device_data(); const int orientation = static_cast(hit_orientation); #pragma omp target device(gpu_id) \ is_device_ptr(d_exclude_primitives, d_surface_hit) { - intersect_surface_tree(volume_tlas_dd, + intersect_surface_tree(volume_group_dd, ray, d_surface_hit, orientation, + ID_NONE, d_exclude_primitives, exclude_count); } @@ -181,21 +181,21 @@ intersect_surface_tree_scalar(const cubql::Context& context, void intersect_surface_tree_batch(const cubql::Context& context, - const CuBQLVolumeTLAS::DD* d_volume_to_tlas, + const CuBQLVolumeGroup::DD* d_volume_to_group, XDGRayHit* d_ray_hits, std::size_t num_rays, HitOrientation hit_orientation) { if (num_rays == 0) return; - if (!d_volume_to_tlas || !d_ray_hits) { + if (!d_volume_to_group || !d_ray_hits) { fatal_error("Invalid cuBQL batch intersection buffers"); } const int gpu_id = context.gpuID; #pragma omp target teams distribute parallel for device(gpu_id) \ - is_device_ptr(d_volume_to_tlas, d_ray_hits) + is_device_ptr(d_volume_to_group, d_ray_hits) for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { const XDGRayHit ray_hit = d_ray_hits[ray_id]; @@ -218,12 +218,13 @@ intersect_surface_tree_batch(const cubql::Context& context, ray.tMax = ray_hit.t_max; ray.volume = ray_hit.volume; - const CuBQLVolumeTLAS::DD volume_tlas = d_volume_to_tlas[ray.volume]; + const CuBQLVolumeGroup::DD volume_group = d_volume_to_group[ray.volume]; - intersect_surface_tree(volume_tlas, + intersect_surface_tree(volume_group, ray, &hit, static_cast(hit_orientation), + ray_hit.last_hit_primitive, nullptr, 0); } @@ -233,7 +234,8 @@ intersect_surface_tree_batch(const cubql::Context& context, d_ray_hits[ray_id].primitive = hit.primitive; d_ray_hits[ray_id].point_in_volume = static_cast(hit.piv); d_ray_hits[ray_id].next_volume = hit.next_volume; - d_ray_hits[ray_id].boundary_condition = static_cast(hit.boundary_condition); + d_ray_hits[ray_id].boundary_condition = + static_cast(hit.boundary_condition); d_ray_hits[ray_id].normal[0] = hit.normal.x; d_ray_hits[ray_id].normal[1] = hit.normal.y; d_ray_hits[ray_id].normal[2] = hit.normal.z; diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index f8f64280..0af9fedf 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -12,6 +12,7 @@ #include "cuBQL/traversal/rayQueries.h" #include +#include namespace xdg { @@ -27,42 +28,43 @@ CuBQLRayTracer::CuBQLRayTracer() CuBQLRayTracer::~CuBQLRayTracer() { - if (d_volume_to_tlas_) { - omp_target_free(d_volume_to_tlas_, context_.gpuID); - d_volume_to_tlas_ = nullptr; + if (d_volume_to_group_) { + omp_target_free(d_volume_to_group_, context_.gpuID); + d_volume_to_group_ = nullptr; } - for (auto& [tree, tlas] : tree_to_volume_tlas_) { - tlas.release(); + for (auto& [tree, group] : tree_to_volume_group_) { + group.release(); } - for (auto& [surface, blas] : surface_to_blas_map_) { - blas.release(); + for (auto& [surface, mesh] : surface_to_mesh_) { + mesh.release(); } } void CuBQLRayTracer::init() { - upload_volume_to_tlas_table_(); + upload_volume_to_group_table_(); initialized_ = true; } -void CuBQLRayTracer::upload_volume_to_tlas_table_() +void CuBQLRayTracer::upload_volume_to_group_table_() { - if (d_volume_to_tlas_) { - omp_target_free(d_volume_to_tlas_, context_.gpuID); - d_volume_to_tlas_ = nullptr; + if (d_volume_to_group_) { + omp_target_free(d_volume_to_group_, context_.gpuID); + d_volume_to_group_ = nullptr; } - if (volume_to_tlas_.empty()) { + if (volume_to_group_.empty()) { return; } - d_volume_to_tlas_ = static_cast - (omp_target_alloc(volume_to_tlas_.size() * sizeof(CuBQLVolumeTLAS::DD), context_.gpuID)); - omp_target_memcpy(d_volume_to_tlas_, - volume_to_tlas_.data(), - volume_to_tlas_.size() * sizeof(CuBQLVolumeTLAS::DD), + d_volume_to_group_ = static_cast + (omp_target_alloc(volume_to_group_.size() * sizeof(CuBQLVolumeGroup::DD), + context_.gpuID)); + omp_target_memcpy(d_volume_to_group_, + volume_to_group_.data(), + volume_to_group_.size() * sizeof(CuBQLVolumeGroup::DD), 0, 0, context_.gpuID, @@ -78,14 +80,14 @@ CuBQLRayTracer::register_volume(const std::shared_ptr& mesh_manager return {surface_tree, element_tree}; } -CuBQLSurfaceBLAS +CuBQLSurfaceMesh CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manager, - MeshID surface_id, - double bounding_box_bump) + MeshID surface_id) { - auto num_faces = mesh_manager->num_surface_faces(surface_id); + const auto num_faces = mesh_manager->num_surface_faces(surface_id); auto vertices = mesh_manager->get_surface_vertices(surface_id); auto indices = mesh_manager->get_surface_connectivity(surface_id); + auto h_primitive_ids = mesh_manager->get_surface_faces(surface_id); std::vector h_vertices; h_vertices.reserve(vertices.size()); @@ -99,8 +101,6 @@ CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manage h_indices.emplace_back(indices[i], indices[i + 1], indices[i + 2]); } - std::vector h_primitive_refs = mesh_manager->get_surface_faces(surface_id); - // TODO- think about how to better handle omp transfer calls. AutoUploadArrays is one option auto* d_vertices = static_cast (omp_target_alloc(h_vertices.size() * sizeof(cuBQL::vec3d), context_.gpuID)); @@ -122,67 +122,26 @@ CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manage context_.gpuID, context_.hostID); - auto* d_primitive_refs = static_cast - (omp_target_alloc(h_primitive_refs.size() * sizeof(MeshID), context_.gpuID)); - omp_target_memcpy(d_primitive_refs, - h_primitive_refs.data(), - h_primitive_refs.size() * sizeof(MeshID), + auto* d_primitive_ids = static_cast + (omp_target_alloc(h_primitive_ids.size() * sizeof(MeshID), context_.gpuID)); + omp_target_memcpy(d_primitive_ids, + h_primitive_ids.data(), + h_primitive_ids.size() * sizeof(MeshID), 0, 0, context_.gpuID, context_.hostID); - auto* d_aabbs = static_cast - (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3f), context_.gpuID)); - const auto num_primitives = static_cast(h_indices.size()); - - // TODO - Abstract this out into its own bounding_box creation function - #pragma omp target device(context_.gpuID) is_device_ptr(d_vertices, d_indices, d_aabbs) \ - firstprivate(bounding_box_bump) - #pragma omp teams distribute parallel for - for (uint32_t primID = 0; primID < num_primitives; ++primID) { - cuBQL::vec3i indices = d_indices[primID]; - - cuBQL::vec3d A = d_vertices[indices.x]; - cuBQL::vec3d B = d_vertices[indices.y]; - cuBQL::vec3d C = d_vertices[indices.z]; - - cuBQL::box3d aabb; - aabb.extend(A); - aabb.extend(B); - aabb.extend(C); - - const cuBQL::vec3d bump(bounding_box_bump); - aabb.lower = aabb.lower - bump; - aabb.upper = aabb.upper + bump; - - d_aabbs[primID] = cuBQL::box3f(aabb); - } - - cuBQL::BuildConfig blasBuildParams; - // TODO - Try setting leaf params to 1 to see what it does - // Check what default is for CUDA - cuBQL::bvh3f bvh; - cuBQL::build_omp_target(bvh, d_aabbs, num_faces, blasBuildParams, context_.gpuID); - - omp_target_free(d_aabbs, context_.gpuID); - CuBQLSurfaceMesh surface_mesh; surface_mesh.surface_id = surface_id; surface_mesh.d_vertices = d_vertices; surface_mesh.d_indices = d_indices; - surface_mesh.d_primitive_refs = d_primitive_refs; - surface_mesh.num_vertices = h_vertices.size(); - surface_mesh.num_triangles = num_faces; + surface_mesh.d_primitive_ids = d_primitive_ids; + surface_mesh.num_vertices = static_cast(h_vertices.size()); + surface_mesh.num_primitives = static_cast(num_faces); surface_mesh.gpu_id = context_.gpuID; - CuBQLSurfaceBLAS surface_blas; - surface_blas.bvh = bvh; - surface_blas.mesh = surface_mesh; - surface_blas.num_prims = num_faces; - surface_blas.gpu_id = context_.gpuID; - - return surface_blas; + return surface_mesh; } TreeID @@ -195,129 +154,157 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man SurfaceTreeID tree = next_surface_tree_id(); surface_trees_.push_back(tree); auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); - std::vector h_tlas_boxes; - std::vector h_surface_instances; - h_tlas_boxes.reserve(volume_surfaces.size()); - h_surface_instances.reserve(volume_surfaces.size()); - - for (const auto &surf : volume_surfaces) { - auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); - const double max_parent_bbox_bump = std::max(bounding_box_bump(mesh_manager, forward_parent), - bounding_box_bump(mesh_manager, reverse_parent)); - - if (!surface_to_blas_map_.count(surf)) { - surface_to_blas_map_[surf] = register_surface(mesh_manager, surf, max_parent_bbox_bump); - } - CuBQLSurfaceBLAS& surface_blas = surface_to_blas_map_.at(surf); + if (volume_surfaces.empty()) { + fatal_error("Volume {} has no surfaces; cannot build cuBQL surface tree", volume_id); + } - // Store BLAS bounding boxes to build TLAS - const auto surface_bounding_box = mesh_manager->surface_bounding_box(surf); - cuBQL::box3d surface_bounds_dp; - surface_bounds_dp.lower = cuBQL::vec3d(surface_bounding_box.min_x, - surface_bounding_box.min_y, - surface_bounding_box.min_z); - surface_bounds_dp.upper = cuBQL::vec3d(surface_bounding_box.max_x, - surface_bounding_box.max_y, - surface_bounding_box.max_z); + std::uint64_t num_volume_primitives = 0; + for (const MeshID surface : volume_surfaces) { + if (!surface_to_mesh_.count(surface)) { + surface_to_mesh_.emplace( + surface, register_surface(mesh_manager, surface)); + } + num_volume_primitives += surface_to_mesh_.at(surface).num_primitives; + } - const cuBQL::vec3d bump(max_parent_bbox_bump); - surface_bounds_dp.lower = surface_bounds_dp.lower - bump; - surface_bounds_dp.upper = surface_bounds_dp.upper + bump; - cuBQL::box3f surface_bounds(surface_bounds_dp); + std::vector h_surfaces; + std::vector h_prim_refs; + std::vector h_surface_bumps; + h_surfaces.reserve(volume_surfaces.size()); + h_surface_bumps.reserve(volume_surfaces.size()); + h_prim_refs.reserve(static_cast(num_volume_primitives)); - CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; - surface_instance.surface_blas = surface_blas.get_device_data(); + for (const MeshID surface : volume_surfaces) { + const auto [forward_parent, reverse_parent] = + mesh_manager->get_parent_volumes(surface); + const CuBQLSurfaceMesh& surface_mesh = surface_to_mesh_.at(surface); + + CuBQLVolumeGroup::SurfaceDD surface_data; + surface_data.mesh = surface_mesh.get_device_data(); - // Store per-instance topology for this volume and surface. if (volume_id == forward_parent) { - surface_instance.reverse_sense = false; - surface_instance.next_volume = reverse_parent; + surface_data.next_volume = reverse_parent; } else if (volume_id == reverse_parent) { - surface_instance.reverse_sense = true; - surface_instance.next_volume = forward_parent; + surface_data.reverse_sense = true; + surface_data.next_volume = forward_parent; } else { - fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); + fatal_error("Volume {} is not a parent of surface {}", volume_id, surface); } - // Store boundary-condition metadata for this surface instance. const auto property = mesh_manager->get_surface_property( - surf, PropertyType::BOUNDARY_CONDITION); - + surface, PropertyType::BOUNDARY_CONDITION); if (property.value == "vacuum") { - surface_instance.boundary_condition = VACUUM; + surface_data.boundary_condition = VACUUM; } else if (property.value == "reflecting" || property.value == "reflective") { - surface_instance.boundary_condition = REFLECTIVE; + surface_data.boundary_condition = REFLECTIVE; } else if (property.value == "transmission") { - surface_instance.boundary_condition = TRANSMISSION; + surface_data.boundary_condition = TRANSMISSION; } else { fatal_error("Unsupported boundary condition '{}' on surface {}", - property.value, surf); + property.value, surface); } - h_tlas_boxes.push_back(surface_bounds); - h_surface_instances.push_back(surface_instance); - } + const auto surface_index = static_cast(h_surfaces.size()); + h_surfaces.push_back(surface_data); + h_surface_bumps.push_back(std::max( + bounding_box_bump(mesh_manager, forward_parent), + bounding_box_bump(mesh_manager, reverse_parent))); - if (h_surface_instances.empty()) { - fatal_error("Volume {} has no surfaces; cannot build cuBQL surface tree", volume_id); + for (std::uint32_t primitive_index = 0; + primitive_index < surface_mesh.num_primitives; + ++primitive_index) { + h_prim_refs.push_back({surface_index, primitive_index}); + } } - auto* d_tlas_boxes = static_cast - (omp_target_alloc(h_tlas_boxes.size() * sizeof(cuBQL::box3f), context_.gpuID)); - omp_target_memcpy(d_tlas_boxes, - h_tlas_boxes.data(), - h_tlas_boxes.size() * sizeof(cuBQL::box3f), + auto* d_aabbs = static_cast + (omp_target_alloc(h_prim_refs.size() * sizeof(cuBQL::box3f), + context_.gpuID)); + auto* d_surfaces = static_cast + (omp_target_alloc(h_surfaces.size() * sizeof(CuBQLVolumeGroup::SurfaceDD), + context_.gpuID)); + omp_target_memcpy(d_surfaces, + h_surfaces.data(), + h_surfaces.size() * sizeof(CuBQLVolumeGroup::SurfaceDD), + 0, + 0, + context_.gpuID, + context_.hostID); + + auto* d_prim_refs = static_cast + (omp_target_alloc(h_prim_refs.size() * sizeof(CuBQLVolumeGroup::PrimRef), + context_.gpuID)); + omp_target_memcpy(d_prim_refs, + h_prim_refs.data(), + h_prim_refs.size() * sizeof(CuBQLVolumeGroup::PrimRef), 0, 0, context_.gpuID, context_.hostID); - auto* d_surface_instances = static_cast - (omp_target_alloc(h_surface_instances.size() * sizeof(CuBQLVolumeTLAS::SurfaceInstanceDD), context_.gpuID)); - omp_target_memcpy(d_surface_instances, - h_surface_instances.data(), - h_surface_instances.size() * sizeof(CuBQLVolumeTLAS::SurfaceInstanceDD), + auto* d_surface_bumps = static_cast + (omp_target_alloc(h_surface_bumps.size() * sizeof(double), context_.gpuID)); + omp_target_memcpy(d_surface_bumps, + h_surface_bumps.data(), + h_surface_bumps.size() * sizeof(double), 0, 0, context_.gpuID, context_.hostID); - cuBQL::BuildConfig tlasBuildParams; - tlasBuildParams.makeLeafThreshold = 1; - tlasBuildParams.maxAllowedLeafSize = 1; - - CuBQLVolumeTLAS volume_tlas; - volume_tlas.volume_id = volume_id; // store meshid in the TLAS object for easier mapping between the two - volume_tlas.num_surface_instances = static_cast(h_surface_instances.size()); - volume_tlas.gpu_id = context_.gpuID; - volume_tlas.d_surface_instances = d_surface_instances; - cuBQL::build_omp_target(volume_tlas.bvh, - d_tlas_boxes, - volume_tlas.num_surface_instances, - tlasBuildParams, + const auto num_primitives = static_cast(h_prim_refs.size()); + const int gpu_id = context_.gpuID; + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(d_aabbs, d_surfaces, d_prim_refs, d_surface_bumps) + for (std::uint32_t primID = 0; primID < num_primitives; ++primID) { + const auto primitive = d_prim_refs[primID]; + const auto mesh = d_surfaces[primitive.surface_index].mesh; + const auto indices = mesh.indices[primitive.primitive_index]; + + cuBQL::box3d aabb; + aabb.extend(mesh.vertices[indices.x]); + aabb.extend(mesh.vertices[indices.y]); + aabb.extend(mesh.vertices[indices.z]); + + // get correct bump for surface + const cuBQL::vec3d bump(d_surface_bumps[primitive.surface_index]); + aabb.lower = aabb.lower - bump; + aabb.upper = aabb.upper + bump; + d_aabbs[primID] = cuBQL::box3f(aabb); + } + + CuBQLVolumeGroup volume_group; + volume_group.d_surfaces = d_surfaces; + volume_group.d_prim_refs = d_prim_refs; + volume_group.num_surfaces = static_cast(h_surfaces.size()); + volume_group.num_primitives = num_primitives; + volume_group.gpu_id = context_.gpuID; + + cuBQL::BuildConfig build_params; + cuBQL::build_omp_target(volume_group.bvh, + d_aabbs, + num_primitives, + build_params, context_.gpuID); - omp_target_free(d_tlas_boxes, context_.gpuID); + omp_target_free(d_surface_bumps, context_.gpuID); + omp_target_free(d_aabbs, context_.gpuID); - // Still required for lifetime and scalar calls which need to resolve TreeID->volume_tlas on CPU side. - auto result = tree_to_volume_tlas_.emplace(tree, std::move(volume_tlas)); + // Retain owning objects for scalar TreeID lookups and allocation lifetime. + auto result = tree_to_volume_group_.emplace(tree, std::move(volume_group)); auto it = result.first; - // Keep a dense host-side MeshID -> TLAS device-data table for prepared queries. - // The TLAS object in tree_to_volume_tlas_ owns the device allocations; this table - // only stores lightweight DD views indexed by volume ID. Upload to device once in - // init(), unless a volume is registered after initialization. - + // Keep a dense host-side MeshID -> group device-data table for batch queries. const auto volume_index = static_cast(volume_id); - if (volume_index >= volume_to_tlas_.size()) { - volume_to_tlas_.resize(volume_index + 1); + if (volume_index >= volume_to_group_.size()) { + volume_to_group_.resize(volume_index + 1); } - volume_to_tlas_[volume_index] = it->second.get_device_data(); + volume_to_group_[volume_index] = it->second.get_device_data(); if (initialized_) { - upload_volume_to_tlas_table_(); + upload_volume_to_group_table_(); } return tree; @@ -359,7 +346,7 @@ bool CuBQLRayTracer::point_in_volume(TreeID tree, const std::vector* exclude_primitives) const { const auto& context = context_; - const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); + const CuBQLVolumeGroup& volume_group = tree_to_volume_group_.at(tree); // Use provided direction or if Direction == nulptr use default direction Direction directionUsed = (direction != nullptr) ? Direction{direction->x, direction->y, direction->z} @@ -375,7 +362,12 @@ bool CuBQLRayTracer::point_in_volume(TreeID tree, CuBQLSurfaceHit surface_hit; // TODO - Maybe we can come up with a better name for this - intersect_surface_tree_scalar(context, volume_tlas, ray, surface_hit, HitOrientation::ANY, exclude_primitives); + intersect_surface_tree_scalar(context, + volume_group, + ray, + surface_hit, + HitOrientation::ANY, + exclude_primitives); // if the ray hit nothing the point must be outside the volume if (surface_hit.primitive == ID_NONE) return false; @@ -392,7 +384,7 @@ CuBQLRayTracer::ray_fire(TreeID tree, std::vector* const exclude_primitives) { const auto& context = context_; - const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); + const CuBQLVolumeGroup& volume_group = tree_to_volume_group_.at(tree); CuBQLRay ray; ray.origin = cuBQL::vec3d(origin.x, origin.y, origin.z); @@ -403,7 +395,12 @@ CuBQLRayTracer::ray_fire(TreeID tree, CuBQLSurfaceHit surface_hit; // TODO - Maybe we can come up with a better name for this - intersect_surface_tree_scalar(context, volume_tlas, ray, surface_hit, hitOrientation, exclude_primitives); + intersect_surface_tree_scalar(context, + volume_group, + ray, + surface_hit, + hitOrientation, + exclude_primitives); if (surface_hit.primitive == ID_NONE) { return {INFTY, ID_NONE}; @@ -454,12 +451,12 @@ CuBQLRayTracer::ray_fire_batch(const XDGRayHitBuffer& ray_hits, fatal_error("Invalid cuBQL XDG ray-hit buffer"); } - if (!d_volume_to_tlas_) { - fatal_error("cuBQL volume TLAS lookup table has not been uploaded"); + if (!d_volume_to_group_) { + fatal_error("cuBQL volume-group lookup table has not been uploaded"); } intersect_surface_tree_batch(context_, - d_volume_to_tlas_, + d_volume_to_group_, ray_hits.data, ray_hits.count, hit_orientation); diff --git a/src/cuBQL/triangles.cpp b/src/cuBQL/triangles.cpp index 73cf1fe3..9a3cf5ac 100644 --- a/src/cuBQL/triangles.cpp +++ b/src/cuBQL/triangles.cpp @@ -14,26 +14,16 @@ void CuBQLSurfaceMesh::release() omp_target_free(d_indices, gpu_id); d_indices = nullptr; } - if (d_primitive_refs) { - omp_target_free(d_primitive_refs, gpu_id); - d_primitive_refs = nullptr; + if (d_primitive_ids) { + omp_target_free(d_primitive_ids, gpu_id); + d_primitive_ids = nullptr; } -} -void CuBQLSurfaceBLAS::release() -{ - if (bvh.primIDs) { - omp_target_free(bvh.primIDs, gpu_id); - bvh.primIDs = nullptr; - } - if (bvh.nodes) { - omp_target_free(bvh.nodes, gpu_id); - bvh.nodes = nullptr; - } - mesh.release(); + num_vertices = 0; + num_primitives = 0; } -void CuBQLVolumeTLAS::release() +void CuBQLVolumeGroup::release() { if (bvh.primIDs) { omp_target_free(bvh.primIDs, gpu_id); @@ -43,10 +33,19 @@ void CuBQLVolumeTLAS::release() omp_target_free(bvh.nodes, gpu_id); bvh.nodes = nullptr; } - if (d_surface_instances) { - omp_target_free(d_surface_instances, gpu_id); - d_surface_instances = nullptr; + if (d_surfaces) { + omp_target_free(d_surfaces, gpu_id); + d_surfaces = nullptr; + } + if (d_prim_refs) { + omp_target_free(d_prim_refs, gpu_id); + d_prim_refs = nullptr; } + + bvh.numNodes = 0; + bvh.numPrims = 0; + num_surfaces = 0; + num_primitives = 0; } } // namespace xdg diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index 002cdbf7..dab27410 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -249,6 +249,7 @@ int main(int argc, char** argv) ray_hit.t_min = 0.0; ray_hit.t_max = INFTY; ray_hit.volume = volume; + ray_hit.last_hit_primitive = ID_NONE; ray_hit.distance = INFTY; ray_hit.surface = ID_NONE; ray_hit.primitive = ID_NONE; From 771984dbfc48a0a54d3559cfea72aea6ca3201f7 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 20 Jul 2026 14:39:49 +0100 Subject: [PATCH 15/23] Track maximum parent-volume bump on shared cuBQL surface meshes --- include/xdg/cuBQL/triangles.h | 3 ++ src/cuBQL/intersection.cpp | 19 ++++-------- src/cuBQL/ray_tracer.cpp | 57 ++++++++++++----------------------- 3 files changed, 28 insertions(+), 51 deletions(-) diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h index 0f190d4d..07016816 100644 --- a/include/xdg/cuBQL/triangles.h +++ b/include/xdg/cuBQL/triangles.h @@ -24,6 +24,7 @@ struct CuBQLSurfaceMesh { struct DD { // Topological metadata MeshID surface_id {ID_NONE}; + double max_parent_volume_bump {0.0}; // Geometric data const cuBQL::vec3d* vertices {nullptr}; @@ -33,6 +34,7 @@ struct CuBQLSurfaceMesh { // Topological metadata MeshID surface_id {ID_NONE}; + double max_parent_volume_bump {0.0}; // Device buffers for primitive data cuBQL::vec3d* d_vertices {nullptr}; @@ -48,6 +50,7 @@ struct CuBQLSurfaceMesh { { return { surface_id, + max_parent_volume_bump, d_vertices, d_indices, d_primitive_ids diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index 49a3a151..e469803f 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -73,8 +73,7 @@ static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, normal_dot_direction = -normal_dot_direction; } - if (orientation_cull(normal_dot_direction, - static_cast(orientation))) { + if (orientation_cull(normal_dot_direction, static_cast(orientation))) { return reject_candidate(traversal_ray); } @@ -104,9 +103,7 @@ static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, }; // Single level traversal call for a shrinking ray query against the flattened BVH of the volume group. - cuBQL::shrinkingRayQuery::forEachPrim(intersect_prim, - volume_group.bvh, - traversal_ray); + cuBQL::shrinkingRayQuery::forEachPrim(intersect_prim, volume_group.bvh, traversal_ray); } #pragma omp end declare target @@ -126,6 +123,7 @@ intersect_surface_tree_scalar(const cubql::Context& context, exclude_count = static_cast(exclude_primitives->size()); d_exclude_primitives = static_cast (omp_target_alloc(exclude_count * sizeof(MeshID), gpu_id)); + omp_target_memcpy(d_exclude_primitives, exclude_primitives->data(), exclude_count * sizeof(MeshID), @@ -208,12 +206,8 @@ intersect_surface_tree_batch(const cubql::Context& context, if (ray_hit.volume != ID_NONE) { CuBQLRay ray; - ray.origin = cuBQL::vec3d(ray_hit.origin[0], - ray_hit.origin[1], - ray_hit.origin[2]); - ray.direction = cuBQL::vec3d(ray_hit.direction[0], - ray_hit.direction[1], - ray_hit.direction[2]); + ray.origin = cuBQL::vec3d(ray_hit.origin[0], ray_hit.origin[1], ray_hit.origin[2]); + ray.direction = cuBQL::vec3d(ray_hit.direction[0], ray_hit.direction[1], ray_hit.direction[2]); ray.tMin = ray_hit.t_min; ray.tMax = ray_hit.t_max; ray.volume = ray_hit.volume; @@ -234,8 +228,7 @@ intersect_surface_tree_batch(const cubql::Context& context, d_ray_hits[ray_id].primitive = hit.primitive; d_ray_hits[ray_id].point_in_volume = static_cast(hit.piv); d_ray_hits[ray_id].next_volume = hit.next_volume; - d_ray_hits[ray_id].boundary_condition = - static_cast(hit.boundary_condition); + d_ray_hits[ray_id].boundary_condition = static_cast(hit.boundary_condition); d_ray_hits[ray_id].normal[0] = hit.normal.x; d_ray_hits[ray_id].normal[1] = hit.normal.y; d_ray_hits[ray_id].normal[2] = hit.normal.z; diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 0af9fedf..fc25677a 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -159,25 +159,26 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man fatal_error("Volume {} has no surfaces; cannot build cuBQL surface tree", volume_id); } + const double volume_bump = bounding_box_bump(mesh_manager, volume_id); + std::uint64_t num_volume_primitives = 0; for (const MeshID surface : volume_surfaces) { if (!surface_to_mesh_.count(surface)) { - surface_to_mesh_.emplace( - surface, register_surface(mesh_manager, surface)); + surface_to_mesh_.emplace(surface, register_surface(mesh_manager, surface)); } - num_volume_primitives += surface_to_mesh_.at(surface).num_primitives; + + auto& surface_mesh = surface_to_mesh_.at(surface); + surface_mesh.max_parent_volume_bump = std::max(surface_mesh.max_parent_volume_bump, volume_bump); + num_volume_primitives += surface_mesh.num_primitives; } std::vector h_surfaces; std::vector h_prim_refs; - std::vector h_surface_bumps; h_surfaces.reserve(volume_surfaces.size()); - h_surface_bumps.reserve(volume_surfaces.size()); h_prim_refs.reserve(static_cast(num_volume_primitives)); for (const MeshID surface : volume_surfaces) { - const auto [forward_parent, reverse_parent] = - mesh_manager->get_parent_volumes(surface); + const auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surface); const CuBQLSurfaceMesh& surface_mesh = surface_to_mesh_.at(surface); CuBQLVolumeGroup::SurfaceDD surface_data; @@ -192,30 +193,22 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man fatal_error("Volume {} is not a parent of surface {}", volume_id, surface); } - const auto property = mesh_manager->get_surface_property( - surface, PropertyType::BOUNDARY_CONDITION); + const auto property = mesh_manager->get_surface_property(surface, PropertyType::BOUNDARY_CONDITION); if (property.value == "vacuum") { surface_data.boundary_condition = VACUUM; - } else if (property.value == "reflecting" || - property.value == "reflective") { + } else if (property.value == "reflecting" || property.value == "reflective") { surface_data.boundary_condition = REFLECTIVE; } else if (property.value == "transmission") { surface_data.boundary_condition = TRANSMISSION; } else { - fatal_error("Unsupported boundary condition '{}' on surface {}", - property.value, surface); + fatal_error("Unsupported boundary condition '{}' on surface {}", property.value, surface); } const auto surface_index = static_cast(h_surfaces.size()); h_surfaces.push_back(surface_data); - h_surface_bumps.push_back(std::max( - bounding_box_bump(mesh_manager, forward_parent), - bounding_box_bump(mesh_manager, reverse_parent))); - - for (std::uint32_t primitive_index = 0; - primitive_index < surface_mesh.num_primitives; - ++primitive_index) { - h_prim_refs.push_back({surface_index, primitive_index}); + + for (std::uint32_t prim_index = 0; prim_index < surface_mesh.num_primitives; ++prim_index) { + h_prim_refs.push_back({surface_index, prim_index}); } } @@ -244,23 +237,14 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man context_.gpuID, context_.hostID); - auto* d_surface_bumps = static_cast - (omp_target_alloc(h_surface_bumps.size() * sizeof(double), context_.gpuID)); - omp_target_memcpy(d_surface_bumps, - h_surface_bumps.data(), - h_surface_bumps.size() * sizeof(double), - 0, - 0, - context_.gpuID, - context_.hostID); - const auto num_primitives = static_cast(h_prim_refs.size()); const int gpu_id = context_.gpuID; #pragma omp target teams distribute parallel for device(gpu_id) \ - is_device_ptr(d_aabbs, d_surfaces, d_prim_refs, d_surface_bumps) + is_device_ptr(d_aabbs, d_surfaces, d_prim_refs) for (std::uint32_t primID = 0; primID < num_primitives; ++primID) { const auto primitive = d_prim_refs[primID]; - const auto mesh = d_surfaces[primitive.surface_index].mesh; + const auto surface = d_surfaces[primitive.surface_index]; + const auto mesh = surface.mesh; const auto indices = mesh.indices[primitive.primitive_index]; cuBQL::box3d aabb; @@ -268,8 +252,7 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man aabb.extend(mesh.vertices[indices.y]); aabb.extend(mesh.vertices[indices.z]); - // get correct bump for surface - const cuBQL::vec3d bump(d_surface_bumps[primitive.surface_index]); + const cuBQL::vec3d bump(mesh.max_parent_volume_bump); aabb.lower = aabb.lower - bump; aabb.upper = aabb.upper + bump; d_aabbs[primID] = cuBQL::box3f(aabb); @@ -289,7 +272,6 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man build_params, context_.gpuID); - omp_target_free(d_surface_bumps, context_.gpuID); omp_target_free(d_aabbs, context_.gpuID); // Retain owning objects for scalar TreeID lookups and allocation lifetime. @@ -351,7 +333,6 @@ bool CuBQLRayTracer::point_in_volume(TreeID tree, // Use provided direction or if Direction == nulptr use default direction Direction directionUsed = (direction != nullptr) ? Direction{direction->x, direction->y, direction->z} : Direction{1. / std::sqrt(2.0), 1. / std::sqrt(2.0), 0.0}; - CuBQLRay ray; ray.origin = cuBQL::vec3d(point.x, point.y, point.z); @@ -360,7 +341,7 @@ bool CuBQLRayTracer::point_in_volume(TreeID tree, ray.tMax = INFTY; CuBQLSurfaceHit surface_hit; - + // TODO - Maybe we can come up with a better name for this intersect_surface_tree_scalar(context, volume_group, From 14652faf4afa625896f4c7e64387f0f2d0f01fc1 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 13 Jul 2026 11:34:59 +0100 Subject: [PATCH 16/23] Made a start on adding an event based particle sim - Basic structure of algorithm in place - Event based simulation data structure - New EventQueue struct for handling device queues --- tools/particle_sim_event.cpp | 121 ++++++++++++ tools/particle_sim_event.h | 347 +++++++++++++++++++++++++++++++++++ 2 files changed, 468 insertions(+) create mode 100644 tools/particle_sim_event.cpp create mode 100644 tools/particle_sim_event.h diff --git a/tools/particle_sim_event.cpp b/tools/particle_sim_event.cpp new file mode 100644 index 00000000..47d1328f --- /dev/null +++ b/tools/particle_sim_event.cpp @@ -0,0 +1,121 @@ +#include +#include +#include + +#include "xdg/error.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/vec3da.h" +#include "xdg/xdg.h" + +#include "argparse/argparse.hpp" + +#include "particle_sim.h" + +using namespace xdg; + +int main(int argc, char** argv) { + +// argument parsing +argparse::ArgumentParser args("XDG Particle Pseudo-Simulation", "1.0", argparse::default_arguments::help); + +args.add_argument("filename") + .help("Path to the input file"); + +args.add_argument("-v", "--verbose") + .default_value(false) + .implicit_value(true) + .help("Enable verbose output of particle events"); + +args.add_argument("-m", "--mfp") + .default_value(1.0) + .help("Mean free path of the particles").scan<'g', double>(); + +args.add_argument("-n", "--n-particles") + .default_value(100u) + .help("Number of particles to simulate").scan<'u', uint32_t>(); + +args.add_argument("-e", "--max-events") + .default_value(1000u) + .help("Maximum number of events per particle").scan<'u', uint32_t>(); + +args.add_argument("-g", "--ipc-graveyard") + .default_value(false) + .implicit_value(true) + .help("Treat the implicit complement as a graveyard (i.e. particles that enter it are killed)"); + +args.add_argument("-m", "--mesh-library") + .help("Mesh library to use. One of (MOAB, LIBMESH)") + .default_value("MOAB"); + +args.add_argument("-r", "--rt-library") + .help("Ray tracing library to use. One of (EMBREE, GPRT)") + .default_value("EMBREE"); +try { + args.parse_args(argc, argv); +} +catch (const std::runtime_error& err) { + std::cout << err.what() << std::endl; + std::cout << args; + exit(0); +} + +// Problem Setup +srand48(42); + +SimulationData sim_data; + +// create a mesh manager +std::string mesh_str = args.get("--mesh-library"); +std::string rt_str = args.get("--rt-library"); + +RTLibrary rt_lib; +if (rt_str == "EMBREE") + rt_lib = RTLibrary::EMBREE; +else if (rt_str == "GPRT") + rt_lib = RTLibrary::GPRT; +else + fatal_error("Invalid ray tracing library '{}' specified", rt_str); + +MeshLibrary mesh_lib; +if (mesh_str == "MOAB") + mesh_lib = MeshLibrary::MOAB; +else if (mesh_str == "LIBMESH") { + mesh_lib = MeshLibrary::LIBMESH; + if (rt_lib == RTLibrary::GPRT) + fatal_error("LibMesh is not currently supported with GPRT"); +} +else + fatal_error("Invalid mesh library '{}' specified", mesh_str); + +// create an XDG instance with the specified mesh and ray tracing library +std::shared_ptr xdg = XDG::create(mesh_lib, rt_lib); +const auto& mm = xdg->mesh_manager(); +mm->load_file(args.get("filename")); +mm->init(); +mm->parse_metadata(); +xdg->prepare_raytracer(); + +sim_data.xdg_ = xdg; + +// update the mean free path +sim_data.mfp_ = args.get("--mfp"); + +sim_data.verbose_particles_ = args.get("--verbose"); +sim_data.implicit_complement_is_graveyard_ = args.get("--ipc-graveyard"); +sim_data.n_particles_ = args.get("--n-particles"); +sim_data.max_events_ = args.get("--max-events"); + +transport_particles(sim_data); +// transport_particles_event_based(sim_data); + +// report distances in each cell in a table +write_message("Cell Track Lengths"); +write_message("-----------"); +for (const auto& [cell, dist] : sim_data.cell_tracks) { + write_message("Cell {}: {}", cell, dist); +} +write_message("-----------"); + + +return 0; +} diff --git a/tools/particle_sim_event.h b/tools/particle_sim_event.h new file mode 100644 index 00000000..11fab9d0 --- /dev/null +++ b/tools/particle_sim_event.h @@ -0,0 +1,347 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "xdg/error.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/vec3da.h" +#include "xdg/xdg.h" + +using namespace xdg; + +// Lightweight append queue used by event queues. +// Making use of the same internal device data view pattern I borrowed from cuBQL/DPRT +template +struct DeviceAppendQueue { + struct DD { + T* data {nullptr}; + int* size {nullptr}; // Number of current valid entries + int capacity {0}; // Current allocated limit on number of entries before container overflow + + int thread_safe_append(const T& value) + { + int idx; + #pragma omp atomic capture + idx = (*size)++; + + if (idx >= capacity) { + #pragma omp atomic write + *size = capacity; + return -1; + } + + data[idx] = value; + return idx; + }; + }; + + T* d_data {nullptr}; + int* d_size {nullptr}; + int h_size {0}; + int capacity {0}; + int gpu_id {0}; + int host_id {omp_get_initial_device()}; + + void allocate(int capacity_, int gpu_id_) + { + release(); + + capacity = capacity_; + gpu_id = gpu_id_; + host_id = omp_get_initial_device(); + h_size = 0; + + if (capacity == 0) return; + + d_data = static_cast(omp_target_alloc(capacity * sizeof(T), gpu_id)); + d_size = static_cast(omp_target_alloc(sizeof(int), gpu_id)); + + if (!d_data || !d_size) { + release(); + fatal_error("Failed to allocate event queue on OpenMP target device."); + } + + resize(0); + } + + void release() + { + if (d_data) { + omp_target_free(d_data, gpu_id); + d_data = nullptr; + } + if (d_size) { + omp_target_free(d_size, gpu_id); + d_size = nullptr; + } + + h_size = 0; + capacity = 0; + } + + void resize(int size) + { + h_size = size; + if (!d_size) return; + + omp_target_memcpy(d_size, + &h_size, + sizeof(int), + 0, + 0, + gpu_id, + host_id); + } + + // Small wrapper to call resize with size 0 + void reset() + { + resize(0); + } + + void sync_size_device_to_host() + { + if (!d_size) return; + + omp_target_memcpy(&h_size, + d_size, + sizeof(int), + 0, + 0, + host_id, + gpu_id); + } + + int size() const { return h_size; } + + DD get_device_data() const + { + return {d_data, d_size, capacity}; + } +}; +struct EventQueueItem { + uint32_t idx; // particle index in event-based particle buffer +}; + +using ParticleEventQueue = DeviceAppendQueue; + +struct Particle { + +Particle(std::shared_ptr xdg, uint32_t id, uint32_t max_events, bool verbose=true, bool ipc_graveyard=false) +: verbose_(verbose), xdg_(xdg), id_(id), max_events_(max_events), ipc_graveyard_(ipc_graveyard) {} + +template +void log (const std::string& msg, const Params&... fmt_args) { + if (!verbose_) return; + write_message(msg, fmt_args...); +} + +void initialize() { + // TODO: replace with sampling + r_ = {0.0, 0.0, 0.0}; + u_ = {1.0, 0.0, 0.0}; + + volume_ = xdg_->find_volume(r_, u_); + log("Particle {} initialized in volume {}", id_, volume_); +} + +void surf_dist() { + surface_intersection_ = xdg_->ray_fire(volume_, r_, u_, INFTY, HitOrientation::EXITING, &history_); + if (surface_intersection_.first == 0.0) { + fatal_error("Particle {} stuck at position ({}, {}, {}) on surfacce {}", id_, r_.x, r_.y, r_.z, surface_intersection_.second); + alive_ = false; + return; + } + if (surface_intersection_.second == ID_NONE) { + fatal_error("Particle {} lost in volume {}", id_, volume_); + alive_ = false; + return; + } + log("Intersected surface {} at distance {} ", surface_intersection_.second, surface_intersection_.first); +} + +void sample_collision_distance(double mfp) { + collision_distance_ = -std::log(1.0 - drand48()) * mfp; +} + +void collide() { + n_events_++; + log("Event {} for particle {}", n_events_, id_); + u_ = rand_dir(); + log("Particle {} collides with material at position ({}, {}, {}), new direction is ({}, {}, {})", id_, r_.x, r_.y, r_.z, u_.z, u_.y, u_.z); + history_.clear(); +} + +void advance(std::unordered_map& cell_tracks) +{ + log("Comparing surface intersection distance {} to collision distance {}", surface_intersection_.first, collision_distance_); + if (collision_distance_ < surface_intersection_.first) { + r_ += collision_distance_ * u_; + cell_tracks[volume_] += collision_distance_; + log("Particle {} collides with material at position ({}, {}, {}) ", id_, r_.x, r_.y, r_.z); + } else { + r_ += surface_intersection_.first * u_; + cell_tracks[volume_] += surface_intersection_.first; + log("Particle {} advances to surface {} at position ({}, {}, {}) ", id_, surface_intersection_.second, r_.x, r_.y, r_.z); + } +} + +void cross_surface() +{ + n_events_++; + log("Event {} for particle {}", n_events_, id_); + auto boundary_condition = xdg_->mesh_manager()->get_surface_property(surface_intersection_.second, PropertyType::BOUNDARY_CONDITION); + // check for the surface boundary condition + if (boundary_condition.value == "reflecting" || boundary_condition.value == "reflective") { + log("Particle {} reflects off surface {}", id_, surface_intersection_.second); + log("Direction before reflection: ({}, {}, {})", u_.x, u_.y, u_.z); + + Direction normal = xdg_->surface_normal(surface_intersection_.second, r_, &history_); + log("Normal to surface: ({}, {}, {})", normal.x, normal.y, normal.z); + + double proj = dot(normal, u_); + double mag = normal.length(); + normal = normal * (2.0 * proj/mag); + u_ = u_ - normal; + u_ = u_.normalize(); + log("Direction after reflection: ({}, {}, {})", u_.x, u_.y, u_.z); + // reset to last intersection + if (history_.size() > 0) { + log("Resetting particle history to last intersection"); + history_ = {history_.back()}; + } + } else if (boundary_condition.value == "vacuum") { + log("Particle {} encounters vacuum boundary at surface {}", id_, surface_intersection_.second); + alive_ = false; + } else { + volume_ = xdg_->mesh_manager()->next_volume(volume_, surface_intersection_.second); + log("Particle {} enters volume {}", id_, volume_); + if (ipc_graveyard_ && volume_ == xdg_->mesh_manager()->implicit_complement()) volume_ = ID_NONE; + if (volume_ == ID_NONE) { + alive_ = false; + return; + } + } +} + +// Data Members +bool verbose_ {true}; +std::shared_ptr xdg_; +uint32_t id_ {0}; +int32_t max_events_ {1000}; +bool ipc_graveyard_ {false}; + +Position r_; +Direction u_; +MeshID volume_ {ID_NONE}; +std::vector history_ {}; +std::pair surface_intersection_ {INFTY, ID_NONE}; +double collision_distance_ {INFTY}; +int32_t n_events_ {0}; +bool alive_ {true}; +}; + +struct EventSimulationData { + std::shared_ptr xdg_; + double mfp_ {1.0}; + uint32_t n_particles_ {100}; + uint32_t max_events_ {1000}; + bool verbose_particles_ {false}; + bool implicit_complement_is_graveyard_ {false}; + std::unordered_map cell_tracks; + + uint32_t max_particles_in_flight_ {100000}; + std::vector particles; + ParticleEventQueue advance_particle_queue; + ParticleEventQueue surface_crossing_queue; + ParticleEventQueue collision_queue; +}; + +void process_init_events(EventSimulationData& sim_data); +void process_advance_particle_events(EventSimulationData& sim_data); +void process_surface_crossing_events(EventSimulationData& sim_data); +void process_collision_events(EventSimulationData& sim_data); + +inline void transport_particle_event_based(EventSimulationData& sim_data) { + // MPI will be needed for multi GPU + // #ifdef OPENMC_MPI + // MPI_Barrier( mpi::intracomm ); + // #endif + + /* + A couple of different concepts from the event based algorithm don't need to be considered here: + - We don't need any fuel/xs lookup related events. + - Death events amount to essentially just setting the flag alive_ = false + - We don't need to worry about secondary particles/a revival bank. This is a consequence of nuclear physics + - We don't worry about in flight particles and batch all particles in one go + + */ + + sim_data.particles.clear(); + sim_data.particles.reserve(sim_data.n_particles_); + + const int gpu_id = omp_get_default_device(); + sim_data.advance_particle_queue.allocate(sim_data.n_particles_, gpu_id); + sim_data.surface_crossing_queue.allocate(sim_data.n_particles_, gpu_id); + sim_data.collision_queue.allocate(sim_data.n_particles_, gpu_id); + + process_init_events(sim_data); + + // Event-based transport loop + while (true) { + // Determine which event kernel has the longest queue + int64_t max = std::max({ + sim_data.advance_particle_queue.size(), + sim_data.surface_crossing_queue.size(), + sim_data.collision_queue.size()}); + + + // Execute event with the longest queue + if (max == 0) { + break; + } else if (max == sim_data.advance_particle_queue.size()) { + process_advance_particle_events(sim_data); + } else if (max == sim_data.surface_crossing_queue.size()) { + process_surface_crossing_events(sim_data); + } else if (max == sim_data.collision_queue.size()) { + process_collision_events(sim_data); + } + } + + + // MPI will be needed for multi gpu + // #ifdef OPENMC_MPI + // MPI_Barrier( mpi::intracomm ); + // #endif +} + + +/* + initialize position/direction + find starting volume (call this globally for all particles [in flight]) + set global particle id + reset event counter + reset alive flag + reset ray history + reset pending surface hit / collision distance + initialize per-particle RNG seed if you stop using global drand48() + enqueue into advance queue + + Essentially set up particle state (for all particles) ready to call advance_particle +*/ + +// void process_death_events(); +/* + Particle death in the pseudo transport app doesn't really mean all that much. + We are essentially just setting the particle.alive_ member to false. + I think a better approach is to actually just handle death as it happens. + Rather than waiting for everything a particle's death state should update when it + occurs. So we don't both with a process_death_events() method. +*/ From 6c0cccd136534747118d3aba7293e24171b71e52 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 13 Jul 2026 12:08:16 +0100 Subject: [PATCH 17/23] Added a new EventParticle class to handle event based particle state - Currently does not share any code with history based Particle class - Updated EventSimulationData to store a device friendly array of particles - Made a new tools::random namespace to share common LCG PRNG logic between XDG miniapps --- tools/particle_sim_event.h | 219 ++++++++++++++++++------------------- tools/random_lcg.h | 43 ++++++++ tools/ray_benchmark.h | 30 +---- 3 files changed, 155 insertions(+), 137 deletions(-) create mode 100644 tools/random_lcg.h diff --git a/tools/particle_sim_event.h b/tools/particle_sim_event.h index 11fab9d0..31347cf5 100644 --- a/tools/particle_sim_event.h +++ b/tools/particle_sim_event.h @@ -14,6 +14,8 @@ #include "xdg/vec3da.h" #include "xdg/xdg.h" +#include "random_lcg.h" + using namespace xdg; // Lightweight append queue used by event queues. @@ -132,133 +134,115 @@ struct EventQueueItem { using ParticleEventQueue = DeviceAppendQueue; -struct Particle { - -Particle(std::shared_ptr xdg, uint32_t id, uint32_t max_events, bool verbose=true, bool ipc_graveyard=false) -: verbose_(verbose), xdg_(xdg), id_(id), max_events_(max_events), ipc_graveyard_(ipc_graveyard) {} +#ifdef _OPENMP +#pragma omp declare target +#endif -template -void log (const std::string& msg, const Params&... fmt_args) { - if (!verbose_) return; - write_message(msg, fmt_args...); -} +struct EventParticle { + void initialize(uint32_t id, + std::uint32_t seed, + Position r, + Direction u, + MeshID volume) + { + id_ = id; + r_.x = r.x; + r_.y = r.y; + r_.z = r.z; + u_.x = u.x; + u_.y = u.y; + u_.z = u.z; + volume_ = volume; + surface_hit_ = ID_NONE; + surface_hit_distance_ = INFTY; + collision_distance_ = INFTY; + last_surface_hit_ = ID_NONE; + rng_state_ = seed ^ id; + n_events_ = 0; + alive_ = true; + } -void initialize() { - // TODO: replace with sampling - r_ = {0.0, 0.0, 0.0}; - u_ = {1.0, 0.0, 0.0}; + void sample_collision_distance(double mfp) + { + collision_distance_ = -std::log(1.0 - xdg::tools::random::rand01(rng_state_)) * mfp; + } - volume_ = xdg_->find_volume(r_, u_); - log("Particle {} initialized in volume {}", id_, volume_); -} + double advance() + { + double distance; -void surf_dist() { - surface_intersection_ = xdg_->ray_fire(volume_, r_, u_, INFTY, HitOrientation::EXITING, &history_); - if (surface_intersection_.first == 0.0) { - fatal_error("Particle {} stuck at position ({}, {}, {}) on surfacce {}", id_, r_.x, r_.y, r_.z, surface_intersection_.second); - alive_ = false; - return; - } - if (surface_intersection_.second == ID_NONE) { - fatal_error("Particle {} lost in volume {}", id_, volume_); - alive_ = false; - return; - } - log("Intersected surface {} at distance {} ", surface_intersection_.second, surface_intersection_.first); -} + if (collision_distance_ < surface_hit_distance_) { + distance = collision_distance_; + } else { + distance = surface_hit_distance_; + } -void sample_collision_distance(double mfp) { - collision_distance_ = -std::log(1.0 - drand48()) * mfp; -} + // explicit scalar to avoid compilation issues with vec3da operator overloads on device + // TODO - Try via vec3da operator overload + r_.x += distance * u_.x; + r_.y += distance * u_.y; + r_.z += distance * u_.z; -void collide() { - n_events_++; - log("Event {} for particle {}", n_events_, id_); - u_ = rand_dir(); - log("Particle {} collides with material at position ({}, {}, {}), new direction is ({}, {}, {})", id_, r_.x, r_.y, r_.z, u_.z, u_.y, u_.z); - history_.clear(); -} + return distance; + } -void advance(std::unordered_map& cell_tracks) -{ - log("Comparing surface intersection distance {} to collision distance {}", surface_intersection_.first, collision_distance_); - if (collision_distance_ < surface_intersection_.first) { - r_ += collision_distance_ * u_; - cell_tracks[volume_] += collision_distance_; - log("Particle {} collides with material at position ({}, {}, {}) ", id_, r_.x, r_.y, r_.z); - } else { - r_ += surface_intersection_.first * u_; - cell_tracks[volume_] += surface_intersection_.first; - log("Particle {} advances to surface {} at position ({}, {}, {}) ", id_, surface_intersection_.second, r_.x, r_.y, r_.z); + void collide() + { + n_events_++; + + double direction[3]; + xdg::tools::random::random_unit_dir_lcg(rng_state_, direction); + u_.x = direction[0]; + u_.y = direction[1]; + u_.z = direction[2]; + + // reset surface-hit-state explicitly + surface_hit_ = ID_NONE; + surface_hit_distance_ = INFTY; + last_surface_hit_ = ID_NONE; } -} -void cross_surface() -{ - n_events_++; - log("Event {} for particle {}", n_events_, id_); - auto boundary_condition = xdg_->mesh_manager()->get_surface_property(surface_intersection_.second, PropertyType::BOUNDARY_CONDITION); - // check for the surface boundary condition - if (boundary_condition.value == "reflecting" || boundary_condition.value == "reflective") { - log("Particle {} reflects off surface {}", id_, surface_intersection_.second); - log("Direction before reflection: ({}, {}, {})", u_.x, u_.y, u_.z); - - Direction normal = xdg_->surface_normal(surface_intersection_.second, r_, &history_); - log("Normal to surface: ({}, {}, {})", normal.x, normal.y, normal.z); - - double proj = dot(normal, u_); - double mag = normal.length(); - normal = normal * (2.0 * proj/mag); - u_ = u_ - normal; - u_ = u_.normalize(); - log("Direction after reflection: ({}, {}, {})", u_.x, u_.y, u_.z); - // reset to last intersection - if (history_.size() > 0) { - log("Resetting particle history to last intersection"); - history_ = {history_.back()}; - } - } else if (boundary_condition.value == "vacuum") { - log("Particle {} encounters vacuum boundary at surface {}", id_, surface_intersection_.second); - alive_ = false; - } else { - volume_ = xdg_->mesh_manager()->next_volume(volume_, surface_intersection_.second); - log("Particle {} enters volume {}", id_, volume_); - if (ipc_graveyard_ && volume_ == xdg_->mesh_manager()->implicit_complement()) volume_ = ID_NONE; - if (volume_ == ID_NONE) { - alive_ = false; - return; - } + void mark_surface_hit(MeshID surface, double distance) + { + surface_hit_ = surface; + surface_hit_distance_ = distance; + last_surface_hit_ = surface; } -} -// Data Members -bool verbose_ {true}; -std::shared_ptr xdg_; -uint32_t id_ {0}; -int32_t max_events_ {1000}; -bool ipc_graveyard_ {false}; - -Position r_; -Direction u_; -MeshID volume_ {ID_NONE}; -std::vector history_ {}; -std::pair surface_intersection_ {INFTY, ID_NONE}; -double collision_distance_ {INFTY}; -int32_t n_events_ {0}; -bool alive_ {true}; + uint32_t id_ {0}; + Position r_ {}; + Direction u_ {}; + MeshID volume_ {ID_NONE}; + + MeshID surface_hit_ {ID_NONE}; + double surface_hit_distance_ {INFTY}; + double collision_distance_ {INFTY}; + MeshID last_surface_hit_ {ID_NONE}; + + std::uint32_t rng_state_ {0}; + int32_t n_events_ {0}; + bool alive_ {true}; }; +#ifdef _OPENMP +#pragma omp end declare target +#endif + struct EventSimulationData { std::shared_ptr xdg_; double mfp_ {1.0}; - uint32_t n_particles_ {100}; + std::uint32_t seed_ {42}; + uint32_t n_particles_ {100000}; uint32_t max_events_ {1000}; bool verbose_particles_ {false}; bool implicit_complement_is_graveyard_ {false}; std::unordered_map cell_tracks; uint32_t max_particles_in_flight_ {100000}; - std::vector particles; + EventParticle* device_particles {nullptr}; + int gpu_id {0}; + int host_id {omp_get_initial_device()}; + XDGRayHitBuffer ray_hits; ParticleEventQueue advance_particle_queue; ParticleEventQueue surface_crossing_queue; ParticleEventQueue collision_queue; @@ -284,13 +268,26 @@ inline void transport_particle_event_based(EventSimulationData& sim_data) { */ - sim_data.particles.clear(); - sim_data.particles.reserve(sim_data.n_particles_); + if (sim_data.n_particles_ == 0) { + fatal_error("Number of event particles must be greater than 0"); + } + + if (sim_data.n_particles_ > sim_data.max_particles_in_flight_) { + fatal_error("Event particle refill is not implemented; n_particles must be <= max_particles_in_flight"); + } + + sim_data.gpu_id = omp_get_default_device(); + sim_data.host_id = omp_get_initial_device(); + + sim_data.device_particles = static_cast( + omp_target_alloc(sim_data.n_particles_ * sizeof(EventParticle), sim_data.gpu_id)); + if (!sim_data.device_particles) { + fatal_error("Failed to allocate event particles on OpenMP target device."); + } - const int gpu_id = omp_get_default_device(); - sim_data.advance_particle_queue.allocate(sim_data.n_particles_, gpu_id); - sim_data.surface_crossing_queue.allocate(sim_data.n_particles_, gpu_id); - sim_data.collision_queue.allocate(sim_data.n_particles_, gpu_id); + sim_data.advance_particle_queue.allocate(sim_data.n_particles_, sim_data.gpu_id); + sim_data.surface_crossing_queue.allocate(sim_data.n_particles_, sim_data.gpu_id); + sim_data.collision_queue.allocate(sim_data.n_particles_, sim_data.gpu_id); process_init_events(sim_data); diff --git a/tools/random_lcg.h b/tools/random_lcg.h new file mode 100644 index 00000000..aed5135d --- /dev/null +++ b/tools/random_lcg.h @@ -0,0 +1,43 @@ +#ifndef _XDG_TOOLS_RANDOM_LCG_H +#define _XDG_TOOLS_RANDOM_LCG_H + +#include +#include + +namespace xdg::tools::random { + +#ifdef _OPENMP +#pragma omp declare target +#endif + +inline double rand01(std::uint32_t& state) +{ + state = state * 1664525u + 1013904223u; + return static_cast(state) * (1.0 / 4294967296.0); +} + +inline void random_unit_dir_lcg(std::uint32_t& state, double direction[3]) +{ + double x1; + double x2; + double s; + + do { + x1 = rand01(state) * 2.0 - 1.0; + x2 = rand01(state) * 2.0 - 1.0; + s = x1 * x1 + x2 * x2; + } while (s <= 0.0 || s >= 1.0); + + const double t = 2.0 * std::sqrt(1.0 - s); + direction[0] = x1 * t; + direction[1] = x2 * t; + direction[2] = 1.0 - 2.0 * s; +} + +#ifdef _OPENMP +#pragma omp end declare target +#endif + +} // namespace xdg::tools::random + +#endif // _XDG_TOOLS_RANDOM_LCG_H diff --git a/tools/ray_benchmark.h b/tools/ray_benchmark.h index 65a48339..3ebf7794 100644 --- a/tools/ray_benchmark.h +++ b/tools/ray_benchmark.h @@ -4,6 +4,8 @@ #include #include +#include "random_lcg.h" + namespace xdg::tools::benchmark { struct SourceSample { @@ -15,30 +17,6 @@ struct SourceSample { #pragma omp declare target #endif -inline double rand01(std::uint32_t& state) -{ - state = state * 1664525u + 1013904223u; - return static_cast(state) * (1.0 / 4294967296.0); -} - -inline void random_unit_dir_lcg(std::uint32_t& state, double direction[3]) -{ - double x1; - double x2; - double s; - - do { - x1 = rand01(state) * 2.0 - 1.0; - x2 = rand01(state) * 2.0 - 1.0; - s = x1 * x1 + x2 * x2; - } while (s <= 0.0 || s >= 1.0); - - const double t = 2.0 * std::sqrt(1.0 - s); - direction[0] = x1 * t; - direction[1] = x2 * t; - direction[2] = 1.0 - 2.0 * s; -} - inline SourceSample random_spherical_source(double origin_x, double origin_y, double origin_z, @@ -46,14 +24,14 @@ inline SourceSample random_spherical_source(double origin_x, double source_radius) { SourceSample sample; - random_unit_dir_lcg(state, sample.direction); + xdg::tools::random::random_unit_dir_lcg(state, sample.direction); sample.position[0] = origin_x; sample.position[1] = origin_y; sample.position[2] = origin_z; if (source_radius > 0.0) { - const double radius = source_radius * std::cbrt(rand01(state)); + const double radius = source_radius * std::cbrt(xdg::tools::random::rand01(state)); sample.position[0] += sample.direction[0] * radius; sample.position[1] += sample.direction[1] * radius; sample.position[2] += sample.direction[2] * radius; From c7d04857efe10359c58bb425fae86bc1de0ad75a Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 13 Jul 2026 15:47:33 +0100 Subject: [PATCH 18/23] Implemented process_init_events --- tools/particle_sim_event.h | 53 ++++++++++++++++++++++++++++---------- 1 file changed, 40 insertions(+), 13 deletions(-) diff --git a/tools/particle_sim_event.h b/tools/particle_sim_event.h index 31347cf5..9e69c54e 100644 --- a/tools/particle_sim_event.h +++ b/tools/particle_sim_event.h @@ -317,22 +317,49 @@ inline void transport_particle_event_based(EventSimulationData& sim_data) { // #ifdef OPENMC_MPI // MPI_Barrier( mpi::intracomm ); // #endif + + sim_data.advance_particle_queue.release(); + sim_data.surface_crossing_queue.release(); + sim_data.collision_queue.release(); + + omp_target_free(sim_data.device_particles, sim_data.gpu_id); + sim_data.device_particles = nullptr; } +inline void process_init_events(EventSimulationData& sim_data) +{ + // Match the history-based miniapp source initialization. + Position r {0.0, 0.0, 0.0}; + Direction u {1.0, 0.0, 0.0}; + MeshID volume = sim_data.xdg_->find_volume(r, u); + + EventParticle* device_particles = sim_data.device_particles; + const int n_particles = static_cast(sim_data.n_particles_); + const std::uint32_t seed = sim_data.seed_; + const int gpu_id = sim_data.gpu_id; + auto advance_queue = sim_data.advance_particle_queue.get_device_data(); + + if (!device_particles) { + fatal_error("Error allocating event particle device storage."); + } + + sim_data.advance_particle_queue.reset(); + sim_data.surface_crossing_queue.reset(); + sim_data.collision_queue.reset(); + + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(device_particles) \ + firstprivate(advance_queue, seed, r, u, volume) + for (int i = 0; i < n_particles; ++i) { + device_particles[i].initialize(static_cast(i), seed, r, u, volume); + advance_queue.thread_safe_append({static_cast(i)}); // no need for xs lookup so we just append particle to queue + } + + sim_data.advance_particle_queue.sync_size_device_to_host(); // ensure host side event scheduler knows the correct queue size +} + + -/* - initialize position/direction - find starting volume (call this globally for all particles [in flight]) - set global particle id - reset event counter - reset alive flag - reset ray history - reset pending surface hit / collision distance - initialize per-particle RNG seed if you stop using global drand48() - enqueue into advance queue - - Essentially set up particle state (for all particles) ready to call advance_particle -*/ // void process_death_events(); /* From 93db09aba2e95b717f521caa899ee030050138e0 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 13 Jul 2026 17:41:50 +0100 Subject: [PATCH 19/23] Implemented process_advance_particle_events() --- tools/particle_sim_event.h | 109 ++++++++++++++++++++++++++++++++++++- 1 file changed, 106 insertions(+), 3 deletions(-) diff --git a/tools/particle_sim_event.h b/tools/particle_sim_event.h index 9e69c54e..c0bbb6e6 100644 --- a/tools/particle_sim_event.h +++ b/tools/particle_sim_event.h @@ -160,6 +160,7 @@ struct EventParticle { rng_state_ = seed ^ id; n_events_ = 0; alive_ = true; + stuck_ = false; } void sample_collision_distance(double mfp) @@ -202,7 +203,7 @@ struct EventParticle { last_surface_hit_ = ID_NONE; } - void mark_surface_hit(MeshID surface, double distance) + void store_surface_hit(MeshID surface, double distance) { surface_hit_ = surface; surface_hit_distance_ = distance; @@ -222,6 +223,7 @@ struct EventParticle { std::uint32_t rng_state_ {0}; int32_t n_events_ {0}; bool alive_ {true}; + bool stuck_ {false}; }; #ifdef _OPENMP @@ -276,7 +278,8 @@ inline void transport_particle_event_based(EventSimulationData& sim_data) { fatal_error("Event particle refill is not implemented; n_particles must be <= max_particles_in_flight"); } - sim_data.gpu_id = omp_get_default_device(); + sim_data.ray_hits = sim_data.xdg_->allocate_ray_hits(sim_data.n_particles_); + sim_data.gpu_id = sim_data.ray_hits.device_id; sim_data.host_id = omp_get_initial_device(); sim_data.device_particles = static_cast( @@ -324,6 +327,8 @@ inline void transport_particle_event_based(EventSimulationData& sim_data) { omp_target_free(sim_data.device_particles, sim_data.gpu_id); sim_data.device_particles = nullptr; + + sim_data.xdg_->free_ray_hits(sim_data.ray_hits); } inline void process_init_events(EventSimulationData& sim_data) @@ -358,7 +363,105 @@ inline void process_init_events(EventSimulationData& sim_data) sim_data.advance_particle_queue.sync_size_device_to_host(); // ensure host side event scheduler knows the correct queue size } +inline void process_advance_particle_events(EventSimulationData& sim_data) +{ + const int n_advance = sim_data.advance_particle_queue.size(); + + if (n_advance == 0) + { + // TODO - Once things are definitely working we can probably remove this warning/check + warning("Advance_particle_events launched with a queue size of 0. Early return called..."); + return; + } + + EventParticle* device_particles = sim_data.device_particles; + XDGRayHit* ray_hits = sim_data.ray_hits.data; + auto advance_queue = sim_data.advance_particle_queue.get_device_data(); + auto surface_crossing_queue = sim_data.surface_crossing_queue.get_device_data(); + auto collision_queue = sim_data.collision_queue.get_device_data(); + const double mfp = sim_data.mfp_; + const int gpu_id = sim_data.gpu_id; + + // rayhit packing kernel + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(device_particles, ray_hits) \ + firstprivate(advance_queue) + for (int i = 0; i < n_advance; ++i) { + const uint32_t particle_idx = advance_queue.data[i].idx; + EventParticle& p = device_particles[particle_idx]; + + ray_hits[i].origin[0] = p.r_.x; + ray_hits[i].origin[1] = p.r_.y; + ray_hits[i].origin[2] = p.r_.z; + ray_hits[i].direction[0] = p.u_.x; + ray_hits[i].direction[1] = p.u_.y; + ray_hits[i].direction[2] = p.u_.z; + ray_hits[i].t_min = 0.0; + ray_hits[i].t_max = INFTY; + ray_hits[i].volume = p.volume_; + ray_hits[i].distance = INFTY; + ray_hits[i].surface = ID_NONE; + ray_hits[i].primitive = ID_NONE; + ray_hits[i].point_in_volume = OUTSIDE; + } + + + /* + XDG ray_fire_batch is host-orchestrated but operates on a device-resident + XDGRayHitBuffer. The particle kernels pack rays into that device buffer, + the host calls xdg->ray_fire_batch(), and a later particle kernel consumes + the hit results from the same device buffer. + + This differs from OpenMC's GPU path, where Particle::event_advance() calls + the device-callable distance-to-boundary logic directly inside one OpenMP + target kernel. Here we instead call a ray packing kernel followed by xdg's + ray_fire kernel and then the particle advance kernel but doing so should + allow us to benefit from GPU accelerated ray tracing against the CAD. + */ + + // Create a view over the actively packed portion of the preallocated ray-hit buffer + XDGRayHitBuffer active_hits {sim_data.ray_hits.data, + static_cast(n_advance), + sim_data.ray_hits.device_id}; + + sim_data.xdg_->ray_fire_batch(active_hits); // Perform GPU-accelerated ray tracing via XDG backend + + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(device_particles, ray_hits) \ + firstprivate(advance_queue, surface_crossing_queue, collision_queue, mfp) + for (int i = 0; i < n_advance; i++) { + const uint32_t particle_idx = advance_queue.data[i].idx; + EventParticle& p = device_particles[particle_idx]; + const XDGRayHit& hit = ray_hits[i]; + + // Set particle state to stuck. Also killed for now but we could tally the number of stuck particles later + if (hit.distance == 0.0) { + p.alive_ = false; + p.stuck_ = true; + continue; + } + + // Set particle state to killed + if (hit.surface == ID_NONE) { + p.alive_ = false; + continue; + } + p.store_surface_hit(hit.surface, hit.distance); + p.sample_collision_distance(mfp); + p.advance(); + + if (p.collision_distance_ < p.surface_hit_distance_) { + collision_queue.thread_safe_append({particle_idx}); + } else { + surface_crossing_queue.thread_safe_append({particle_idx}); + } + } + + sim_data.surface_crossing_queue.sync_size_device_to_host(); + sim_data.collision_queue.sync_size_device_to_host(); + sim_data.advance_particle_queue.reset(); +} // void process_death_events(); @@ -367,5 +470,5 @@ inline void process_init_events(EventSimulationData& sim_data) We are essentially just setting the particle.alive_ member to false. I think a better approach is to actually just handle death as it happens. Rather than waiting for everything a particle's death state should update when it - occurs. So we don't both with a process_death_events() method. + occurs. So we don't bother with a process_death_events() method. */ From ce850705a6e038f489d64a9b451cef1a32b7edf9 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 13 Jul 2026 17:59:48 +0100 Subject: [PATCH 20/23] Implemented process_collision_events() --- tools/particle_sim_event.h | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/tools/particle_sim_event.h b/tools/particle_sim_event.h index c0bbb6e6..e8df7726 100644 --- a/tools/particle_sim_event.h +++ b/tools/particle_sim_event.h @@ -463,6 +463,41 @@ inline void process_advance_particle_events(EventSimulationData& sim_data) sim_data.advance_particle_queue.reset(); } +inline void process_collision_events(EventSimulationData& sim_data) +{ + EventParticle* device_particles = sim_data.device_particles; + const int n_collisions = sim_data.collision_queue.size(); + + if (n_collisions == 0) { + return; + } + + auto advance_queue = sim_data.advance_particle_queue.get_device_data(); + auto collision_queue = sim_data.collision_queue.get_device_data(); + const int gpu_id = sim_data.gpu_id; + const int max_events = sim_data.max_events_; + + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(device_particles) \ + firstprivate(advance_queue, collision_queue, max_events) + for (int i = 0; i < n_collisions; i++) + { + const uint32_t particle_idx = collision_queue.data[i].idx; + EventParticle& p = device_particles[particle_idx]; + + p.collide(); + + if (p.n_events_ >= max_events) { + p.alive_ = false; + continue; + } + + advance_queue.thread_safe_append({particle_idx}); + } + + sim_data.advance_particle_queue.sync_size_device_to_host(); + sim_data.collision_queue.reset(); +} // void process_death_events(); /* From dcd15727f1be4014bda1d1442ec4659c6b4c4e21 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 14 Jul 2026 15:05:01 +0100 Subject: [PATCH 21/23] Implemented process_surface_crossing_events() --- tools/particle_sim_event.h | 128 +++++++++++++++++++++++++++++++++++-- 1 file changed, 123 insertions(+), 5 deletions(-) diff --git a/tools/particle_sim_event.h b/tools/particle_sim_event.h index e8df7726..f7c603c2 100644 --- a/tools/particle_sim_event.h +++ b/tools/particle_sim_event.h @@ -157,6 +157,11 @@ struct EventParticle { surface_hit_distance_ = INFTY; collision_distance_ = INFTY; last_surface_hit_ = ID_NONE; + next_volume_ = ID_NONE; + boundary_condition_ = UNSET; + surface_normal_.x = 0.0; + surface_normal_.y = 0.0; + surface_normal_.z = 0.0; rng_state_ = seed ^ id; n_events_ = 0; alive_ = true; @@ -201,15 +206,79 @@ struct EventParticle { surface_hit_ = ID_NONE; surface_hit_distance_ = INFTY; last_surface_hit_ = ID_NONE; + next_volume_ = ID_NONE; + boundary_condition_ = UNSET; + surface_normal_.x = 0.0; + surface_normal_.y = 0.0; + surface_normal_.z = 0.0; } - void store_surface_hit(MeshID surface, double distance) + void store_surface_crossing(MeshID surface, + double distance, + MeshID next_volume, + SurfaceBoundaryCondition boundary_condition, + const double normal[3]) { surface_hit_ = surface; surface_hit_distance_ = distance; last_surface_hit_ = surface; + next_volume_ = next_volume; + boundary_condition_ = boundary_condition; + surface_normal_.x = normal[0]; + surface_normal_.y = normal[1]; + surface_normal_.z = normal[2]; } + void surface_cross() + { + n_events_++; + + switch (boundary_condition_) { + case SurfaceBoundaryCondition::TRANSMISSION: + volume_ = next_volume_; + + // TODO: Restore optional implicit-complement graveyard handling. + if (volume_ == ID_NONE) { + alive_ = false; + } + break; + + case SurfaceBoundaryCondition::VACUUM: + alive_ = false; + break; + + case SurfaceBoundaryCondition::REFLECTIVE: { + const double nx = surface_normal_.x; + const double ny = surface_normal_.y; + const double nz = surface_normal_.z; + + // cuBQL returns the raw triangle normal. Dividing by n dot n applies + // the reflection formula without first normalising meaning we skip a square root. + const double normal_squared = nx * nx + ny * ny + nz * nz; + const double projection = u_.x * nx + u_.y * ny + u_.z * nz; + const double scale = 2.0 * projection / normal_squared; + + u_.x -= scale * nx; + u_.y -= scale * ny; + u_.z -= scale * nz; + + // Normalize the reflected particle direction. + const double direction_squared = u_.x * u_.x + u_.y * u_.y + u_.z * u_.z; + const double inverse_direction_length = 1.0 / std::sqrt(direction_squared); + u_.x *= inverse_direction_length; + u_.y *= inverse_direction_length; + u_.z *= inverse_direction_length; + break; + } + + case SurfaceBoundaryCondition::UNSET: + default: + alive_ = false; + break; + } + } + + // Data members uint32_t id_ {0}; Position r_ {}; Direction u_ {}; @@ -219,6 +288,9 @@ struct EventParticle { double surface_hit_distance_ {INFTY}; double collision_distance_ {INFTY}; MeshID last_surface_hit_ {ID_NONE}; + MeshID next_volume_ {ID_NONE}; + SurfaceBoundaryCondition boundary_condition_ {UNSET}; + Direction surface_normal_ {0.0}; std::uint32_t rng_state_ {0}; int32_t n_events_ {0}; @@ -234,13 +306,13 @@ struct EventSimulationData { std::shared_ptr xdg_; double mfp_ {1.0}; std::uint32_t seed_ {42}; - uint32_t n_particles_ {100000}; + uint32_t n_particles_ {1000000}; uint32_t max_events_ {1000}; bool verbose_particles_ {false}; bool implicit_complement_is_graveyard_ {false}; std::unordered_map cell_tracks; - uint32_t max_particles_in_flight_ {100000}; + uint32_t max_particles_in_flight_ {1000000}; EventParticle* device_particles {nullptr}; int gpu_id {0}; int host_id {omp_get_initial_device()}; @@ -396,13 +468,19 @@ inline void process_advance_particle_events(EventSimulationData& sim_data) ray_hits[i].direction[0] = p.u_.x; ray_hits[i].direction[1] = p.u_.y; ray_hits[i].direction[2] = p.u_.z; - ray_hits[i].t_min = 0.0; + // Avoid immediately re-hitting the surface from which the particle starts. + ray_hits[i].t_min = TINY_BIT; ray_hits[i].t_max = INFTY; ray_hits[i].volume = p.volume_; ray_hits[i].distance = INFTY; ray_hits[i].surface = ID_NONE; ray_hits[i].primitive = ID_NONE; ray_hits[i].point_in_volume = OUTSIDE; + ray_hits[i].next_volume = ID_NONE; + ray_hits[i].boundary_condition = static_cast(UNSET); + ray_hits[i].normal[0] = 0.0; + ray_hits[i].normal[1] = 0.0; + ray_hits[i].normal[2] = 0.0; } @@ -447,7 +525,12 @@ inline void process_advance_particle_events(EventSimulationData& sim_data) continue; } - p.store_surface_hit(hit.surface, hit.distance); + p.store_surface_crossing( + hit.surface, + hit.distance, + hit.next_volume, + static_cast(hit.boundary_condition), + hit.normal); p.sample_collision_distance(mfp); p.advance(); @@ -499,6 +582,41 @@ inline void process_collision_events(EventSimulationData& sim_data) sim_data.collision_queue.reset(); } +inline void process_surface_crossing_events(EventSimulationData& sim_data) +{ + EventParticle* device_particles = sim_data.device_particles; + const int n_surface_crossings = sim_data.surface_crossing_queue.size(); + + if (n_surface_crossings == 0) { + return; + } + + auto advance_queue = sim_data.advance_particle_queue.get_device_data(); + auto surface_crossing_queue = sim_data.surface_crossing_queue.get_device_data(); + const int gpu_id = sim_data.gpu_id; + const int max_events = sim_data.max_events_; + + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(device_particles) \ + firstprivate(advance_queue, surface_crossing_queue, max_events) + for (int i = 0; i < n_surface_crossings; i++) + { + const uint32_t particle_idx = surface_crossing_queue.data[i].idx; + EventParticle& p = device_particles[particle_idx]; + p.surface_cross(); + + if (!p.alive_ || p.n_events_ >= max_events) { + p.alive_ = false; + continue; + } + + advance_queue.thread_safe_append({particle_idx}); + } + + sim_data.advance_particle_queue.sync_size_device_to_host(); + sim_data.surface_crossing_queue.reset(); +} + // void process_death_events(); /* Particle death in the pseudo transport app doesn't really mean all that much. From 6fae809be25fecd491fe34bc3dd7d18b96395dd9 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 14 Jul 2026 15:05:20 +0100 Subject: [PATCH 22/23] Wired up new miniapp in CMake --- tools/CMakeLists.txt | 1 + tools/particle_sim_event.cpp | 18 ++++++++---------- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 883cbe84..a2acd610 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -1,5 +1,6 @@ set(TOOL_NAMES particle_sim +particle_sim_event ray_fire ray_benchmark find_volume diff --git a/tools/particle_sim_event.cpp b/tools/particle_sim_event.cpp index 47d1328f..b2f15308 100644 --- a/tools/particle_sim_event.cpp +++ b/tools/particle_sim_event.cpp @@ -9,7 +9,7 @@ #include "argparse/argparse.hpp" -#include "particle_sim.h" +#include "particle_sim_event.h" using namespace xdg; @@ -31,7 +31,7 @@ args.add_argument("-m", "--mfp") .help("Mean free path of the particles").scan<'g', double>(); args.add_argument("-n", "--n-particles") - .default_value(100u) + .default_value(1000000u) .help("Number of particles to simulate").scan<'u', uint32_t>(); args.add_argument("-e", "--max-events") @@ -48,8 +48,8 @@ args.add_argument("-m", "--mesh-library") .default_value("MOAB"); args.add_argument("-r", "--rt-library") - .help("Ray tracing library to use. One of (EMBREE, GPRT)") - .default_value("EMBREE"); + .help("Ray tracing library to use. Event transport currently requires CUBQL") + .default_value("CUBQL"); try { args.parse_args(argc, argv); } @@ -59,10 +59,7 @@ catch (const std::runtime_error& err) { exit(0); } -// Problem Setup -srand48(42); - -SimulationData sim_data; +EventSimulationData sim_data; // create a mesh manager std::string mesh_str = args.get("--mesh-library"); @@ -73,6 +70,8 @@ if (rt_str == "EMBREE") rt_lib = RTLibrary::EMBREE; else if (rt_str == "GPRT") rt_lib = RTLibrary::GPRT; +else if (rt_str == "CUBQL") + rt_lib = RTLibrary::CUBQL; else fatal_error("Invalid ray tracing library '{}' specified", rt_str); @@ -105,8 +104,7 @@ sim_data.implicit_complement_is_graveyard_ = args.get("--ipc-graveyard"); sim_data.n_particles_ = args.get("--n-particles"); sim_data.max_events_ = args.get("--max-events"); -transport_particles(sim_data); -// transport_particles_event_based(sim_data); +transport_particle_event_based(sim_data); // report distances in each cell in a table write_message("Cell Track Lengths"); From c01d47eba3164ae4b887ab90cada290d6aeb8db1 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 20 Jul 2026 16:37:38 +0100 Subject: [PATCH 23/23] Added profiling/timing statistics --- tools/particle_sim_event.cpp | 149 +++++++++++++++++++++++++++++++++-- tools/particle_sim_event.h | 60 ++++++++++++++ 2 files changed, 202 insertions(+), 7 deletions(-) diff --git a/tools/particle_sim_event.cpp b/tools/particle_sim_event.cpp index b2f15308..4521a0d9 100644 --- a/tools/particle_sim_event.cpp +++ b/tools/particle_sim_event.cpp @@ -1,8 +1,13 @@ +#include #include #include #include +#include + +#include #include "xdg/error.h" +#include "xdg/timer.h" #include "xdg/mesh_manager_interface.h" #include "xdg/vec3da.h" #include "xdg/xdg.h" @@ -50,6 +55,12 @@ args.add_argument("-m", "--mesh-library") args.add_argument("-r", "--rt-library") .help("Ray tracing library to use. Event transport currently requires CUBQL") .default_value("CUBQL"); + +args.add_argument("--format") + .default_value("human") + .choices("human", "csv") + .help("stdout format. Human readable (default) or csv"); + try { args.parse_args(argc, argv); } @@ -59,11 +70,17 @@ catch (const std::runtime_error& err) { exit(0); } +Timer wall_timer; +wall_timer.start(); + EventSimulationData sim_data; // create a mesh manager std::string mesh_str = args.get("--mesh-library"); std::string rt_str = args.get("--rt-library"); +const std::string model_filename = args.get("filename"); +const std::string model_name = std::filesystem::path(model_filename).filename().string(); +const std::string output_format = args.get("--format"); RTLibrary rt_lib; if (rt_str == "EMBREE") @@ -87,15 +104,33 @@ else fatal_error("Invalid mesh library '{}' specified", mesh_str); // create an XDG instance with the specified mesh and ray tracing library +Timer xdg_setup_timer; +xdg_setup_timer.start(); std::shared_ptr xdg = XDG::create(mesh_lib, rt_lib); const auto& mm = xdg->mesh_manager(); -mm->load_file(args.get("filename")); +mm->load_file(model_filename); mm->init(); mm->parse_metadata(); xdg->prepare_raytracer(); +xdg_setup_timer.stop(); +sim_data.profiling.xdg_setup_s += xdg_setup_timer.elapsed(); sim_data.xdg_ = xdg; +const int model_num_volumes = mm->num_volumes(); +const int model_num_surfaces = mm->num_surfaces(); +const int model_num_volume_elements = mm->num_volume_elements(); +const int model_num_vertices = mm->num_vertices(); +std::uint64_t model_num_surface_primitives = 0; +for (MeshID surface : mm->surfaces()) { + model_num_surface_primitives += static_cast(mm->num_surface_faces(surface)); +} + +std::uint64_t model_num_volume_surface_refs = 0; +for (MeshID volume : mm->volumes()) { + model_num_volume_surface_refs += static_cast(mm->num_volume_faces(volume)); +} + // update the mean free path sim_data.mfp_ = args.get("--mfp"); @@ -106,13 +141,113 @@ sim_data.max_events_ = args.get("--max-events"); transport_particle_event_based(sim_data); -// report distances in each cell in a table -write_message("Cell Track Lengths"); -write_message("-----------"); -for (const auto& [cell, dist] : sim_data.cell_tracks) { - write_message("Cell {}: {}", cell, dist); +wall_timer.stop(); +const double wall_time = wall_timer.elapsed(); +const auto& profiling = sim_data.profiling; +const double ray_throughput = + profiling.advance_ray_trace_s > 0.0 + ? static_cast(profiling.rays_traced) / profiling.advance_ray_trace_s + : 0.0; + +const std::vector csv_columns { + "model", + "mesh_library", + "rt_library", + "model_num_volumes", + "model_num_surfaces", + "model_num_volume_elements", + "model_num_vertices", + "model_num_surface_primitives", + "model_num_volume_surface_refs", + "n_particles", + "max_events", + "mean_free_path", + "wall_time_s", + "profile_xdg_setup_s", + "profile_transport_s", + "profile_advance_total_s", + "profile_advance_pack_rays_s", + "profile_advance_ray_trace_s", + "profile_advance_update_particles_s", + "profile_collision_s", + "profile_surface_crossing_s", + "profile_ray_batches", + "profile_collision_calls", + "profile_surface_crossing_calls", + "profile_rays_traced", + "profile_ray_throughput_rays_per_s" +}; + +const std::vector csv_values { + model_name, + mesh_str, + rt_str, + fmt::format("{}", model_num_volumes), + fmt::format("{}", model_num_surfaces), + fmt::format("{}", model_num_volume_elements), + fmt::format("{}", model_num_vertices), + fmt::format("{}", model_num_surface_primitives), + fmt::format("{}", model_num_volume_surface_refs), + fmt::format("{}", sim_data.n_particles_), + fmt::format("{}", sim_data.max_events_), + fmt::format("{}", sim_data.mfp_), + fmt::format("{}", wall_time), + fmt::format("{}", profiling.xdg_setup_s), + fmt::format("{}", profiling.transport_s), + fmt::format("{}", profiling.advance_total_s), + fmt::format("{}", profiling.advance_pack_rays_s), + fmt::format("{}", profiling.advance_ray_trace_s), + fmt::format("{}", profiling.advance_update_particles_s), + fmt::format("{}", profiling.collision_s), + fmt::format("{}", profiling.surface_crossing_s), + fmt::format("{}", profiling.advance_calls), + fmt::format("{}", profiling.collision_calls), + fmt::format("{}", profiling.surface_crossing_calls), + fmt::format("{}", profiling.rays_traced), + fmt::format("{}", ray_throughput) +}; + +if (output_format == "csv") { + std::cout << fmt::format("{}\n", fmt::join(csv_columns, ",")); + std::cout << fmt::format("{}\n", fmt::join(csv_values, ",")); +} else { + std::cout << "\nXDG event-based particle pseudo-simulation\n"; + std::cout << "----------------------------------------\n"; + std::cout << "Model : " << model_name << "\n"; + std::cout << "Mesh library : " << mesh_str << "\n"; + std::cout << "Ray tracing library : " << rt_str << "\n"; + std::cout << "Volumes : " << model_num_volumes << "\n"; + std::cout << "Surfaces : " << model_num_surfaces << "\n"; + std::cout << "Volume elements : " << model_num_volume_elements << "\n"; + std::cout << "Vertices : " << model_num_vertices << "\n"; + std::cout << "Surface primitives : " << model_num_surface_primitives << "\n"; + std::cout << "Volume surface refs : " << model_num_volume_surface_refs << "\n"; + std::cout << "Particles : " << sim_data.n_particles_ << "\n"; + std::cout << "Max events/particle : " << sim_data.max_events_ << "\n"; + std::cout << "Mean free path : " << sim_data.mfp_ << "\n"; + std::cout << "----------------------------------------\n"; + std::cout << "Full wall-clock time : " << wall_time << " s\n"; + std::cout << "XDG setup : " << profiling.xdg_setup_s << " s\n"; + std::cout << "Transport time : " << profiling.transport_s << " s\n"; + std::cout << "Advance total : " << profiling.advance_total_s + << " s (" << profiling.advance_calls << " calls)\n"; + std::cout << " Ray batches : " << profiling.advance_calls << "\n"; + std::cout << " Pack rays : " << profiling.advance_pack_rays_s << " s\n"; + std::cout << " Ray trace : " << profiling.advance_ray_trace_s << " s\n"; + std::cout << " Rays traced : " + << fmt::format("{:.6e}", static_cast(profiling.rays_traced)) << "\n"; + std::cout << " Ray throughput : " << ray_throughput << " rays/s\n"; + std::cout << " Update particles : " << profiling.advance_update_particles_s << " s\n"; + std::cout << "Collision events : " << profiling.collision_s + << " s (" << profiling.collision_calls << " calls)\n"; + std::cout << "Surface crossings : " << profiling.surface_crossing_s + << " s (" << profiling.surface_crossing_calls << " calls)\n"; + std::cout << "----------------------------------------\n"; + std::cout << "Cell track lengths\n"; + for (const auto& [cell, dist] : sim_data.cell_tracks) { + std::cout << "Cell " << cell << " : " << dist << "\n"; + } } -write_message("-----------"); return 0; diff --git a/tools/particle_sim_event.h b/tools/particle_sim_event.h index f7c603c2..ba48a55d 100644 --- a/tools/particle_sim_event.h +++ b/tools/particle_sim_event.h @@ -11,6 +11,7 @@ #include "xdg/error.h" #include "xdg/mesh_manager_interface.h" +#include "xdg/timer.h" #include "xdg/vec3da.h" #include "xdg/xdg.h" @@ -303,6 +304,21 @@ struct EventParticle { #endif struct EventSimulationData { + struct Profiling { + double xdg_setup_s {0.0}; + double transport_s {0.0}; + double advance_total_s {0.0}; + double advance_pack_rays_s {0.0}; + double advance_ray_trace_s {0.0}; + double advance_update_particles_s {0.0}; + double collision_s {0.0}; + double surface_crossing_s {0.0}; + std::uint64_t advance_calls {0}; + std::uint64_t collision_calls {0}; + std::uint64_t surface_crossing_calls {0}; + std::uint64_t rays_traced {0}; + }; + std::shared_ptr xdg_; double mfp_ {1.0}; std::uint32_t seed_ {42}; @@ -320,6 +336,7 @@ struct EventSimulationData { ParticleEventQueue advance_particle_queue; ParticleEventQueue surface_crossing_queue; ParticleEventQueue collision_queue; + Profiling profiling; }; void process_init_events(EventSimulationData& sim_data); @@ -328,6 +345,13 @@ void process_surface_crossing_events(EventSimulationData& sim_data); void process_collision_events(EventSimulationData& sim_data); inline void transport_particle_event_based(EventSimulationData& sim_data) { + const double xdg_setup_s = sim_data.profiling.xdg_setup_s; + sim_data.profiling = {}; + sim_data.profiling.xdg_setup_s = xdg_setup_s; + + Timer transport_timer; + transport_timer.start(); + // MPI will be needed for multi GPU // #ifdef OPENMC_MPI // MPI_Barrier( mpi::intracomm ); @@ -401,6 +425,9 @@ inline void transport_particle_event_based(EventSimulationData& sim_data) { sim_data.device_particles = nullptr; sim_data.xdg_->free_ray_hits(sim_data.ray_hits); + + transport_timer.stop(); + sim_data.profiling.transport_s += transport_timer.elapsed(); } inline void process_init_events(EventSimulationData& sim_data) @@ -446,6 +473,11 @@ inline void process_advance_particle_events(EventSimulationData& sim_data) return; } + Timer total_timer; + total_timer.start(); + sim_data.profiling.advance_calls++; + sim_data.profiling.rays_traced += static_cast(n_advance); + EventParticle* device_particles = sim_data.device_particles; XDGRayHit* ray_hits = sim_data.ray_hits.data; auto advance_queue = sim_data.advance_particle_queue.get_device_data(); @@ -454,6 +486,9 @@ inline void process_advance_particle_events(EventSimulationData& sim_data) const double mfp = sim_data.mfp_; const int gpu_id = sim_data.gpu_id; + Timer timer; + timer.start(); + // rayhit packing kernel #pragma omp target teams distribute parallel for device(gpu_id) \ is_device_ptr(device_particles, ray_hits) \ @@ -482,6 +517,8 @@ inline void process_advance_particle_events(EventSimulationData& sim_data) ray_hits[i].normal[1] = 0.0; ray_hits[i].normal[2] = 0.0; } + timer.stop(); + sim_data.profiling.advance_pack_rays_s += timer.elapsed(); /* @@ -502,8 +539,14 @@ inline void process_advance_particle_events(EventSimulationData& sim_data) static_cast(n_advance), sim_data.ray_hits.device_id}; + timer.reset(); + timer.start(); sim_data.xdg_->ray_fire_batch(active_hits); // Perform GPU-accelerated ray tracing via XDG backend + timer.stop(); + sim_data.profiling.advance_ray_trace_s += timer.elapsed(); + timer.reset(); + timer.start(); #pragma omp target teams distribute parallel for device(gpu_id) \ is_device_ptr(device_particles, ray_hits) \ firstprivate(advance_queue, surface_crossing_queue, collision_queue, mfp) @@ -540,10 +583,15 @@ inline void process_advance_particle_events(EventSimulationData& sim_data) surface_crossing_queue.thread_safe_append({particle_idx}); } } + timer.stop(); + sim_data.profiling.advance_update_particles_s += timer.elapsed(); sim_data.surface_crossing_queue.sync_size_device_to_host(); sim_data.collision_queue.sync_size_device_to_host(); sim_data.advance_particle_queue.reset(); + + total_timer.stop(); + sim_data.profiling.advance_total_s += total_timer.elapsed(); } inline void process_collision_events(EventSimulationData& sim_data) @@ -555,6 +603,10 @@ inline void process_collision_events(EventSimulationData& sim_data) return; } + Timer timer; + timer.start(); + sim_data.profiling.collision_calls++; + auto advance_queue = sim_data.advance_particle_queue.get_device_data(); auto collision_queue = sim_data.collision_queue.get_device_data(); const int gpu_id = sim_data.gpu_id; @@ -580,6 +632,8 @@ inline void process_collision_events(EventSimulationData& sim_data) sim_data.advance_particle_queue.sync_size_device_to_host(); sim_data.collision_queue.reset(); + timer.stop(); + sim_data.profiling.collision_s += timer.elapsed(); } inline void process_surface_crossing_events(EventSimulationData& sim_data) @@ -591,6 +645,10 @@ inline void process_surface_crossing_events(EventSimulationData& sim_data) return; } + Timer timer; + timer.start(); + sim_data.profiling.surface_crossing_calls++; + auto advance_queue = sim_data.advance_particle_queue.get_device_data(); auto surface_crossing_queue = sim_data.surface_crossing_queue.get_device_data(); const int gpu_id = sim_data.gpu_id; @@ -615,6 +673,8 @@ inline void process_surface_crossing_events(EventSimulationData& sim_data) sim_data.advance_particle_queue.sync_size_device_to_host(); sim_data.surface_crossing_queue.reset(); + timer.stop(); + sim_data.profiling.surface_crossing_s += timer.elapsed(); } // void process_death_events();