From 24dccf89933bbcf9ac05a6f55b23a96cedfca48a Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 14:44:10 +0100 Subject: [PATCH 01/21] 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 0c185556..62124d46 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 @@ -225,6 +242,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 @@ -285,7 +319,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) @@ -323,6 +362,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 dfba5d7d04ffcebe93aeed85e37df42e0390c27f Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 14:52:51 +0100 Subject: [PATCH 02/21] 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 caf43fd3..ba652afb 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 @@ -181,4 +183,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 75face29..603b8f2b 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 32bbd86c82daeb5475b213fa7cfea427d18ca512 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 14:57:42 +0100 Subject: [PATCH 03/21] 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 b2e6776d..b979d6d8 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 fde3538a00b0f9a62dbec1cc00d63d6039747686 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:07:32 +0100 Subject: [PATCH 04/21] 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 0b9ee4b339ba25544bb91b2a16e3eee841dfa1bc Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:08:50 +0100 Subject: [PATCH 05/21] 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 23b95f1d..11f8b668 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -55,6 +55,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 bb3b1451..b11a25d8 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 MockedTriTetMesh", "[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 501fa5f6..de6c0106 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 MockedTriTetMesh (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 b28c9ddbfaa0c38026ee42cfc1607a61342ffd49 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:09:38 +0100 Subject: [PATCH 06/21] 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 da5ffed717232eb592d517f6daed17ce229f443d Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:10:13 +0100 Subject: [PATCH 07/21] 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 4f71cdcdd14ed750b201620ee4449aad584a85cc Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:10:33 +0100 Subject: [PATCH 08/21] 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 88b78a6b2c2ef8484f5d594d79eee4d9c2d544cc Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 15:15:25 +0100 Subject: [PATCH 09/21] 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 cff9b8818e649339f898c3bc19b4d98bf7d8bfb2 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 22 Jun 2026 16:28:38 +0100 Subject: [PATCH 10/21] 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 3b5ccab96687fc790f8c53a176f46d1c4981a460 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 25 Jun 2026 13:15:33 +0100 Subject: [PATCH 11/21] 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 603b8f2b..30a84710 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 3f8735bc7bf430e6088a6010ad106681d548ed65 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 25 Jun 2026 17:19:47 +0100 Subject: [PATCH 12/21] 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 f7909f6961a7a083c2a0390e36d92c0efde2eed9 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 14 Jul 2026 11:51:35 +0100 Subject: [PATCH 13/21] 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 d0228d7e6cf3280a57eda55038555d6bc4ce0f14 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 17 Jul 2026 12:32:55 +0100 Subject: [PATCH 14/21] 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 154422506ec9c3554a935980461c680d124ef2ad Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 20 Jul 2026 14:39:49 +0100 Subject: [PATCH 15/21] 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 84cd932fae4d9105ff337ba586cb094effadd4c1 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 31 Jul 2026 12:32:26 +0100 Subject: [PATCH 16/21] Added test for batch ray fire --- include/xdg/util/rng.h | 35 +++++++++ tests/CMakeLists.txt | 4 ++ tests/test_xdg_batch_api.cpp | 133 +++++++++++++++++++++++++++++++++++ tools/ray_benchmark.h | 30 ++------ 4 files changed, 176 insertions(+), 26 deletions(-) create mode 100644 tests/test_xdg_batch_api.cpp diff --git a/include/xdg/util/rng.h b/include/xdg/util/rng.h index 539fba9c..e2069327 100644 --- a/include/xdg/util/rng.h +++ b/include/xdg/util/rng.h @@ -1,6 +1,8 @@ #ifndef XDG_UTIL_RNG_H #define XDG_UTIL_RNG_H +#include +#include #include namespace xdg { @@ -20,6 +22,39 @@ inline double rand_double(double min=0.0, double max=1.0) return min + (max - min) * dis(gen); } +#ifdef _OPENMP +#pragma omp declare target +#endif + +inline double lcg_rand01(std::uint32_t& state) +{ + state = state * 1664525u + 1013904223u; + return static_cast(state) * (1.0 / 4294967296.0); +} + +inline void random_unit_direction_lcg(std::uint32_t& state, + double direction[3]) +{ + double x1; + double x2; + double s; + + do { + x1 = lcg_rand01(state) * 2.0 - 1.0; + x2 = lcg_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 #endif // XDG_UTIL_RNG_H diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 11f8b668..fc7e0837 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -39,6 +39,10 @@ if (XDG_ENABLE_MOAB) list(APPEND TEST_NAMES test_ray_tracer_cross_check) endif() +if (XDG_ENABLE_CUBQL) + list(APPEND TEST_NAMES test_xdg_batch_api) +endif() + if (XDG_ENABLE_MOAB AND XDG_BUILD_TOOLS) list(APPEND TEST_NAMES test_overlap_check) endif() diff --git a/tests/test_xdg_batch_api.cpp b/tests/test_xdg_batch_api.cpp new file mode 100644 index 00000000..b4b9ae0e --- /dev/null +++ b/tests/test_xdg_batch_api.cpp @@ -0,0 +1,133 @@ +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "mesh_mock.h" +#include "util.h" +#include "xdg/constants.h" +#include "xdg/device_ray.h" +#include "xdg/util/rng.h" +#include "xdg/xdg.h" + +using namespace xdg; +using namespace xdg::test; + +constexpr double tolerance = 1e-6; + +TEST_CASE("XDG batch ray-hit buffer lifecycle", "[batch][cubql]") +{ + check_ray_tracer_supported(RTLibrary::CUBQL); + + auto mesh_manager = std::make_shared(false); + auto xdg = std::make_shared(mesh_manager, RTLibrary::CUBQL); + + const XDGRayHitBuffer empty_buffer = xdg->allocate_ray_hits(0); + REQUIRE(empty_buffer.data == nullptr); + REQUIRE(empty_buffer.count == 0); + REQUIRE(empty_buffer.device_id == -1); + + XDGRayHitBuffer ray_hits = xdg->allocate_ray_hits(64); + REQUIRE(ray_hits.data != nullptr); + REQUIRE(ray_hits.count == 64); + REQUIRE(ray_hits.device_id >= 0); + xdg->free_ray_hits(ray_hits); + REQUIRE(ray_hits.data == nullptr); + REQUIRE(ray_hits.count == 0); + REQUIRE(ray_hits.device_id == -1); +} + +TEST_CASE("XDG batch ray fire accepts an empty buffer", + "[rayfire][batch][cubql]") +{ + check_ray_tracer_supported(RTLibrary::CUBQL); + + auto mesh_manager = std::make_shared(false); + auto xdg = std::make_shared(mesh_manager, RTLibrary::CUBQL); + + REQUIRE_NOTHROW(xdg->ray_fire_batch({})); +} + +// The batch ray-fire API is currently implemented only by cuBQL. +// This test also populates its input buffer using OpenMP target offload, +// so support for another backend may require a different population path. +TEST_CASE("XDG batch ray fire matches scalar queries on MeshMock", + "[rayfire][batch][cubql]") +{ + check_ray_tracer_supported(RTLibrary::CUBQL); + + auto mesh_manager = std::make_shared(false); + mesh_manager->init(); + + auto xdg = std::make_shared(mesh_manager, RTLibrary::CUBQL); + xdg->prepare_raytracer(); + + const std::size_t batch_sizes[] = {1, 64}; + for (std::size_t num_rays : batch_sizes) { + DYNAMIC_SECTION("N=" << num_rays << " matches scalar") + { + XDGRayHitBuffer ray_hits = xdg->allocate_ray_hits(num_rays); + + REQUIRE(ray_hits.data != nullptr); + + XDGRayHit* d_ray_hits = ray_hits.data; + const int device_id = ray_hits.device_id; + const std::uint32_t seed = 12345u; + const MeshID volume = mesh_manager->volumes()[0]; + + #pragma omp target teams distribute parallel for device(device_id) \ + is_device_ptr(d_ray_hits) + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + std::uint32_t state = seed ^ static_cast(ray_id); + + XDGRayHit ray_hit {}; + ray_hit.origin[0] = 0.0; + ray_hit.origin[1] = 0.0; + ray_hit.origin[2] = 0.0; + random_unit_direction_lcg(state, ray_hit.direction); + ray_hit.t_min = 0.0; + ray_hit.t_max = INFTY; + ray_hit.volume = volume; + ray_hit.last_hit_primitive = ID_NONE; + d_ray_hits[ray_id] = ray_hit; + } + + xdg->ray_fire_batch(ray_hits); // Perform device side ray tracing + + std::vector batch_results(num_rays); + const int copy_status = omp_target_memcpy(batch_results.data(), + ray_hits.data, + num_rays * sizeof(XDGRayHit), + 0, + 0, + omp_get_initial_device(), + ray_hits.device_id); + REQUIRE(copy_status == 0); // Ensure omp memcpy was successful before checking results + + xdg->free_ray_hits(ray_hits); + + // Compare cuBQL batch results with cuBQL scalar ray_fire results + // Note that cuBQL scalar ray_fire is cross-checked against embree in test_ray_tracer_cross_check.cpp + const Position origin {0.0, 0.0, 0.0}; + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + std::uint32_t state = seed ^ static_cast(ray_id); + double direction_values[3]; + random_unit_direction_lcg(state, direction_values); + const Direction direction {direction_values[0], direction_values[1], direction_values[2]}; + const auto scalar_result = xdg->ray_fire(volume, origin, direction); + + INFO("ray index: " << ray_id); + INFO("direction: " << direction); + INFO("batch result: surface=" << batch_results[ray_id].surface << ", distance=" << std::setprecision(17) << batch_results[ray_id].distance); + INFO("scalar result: surface=" << scalar_result.second << ", distance=" << std::setprecision(17) << scalar_result.first); + REQUIRE(batch_results[ray_id].surface == scalar_result.second); + REQUIRE_THAT(batch_results[ray_id].distance, Catch::Matchers::WithinAbs(scalar_result.first, tolerance)); + } + } + } +} diff --git a/tools/ray_benchmark.h b/tools/ray_benchmark.h index 65a48339..c255935d 100644 --- a/tools/ray_benchmark.h +++ b/tools/ray_benchmark.h @@ -4,6 +4,8 @@ #include #include +#include "xdg/util/rng.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); + random_unit_direction_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(lcg_rand01(state)); sample.position[0] += sample.direction[0] * radius; sample.position[1] += sample.direction[1] * radius; sample.position[2] += sample.direction[2] * radius; From f8953440611f7f17d936c6e4c4954021aca43341 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 31 Jul 2026 16:55:19 +0100 Subject: [PATCH 17/21] Extended batch api tests to include a mult-volume test with full payload checking --- tests/test_xdg_batch_api.cpp | 132 ++++++++++++++++++++++++++++++++++- 1 file changed, 130 insertions(+), 2 deletions(-) diff --git a/tests/test_xdg_batch_api.cpp b/tests/test_xdg_batch_api.cpp index b4b9ae0e..99f937df 100644 --- a/tests/test_xdg_batch_api.cpp +++ b/tests/test_xdg_batch_api.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -67,9 +68,9 @@ TEST_CASE("XDG batch ray fire matches scalar queries on MeshMock", auto xdg = std::make_shared(mesh_manager, RTLibrary::CUBQL); xdg->prepare_raytracer(); - const std::size_t batch_sizes[] = {1, 64}; + const std::size_t batch_sizes[] = {1, 1024}; for (std::size_t num_rays : batch_sizes) { - DYNAMIC_SECTION("N=" << num_rays << " matches scalar") + DYNAMIC_SECTION("N=" << num_rays) { XDGRayHitBuffer ray_hits = xdg->allocate_ray_hits(num_rays); @@ -131,3 +132,130 @@ TEST_CASE("XDG batch ray fire matches scalar queries on MeshMock", } } } + +TEST_CASE("XDG batch ray fire multi-volume", "[rayfire][batch][cubql]") +{ + check_ray_tracer_supported(RTLibrary::CUBQL); + + auto xdg = XDG::create(MeshLibrary::MOAB, RTLibrary::CUBQL); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("pwr_pincell.h5m"); + mesh_manager->init(); + mesh_manager->parse_metadata(); + xdg->prepare_raytracer(); + + struct VolumeCase { + MeshID volume; + Position origin; + }; + + const std::array volume_cases {{ + {1, {0.0, 0.01, 0.1}}, + {2, {0.42, 0.01, 0.1}}, + {3, {0.50, 0.01, 0.1}} + }}; + + for (const auto& volume_case : volume_cases) { + CAPTURE(volume_case.volume, volume_case.origin); + REQUIRE(xdg->point_in_volume(volume_case.volume, volume_case.origin)); + } + + const std::array axis_aligned_directions {{ + {1.0, 0.0, 0.0}, + {-1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, -1.0, 0.0}, + {0.0, 0.0, 1.0}, + {0.0, 0.0, -1.0} + }}; + + // Build host-side ray-hit buffer and copy to device to handle every ray in one batch + std::vector host_ray_hits; + host_ray_hits.reserve(volume_cases.size() * axis_aligned_directions.size()); + + for (const auto& direction : axis_aligned_directions) { + for (const auto& volume_case : volume_cases) { + XDGRayHit ray_hit {}; + ray_hit.origin[0] = volume_case.origin.x; + ray_hit.origin[1] = volume_case.origin.y; + ray_hit.origin[2] = volume_case.origin.z; + ray_hit.direction[0] = direction.x; + ray_hit.direction[1] = direction.y; + ray_hit.direction[2] = direction.z; + ray_hit.t_min = 0.0; + ray_hit.t_max = INFTY; + ray_hit.volume = volume_case.volume; + ray_hit.last_hit_primitive = ID_NONE; + host_ray_hits.push_back(ray_hit); + } + } + + XDGRayHitBuffer ray_hits = xdg->allocate_ray_hits(host_ray_hits.size()); + const std::size_t buffer_size = host_ray_hits.size() * sizeof(XDGRayHit); + REQUIRE(omp_target_memcpy(ray_hits.data, + host_ray_hits.data(), + buffer_size, + 0, + 0, + ray_hits.device_id, + omp_get_initial_device()) == 0); + + xdg->ray_fire_batch(ray_hits); + + REQUIRE(omp_target_memcpy(host_ray_hits.data(), + ray_hits.data, + buffer_size, + 0, + 0, + omp_get_initial_device(), + ray_hits.device_id) == 0); + + xdg->free_ray_hits(ray_hits); + + for (std::size_t ray_id = 0; ray_id < host_ray_hits.size(); ++ray_id) { + const auto& batch_ray_hit = host_ray_hits[ray_id]; // recover rayhit for this ray from batch results + + const Position origin {batch_ray_hit.origin[0], batch_ray_hit.origin[1], batch_ray_hit.origin[2]}; + const Direction direction {batch_ray_hit.direction[0], batch_ray_hit.direction[1], batch_ray_hit.direction[2]}; + std::vector last_hit_prim; + const auto scalar_hit = xdg->ray_fire(batch_ray_hit.volume, origin, direction, INFTY, HitOrientation::EXITING, &last_hit_prim); + + REQUIRE(scalar_hit.second != ID_NONE); + REQUIRE(last_hit_prim.size() == 1); + + // Get expected values for this ray from scalar ray_fire and mesh manager for direct comparison against batch ray_fire results + const MeshID expected_surface = scalar_hit.second; + const double expected_distance = scalar_hit.first; + const MeshID expected_primitive = last_hit_prim.back(); + const PointInVolume expected_point_in_volume = xdg->point_in_volume(batch_ray_hit.volume, origin, &direction) ? INSIDE : OUTSIDE; + const MeshID expected_next_volume = mesh_manager->next_volume(batch_ray_hit.volume, expected_surface); + const auto boundary_property = mesh_manager->get_surface_property(expected_surface, PropertyType::BOUNDARY_CONDITION); + SurfaceBoundaryCondition expected_boundary_condition = UNSET; + if (boundary_property.value == "transmission") { + expected_boundary_condition = TRANSMISSION; + } else if (boundary_property.value == "reflecting") { + expected_boundary_condition = REFLECTIVE; + } + const auto vertices = mesh_manager->face_vertices(expected_primitive); + const Direction expected_normal = (vertices[1] - vertices[0]).cross(vertices[2] - vertices[0]); + + INFO("ray index: " << ray_id); + INFO("volume: " << batch_ray_hit.volume); + INFO("origin: " << origin); + INFO("direction: " << direction); + INFO("batch ray hit: surface=" << batch_ray_hit.surface << ", distance=" << std::setprecision(17) << batch_ray_hit.distance); + INFO("scalar hit: surface=" << expected_surface << ", distance=" << std::setprecision(17) << expected_distance); + INFO("boundary property: " << boundary_property.value); + + REQUIRE(expected_boundary_condition != UNSET); + REQUIRE(batch_ray_hit.surface == expected_surface); + REQUIRE_THAT(batch_ray_hit.distance, Catch::Matchers::WithinAbs(expected_distance, tolerance)); + REQUIRE(batch_ray_hit.primitive == expected_primitive); + REQUIRE(batch_ray_hit.point_in_volume == expected_point_in_volume); + REQUIRE(batch_ray_hit.next_volume == expected_next_volume); + REQUIRE(batch_ray_hit.boundary_condition == expected_boundary_condition); + REQUIRE_THAT(batch_ray_hit.normal[0], Catch::Matchers::WithinAbs(expected_normal.x, tolerance)); + REQUIRE_THAT(batch_ray_hit.normal[1], Catch::Matchers::WithinAbs(expected_normal.y, tolerance)); + REQUIRE_THAT(batch_ray_hit.normal[2], Catch::Matchers::WithinAbs(expected_normal.z, tolerance)); + } +} From e0a4bdfd30347629f5d2da75ea386bd1525ec68a Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 31 Jul 2026 17:15:23 +0100 Subject: [PATCH 18/21] Updated output on test failiure for batch ray fire --- tests/test_xdg_batch_api.cpp | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/tests/test_xdg_batch_api.cpp b/tests/test_xdg_batch_api.cpp index 99f937df..562d862d 100644 --- a/tests/test_xdg_batch_api.cpp +++ b/tests/test_xdg_batch_api.cpp @@ -133,7 +133,8 @@ TEST_CASE("XDG batch ray fire matches scalar queries on MeshMock", } } -TEST_CASE("XDG batch ray fire multi-volume", "[rayfire][batch][cubql]") +TEST_CASE("XDG batch ray fire multi-volume check full rayhit results", + "[rayfire][batch][cubql]") { check_ray_tracer_supported(RTLibrary::CUBQL); @@ -211,10 +212,10 @@ TEST_CASE("XDG batch ray fire multi-volume", "[rayfire][batch][cubql]") ray_hits.device_id) == 0); xdg->free_ray_hits(ray_hits); - + + // Loop over each ray and perform a scalar ray_fire to compare against the batch results for (std::size_t ray_id = 0; ray_id < host_ray_hits.size(); ++ray_id) { const auto& batch_ray_hit = host_ray_hits[ray_id]; // recover rayhit for this ray from batch results - const Position origin {batch_ray_hit.origin[0], batch_ray_hit.origin[1], batch_ray_hit.origin[2]}; const Direction direction {batch_ray_hit.direction[0], batch_ray_hit.direction[1], batch_ray_hit.direction[2]}; std::vector last_hit_prim; @@ -223,7 +224,6 @@ TEST_CASE("XDG batch ray fire multi-volume", "[rayfire][batch][cubql]") REQUIRE(scalar_hit.second != ID_NONE); REQUIRE(last_hit_prim.size() == 1); - // Get expected values for this ray from scalar ray_fire and mesh manager for direct comparison against batch ray_fire results const MeshID expected_surface = scalar_hit.second; const double expected_distance = scalar_hit.first; const MeshID expected_primitive = last_hit_prim.back(); @@ -239,14 +239,20 @@ TEST_CASE("XDG batch ray fire multi-volume", "[rayfire][batch][cubql]") const auto vertices = mesh_manager->face_vertices(expected_primitive); const Direction expected_normal = (vertices[1] - vertices[0]).cross(vertices[2] - vertices[0]); + // Outputting full rayhit state on failure INFO("ray index: " << ray_id); - INFO("volume: " << batch_ray_hit.volume); + INFO("input volume: " << batch_ray_hit.volume); INFO("origin: " << origin); INFO("direction: " << direction); - INFO("batch ray hit: surface=" << batch_ray_hit.surface << ", distance=" << std::setprecision(17) << batch_ray_hit.distance); - INFO("scalar hit: surface=" << expected_surface << ", distance=" << std::setprecision(17) << expected_distance); - INFO("boundary property: " << boundary_property.value); - + INFO("surface: batch=" << batch_ray_hit.surface << ", scalar=" << expected_surface); + INFO("distance: batch=" << std::setprecision(17) << batch_ray_hit.distance << ", scalar=" << expected_distance); + INFO("primitive: batch=" << batch_ray_hit.primitive << ", scalar=" << expected_primitive); + INFO("point in volume: batch=" << batch_ray_hit.point_in_volume << ", scalar=" << expected_point_in_volume); + INFO("next volume: batch=" << batch_ray_hit.next_volume << ", scalar=" << expected_next_volume); + INFO("boundary condition: batch=" << batch_ray_hit.boundary_condition << ", scalar=" << expected_boundary_condition << " (" << boundary_property.value << ")"); + INFO("normal: batch=[" << batch_ray_hit.normal[0] << " " << batch_ray_hit.normal[1] << " " << batch_ray_hit.normal[2] << "], scalar=" << expected_normal); + + // Perform assertions per ray REQUIRE(expected_boundary_condition != UNSET); REQUIRE(batch_ray_hit.surface == expected_surface); REQUIRE_THAT(batch_ray_hit.distance, Catch::Matchers::WithinAbs(expected_distance, tolerance)); From de86dbd556f6034553b0dee9f06d26102f48bd23 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 3 Aug 2026 12:56:35 +0100 Subject: [PATCH 19/21] Fix header file include of mesh_mocks after quad+hex rebase --- tests/test_files | 2 +- tests/test_xdg_batch_api.cpp | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/tests/test_files b/tests/test_files index aa41be5e..323951e2 160000 --- a/tests/test_files +++ b/tests/test_files @@ -1 +1 @@ -Subproject commit aa41be5ea237951d3e9787e61c55c10984efea65 +Subproject commit 323951e29d44277bb74fab98d8632dc32635883b diff --git a/tests/test_xdg_batch_api.cpp b/tests/test_xdg_batch_api.cpp index 562d862d..8c56599e 100644 --- a/tests/test_xdg_batch_api.cpp +++ b/tests/test_xdg_batch_api.cpp @@ -9,7 +9,7 @@ #include #include -#include "mesh_mock.h" +#include "mesh_mocks.h" #include "util.h" #include "xdg/constants.h" #include "xdg/device_ray.h" @@ -74,8 +74,6 @@ TEST_CASE("XDG batch ray fire matches scalar queries on MeshMock", { XDGRayHitBuffer ray_hits = xdg->allocate_ray_hits(num_rays); - REQUIRE(ray_hits.data != nullptr); - XDGRayHit* d_ray_hits = ray_hits.data; const int device_id = ray_hits.device_id; const std::uint32_t seed = 12345u; From 4703daeeea00974ae3991cef15b626dc420b07fa Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 3 Aug 2026 12:58:59 +0100 Subject: [PATCH 20/21] Determine expected normals in batch_ray_fire_test via mesh_manager --- tests/test_xdg_batch_api.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/test_xdg_batch_api.cpp b/tests/test_xdg_batch_api.cpp index 8c56599e..b2322693 100644 --- a/tests/test_xdg_batch_api.cpp +++ b/tests/test_xdg_batch_api.cpp @@ -234,8 +234,7 @@ TEST_CASE("XDG batch ray fire multi-volume check full rayhit results", } else if (boundary_property.value == "reflecting") { expected_boundary_condition = REFLECTIVE; } - const auto vertices = mesh_manager->face_vertices(expected_primitive); - const Direction expected_normal = (vertices[1] - vertices[0]).cross(vertices[2] - vertices[0]); + const Direction expected_normal = mesh_manager->face_normal(expected_primitive); // Outputting full rayhit state on failure INFO("ray index: " << ray_id); From b857939d316fdb9ea005eedb49af9bac3e04356f Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 4 Aug 2026 14:28:40 +0100 Subject: [PATCH 21/21] Update test_files submodule --- tests/test_files | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_files b/tests/test_files index 323951e2..aa41be5e 160000 --- a/tests/test_files +++ b/tests/test_files @@ -1 +1 @@ -Subproject commit 323951e29d44277bb74fab98d8632dc32635883b +Subproject commit aa41be5ea237951d3e9787e61c55c10984efea65