From 565677f6af4f2695e24d69e0689387e292f31cc1 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 3 Jun 2026 16:23:33 +0100 Subject: [PATCH 01/44] Working ray_benchmark miniapp --- tools/ray_benchmark.cpp | 124 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index dea28351..6c49d063 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -1,16 +1,28 @@ #include #include +<<<<<<< HEAD #include +======= +>>>>>>> 1003ecb (Working ray_benchmark miniapp) #include #include #include #include #include +<<<<<<< HEAD #include "argparse/argparse.hpp" #include #include "xdg/config.h" +======= +#ifdef XDG_OPENMP +#include +#endif + +#include "argparse/argparse.hpp" + +>>>>>>> 1003ecb (Working ray_benchmark miniapp) #include "xdg/constants.h" #include "xdg/error.h" #include "xdg/timer.h" @@ -19,12 +31,19 @@ #include "ray_benchmark.h" +<<<<<<< HEAD +======= +>>>>>>> 1003ecb (Working ray_benchmark miniapp) using namespace xdg; int main(int argc, char** argv) { +<<<<<<< HEAD argparse::ArgumentParser args("XDG Raytracing throughput benchmarking tool", +======= + argparse::ArgumentParser args("XDG Ray Tracing throughput benchmarking tool", +>>>>>>> 1003ecb (Working ray_benchmark miniapp) "1.0", argparse::default_arguments::help); @@ -46,6 +65,7 @@ int main(int argc, char** argv) .scan<'u', std::uint32_t>(); args.add_argument("-o", "-p", "--origin", "--position") +<<<<<<< HEAD .help("Ray origin/position. Defaults to the center of the model bounding box") .scan<'g', double>() .nargs(3); @@ -55,6 +75,13 @@ int main(int argc, char** argv) .implicit_value(true) .help("Use the queried volume bounding box center as the ray origin"); +======= + .default_value(std::vector{0.0, 0.0, 0.0}) + .help("Ray origin/position") + .scan<'g', double>() + .nargs(3); + +>>>>>>> 1003ecb (Working ray_benchmark miniapp) args.add_argument("-m", "--mesh-library") .help("Mesh library to use. One of (MOAB, LIBMESH)") .default_value("MOAB"); @@ -73,11 +100,14 @@ int main(int argc, char** argv) .help("Radius of a scattered source around the origin") .scan<'g', double>(); +<<<<<<< HEAD args.add_argument("--format") .default_value("human") .choices("human", "csv") .help("stdout format. Human readable (default) or csv"); +======= +>>>>>>> 1003ecb (Working ray_benchmark miniapp) args.add_description( "Benchmarks ray-fire throughput for a selected mesh volume. A source " "position is provided and ray directions are randomly generated from it."); @@ -112,12 +142,19 @@ int main(int argc, char** argv) } const MeshID volume = args.get("volume"); +<<<<<<< HEAD const std::string model_filename = args.get("filename"); const std::string model_name = std::filesystem::path(model_filename).filename().string(); const std::size_t num_rays = args.get("--num-rays"); const std::uint32_t seed = args.get("--seed"); const double source_radius = args.get("--source-radius"); const std::string output_format = args.get("--format"); +======= + const std::size_t num_rays = args.get("--num-rays"); + const std::uint32_t seed = args.get("--seed"); + const Position origin = args.get>("--origin"); + const double source_radius = args.get("--source-radius"); +>>>>>>> 1003ecb (Working ray_benchmark miniapp) Timer wall_timer; Timer setup_timer; @@ -130,6 +167,7 @@ int main(int argc, char** argv) setup_timer.start(); std::shared_ptr xdg = XDG::create(mesh_lib, rt_lib); const auto& mesh_manager = xdg->mesh_manager(); +<<<<<<< HEAD mesh_manager->load_file(model_filename); mesh_manager->init(); @@ -222,6 +260,73 @@ int main(int argc, char** argv) const double trace_time = trace_timer.elapsed(); const double end_to_end_time = generation_time + trace_time; const double setup_time = setup_timer.elapsed(); +======= + mesh_manager->load_file(args.get("filename")); + mesh_manager->init(); + + if (args.get("--list")) { + for (auto mesh_volume : mesh_manager->volumes()) { + std::cout << mesh_volume << std::endl; + } + return 0; + } + + xdg->prepare_volume_for_raytracing(volume); + xdg->ray_tracing_interface()->init(); + setup_timer.stop(); + + if (rt_lib == RTLibrary::EMBREE) { + #ifdef XDG_OPENMP + rt_label += " (" + std::to_string(omp_get_max_threads()) + " CPU threads)"; + #endif + } + + std::cout << "Volume ID: " << volume + << " with: " << mesh_manager->num_volume_faces(volume) + << " faces" << std::endl; + std::cout << "Starting ray fire benchmark with " << num_rays + << " rays using " << rt_label << "\n" << std::endl; + std::cout << "XDG initialisation time = " + << setup_timer.elapsed() << "s" << std::endl; + + if (rt_lib == RTLibrary::EMBREE) { + // Generate random rays from source + generation_timer.start(); + std::vector origins(num_rays); + std::vector directions(num_rays); + + #pragma omp parallel for schedule(static) + 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(); + + #pragma omp parallel for schedule(static) + for (std::size_t i = 0; i < num_rays; ++i) { + (void) xdg->ray_fire(volume, origins[i], directions[i]); + } + + trace_timer.stop(); + } + + const double generation_time = generation_timer.elapsed(); + const double trace_time = trace_timer.elapsed(); + const double end_to_end_time = generation_time + trace_time; +>>>>>>> 1003ecb (Working ray_benchmark miniapp) const double trace_only_rps = trace_time > 0.0 ? static_cast(num_rays) / trace_time : 0.0; @@ -230,6 +335,7 @@ int main(int argc, char** argv) : 0.0; wall_timer.stop(); +<<<<<<< HEAD const double wall_time = wall_timer.elapsed(); const std::vector csv_columns { @@ -317,6 +423,24 @@ int main(int argc, char** argv) std::cout << "End-to-end throughput : " << end_to_end_rps << " rays/s\n"; std::cout << "Trace-only throughput : " << trace_only_rps << " rays/s\n"; } +======= + + std::cout << "Random ray generation time = " + << generation_time << "s" << std::endl; + std::cout << "Generation + tracing time = " + << end_to_end_time << "s" << std::endl; + std::cout << "End-to-end throughput = " + << end_to_end_rps << " rays/s" << std::endl; + std::cout << "Full wall-clock time = " + << wall_timer.elapsed() << "s (post-argparse)" << std::endl; + + std::cout << "----------------------------------------" << std::endl; + std::cout << "Ray tracing time (trace-only)= " + << trace_time << "s for " << num_rays << " rays" << std::endl; + std::cout << "Trace-only throughput = " + << trace_only_rps << " rays/s" << std::endl; + std::cout << "----------------------------------------" << std::endl; +>>>>>>> 1003ecb (Working ray_benchmark miniapp) return 0; } From 8305ab3e1ef9292d2d7b88873e3489b3ba6eb358 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 10 Jun 2026 10:12:21 +0100 Subject: [PATCH 02/44] Default origin to model bbox + add flag to use volume bbox --- tools/ray_benchmark.cpp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index 6c49d063..1b0a8879 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -152,7 +152,6 @@ int main(int argc, char** argv) ======= const std::size_t num_rays = args.get("--num-rays"); const std::uint32_t seed = args.get("--seed"); - const Position origin = args.get>("--origin"); const double source_radius = args.get("--source-radius"); >>>>>>> 1003ecb (Working ray_benchmark miniapp) @@ -271,6 +270,20 @@ int main(int argc, char** argv) return 0; } + const auto origin_arg = args.present>("--origin"); + bool use_volume_center = args.get("--volume-center"); + if (origin_arg && use_volume_center) { + warning("--volume-center enabled but an explicit origin was also provided. The explicit origin will be used and the volume center will be ignored."); + use_volume_center = false; + } + + Position origin = mesh_manager->global_bounding_box().center(); + if (origin_arg) { + origin = Position(origin_arg.value()); + } else if (use_volume_center) { + origin = mesh_manager->volume_bounding_box(volume).center(); + } + xdg->prepare_volume_for_raytracing(volume); xdg->ray_tracing_interface()->init(); setup_timer.stop(); From 71e9597d52d9efd7831c39780e5588fc39302523 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 10 Jun 2026 10:49:18 +0100 Subject: [PATCH 03/44] Added fmt delimitting + spacing for volume listing --- tools/ray_benchmark.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index 1b0a8879..1d7615cd 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -30,6 +30,8 @@ #include "xdg/xdg.h" #include "ray_benchmark.h" +#include + <<<<<<< HEAD @@ -264,9 +266,7 @@ int main(int argc, char** argv) mesh_manager->init(); if (args.get("--list")) { - for (auto mesh_volume : mesh_manager->volumes()) { - std::cout << mesh_volume << std::endl; - } + std::cout << "[" << fmt::format("{}", fmt::join(mesh_manager->volumes(), ", ")) << "]\n"; return 0; } From 817a554de41026c53854052e4a99da9a140cb650 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 10 Jun 2026 11:21:56 +0100 Subject: [PATCH 04/44] Attempting to make use of XDGConfig and failing... --- tools/ray_benchmark.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index 1d7615cd..f49b6fd0 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -21,6 +21,7 @@ #endif #include "argparse/argparse.hpp" +#include >>>>>>> 1003ecb (Working ray_benchmark miniapp) #include "xdg/constants.h" @@ -30,7 +31,6 @@ #include "xdg/xdg.h" #include "ray_benchmark.h" -#include <<<<<<< HEAD @@ -287,11 +287,9 @@ int main(int argc, char** argv) xdg->prepare_volume_for_raytracing(volume); xdg->ray_tracing_interface()->init(); setup_timer.stop(); - if (rt_lib == RTLibrary::EMBREE) { - #ifdef XDG_OPENMP - rt_label += " (" + std::to_string(omp_get_max_threads()) + " CPU threads)"; - #endif + rt_label += " (" + std::to_string(XDGConfig::config().n_threads()) + + " CPU threads)"; } std::cout << "Volume ID: " << volume From 4eb7af91e0b882937dd2b12bfc64269915599305 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 10 Jun 2026 13:34:30 +0100 Subject: [PATCH 05/44] Added a flag to specify human readable stdout or csv --- tools/ray_benchmark.cpp | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index f49b6fd0..a06f598a 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -292,13 +292,7 @@ int main(int argc, char** argv) + " CPU threads)"; } - std::cout << "Volume ID: " << volume - << " with: " << mesh_manager->num_volume_faces(volume) - << " faces" << std::endl; - std::cout << "Starting ray fire benchmark with " << num_rays - << " rays using " << rt_label << "\n" << std::endl; - std::cout << "XDG initialisation time = " - << setup_timer.elapsed() << "s" << std::endl; + const auto num_faces = mesh_manager->num_volume_faces(volume); if (rt_lib == RTLibrary::EMBREE) { // Generate random rays from source From a6842904e8938ad9320b691cf160a89a423bdf1c Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 10 Jun 2026 13:38:25 +0100 Subject: [PATCH 06/44] Set omp pragma schedule to runtime --- tools/ray_benchmark.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index a06f598a..3b1ae8d7 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -300,7 +300,7 @@ int main(int argc, char** argv) std::vector origins(num_rays); std::vector directions(num_rays); - #pragma omp parallel for schedule(static) + #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, @@ -320,7 +320,7 @@ int main(int argc, char** argv) // Trace rays trace_timer.start(); - #pragma omp parallel for schedule(static) + #pragma omp parallel for schedule(runtime) for (std::size_t i = 0; i < num_rays; ++i) { (void) xdg->ray_fire(volume, origins[i], directions[i]); } From c6d71fe36ae6ee98514c7956f1b59b32faacb81c Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 27 Apr 2026 13:09:07 +0100 Subject: [PATCH 07/44] Added cuBQL as a submodule --- .gitmodules | 3 +++ vendor/cuBQL | 1 + 2 files changed, 4 insertions(+) create mode 160000 vendor/cuBQL 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/vendor/cuBQL b/vendor/cuBQL new file mode 160000 index 00000000..5a651b37 --- /dev/null +++ b/vendor/cuBQL @@ -0,0 +1 @@ +Subproject commit 5a651b3787e0cfae123eb427088fa73921e5a5fe From ca62ec5ff9f13f3b98a5229c0c07bb45dd225a34 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 28 Apr 2026 15:05:36 +0100 Subject: [PATCH 08/44] Added required CMake wiring + a CMakePresets.json for offload flags --- CMakeLists.txt | 34 +++++++++++++++++++++++++++++++--- CMakePresets.json | 45 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+), 3 deletions(-) create mode 100644 CMakePresets.json diff --git a/CMakeLists.txt b/CMakeLists.txt index df2d6cb3..33209afd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,6 +7,7 @@ option(XDG_ENABLE_LIBMESH "Enable support for the libMesh mesh library" OFF option(XDG_LINK_MPI "Link with MPI (for dependency compatibility)" OFF) option(XDG_ENABLE_EMBREE "Enable support for the Embree ray tracing library" ON) option(XDG_ENABLE_GPRT "Enable support for the GPRT ray tracing library" OFF) +option(XDG_ENABLE_CUBQL "Enable support for the cuBQL ray tracing library" OFF) option(XDG_BUILD_TESTS "Enable C++ unit testing" ON) option(XDG_BUILD_TOOLS "Enable tools and miniapps" ON) @@ -20,6 +21,10 @@ if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose build type" FORCE) endif() +if(DEFINED XDG_CMAKE_PRESET) + message(STATUS "XDG CMake preset: ${XDG_CMAKE_PRESET}") +endif() + # Compiler options (things in this section may not be platform-portable) set(CMAKE_POSITION_INDEPENDENT_CODE ON) @@ -111,6 +116,12 @@ if(XDG_ENABLE_GPRT) ) endif() +if(XDG_ENABLE_CUBQL) + list(APPEND VENDOR_PATHS + vendor/cuBQL + ) +endif() + if(GIT_FOUND AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/.git") option(XDG_GIT_SUBMODULE "Check submodules during build" ON) if(XDG_GIT_SUBMODULE) @@ -175,11 +186,12 @@ if (NOT XDG_ENABLE_MOAB AND NOT XDG_ENABLE_LIBMESH) endif() # Ensure at least one ray tracing backend is enabled -if (NOT XDG_ENABLE_EMBREE AND NOT XDG_ENABLE_GPRT) +if (NOT XDG_ENABLE_EMBREE AND NOT XDG_ENABLE_GPRT AND NOT XDG_ENABLE_CUBQL) message(FATAL_ERROR "No ray tracing backend enabled. Enable at least one of:\n" " -DXDG_ENABLE_EMBREE=ON\n" - " -DXDG_ENABLE_GPRT=ON") + " -DXDG_ENABLE_GPRT=ON\n" + " -DXDG_ENABLE_CUBQL=ON") endif() # GPRT @@ -188,6 +200,11 @@ if (XDG_ENABLE_GPRT) add_subdirectory(vendor/GPRT) endif() +if (XDG_ENABLE_CUBQL) + set(CUBQL_OMP ON CACHE BOOL "Build cuBQL with OpenMP target offload" FORCE) + add_subdirectory(vendor/cuBQL) +endif() + list(APPEND xdg_sources src/geometry/measure.cpp src/geometry/plucker.cpp @@ -221,6 +238,12 @@ dbl_deviceCode endif() +if (XDG_ENABLE_CUBQL) +list(APPEND xdg_sources +src/cuBQL/ray_tracer.cpp +) +endif() + if (XDG_ENABLE_LIBMESH) list(APPEND xdg_sources src/libmesh/mesh_manager.cpp @@ -281,7 +304,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) diff --git a/CMakePresets.json b/CMakePresets.json new file mode 100644 index 00000000..d1e26cb8 --- /dev/null +++ b/CMakePresets.json @@ -0,0 +1,45 @@ +{ + "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" + } + } + ], + "buildPresets": [ + { + "name": "cubql_llvm_ada", + "configurePreset": "cubql_llvm_ada" + } + ] +} From a7f5f71966c3a486c8c1195afc807b036bade74a Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 28 Apr 2026 15:06:37 +0100 Subject: [PATCH 09/44] Added constructor setup for a cubql backend --- include/xdg/constants.h | 8 +++++--- include/xdg/ray_tracers.h | 6 +++++- src/xdg.cpp | 14 ++++++++++++++ 3 files changed, 24 insertions(+), 4 deletions(-) diff --git a/include/xdg/constants.h b/include/xdg/constants.h index 38df9db1..e143d537 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -54,7 +54,8 @@ enum class MeshLibrary { // Ray Tracing library identifier enum class RTLibrary { EMBREE, - GPRT + GPRT, + CUBQL }; static const std::map MESH_LIB_TO_STR = @@ -67,7 +68,8 @@ static const std::map MESH_LIB_TO_STR = static const std::map RT_LIB_TO_STR = { {RTLibrary::EMBREE, "EMBREE"}, - {RTLibrary::GPRT, "GPRT"} + {RTLibrary::GPRT, "GPRT"}, + {RTLibrary::CUBQL, "CUBQL"} }; // Mesh identifer type @@ -148,4 +150,4 @@ struct formatter : fmt::formatter { } -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/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/xdg.cpp b/src/xdg.cpp index cc738008..18737bf2 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -31,6 +31,14 @@ XDG::XDG(std::shared_ptr mesh_manager, RTLibrary ray_tracing_lib) #else fatal_error("This build was not compiled with GPRT support (XDG_ENABLE_GPRT=OFF)."); #endif + + case RTLibrary::CUBQL: + #ifdef XDG_ENABLE_CUBQL + set_ray_tracing_interface(std::make_shared()); + break; + #else + fatal_error("This build was not compiled with cuBQL support (XDG_ENABLE_CUBQL=OFF)."); + #endif } } @@ -84,6 +92,9 @@ std::shared_ptr XDG::create(MeshLibrary mesh_lib, RTLibrary ray_tracing_lib #ifdef XDG_ENABLE_GPRT if (ray_tracing_lib == RTLibrary::GPRT) return std::make_shared(); #endif + #ifdef XDG_ENABLE_CUBQL + if (ray_tracing_lib == RTLibrary::CUBQL) return std::make_shared(); + #endif // If no supported ray tracing library throw an error std::string msg = fmt::format("Invalid ray tracing library '{}'. Supported:", RT_LIB_TO_STR.at(ray_tracing_lib)); @@ -93,6 +104,9 @@ std::shared_ptr XDG::create(MeshLibrary mesh_lib, RTLibrary ray_tracing_lib #ifdef XDG_ENABLE_GPRT msg += " GPRT"; #endif + #ifdef XDG_ENABLE_CUBQL + msg += " CUBQL"; + #endif fatal_error(msg); }; From 6b6f0709a692fadf26260593703dc2b0409d6b5a Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 28 Apr 2026 15:07:12 +0100 Subject: [PATCH 10/44] Added cubql ray tracer header + implmentation stubs --- include/xdg/cuBQL/ray_tracer.h | 68 +++++++++++++++++++++ src/cuBQL/ray_tracer.cpp | 104 +++++++++++++++++++++++++++++++++ 2 files changed, 172 insertions(+) create mode 100644 include/xdg/cuBQL/ray_tracer.h create mode 100644 src/cuBQL/ray_tracer.cpp diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h new file mode 100644 index 00000000..eae80009 --- /dev/null +++ b/include/xdg/cuBQL/ray_tracer.h @@ -0,0 +1,68 @@ +#ifndef _XDG_CUBQL_RAY_TRACING_INTERFACE_H +#define _XDG_CUBQL_RAY_TRACING_INTERFACE_H + +#include +#include +#include + +#include "xdg/constants.h" +#include "xdg/geometry_data.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/ray.h" +#include "xdg/ray_tracing_interface.h" + + + +namespace xdg { + +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; + + std::pair closest(TreeID tree, + const Position& origin) override; + + bool occluded(TreeID tree, + const Position& origin, + const Direction& direction, + double& dist) const override; +}; + +} // namespace xdg + +#endif // include guard diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp new file mode 100644 index 00000000..05de9ea7 --- /dev/null +++ b/src/cuBQL/ray_tracer.cpp @@ -0,0 +1,104 @@ +#include "xdg/cuBQL/ray_tracer.h" +#include "xdg/error.h" + +#include "cuBQL/bvh.h" +#include "cuBQL/builder/omp.h" + +namespace xdg { + +CuBQLRayTracer::CuBQLRayTracer() = default; + +CuBQLRayTracer::~CuBQLRayTracer() = default; + +void CuBQLRayTracer::init() +{ +} + +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}; +} + +TreeID +CuBQLRayTracer::create_surface_tree(const std::shared_ptr&, + MeshID) +{ + // Testing cubql compilation + cuBQL::bvh_t bvh; + + + return TREE_NONE; +} + +TreeID +CuBQLRayTracer::create_element_tree(const std::shared_ptr&, + MeshID) +{ + fatal_error("Element trees not currently supported with cuBQL ray tracer"); + return TREE_NONE; +} + +void CuBQLRayTracer::create_global_surface_tree() +{ + fatal_error("Global surface trees not currently supported with cuBQL ray tracer"); +} + +void CuBQLRayTracer::create_global_element_tree() +{ + fatal_error("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, + const Position&, + const Direction*, + const std::vector*) const +{ + fatal_error("Point-in-volume queries not currently supported with cuBQL ray tracer"); + return false; +} + +std::pair +CuBQLRayTracer::ray_fire(TreeID, + const Position&, + const Direction&, + const double, + HitOrientation, + std::vector* const) +{ + fatal_error("Ray-fire queries not currently supported with cuBQL ray tracer"); + return {INFTY, ID_NONE}; +} + +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 From bd725f552213502189462d86a1fe9c59825e381b Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 28 Apr 2026 15:36:30 +0100 Subject: [PATCH 11/44] Added a nvhpc preset --- CMakePresets.json | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/CMakePresets.json b/CMakePresets.json index d1e26cb8..09b86db1 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -34,12 +34,35 @@ "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 e569649dfcb328c7fd63474ea724baea07ed0f5d Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 28 Apr 2026 17:27:44 +0100 Subject: [PATCH 12/44] Added some guards to ensure that cubql path stays on OMP and doesnt switch to cuda path --- src/cuBQL/ray_tracer.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 05de9ea7..fa9c3a02 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -1,6 +1,11 @@ #include "xdg/cuBQL/ray_tracer.h" #include "xdg/error.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 "cuBQL/bvh.h" #include "cuBQL/builder/omp.h" From f16ffe09b9a5a2da396fe9de6525bee96bd7faeb Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 29 Apr 2026 16:13:17 +0100 Subject: [PATCH 13/44] Added a TODO comment to fix the loading of libomptarget.so --- CMakeLists.txt | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 33209afd..046dcaa5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -348,6 +348,22 @@ if (XDG_ENABLE_GPRT) endif() target_link_libraries(xdg PRIVATE fmt::fmt) +if (XDG_ENABLE_CUBQL) + target_compile_definitions(xdg PUBLIC XDG_ENABLE_CUBQL) + target_link_libraries(xdg $) + # 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 fmt::fmt) + +#TODO: refactor to only be required when tools are enabled +target_link_libraries(xdg indicators::indicators) # ========================== # Link ray tracing libraries From ddc5cfd125588aca8ad34211d3ce68f2cdb59c2f Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 30 Apr 2026 12:53:01 +0100 Subject: [PATCH 14/44] Working towards creation of a surface BLAS from mesh manager data --- include/xdg/cuBQL/ray_tracer.h | 4 ++ src/cuBQL/ray_tracer.cpp | 75 +++++++++++++++++++++++++++++++--- 2 files changed, 74 insertions(+), 5 deletions(-) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index eae80009..b2950f51 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -2,6 +2,7 @@ #define _XDG_CUBQL_RAY_TRACING_INTERFACE_H #include +#include #include #include @@ -61,6 +62,9 @@ class CuBQLRayTracer : public RayTracer { const Position& origin, const Direction& direction, double& dist) const override; + +private: + }; } // namespace xdg diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index fa9c3a02..91f7e281 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -6,6 +6,7 @@ #undef __CUDA_ARCH__ #endif +#include #include "cuBQL/bvh.h" #include "cuBQL/builder/omp.h" @@ -29,13 +30,77 @@ CuBQLRayTracer::register_volume(const std::shared_ptr& mesh_manager } TreeID -CuBQLRayTracer::create_surface_tree(const std::shared_ptr&, - MeshID) +CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_manager, + MeshID volume_id) { - // Testing cubql compilation - cuBQL::bvh_t bvh; - + int gpu_id = 0; // TODO - how to manage GPU IDs in a multi-GPU system? + + cuBQL::box3d *d_boxes = nullptr; // box3d is an alias for box_t + int num_boxes = 0; + + cuBQL::bvh3d bvh; // bvh3d is an alias for BinaryBVH + // bvh_t is an alias for BinaryBVH, so bvh_t is also BinaryBVH + + cuBQL::BuildConfig buildParams; // It looks like spatial median is the only supported method for omp builder + // Looks like the cuBQL::gpuBuilder is only pulled in when cubql is built with CUDA support + + cuBQL::build_omp_target(bvh, d_boxes, num_boxes, buildParams, gpu_id); // build bvh on GPU 0 using OpenMP target offloading + + SurfaceTreeID tree = next_surface_tree_id(); + surface_trees_.push_back(tree); + auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); + + for (const auto &surf : volume_surfaces) { + auto num_faces = mesh_manager->num_surface_faces(surf); + auto vertices = mesh_manager->get_surface_vertices(surf); + auto indices = mesh_manager->get_surface_connectivity(surf); + + // Get host side storage for vertices using cubql friendly types + std::vector h_vertices; + h_vertices.reserve(vertices.size()); + for (const auto& vertex : vertices) { + h_vertices.emplace_back(vertex.x, vertex.y, vertex.z); + } + + // Get host side storage for indices using cubql friendly types + 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]); + } + + // copy vertices and indices to device + auto* d_vertices = static_cast + (omp_target_alloc(h_vertices.size() * sizeof(cuBQL::vec3d), gpu_id)); + + auto* d_indices = static_cast + (omp_target_alloc(h_indices.size() * sizeof(cuBQL::vec3i), gpu_id)); + + // Create device storage for triangle AABBs to be computed in parallel on GPU + auto* d_aabbs = static_cast + (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3d), gpu_id)); + + // Create AABBs for each triangle in parallel on GPU + // TODO - Should this be its own function? + // TODO - How can this be extended for tets and other element types? + #pragma omp target device(gpu_id) is_device_ptr(d_vertices, d_indices, d_aabbs) + #pragma omp teams distribute parallel for + for (uint32_t primID = 0; primID < num_faces; ++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); + + d_aabbs[primID] = aabb; + } + } return TREE_NONE; } From 3c832b003d122148704c07a5f5d0f563b2d62ada Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 30 Apr 2026 15:08:23 +0100 Subject: [PATCH 15/44] Copy vertices+indices to device and build bvh --- src/cuBQL/ray_tracer.cpp | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 91f7e281..aa8e724d 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -35,23 +35,27 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man { int gpu_id = 0; // TODO - how to manage GPU IDs in a multi-GPU system? + int host_id = omp_get_initial_device(); - cuBQL::box3d *d_boxes = nullptr; // box3d is an alias for box_t - int num_boxes = 0; - cuBQL::bvh3d bvh; // bvh3d is an alias for BinaryBVH - // bvh_t is an alias for BinaryBVH, so bvh_t is also BinaryBVH + // cuBQL::box3d *d_boxes = nullptr; // box3d is an alias for box_t + // int num_boxes = 0; + + // cuBQL::bvh3d bvh; // bvh3d is an alias for BinaryBVH + // // bvh_t is an alias for BinaryBVH, so bvh_t is also BinaryBVH cuBQL::BuildConfig buildParams; // It looks like spatial median is the only supported method for omp builder - // Looks like the cuBQL::gpuBuilder is only pulled in when cubql is built with CUDA support + // // Looks like the cuBQL::gpuBuilder is only pulled in when cubql is built with CUDA support - cuBQL::build_omp_target(bvh, d_boxes, num_boxes, buildParams, gpu_id); // build bvh on GPU 0 using OpenMP target offloading + // cuBQL::build_omp_target(bvh, d_boxes, num_boxes, buildParams, gpu_id); // build bvh on GPU 0 using OpenMP target offloading SurfaceTreeID tree = next_surface_tree_id(); surface_trees_.push_back(tree); auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); for (const auto &surf : volume_surfaces) { + cuBQL::bvh3d bvh; // bvh3d is an alias for BinaryBVH + auto num_faces = mesh_manager->num_surface_faces(surf); auto vertices = mesh_manager->get_surface_vertices(surf); auto indices = mesh_manager->get_surface_connectivity(surf); @@ -73,9 +77,13 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man // copy vertices and indices to device auto* d_vertices = static_cast (omp_target_alloc(h_vertices.size() * sizeof(cuBQL::vec3d), gpu_id)); + omp_target_memcpy(d_vertices, h_vertices.data(), + h_vertices.size() * sizeof(cuBQL::vec3d), 0, 0, gpu_id, host_id); auto* d_indices = static_cast (omp_target_alloc(h_indices.size() * sizeof(cuBQL::vec3i), gpu_id)); + omp_target_memcpy(d_indices, h_indices.data(), + h_indices.size() * sizeof(cuBQL::vec3i), 0, 0, gpu_id, host_id); // Create device storage for triangle AABBs to be computed in parallel on GPU auto* d_aabbs = static_cast @@ -86,7 +94,7 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man // TODO - How can this be extended for tets and other element types? #pragma omp target device(gpu_id) is_device_ptr(d_vertices, d_indices, d_aabbs) #pragma omp teams distribute parallel for - for (uint32_t primID = 0; primID < num_faces; ++primID) { + for (uint32_t primID = 0; primID < h_indices.size(); ++primID) { cuBQL::vec3i indices = d_indices[primID]; cuBQL::vec3d A = d_vertices[indices.x]; @@ -99,7 +107,11 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man aabb.extend(C); d_aabbs[primID] = aabb; - } + } // This is our AABB population kernel written with OpenMP target offloading pragmas to run in parallel on the GPU + + // Construct the bvh on the gpu using the AABBs with openmp pathway + cuBQL::build_omp_target(bvh, d_aabbs, num_faces, buildParams, gpu_id); + } return TREE_NONE; } From 7983ac078b5a0928c5d96629f746e896ac16e818 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 1 May 2026 15:32:08 +0100 Subject: [PATCH 16/44] Added objects for storage of bvhs --- include/xdg/cuBQL/ray_tracer.h | 10 ++++++++++ src/cuBQL/ray_tracer.cpp | 22 ++++++++++------------ 2 files changed, 20 insertions(+), 12 deletions(-) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index b2950f51..a44c57b6 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -12,6 +12,13 @@ #include "xdg/ray.h" #include "xdg/ray_tracing_interface.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 "cuBQL/bvh.h" + namespace xdg { @@ -64,6 +71,9 @@ class CuBQLRayTracer : public RayTracer { double& dist) const override; private: + std::vector surface_bvhs_; // BVH for each surface tree +// std::unordered_map surface_volume_tree_to_bvh_map; // Map from surface tree IDs to their corresponding cuBQL BVH structures + }; diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index aa8e724d..dc71cde8 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -1,15 +1,9 @@ #include "xdg/cuBQL/ray_tracer.h" #include "xdg/error.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 "cuBQL/bvh.h" #include "cuBQL/builder/omp.h" - +#include "cuBQL/traversal/shrinkingRadiusQuery.h" namespace xdg { CuBQLRayTracer::CuBQLRayTracer() = default; @@ -111,7 +105,7 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man // Construct the bvh on the gpu using the AABBs with openmp pathway cuBQL::build_omp_target(bvh, d_aabbs, num_faces, buildParams, gpu_id); - + surface_bvhs_.push_back(bvh); // Store the BVH for this surface tree } return TREE_NONE; } @@ -157,12 +151,16 @@ bool CuBQLRayTracer::point_in_volume(TreeID, std::pair CuBQLRayTracer::ray_fire(TreeID, - const Position&, - const Direction&, - const double, - HitOrientation, + const Position& origin, + const Direction& direction, + const double tmax, + HitOrientation hitOrientation, std::vector* const) { + + // For a closest hit query this is the traversal template we want to make use of + cuBQL::shrinkingRadiusQuery::forEachPrim(prim_lambda, node_lambda, surface_bvhs_[0], tmax) + fatal_error("Ray-fire queries not currently supported with cuBQL ray tracer"); return {INFTY, ID_NONE}; } From 49a400e1b999c718b054138aef4bb0691e3f0b56 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 1 May 2026 16:51:37 +0100 Subject: [PATCH 17/44] Prototyping with codex. Working ray-fire. MAKE SURE TO RESET TO PREVIOUS COMMIT --- include/xdg/cuBQL/ray_tracer.h | 26 ++- src/cuBQL/ray_tracer.cpp | 345 +++++++++++++++++++++++++++++++-- tests/test_ray_fire.cpp | 3 +- tests/util.h | 20 +- tools/ray_fire.cpp | 6 +- 5 files changed, 378 insertions(+), 22 deletions(-) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index a44c57b6..b12681f6 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -71,10 +71,28 @@ class CuBQLRayTracer : public RayTracer { double& dist) const override; private: - std::vector surface_bvhs_; // BVH for each surface tree -// std::unordered_map surface_volume_tree_to_bvh_map; // Map from surface tree IDs to their corresponding cuBQL BVH structures - - + struct CuBQLSurfaceBVH { + MeshID surface {ID_NONE}; + cuBQL::bvh3d bvh; + + cuBQL::vec3d* d_vertices {nullptr}; + cuBQL::vec3i* d_indices {nullptr}; + cuBQL::vec3d* d_normals {nullptr}; + MeshID* d_primitive_refs {nullptr}; + + uint32_t num_vertices {0}; + uint32_t num_faces {0}; + int gpu_id {0}; + }; + + struct CuBQLRayHit { + double distance {INFTY}; + MeshID surface {ID_NONE}; + MeshID primitive {ID_NONE}; + }; + + std::vector surface_bvhs_; + std::unordered_map> tree_to_surface_bvh_indices_; }; } // namespace xdg diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index dc71cde8..9da83d6a 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -3,12 +3,154 @@ #include #include "cuBQL/builder/omp.h" -#include "cuBQL/traversal/shrinkingRadiusQuery.h" +#include "cuBQL/math/Ray.h" +#include "cuBQL/queries/triangleData/Triangle.h" +#include "cuBQL/queries/triangleData/math/rayTriangleIntersections.h" +#include "cuBQL/traversal/rayQueries.h" + namespace xdg { +struct CuBQLPluckerIntersectionResult { + bool hit {false}; + double t {0.0}; +}; + +static inline __cubql_both double cubql_ray_hit_tolerance(double t) +{ + constexpr double tolerance = 64.0 * 2.2204460492503131e-16; + return tolerance * (1.0 + cuBQL::abst(t)); +} + +static inline __cubql_both bool cubql_plucker_first(cuBQL::vec3d a, + cuBQL::vec3d b) +{ + if (a.x < b.x) return true; + if (a.x > b.x) return false; + + if (a.y < b.y) return true; + if (a.y > b.y) return false; + + return a.z < b.z; +} + +static inline __cubql_both double cubql_plucker_edge_test(cuBQL::vec3d vertex_a, + cuBQL::vec3d vertex_b, + cuBQL::vec3d ray, + cuBQL::vec3d ray_normal) +{ + double pip; + if (cubql_plucker_first(vertex_a, vertex_b)) { + const cuBQL::vec3d edge = vertex_b - vertex_a; + const cuBQL::vec3d edge_normal = cross(edge, vertex_a); + pip = dot(ray, edge_normal) + dot(ray_normal, edge); + } else { + const cuBQL::vec3d edge = vertex_a - vertex_b; + const cuBQL::vec3d edge_normal = cross(edge, vertex_b); + pip = dot(ray, edge_normal) + dot(ray_normal, edge); + pip = -pip; + } + + constexpr double dbl_zero_tol = 20.0 * 2.2204460492503131e-16; + if (cuBQL::abst(pip) < dbl_zero_tol) { + pip = 0.0; + } + return pip; +} + +static inline __cubql_both CuBQLPluckerIntersectionResult +cubql_plucker_ray_tri_intersect(cuBQL::vec3d vertices[3], + cuBQL::vec3d origin, + cuBQL::vec3d direction, + double t_max, + double t_min) +{ + const cuBQL::vec3d ray_a = direction; + const cuBQL::vec3d ray_b = cross(direction, origin); + + const double plucker_coord_0 = + cubql_plucker_edge_test(vertices[0], vertices[1], ray_a, ray_b); + + const double plucker_coord_1 = + cubql_plucker_edge_test(vertices[1], vertices[2], ray_a, ray_b); + + if ((0.0 < plucker_coord_0 && 0.0 > plucker_coord_1) || + (0.0 > plucker_coord_0 && 0.0 < plucker_coord_1)) { + return {}; + } + + const double plucker_coord_2 = + cubql_plucker_edge_test(vertices[2], vertices[0], ray_a, ray_b); + + if ((0.0 < plucker_coord_1 && 0.0 > plucker_coord_2) || + (0.0 > plucker_coord_1 && 0.0 < plucker_coord_2) || + (0.0 < plucker_coord_0 && 0.0 > plucker_coord_2) || + (0.0 > plucker_coord_0 && 0.0 < plucker_coord_2)) { + return {}; + } + + if (plucker_coord_0 == 0.0 && + plucker_coord_1 == 0.0 && + plucker_coord_2 == 0.0) { + return {}; + } + + const double inverse_sum = + 1.0 / (plucker_coord_0 + plucker_coord_1 + plucker_coord_2); + + const cuBQL::vec3d intersection = + plucker_coord_0 * inverse_sum * vertices[2] + + plucker_coord_1 * inverse_sum * vertices[0] + + plucker_coord_2 * inverse_sum * vertices[1]; + + int idx = 0; + double max_abs_dir = 0.0; + for (int i = 0; i < 3; ++i) { + if (cuBQL::abst(direction[i]) > max_abs_dir) { + idx = i; + max_abs_dir = cuBQL::abst(direction[i]); + } + } + + double dist_out = (intersection[idx] - origin[idx]) / direction[idx]; + + const double u = plucker_coord_2 * inverse_sum; + const double v = plucker_coord_0 * inverse_sum; + if (u < 0.0 || v < 0.0 || (u + v) > 1.0) { + return {}; + } + + if (dist_out < t_min || dist_out > t_max) { + return {}; + } + + return {true, dist_out}; +} + CuBQLRayTracer::CuBQLRayTracer() = default; -CuBQLRayTracer::~CuBQLRayTracer() = default; +CuBQLRayTracer::~CuBQLRayTracer() +{ + for (auto& surface_bvh : surface_bvhs_) { + cuBQL::omp::Context context(surface_bvh.gpu_id); + + if (surface_bvh.bvh.nodes || surface_bvh.bvh.primIDs) { + cuBQL::omp::freeBVH(surface_bvh.bvh, &context); + } + + if (surface_bvh.d_vertices) { + omp_target_free(surface_bvh.d_vertices, surface_bvh.gpu_id); + } + if (surface_bvh.d_indices) { + omp_target_free(surface_bvh.d_indices, surface_bvh.gpu_id); + } + if (surface_bvh.d_normals) { + omp_target_free(surface_bvh.d_normals, surface_bvh.gpu_id); + } + if (surface_bvh.d_primitive_refs) { + omp_target_free(surface_bvh.d_primitive_refs, surface_bvh.gpu_id); + } + } +} void CuBQLRayTracer::init() { @@ -45,11 +187,10 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man SurfaceTreeID tree = next_surface_tree_id(); surface_trees_.push_back(tree); + auto& surface_bvh_indices = tree_to_surface_bvh_indices_[tree]; auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); for (const auto &surf : volume_surfaces) { - cuBQL::bvh3d bvh; // bvh3d is an alias for BinaryBVH - auto num_faces = mesh_manager->num_surface_faces(surf); auto vertices = mesh_manager->get_surface_vertices(surf); auto indices = mesh_manager->get_surface_connectivity(surf); @@ -68,6 +209,26 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man h_indices.emplace_back(indices[i], indices[i + 1], indices[i + 2]); } + // Get storage for primitive refs so a hit can be mapped back to a mesh face. + std::vector h_primitive_refs; + h_primitive_refs.reserve(num_faces); + + std::vector h_normals; + h_normals.reserve(num_faces); + + auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); + for (const auto& face : mesh_manager->get_surface_faces(surf)) { + h_primitive_refs.push_back(face); + + auto normal = mesh_manager->face_normal(face); + if (volume_id == reverse_parent) { + normal = -normal; + } else if (volume_id != forward_parent) { + fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); + } + h_normals.emplace_back(normal.x, normal.y, normal.z); + } + // copy vertices and indices to device auto* d_vertices = static_cast (omp_target_alloc(h_vertices.size() * sizeof(cuBQL::vec3d), gpu_id)); @@ -79,16 +240,27 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man omp_target_memcpy(d_indices, h_indices.data(), h_indices.size() * sizeof(cuBQL::vec3i), 0, 0, gpu_id, host_id); + auto* d_normals = static_cast + (omp_target_alloc(h_normals.size() * sizeof(cuBQL::vec3d), gpu_id)); + omp_target_memcpy(d_normals, h_normals.data(), + h_normals.size() * sizeof(cuBQL::vec3d), 0, 0, gpu_id, host_id); + + auto* d_primitive_refs = static_cast + (omp_target_alloc(h_primitive_refs.size() * sizeof(MeshID), gpu_id)); + omp_target_memcpy(d_primitive_refs, h_primitive_refs.data(), + h_primitive_refs.size() * sizeof(MeshID), 0, 0, gpu_id, host_id); + // Create device storage for triangle AABBs to be computed in parallel on GPU auto* d_aabbs = static_cast (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3d), gpu_id)); + const auto num_primitives = static_cast(h_indices.size()); // Create AABBs for each triangle in parallel on GPU // TODO - Should this be its own function? // TODO - How can this be extended for tets and other element types? #pragma omp target device(gpu_id) is_device_ptr(d_vertices, d_indices, d_aabbs) #pragma omp teams distribute parallel for - for (uint32_t primID = 0; primID < h_indices.size(); ++primID) { + for (uint32_t primID = 0; primID < num_primitives; ++primID) { cuBQL::vec3i indices = d_indices[primID]; cuBQL::vec3d A = d_vertices[indices.x]; @@ -104,17 +276,33 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man } // This is our AABB population kernel written with OpenMP target offloading pragmas to run in parallel on the GPU // Construct the bvh on the gpu using the AABBs with openmp pathway + cuBQL::bvh3d bvh; // bvh3d is an alias for BinaryBVH cuBQL::build_omp_target(bvh, d_aabbs, num_faces, buildParams, gpu_id); - surface_bvhs_.push_back(bvh); // Store the BVH for this surface tree + + omp_target_free(d_aabbs, gpu_id); + + CuBQLSurfaceBVH surface_bvh; + surface_bvh.surface = surf; + surface_bvh.bvh = bvh; + surface_bvh.d_vertices = d_vertices; + surface_bvh.d_indices = d_indices; + surface_bvh.d_normals = d_normals; + surface_bvh.d_primitive_refs = d_primitive_refs; + surface_bvh.num_vertices = h_vertices.size(); + surface_bvh.num_faces = num_faces; + surface_bvh.gpu_id = gpu_id; + + surface_bvh_indices.push_back(surface_bvhs_.size()); + surface_bvhs_.push_back(surface_bvh); } - return TREE_NONE; + return tree; } TreeID CuBQLRayTracer::create_element_tree(const std::shared_ptr&, MeshID) { - fatal_error("Element trees not currently supported with cuBQL ray tracer"); + warning("Element trees not currently supported with cuBQL ray tracer"); return TREE_NONE; } @@ -150,19 +338,148 @@ bool CuBQLRayTracer::point_in_volume(TreeID, } std::pair -CuBQLRayTracer::ray_fire(TreeID, +CuBQLRayTracer::ray_fire(TreeID tree, const Position& origin, const Direction& direction, const double tmax, HitOrientation hitOrientation, - std::vector* const) + std::vector* const exclude_primitives) { + int gpu_id = 0; // TODO - how to manage GPU IDs in a multi-GPU system? + int host_id = omp_get_initial_device(); - // For a closest hit query this is the traversal template we want to make use of - cuBQL::shrinkingRadiusQuery::forEachPrim(prim_lambda, node_lambda, surface_bvhs_[0], tmax) + 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, + host_id); + } - fatal_error("Ray-fire queries not currently supported with cuBQL ray tracer"); - return {INFTY, ID_NONE}; + double closest_distance = tmax; + MeshID closest_surface = ID_NONE; + MeshID closest_primitive = ID_NONE; + auto* d_surface_hit = static_cast + (omp_target_alloc(sizeof(CuBQLRayHit), gpu_id)); + + const auto& surface_bvh_indices = tree_to_surface_bvh_indices_.at(tree); + for (const auto surface_bvh_index : surface_bvh_indices) { + const auto& surface_bvh = surface_bvhs_.at(surface_bvh_index); + + cuBQL::bvh3d bvh = surface_bvh.bvh; + const cuBQL::vec3d* d_vertices = surface_bvh.d_vertices; + const cuBQL::vec3i* d_indices = surface_bvh.d_indices; + const cuBQL::vec3d* d_normals = surface_bvh.d_normals; + const MeshID* d_primitive_refs = surface_bvh.d_primitive_refs; + + CuBQLRayHit surface_hit; + surface_hit.distance = closest_distance + cubql_ray_hit_tolerance(closest_distance); + omp_target_memcpy(d_surface_hit, + &surface_hit, + sizeof(CuBQLRayHit), + 0, + 0, + gpu_id, + host_id); + + const int orientation = static_cast(hitOrientation); + const cuBQL::vec3d ray_origin(origin.x, origin.y, origin.z); + const cuBQL::vec3d ray_direction(direction.x, direction.y, direction.z); + + #pragma omp target device(gpu_id) \ + is_device_ptr(d_vertices, d_indices, d_normals, d_primitive_refs, d_exclude_primitives, d_surface_hit) + { + cuBQL::ray3d ray(ray_origin, ray_direction, 0.0, d_surface_hit->distance); + + auto intersect_prim = [=, &ray] + (uint32_t prim_id) -> double + { + const MeshID primitive_ref = d_primitive_refs[prim_id]; + + for (int i = 0; i < exclude_count; ++i) { + if (d_exclude_primitives[i] == primitive_ref) { + return ray.tMax; + } + } + + const cuBQL::vec3i index = d_indices[prim_id]; + cuBQL::vec3d vertices[3] = { + d_vertices[index.x], + d_vertices[index.y], + d_vertices[index.z] + }; + + const double normal_dot_direction = dot(d_normals[prim_id], ray.direction); + bool culled = false; + if (orientation == static_cast(HitOrientation::EXITING)) { + culled = normal_dot_direction < 0.0; + } else if (orientation == static_cast(HitOrientation::ENTERING)) { + culled = normal_dot_direction >= 0.0; + } + if (culled) { + return ray.tMax; + } + + auto intersection = cubql_plucker_ray_tri_intersect(vertices, + ray.origin, + ray.direction, + ray.tMax, + ray.tMin); + if (intersection.hit) { + d_surface_hit->distance = intersection.t; + d_surface_hit->primitive = primitive_ref; + ray.tMax = intersection.t; + } + + return ray.tMax; + }; + + cuBQL::shrinkingRayQuery::forEachPrim(intersect_prim, bvh, ray); + } + + omp_target_memcpy(&surface_hit, + d_surface_hit, + sizeof(CuBQLRayHit), + 0, + 0, + host_id, + gpu_id); + + const double hit_tolerance = cubql_ray_hit_tolerance(closest_distance); + const bool closer_hit = surface_hit.distance < closest_distance - hit_tolerance; + const bool tied_hit = cuBQL::abst(surface_hit.distance - closest_distance) <= hit_tolerance; + const bool lower_surface_tie = closest_surface == ID_NONE || + surface_bvh.surface < closest_surface; + + if (surface_hit.primitive != ID_NONE && (closer_hit || (tied_hit && lower_surface_tie))) { + closest_distance = surface_hit.distance; + closest_surface = surface_bvh.surface; + closest_primitive = surface_hit.primitive; + } + } + + omp_target_free(d_surface_hit, gpu_id); + + if (d_exclude_primitives) { + omp_target_free(d_exclude_primitives, gpu_id); + } + + if (closest_surface == ID_NONE) { + return {INFTY, ID_NONE}; + } + + if (exclude_primitives) { + exclude_primitives->push_back(closest_primitive); + } + + return {closest_distance, closest_surface}; } std::pair diff --git a/tests/test_ray_fire.cpp b/tests/test_ray_fire.cpp index be816c36..70a5353f 100644 --- a/tests/test_ray_fire.cpp +++ b/tests/test_ray_fire.cpp @@ -18,7 +18,8 @@ using namespace xdg::test; TEMPLATE_TEST_CASE("Ray Fire on MeshMock (per-backend sections)", "[rayfire][mock]", Embree_Raytracer, - GPRT_Raytracer) + GPRT_Raytracer, + CuBQL_Raytracer) { // Generate one test run per enabled backend constexpr auto rt_backend = TestType::value; diff --git a/tests/util.h b/tests/util.h index a9841086..86657d9d 100644 --- a/tests/util.h +++ b/tests/util.h @@ -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 true; + #else + return false; + #endif } return false; @@ -81,6 +88,12 @@ inline void check_ray_tracer_supported(xdg::RTLibrary rt) { if (!ray_tracer_available(rt)) { SKIP(fmt::format("{} backend unavailable; skipping.", xdg::RT_LIB_TO_STR.at(rt))); } + #endif + #ifndef XDG_ENABLE_CUBQL + if(!ray_tracer_available(xdg::RTLibrary::CUBQL)) { + SKIP(fmt::format("{} backend unavailable; skipping.", xdg::RT_LIB_TO_STR.at(rt))); + } + #endif } inline void check_mesh_library_supported(xdg::MeshLibrary mesh) { @@ -118,5 +131,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; } 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 812b0f20a72c242c17080d154b2687bceab3df47 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 6 May 2026 14:42:19 +0100 Subject: [PATCH 18/44] Added some additional comments and cleanups --- include/xdg/cuBQL/ray_tracer.h | 2 +- src/cuBQL/ray_tracer.cpp | 37 ++++++++++++++++++++++++---------- 2 files changed, 27 insertions(+), 12 deletions(-) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index b12681f6..f0a59b04 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -85,7 +85,7 @@ class CuBQLRayTracer : public RayTracer { int gpu_id {0}; }; - struct CuBQLRayHit { + struct CubqlHit { double distance {INFTY}; MeshID surface {ID_NONE}; MeshID primitive {ID_NONE}; diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 9da83d6a..bc9680ac 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -220,6 +220,9 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man for (const auto& face : mesh_manager->get_surface_faces(surf)) { h_primitive_refs.push_back(face); + // TODO: This pre-orients normals for the current volume tree. If cuBQL + // starts reusing a single surface BVH/BLAS across both parent volumes, + // normals must be flipped at intersection time based on the queried volume. auto normal = mesh_manager->face_normal(face); if (volume_id == reverse_parent) { normal = -normal; @@ -366,8 +369,8 @@ CuBQLRayTracer::ray_fire(TreeID tree, double closest_distance = tmax; MeshID closest_surface = ID_NONE; MeshID closest_primitive = ID_NONE; - auto* d_surface_hit = static_cast - (omp_target_alloc(sizeof(CuBQLRayHit), gpu_id)); + auto* d_surface_hit = static_cast + (omp_target_alloc(sizeof(CubqlHit), gpu_id)); const auto& surface_bvh_indices = tree_to_surface_bvh_indices_.at(tree); for (const auto surface_bvh_index : surface_bvh_indices) { @@ -379,11 +382,11 @@ CuBQLRayTracer::ray_fire(TreeID tree, const cuBQL::vec3d* d_normals = surface_bvh.d_normals; const MeshID* d_primitive_refs = surface_bvh.d_primitive_refs; - CuBQLRayHit surface_hit; + CubqlHit surface_hit; surface_hit.distance = closest_distance + cubql_ray_hit_tolerance(closest_distance); omp_target_memcpy(d_surface_hit, &surface_hit, - sizeof(CuBQLRayHit), + sizeof(CubqlHit), 0, 0, gpu_id, @@ -396,9 +399,17 @@ CuBQLRayTracer::ray_fire(TreeID tree, #pragma omp target device(gpu_id) \ is_device_ptr(d_vertices, d_indices, d_normals, d_primitive_refs, d_exclude_primitives, d_surface_hit) { - cuBQL::ray3d ray(ray_origin, ray_direction, 0.0, d_surface_hit->distance); - - auto intersect_prim = [=, &ray] + cuBQL::ray3d ray; + ray.origin = ray_origin; + ray.direction = ray_direction; + ray.tMin = 0.0; + ray.tMax = d_surface_hit->distance; + + // lambda for primitive intersection + // - culls excluded primitives + // - culls backfacing or frontfacing hits based on hitOrientation + // - calls plucker intersection + auto xdg_plucker_intersect_prim = [=, &ray] (uint32_t prim_id) -> double { const MeshID primitive_ref = d_primitive_refs[prim_id]; @@ -441,23 +452,27 @@ CuBQLRayTracer::ray_fire(TreeID tree, return ray.tMax; }; - cuBQL::shrinkingRayQuery::forEachPrim(intersect_prim, bvh, ray); + // Traversal template from cuBQL + cuBQL::shrinkingRayQuery::forEachPrim(xdg_plucker_intersect_prim, bvh, ray); } omp_target_memcpy(&surface_hit, d_surface_hit, - sizeof(CuBQLRayHit), + sizeof(CubqlHit), 0, 0, host_id, gpu_id); - + + // Temporary cuBQL-only closest-hit reduction while each surface BVH is queried + // independently. Embree/GPRT delegate this to a single scene/TLAS traversal; + // once cuBQL does the same, this tolerance/tie-break logic should go away. const double hit_tolerance = cubql_ray_hit_tolerance(closest_distance); const bool closer_hit = surface_hit.distance < closest_distance - hit_tolerance; const bool tied_hit = cuBQL::abst(surface_hit.distance - closest_distance) <= hit_tolerance; const bool lower_surface_tie = closest_surface == ID_NONE || surface_bvh.surface < closest_surface; - + if (surface_hit.primitive != ID_NONE && (closer_hit || (tied_hit && lower_surface_tie))) { closest_distance = surface_hit.distance; closest_surface = surface_bvh.surface; From cf4346fd16673559d88bcc6ae7959b9e7998df9c Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 7 May 2026 13:50:01 +0100 Subject: [PATCH 19/44] Cleaned up the plucker intersection function for cuBQL --- CMakeLists.txt | 8 +++ include/xdg/geometry/dp_math.h | 22 +++++- src/cuBQL/ray_tracer.cpp | 124 +++------------------------------ 3 files changed, 38 insertions(+), 116 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 046dcaa5..ae99db2a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -242,6 +242,14 @@ if (XDG_ENABLE_CUBQL) list(APPEND xdg_sources 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 + PROPERTIES COMPILE_DEFINITIONS XDG_DP_MATH_CUBQL +) + endif() if (XDG_ENABLE_LIBMESH) 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/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index bc9680ac..ea0e2216 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -1,5 +1,6 @@ #include "xdg/cuBQL/ray_tracer.h" #include "xdg/error.h" +#include "xdg/geometry/plucker.h" #include #include "cuBQL/builder/omp.h" @@ -10,122 +11,12 @@ namespace xdg { -struct CuBQLPluckerIntersectionResult { - bool hit {false}; - double t {0.0}; -}; - static inline __cubql_both double cubql_ray_hit_tolerance(double t) { constexpr double tolerance = 64.0 * 2.2204460492503131e-16; return tolerance * (1.0 + cuBQL::abst(t)); } -static inline __cubql_both bool cubql_plucker_first(cuBQL::vec3d a, - cuBQL::vec3d b) -{ - if (a.x < b.x) return true; - if (a.x > b.x) return false; - - if (a.y < b.y) return true; - if (a.y > b.y) return false; - - return a.z < b.z; -} - -static inline __cubql_both double cubql_plucker_edge_test(cuBQL::vec3d vertex_a, - cuBQL::vec3d vertex_b, - cuBQL::vec3d ray, - cuBQL::vec3d ray_normal) -{ - double pip; - if (cubql_plucker_first(vertex_a, vertex_b)) { - const cuBQL::vec3d edge = vertex_b - vertex_a; - const cuBQL::vec3d edge_normal = cross(edge, vertex_a); - pip = dot(ray, edge_normal) + dot(ray_normal, edge); - } else { - const cuBQL::vec3d edge = vertex_a - vertex_b; - const cuBQL::vec3d edge_normal = cross(edge, vertex_b); - pip = dot(ray, edge_normal) + dot(ray_normal, edge); - pip = -pip; - } - - constexpr double dbl_zero_tol = 20.0 * 2.2204460492503131e-16; - if (cuBQL::abst(pip) < dbl_zero_tol) { - pip = 0.0; - } - return pip; -} - -static inline __cubql_both CuBQLPluckerIntersectionResult -cubql_plucker_ray_tri_intersect(cuBQL::vec3d vertices[3], - cuBQL::vec3d origin, - cuBQL::vec3d direction, - double t_max, - double t_min) -{ - const cuBQL::vec3d ray_a = direction; - const cuBQL::vec3d ray_b = cross(direction, origin); - - const double plucker_coord_0 = - cubql_plucker_edge_test(vertices[0], vertices[1], ray_a, ray_b); - - const double plucker_coord_1 = - cubql_plucker_edge_test(vertices[1], vertices[2], ray_a, ray_b); - - if ((0.0 < plucker_coord_0 && 0.0 > plucker_coord_1) || - (0.0 > plucker_coord_0 && 0.0 < plucker_coord_1)) { - return {}; - } - - const double plucker_coord_2 = - cubql_plucker_edge_test(vertices[2], vertices[0], ray_a, ray_b); - - if ((0.0 < plucker_coord_1 && 0.0 > plucker_coord_2) || - (0.0 > plucker_coord_1 && 0.0 < plucker_coord_2) || - (0.0 < plucker_coord_0 && 0.0 > plucker_coord_2) || - (0.0 > plucker_coord_0 && 0.0 < plucker_coord_2)) { - return {}; - } - - if (plucker_coord_0 == 0.0 && - plucker_coord_1 == 0.0 && - plucker_coord_2 == 0.0) { - return {}; - } - - const double inverse_sum = - 1.0 / (plucker_coord_0 + plucker_coord_1 + plucker_coord_2); - - const cuBQL::vec3d intersection = - plucker_coord_0 * inverse_sum * vertices[2] + - plucker_coord_1 * inverse_sum * vertices[0] + - plucker_coord_2 * inverse_sum * vertices[1]; - - int idx = 0; - double max_abs_dir = 0.0; - for (int i = 0; i < 3; ++i) { - if (cuBQL::abst(direction[i]) > max_abs_dir) { - idx = i; - max_abs_dir = cuBQL::abst(direction[i]); - } - } - - double dist_out = (intersection[idx] - origin[idx]) / direction[idx]; - - const double u = plucker_coord_2 * inverse_sum; - const double v = plucker_coord_0 * inverse_sum; - if (u < 0.0 || v < 0.0 || (u + v) > 1.0) { - return {}; - } - - if (dist_out < t_min || dist_out > t_max) { - return {}; - } - - return {true, dist_out}; -} - CuBQLRayTracer::CuBQLRayTracer() = default; CuBQLRayTracer::~CuBQLRayTracer() @@ -438,11 +329,14 @@ CuBQLRayTracer::ray_fire(TreeID tree, return ray.tMax; } - auto intersection = cubql_plucker_ray_tri_intersect(vertices, - ray.origin, - ray.direction, - ray.tMax, - ray.tMin); + // Reuse same plucker intersection function applied to other ray tracers + auto intersection = plucker_ray_tri_intersect(vertices, + ray.origin, + ray.direction, + ray.tMax, + ray.tMin, + false, + 0); if (intersection.hit) { d_surface_hit->distance = intersection.t; d_surface_hit->primitive = primitive_ref; From b8a5c704c6d77cd0e033d85b2e590f1c558245e9 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 7 May 2026 13:54:27 +0100 Subject: [PATCH 20/44] Removed constexpr EARLY_EXIT as nvlink couldn't resolve it --- include/xdg/geometry/plucker.h | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) 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}; } From 7b7937c3b5d7cf09f8a69cea843f086f2b50c9a7 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 14 May 2026 14:01:19 +0100 Subject: [PATCH 21/44] Fixes for rebase from main --- CMakeLists.txt | 6 +----- vendor/GPRT | 2 +- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ae99db2a..697a20e4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -355,7 +355,6 @@ if (XDG_ENABLE_GPRT) target_link_options(xdg PRIVATE -Wl,--unresolved-symbols=ignore-in-shared-libs) endif() -target_link_libraries(xdg PRIVATE fmt::fmt) if (XDG_ENABLE_CUBQL) target_compile_definitions(xdg PUBLIC XDG_ENABLE_CUBQL) target_link_libraries(xdg $) @@ -368,10 +367,7 @@ if (XDG_ENABLE_CUBQL) # whether install RPATH should remain environment-module based or be opt-in. endif() -target_link_libraries(xdg fmt::fmt) - -#TODO: refactor to only be required when tools are enabled -target_link_libraries(xdg indicators::indicators) +target_link_libraries(xdg PRIVATE fmt::fmt) # ========================== # Link ray tracing libraries diff --git a/vendor/GPRT b/vendor/GPRT index f1e95e41..7b55a1ae 160000 --- a/vendor/GPRT +++ b/vendor/GPRT @@ -1 +1 @@ -Subproject commit f1e95e4188cde591547d6b4a33a70bf2afaeec59 +Subproject commit 7b55a1ae432b226be18592f6f2e7f8980c8df9ce From 615ec1577e98aef714843e2cf8c8392b7615ee33 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 14 May 2026 16:35:19 +0100 Subject: [PATCH 22/44] Updated CMakeLists to make proper use of PRIVATE after rebasing recent changes from main --- CMakeLists.txt | 2 +- tests/CMakeLists.txt | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 697a20e4..484e71f4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -357,7 +357,7 @@ endif() if (XDG_ENABLE_CUBQL) target_compile_definitions(xdg PUBLIC XDG_ENABLE_CUBQL) - target_link_libraries(xdg $) + 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. diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 21691467..bdb1f7e8 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -54,6 +54,9 @@ foreach(test ${TEST_NAMES}) if (XDG_ENABLE_MOAB) target_link_libraries(${test} PRIVATE MOAB) endif() + if (XDG_ENABLE_CUBQL) + target_link_libraries(${test} PRIVATE $) + endif() set_target_properties(${test} PROPERTIES BUILD_RPATH "$") catch_discover_tests(${test} From 9fa439ba7d81a1464a039e0a2d909b4af3f80a4e Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 15 May 2026 13:42:54 +0100 Subject: [PATCH 23/44] Fix runtime check for available rt libraries --- tests/util.h | 6 ------ 1 file changed, 6 deletions(-) diff --git a/tests/util.h b/tests/util.h index 86657d9d..074ff899 100644 --- a/tests/util.h +++ b/tests/util.h @@ -88,12 +88,6 @@ inline void check_ray_tracer_supported(xdg::RTLibrary rt) { if (!ray_tracer_available(rt)) { SKIP(fmt::format("{} backend unavailable; skipping.", xdg::RT_LIB_TO_STR.at(rt))); } - #endif - #ifndef XDG_ENABLE_CUBQL - if(!ray_tracer_available(xdg::RTLibrary::CUBQL)) { - SKIP(fmt::format("{} backend unavailable; skipping.", xdg::RT_LIB_TO_STR.at(rt))); - } - #endif } inline void check_mesh_library_supported(xdg::MeshLibrary mesh) { From 6c244a2256560784b23d10fbe5496ee14e932d54 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 15 May 2026 14:21:24 +0100 Subject: [PATCH 24/44] Reducing the number of queries in cross check for now since cubql pathway a little slow --- tests/test_ray_tracer_cross_check.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_ray_tracer_cross_check.cpp b/tests/test_ray_tracer_cross_check.cpp index c75b0273..3cde9f73 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(); } From ac205984b2459995a28f105a006968481fc1cc55 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 15 May 2026 14:24:03 +0100 Subject: [PATCH 25/44] Reworked normal/sense handling so normals are calculated and sense handling is performed at intersection time Also had to change the create_global_*_trees stubs for cuBQL to return a warning rather than fatal_error Cross check verified to produce consistent results across all three ray tracers --- include/xdg/cuBQL/ray_tracer.h | 4 ++- src/cuBQL/ray_tracer.cpp | 49 ++++++++++++++-------------------- 2 files changed, 23 insertions(+), 30 deletions(-) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index f0a59b04..528e16ce 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -73,11 +73,12 @@ class CuBQLRayTracer : public RayTracer { private: struct CuBQLSurfaceBVH { MeshID surface {ID_NONE}; + MeshID forward_parent {ID_NONE}; + MeshID reverse_parent {ID_NONE}; cuBQL::bvh3d bvh; cuBQL::vec3d* d_vertices {nullptr}; cuBQL::vec3i* d_indices {nullptr}; - cuBQL::vec3d* d_normals {nullptr}; MeshID* d_primitive_refs {nullptr}; uint32_t num_vertices {0}; @@ -93,6 +94,7 @@ class CuBQLRayTracer : public RayTracer { std::vector surface_bvhs_; std::unordered_map> tree_to_surface_bvh_indices_; + std::unordered_map surface_tree_to_volume_; }; } // namespace xdg diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index ea0e2216..22406943 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -34,9 +34,6 @@ CuBQLRayTracer::~CuBQLRayTracer() if (surface_bvh.d_indices) { omp_target_free(surface_bvh.d_indices, surface_bvh.gpu_id); } - if (surface_bvh.d_normals) { - omp_target_free(surface_bvh.d_normals, surface_bvh.gpu_id); - } if (surface_bvh.d_primitive_refs) { omp_target_free(surface_bvh.d_primitive_refs, surface_bvh.gpu_id); } @@ -78,6 +75,7 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man SurfaceTreeID tree = next_surface_tree_id(); surface_trees_.push_back(tree); + surface_tree_to_volume_[tree] = volume_id; auto& surface_bvh_indices = tree_to_surface_bvh_indices_[tree]; auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); @@ -104,23 +102,8 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man std::vector h_primitive_refs; h_primitive_refs.reserve(num_faces); - std::vector h_normals; - h_normals.reserve(num_faces); - - auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); for (const auto& face : mesh_manager->get_surface_faces(surf)) { h_primitive_refs.push_back(face); - - // TODO: This pre-orients normals for the current volume tree. If cuBQL - // starts reusing a single surface BVH/BLAS across both parent volumes, - // normals must be flipped at intersection time based on the queried volume. - auto normal = mesh_manager->face_normal(face); - if (volume_id == reverse_parent) { - normal = -normal; - } else if (volume_id != forward_parent) { - fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); - } - h_normals.emplace_back(normal.x, normal.y, normal.z); } // copy vertices and indices to device @@ -134,11 +117,6 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man omp_target_memcpy(d_indices, h_indices.data(), h_indices.size() * sizeof(cuBQL::vec3i), 0, 0, gpu_id, host_id); - auto* d_normals = static_cast - (omp_target_alloc(h_normals.size() * sizeof(cuBQL::vec3d), gpu_id)); - omp_target_memcpy(d_normals, h_normals.data(), - h_normals.size() * sizeof(cuBQL::vec3d), 0, 0, gpu_id, host_id); - auto* d_primitive_refs = static_cast (omp_target_alloc(h_primitive_refs.size() * sizeof(MeshID), gpu_id)); omp_target_memcpy(d_primitive_refs, h_primitive_refs.data(), @@ -175,12 +153,19 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man omp_target_free(d_aabbs, gpu_id); + auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); + CuBQLSurfaceBVH surface_bvh; + surface_bvh.surface = surf; + surface_bvh.forward_parent = forward_parent; + surface_bvh.reverse_parent = reverse_parent; + if (volume_id != forward_parent && volume_id != reverse_parent) { + fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); + } surface_bvh.bvh = bvh; surface_bvh.d_vertices = d_vertices; surface_bvh.d_indices = d_indices; - surface_bvh.d_normals = d_normals; surface_bvh.d_primitive_refs = d_primitive_refs; surface_bvh.num_vertices = h_vertices.size(); surface_bvh.num_faces = num_faces; @@ -202,12 +187,12 @@ CuBQLRayTracer::create_element_tree(const std::shared_ptr&, void CuBQLRayTracer::create_global_surface_tree() { - fatal_error("Global surface trees not currently supported with cuBQL ray tracer"); + warning("Global surface trees not currently supported with cuBQL ray tracer"); } void CuBQLRayTracer::create_global_element_tree() { - fatal_error("Global element trees not currently supported with cuBQL ray tracer"); + warning("Global element trees not currently supported with cuBQL ray tracer"); } MeshID CuBQLRayTracer::find_element(const Position&) const @@ -263,6 +248,7 @@ CuBQLRayTracer::ray_fire(TreeID tree, auto* d_surface_hit = static_cast (omp_target_alloc(sizeof(CubqlHit), gpu_id)); + const MeshID query_volume = surface_tree_to_volume_.at(tree); const auto& surface_bvh_indices = tree_to_surface_bvh_indices_.at(tree); for (const auto surface_bvh_index : surface_bvh_indices) { const auto& surface_bvh = surface_bvhs_.at(surface_bvh_index); @@ -270,8 +256,8 @@ CuBQLRayTracer::ray_fire(TreeID tree, cuBQL::bvh3d bvh = surface_bvh.bvh; const cuBQL::vec3d* d_vertices = surface_bvh.d_vertices; const cuBQL::vec3i* d_indices = surface_bvh.d_indices; - const cuBQL::vec3d* d_normals = surface_bvh.d_normals; const MeshID* d_primitive_refs = surface_bvh.d_primitive_refs; + const bool reverse_sense = query_volume == surface_bvh.reverse_parent; CubqlHit surface_hit; surface_hit.distance = closest_distance + cubql_ray_hit_tolerance(closest_distance); @@ -288,7 +274,7 @@ CuBQLRayTracer::ray_fire(TreeID tree, const cuBQL::vec3d ray_direction(direction.x, direction.y, direction.z); #pragma omp target device(gpu_id) \ - is_device_ptr(d_vertices, d_indices, d_normals, d_primitive_refs, d_exclude_primitives, d_surface_hit) + is_device_ptr(d_vertices, d_indices, d_primitive_refs, d_exclude_primitives, d_surface_hit) { cuBQL::ray3d ray; ray.origin = ray_origin; @@ -318,7 +304,12 @@ CuBQLRayTracer::ray_fire(TreeID tree, d_vertices[index.z] }; - const double normal_dot_direction = dot(d_normals[prim_id], ray.direction); + cuBQL::vec3d normal = cuBQL::cross(vertices[1] - vertices[0], vertices[2] - vertices[0]); + if (reverse_sense) { + normal = -normal; + } + + const double normal_dot_direction = dot(normal, ray.direction); bool culled = false; if (orientation == static_cast(HitOrientation::EXITING)) { culled = normal_dot_direction < 0.0; From d130c41e9c0c403f792a120591069372dc502e71 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 15 May 2026 14:39:48 +0100 Subject: [PATCH 26/44] Added an omp_target_device_probe and put it into the same header as the vk_device probe --- ...ulkan_probe.h => available_device_probe.h} | 48 ++++++++++++++++++- src/cuBQL/ray_tracer.cpp | 8 +++- src/gprt/ray_tracer.cpp | 2 +- tests/util.h | 4 +- 4 files changed, 57 insertions(+), 5 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/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 22406943..4e3e7add 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -1,6 +1,7 @@ #include "xdg/cuBQL/ray_tracer.h" #include "xdg/error.h" #include "xdg/geometry/plucker.h" +#include "xdg/available_device_probe.h" #include #include "cuBQL/builder/omp.h" @@ -17,7 +18,12 @@ static inline __cubql_both double cubql_ray_hit_tolerance(double t) return tolerance * (1.0 + cuBQL::abst(t)); } -CuBQLRayTracer::CuBQLRayTracer() = default; +CuBQLRayTracer::CuBQLRayTracer() +{ + if (!system_has_omp_target_device()) { + fatal_error("No OpenMP target capable device found; cannot initialize cuBQL ray tracer."); + } +} CuBQLRayTracer::~CuBQLRayTracer() { diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index 15d21a7b..9db035ec 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -1,6 +1,6 @@ #include "xdg/gprt/ray_tracer.h" #include "gprt/gprt.h" -#include "xdg/gprt/vulkan_probe.h" +#include "xdg/available_device_probe.h" namespace xdg { diff --git a/tests/util.h b/tests/util.h index 074ff899..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 { @@ -52,7 +52,7 @@ inline bool ray_tracer_available(xdg::RTLibrary rt) { case xdg::RTLibrary::CUBQL: #ifdef XDG_ENABLE_CUBQL - return true; + return system_has_omp_target_device(); #else return false; #endif From d0d7366276428a546a0495cdb172f44923c52b33 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 15 May 2026 16:38:16 +0100 Subject: [PATCH 27/44] Introduced a cubql_backend Context object to store host_id and device_id --- include/xdg/cuBQL/cuBQL_backend.h | 30 +++++++++++++++++++++++++++ include/xdg/cuBQL/ray_tracer.h | 17 ++++++---------- src/cuBQL/ray_tracer.cpp | 34 +++++++++++++++++-------------- 3 files changed, 55 insertions(+), 26 deletions(-) create mode 100644 include/xdg/cuBQL/cuBQL_backend.h 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/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index 528e16ce..c737ded6 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -8,19 +8,11 @@ #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" -// 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" - - - namespace xdg { class CuBQLRayTracer : public RayTracer { @@ -71,7 +63,9 @@ class CuBQLRayTracer : public RayTracer { double& dist) const override; private: - struct CuBQLSurfaceBVH { + cubql::Context context_; + + struct CuBQLSurfaceBLAS { MeshID surface {ID_NONE}; MeshID forward_parent {ID_NONE}; MeshID reverse_parent {ID_NONE}; @@ -85,6 +79,7 @@ class CuBQLRayTracer : public RayTracer { uint32_t num_faces {0}; int gpu_id {0}; }; + struct CubqlHit { double distance {INFTY}; @@ -92,7 +87,7 @@ class CuBQLRayTracer : public RayTracer { MeshID primitive {ID_NONE}; }; - std::vector surface_bvhs_; + std::vector surface_bvhs_; std::unordered_map> tree_to_surface_bvh_indices_; std::unordered_map surface_tree_to_volume_; }; diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 4e3e7add..96705041 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -23,6 +23,9 @@ 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() @@ -64,8 +67,8 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man MeshID volume_id) { - int gpu_id = 0; // TODO - how to manage GPU IDs in a multi-GPU system? - int host_id = omp_get_initial_device(); + const auto& context = context_; + const int gpu_id = context.gpuID; // cuBQL::box3d *d_boxes = nullptr; // box3d is an alias for box_t @@ -116,17 +119,17 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man auto* d_vertices = static_cast (omp_target_alloc(h_vertices.size() * sizeof(cuBQL::vec3d), gpu_id)); omp_target_memcpy(d_vertices, h_vertices.data(), - h_vertices.size() * sizeof(cuBQL::vec3d), 0, 0, gpu_id, host_id); + h_vertices.size() * sizeof(cuBQL::vec3d), 0, 0, gpu_id, context.hostID); auto* d_indices = static_cast (omp_target_alloc(h_indices.size() * sizeof(cuBQL::vec3i), gpu_id)); omp_target_memcpy(d_indices, h_indices.data(), - h_indices.size() * sizeof(cuBQL::vec3i), 0, 0, gpu_id, host_id); + h_indices.size() * sizeof(cuBQL::vec3i), 0, 0, gpu_id, context.hostID); auto* d_primitive_refs = static_cast (omp_target_alloc(h_primitive_refs.size() * sizeof(MeshID), gpu_id)); omp_target_memcpy(d_primitive_refs, h_primitive_refs.data(), - h_primitive_refs.size() * sizeof(MeshID), 0, 0, gpu_id, host_id); + h_primitive_refs.size() * sizeof(MeshID), 0, 0, gpu_id, context.hostID); // Create device storage for triangle AABBs to be computed in parallel on GPU auto* d_aabbs = static_cast @@ -160,15 +163,14 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man omp_target_free(d_aabbs, gpu_id); auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); + if (volume_id != forward_parent && volume_id != reverse_parent) { + fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); + } - CuBQLSurfaceBVH surface_bvh; - + CuBQLSurfaceBLAS surface_bvh; surface_bvh.surface = surf; surface_bvh.forward_parent = forward_parent; surface_bvh.reverse_parent = reverse_parent; - if (volume_id != forward_parent && volume_id != reverse_parent) { - fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); - } surface_bvh.bvh = bvh; surface_bvh.d_vertices = d_vertices; surface_bvh.d_indices = d_indices; @@ -180,6 +182,8 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man surface_bvh_indices.push_back(surface_bvhs_.size()); surface_bvhs_.push_back(surface_bvh); } + + return tree; } @@ -230,8 +234,8 @@ CuBQLRayTracer::ray_fire(TreeID tree, HitOrientation hitOrientation, std::vector* const exclude_primitives) { - int gpu_id = 0; // TODO - how to manage GPU IDs in a multi-GPU system? - int host_id = omp_get_initial_device(); + const auto& context = context_; + const int gpu_id = context.gpuID; MeshID* d_exclude_primitives = nullptr; int exclude_count = 0; @@ -245,7 +249,7 @@ CuBQLRayTracer::ray_fire(TreeID tree, 0, 0, gpu_id, - host_id); + context.hostID); } double closest_distance = tmax; @@ -273,7 +277,7 @@ CuBQLRayTracer::ray_fire(TreeID tree, 0, 0, gpu_id, - host_id); + context.hostID); const int orientation = static_cast(hitOrientation); const cuBQL::vec3d ray_origin(origin.x, origin.y, origin.z); @@ -352,7 +356,7 @@ CuBQLRayTracer::ray_fire(TreeID tree, sizeof(CubqlHit), 0, 0, - host_id, + context.hostID, gpu_id); // Temporary cuBQL-only closest-hit reduction while each surface BVH is queried From 5a8faf327516dd7db2608310ec15ddb213c2ba37 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 15 May 2026 18:03:37 +0100 Subject: [PATCH 28/44] Started adding some more abstraction towards two-level traversal scheme --- include/xdg/cuBQL/ray_tracer.h | 21 +----- src/cuBQL/ray_tracer.cpp | 128 ++++++++++++++++++++++----------- 2 files changed, 90 insertions(+), 59 deletions(-) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index c737ded6..e9cc2308 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -65,30 +65,15 @@ class CuBQLRayTracer : public RayTracer { private: cubql::Context context_; - struct CuBQLSurfaceBLAS { - MeshID surface {ID_NONE}; - MeshID forward_parent {ID_NONE}; - MeshID reverse_parent {ID_NONE}; - cuBQL::bvh3d bvh; - - cuBQL::vec3d* d_vertices {nullptr}; - cuBQL::vec3i* d_indices {nullptr}; - MeshID* d_primitive_refs {nullptr}; - - uint32_t num_vertices {0}; - uint32_t num_faces {0}; - int gpu_id {0}; - }; - - struct CubqlHit { double distance {INFTY}; MeshID surface {ID_NONE}; MeshID primitive {ID_NONE}; }; - std::vector surface_bvhs_; - std::unordered_map> tree_to_surface_bvh_indices_; + std::vector surface_meshes_; + std::vector surface_blases_; + std::unordered_map> tree_to_surface_blas_indices_; std::unordered_map surface_tree_to_volume_; }; diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 96705041..81d69181 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -30,21 +30,30 @@ CuBQLRayTracer::CuBQLRayTracer() CuBQLRayTracer::~CuBQLRayTracer() { - for (auto& surface_bvh : surface_bvhs_) { - cuBQL::omp::Context context(surface_bvh.gpu_id); + for (auto& surface_blas : surface_blases_) { + cuBQL::omp::Context context(surface_blas.gpu_id); - if (surface_bvh.bvh.nodes || surface_bvh.bvh.primIDs) { - cuBQL::omp::freeBVH(surface_bvh.bvh, &context); + if (surface_blas.bvh.nodes || surface_blas.bvh.primIDs) { + cuBQL::omp::freeBVH(surface_blas.bvh, &context); } - if (surface_bvh.d_vertices) { - omp_target_free(surface_bvh.d_vertices, surface_bvh.gpu_id); + if (surface_blas.d_meshDDs) { + omp_target_free(surface_blas.d_meshDDs, surface_blas.gpu_id); } - if (surface_bvh.d_indices) { - omp_target_free(surface_bvh.d_indices, surface_bvh.gpu_id); + if (surface_blas.d_primRefs) { + omp_target_free(surface_blas.d_primRefs, surface_blas.gpu_id); } - if (surface_bvh.d_primitive_refs) { - omp_target_free(surface_bvh.d_primitive_refs, surface_bvh.gpu_id); + } + + for (auto& surface_mesh : surface_meshes_) { + if (surface_mesh.d_vertices) { + omp_target_free(surface_mesh.d_vertices, surface_mesh.gpu_id); + } + if (surface_mesh.d_indices) { + omp_target_free(surface_mesh.d_indices, surface_mesh.gpu_id); + } + if (surface_mesh.d_primitive_refs) { + omp_target_free(surface_mesh.d_primitive_refs, surface_mesh.gpu_id); } } } @@ -85,7 +94,7 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man SurfaceTreeID tree = next_surface_tree_id(); surface_trees_.push_back(tree); surface_tree_to_volume_[tree] = volume_id; - auto& surface_bvh_indices = tree_to_surface_bvh_indices_[tree]; + auto& surface_blas_indices = tree_to_surface_blas_indices_[tree]; auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); for (const auto &surf : volume_surfaces) { @@ -167,20 +176,55 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); } - CuBQLSurfaceBLAS surface_bvh; - surface_bvh.surface = surf; - surface_bvh.forward_parent = forward_parent; - surface_bvh.reverse_parent = reverse_parent; - surface_bvh.bvh = bvh; - surface_bvh.d_vertices = d_vertices; - surface_bvh.d_indices = d_indices; - surface_bvh.d_primitive_refs = d_primitive_refs; - surface_bvh.num_vertices = h_vertices.size(); - surface_bvh.num_faces = num_faces; - surface_bvh.gpu_id = gpu_id; - - surface_bvh_indices.push_back(surface_bvhs_.size()); - surface_bvhs_.push_back(surface_bvh); + CuBQLSurfaceMesh surface_mesh; + surface_mesh.dd.surface_id = surf; + surface_mesh.dd.forward_parent = forward_parent; + surface_mesh.dd.reverse_parent = reverse_parent; + 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 = gpu_id; + + std::vector h_meshDDs {surface_mesh.get_device_data()}; + auto* d_meshDDs = static_cast + (omp_target_alloc(h_meshDDs.size() * sizeof(CuBQLSurfaceMesh::DD), gpu_id)); + omp_target_memcpy(d_meshDDs, + h_meshDDs.data(), + h_meshDDs.size() * sizeof(CuBQLSurfaceMesh::DD), + 0, + 0, + gpu_id, + context.hostID); + + std::vector h_primRefs; + h_primRefs.reserve(num_faces); + for (uint32_t primID = 0; primID < num_faces; ++primID) { + h_primRefs.push_back({0, static_cast(primID)}); + } + + auto* d_primRefs = static_cast + (omp_target_alloc(h_primRefs.size() * sizeof(CuBQLSurfaceBLAS::PrimRef), gpu_id)); + omp_target_memcpy(d_primRefs, + h_primRefs.data(), + h_primRefs.size() * sizeof(CuBQLSurfaceBLAS::PrimRef), + 0, + 0, + gpu_id, + context.hostID); + + CuBQLSurfaceBLAS surface_blas; + surface_blas.bvh = bvh; + surface_blas.d_meshDDs = d_meshDDs; + surface_blas.d_primRefs = d_primRefs; + surface_blas.num_meshes = h_meshDDs.size(); + surface_blas.num_prims = h_primRefs.size(); + surface_blas.gpu_id = gpu_id; + + surface_blas_indices.push_back(surface_blases_.size()); + surface_meshes_.push_back(surface_mesh); + surface_blases_.push_back(surface_blas); } @@ -259,15 +303,13 @@ CuBQLRayTracer::ray_fire(TreeID tree, (omp_target_alloc(sizeof(CubqlHit), gpu_id)); const MeshID query_volume = surface_tree_to_volume_.at(tree); - const auto& surface_bvh_indices = tree_to_surface_bvh_indices_.at(tree); - for (const auto surface_bvh_index : surface_bvh_indices) { - const auto& surface_bvh = surface_bvhs_.at(surface_bvh_index); + const auto& surface_blas_indices = tree_to_surface_blas_indices_.at(tree); + for (const auto surface_blas_index : surface_blas_indices) { + const auto& surface_blas = surface_blases_.at(surface_blas_index); - cuBQL::bvh3d bvh = surface_bvh.bvh; - const cuBQL::vec3d* d_vertices = surface_bvh.d_vertices; - const cuBQL::vec3i* d_indices = surface_bvh.d_indices; - const MeshID* d_primitive_refs = surface_bvh.d_primitive_refs; - const bool reverse_sense = query_volume == surface_bvh.reverse_parent; + cuBQL::bvh3d bvh = surface_blas.bvh; + const CuBQLSurfaceMesh::DD* d_meshDDs = surface_blas.d_meshDDs; + const CuBQLSurfaceBLAS::PrimRef* d_primRefs = surface_blas.d_primRefs; CubqlHit surface_hit; surface_hit.distance = closest_distance + cubql_ray_hit_tolerance(closest_distance); @@ -284,7 +326,7 @@ CuBQLRayTracer::ray_fire(TreeID tree, const cuBQL::vec3d ray_direction(direction.x, direction.y, direction.z); #pragma omp target device(gpu_id) \ - is_device_ptr(d_vertices, d_indices, d_primitive_refs, d_exclude_primitives, d_surface_hit) + is_device_ptr(d_meshDDs, d_primRefs, d_exclude_primitives, d_surface_hit) { cuBQL::ray3d ray; ray.origin = ray_origin; @@ -299,7 +341,9 @@ CuBQLRayTracer::ray_fire(TreeID tree, auto xdg_plucker_intersect_prim = [=, &ray] (uint32_t prim_id) -> double { - const MeshID primitive_ref = d_primitive_refs[prim_id]; + const CuBQLSurfaceBLAS::PrimRef prim_ref = d_primRefs[prim_id]; + const CuBQLSurfaceMesh::DD mesh = d_meshDDs[prim_ref.meshID]; + const MeshID primitive_ref = mesh.primitive_refs[prim_ref.primID]; for (int i = 0; i < exclude_count; ++i) { if (d_exclude_primitives[i] == primitive_ref) { @@ -307,14 +351,15 @@ CuBQLRayTracer::ray_fire(TreeID tree, } } - const cuBQL::vec3i index = d_indices[prim_id]; + const cuBQL::vec3i index = mesh.indices[prim_ref.primID]; cuBQL::vec3d vertices[3] = { - d_vertices[index.x], - d_vertices[index.y], - d_vertices[index.z] + 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]); + const bool reverse_sense = query_volume == mesh.reverse_parent; if (reverse_sense) { normal = -normal; } @@ -340,6 +385,7 @@ CuBQLRayTracer::ray_fire(TreeID tree, 0); if (intersection.hit) { d_surface_hit->distance = intersection.t; + d_surface_hit->surface = mesh.surface_id; d_surface_hit->primitive = primitive_ref; ray.tMax = intersection.t; } @@ -366,11 +412,11 @@ CuBQLRayTracer::ray_fire(TreeID tree, const bool closer_hit = surface_hit.distance < closest_distance - hit_tolerance; const bool tied_hit = cuBQL::abst(surface_hit.distance - closest_distance) <= hit_tolerance; const bool lower_surface_tie = closest_surface == ID_NONE || - surface_bvh.surface < closest_surface; + surface_hit.surface < closest_surface; if (surface_hit.primitive != ID_NONE && (closer_hit || (tied_hit && lower_surface_tie))) { closest_distance = surface_hit.distance; - closest_surface = surface_bvh.surface; + closest_surface = surface_hit.surface; closest_primitive = surface_hit.primitive; } } From 2b218aa7fc0a29a2ac8be7d777a611f80ae69fb9 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 15 May 2026 18:05:25 +0100 Subject: [PATCH 29/44] Forgot to commit the new file which holds the cuBQL structs --- include/xdg/cuBQL/triangles.h | 92 +++++++++++++++++++++++++++++++++++ 1 file changed, 92 insertions(+) create mode 100644 include/xdg/cuBQL/triangles.h diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h new file mode 100644 index 00000000..df4cc270 --- /dev/null +++ b/include/xdg/cuBQL/triangles.h @@ -0,0 +1,92 @@ +#ifndef _XDG_CUBQL_TRIANGLES_H +#define _XDG_CUBQL_TRIANGLES_H + +#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 { + +// Triangle mesh data structure for storing triangle vertex/index data. +// Including a device-data struct (DD) +struct CuBQLSurfaceMesh { + struct DD { + // Topological metadata + MeshID surface_id {ID_NONE}; + MeshID forward_parent {ID_NONE}; + MeshID reverse_parent {ID_NONE}; + + // Geometric data + const cuBQL::vec3d* vertices {nullptr}; + const cuBQL::vec3i* indices {nullptr}; + const MeshID* primitive_refs {nullptr}; + }; + + DD dd; + + // Device buffers for triangle data + cuBQL::vec3d* d_vertices {nullptr}; + cuBQL::vec3i* d_indices {nullptr}; + MeshID* d_primitive_refs {nullptr}; + + // Metadata of the surface mesh TODO - do we really need these stored here? + 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 + { + DD device_data; + device_data.surface_id = dd.surface_id; + device_data.forward_parent = dd.forward_parent; + device_data.reverse_parent = dd.reverse_parent; + device_data.vertices = d_vertices; + device_data.indices = d_indices; + device_data.primitive_refs = d_primitive_refs; + return device_data; + } +}; + +struct CuBQLSurfaceBLAS { + struct PrimRef { + int meshID {0}; + int primID {0}; + }; + + struct DD { + CuBQLSurfaceMesh::DD* meshes {nullptr}; + PrimRef* primRefs {nullptr}; + cuBQL::bvh3d bvh; + }; + + cuBQL::bvh3d bvh; + CuBQLSurfaceMesh::DD* d_meshDDs {nullptr}; + PrimRef* d_primRefs {nullptr}; + + uint32_t num_meshes {0}; + uint32_t num_prims {0}; + int gpu_id {0}; + + DD get_device_data() const + { + DD device_data; + device_data.meshes = d_meshDDs; + device_data.primRefs = d_primRefs; + device_data.bvh = bvh; + return device_data; + } +}; + +} // namespace xdg + +#endif // include guard From afe30249d2da81a008d65b0d534ca14f23a8a14f Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 18 May 2026 17:11:22 +0100 Subject: [PATCH 30/44] Two level traversal scheme in place and working - New structs - Corrected lambdas - New mappings --- CMakeLists.txt | 1 + include/xdg/cuBQL/ray_tracer.h | 4 +- include/xdg/cuBQL/triangles.h | 98 ++++++--- src/cuBQL/ray_tracer.cpp | 364 ++++++++++++++++----------------- src/cuBQL/triangles.cpp | 55 +++++ 5 files changed, 304 insertions(+), 218 deletions(-) create mode 100644 src/cuBQL/triangles.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 484e71f4..abb172d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -240,6 +240,7 @@ endif() if (XDG_ENABLE_CUBQL) list(APPEND xdg_sources +src/cuBQL/triangles.cpp src/cuBQL/ray_tracer.cpp ) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index e9cc2308..960a5783 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -71,9 +71,7 @@ class CuBQLRayTracer : public RayTracer { MeshID primitive {ID_NONE}; }; - std::vector surface_meshes_; - std::vector surface_blases_; - std::unordered_map> tree_to_surface_blas_indices_; + std::unordered_map tree_to_volume_tlas_; std::unordered_map surface_tree_to_volume_; }; diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h index df4cc270..51399fed 100644 --- a/include/xdg/cuBQL/triangles.h +++ b/include/xdg/cuBQL/triangles.h @@ -2,6 +2,7 @@ #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. @@ -16,14 +17,15 @@ namespace xdg { -// Triangle mesh data structure for storing triangle vertex/index data. -// Including a device-data struct (DD) +/* + 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}; - MeshID forward_parent {ID_NONE}; - MeshID reverse_parent {ID_NONE}; // Geometric data const cuBQL::vec3d* vertices {nullptr}; @@ -31,14 +33,14 @@ struct CuBQLSurfaceMesh { const MeshID* primitive_refs {nullptr}; }; - DD dd; + // 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}; - // Metadata of the surface mesh TODO - do we really need these stored here? uint32_t num_vertices {0}; uint32_t num_triangles {0}; int gpu_id {0}; @@ -46,45 +48,81 @@ struct CuBQLSurfaceMesh { // Accessor for Device Data struct, which is passed to cuBQL BVH traversal/intersection functions DD get_device_data() const { - DD device_data; - device_data.surface_id = dd.surface_id; - device_data.forward_parent = dd.forward_parent; - device_data.reverse_parent = dd.reverse_parent; - device_data.vertices = d_vertices; - device_data.indices = d_indices; - device_data.primitive_refs = d_primitive_refs; - return device_data; + return { + surface_id, + d_vertices, + d_indices, + d_primitive_refs + }; } + + void release(); }; +/* + Owns a cuBQL bottom-level acceleration structure over surface triangles. + The nested DD type is the compact device-data view used during traversal. + + TODO - DPRT also supports grouping multiple triangle meshes into one BLAS. + xdg currently uses a simpler one topological surface mesh per BLAS layout + but maybe it would be useful for us to also support multiple surfaces being + tied to the same BLAS? +*/ struct CuBQLSurfaceBLAS { - struct PrimRef { - int meshID {0}; - int primID {0}; + struct DD { + CuBQLSurfaceMesh::DD mesh; // Mesh data device handle + cuBQL::bvh3d bvh; // B:AS device handle + }; + + cuBQL::bvh3d 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 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}; }; struct DD { - CuBQLSurfaceMesh::DD* meshes {nullptr}; - PrimRef* primRefs {nullptr}; - cuBQL::bvh3d bvh; + const SurfaceInstanceDD* surface_instances {nullptr}; + cuBQL::bvh3d bvh; // TLAS device handle }; - cuBQL::bvh3d bvh; - CuBQLSurfaceMesh::DD* d_meshDDs {nullptr}; - PrimRef* d_primRefs {nullptr}; + cuBQL::bvh3d bvh; // TLAS host handle + SurfaceInstanceDD* d_surface_instances {nullptr}; + std::vector surface_blases; - uint32_t num_meshes {0}; - uint32_t num_prims {0}; + uint32_t num_surface_instances {0}; int gpu_id {0}; DD get_device_data() const { - DD device_data; - device_data.meshes = d_meshDDs; - device_data.primRefs = d_primRefs; - device_data.bvh = bvh; - return device_data; + return {d_surface_instances, bvh}; } + + void release(); }; } // namespace xdg diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 81d69181..8063efe2 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -12,12 +12,6 @@ namespace xdg { -static inline __cubql_both double cubql_ray_hit_tolerance(double t) -{ - constexpr double tolerance = 64.0 * 2.2204460492503131e-16; - return tolerance * (1.0 + cuBQL::abst(t)); -} - CuBQLRayTracer::CuBQLRayTracer() { if (!system_has_omp_target_device()) { @@ -30,31 +24,8 @@ CuBQLRayTracer::CuBQLRayTracer() CuBQLRayTracer::~CuBQLRayTracer() { - for (auto& surface_blas : surface_blases_) { - cuBQL::omp::Context context(surface_blas.gpu_id); - - if (surface_blas.bvh.nodes || surface_blas.bvh.primIDs) { - cuBQL::omp::freeBVH(surface_blas.bvh, &context); - } - - if (surface_blas.d_meshDDs) { - omp_target_free(surface_blas.d_meshDDs, surface_blas.gpu_id); - } - if (surface_blas.d_primRefs) { - omp_target_free(surface_blas.d_primRefs, surface_blas.gpu_id); - } - } - - for (auto& surface_mesh : surface_meshes_) { - if (surface_mesh.d_vertices) { - omp_target_free(surface_mesh.d_vertices, surface_mesh.gpu_id); - } - if (surface_mesh.d_indices) { - omp_target_free(surface_mesh.d_indices, surface_mesh.gpu_id); - } - if (surface_mesh.d_primitive_refs) { - omp_target_free(surface_mesh.d_primitive_refs, surface_mesh.gpu_id); - } + for (auto& [tree, tlas] : tree_to_volume_tlas_) { + tlas.release(); } } @@ -94,8 +65,13 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man SurfaceTreeID tree = next_surface_tree_id(); surface_trees_.push_back(tree); surface_tree_to_volume_[tree] = volume_id; - auto& surface_blas_indices = tree_to_surface_blas_indices_[tree]; auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); + std::vector h_tlas_boxes; + std::vector h_surface_instances; + std::vector h_surface_blases; + h_tlas_boxes.reserve(volume_surfaces.size()); + h_surface_instances.reserve(volume_surfaces.size()); + h_surface_blases.reserve(volume_surfaces.size()); for (const auto &surf : volume_surfaces) { auto num_faces = mesh_manager->num_surface_faces(surf); @@ -117,12 +93,7 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man } // Get storage for primitive refs so a hit can be mapped back to a mesh face. - std::vector h_primitive_refs; - h_primitive_refs.reserve(num_faces); - - for (const auto& face : mesh_manager->get_surface_faces(surf)) { - h_primitive_refs.push_back(face); - } + std::vector h_primitive_refs = mesh_manager->get_surface_faces(surf); // copy vertices and indices to device auto* d_vertices = static_cast @@ -177,9 +148,7 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man } CuBQLSurfaceMesh surface_mesh; - surface_mesh.dd.surface_id = surf; - surface_mesh.dd.forward_parent = forward_parent; - surface_mesh.dd.reverse_parent = reverse_parent; + surface_mesh.surface_id = surf; surface_mesh.d_vertices = d_vertices; surface_mesh.d_indices = d_indices; surface_mesh.d_primitive_refs = d_primitive_refs; @@ -187,46 +156,68 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man surface_mesh.num_triangles = num_faces; surface_mesh.gpu_id = gpu_id; - std::vector h_meshDDs {surface_mesh.get_device_data()}; - auto* d_meshDDs = static_cast - (omp_target_alloc(h_meshDDs.size() * sizeof(CuBQLSurfaceMesh::DD), gpu_id)); - omp_target_memcpy(d_meshDDs, - h_meshDDs.data(), - h_meshDDs.size() * sizeof(CuBQLSurfaceMesh::DD), - 0, - 0, - gpu_id, - context.hostID); - - std::vector h_primRefs; - h_primRefs.reserve(num_faces); - for (uint32_t primID = 0; primID < num_faces; ++primID) { - h_primRefs.push_back({0, static_cast(primID)}); - } - - auto* d_primRefs = static_cast - (omp_target_alloc(h_primRefs.size() * sizeof(CuBQLSurfaceBLAS::PrimRef), gpu_id)); - omp_target_memcpy(d_primRefs, - h_primRefs.data(), - h_primRefs.size() * sizeof(CuBQLSurfaceBLAS::PrimRef), - 0, - 0, - gpu_id, - context.hostID); - CuBQLSurfaceBLAS surface_blas; surface_blas.bvh = bvh; - surface_blas.d_meshDDs = d_meshDDs; - surface_blas.d_primRefs = d_primRefs; - surface_blas.num_meshes = h_meshDDs.size(); - surface_blas.num_prims = h_primRefs.size(); + surface_blas.mesh = surface_mesh; + surface_blas.num_prims = num_faces; surface_blas.gpu_id = gpu_id; - surface_blas_indices.push_back(surface_blases_.size()); - surface_meshes_.push_back(surface_mesh); - surface_blases_.push_back(surface_blas); + // Store BLAS bounding boxes to build TLAS + const auto surface_bounding_box = mesh_manager->surface_bounding_box(surf); + cuBQL::box3d surface_bounds; + surface_bounds.lower = cuBQL::vec3d(surface_bounding_box.min_x, + surface_bounding_box.min_y, + surface_bounding_box.min_z); + surface_bounds.upper = cuBQL::vec3d(surface_bounding_box.max_x, + surface_bounding_box.max_y, + surface_bounding_box.max_z); + + h_tlas_boxes.push_back(surface_bounds); + h_surface_instances.push_back({surface_blas.get_device_data()}); + h_surface_blases.push_back(surface_blas); } + 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::box3d), gpu_id)); + omp_target_memcpy(d_tlas_boxes, + h_tlas_boxes.data(), + h_tlas_boxes.size() * sizeof(cuBQL::box3d), + 0, + 0, + gpu_id, + context.hostID); + + auto* d_surface_instances = static_cast + (omp_target_alloc(h_surface_instances.size() * sizeof(CuBQLVolumeTLAS::SurfaceInstanceDD), gpu_id)); + omp_target_memcpy(d_surface_instances, + h_surface_instances.data(), + h_surface_instances.size() * sizeof(CuBQLVolumeTLAS::SurfaceInstanceDD), + 0, + 0, + gpu_id, + context.hostID); + + cuBQL::BuildConfig tlasBuildParams(1); + tlasBuildParams.maxAllowedLeafSize = 1; + + CuBQLVolumeTLAS volume_tlas; + volume_tlas.num_surface_instances = static_cast(h_surface_instances.size()); + volume_tlas.gpu_id = gpu_id; + volume_tlas.d_surface_instances = d_surface_instances; + volume_tlas.surface_blases = std::move(h_surface_blases); + cuBQL::build_omp_target(volume_tlas.bvh, + d_tlas_boxes, + volume_tlas.num_surface_instances, + tlasBuildParams, + gpu_id); + + omp_target_free(d_tlas_boxes, gpu_id); + + tree_to_volume_tlas_.emplace(tree, std::move(volume_tlas)); return tree; } @@ -302,123 +293,126 @@ CuBQLRayTracer::ray_fire(TreeID tree, auto* d_surface_hit = static_cast (omp_target_alloc(sizeof(CubqlHit), gpu_id)); - const MeshID query_volume = surface_tree_to_volume_.at(tree); - const auto& surface_blas_indices = tree_to_surface_blas_indices_.at(tree); - for (const auto surface_blas_index : surface_blas_indices) { - const auto& surface_blas = surface_blases_.at(surface_blas_index); - - cuBQL::bvh3d bvh = surface_blas.bvh; - const CuBQLSurfaceMesh::DD* d_meshDDs = surface_blas.d_meshDDs; - const CuBQLSurfaceBLAS::PrimRef* d_primRefs = surface_blas.d_primRefs; - - CubqlHit surface_hit; - surface_hit.distance = closest_distance + cubql_ray_hit_tolerance(closest_distance); - omp_target_memcpy(d_surface_hit, - &surface_hit, - sizeof(CubqlHit), - 0, - 0, - gpu_id, - context.hostID); - - const int orientation = static_cast(hitOrientation); - const cuBQL::vec3d ray_origin(origin.x, origin.y, origin.z); - const cuBQL::vec3d ray_direction(direction.x, direction.y, direction.z); - - #pragma omp target device(gpu_id) \ - is_device_ptr(d_meshDDs, d_primRefs, d_exclude_primitives, d_surface_hit) + const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); + const auto volume_tlas_dd = volume_tlas.get_device_data(); + + CubqlHit surface_hit; + surface_hit.distance = tmax; + omp_target_memcpy(d_surface_hit, + &surface_hit, + sizeof(CubqlHit), + 0, + 0, + gpu_id, + context.hostID); + + const int orientation = static_cast(hitOrientation); + const cuBQL::vec3d ray_origin(origin.x, origin.y, origin.z); + const cuBQL::vec3d ray_direction(direction.x, direction.y, direction.z); + + #pragma omp target device(gpu_id) \ + is_device_ptr(d_exclude_primitives, d_surface_hit) + { + cuBQL::ray3d world_ray; + world_ray.origin = ray_origin; + world_ray.direction = ray_direction; + world_ray.tMin = 0.0; + world_ray.tMax = d_surface_hit->distance; + + CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; + + // lambda to go from tlas->blas + auto enter_blas = [=, &surface_instance, &world_ray] + (cuBQL::ray3d& out_ray, cuBQL::bvh3d& out_bvh, int instance_id) { - cuBQL::ray3d ray; - ray.origin = ray_origin; - ray.direction = ray_direction; - ray.tMin = 0.0; - ray.tMax = d_surface_hit->distance; - - // lambda for primitive intersection - // - culls excluded primitives - // - culls backfacing or frontfacing hits based on hitOrientation - // - calls plucker intersection - auto xdg_plucker_intersect_prim = [=, &ray] - (uint32_t prim_id) -> double - { - const CuBQLSurfaceBLAS::PrimRef prim_ref = d_primRefs[prim_id]; - const CuBQLSurfaceMesh::DD mesh = d_meshDDs[prim_ref.meshID]; - const MeshID primitive_ref = mesh.primitive_refs[prim_ref.primID]; - - for (int i = 0; i < exclude_count; ++i) { - if (d_exclude_primitives[i] == primitive_ref) { - return ray.tMax; - } - } - - const cuBQL::vec3i index = mesh.indices[prim_ref.primID]; - 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]); - const bool reverse_sense = query_volume == mesh.reverse_parent; - if (reverse_sense) { - normal = -normal; - } - - const double normal_dot_direction = dot(normal, ray.direction); - bool culled = false; - if (orientation == static_cast(HitOrientation::EXITING)) { - culled = normal_dot_direction < 0.0; - } else if (orientation == static_cast(HitOrientation::ENTERING)) { - culled = normal_dot_direction >= 0.0; - } - if (culled) { - return ray.tMax; - } + surface_instance = volume_tlas_dd.surface_instances[instance_id]; + out_ray = world_ray; // No ray transformation because everything is in world coords with xdg + out_bvh = surface_instance.surface_blas.bvh; // Get the BLAS for the surface instance we hit in the TLAS + }; + + + // lambda for primitive intersection + // - culls excluded primitives + // - culls backfacing or frontfacing hits based on hitOrientation + // - calls plucker intersection + auto xdg_plucker_intersect_prim = [=, &world_ray, &surface_instance] + (uint32_t prim_id) -> double + { + const CuBQLSurfaceMesh::DD mesh = surface_instance.surface_blas.mesh; + const MeshID primitive_ref = mesh.primitive_refs[prim_id]; - // Reuse same plucker intersection function applied to other ray tracers - auto intersection = plucker_ray_tri_intersect(vertices, - ray.origin, - ray.direction, - ray.tMax, - ray.tMin, - false, - 0); - if (intersection.hit) { - d_surface_hit->distance = intersection.t; - d_surface_hit->surface = mesh.surface_id; - d_surface_hit->primitive = primitive_ref; - ray.tMax = intersection.t; + for (int i = 0; i < exclude_count; ++i) { + if (d_exclude_primitives[i] == primitive_ref) { + return world_ray.tMax; } + } - return 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] }; - // Traversal template from cuBQL - cuBQL::shrinkingRayQuery::forEachPrim(xdg_plucker_intersect_prim, bvh, ray); - } - - omp_target_memcpy(&surface_hit, - d_surface_hit, - sizeof(CubqlHit), - 0, - 0, - context.hostID, - gpu_id); - - // Temporary cuBQL-only closest-hit reduction while each surface BVH is queried - // independently. Embree/GPRT delegate this to a single scene/TLAS traversal; - // once cuBQL does the same, this tolerance/tie-break logic should go away. - const double hit_tolerance = cubql_ray_hit_tolerance(closest_distance); - const bool closer_hit = surface_hit.distance < closest_distance - hit_tolerance; - const bool tied_hit = cuBQL::abst(surface_hit.distance - closest_distance) <= hit_tolerance; - const bool lower_surface_tie = closest_surface == ID_NONE || - surface_hit.surface < closest_surface; + 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, world_ray.direction); + bool culled = false; + if (orientation == static_cast(HitOrientation::EXITING)) { + culled = normal_dot_direction < 0.0; + } else if (orientation == static_cast(HitOrientation::ENTERING)) { + culled = normal_dot_direction >= 0.0; + } + if (culled) { + return world_ray.tMax; + } + + // Reuse same plucker intersection function applied to other ray tracers + auto intersection = plucker_ray_tri_intersect(vertices, + world_ray.origin, + world_ray.direction, + world_ray.tMax, + world_ray.tMin, + false, + 0); + if (intersection.hit) { + d_surface_hit->distance = intersection.t; + d_surface_hit->surface = mesh.surface_id; + d_surface_hit->primitive = primitive_ref; + world_ray.tMax = intersection.t; + } + + return world_ray.tMax; + }; + + // We must provide the lambda to run on exit but it is empty for our use case + // since there is no need for cleanup or transforming the ray back to world space + auto leave_blas = []() -> void {}; + + // Two level Traversal template from cuBQL + cuBQL::shrinkingRayQuery::twoLevel::forEachPrim(enter_blas, + leave_blas, + xdg_plucker_intersect_prim, + volume_tlas_dd.bvh, + world_ray); - if (surface_hit.primitive != ID_NONE && (closer_hit || (tied_hit && lower_surface_tie))) { - closest_distance = surface_hit.distance; - closest_surface = surface_hit.surface; - closest_primitive = surface_hit.primitive; - } + } + + omp_target_memcpy(&surface_hit, + d_surface_hit, + sizeof(CubqlHit), + 0, + 0, + context.hostID, + gpu_id); + + if (surface_hit.primitive != ID_NONE) { + closest_distance = surface_hit.distance; + closest_surface = surface_hit.surface; + closest_primitive = surface_hit.primitive; } omp_target_free(d_surface_hit, gpu_id); diff --git a/src/cuBQL/triangles.cpp b/src/cuBQL/triangles.cpp new file mode 100644 index 00000000..7bf202cb --- /dev/null +++ b/src/cuBQL/triangles.cpp @@ -0,0 +1,55 @@ +#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; + } + for (auto& surface_blas : surface_blases) { + surface_blas.release(); + } +} + +} // namespace xdg From 089ba15948b63d91be81fd9182f4cc3cc12d8a4c Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 18 May 2026 17:21:13 +0100 Subject: [PATCH 31/44] Fixed sense handling for surface second visit --- src/cuBQL/ray_tracer.cpp | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 8063efe2..5f5fe271 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -143,9 +143,6 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man omp_target_free(d_aabbs, gpu_id); auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); - if (volume_id != forward_parent && volume_id != reverse_parent) { - fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); - } CuBQLSurfaceMesh surface_mesh; surface_mesh.surface_id = surf; @@ -172,8 +169,20 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man surface_bounding_box.max_y, surface_bounding_box.max_z); + 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_blas.get_device_data()}); + h_surface_instances.push_back(surface_instance); h_surface_blases.push_back(surface_blas); } From b18eeb9ffc872ac9df516eb47f108cfeb2f13162 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 20 May 2026 12:23:55 +0100 Subject: [PATCH 32/44] Refactor to build only a single CuBQLSurfaceMesh and CuBQLSurfaceBLAS per topological surface --- include/xdg/cuBQL/ray_tracer.h | 5 +- include/xdg/cuBQL/triangles.h | 1 - src/cuBQL/ray_tracer.cpp | 210 ++++++++++++++++++--------------- src/cuBQL/triangles.cpp | 3 - vendor/GPRT | 2 +- 5 files changed, 123 insertions(+), 98 deletions(-) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index 960a5783..cb373623 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -63,6 +63,9 @@ class CuBQLRayTracer : public RayTracer { double& dist) const override; private: + CuBQLSurfaceBLAS + register_surface(const std::shared_ptr& mesh_manager, MeshID surface_id); + cubql::Context context_; struct CubqlHit { @@ -72,7 +75,7 @@ class CuBQLRayTracer : public RayTracer { }; std::unordered_map tree_to_volume_tlas_; - std::unordered_map surface_tree_to_volume_; + std::unordered_map surface_to_blas_map_; }; } // namespace xdg diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h index 51399fed..dca264ad 100644 --- a/include/xdg/cuBQL/triangles.h +++ b/include/xdg/cuBQL/triangles.h @@ -112,7 +112,6 @@ struct CuBQLVolumeTLAS { cuBQL::bvh3d bvh; // TLAS host handle SurfaceInstanceDD* d_surface_instances {nullptr}; - std::vector surface_blases; uint32_t num_surface_instances {0}; int gpu_id {0}; diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 5f5fe271..941ab752 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -27,6 +27,10 @@ CuBQLRayTracer::~CuBQLRayTracer() for (auto& [tree, tlas] : tree_to_volume_tlas_) { tlas.release(); } + + for (auto& [surface, blas] : surface_to_blas_map_) { + blas.release(); + } } void CuBQLRayTracer::init() @@ -42,9 +46,112 @@ CuBQLRayTracer::register_volume(const std::shared_ptr& mesh_manager return {surface_tree, element_tree}; } +CuBQLSurfaceBLAS +CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manager, + MeshID surface_id) +{ + const auto& context = context_; + const int gpu_id = context.gpuID; + 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), gpu_id)); + omp_target_memcpy(d_vertices, + h_vertices.data(), + h_vertices.size() * sizeof(cuBQL::vec3d), + 0, + 0, + gpu_id, + context.hostID); + + auto* d_indices = static_cast + (omp_target_alloc(h_indices.size() * sizeof(cuBQL::vec3i), gpu_id)); + omp_target_memcpy(d_indices, + h_indices.data(), + h_indices.size() * sizeof(cuBQL::vec3i), + 0, + 0, + gpu_id, + context.hostID); + + auto* d_primitive_refs = static_cast + (omp_target_alloc(h_primitive_refs.size() * sizeof(MeshID), gpu_id)); + omp_target_memcpy(d_primitive_refs, + h_primitive_refs.data(), + h_primitive_refs.size() * sizeof(MeshID), + 0, + 0, + gpu_id, + context.hostID); + + auto* d_aabbs = static_cast + (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3d), gpu_id)); + const auto num_primitives = static_cast(h_indices.size()); + + // TODO - Abstract this out into its own bounding_box creation function + #pragma omp target device(gpu_id) is_device_ptr(d_vertices, d_indices, d_aabbs) + #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); + + d_aabbs[primID] = aabb; + } + + cuBQL::BuildConfig blasBuildParams; + // TODO - Try setting leaf params to 1 to see what it does + // Check what default is for CUDA + cuBQL::bvh3d bvh; + cuBQL::build_omp_target(bvh, d_aabbs, num_faces, blasBuildParams, gpu_id); + + omp_target_free(d_aabbs, gpu_id); + + 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 = gpu_id; + + CuBQLSurfaceBLAS surface_blas; + surface_blas.bvh = bvh; + surface_blas.mesh = surface_mesh; + surface_blas.num_prims = num_faces; + surface_blas.gpu_id = gpu_id; + + return surface_blas; +} + TreeID CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_manager, - MeshID volume_id) + MeshID volume_id) { const auto& context = context_; @@ -55,110 +162,30 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man // int num_boxes = 0; // cuBQL::bvh3d bvh; // bvh3d is an alias for BinaryBVH - // // bvh_t is an alias for BinaryBVH, so bvh_t is also BinaryBVH + // bvh_t is an alias for BinaryBVH, so bvh_t is also BinaryBVH + + // Looks like the cuBQL::gpuBuilder is only pulled in when cubql is built with CUDA support + + // It looks like spatial median is the only supported method for omp builder - cuBQL::BuildConfig buildParams; // It looks like spatial median is the only supported method for omp builder - // // Looks like the cuBQL::gpuBuilder is only pulled in when cubql is built with CUDA support - // cuBQL::build_omp_target(bvh, d_boxes, num_boxes, buildParams, gpu_id); // build bvh on GPU 0 using OpenMP target offloading SurfaceTreeID tree = next_surface_tree_id(); surface_trees_.push_back(tree); - surface_tree_to_volume_[tree] = volume_id; auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); std::vector h_tlas_boxes; std::vector h_surface_instances; - std::vector h_surface_blases; h_tlas_boxes.reserve(volume_surfaces.size()); h_surface_instances.reserve(volume_surfaces.size()); - h_surface_blases.reserve(volume_surfaces.size()); for (const auto &surf : volume_surfaces) { - auto num_faces = mesh_manager->num_surface_faces(surf); - auto vertices = mesh_manager->get_surface_vertices(surf); - auto indices = mesh_manager->get_surface_connectivity(surf); - - // Get host side storage for vertices using cubql friendly types - std::vector h_vertices; - h_vertices.reserve(vertices.size()); - for (const auto& vertex : vertices) { - h_vertices.emplace_back(vertex.x, vertex.y, vertex.z); - } - - // Get host side storage for indices using cubql friendly types - 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]); + if (!surface_to_blas_map_.count(surf)) { + surface_to_blas_map_[surf] = register_surface(mesh_manager, surf); } - // Get storage for primitive refs so a hit can be mapped back to a mesh face. - std::vector h_primitive_refs = mesh_manager->get_surface_faces(surf); - - // copy vertices and indices to device - auto* d_vertices = static_cast - (omp_target_alloc(h_vertices.size() * sizeof(cuBQL::vec3d), gpu_id)); - omp_target_memcpy(d_vertices, h_vertices.data(), - h_vertices.size() * sizeof(cuBQL::vec3d), 0, 0, gpu_id, context.hostID); - - auto* d_indices = static_cast - (omp_target_alloc(h_indices.size() * sizeof(cuBQL::vec3i), gpu_id)); - omp_target_memcpy(d_indices, h_indices.data(), - h_indices.size() * sizeof(cuBQL::vec3i), 0, 0, gpu_id, context.hostID); - - auto* d_primitive_refs = static_cast - (omp_target_alloc(h_primitive_refs.size() * sizeof(MeshID), gpu_id)); - omp_target_memcpy(d_primitive_refs, h_primitive_refs.data(), - h_primitive_refs.size() * sizeof(MeshID), 0, 0, gpu_id, context.hostID); - - // Create device storage for triangle AABBs to be computed in parallel on GPU - auto* d_aabbs = static_cast - (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3d), gpu_id)); - const auto num_primitives = static_cast(h_indices.size()); - - // Create AABBs for each triangle in parallel on GPU - // TODO - Should this be its own function? - // TODO - How can this be extended for tets and other element types? - #pragma omp target device(gpu_id) is_device_ptr(d_vertices, d_indices, d_aabbs) - #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); - - d_aabbs[primID] = aabb; - } // This is our AABB population kernel written with OpenMP target offloading pragmas to run in parallel on the GPU - - // Construct the bvh on the gpu using the AABBs with openmp pathway - cuBQL::bvh3d bvh; // bvh3d is an alias for BinaryBVH - cuBQL::build_omp_target(bvh, d_aabbs, num_faces, buildParams, gpu_id); - - omp_target_free(d_aabbs, gpu_id); - + CuBQLSurfaceBLAS& surface_blas = surface_to_blas_map_.at(surf); auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); - CuBQLSurfaceMesh surface_mesh; - surface_mesh.surface_id = surf; - 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 = gpu_id; - - CuBQLSurfaceBLAS surface_blas; - surface_blas.bvh = bvh; - surface_blas.mesh = surface_mesh; - surface_blas.num_prims = num_faces; - surface_blas.gpu_id = gpu_id; - // Store BLAS bounding boxes to build TLAS const auto surface_bounding_box = mesh_manager->surface_bounding_box(surf); cuBQL::box3d surface_bounds; @@ -183,7 +210,6 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man h_tlas_boxes.push_back(surface_bounds); h_surface_instances.push_back(surface_instance); - h_surface_blases.push_back(surface_blas); } if (h_surface_instances.empty()) { @@ -210,14 +236,14 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man gpu_id, context.hostID); - cuBQL::BuildConfig tlasBuildParams(1); + cuBQL::BuildConfig tlasBuildParams; + tlasBuildParams.makeLeafThreshold = 1; tlasBuildParams.maxAllowedLeafSize = 1; CuBQLVolumeTLAS volume_tlas; volume_tlas.num_surface_instances = static_cast(h_surface_instances.size()); volume_tlas.gpu_id = gpu_id; volume_tlas.d_surface_instances = d_surface_instances; - volume_tlas.surface_blases = std::move(h_surface_blases); cuBQL::build_omp_target(volume_tlas.bvh, d_tlas_boxes, volume_tlas.num_surface_instances, diff --git a/src/cuBQL/triangles.cpp b/src/cuBQL/triangles.cpp index 7bf202cb..73cf1fe3 100644 --- a/src/cuBQL/triangles.cpp +++ b/src/cuBQL/triangles.cpp @@ -47,9 +47,6 @@ void CuBQLVolumeTLAS::release() omp_target_free(d_surface_instances, gpu_id); d_surface_instances = nullptr; } - for (auto& surface_blas : surface_blases) { - surface_blas.release(); - } } } // namespace xdg diff --git a/vendor/GPRT b/vendor/GPRT index 7b55a1ae..f1e95e41 160000 --- a/vendor/GPRT +++ b/vendor/GPRT @@ -1 +1 @@ -Subproject commit 7b55a1ae432b226be18592f6f2e7f8980c8df9ce +Subproject commit f1e95e4188cde591547d6b4a33a70bf2afaeec59 From bd12a16b972b09401d748e694d9a8d2c034e1d90 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 22 May 2026 14:08:47 +0100 Subject: [PATCH 33/44] Added point_in_volume query --- include/xdg/cuBQL/ray_tracer.h | 6 - include/xdg/cuBQL/triangles.h | 11 +- src/cuBQL/ray_tracer.cpp | 216 ++++++++++++++++++++++++++------- tests/test_point_in_volume.cpp | 3 +- 4 files changed, 178 insertions(+), 58 deletions(-) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index cb373623..89829606 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -68,12 +68,6 @@ class CuBQLRayTracer : public RayTracer { cubql::Context context_; - struct CubqlHit { - double distance {INFTY}; - MeshID surface {ID_NONE}; - MeshID primitive {ID_NONE}; - }; - std::unordered_map tree_to_volume_tlas_; std::unordered_map surface_to_blas_map_; }; diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h index dca264ad..67a73c8e 100644 --- a/include/xdg/cuBQL/triangles.h +++ b/include/xdg/cuBQL/triangles.h @@ -60,13 +60,8 @@ struct CuBQLSurfaceMesh { }; /* - Owns a cuBQL bottom-level acceleration structure over surface triangles. + 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. - - TODO - DPRT also supports grouping multiple triangle meshes into one BLAS. - xdg currently uses a simpler one topological surface mesh per BLAS layout - but maybe it would be useful for us to also support multiple surfaces being - tied to the same BLAS? */ struct CuBQLSurfaceBLAS { struct DD { @@ -89,7 +84,7 @@ struct CuBQLSurfaceBLAS { }; /* - Owns a cuBQL top-level acceleration structure for one topological volume. + 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. */ @@ -102,7 +97,7 @@ struct CuBQLVolumeTLAS { */ struct SurfaceInstanceDD { CuBQLSurfaceBLAS::DD surface_blas; - bool reverse_sense {false}; + bool reverse_sense {false}; // value set in create_surface_tree based on parent vols }; struct DD { diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 941ab752..e584d9c4 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -1,4 +1,5 @@ #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" @@ -287,13 +288,154 @@ MeshID CuBQLRayTracer::find_element(TreeID, const Position&) const return ID_NONE; } -bool CuBQLRayTracer::point_in_volume(TreeID, - const Position&, - const Direction*, - const std::vector*) const +bool CuBQLRayTracer::point_in_volume(TreeID tree, + const Position& point, + const Direction* direction, + const std::vector* exclude_primitives) const { - fatal_error("Point-in-volume queries not currently supported with cuBQL ray tracer"); - return false; + const auto& context = context_; + const int gpu_id = context.gpuID; + const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); + const auto volume_tlas_dd = volume_tlas.get_device_data(); + + 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_hit = static_cast + (omp_target_alloc(sizeof(CuBQLSurfaceHit), gpu_id)); + + CuBQLSurfaceHit hit; + omp_target_memcpy(d_hit, + &hit, + sizeof(CuBQLSurfaceHit), + 0, + 0, + gpu_id, + context.hostID); + + // 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}; + + const cuBQL::vec3d ray_origin(point.x, point.y, point.z); + const cuBQL::vec3d ray_direction(directionUsed.x, directionUsed.y, directionUsed.z); + + #pragma omp target device(gpu_id) \ + is_device_ptr(d_exclude_primitives, d_hit) + { + cuBQL::ray3d world_ray; + world_ray.origin = ray_origin; + world_ray.direction = ray_direction; + world_ray.tMin = 0.0; + world_ray.tMax = d_hit->distance; + + CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; + + // lambda to go from tlas->blas + auto enter_blas = [=, &surface_instance, &world_ray] + (cuBQL::ray3d& out_ray, cuBQL::bvh3d& out_bvh, int instance_id) + { + surface_instance = volume_tlas_dd.surface_instances[instance_id]; + out_ray = world_ray; // No ray transformation because everything is in world coords with xdg + out_bvh = surface_instance.surface_blas.bvh; // Get the BLAS for the surface instance we hit in the TLAS + }; + + + // lambda for primitive intersection + // - culls excluded primitives + // - culls backfacing or frontfacing hits based on hitOrientation + // - calls plucker intersection + auto xdg_plucker_intersect_prim = [=, &world_ray, &surface_instance] + (uint32_t prim_id) -> double + { + 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 (d_exclude_primitives[i] == primitive_ref) { + return world_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, world_ray.direction); + if (orientation_cull(normal_dot_direction, HitOrientation::ANY)) { + return world_ray.tMax; + } + + // Reuse same plucker intersection function applied to other ray tracers + auto intersection = plucker_ray_tri_intersect(vertices, + world_ray.origin, + world_ray.direction, + world_ray.tMax, + world_ray.tMin, + false, + 0); + if (intersection.hit) { + d_hit->distance = intersection.t; + d_hit->surface = mesh.surface_id; + d_hit->primitive = primitive_ref; + d_hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; + world_ray.tMax = intersection.t; + } + + return world_ray.tMax; + }; + + // We must provide the lambda to run on exit but it is empty for our use case + // since there is no need for cleanup or transforming the ray back to world space + auto leave_blas = []() -> void {}; + + // Two level Traversal template from cuBQL + cuBQL::shrinkingRayQuery::twoLevel::forEachPrim(enter_blas, + leave_blas, + xdg_plucker_intersect_prim, + volume_tlas_dd.bvh, + world_ray); + + } + + omp_target_memcpy(&hit, + d_hit, + sizeof(CuBQLSurfaceHit), + 0, + 0, + context.hostID, + gpu_id); + + omp_target_free(d_hit, gpu_id); + + if (d_exclude_primitives) { + omp_target_free(d_exclude_primitives, gpu_id); + } + + // if the ray hit nothing the point must be outside the volume + if (hit.primitive == ID_NONE) return false; + + return hit.piv == INSIDE; } std::pair @@ -322,37 +464,35 @@ CuBQLRayTracer::ray_fire(TreeID tree, context.hostID); } - double closest_distance = tmax; - MeshID closest_surface = ID_NONE; - MeshID closest_primitive = ID_NONE; - auto* d_surface_hit = static_cast - (omp_target_alloc(sizeof(CubqlHit), gpu_id)); + auto* d_hit = static_cast + (omp_target_alloc(sizeof(CuBQLSurfaceHit), gpu_id)); const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); const auto volume_tlas_dd = volume_tlas.get_device_data(); - CubqlHit surface_hit; - surface_hit.distance = tmax; - omp_target_memcpy(d_surface_hit, - &surface_hit, - sizeof(CubqlHit), + CuBQLSurfaceHit hit; + hit.distance = tmax; + omp_target_memcpy(d_hit, + &hit, + sizeof(CuBQLSurfaceHit), 0, 0, gpu_id, context.hostID); - - const int orientation = static_cast(hitOrientation); + + // These are implicitly mapped to our openmp target region when we try to use them + const int ray_orientation = static_cast(hitOrientation); const cuBQL::vec3d ray_origin(origin.x, origin.y, origin.z); const cuBQL::vec3d ray_direction(direction.x, direction.y, direction.z); #pragma omp target device(gpu_id) \ - is_device_ptr(d_exclude_primitives, d_surface_hit) + is_device_ptr(d_exclude_primitives, d_hit) { cuBQL::ray3d world_ray; world_ray.origin = ray_origin; world_ray.direction = ray_direction; world_ray.tMin = 0.0; - world_ray.tMax = d_surface_hit->distance; + world_ray.tMax = d_hit->distance; CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; @@ -395,13 +535,8 @@ CuBQLRayTracer::ray_fire(TreeID tree, } const double normal_dot_direction = dot(normal, world_ray.direction); - bool culled = false; - if (orientation == static_cast(HitOrientation::EXITING)) { - culled = normal_dot_direction < 0.0; - } else if (orientation == static_cast(HitOrientation::ENTERING)) { - culled = normal_dot_direction >= 0.0; - } - if (culled) { + if (orientation_cull(normal_dot_direction, + static_cast(ray_orientation))) { return world_ray.tMax; } @@ -414,9 +549,10 @@ CuBQLRayTracer::ray_fire(TreeID tree, false, 0); if (intersection.hit) { - d_surface_hit->distance = intersection.t; - d_surface_hit->surface = mesh.surface_id; - d_surface_hit->primitive = primitive_ref; + d_hit->distance = intersection.t; + d_hit->surface = mesh.surface_id; + d_hit->primitive = primitive_ref; + d_hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; world_ray.tMax = intersection.t; } @@ -436,35 +572,29 @@ CuBQLRayTracer::ray_fire(TreeID tree, } - omp_target_memcpy(&surface_hit, - d_surface_hit, - sizeof(CubqlHit), + omp_target_memcpy(&hit, + d_hit, + sizeof(CuBQLSurfaceHit), 0, 0, context.hostID, gpu_id); - if (surface_hit.primitive != ID_NONE) { - closest_distance = surface_hit.distance; - closest_surface = surface_hit.surface; - closest_primitive = surface_hit.primitive; - } - - omp_target_free(d_surface_hit, gpu_id); + omp_target_free(d_hit, gpu_id); if (d_exclude_primitives) { omp_target_free(d_exclude_primitives, gpu_id); } - if (closest_surface == ID_NONE) { + if (hit.primitive == ID_NONE) { return {INFTY, ID_NONE}; } if (exclude_primitives) { - exclude_primitives->push_back(closest_primitive); + exclude_primitives->push_back(hit.primitive); } - return {closest_distance, closest_surface}; + return {hit.distance, hit.surface}; } std::pair diff --git a/tests/test_point_in_volume.cpp b/tests/test_point_in_volume.cpp index ae34e823..65874f9f 100644 --- a/tests/test_point_in_volume.cpp +++ b/tests/test_point_in_volume.cpp @@ -16,7 +16,8 @@ using namespace xdg::test; TEMPLATE_TEST_CASE("Point-in-volume on MeshMock", "[piv][mock]", Embree_Raytracer, - GPRT_Raytracer) + GPRT_Raytracer, + CuBQL_Raytracer) { constexpr auto rt_backend = TestType::value; From e7e260df5ace5d54ca2e56597d2e6cb1394a9f9d Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 22 May 2026 14:34:31 +0100 Subject: [PATCH 34/44] Extended ray tracer cross checks to also do point_in_volume for pincell --- tests/test_ray_tracer_cross_check.cpp | 55 +++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/tests/test_ray_tracer_cross_check.cpp b/tests/test_ray_tracer_cross_check.cpp index 3cde9f73..d2407edf 100644 --- a/tests/test_ray_tracer_cross_check.cpp +++ b/tests/test_ray_tracer_cross_check.cpp @@ -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 From d16f4813794e62a83186123c13f5f2447685fbae Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 22 May 2026 14:42:17 +0100 Subject: [PATCH 35/44] Added new intersection.h and intersection.cpp to contain helpers/structs for the intersection routine used with the cubql backend --- include/xdg/cuBQL/intersection.h | 63 ++++++++++++++ src/cuBQL/intersection.cpp | 145 +++++++++++++++++++++++++++++++ 2 files changed, 208 insertions(+) create mode 100644 include/xdg/cuBQL/intersection.h create mode 100644 src/cuBQL/intersection.cpp diff --git a/include/xdg/cuBQL/intersection.h b/include/xdg/cuBQL/intersection.h new file mode 100644 index 00000000..b2fc2e23 --- /dev/null +++ b/include/xdg/cuBQL/intersection.h @@ -0,0 +1,63 @@ +#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 "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}; +}; + +// 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; +} + +CuBQLSurfaceHit +intersect_surface_tree(const cubql::Context& context, + const CuBQLVolumeTLAS& volume_tlas, + const CuBQLRay& ray, + HitOrientation hit_orientation, + const std::vector* exclude_primitives); + +} // namespace xdg + +#endif // include guard diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp new file mode 100644 index 00000000..b08260c1 --- /dev/null +++ b/src/cuBQL/intersection.cpp @@ -0,0 +1,145 @@ +#include "xdg/cuBQL/intersection.h" +#include "xdg/geometry/plucker.h" + +#include + +#include "cuBQL/math/Ray.h" +#include "cuBQL/traversal/rayQueries.h" + +namespace xdg { + +CuBQLSurfaceHit +intersect_surface_tree(const cubql::Context& context, + const CuBQLVolumeTLAS& volume_tlas, + const CuBQLRay& ray, + 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)); + + CuBQLSurfaceHit surface_hit; + 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) + { + cuBQL::ray3d world_ray; + world_ray.origin = ray.origin; + world_ray.direction = ray.direction; + world_ray.tMin = ray.tmin; + world_ray.tMax = d_surface_hit->distance; + + CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; + + auto enter_blas = [=, &surface_instance, &world_ray] + (cuBQL::ray3d& out_ray, cuBQL::bvh3d& out_bvh, int instance_id) + { + surface_instance = volume_tlas_dd.surface_instances[instance_id]; + out_ray = world_ray; + out_bvh = surface_instance.surface_blas.bvh; + }; + + auto xdg_plucker_intersect_prim = [=, &world_ray, &surface_instance] + (uint32_t prim_id) -> double + { + 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 (d_exclude_primitives[i] == primitive_ref) { + return world_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, world_ray.direction); + if (orientation_cull(normal_dot_direction, + static_cast(orientation))) { + return world_ray.tMax; + } + + auto intersection = plucker_ray_tri_intersect(vertices, + world_ray.origin, + world_ray.direction, + world_ray.tMax, + world_ray.tMin, + false, + 0); + if (intersection.hit) { + d_surface_hit->distance = intersection.t; + d_surface_hit->surface = mesh.surface_id; + d_surface_hit->primitive = primitive_ref; + d_surface_hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; + world_ray.tMax = intersection.t; + } + + return world_ray.tMax; + }; + + auto leave_blas = []() -> void {}; + + cuBQL::shrinkingRayQuery::twoLevel::forEachPrim(enter_blas, + leave_blas, + xdg_plucker_intersect_prim, + volume_tlas_dd.bvh, + world_ray); + } + + 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 surface_hit; +} + +} // namespace xdg From 2794efd8abce15593f446fd48898fa40a0be122d Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 22 May 2026 16:14:55 +0100 Subject: [PATCH 36/44] Added cuBQL support to particle sim --- 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 84e2200a43cd8bdef6df6cc96824d8eb986c97b9 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 1 Jun 2026 17:08:42 +0100 Subject: [PATCH 37/44] Abstracted shared intersection logic between ray_fire and point_in_volume into its own function I am essentially trying to make this look more like our embree interface --- CMakeLists.txt | 2 + include/xdg/cuBQL/intersection.h | 17 +- src/cuBQL/intersection.cpp | 10 +- src/cuBQL/ray_tracer.cpp | 281 +++---------------------------- 4 files changed, 40 insertions(+), 270 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index abb172d5..7196b384 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -241,6 +241,7 @@ endif() if (XDG_ENABLE_CUBQL) list(APPEND xdg_sources src/cuBQL/triangles.cpp +src/cuBQL/intersection.cpp src/cuBQL/ray_tracer.cpp ) @@ -248,6 +249,7 @@ src/cuBQL/ray_tracer.cpp # 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 ) diff --git a/include/xdg/cuBQL/intersection.h b/include/xdg/cuBQL/intersection.h index b2fc2e23..09f6283c 100644 --- a/include/xdg/cuBQL/intersection.h +++ b/include/xdg/cuBQL/intersection.h @@ -18,10 +18,20 @@ namespace xdg { struct CuBQLRay { cuBQL::vec3d origin; cuBQL::vec3d direction; - double tmin {0.0}; - double tmax {INFTY}; + double tMin {0.0}; + double tMax {INFTY}; }; +/* 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}; @@ -51,10 +61,11 @@ inline bool orientation_cull(double normal_dot_direction, return false; } -CuBQLSurfaceHit +void intersect_surface_tree(const cubql::Context& context, const CuBQLVolumeTLAS& volume_tlas, const CuBQLRay& ray, + CuBQLSurfaceHit& hit, HitOrientation hit_orientation, const std::vector* exclude_primitives); diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index b08260c1..77501ef9 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -8,10 +8,11 @@ namespace xdg { -CuBQLSurfaceHit +void intersect_surface_tree(const cubql::Context& context, const CuBQLVolumeTLAS& volume_tlas, const CuBQLRay& ray, + CuBQLSurfaceHit& surface_hit, HitOrientation hit_orientation, const std::vector* exclude_primitives) { @@ -35,8 +36,7 @@ intersect_surface_tree(const cubql::Context& context, auto* d_surface_hit = static_cast (omp_target_alloc(sizeof(CuBQLSurfaceHit), gpu_id)); - CuBQLSurfaceHit surface_hit; - surface_hit.distance = ray.tmax; + surface_hit.distance = ray.tMax; omp_target_memcpy(d_surface_hit, &surface_hit, sizeof(CuBQLSurfaceHit), @@ -54,7 +54,7 @@ intersect_surface_tree(const cubql::Context& context, cuBQL::ray3d world_ray; world_ray.origin = ray.origin; world_ray.direction = ray.direction; - world_ray.tMin = ray.tmin; + world_ray.tMin = ray.tMin; world_ray.tMax = d_surface_hit->distance; CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; @@ -139,7 +139,7 @@ intersect_surface_tree(const cubql::Context& context, omp_target_free(d_exclude_primitives, gpu_id); } - return surface_hit; + return; } } // namespace xdg diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index e584d9c4..ac08340b 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -296,146 +296,26 @@ bool CuBQLRayTracer::point_in_volume(TreeID tree, const auto& context = context_; const int gpu_id = context.gpuID; const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); - const auto volume_tlas_dd = volume_tlas.get_device_data(); - - 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_hit = static_cast - (omp_target_alloc(sizeof(CuBQLSurfaceHit), gpu_id)); - - CuBQLSurfaceHit hit; - omp_target_memcpy(d_hit, - &hit, - sizeof(CuBQLSurfaceHit), - 0, - 0, - gpu_id, - context.hostID); // 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}; - const cuBQL::vec3d ray_origin(point.x, point.y, point.z); - const cuBQL::vec3d ray_direction(directionUsed.x, directionUsed.y, directionUsed.z); - #pragma omp target device(gpu_id) \ - is_device_ptr(d_exclude_primitives, d_hit) - { - cuBQL::ray3d world_ray; - world_ray.origin = ray_origin; - world_ray.direction = ray_direction; - world_ray.tMin = 0.0; - world_ray.tMax = d_hit->distance; - - CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; - - // lambda to go from tlas->blas - auto enter_blas = [=, &surface_instance, &world_ray] - (cuBQL::ray3d& out_ray, cuBQL::bvh3d& out_bvh, int instance_id) - { - surface_instance = volume_tlas_dd.surface_instances[instance_id]; - out_ray = world_ray; // No ray transformation because everything is in world coords with xdg - out_bvh = surface_instance.surface_blas.bvh; // Get the BLAS for the surface instance we hit in the TLAS - }; - - - // lambda for primitive intersection - // - culls excluded primitives - // - culls backfacing or frontfacing hits based on hitOrientation - // - calls plucker intersection - auto xdg_plucker_intersect_prim = [=, &world_ray, &surface_instance] - (uint32_t prim_id) -> double - { - 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 (d_exclude_primitives[i] == primitive_ref) { - return world_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, world_ray.direction); - if (orientation_cull(normal_dot_direction, HitOrientation::ANY)) { - return world_ray.tMax; - } - - // Reuse same plucker intersection function applied to other ray tracers - auto intersection = plucker_ray_tri_intersect(vertices, - world_ray.origin, - world_ray.direction, - world_ray.tMax, - world_ray.tMin, - false, - 0); - if (intersection.hit) { - d_hit->distance = intersection.t; - d_hit->surface = mesh.surface_id; - d_hit->primitive = primitive_ref; - d_hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; - world_ray.tMax = intersection.t; - } - - return world_ray.tMax; - }; - - // We must provide the lambda to run on exit but it is empty for our use case - // since there is no need for cleanup or transforming the ray back to world space - auto leave_blas = []() -> void {}; - - // Two level Traversal template from cuBQL - cuBQL::shrinkingRayQuery::twoLevel::forEachPrim(enter_blas, - leave_blas, - xdg_plucker_intersect_prim, - volume_tlas_dd.bvh, - world_ray); - - } + 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; - omp_target_memcpy(&hit, - d_hit, - sizeof(CuBQLSurfaceHit), - 0, - 0, - context.hostID, - gpu_id); - - omp_target_free(d_hit, gpu_id); + CuBQLSurfaceHit surface_hit; - if (d_exclude_primitives) { - omp_target_free(d_exclude_primitives, gpu_id); - } + intersect_surface_tree(context, volume_tlas, ray, surface_hit, HitOrientation::ANY, exclude_primitives); // if the ray hit nothing the point must be outside the volume - if (hit.primitive == ID_NONE) return false; + if (surface_hit.primitive == ID_NONE) return false; - return hit.piv == INSIDE; + return surface_hit.piv == INSIDE; } std::pair @@ -448,153 +328,30 @@ CuBQLRayTracer::ray_fire(TreeID tree, { const auto& context = context_; 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_hit = static_cast - (omp_target_alloc(sizeof(CuBQLSurfaceHit), gpu_id)); - const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); - const auto volume_tlas_dd = volume_tlas.get_device_data(); - CuBQLSurfaceHit hit; - hit.distance = tmax; - omp_target_memcpy(d_hit, - &hit, - sizeof(CuBQLSurfaceHit), - 0, - 0, - gpu_id, - context.hostID); - // These are implicitly mapped to our openmp target region when we try to use them const int ray_orientation = static_cast(hitOrientation); - const cuBQL::vec3d ray_origin(origin.x, origin.y, origin.z); - const cuBQL::vec3d ray_direction(direction.x, direction.y, direction.z); - - #pragma omp target device(gpu_id) \ - is_device_ptr(d_exclude_primitives, d_hit) - { - cuBQL::ray3d world_ray; - world_ray.origin = ray_origin; - world_ray.direction = ray_direction; - world_ray.tMin = 0.0; - world_ray.tMax = d_hit->distance; - CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; + 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; - // lambda to go from tlas->blas - auto enter_blas = [=, &surface_instance, &world_ray] - (cuBQL::ray3d& out_ray, cuBQL::bvh3d& out_bvh, int instance_id) - { - surface_instance = volume_tlas_dd.surface_instances[instance_id]; - out_ray = world_ray; // No ray transformation because everything is in world coords with xdg - out_bvh = surface_instance.surface_blas.bvh; // Get the BLAS for the surface instance we hit in the TLAS - }; - - - // lambda for primitive intersection - // - culls excluded primitives - // - culls backfacing or frontfacing hits based on hitOrientation - // - calls plucker intersection - auto xdg_plucker_intersect_prim = [=, &world_ray, &surface_instance] - (uint32_t prim_id) -> double - { - 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 (d_exclude_primitives[i] == primitive_ref) { - return world_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, world_ray.direction); - if (orientation_cull(normal_dot_direction, - static_cast(ray_orientation))) { - return world_ray.tMax; - } - - // Reuse same plucker intersection function applied to other ray tracers - auto intersection = plucker_ray_tri_intersect(vertices, - world_ray.origin, - world_ray.direction, - world_ray.tMax, - world_ray.tMin, - false, - 0); - if (intersection.hit) { - d_hit->distance = intersection.t; - d_hit->surface = mesh.surface_id; - d_hit->primitive = primitive_ref; - d_hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; - world_ray.tMax = intersection.t; - } - - return world_ray.tMax; - }; - - // We must provide the lambda to run on exit but it is empty for our use case - // since there is no need for cleanup or transforming the ray back to world space - auto leave_blas = []() -> void {}; - - // Two level Traversal template from cuBQL - cuBQL::shrinkingRayQuery::twoLevel::forEachPrim(enter_blas, - leave_blas, - xdg_plucker_intersect_prim, - volume_tlas_dd.bvh, - world_ray); - - } + CuBQLSurfaceHit surface_hit; - omp_target_memcpy(&hit, - d_hit, - sizeof(CuBQLSurfaceHit), - 0, - 0, - context.hostID, - gpu_id); - - omp_target_free(d_hit, gpu_id); - - if (d_exclude_primitives) { - omp_target_free(d_exclude_primitives, gpu_id); - } + intersect_surface_tree(context, volume_tlas, ray, surface_hit, hitOrientation, exclude_primitives); - if (hit.primitive == ID_NONE) { + if (surface_hit.primitive == ID_NONE) { return {INFTY, ID_NONE}; } if (exclude_primitives) { - exclude_primitives->push_back(hit.primitive); + exclude_primitives->push_back(surface_hit.primitive); } - return {hit.distance, hit.surface}; + return {surface_hit.distance, surface_hit.surface}; } std::pair From 7959a44f82c2c41ea95dff353f4d178528282b4c Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 2 Jun 2026 13:33:34 +0100 Subject: [PATCH 38/44] Added setup for dynamic-ray volume querying - Volume to query against carried on ray - Lookup table for volume mesh id to BVH uploaded to device --- include/xdg/cuBQL/ray_tracer.h | 6 ++++ include/xdg/cuBQL/triangles.h | 4 ++- src/cuBQL/ray_tracer.cpp | 64 +++++++++++++++++++++++++++++----- 3 files changed, 64 insertions(+), 10 deletions(-) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index 89829606..b04c990e 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -66,10 +66,16 @@ class CuBQLRayTracer : public RayTracer { CuBQLSurfaceBLAS register_surface(const std::shared_ptr& mesh_manager, MeshID surface_id); + 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 diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h index 67a73c8e..2b9f62cb 100644 --- a/include/xdg/cuBQL/triangles.h +++ b/include/xdg/cuBQL/triangles.h @@ -101,10 +101,12 @@ struct CuBQLVolumeTLAS { }; struct DD { + MeshID volume_id {ID_NONE}; const SurfaceInstanceDD* surface_instances {nullptr}; cuBQL::bvh3d bvh; // TLAS device handle }; + MeshID volume_id {ID_NONE}; cuBQL::bvh3d bvh; // TLAS host handle SurfaceInstanceDD* d_surface_instances {nullptr}; @@ -113,7 +115,7 @@ struct CuBQLVolumeTLAS { DD get_device_data() const { - return {d_surface_instances, bvh}; + return {volume_id, d_surface_instances, bvh}; } void release(); diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index ac08340b..b38eeeb8 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -25,6 +25,11 @@ CuBQLRayTracer::CuBQLRayTracer() 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(); } @@ -36,6 +41,32 @@ CuBQLRayTracer::~CuBQLRayTracer() void CuBQLRayTracer::init() { + upload_volume_to_tlas_table_(); + initialized_ = true; +} + +void CuBQLRayTracer::upload_volume_to_tlas_table_() +{ + const int gpu_id = context_.gpuID; + + if (d_volume_to_tlas_) { + omp_target_free(d_volume_to_tlas_, gpu_id); + 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), gpu_id)); + omp_target_memcpy(d_volume_to_tlas_, + volume_to_tlas_.data(), + volume_to_tlas_.size() * sizeof(CuBQLVolumeTLAS::DD), + 0, + 0, + gpu_id, + context_.hostID); } std::pair @@ -105,7 +136,7 @@ CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manage auto* d_aabbs = static_cast (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3d), gpu_id)); const auto num_primitives = static_cast(h_indices.size()); - + // TODO - Abstract this out into its own bounding_box creation function #pragma omp target device(gpu_id) is_device_ptr(d_vertices, d_indices, d_aabbs) #pragma omp teams distribute parallel for @@ -242,6 +273,7 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man 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 = gpu_id; volume_tlas.d_surface_instances = d_surface_instances; @@ -253,7 +285,24 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man omp_target_free(d_tlas_boxes, gpu_id); - tree_to_volume_tlas_.emplace(tree, std::move(volume_tlas)); + // 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; } @@ -294,7 +343,6 @@ bool CuBQLRayTracer::point_in_volume(TreeID tree, const std::vector* exclude_primitives) const { const auto& context = context_; - const int gpu_id = context.gpuID; const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); // Use provided direction or if Direction == nulptr use default direction @@ -310,7 +358,8 @@ bool CuBQLRayTracer::point_in_volume(TreeID tree, CuBQLSurfaceHit surface_hit; - intersect_surface_tree(context, volume_tlas, ray, surface_hit, HitOrientation::ANY, exclude_primitives); + // 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; @@ -327,12 +376,8 @@ CuBQLRayTracer::ray_fire(TreeID tree, std::vector* const exclude_primitives) { const auto& context = context_; - const int gpu_id = context.gpuID; const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); - // These are implicitly mapped to our openmp target region when we try to use them - const int ray_orientation = static_cast(hitOrientation); - CuBQLRay ray; ray.origin = cuBQL::vec3d(origin.x, origin.y, origin.z); ray.direction = cuBQL::vec3d(direction.x, direction.y, direction.z); @@ -341,7 +386,8 @@ CuBQLRayTracer::ray_fire(TreeID tree, CuBQLSurfaceHit surface_hit; - intersect_surface_tree(context, volume_tlas, ray, surface_hit, hitOrientation, exclude_primitives); + // 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}; From 84ab66ce722b73f68e84d9d7517f92594bb7f394 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 2 Jun 2026 14:45:27 +0100 Subject: [PATCH 39/44] Abstracted the ray traversal into its own inline function which - intersect_surface_tree() intersects a single ray with BVH - Can be used from within different intersection algorithms, i.e scalar and batch (when implemented) --- include/xdg/cuBQL/intersection.h | 15 +-- src/cuBQL/intersection.cpp | 180 +++++++++++++++++-------------- 2 files changed, 111 insertions(+), 84 deletions(-) diff --git a/include/xdg/cuBQL/intersection.h b/include/xdg/cuBQL/intersection.h index 09f6283c..06f31ee6 100644 --- a/include/xdg/cuBQL/intersection.h +++ b/include/xdg/cuBQL/intersection.h @@ -20,6 +20,7 @@ struct CuBQLRay { 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*/ @@ -61,13 +62,15 @@ inline bool orientation_cull(double normal_dot_direction, return false; } +// Wrapper for launching a single ray intersection query against the surface tree, with Host<->Device staging of ray and hit data void -intersect_surface_tree(const cubql::Context& context, - const CuBQLVolumeTLAS& volume_tlas, - const CuBQLRay& ray, - CuBQLSurfaceHit& hit, - HitOrientation hit_orientation, - const std::vector* exclude_primitives); +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); + } // namespace xdg diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index 77501ef9..5fbc00a8 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -8,13 +8,103 @@ 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 ray, + CuBQLSurfaceHit* hit, + int orientation, + const MeshID* exclude_primitives, + int exclude_count) +{ + cuBQL::ray3d world_ray; + world_ray.origin = ray.origin; + world_ray.direction = ray.direction; + world_ray.tMin = ray.tMin; + world_ray.tMax = hit->distance; + + CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; + + auto enter_blas = [=, &surface_instance, &world_ray] + (cuBQL::ray3d& out_ray, cuBQL::bvh3d& out_bvh, int instance_id) + { + surface_instance = volume_tlas.surface_instances[instance_id]; + out_ray = world_ray; + out_bvh = surface_instance.surface_blas.bvh; + }; + + auto intersect_prim = [=, &world_ray, &surface_instance] + (uint32_t prim_id) -> double + { + 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 world_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, world_ray.direction); + + if (orientation_cull(normal_dot_direction, + static_cast(orientation))) { + return world_ray.tMax; + } + + auto intersection = plucker_ray_tri_intersect(vertices, + world_ray.origin, + world_ray.direction, + world_ray.tMax, + world_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; + world_ray.tMax = intersection.t; + } + + return world_ray.tMax; + }; + + auto leave_blas = []() -> void {}; + + cuBQL::shrinkingRayQuery::twoLevel::forEachPrim(enter_blas, + leave_blas, + intersect_prim, + volume_tlas.bvh, + world_ray); +} +#pragma omp end declare target + +// Wrapper for launching a single ray intersection query against the surface tree, with Host<->Device staging of ray and hit data void -intersect_surface_tree(const cubql::Context& context, - const CuBQLVolumeTLAS& volume_tlas, - const CuBQLRay& ray, - CuBQLSurfaceHit& surface_hit, - HitOrientation hit_orientation, - const std::vector* exclude_primitives) +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; @@ -51,78 +141,12 @@ intersect_surface_tree(const cubql::Context& context, #pragma omp target device(gpu_id) \ is_device_ptr(d_exclude_primitives, d_surface_hit) { - cuBQL::ray3d world_ray; - world_ray.origin = ray.origin; - world_ray.direction = ray.direction; - world_ray.tMin = ray.tMin; - world_ray.tMax = d_surface_hit->distance; - - CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; - - auto enter_blas = [=, &surface_instance, &world_ray] - (cuBQL::ray3d& out_ray, cuBQL::bvh3d& out_bvh, int instance_id) - { - surface_instance = volume_tlas_dd.surface_instances[instance_id]; - out_ray = world_ray; - out_bvh = surface_instance.surface_blas.bvh; - }; - - auto xdg_plucker_intersect_prim = [=, &world_ray, &surface_instance] - (uint32_t prim_id) -> double - { - 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 (d_exclude_primitives[i] == primitive_ref) { - return world_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, world_ray.direction); - if (orientation_cull(normal_dot_direction, - static_cast(orientation))) { - return world_ray.tMax; - } - - auto intersection = plucker_ray_tri_intersect(vertices, - world_ray.origin, - world_ray.direction, - world_ray.tMax, - world_ray.tMin, - false, - 0); - if (intersection.hit) { - d_surface_hit->distance = intersection.t; - d_surface_hit->surface = mesh.surface_id; - d_surface_hit->primitive = primitive_ref; - d_surface_hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; - world_ray.tMax = intersection.t; - } - - return world_ray.tMax; - }; - - auto leave_blas = []() -> void {}; - - cuBQL::shrinkingRayQuery::twoLevel::forEachPrim(enter_blas, - leave_blas, - xdg_plucker_intersect_prim, - volume_tlas_dd.bvh, - world_ray); + intersect_surface_tree(volume_tlas_dd, + ray, + d_surface_hit, + orientation, + d_exclude_primitives, + exclude_count); } omp_target_memcpy(&surface_hit, From 9cb25171f536d7406f12751d86aceb321ceafa5d Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 3 Jun 2026 14:54:46 +0100 Subject: [PATCH 40/44] Added a basic batch ray_fire method - Expects populated device pointer of d_rays and populates d_hits --- include/xdg/cuBQL/intersection.h | 14 +++++++- include/xdg/cuBQL/ray_tracer.h | 9 +++++ src/cuBQL/intersection.cpp | 56 ++++++++++++++++++++++++++------ src/cuBQL/ray_tracer.cpp | 20 ++++++++++++ 4 files changed, 88 insertions(+), 11 deletions(-) diff --git a/include/xdg/cuBQL/intersection.h b/include/xdg/cuBQL/intersection.h index 06f31ee6..50d95141 100644 --- a/include/xdg/cuBQL/intersection.h +++ b/include/xdg/cuBQL/intersection.h @@ -7,6 +7,7 @@ #undef __CUDA_ARCH__ #endif +#include #include #include "xdg/constants.h" @@ -62,7 +63,10 @@ inline bool orientation_cull(double normal_dot_direction, return false; } -// Wrapper for launching a single ray intersection query against the surface tree, with Host<->Device staging of ray and hit data +/* +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, @@ -72,6 +76,14 @@ intersect_surface_tree_scalar(const cubql::Context& context, 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 index b04c990e..388ef29f 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -1,6 +1,7 @@ #ifndef _XDG_CUBQL_RAY_TRACING_INTERFACE_H #define _XDG_CUBQL_RAY_TRACING_INTERFACE_H +#include #include #include #include @@ -15,6 +16,9 @@ namespace xdg { +struct CuBQLRay; +struct CuBQLSurfaceHit; + class CuBQLRayTracer : public RayTracer { public: CuBQLRayTracer(); @@ -54,6 +58,11 @@ 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); + std::pair closest(TreeID tree, const Position& origin) override; diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index 5fbc00a8..eccb2f91 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -1,7 +1,8 @@ +#include + #include "xdg/cuBQL/intersection.h" #include "xdg/geometry/plucker.h" - -#include +#include "xdg/error.h" #include "cuBQL/math/Ray.h" #include "cuBQL/traversal/rayQueries.h" @@ -10,13 +11,12 @@ 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 ray, - CuBQLSurfaceHit* hit, - int orientation, - const MeshID* exclude_primitives, - int exclude_count) +static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, + CuBQLRay ray, + CuBQLSurfaceHit* hit, + int orientation, + const MeshID* exclude_primitives, + int exclude_count) { cuBQL::ray3d world_ray; world_ray.origin = ray.origin; @@ -97,7 +97,6 @@ static inline void intersect_surface_tree( } #pragma omp end declare target -// Wrapper for launching a single ray intersection query against the surface tree, with Host<->Device staging of ray and hit data void intersect_surface_tree_scalar(const cubql::Context& context, const CuBQLVolumeTLAS& volume_tlas, @@ -166,4 +165,41 @@ intersect_surface_tree_scalar(const cubql::Context& context, 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 index b38eeeb8..2ee22d6b 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -400,6 +400,26 @@ CuBQLRayTracer::ray_fire(TreeID tree, 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&) { From 4876116b710b7d5d565829e9e263e96b978862d9 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 4 Jun 2026 17:22:50 +0100 Subject: [PATCH 41/44] Working ray benchmark with cubql --- tools/ray_benchmark.cpp | 71 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index 3b1ae8d7..76f63428 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -32,6 +32,10 @@ #include "ray_benchmark.h" +#ifdef XDG_ENABLE_CUBQL +#include "xdg/cuBQL/intersection.h" +#include "xdg/cuBQL/ray_tracer.h" +#endif <<<<<<< HEAD @@ -130,6 +134,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); } @@ -327,6 +333,71 @@ int main(int argc, char** argv) trace_timer.stop(); } + 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 and trace in one step + 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 + trace_timer.start(); + rti->ray_fire_batch(d_rays, d_hits, num_rays); + trace_timer.stop(); + + omp_target_free(d_hits, gpu_id); + omp_target_free(d_rays, gpu_id); + #endif + } const double generation_time = generation_timer.elapsed(); const double trace_time = trace_timer.elapsed(); From cb7fee9989f10fd235447b2b0e1890b6a9d5055d Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 15 Jun 2026 17:21:53 +0100 Subject: [PATCH 42/44] Some minor cleanup --- src/cuBQL/ray_tracer.cpp | 76 +++++++++++++++------------------------- tools/CMakeLists.txt | 4 +++ 2 files changed, 33 insertions(+), 47 deletions(-) diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 2ee22d6b..94e87b2e 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -47,10 +47,8 @@ void CuBQLRayTracer::init() void CuBQLRayTracer::upload_volume_to_tlas_table_() { - const int gpu_id = context_.gpuID; - if (d_volume_to_tlas_) { - omp_target_free(d_volume_to_tlas_, gpu_id); + omp_target_free(d_volume_to_tlas_, context_.gpuID); d_volume_to_tlas_ = nullptr; } @@ -59,13 +57,13 @@ void CuBQLRayTracer::upload_volume_to_tlas_table_() } d_volume_to_tlas_ = static_cast - (omp_target_alloc(volume_to_tlas_.size() * sizeof(CuBQLVolumeTLAS::DD), gpu_id)); + (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, - gpu_id, + context_.gpuID, context_.hostID); } @@ -82,8 +80,6 @@ CuBQLSurfaceBLAS CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manager, MeshID surface_id) { - const auto& context = context_; - const int gpu_id = context.gpuID; 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); @@ -104,41 +100,41 @@ CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manage // 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), gpu_id)); + (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, - gpu_id, - context.hostID); + context_.gpuID, + context_.hostID); auto* d_indices = static_cast - (omp_target_alloc(h_indices.size() * sizeof(cuBQL::vec3i), gpu_id)); + (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, - gpu_id, - context.hostID); + context_.gpuID, + context_.hostID); auto* d_primitive_refs = static_cast - (omp_target_alloc(h_primitive_refs.size() * sizeof(MeshID), gpu_id)); + (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, - gpu_id, - context.hostID); + context_.gpuID, + context_.hostID); auto* d_aabbs = static_cast - (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3d), gpu_id)); + (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3d), 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(gpu_id) is_device_ptr(d_vertices, d_indices, d_aabbs) + #pragma omp target device(context_.gpuID) is_device_ptr(d_vertices, d_indices, d_aabbs) #pragma omp teams distribute parallel for for (uint32_t primID = 0; primID < num_primitives; ++primID) { cuBQL::vec3i indices = d_indices[primID]; @@ -159,9 +155,9 @@ CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manage // TODO - Try setting leaf params to 1 to see what it does // Check what default is for CUDA cuBQL::bvh3d bvh; - cuBQL::build_omp_target(bvh, d_aabbs, num_faces, blasBuildParams, gpu_id); + cuBQL::build_omp_target(bvh, d_aabbs, num_faces, blasBuildParams, context_.gpuID); - omp_target_free(d_aabbs, gpu_id); + omp_target_free(d_aabbs, context_.gpuID); CuBQLSurfaceMesh surface_mesh; surface_mesh.surface_id = surface_id; @@ -170,13 +166,13 @@ CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manage 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 = gpu_id; + 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 = gpu_id; + surface_blas.gpu_id = context_.gpuID; return surface_blas; } @@ -185,22 +181,8 @@ TreeID CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_manager, MeshID volume_id) { - - const auto& context = context_; - const int gpu_id = context.gpuID; - - - // cuBQL::box3d *d_boxes = nullptr; // box3d is an alias for box_t - // int num_boxes = 0; - - // cuBQL::bvh3d bvh; // bvh3d is an alias for BinaryBVH - // bvh_t is an alias for BinaryBVH, so bvh_t is also BinaryBVH - - // Looks like the cuBQL::gpuBuilder is only pulled in when cubql is built with CUDA support - - // It looks like spatial median is the only supported method for omp builder - - // cuBQL::build_omp_target(bvh, d_boxes, num_boxes, buildParams, gpu_id); // build bvh on GPU 0 using OpenMP target offloading + // 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); @@ -249,24 +231,24 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man } auto* d_tlas_boxes = static_cast - (omp_target_alloc(h_tlas_boxes.size() * sizeof(cuBQL::box3d), gpu_id)); + (omp_target_alloc(h_tlas_boxes.size() * sizeof(cuBQL::box3d), context_.gpuID)); omp_target_memcpy(d_tlas_boxes, h_tlas_boxes.data(), h_tlas_boxes.size() * sizeof(cuBQL::box3d), 0, 0, - gpu_id, - context.hostID); + context_.gpuID, + context_.hostID); auto* d_surface_instances = static_cast - (omp_target_alloc(h_surface_instances.size() * sizeof(CuBQLVolumeTLAS::SurfaceInstanceDD), gpu_id)); + (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, - gpu_id, - context.hostID); + context_.gpuID, + context_.hostID); cuBQL::BuildConfig tlasBuildParams; tlasBuildParams.makeLeafThreshold = 1; @@ -275,15 +257,15 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man 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 = gpu_id; + 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, - gpu_id); + context_.gpuID); - omp_target_free(d_tlas_boxes, gpu_id); + 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)); 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() From 7f860dbc1a21b6d0d2784eaf26751a5c481efaf9 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 16 Jun 2026 13:31:46 +0100 Subject: [PATCH 43/44] Implemented concrete mixed precision traversal --- include/xdg/cuBQL/ray_tracer.h | 4 ++- include/xdg/cuBQL/triangles.h | 8 ++--- src/cuBQL/intersection.cpp | 47 ++++++++++++++++------------- src/cuBQL/ray_tracer.cpp | 54 ++++++++++++++++++++++------------ 4 files changed, 68 insertions(+), 45 deletions(-) diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h index 388ef29f..c422ecb3 100644 --- a/include/xdg/cuBQL/ray_tracer.h +++ b/include/xdg/cuBQL/ray_tracer.h @@ -73,7 +73,9 @@ class CuBQLRayTracer : public RayTracer { private: CuBQLSurfaceBLAS - register_surface(const std::shared_ptr& mesh_manager, MeshID surface_id); + register_surface(const std::shared_ptr& mesh_manager, + MeshID surface_id, + double bounding_box_bump); void upload_volume_to_tlas_table_(); diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h index 2b9f62cb..20e53c91 100644 --- a/include/xdg/cuBQL/triangles.h +++ b/include/xdg/cuBQL/triangles.h @@ -66,10 +66,10 @@ struct CuBQLSurfaceMesh { struct CuBQLSurfaceBLAS { struct DD { CuBQLSurfaceMesh::DD mesh; // Mesh data device handle - cuBQL::bvh3d bvh; // B:AS device handle + cuBQL::bvh3f bvh; // BLAS device handle }; - cuBQL::bvh3d bvh; // BLAS host handle + cuBQL::bvh3f bvh; // BLAS host handle CuBQLSurfaceMesh mesh; // Surface mesh host owner uint32_t num_prims {0}; @@ -103,11 +103,11 @@ struct CuBQLVolumeTLAS { struct DD { MeshID volume_id {ID_NONE}; const SurfaceInstanceDD* surface_instances {nullptr}; - cuBQL::bvh3d bvh; // TLAS device handle + cuBQL::bvh3f bvh; // TLAS device handle }; MeshID volume_id {ID_NONE}; - cuBQL::bvh3d bvh; // TLAS host handle + cuBQL::bvh3f bvh; // TLAS host handle SurfaceInstanceDD* d_surface_instances {nullptr}; uint32_t num_surface_instances {0}; diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp index eccb2f91..cc5b1bdb 100644 --- a/src/cuBQL/intersection.cpp +++ b/src/cuBQL/intersection.cpp @@ -12,37 +12,40 @@ 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 ray, + CuBQLRay intersection_ray, CuBQLSurfaceHit* hit, int orientation, const MeshID* exclude_primitives, int exclude_count) { - cuBQL::ray3d world_ray; - world_ray.origin = ray.origin; - world_ray.direction = ray.direction; - world_ray.tMin = ray.tMin; - world_ray.tMax = hit->distance; + // 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, &world_ray] - (cuBQL::ray3d& out_ray, cuBQL::bvh3d& out_bvh, int instance_id) + 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 = world_ray; + out_ray = traversal_ray; out_bvh = surface_instance.surface_blas.bvh; }; - auto intersect_prim = [=, &world_ray, &surface_instance] - (uint32_t prim_id) -> double + 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 world_ray.tMax; + return traversal_ray.tMax; } } @@ -61,18 +64,18 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, normal = -normal; } - const double normal_dot_direction = dot(normal, world_ray.direction); + const double normal_dot_direction = dot(normal, intersection_ray.direction); if (orientation_cull(normal_dot_direction, static_cast(orientation))) { - return world_ray.tMax; + return traversal_ray.tMax; } auto intersection = plucker_ray_tri_intersect(vertices, - world_ray.origin, - world_ray.direction, - world_ray.tMax, - world_ray.tMin, + intersection_ray.origin, + intersection_ray.direction, + hit->distance, + intersection_ray.tMin, false, 0); @@ -81,10 +84,12 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, hit->surface = mesh.surface_id; hit->primitive = primitive_ref; hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; - world_ray.tMax = intersection.t; + traversal_ray.tMax = static_cast(intersection.t); } - return world_ray.tMax; + // 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 {}; @@ -93,7 +98,7 @@ static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, leave_blas, intersect_prim, volume_tlas.bvh, - world_ray); + traversal_ray); } #pragma omp end declare target diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp index 94e87b2e..c41bc4f2 100644 --- a/src/cuBQL/ray_tracer.cpp +++ b/src/cuBQL/ray_tracer.cpp @@ -11,6 +11,8 @@ #include "cuBQL/queries/triangleData/math/rayTriangleIntersections.h" #include "cuBQL/traversal/rayQueries.h" +#include + namespace xdg { CuBQLRayTracer::CuBQLRayTracer() @@ -78,7 +80,8 @@ CuBQLRayTracer::register_volume(const std::shared_ptr& mesh_manager CuBQLSurfaceBLAS CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manager, - MeshID surface_id) + 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); @@ -129,12 +132,13 @@ CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manage context_.gpuID, context_.hostID); - auto* d_aabbs = static_cast - (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3d), context_.gpuID)); + 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) + #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]; @@ -148,13 +152,17 @@ CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manage aabb.extend(B); aabb.extend(C); - d_aabbs[primID] = aabb; + 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::bvh3d bvh; + cuBQL::bvh3f bvh; cuBQL::build_omp_target(bvh, d_aabbs, num_faces, blasBuildParams, context_.gpuID); omp_target_free(d_aabbs, context_.gpuID); @@ -187,28 +195,36 @@ 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_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); + surface_to_blas_map_[surf] = register_surface(mesh_manager, surf, max_parent_bbox_bump); } CuBQLSurfaceBLAS& surface_blas = surface_to_blas_map_.at(surf); - auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); // Store BLAS bounding boxes to build TLAS const auto surface_bounding_box = mesh_manager->surface_bounding_box(surf); - cuBQL::box3d surface_bounds; - surface_bounds.lower = cuBQL::vec3d(surface_bounding_box.min_x, - surface_bounding_box.min_y, - surface_bounding_box.min_z); - surface_bounds.upper = cuBQL::vec3d(surface_bounding_box.max_x, - surface_bounding_box.max_y, - surface_bounding_box.max_z); + 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(); @@ -230,11 +246,11 @@ CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_man 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::box3d), context_.gpuID)); + 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::box3d), + h_tlas_boxes.size() * sizeof(cuBQL::box3f), 0, 0, context_.gpuID, From 59d99f49af9e6b8e49d4bc71b87c869f23d42e2e Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 18 Jun 2026 14:34:54 +0100 Subject: [PATCH 44/44] Re-factored ray_benchmark to work with new changes after rebase from main --- tools/ray_benchmark.cpp | 257 +++++++++++----------------------------- 1 file changed, 69 insertions(+), 188 deletions(-) diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index 76f63428..847389d7 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -1,29 +1,16 @@ #include #include -<<<<<<< HEAD #include -======= ->>>>>>> 1003ecb (Working ray_benchmark miniapp) #include #include #include #include #include -<<<<<<< HEAD #include "argparse/argparse.hpp" #include #include "xdg/config.h" -======= -#ifdef XDG_OPENMP -#include -#endif - -#include "argparse/argparse.hpp" -#include - ->>>>>>> 1003ecb (Working ray_benchmark miniapp) #include "xdg/constants.h" #include "xdg/error.h" #include "xdg/timer.h" @@ -33,23 +20,16 @@ #include "ray_benchmark.h" #ifdef XDG_ENABLE_CUBQL +#include #include "xdg/cuBQL/intersection.h" #include "xdg/cuBQL/ray_tracer.h" #endif -<<<<<<< HEAD - -======= ->>>>>>> 1003ecb (Working ray_benchmark miniapp) using namespace xdg; int main(int argc, char** argv) { -<<<<<<< HEAD argparse::ArgumentParser args("XDG Raytracing throughput benchmarking tool", -======= - argparse::ArgumentParser args("XDG Ray Tracing throughput benchmarking tool", ->>>>>>> 1003ecb (Working ray_benchmark miniapp) "1.0", argparse::default_arguments::help); @@ -71,7 +51,6 @@ int main(int argc, char** argv) .scan<'u', std::uint32_t>(); args.add_argument("-o", "-p", "--origin", "--position") -<<<<<<< HEAD .help("Ray origin/position. Defaults to the center of the model bounding box") .scan<'g', double>() .nargs(3); @@ -81,19 +60,12 @@ int main(int argc, char** argv) .implicit_value(true) .help("Use the queried volume bounding box center as the ray origin"); -======= - .default_value(std::vector{0.0, 0.0, 0.0}) - .help("Ray origin/position") - .scan<'g', double>() - .nargs(3); - ->>>>>>> 1003ecb (Working ray_benchmark miniapp) args.add_argument("-m", "--mesh-library") .help("Mesh library to use. One of (MOAB, LIBMESH)") .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") @@ -106,14 +78,11 @@ int main(int argc, char** argv) .help("Radius of a scattered source around the origin") .scan<'g', double>(); -<<<<<<< HEAD args.add_argument("--format") .default_value("human") .choices("human", "csv") .help("stdout format. Human readable (default) or csv"); -======= ->>>>>>> 1003ecb (Working ray_benchmark miniapp) args.add_description( "Benchmarks ray-fire throughput for a selected mesh volume. A source " "position is provided and ray directions are randomly generated from it."); @@ -150,18 +119,12 @@ int main(int argc, char** argv) } const MeshID volume = args.get("volume"); -<<<<<<< HEAD const std::string model_filename = args.get("filename"); const std::string model_name = std::filesystem::path(model_filename).filename().string(); const std::size_t num_rays = args.get("--num-rays"); const std::uint32_t seed = args.get("--seed"); const double source_radius = args.get("--source-radius"); const std::string output_format = args.get("--format"); -======= - const std::size_t num_rays = args.get("--num-rays"); - const std::uint32_t seed = args.get("--seed"); - const double source_radius = args.get("--source-radius"); ->>>>>>> 1003ecb (Working ray_benchmark miniapp) Timer wall_timer; Timer setup_timer; @@ -174,7 +137,6 @@ int main(int argc, char** argv) setup_timer.start(); std::shared_ptr xdg = XDG::create(mesh_lib, rt_lib); const auto& mesh_manager = xdg->mesh_manager(); -<<<<<<< HEAD mesh_manager->load_file(model_filename); mesh_manager->init(); @@ -216,91 +178,14 @@ int main(int argc, char** argv) setup_timer.stop(); - 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(); - - // Trace rays - trace_timer.start(); - std::size_t num_hits = 0; + if (num_rays < 1) fatal_error("Number of rays must be greater than 0"); - #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++; - } - - 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; - const double setup_time = setup_timer.elapsed(); -======= - mesh_manager->load_file(args.get("filename")); - mesh_manager->init(); - - if (args.get("--list")) { - std::cout << "[" << fmt::format("{}", fmt::join(mesh_manager->volumes(), ", ")) << "]\n"; - return 0; - } - - const auto origin_arg = args.present>("--origin"); - bool use_volume_center = args.get("--volume-center"); - if (origin_arg && use_volume_center) { - warning("--volume-center enabled but an explicit origin was also provided. The explicit origin will be used and the volume center will be ignored."); - use_volume_center = false; - } - - Position origin = mesh_manager->global_bounding_box().center(); - if (origin_arg) { - origin = Position(origin_arg.value()); - } else if (use_volume_center) { - origin = mesh_manager->volume_bounding_box(volume).center(); - } - - xdg->prepare_volume_for_raytracing(volume); - xdg->ray_tracing_interface()->init(); - setup_timer.stop(); 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); - - if (rt_lib == RTLibrary::EMBREE) { // Generate random rays from source generation_timer.start(); std::vector origins(num_rays); @@ -323,86 +208,101 @@ int main(int argc, char** argv) } generation_timer.stop(); - // Trace rays - trace_timer.start(); + std::vector hit_surfaces(num_rays, ID_NONE); + trace_timer.start(); #pragma omp parallel for schedule(runtime) for (std::size_t i = 0; i < num_rays; ++i) { - (void) xdg->ray_fire(volume, origins[i], directions[i]); + hit_surfaces[i] = xdg->ray_fire(volume, origins[i], directions[i]).second; // just return surface id of hit } - trace_timer.stop(); + + // 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)."); + 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 and trace in one step - generation_timer.start(); + // Generate random rays directly on the target device. + generation_timer.start(); - const int gpu_id = omp_get_default_device(); + const int gpu_id = omp_get_default_device(); - CuBQLRay* d_rays = static_cast( - omp_target_alloc(num_rays * sizeof(CuBQLRay), gpu_id)); + 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"); - } + 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; + 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); + #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); + 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; - } + 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; - CuBQLSurfaceHit* d_hits = static_cast( - omp_target_alloc(num_rays * sizeof(CuBQLSurfaceHit), gpu_id)); + d_rays[ray_id] = ray; + } - if (!d_hits) { - omp_target_free(d_rays, gpu_id); - fatal_error("Failed to allocate cuBQL hit buffer"); - } + CuBQLSurfaceHit* d_hits = static_cast( + omp_target_alloc(num_rays * sizeof(CuBQLSurfaceHit), gpu_id)); - generation_timer.stop(); + if (!d_hits) { + omp_target_free(d_rays, gpu_id); + fatal_error("Failed to allocate cuBQL hit buffer"); + } - // Trace rays - trace_timer.start(); - rti->ray_fire_batch(d_rays, d_hits, num_rays); - trace_timer.stop(); + 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(); - omp_target_free(d_hits, gpu_id); - omp_target_free(d_rays, gpu_id); + #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 } + 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; ->>>>>>> 1003ecb (Working ray_benchmark miniapp) + const double setup_time = setup_timer.elapsed(); const double trace_only_rps = trace_time > 0.0 ? static_cast(num_rays) / trace_time : 0.0; @@ -411,7 +311,6 @@ int main(int argc, char** argv) : 0.0; wall_timer.stop(); -<<<<<<< HEAD const double wall_time = wall_timer.elapsed(); const std::vector csv_columns { @@ -499,24 +398,6 @@ int main(int argc, char** argv) std::cout << "End-to-end throughput : " << end_to_end_rps << " rays/s\n"; std::cout << "Trace-only throughput : " << trace_only_rps << " rays/s\n"; } -======= - - std::cout << "Random ray generation time = " - << generation_time << "s" << std::endl; - std::cout << "Generation + tracing time = " - << end_to_end_time << "s" << std::endl; - std::cout << "End-to-end throughput = " - << end_to_end_rps << " rays/s" << std::endl; - std::cout << "Full wall-clock time = " - << wall_timer.elapsed() << "s (post-argparse)" << std::endl; - - std::cout << "----------------------------------------" << std::endl; - std::cout << "Ray tracing time (trace-only)= " - << trace_time << "s for " << num_rays << " rays" << std::endl; - std::cout << "Trace-only throughput = " - << trace_only_rps << " rays/s" << std::endl; - std::cout << "----------------------------------------" << std::endl; ->>>>>>> 1003ecb (Working ray_benchmark miniapp) return 0; }