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 fb762521..d05306a1 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,20 @@ dbl_deviceCode endif() +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) list(APPEND xdg_sources src/libmesh/mesh_manager.cpp @@ -281,7 +312,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 +355,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..6a6c8111 --- /dev/null +++ b/include/xdg/cuBQL/cuBQL_backend.h @@ -0,0 +1,123 @@ +#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_backend { + +struct Context { + int gpuID {0}; + int hostID {omp_get_initial_device()}; +}; + +template +struct AutoUploadArray { + AutoUploadArray() = default; + + AutoUploadArray(Context context, const InputT* input, std::size_t count) + : context(context), count(count) + { + if (count == 0) return; + + elements = static_cast(omp_target_alloc(count * sizeof(T), context.gpuID)); + if (!elements) { + fatal_error("Failed to allocate cuBQL device array."); + } + needsFree = true; + + int result = 0; + if constexpr (std::is_same_v) { + result = omp_target_memcpy(elements, + input, + count * sizeof(T), + 0, + 0, + context.gpuID, + context.hostID); + } else { + std::vector converted(count); + for (std::size_t i = 0; i < count; ++i) { + converted[i] = T(input[i]); + } + + result = omp_target_memcpy(elements, + converted.data(), + count * sizeof(T), + 0, + 0, + context.gpuID, + context.hostID); + } + + if (result != 0) { + reset(); + fatal_error("Failed to upload cuBQL device array."); + } + } + + ~AutoUploadArray() + { + reset(); + } + + AutoUploadArray(const AutoUploadArray&) = delete; + AutoUploadArray& operator=(const AutoUploadArray&) = delete; + + AutoUploadArray(AutoUploadArray&& other) noexcept + { + *this = std::move(other); + } + + AutoUploadArray& operator=(AutoUploadArray&& other) noexcept + { + if (this != &other) { + reset(); + + context = other.context; + elements = other.elements; + count = other.count; + needsFree = other.needsFree; + + other.elements = nullptr; + other.count = 0; + other.needsFree = false; + } + + return *this; + } + + void reset() + { + if (needsFree && elements) { + omp_target_free(elements, context.gpuID); + } + + elements = nullptr; + count = 0; + needsFree = false; + } + + T* get() const { return elements; } + + T* elements {nullptr}; + std::size_t count {0}; + bool needsFree {false}; + Context context {}; +}; + +} // namespace xdg_cubql_backend + +#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..6c8a6501 --- /dev/null +++ b/include/xdg/cuBQL/ray_tracer.h @@ -0,0 +1,97 @@ +#ifndef _XDG_CUBQL_RAY_TRACING_INTERFACE_H +#define _XDG_CUBQL_RAY_TRACING_INTERFACE_H + +#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 { + +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; + +private: + xdg_cubql_backend::Context context_; + + struct CuBQLSurfaceBLAS { + MeshID surface {ID_NONE}; + MeshID forward_parent {ID_NONE}; + MeshID reverse_parent {ID_NONE}; + cuBQL::bvh3d bvh; + + xdg_cubql_backend::AutoUploadArray vertices; + xdg_cubql_backend::AutoUploadArray indices; + xdg_cubql_backend::AutoUploadArray primitive_refs; + + 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::unordered_map surface_tree_to_volume_; +}; + +} // 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/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp new file mode 100644 index 00000000..3660fb44 --- /dev/null +++ b/src/cuBQL/ray_tracer.cpp @@ -0,0 +1,383 @@ +#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" +#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 { + +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()) { + 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() +{ + 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); + } + + } +} + +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& 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 + + 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& 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) { + 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]); + } + + // 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); + } + + xdg_cubql_backend::AutoUploadArray d_vertices(context, h_vertices.data(), h_vertices.size()); + xdg_cubql_backend::AutoUploadArray d_indices(context, h_indices.data(), h_indices.size()); + xdg_cubql_backend::AutoUploadArray d_primitive_refs(context, h_primitive_refs.data(), h_primitive_refs.size()); + + // 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? + auto* vertices_ptr = d_vertices.get(); + auto* indices_ptr = d_indices.get(); + #pragma omp target device(gpu_id) is_device_ptr(vertices_ptr, indices_ptr, d_aabbs) + #pragma omp teams distribute parallel for + for (uint32_t primID = 0; primID < num_primitives; ++primID) { + cuBQL::vec3i indices = indices_ptr[primID]; + + cuBQL::vec3d A = vertices_ptr[indices.x]; + cuBQL::vec3d B = vertices_ptr[indices.y]; + cuBQL::vec3d C = vertices_ptr[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); + + 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); + } + + 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.vertices = std::move(d_vertices); + surface_bvh.indices = std::move(d_indices); + surface_bvh.primitive_refs = std::move(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(std::move(surface_bvh)); + } + + + 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, + 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 tree, + const Position& origin, + const Direction& direction, + const double tmax, + HitOrientation hitOrientation, + std::vector* const exclude_primitives) +{ + const auto& context = context_; + const int gpu_id = context.gpuID; + + int exclude_count = 0; + xdg_cubql_backend::AutoUploadArray d_exclude_primitives; + if (exclude_primitives && !exclude_primitives->empty()) { + exclude_count = static_cast(exclude_primitives->size()); + d_exclude_primitives = xdg_cubql_backend::AutoUploadArray(context, + exclude_primitives->data(), + exclude_count); + } + + 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)); + + 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); + + cuBQL::bvh3d bvh = surface_bvh.bvh; + const cuBQL::vec3d* d_vertices = surface_bvh.vertices.get(); + const cuBQL::vec3i* d_indices = surface_bvh.indices.get(); + const MeshID* d_primitive_refs = surface_bvh.primitive_refs.get(); + const MeshID* d_exclude_primitive_ptr = d_exclude_primitives.get(); + const bool reverse_sense = query_volume == surface_bvh.reverse_parent; + + 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_vertices, d_indices, d_primitive_refs, d_exclude_primitive_ptr, d_surface_hit) + { + 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]; + + for (int i = 0; i < exclude_count; ++i) { + if (d_exclude_primitive_ptr[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] + }; + + 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; + } else if (orientation == static_cast(HitOrientation::ENTERING)) { + culled = normal_dot_direction >= 0.0; + } + if (culled) { + return ray.tMax; + } + + // 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; + ray.tMax = intersection.t; + } + + return ray.tMax; + }; + + // 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_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 (closest_surface == ID_NONE) { + return {INFTY, ID_NONE}; + } + + if (exclude_primitives) { + exclude_primitives->push_back(closest_primitive); + } + + return {closest_distance, closest_surface}; +} + +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/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 371a9a1c..83816e6c 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); }; @@ -325,4 +339,4 @@ double XDG::measure_volume_area(MeshID volume) const return area; } -} // namespace xdg \ No newline at end of file +} // namespace xdg diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 222626f2..0009dbfb 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -53,6 +53,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_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..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(); } 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/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/GPRT b/vendor/GPRT index f1e95e41..7b55a1ae 160000 --- a/vendor/GPRT +++ b/vendor/GPRT @@ -1 +1 @@ -Subproject commit f1e95e4188cde591547d6b4a33a70bf2afaeec59 +Subproject commit 7b55a1ae432b226be18592f6f2e7f8980c8df9ce 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