diff --git a/.gitmodules b/.gitmodules index 77e2a79a..dd0e33bd 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,3 +20,6 @@ [submodule "vendor/GPRT"] path = vendor/GPRT url = https://github.com/gprt-org/GPRT.git +[submodule "vendor/cuBQL"] + path = vendor/cuBQL + url = https://github.com/NVIDIA/cuBQL diff --git a/CMakeLists.txt b/CMakeLists.txt index df2d6cb3..7196b384 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,6 +7,7 @@ option(XDG_ENABLE_LIBMESH "Enable support for the libMesh mesh library" OFF option(XDG_LINK_MPI "Link with MPI (for dependency compatibility)" OFF) option(XDG_ENABLE_EMBREE "Enable support for the Embree ray tracing library" ON) option(XDG_ENABLE_GPRT "Enable support for the GPRT ray tracing library" OFF) +option(XDG_ENABLE_CUBQL "Enable support for the cuBQL ray tracing library" OFF) option(XDG_BUILD_TESTS "Enable C++ unit testing" ON) option(XDG_BUILD_TOOLS "Enable tools and miniapps" ON) @@ -20,6 +21,10 @@ if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose build type" FORCE) endif() +if(DEFINED XDG_CMAKE_PRESET) + message(STATUS "XDG CMake preset: ${XDG_CMAKE_PRESET}") +endif() + # Compiler options (things in this section may not be platform-portable) set(CMAKE_POSITION_INDEPENDENT_CODE ON) @@ -111,6 +116,12 @@ if(XDG_ENABLE_GPRT) ) endif() +if(XDG_ENABLE_CUBQL) + list(APPEND VENDOR_PATHS + vendor/cuBQL + ) +endif() + if(GIT_FOUND AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/.git") option(XDG_GIT_SUBMODULE "Check submodules during build" ON) if(XDG_GIT_SUBMODULE) @@ -175,11 +186,12 @@ if (NOT XDG_ENABLE_MOAB AND NOT XDG_ENABLE_LIBMESH) endif() # Ensure at least one ray tracing backend is enabled -if (NOT XDG_ENABLE_EMBREE AND NOT XDG_ENABLE_GPRT) +if (NOT XDG_ENABLE_EMBREE AND NOT XDG_ENABLE_GPRT AND NOT XDG_ENABLE_CUBQL) message(FATAL_ERROR "No ray tracing backend enabled. Enable at least one of:\n" " -DXDG_ENABLE_EMBREE=ON\n" - " -DXDG_ENABLE_GPRT=ON") + " -DXDG_ENABLE_GPRT=ON\n" + " -DXDG_ENABLE_CUBQL=ON") endif() # GPRT @@ -188,6 +200,11 @@ if (XDG_ENABLE_GPRT) add_subdirectory(vendor/GPRT) endif() +if (XDG_ENABLE_CUBQL) + set(CUBQL_OMP ON CACHE BOOL "Build cuBQL with OpenMP target offload" FORCE) + add_subdirectory(vendor/cuBQL) +endif() + list(APPEND xdg_sources src/geometry/measure.cpp src/geometry/plucker.cpp @@ -221,6 +238,23 @@ dbl_deviceCode endif() +if (XDG_ENABLE_CUBQL) +list(APPEND xdg_sources +src/cuBQL/triangles.cpp +src/cuBQL/intersection.cpp +src/cuBQL/ray_tracer.cpp +) + +# We need a precompile definition to switch to using the cuBQL math types in the shared +# plucker intersection code. The compile definition is used in dp__math.h +set_source_files_properties( + src/cuBQL/ray_tracer.cpp + src/cuBQL/intersection.cpp + PROPERTIES COMPILE_DEFINITIONS XDG_DP_MATH_CUBQL +) + +endif() + if (XDG_ENABLE_LIBMESH) list(APPEND xdg_sources src/libmesh/mesh_manager.cpp @@ -281,7 +315,12 @@ if (${CMAKE_BUILD_TYPE} MATCHES "Debug") endif() # attempt to find OpenMP and include it if found -find_package(OpenMP) +if (XDG_ENABLE_CUBQL) + find_package(OpenMP REQUIRED) +else() + find_package(OpenMP) +endif() + if (OpenMP_CXX_FOUND) target_link_libraries(xdg PRIVATE OpenMP::OpenMP_CXX) target_compile_definitions(xdg PRIVATE XDG_HAVE_OPENMP) @@ -319,6 +358,18 @@ if (XDG_ENABLE_GPRT) target_link_options(xdg PRIVATE -Wl,--unresolved-symbols=ignore-in-shared-libs) endif() +if (XDG_ENABLE_CUBQL) + target_compile_definitions(xdg PUBLIC XDG_ENABLE_CUBQL) + target_link_libraries(xdg PRIVATE $) + # TODO: Stop relying on LD_LIBRARY_PATH for LLVM OpenMP offload runtimes. + # As a temporary measure whilst figuring out our way around cuBQL this is + # okay but in the long run we should aim for a more robust solutions here. + # Clang injects libomptarget when offload flags are supplied by presets, but + # CMake does not currently add that compiler runtime directory to RPATH. + # Add targeted BUILD_RPATH handling for libomptarget/libomp, and decide + # whether install RPATH should remain environment-module based or be opt-in. +endif() + target_link_libraries(xdg PRIVATE fmt::fmt) # ========================== diff --git a/CMakePresets.json b/CMakePresets.json new file mode 100644 index 00000000..09b86db1 --- /dev/null +++ b/CMakePresets.json @@ -0,0 +1,68 @@ +{ + "version": 3, + "configurePresets": [ + { + "name": "base", + "generator": "Unix Makefiles", + "binaryDir": "${sourceDir}/build/${presetName}", + "cacheVariables": { + "CMAKE_CXX_FLAGS": "$env{PRESET_CXX_FLAGS} $env{LLVM_CXX_FLAGS} $env{COMMON_CXX_FLAGS}", + "XDG_CMAKE_PRESET": "${presetName}" + }, + "environment": { + "COMMON_CXX_FLAGS": "" + } + }, + { + "name": "llvm", + "inherits": ["base"], + "cacheVariables": { + "CMAKE_C_COMPILER": "clang", + "CMAKE_CXX_COMPILER": "clang++" + }, + "environment": { + "LLVM_CXX_FLAGS": "-fopenmp -fopenmp-cuda-mode" + } + }, + { + "name": "cubql_llvm_ada", + "inherits": ["llvm"], + "displayName": "cuBQL LLVM OpenMP offload RTX 2000 Ada", + "cacheVariables": { + "XDG_ENABLE_CUBQL": "ON" + }, + "environment": { + "PRESET_CXX_FLAGS": "-fopenmp-targets=nvptx64 -Xopenmp-target -march=sm_80" + } + }, + { + "name": "nvhpc", + "inherits": ["base"], + "cacheVariables": { + "CMAKE_C_COMPILER": "nvc", + "CMAKE_CXX_COMPILER": "nvc++" + } + }, + { + "name": "cubql_nvhpc_ada", + "inherits": ["nvhpc"], + "displayName": "cuBQL NVHPC OpenMP offload RTX 2000 Ada", + "cacheVariables": { + "XDG_ENABLE_CUBQL": "ON" + }, + "environment": { + "PRESET_CXX_FLAGS": "-mp=gpu -Minfo=mp -gpu=cc89" + } + } + ], + "buildPresets": [ + { + "name": "cubql_llvm_ada", + "configurePreset": "cubql_llvm_ada" + }, + { + "name": "cubql_nvhpc_ada", + "configurePreset": "cubql_nvhpc_ada" + } + ] +} 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/include/xdg/constants.h b/include/xdg/constants.h index 38df9db1..e143d537 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -54,7 +54,8 @@ enum class MeshLibrary { // Ray Tracing library identifier enum class RTLibrary { EMBREE, - GPRT + GPRT, + CUBQL }; static const std::map MESH_LIB_TO_STR = @@ -67,7 +68,8 @@ static const std::map MESH_LIB_TO_STR = static const std::map RT_LIB_TO_STR = { {RTLibrary::EMBREE, "EMBREE"}, - {RTLibrary::GPRT, "GPRT"} + {RTLibrary::GPRT, "GPRT"}, + {RTLibrary::CUBQL, "CUBQL"} }; // Mesh identifer type @@ -148,4 +150,4 @@ struct formatter : fmt::formatter { } -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/cuBQL/cuBQL_backend.h b/include/xdg/cuBQL/cuBQL_backend.h new file mode 100644 index 00000000..b348802f --- /dev/null +++ b/include/xdg/cuBQL/cuBQL_backend.h @@ -0,0 +1,30 @@ +#ifndef _XDG_CUBQL_BACKEND_H +#define _XDG_CUBQL_BACKEND_H + +#include +#include +#include +#include + +#include + +// Guards to prevent CUDA headers from being included in host code, which causes +// failed compilation with LLVM-clang. +#if defined(__CUDA_ARCH__) && !defined(__CUDACC__) +#undef __CUDA_ARCH__ +#endif + +#include "xdg/error.h" + +namespace xdg::cubql { + +struct Context { + int gpuID {0}; + int hostID {omp_get_initial_device()}; +}; + + + +} // namespace xdg::cubql + +#endif // include guard diff --git a/include/xdg/cuBQL/intersection.h b/include/xdg/cuBQL/intersection.h new file mode 100644 index 00000000..50d95141 --- /dev/null +++ b/include/xdg/cuBQL/intersection.h @@ -0,0 +1,89 @@ +#ifndef _XDG_CUBQL_INTERSECTION_H +#define _XDG_CUBQL_INTERSECTION_H + +// Guards to prevent CUDA headers from being included in host code, which causes +// failed compilation with LLVM-clang. +#if defined(__CUDA_ARCH__) && !defined(__CUDACC__) +#undef __CUDA_ARCH__ +#endif + +#include +#include + +#include "xdg/constants.h" +#include "xdg/cuBQL/triangles.h" +#include "cuBQL/math/vec.h" + +namespace xdg { + +struct CuBQLRay { + cuBQL::vec3d origin; + cuBQL::vec3d direction; + double tMin {0.0}; + double tMax {INFTY}; + MeshID volume {ID_NONE}; // volume we are tracing ray against +}; + +/* POD SurfaceRay struct for external population*/ +// struct CuBQLSurfaceRay { +// double origin[3]; +// double direction[3]; +// uint32_t volume_slot; +// uint32_t enabled; +// const MeshID* exclude_primitives; +// int32_t exclude_count; +// }; + +// TODO - Consider whether this is useful/necessary as its own struct +// struct CuBQLExcludeList { +// const MeshID* primitives {nullptr}; +// int count {0}; +// }; + +struct CuBQLSurfaceHit { + double distance {INFTY}; + MeshID surface {ID_NONE}; + MeshID primitive {ID_NONE}; + PointInVolume piv {OUTSIDE}; + + bool hit_found() const { return primitive != ID_NONE; } +}; + +inline bool orientation_cull(double normal_dot_direction, + HitOrientation orientation) +{ + if (orientation == HitOrientation::ANY) return false; + + if (orientation == HitOrientation::EXITING && normal_dot_direction < 0.0) { + return true; + } else if (orientation == HitOrientation::ENTERING && normal_dot_direction >= 0.0) { + return true; + } + + return false; +} + +/* +Wrapper for launching a single ray intersection query against the surface tree, with Host<->Device staging of ray and hit data +Performs host side staging and transfer hit data back to host after device side traversal +*/ +void +intersect_surface_tree_scalar(const cubql::Context& context, + const CuBQLVolumeTLAS& volume_tlas, + const CuBQLRay& ray, + CuBQLSurfaceHit& hit, + HitOrientation hit_orientation, + const std::vector* exclude_primitives); + + +void +intersect_surface_tree_batch(const cubql::Context& context, + const CuBQLVolumeTLAS::DD* d_volume_to_tlas, + const CuBQLRay* d_rays, + CuBQLSurfaceHit* d_hits, + std::size_t num_rays, + HitOrientation hit_orientation); + +} // namespace xdg + +#endif // include guard diff --git a/include/xdg/cuBQL/ray_tracer.h b/include/xdg/cuBQL/ray_tracer.h new file mode 100644 index 00000000..c422ecb3 --- /dev/null +++ b/include/xdg/cuBQL/ray_tracer.h @@ -0,0 +1,94 @@ +#ifndef _XDG_CUBQL_RAY_TRACING_INTERFACE_H +#define _XDG_CUBQL_RAY_TRACING_INTERFACE_H + +#include +#include +#include +#include +#include + +#include "xdg/constants.h" +#include "xdg/geometry_data.h" +#include "xdg/cuBQL/triangles.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/ray.h" +#include "xdg/ray_tracing_interface.h" + +namespace xdg { + +struct CuBQLRay; +struct CuBQLSurfaceHit; + +class CuBQLRayTracer : public RayTracer { +public: + CuBQLRayTracer(); + ~CuBQLRayTracer() override; + + RTLibrary library() const override { return RTLibrary::CUBQL; } + + void init() override; + + std::pair + register_volume(const std::shared_ptr& mesh_manager, + MeshID volume) override; + + TreeID create_surface_tree(const std::shared_ptr& mesh_manager, + MeshID volume) override; + + TreeID create_element_tree(const std::shared_ptr& mesh_manager, + MeshID volume) override; + + void create_global_surface_tree() override; + + void create_global_element_tree() override; + + MeshID find_element(const Position& point) const override; + + MeshID find_element(TreeID tree, const Position& point) const override; + + bool point_in_volume(TreeID tree, + const Position& point, + const Direction* direction = nullptr, + const std::vector* exclude_primitives = nullptr) const override; + + std::pair ray_fire(TreeID tree, + const Position& origin, + const Direction& direction, + const double dist_limit = INFTY, + HitOrientation orientation = HitOrientation::EXITING, + std::vector* const exclude_primitives = nullptr) override; + + void ray_fire_batch(const CuBQLRay* d_rays, + CuBQLSurfaceHit* d_hits, + std::size_t num_rays, + HitOrientation orientation = HitOrientation::EXITING); + + std::pair closest(TreeID tree, + const Position& origin) override; + + bool occluded(TreeID tree, + const Position& origin, + const Direction& direction, + double& dist) const override; + +private: + CuBQLSurfaceBLAS + register_surface(const std::shared_ptr& mesh_manager, + MeshID surface_id, + double bounding_box_bump); + + void upload_volume_to_tlas_table_(); + + cubql::Context context_; + + std::unordered_map tree_to_volume_tlas_; + std::unordered_map surface_to_blas_map_; + + std::vector volume_to_tlas_; + CuBQLVolumeTLAS::DD* d_volume_to_tlas_ {nullptr}; + bool initialized_ {false}; +}; + +} // namespace xdg + +#endif // include guard diff --git a/include/xdg/cuBQL/triangles.h b/include/xdg/cuBQL/triangles.h new file mode 100644 index 00000000..20e53c91 --- /dev/null +++ b/include/xdg/cuBQL/triangles.h @@ -0,0 +1,126 @@ +#ifndef _XDG_CUBQL_TRIANGLES_H +#define _XDG_CUBQL_TRIANGLES_H + +#include +#include + +// Guards to prevent CUDA headers from being included in host code, which causes +// failed compilation with LLVM-clang. +#if defined(__CUDA_ARCH__) && !defined(__CUDACC__) +#undef __CUDA_ARCH__ +#endif + +#include "cuBQL/bvh.h" +#include "cuBQL/math/vec.h" +#include "xdg/constants.h" +#include "xdg/cuBQL/cuBQL_backend.h" + +namespace xdg { + +/* + Owns the triangle buffers for one topological surface. The nested DD type is + the compact device-data view copied into OpenMP target regions instead of the + full host-side owner. +*/ +struct CuBQLSurfaceMesh { + struct DD { + // Topological metadata + MeshID surface_id {ID_NONE}; + + // Geometric data + const cuBQL::vec3d* vertices {nullptr}; + const cuBQL::vec3i* indices {nullptr}; + const MeshID* primitive_refs {nullptr}; + }; + + // Topological metadata + MeshID surface_id {ID_NONE}; + + // Device buffers for triangle data + cuBQL::vec3d* d_vertices {nullptr}; + cuBQL::vec3i* d_indices {nullptr}; + MeshID* d_primitive_refs {nullptr}; + + uint32_t num_vertices {0}; + uint32_t num_triangles {0}; + int gpu_id {0}; + + // Accessor for Device Data struct, which is passed to cuBQL BVH traversal/intersection functions + DD get_device_data() const + { + return { + surface_id, + d_vertices, + d_indices, + d_primitive_refs + }; + } + + void release(); +}; + +/* + Owns a cuBQL BVH used as a Bottom Level Acceleration Structure over surface triangles. + The nested DD type is the compact device-data view used during traversal. +*/ +struct CuBQLSurfaceBLAS { + struct DD { + CuBQLSurfaceMesh::DD mesh; // Mesh data device handle + cuBQL::bvh3f bvh; // BLAS device handle + }; + + cuBQL::bvh3f bvh; // BLAS host handle + CuBQLSurfaceMesh mesh; // Surface mesh host owner + + uint32_t num_prims {0}; + int gpu_id {0}; + + DD get_device_data() const + { + return {mesh.get_device_data(), bvh}; + } + + void release(); +}; + +/* + Owns a cuBQL BVH used as a Top-Level Acceleration Structure for one topological volume. + The TLAS groups the surface BLASes that bound that volume and stores + per-volume relationship metadata for each surface instance. +*/ +struct CuBQLVolumeTLAS { + /* + TLAS-local instance payload. The same surface BLAS can participate in + different volume TLASes with different sense, so reverse_sense belongs on + the volume-surface relationship rather than on the reusable surface mesh + or BLAS geometry. + */ + struct SurfaceInstanceDD { + CuBQLSurfaceBLAS::DD surface_blas; + bool reverse_sense {false}; // value set in create_surface_tree based on parent vols + }; + + struct DD { + MeshID volume_id {ID_NONE}; + const SurfaceInstanceDD* surface_instances {nullptr}; + cuBQL::bvh3f bvh; // TLAS device handle + }; + + MeshID volume_id {ID_NONE}; + cuBQL::bvh3f bvh; // TLAS host handle + SurfaceInstanceDD* d_surface_instances {nullptr}; + + uint32_t num_surface_instances {0}; + int gpu_id {0}; + + DD get_device_data() const + { + return {volume_id, d_surface_instances, bvh}; + } + + void release(); +}; + +} // namespace xdg + +#endif // include guard diff --git a/include/xdg/geometry/dp_math.h b/include/xdg/geometry/dp_math.h index 5b162d3f..5dfedffd 100644 --- a/include/xdg/geometry/dp_math.h +++ b/include/xdg/geometry/dp_math.h @@ -22,7 +22,27 @@ namespace dp { static const double INFTY = 1.7976931348623157e+308; // std::numeric_limits::max() is not available in slang } -#else +// TODO - Is this the right way to handle this for cubql openmp target offload compilation? +// In theory, we can compile the C++ pathway that Embree uses but will the vec3da types be +// omptarget friendly? +// For now I have defined this separate compilation pathway which is enabled with a new +// precompile definition only set when when compiling the CuBQLRayTracer that maps to cuBQL's math types +#elif defined(XDG_DP_MATH_CUBQL) +#include "cuBQL/math/common.h" + +// C++ compilation with cuBQL math types, map dp::vec3 -> cuBQL::vec3d +namespace dp { + typedef cuBQL::vec3d vec3; + + inline double dot(vec3 a, vec3 b) { return cuBQL::dot(a, b); } + inline vec3 cross(vec3 a, vec3 b) { return cuBQL::cross(a, b); } + inline double abs(double a) { return cuBQL::abst(a); } + + static constexpr double DBL_ZERO_TOL = 20.0 * 2.2204460492503131e-16; // same as 20 * std::numeric_limits::epsilon() + constexpr double INFTY {std::numeric_limits::max()}; +} + +#else #include "xdg/vec3da.h" // C++ compilation map dp::vec3 -> xdg::Vec3da diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index acf0ebfa..8b7581c9 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -25,8 +25,6 @@ struct PluckerIntersectionResult { double t = 0.0; // Distance along the ray to the intersection point }; -static constexpr PluckerIntersectionResult EXIT_EARLY = {false, 0.0}; - /* Function to return the vertex with the lowest coordinates. To force the same ray-edge computation, the Plücker test needs to use consistent edge representation. This would be more simple with MOAB handles instead of @@ -81,7 +79,7 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], // If orientation is set, confirm that sign of plucker_coordinate indicate // correct orientation of intersection if (useOrientation && orientation * plucker_coord0 > 0) { - return EXIT_EARLY; + return {false, 0.0}; } // Determine the value of the second Plucker coordinate from edge 1 @@ -92,13 +90,13 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], // correct orientation of intersection if (useOrientation) { if (orientation * plucker_coord1 > 0) { - return EXIT_EARLY; + return {false, 0.0}; } // If the orientation is not specified, all plucker_coords must be the same // sign or zero. } else if ((0.0 < plucker_coord0 && 0.0 > plucker_coord1) || (0.0 > plucker_coord0 && 0.0 < plucker_coord1)) { - return EXIT_EARLY; + return {false, 0.0}; } // Determine the value of the third Plucker coordinate from edge 2 @@ -109,7 +107,7 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], // correct orientation of intersection if (useOrientation) { if (orientation * plucker_coord2 > 0) { - return EXIT_EARLY; + return {false, 0.0}; } // If the orientation is not specified, all plucker_coords must be the same // sign or zero. @@ -117,12 +115,12 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], (0.0 > plucker_coord1 && 0.0 < plucker_coord2) || (0.0 < plucker_coord0 && 0.0 > plucker_coord2) || (0.0 > plucker_coord0 && 0.0 < plucker_coord2)) { - return EXIT_EARLY; + return {false, 0.0}; } // check for coplanar case to avoid dividing by zero if (0.0 == plucker_coord0 && 0.0 == plucker_coord1 && 0.0 == plucker_coord2) { - return EXIT_EARLY; + return {false, 0.0}; } // get the distance to intersection @@ -155,7 +153,7 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], } // is the intersection within distance limits? - if (dist_out < tMin || dist_out > tMax) return EXIT_EARLY; + if (dist_out < tMin || dist_out > tMax) return {false, 0.0}; return {true, dist_out}; } diff --git a/include/xdg/ray_tracers.h b/include/xdg/ray_tracers.h index b4816571..2031baf0 100644 --- a/include/xdg/ray_tracers.h +++ b/include/xdg/ray_tracers.h @@ -5,4 +5,8 @@ #ifdef XDG_ENABLE_GPRT #include "xdg/gprt/ray_tracer.h" -#endif \ No newline at end of file +#endif + +#ifdef XDG_ENABLE_CUBQL +#include "xdg/cuBQL/ray_tracer.h" +#endif diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp new file mode 100644 index 00000000..cc5b1bdb --- /dev/null +++ b/src/cuBQL/intersection.cpp @@ -0,0 +1,210 @@ +#include + +#include "xdg/cuBQL/intersection.h" +#include "xdg/geometry/plucker.h" +#include "xdg/error.h" + +#include "cuBQL/math/Ray.h" +#include "cuBQL/traversal/rayQueries.h" + +namespace xdg { + +// Core traversal and intersection routine for a single ray against a given volume tlas +#pragma omp declare target +static inline void intersect_surface_tree(CuBQLVolumeTLAS::DD volume_tlas, + CuBQLRay intersection_ray, + CuBQLSurfaceHit* hit, + int orientation, + const MeshID* exclude_primitives, + int exclude_count) +{ + // cuBQL traverses the FP32 BVH with an FP32 ray; the original CuBQLRay + // remains the FP64 source of truth for the final triangle intersection. + cuBQL::ray3f traversal_ray; + traversal_ray.origin = cuBQL::vec3f(intersection_ray.origin); + traversal_ray.direction = cuBQL::vec3f(intersection_ray.direction); + // TODO: Is this truncation safe enough for tmin and tmax? Pretty sure embree/gprt does a similar truncation for ray bounds + traversal_ray.tMin = static_cast(intersection_ray.tMin); + traversal_ray.tMax = static_cast(hit->distance); + + CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; + + auto enter_blas = [=, &surface_instance, &traversal_ray] + (cuBQL::ray3f& out_ray, cuBQL::bvh3f& out_bvh, int instance_id) + { + surface_instance = volume_tlas.surface_instances[instance_id]; + out_ray = traversal_ray; + out_bvh = surface_instance.surface_blas.bvh; + }; + + auto intersect_prim = [=, &traversal_ray, &surface_instance] + (uint32_t prim_id) -> float + { + const CuBQLSurfaceMesh::DD mesh = surface_instance.surface_blas.mesh; + const MeshID primitive_ref = mesh.primitive_refs[prim_id]; + + for (int i = 0; i < exclude_count; ++i) { + if (exclude_primitives[i] == primitive_ref) { + return traversal_ray.tMax; + } + } + + const cuBQL::vec3i index = mesh.indices[prim_id]; + + cuBQL::vec3d vertices[3] = { + mesh.vertices[index.x], + mesh.vertices[index.y], + mesh.vertices[index.z] + }; + + cuBQL::vec3d normal = cuBQL::cross(vertices[1] - vertices[0], + vertices[2] - vertices[0]); + + if (surface_instance.reverse_sense) { + normal = -normal; + } + + const double normal_dot_direction = dot(normal, intersection_ray.direction); + + if (orientation_cull(normal_dot_direction, + static_cast(orientation))) { + return traversal_ray.tMax; + } + + auto intersection = plucker_ray_tri_intersect(vertices, + intersection_ray.origin, + intersection_ray.direction, + hit->distance, + intersection_ray.tMin, + false, + 0); + + if (intersection.hit) { + hit->distance = intersection.t; + hit->surface = mesh.surface_id; + hit->primitive = primitive_ref; + hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; + traversal_ray.tMax = static_cast(intersection.t); + } + + // Return value is only the FP32 traversal shrink distance. The accepted hit + // distance stored above remains the FP64 Plucker result. + return traversal_ray.tMax; + }; + + auto leave_blas = []() -> void {}; + + cuBQL::shrinkingRayQuery::twoLevel::forEachPrim(enter_blas, + leave_blas, + intersect_prim, + volume_tlas.bvh, + traversal_ray); +} +#pragma omp end declare target + +void +intersect_surface_tree_scalar(const cubql::Context& context, + const CuBQLVolumeTLAS& volume_tlas, + const CuBQLRay& ray, + CuBQLSurfaceHit& surface_hit, + HitOrientation hit_orientation, + const std::vector* exclude_primitives) +{ + const int gpu_id = context.gpuID; + + MeshID* d_exclude_primitives = nullptr; + int exclude_count = 0; + if (exclude_primitives && !exclude_primitives->empty()) { + exclude_count = static_cast(exclude_primitives->size()); + d_exclude_primitives = static_cast + (omp_target_alloc(exclude_count * sizeof(MeshID), gpu_id)); + omp_target_memcpy(d_exclude_primitives, + exclude_primitives->data(), + exclude_count * sizeof(MeshID), + 0, + 0, + gpu_id, + context.hostID); + } + + auto* d_surface_hit = static_cast + (omp_target_alloc(sizeof(CuBQLSurfaceHit), gpu_id)); + + surface_hit.distance = ray.tMax; + omp_target_memcpy(d_surface_hit, + &surface_hit, + sizeof(CuBQLSurfaceHit), + 0, + 0, + gpu_id, + context.hostID); + + const auto volume_tlas_dd = volume_tlas.get_device_data(); + const int orientation = static_cast(hit_orientation); + + #pragma omp target device(gpu_id) \ + is_device_ptr(d_exclude_primitives, d_surface_hit) + { + intersect_surface_tree(volume_tlas_dd, + ray, + d_surface_hit, + orientation, + d_exclude_primitives, + exclude_count); + } + + omp_target_memcpy(&surface_hit, + d_surface_hit, + sizeof(CuBQLSurfaceHit), + 0, + 0, + context.hostID, + gpu_id); + + omp_target_free(d_surface_hit, gpu_id); + + if (d_exclude_primitives) { + omp_target_free(d_exclude_primitives, gpu_id); + } + + return; +} + +void +intersect_surface_tree_batch(const cubql::Context& context, + const CuBQLVolumeTLAS::DD* d_volume_to_tlas, + const CuBQLRay* d_rays, + CuBQLSurfaceHit* d_hits, + std::size_t num_rays, + HitOrientation hit_orientation) +{ + + if (num_rays == 0) return; + + if (!d_volume_to_tlas || !d_rays || !d_hits) { + fatal_error("Invalid cuBQL batch intersection buffers"); + } + + const int gpu_id = context.gpuID; + + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(d_volume_to_tlas, d_rays, d_hits) + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + const CuBQLRay ray = d_rays[ray_id]; + const CuBQLVolumeTLAS::DD volume_tlas = d_volume_to_tlas[ray.volume]; + + CuBQLSurfaceHit hit; + hit.distance = ray.tMax; + + intersect_surface_tree(volume_tlas, + ray, + &hit, + static_cast(hit_orientation), + nullptr, + 0); + + d_hits[ray_id] = hit; + } +} + +} // namespace xdg diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp new file mode 100644 index 00000000..c41bc4f2 --- /dev/null +++ b/src/cuBQL/ray_tracer.cpp @@ -0,0 +1,437 @@ +#include "xdg/cuBQL/ray_tracer.h" +#include "xdg/cuBQL/intersection.h" +#include "xdg/error.h" +#include "xdg/geometry/plucker.h" +#include "xdg/available_device_probe.h" + +#include +#include "cuBQL/builder/omp.h" +#include "cuBQL/math/Ray.h" +#include "cuBQL/queries/triangleData/Triangle.h" +#include "cuBQL/queries/triangleData/math/rayTriangleIntersections.h" +#include "cuBQL/traversal/rayQueries.h" + +#include + +namespace xdg { + +CuBQLRayTracer::CuBQLRayTracer() +{ + if (!system_has_omp_target_device()) { + fatal_error("No OpenMP target capable device found; cannot initialize cuBQL ray tracer."); + } + + context_.gpuID = 0; // TODO - support selecting among multiple OpenMP target devices. + context_.hostID = omp_get_initial_device(); +} + +CuBQLRayTracer::~CuBQLRayTracer() +{ + if (d_volume_to_tlas_) { + omp_target_free(d_volume_to_tlas_, context_.gpuID); + d_volume_to_tlas_ = nullptr; + } + + for (auto& [tree, tlas] : tree_to_volume_tlas_) { + tlas.release(); + } + + for (auto& [surface, blas] : surface_to_blas_map_) { + blas.release(); + } +} + +void CuBQLRayTracer::init() +{ + upload_volume_to_tlas_table_(); + initialized_ = true; +} + +void CuBQLRayTracer::upload_volume_to_tlas_table_() +{ + if (d_volume_to_tlas_) { + omp_target_free(d_volume_to_tlas_, context_.gpuID); + d_volume_to_tlas_ = nullptr; + } + + if (volume_to_tlas_.empty()) { + return; + } + + d_volume_to_tlas_ = static_cast + (omp_target_alloc(volume_to_tlas_.size() * sizeof(CuBQLVolumeTLAS::DD), context_.gpuID)); + omp_target_memcpy(d_volume_to_tlas_, + volume_to_tlas_.data(), + volume_to_tlas_.size() * sizeof(CuBQLVolumeTLAS::DD), + 0, + 0, + context_.gpuID, + context_.hostID); +} + +std::pair +CuBQLRayTracer::register_volume(const std::shared_ptr& mesh_manager, + MeshID volume) +{ + TreeID surface_tree = create_surface_tree(mesh_manager, volume); + TreeID element_tree = create_element_tree(mesh_manager, volume); + return {surface_tree, element_tree}; +} + +CuBQLSurfaceBLAS +CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manager, + MeshID surface_id, + double bounding_box_bump) +{ + auto num_faces = mesh_manager->num_surface_faces(surface_id); + auto vertices = mesh_manager->get_surface_vertices(surface_id); + auto indices = mesh_manager->get_surface_connectivity(surface_id); + + std::vector h_vertices; + h_vertices.reserve(vertices.size()); + for (const auto& vertex : vertices) { + h_vertices.emplace_back(vertex.x, vertex.y, vertex.z); + } + + std::vector h_indices; + h_indices.reserve(indices.size() / 3); + for (size_t i = 0; i < indices.size(); i += 3) { + h_indices.emplace_back(indices[i], indices[i + 1], indices[i + 2]); + } + + std::vector h_primitive_refs = mesh_manager->get_surface_faces(surface_id); + + // TODO- think about how to better handle omp transfer calls. AutoUploadArrays is one option + auto* d_vertices = static_cast + (omp_target_alloc(h_vertices.size() * sizeof(cuBQL::vec3d), context_.gpuID)); + omp_target_memcpy(d_vertices, + h_vertices.data(), + h_vertices.size() * sizeof(cuBQL::vec3d), + 0, + 0, + context_.gpuID, + context_.hostID); + + auto* d_indices = static_cast + (omp_target_alloc(h_indices.size() * sizeof(cuBQL::vec3i), context_.gpuID)); + omp_target_memcpy(d_indices, + h_indices.data(), + h_indices.size() * sizeof(cuBQL::vec3i), + 0, + 0, + context_.gpuID, + context_.hostID); + + auto* d_primitive_refs = static_cast + (omp_target_alloc(h_primitive_refs.size() * sizeof(MeshID), context_.gpuID)); + omp_target_memcpy(d_primitive_refs, + h_primitive_refs.data(), + h_primitive_refs.size() * sizeof(MeshID), + 0, + 0, + context_.gpuID, + context_.hostID); + + auto* d_aabbs = static_cast + (omp_target_alloc(h_indices.size() * sizeof(cuBQL::box3f), context_.gpuID)); + const auto num_primitives = static_cast(h_indices.size()); + + // TODO - Abstract this out into its own bounding_box creation function + #pragma omp target device(context_.gpuID) is_device_ptr(d_vertices, d_indices, d_aabbs) \ + firstprivate(bounding_box_bump) + #pragma omp teams distribute parallel for + for (uint32_t primID = 0; primID < num_primitives; ++primID) { + cuBQL::vec3i indices = d_indices[primID]; + + cuBQL::vec3d A = d_vertices[indices.x]; + cuBQL::vec3d B = d_vertices[indices.y]; + cuBQL::vec3d C = d_vertices[indices.z]; + + cuBQL::box3d aabb; + aabb.extend(A); + aabb.extend(B); + aabb.extend(C); + + const cuBQL::vec3d bump(bounding_box_bump); + aabb.lower = aabb.lower - bump; + aabb.upper = aabb.upper + bump; + + d_aabbs[primID] = cuBQL::box3f(aabb); + } + + cuBQL::BuildConfig blasBuildParams; + // TODO - Try setting leaf params to 1 to see what it does + // Check what default is for CUDA + cuBQL::bvh3f bvh; + cuBQL::build_omp_target(bvh, d_aabbs, num_faces, blasBuildParams, context_.gpuID); + + omp_target_free(d_aabbs, context_.gpuID); + + CuBQLSurfaceMesh surface_mesh; + surface_mesh.surface_id = surface_id; + surface_mesh.d_vertices = d_vertices; + surface_mesh.d_indices = d_indices; + surface_mesh.d_primitive_refs = d_primitive_refs; + surface_mesh.num_vertices = h_vertices.size(); + surface_mesh.num_triangles = num_faces; + surface_mesh.gpu_id = context_.gpuID; + + CuBQLSurfaceBLAS surface_blas; + surface_blas.bvh = bvh; + surface_blas.mesh = surface_mesh; + surface_blas.num_prims = num_faces; + surface_blas.gpu_id = context_.gpuID; + + return surface_blas; +} + +TreeID +CuBQLRayTracer::create_surface_tree(const std::shared_ptr& mesh_manager, + MeshID volume_id) +{ + // TODO - Right now each CuBQLRayTracer instance has a single "Context" which holds a single GPU_ID + // so this will need to be reworked in the future to handle multi-gpus + + SurfaceTreeID tree = next_surface_tree_id(); + surface_trees_.push_back(tree); + auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); + std::vector h_tlas_boxes; + std::vector h_surface_instances; + h_tlas_boxes.reserve(volume_surfaces.size()); + h_surface_instances.reserve(volume_surfaces.size()); + + for (const auto &surf : volume_surfaces) { + auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); + const double max_parent_bbox_bump = std::max(bounding_box_bump(mesh_manager, forward_parent), + bounding_box_bump(mesh_manager, reverse_parent)); + + if (!surface_to_blas_map_.count(surf)) { + surface_to_blas_map_[surf] = register_surface(mesh_manager, surf, max_parent_bbox_bump); + } + + CuBQLSurfaceBLAS& surface_blas = surface_to_blas_map_.at(surf); + + // Store BLAS bounding boxes to build TLAS + const auto surface_bounding_box = mesh_manager->surface_bounding_box(surf); + cuBQL::box3d surface_bounds_dp; + surface_bounds_dp.lower = cuBQL::vec3d(surface_bounding_box.min_x, + surface_bounding_box.min_y, + surface_bounding_box.min_z); + surface_bounds_dp.upper = cuBQL::vec3d(surface_bounding_box.max_x, + surface_bounding_box.max_y, + surface_bounding_box.max_z); + + const cuBQL::vec3d bump(max_parent_bbox_bump); + surface_bounds_dp.lower = surface_bounds_dp.lower - bump; + surface_bounds_dp.upper = surface_bounds_dp.upper + bump; + cuBQL::box3f surface_bounds(surface_bounds_dp); + + CuBQLVolumeTLAS::SurfaceInstanceDD surface_instance; + surface_instance.surface_blas = surface_blas.get_device_data(); + + // Sense setting for each surface instance in the TLAS + if (volume_id == forward_parent) { + surface_instance.reverse_sense = false; + } else if (volume_id == reverse_parent) { + surface_instance.reverse_sense = true; + } else { + fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); + } + + h_tlas_boxes.push_back(surface_bounds); + h_surface_instances.push_back(surface_instance); + } + + if (h_surface_instances.empty()) { + fatal_error("Volume {} has no surfaces; cannot build cuBQL surface tree", volume_id); + } + + auto* d_tlas_boxes = static_cast + (omp_target_alloc(h_tlas_boxes.size() * sizeof(cuBQL::box3f), context_.gpuID)); + omp_target_memcpy(d_tlas_boxes, + h_tlas_boxes.data(), + h_tlas_boxes.size() * sizeof(cuBQL::box3f), + 0, + 0, + context_.gpuID, + context_.hostID); + + auto* d_surface_instances = static_cast + (omp_target_alloc(h_surface_instances.size() * sizeof(CuBQLVolumeTLAS::SurfaceInstanceDD), context_.gpuID)); + omp_target_memcpy(d_surface_instances, + h_surface_instances.data(), + h_surface_instances.size() * sizeof(CuBQLVolumeTLAS::SurfaceInstanceDD), + 0, + 0, + context_.gpuID, + context_.hostID); + + cuBQL::BuildConfig tlasBuildParams; + tlasBuildParams.makeLeafThreshold = 1; + tlasBuildParams.maxAllowedLeafSize = 1; + + CuBQLVolumeTLAS volume_tlas; + volume_tlas.volume_id = volume_id; // store meshid in the TLAS object for easier mapping between the two + volume_tlas.num_surface_instances = static_cast(h_surface_instances.size()); + volume_tlas.gpu_id = context_.gpuID; + volume_tlas.d_surface_instances = d_surface_instances; + cuBQL::build_omp_target(volume_tlas.bvh, + d_tlas_boxes, + volume_tlas.num_surface_instances, + tlasBuildParams, + context_.gpuID); + + omp_target_free(d_tlas_boxes, context_.gpuID); + + // Still required for lifetime and scalar calls which need to resolve TreeID->volume_tlas on CPU side. + auto result = tree_to_volume_tlas_.emplace(tree, std::move(volume_tlas)); + auto it = result.first; + + // Keep a dense host-side MeshID -> TLAS device-data table for prepared queries. + // The TLAS object in tree_to_volume_tlas_ owns the device allocations; this table + // only stores lightweight DD views indexed by volume ID. Upload to device once in + // init(), unless a volume is registered after initialization. + + const auto volume_index = static_cast(volume_id); + if (volume_index >= volume_to_tlas_.size()) { + volume_to_tlas_.resize(volume_index + 1); + } + volume_to_tlas_[volume_index] = it->second.get_device_data(); + + if (initialized_) { + upload_volume_to_tlas_table_(); + } + + return tree; +} + +TreeID +CuBQLRayTracer::create_element_tree(const std::shared_ptr&, + MeshID) +{ + warning("Element trees not currently supported with cuBQL ray tracer"); + return TREE_NONE; +} + +void CuBQLRayTracer::create_global_surface_tree() +{ + warning("Global surface trees not currently supported with cuBQL ray tracer"); +} + +void CuBQLRayTracer::create_global_element_tree() +{ + warning("Global element trees not currently supported with cuBQL ray tracer"); +} + +MeshID CuBQLRayTracer::find_element(const Position&) const +{ + fatal_error("Element queries not currently supported with cuBQL ray tracer"); + return ID_NONE; +} + +MeshID CuBQLRayTracer::find_element(TreeID, const Position&) const +{ + fatal_error("Element queries not currently supported with cuBQL ray tracer"); + return ID_NONE; +} + +bool CuBQLRayTracer::point_in_volume(TreeID tree, + const Position& point, + const Direction* direction, + const std::vector* exclude_primitives) const +{ + const auto& context = context_; + const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); + + // Use provided direction or if Direction == nulptr use default direction + Direction directionUsed = (direction != nullptr) ? Direction{direction->x, direction->y, direction->z} + : Direction{1. / std::sqrt(2.0), 1. / std::sqrt(2.0), 0.0}; + + + CuBQLRay ray; + ray.origin = cuBQL::vec3d(point.x, point.y, point.z); + ray.direction = cuBQL::vec3d(directionUsed.x, directionUsed.y, directionUsed.z); + ray.tMin = 0.0; + ray.tMax = INFTY; + + CuBQLSurfaceHit surface_hit; + + // TODO - Maybe we can come up with a better name for this + intersect_surface_tree_scalar(context, volume_tlas, ray, surface_hit, HitOrientation::ANY, exclude_primitives); + + // if the ray hit nothing the point must be outside the volume + if (surface_hit.primitive == ID_NONE) return false; + + return surface_hit.piv == INSIDE; +} + +std::pair +CuBQLRayTracer::ray_fire(TreeID tree, + const Position& origin, + const Direction& direction, + const double tmax, + HitOrientation hitOrientation, + std::vector* const exclude_primitives) +{ + const auto& context = context_; + const CuBQLVolumeTLAS& volume_tlas = tree_to_volume_tlas_.at(tree); + + CuBQLRay ray; + ray.origin = cuBQL::vec3d(origin.x, origin.y, origin.z); + ray.direction = cuBQL::vec3d(direction.x, direction.y, direction.z); + ray.tMin = 0.0; + ray.tMax = tmax; + + CuBQLSurfaceHit surface_hit; + + // TODO - Maybe we can come up with a better name for this + intersect_surface_tree_scalar(context, volume_tlas, ray, surface_hit, hitOrientation, exclude_primitives); + + if (surface_hit.primitive == ID_NONE) { + return {INFTY, ID_NONE}; + } + + if (exclude_primitives) { + exclude_primitives->push_back(surface_hit.primitive); + } + + return {surface_hit.distance, surface_hit.surface}; +} + +void +CuBQLRayTracer::ray_fire_batch(const CuBQLRay* d_rays, + CuBQLSurfaceHit* d_hits, + std::size_t num_rays, + HitOrientation orientation) +{ + if (num_rays == 0) return; + + if (!d_volume_to_tlas_) { + fatal_error("cuBQL volume TLAS lookup table has not been uploaded"); + } + + intersect_surface_tree_batch(context_, + d_volume_to_tlas_, + d_rays, + d_hits, + num_rays, + orientation); +} + +std::pair +CuBQLRayTracer::closest(TreeID, const Position&) +{ + fatal_error("Closest queries not currently supported with cuBQL ray tracer"); + return {INFTY, ID_NONE}; +} + +bool CuBQLRayTracer::occluded(TreeID, + const Position&, + const Direction&, + double&) const +{ + fatal_error("Occlusion queries not currently supported with cuBQL ray tracer"); + return false; +} + +} // namespace xdg diff --git a/src/cuBQL/triangles.cpp b/src/cuBQL/triangles.cpp new file mode 100644 index 00000000..73cf1fe3 --- /dev/null +++ b/src/cuBQL/triangles.cpp @@ -0,0 +1,52 @@ +#include "xdg/cuBQL/triangles.h" + +#include + +namespace xdg { + +void CuBQLSurfaceMesh::release() +{ + if (d_vertices) { + omp_target_free(d_vertices, gpu_id); + d_vertices = nullptr; + } + if (d_indices) { + omp_target_free(d_indices, gpu_id); + d_indices = nullptr; + } + if (d_primitive_refs) { + omp_target_free(d_primitive_refs, gpu_id); + d_primitive_refs = nullptr; + } +} + +void CuBQLSurfaceBLAS::release() +{ + if (bvh.primIDs) { + omp_target_free(bvh.primIDs, gpu_id); + bvh.primIDs = nullptr; + } + if (bvh.nodes) { + omp_target_free(bvh.nodes, gpu_id); + bvh.nodes = nullptr; + } + mesh.release(); +} + +void CuBQLVolumeTLAS::release() +{ + if (bvh.primIDs) { + omp_target_free(bvh.primIDs, gpu_id); + bvh.primIDs = nullptr; + } + if (bvh.nodes) { + omp_target_free(bvh.nodes, gpu_id); + bvh.nodes = nullptr; + } + if (d_surface_instances) { + omp_target_free(d_surface_instances, gpu_id); + d_surface_instances = nullptr; + } +} + +} // namespace xdg diff --git a/src/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/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); }; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 21691467..bdb1f7e8 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -54,6 +54,9 @@ foreach(test ${TEST_NAMES}) if (XDG_ENABLE_MOAB) target_link_libraries(${test} PRIVATE MOAB) endif() + if (XDG_ENABLE_CUBQL) + target_link_libraries(${test} PRIVATE $) + endif() set_target_properties(${test} PROPERTIES BUILD_RPATH "$") catch_discover_tests(${test} diff --git a/tests/test_point_in_volume.cpp b/tests/test_point_in_volume.cpp index ae34e823..65874f9f 100644 --- a/tests/test_point_in_volume.cpp +++ b/tests/test_point_in_volume.cpp @@ -16,7 +16,8 @@ using namespace xdg::test; TEMPLATE_TEST_CASE("Point-in-volume on MeshMock", "[piv][mock]", Embree_Raytracer, - GPRT_Raytracer) + GPRT_Raytracer, + CuBQL_Raytracer) { constexpr auto rt_backend = TestType::value; diff --git a/tests/test_ray_fire.cpp b/tests/test_ray_fire.cpp index be816c36..70a5353f 100644 --- a/tests/test_ray_fire.cpp +++ b/tests/test_ray_fire.cpp @@ -18,7 +18,8 @@ using namespace xdg::test; TEMPLATE_TEST_CASE("Ray Fire on MeshMock (per-backend sections)", "[rayfire][mock]", Embree_Raytracer, - GPRT_Raytracer) + GPRT_Raytracer, + CuBQL_Raytracer) { // Generate one test run per enabled backend constexpr auto rt_backend = TestType::value; diff --git a/tests/test_ray_tracer_cross_check.cpp b/tests/test_ray_tracer_cross_check.cpp index c75b0273..d2407edf 100644 --- a/tests/test_ray_tracer_cross_check.cpp +++ b/tests/test_ray_tracer_cross_check.cpp @@ -57,7 +57,7 @@ TEST_CASE("Test Pincell RT libraries Cross-Check ray_fire queries", "[moab][rayf } srand48(12345); // set fixed seed for rng - std::vector directions(1000); + std::vector directions(100); for (auto &dir : directions) { dir = rand_dir(); } @@ -86,4 +86,59 @@ TEST_CASE("Test Pincell RT libraries Cross-Check ray_fire queries", "[moab][rayf } } +TEST_CASE("Test Pincell RT libraries Cross-Check point_in_volume queries", "[moab][piv][cross-check]") +{ + const auto rt_cases = make_rt_cases("pincell.h5m"); + if (rt_cases.size() < 2) { + SKIP("Fewer than two ray tracing backends are available; skipping cross-check."); + } + + const std::array directions { + Direction {1.0, 0.0, 0.0}, // axis-aligned x ray + Direction {0.0, 1.0, 0.0}, // axis-aligned y ray + Direction {0.0, 0.0, 1.0}, // axis-aligned z ray + Direction {0.371390676, 0.557086014, 0.742781353} // non-axis ray + }; + + // Points generated by codex with a h5dump to test various edge cases around the pincell geometry + const std::array points { + Position {0.0, 0.0, 0.0}, // pincell center + Position {8.5, 0.25, 0.0}, // inside inner cylinder, r < 9 + Position {9.5, 0.25, 0.0}, // between cylinder radii, 9 < r < 10 + Position {10.5, 0.25, 0.0}, // outside outer cylinder, r > 10 + Position {20.0, 0.25, 0.0}, // inside square cell away from cylinder + Position {24.5, 0.25, 0.0}, // just inside x = 25 square boundary + Position {25.5, 0.25, 0.0}, // between x = 25 and outer x = 27.5 + Position {27.0, 0.25, 0.0}, // just inside outer x = 27.5 boundary + Position {28.0, 0.25, 0.0}, // outside outer x boundary + Position {0.25, 0.0, 19.5}, // just inside positive z = 20 cap + Position {0.25, 0.0, 20.5}, // just outside positive z = 20 cap + Position {0.25, 0.0, -19.5}, // just inside negative z = -20 cap + Position {0.25, 0.0, -20.5}, // just outside negative z = -20 cap + Position {0.25, 24.5, 0.0}, // just inside y = 25 square boundary + Position {0.25, 28.0, 0.0} // outside outer y boundary + }; + + const auto& reference_case = rt_cases.front(); + + for (const auto& volume : reference_case.xdg->mesh_manager()->volumes()) { + for (const auto& point : points) { + for (const auto& direction : directions) { + const auto reference_result = + reference_case.xdg->point_in_volume(volume, point, &direction); + + for (size_t i = 1; i < rt_cases.size(); ++i) { + const auto& candidate = rt_cases[i]; + const auto candidate_result = + candidate.xdg->point_in_volume(volume, point, &direction); + + CAPTURE(volume, point, direction, reference_case.name, candidate.name, + reference_result, candidate_result); + REQUIRE(candidate_result == reference_result); + } + } + } + } +} + // TODO - Add all of the other queries diff --git a/tests/util.h b/tests/util.h index a9841086..96946aa2 100644 --- a/tests/util.h +++ b/tests/util.h @@ -8,7 +8,7 @@ #include "xdg/constants.h" #include "xdg/ray_tracers.h" #include "xdg/mesh_managers.h" -#include "xdg/gprt/vulkan_probe.h" +#include "xdg/available_device_probe.h" namespace xdg::test { @@ -17,7 +17,7 @@ using LibMesh_Interface = std::integral_constant; using GPRT_Raytracer = std::integral_constant; - +using CuBQL_Raytracer = std::integral_constant; } // namespace xdg::test namespace Catch { @@ -49,6 +49,13 @@ inline bool ray_tracer_available(xdg::RTLibrary rt) { #else return false; #endif + + case xdg::RTLibrary::CUBQL: + #ifdef XDG_ENABLE_CUBQL + return system_has_omp_target_device(); + #else + return false; + #endif } return false; @@ -118,5 +125,10 @@ create_raytracer(xdg::RTLibrary rt) { return std::make_shared(); #endif + #ifdef XDG_ENABLE_CUBQL + if (rt == xdg::RTLibrary::CUBQL) + return std::make_shared(); + #endif + return nullptr; } diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 50a0c650..883cbe84 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -38,3 +38,7 @@ if (TARGET overlap-check) ${CMAKE_CURRENT_SOURCE_DIR}/overlap.cpp ) endif() + +if (XDG_ENABLE_CUBQL AND TARGET ray-benchmark) + target_link_libraries(ray-benchmark PRIVATE $) +endif() diff --git a/tools/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); diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp index dea28351..847389d7 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -19,6 +19,11 @@ #include "ray_benchmark.h" +#ifdef XDG_ENABLE_CUBQL +#include +#include "xdg/cuBQL/intersection.h" +#include "xdg/cuBQL/ray_tracer.h" +#endif using namespace xdg; @@ -60,7 +65,7 @@ int main(int argc, char** argv) .default_value("MOAB"); args.add_argument("-rt", "--rt-library") - .help("Ray tracing library to use. Currently implemented: EMBREE") + .help("Ray tracing library to use. Currently implemented: EMBREE, CUBQL") .default_value("EMBREE"); args.add_argument("-l", "--list") @@ -98,6 +103,8 @@ int main(int argc, char** argv) RTLibrary rt_lib; if (rt_str == "EMBREE") { rt_lib = RTLibrary::EMBREE; + } else if (rt_str == "CUBQL") { + rt_lib = RTLibrary::CUBQL; } else { fatal_error("Ray tracing library '{}' is not implemented in this benchmark tool yet", rt_str); } @@ -171,53 +178,127 @@ int main(int argc, char** argv) setup_timer.stop(); + const auto num_faces = mesh_manager->num_volume_faces(volume); + std::size_t num_hits = 0; + if (num_rays < 1) fatal_error("Number of rays must be greater than 0"); + if (rt_lib == RTLibrary::EMBREE) { rt_label += " (" + std::to_string(XDGConfig::config().n_threads()) + " CPU threads)"; - } - const auto num_faces = mesh_manager->num_volume_faces(volume); - - - // Generate random rays from source - generation_timer.start(); - std::vector origins(num_rays); - std::vector directions(num_rays); - - #pragma omp parallel for schedule(runtime) - for (std::size_t i = 0; i < num_rays; ++i) { - std::uint32_t state = seed ^ static_cast(i); - auto sample = tools::benchmark::random_spherical_source(origin.x, - origin.y, - origin.z, - state, - source_radius); - origins[i] = Position(sample.position[0], - sample.position[1], - sample.position[2]); - directions[i] = Direction(sample.direction[0], - sample.direction[1], - sample.direction[2]); - } - generation_timer.stop(); + // Generate random rays from source + generation_timer.start(); + std::vector origins(num_rays); + std::vector directions(num_rays); + + #pragma omp parallel for schedule(runtime) + for (std::size_t i = 0; i < num_rays; ++i) { + std::uint32_t state = seed ^ static_cast(i); + auto sample = tools::benchmark::random_spherical_source(origin.x, + origin.y, + origin.z, + state, + source_radius); + origins[i] = Position(sample.position[0], + sample.position[1], + sample.position[2]); + directions[i] = Direction(sample.direction[0], + sample.direction[1], + sample.direction[2]); + } + generation_timer.stop(); - // Trace rays - trace_timer.start(); + std::vector hit_surfaces(num_rays, ID_NONE); - std::size_t num_hits = 0; + trace_timer.start(); + #pragma omp parallel for schedule(runtime) + for (std::size_t i = 0; i < num_rays; ++i) { + hit_surfaces[i] = xdg->ray_fire(volume, origins[i], directions[i]).second; // just return surface id of hit + } + trace_timer.stop(); - #pragma omp parallel for schedule(runtime) reduction(+:num_hits) - for (std::size_t i = 0; i < num_rays; ++i) { - const auto hit = xdg->ray_fire(volume, origins[i], directions[i]); - if (hit.second != ID_NONE) num_hits++; + // Count hits outside of timing region + for (std::size_t i = 0; i < num_rays; ++i) { + if (hit_surfaces[i] != ID_NONE) num_hits++; + } + } + else if (rt_lib == RTLibrary::CUBQL) { + #ifndef XDG_ENABLE_CUBQL + fatal_error("This build was not compiled with cuBQL support (XDG_ENABLE_CUBQL=OFF)."); + #else + auto rti = std::dynamic_pointer_cast(xdg->ray_tracing_interface()); + + // Generate random rays directly on the target device. + generation_timer.start(); + + const int gpu_id = omp_get_default_device(); + + CuBQLRay* d_rays = static_cast( + omp_target_alloc(num_rays * sizeof(CuBQLRay), gpu_id)); + + if (!d_rays) { + fatal_error("Failed to allocate cuBQL ray buffer"); + } + + const double origin_x = origin.x; + const double origin_y = origin.y; + const double origin_z = origin.z; + + #pragma omp target teams distribute parallel for device(gpu_id) is_device_ptr(d_rays) + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + std::uint32_t state = seed ^ static_cast(ray_id); + + auto sample = tools::benchmark::random_spherical_source(origin_x, + origin_y, + origin_z, + state, + source_radius); + + CuBQLRay ray; + ray.origin = cuBQL::vec3d(sample.position[0], + sample.position[1], + sample.position[2]); + ray.direction = cuBQL::vec3d(sample.direction[0], + sample.direction[1], + sample.direction[2]); + ray.tMin = 0.0; + ray.tMax = INFTY; + ray.volume = volume; + + d_rays[ray_id] = ray; + } + + CuBQLSurfaceHit* d_hits = static_cast( + omp_target_alloc(num_rays * sizeof(CuBQLSurfaceHit), gpu_id)); + + if (!d_hits) { + omp_target_free(d_rays, gpu_id); + fatal_error("Failed to allocate cuBQL hit buffer"); + } + + generation_timer.stop(); + + // Trace rays and count hits on the target device. + trace_timer.start(); + rti->ray_fire_batch(d_rays, d_hits, num_rays); + trace_timer.stop(); + + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(d_hits) reduction(+:num_hits) + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + if (d_hits[ray_id].primitive != ID_NONE) num_hits++; + } + + omp_target_free(d_hits, gpu_id); + omp_target_free(d_rays, gpu_id); + #endif } - - trace_timer.stop(); const std::size_t num_misses = num_rays - num_hits; const double hit_fraction = num_rays > 0 ? static_cast(num_hits) / static_cast(num_rays) : 0.0; + const double generation_time = generation_timer.elapsed(); const double trace_time = trace_timer.elapsed(); const double end_to_end_time = generation_time + trace_time; 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 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