diff --git a/.gitmodules b/.gitmodules index 77e2a79a..dd0e33bd 100644 --- a/.gitmodules +++ b/.gitmodules @@ -20,3 +20,6 @@ [submodule "vendor/GPRT"] path = vendor/GPRT url = https://github.com/gprt-org/GPRT.git +[submodule "vendor/cuBQL"] + path = vendor/cuBQL + url = https://github.com/NVIDIA/cuBQL diff --git a/CMakeLists.txt b/CMakeLists.txt index 0c185556..62124d46 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,6 +7,7 @@ option(XDG_ENABLE_LIBMESH "Enable support for the libMesh mesh library" OFF option(XDG_LINK_MPI "Link with MPI (for dependency compatibility)" OFF) option(XDG_ENABLE_EMBREE "Enable support for the Embree ray tracing library" ON) option(XDG_ENABLE_GPRT "Enable support for the GPRT ray tracing library" OFF) +option(XDG_ENABLE_CUBQL "Enable support for the cuBQL ray tracing library" OFF) option(XDG_BUILD_TESTS "Enable C++ unit testing" ON) option(XDG_BUILD_TOOLS "Enable tools and miniapps" ON) @@ -20,6 +21,10 @@ if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose build type" FORCE) endif() +if(DEFINED XDG_CMAKE_PRESET) + message(STATUS "XDG CMake preset: ${XDG_CMAKE_PRESET}") +endif() + # Compiler options (things in this section may not be platform-portable) set(CMAKE_POSITION_INDEPENDENT_CODE ON) @@ -111,6 +116,12 @@ if(XDG_ENABLE_GPRT) ) endif() +if(XDG_ENABLE_CUBQL) + list(APPEND VENDOR_PATHS + vendor/cuBQL + ) +endif() + if(GIT_FOUND AND EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/.git") option(XDG_GIT_SUBMODULE "Check submodules during build" ON) if(XDG_GIT_SUBMODULE) @@ -175,11 +186,12 @@ if (NOT XDG_ENABLE_MOAB AND NOT XDG_ENABLE_LIBMESH) endif() # Ensure at least one ray tracing backend is enabled -if (NOT XDG_ENABLE_EMBREE AND NOT XDG_ENABLE_GPRT) +if (NOT XDG_ENABLE_EMBREE AND NOT XDG_ENABLE_GPRT AND NOT XDG_ENABLE_CUBQL) message(FATAL_ERROR "No ray tracing backend enabled. Enable at least one of:\n" " -DXDG_ENABLE_EMBREE=ON\n" - " -DXDG_ENABLE_GPRT=ON") + " -DXDG_ENABLE_GPRT=ON\n" + " -DXDG_ENABLE_CUBQL=ON") endif() # GPRT @@ -188,6 +200,11 @@ if (XDG_ENABLE_GPRT) add_subdirectory(vendor/GPRT) endif() +if (XDG_ENABLE_CUBQL) + set(CUBQL_OMP ON CACHE BOOL "Build cuBQL with OpenMP target offload" FORCE) + add_subdirectory(vendor/cuBQL) +endif() + list(APPEND xdg_sources src/geometry/measure.cpp src/geometry/plucker.cpp @@ -225,6 +242,23 @@ dbl_deviceCode endif() +if (XDG_ENABLE_CUBQL) +list(APPEND xdg_sources +src/cuBQL/triangles.cpp +src/cuBQL/intersection.cpp +src/cuBQL/ray_tracer.cpp +) + +# We need a precompile definition to switch to using the cuBQL math types in the shared +# plucker intersection code. The compile definition is used in dp__math.h +set_source_files_properties( + src/cuBQL/ray_tracer.cpp + src/cuBQL/intersection.cpp + PROPERTIES COMPILE_DEFINITIONS XDG_DP_MATH_CUBQL +) + +endif() + if (XDG_ENABLE_LIBMESH) list(APPEND xdg_sources src/libmesh/mesh_manager.cpp @@ -285,7 +319,12 @@ if (${CMAKE_BUILD_TYPE} MATCHES "Debug") endif() # attempt to find OpenMP and include it if found -find_package(OpenMP) +if (XDG_ENABLE_CUBQL) + find_package(OpenMP REQUIRED) +else() + find_package(OpenMP) +endif() + if (OpenMP_CXX_FOUND) target_link_libraries(xdg PRIVATE OpenMP::OpenMP_CXX) target_compile_definitions(xdg PRIVATE XDG_HAVE_OPENMP) @@ -323,6 +362,18 @@ if (XDG_ENABLE_GPRT) target_link_options(xdg PRIVATE -Wl,--unresolved-symbols=ignore-in-shared-libs) endif() +if (XDG_ENABLE_CUBQL) + target_compile_definitions(xdg PUBLIC XDG_ENABLE_CUBQL) + target_link_libraries(xdg PRIVATE $) + # TODO: Stop relying on LD_LIBRARY_PATH for LLVM OpenMP offload runtimes. + # As a temporary measure whilst figuring out our way around cuBQL this is + # okay but in the long run we should aim for a more robust solutions here. + # Clang injects libomptarget when offload flags are supplied by presets, but + # CMake does not currently add that compiler runtime directory to RPATH. + # Add targeted BUILD_RPATH handling for libomptarget/libomp, and decide + # whether install RPATH should remain environment-module based or be opt-in. +endif() + target_link_libraries(xdg PRIVATE fmt::fmt) # ========================== diff --git a/CMakePresets.json b/CMakePresets.json new file mode 100644 index 00000000..ad89fa82 --- /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_89" + } + }, + { + "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 caf43fd3..ba652afb 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -54,7 +54,8 @@ enum class MeshLibrary { // Ray Tracing library identifier enum class RTLibrary { EMBREE, - GPRT + GPRT, + CUBQL }; static const std::map MESH_LIB_TO_STR = @@ -67,7 +68,8 @@ static const std::map MESH_LIB_TO_STR = static const std::map RT_LIB_TO_STR = { {RTLibrary::EMBREE, "EMBREE"}, - {RTLibrary::GPRT, "GPRT"} + {RTLibrary::GPRT, "GPRT"}, + {RTLibrary::CUBQL, "CUBQL"} }; // Mesh identifer type @@ -181,4 +183,4 @@ struct formatter : fmt::formatter { } -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/cuBQL/cuBQL_backend.h b/include/xdg/cuBQL/cuBQL_backend.h new file mode 100644 index 00000000..b348802f --- /dev/null +++ b/include/xdg/cuBQL/cuBQL_backend.h @@ -0,0 +1,30 @@ +#ifndef _XDG_CUBQL_BACKEND_H +#define _XDG_CUBQL_BACKEND_H + +#include +#include +#include +#include + +#include + +// Guards to prevent CUDA headers from being included in host code, which causes +// failed compilation with LLVM-clang. +#if defined(__CUDA_ARCH__) && !defined(__CUDACC__) +#undef __CUDA_ARCH__ +#endif + +#include "xdg/error.h" + +namespace xdg::cubql { + +struct Context { + int gpuID {0}; + int hostID {omp_get_initial_device()}; +}; + + + +} // namespace xdg::cubql + +#endif // include guard diff --git a/include/xdg/cuBQL/intersection.h b/include/xdg/cuBQL/intersection.h new file mode 100644 index 00000000..9a941a4f --- /dev/null +++ b/include/xdg/cuBQL/intersection.h @@ -0,0 +1,82 @@ +#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/device_ray.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 +}; + +// 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}; + MeshID next_volume {ID_NONE}; + SurfaceBoundaryCondition boundary_condition {SurfaceBoundaryCondition::UNSET}; + cuBQL::vec3d normal {0.0}; + + 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 CuBQLVolumeGroup& volume_group, + const CuBQLRay& ray, + CuBQLSurfaceHit& hit, + HitOrientation hit_orientation, + const std::vector* exclude_primitives); + + +void +intersect_surface_tree_batch(const cubql::Context& context, + const CuBQLVolumeGroup::DD* d_volume_to_group, + XDGRayHit* d_ray_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..ebd4bb88 --- /dev/null +++ b/include/xdg/cuBQL/ray_tracer.h @@ -0,0 +1,95 @@ +#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 XDGRayHitBuffer& ray_hits, + HitOrientation hit_orientation = HitOrientation::EXITING) const override; + + XDGRayHitBuffer allocate_ray_hits(std::size_t count) const override; + + void free_ray_hits(XDGRayHitBuffer& ray_hits) const override; + + std::pair closest(TreeID tree, + const Position& origin) override; + + bool occluded(TreeID tree, + const Position& origin, + const Direction& direction, + double& dist) const override; + +private: + CuBQLSurfaceMesh + register_surface(const std::shared_ptr& mesh_manager, + MeshID surface_id); + + void upload_volume_to_group_table_(); + + cubql::Context context_; + + std::unordered_map tree_to_volume_group_; + std::unordered_map surface_to_mesh_; + + std::vector volume_to_group_; + CuBQLVolumeGroup::DD* d_volume_to_group_ {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..07016816 --- /dev/null +++ b/include/xdg/cuBQL/triangles.h @@ -0,0 +1,114 @@ +#ifndef _XDG_CUBQL_TRIANGLES_H +#define _XDG_CUBQL_TRIANGLES_H + +#include + +// Guards to prevent CUDA headers from being included in host code, which causes +// failed compilation with LLVM-clang. +#if defined(__CUDA_ARCH__) && !defined(__CUDACC__) +#undef __CUDA_ARCH__ +#endif + +#include "cuBQL/bvh.h" +#include "cuBQL/math/vec.h" +#include "xdg/constants.h" +#include "xdg/cuBQL/cuBQL_backend.h" + +namespace xdg { + +/** + Owns the primitive buffers for one topological surface. The nested DD type is + the compact, non-owning device-data view used in OpenMP target regions. +*/ +struct CuBQLSurfaceMesh { + struct DD { + // Topological metadata + MeshID surface_id {ID_NONE}; + double max_parent_volume_bump {0.0}; + + // Geometric data + const cuBQL::vec3d* vertices {nullptr}; + const cuBQL::vec3i* indices {nullptr}; + const MeshID* primitive_ids {nullptr}; + }; + + // Topological metadata + MeshID surface_id {ID_NONE}; + double max_parent_volume_bump {0.0}; + + // Device buffers for primitive data + cuBQL::vec3d* d_vertices {nullptr}; + cuBQL::vec3i* d_indices {nullptr}; + MeshID* d_primitive_ids {nullptr}; + + std::uint32_t num_vertices {0}; + std::uint32_t num_primitives {0}; + int gpu_id {0}; + + // Return the non-owning view used by device traversal and intersection code. + DD get_device_data() const + { + return { + surface_id, + max_parent_volume_bump, + d_vertices, + d_indices, + d_primitive_ids + }; + } + + void release(); +}; + +/** + Owns the flattened cuBQL BVH for one topological volume. The BVH is built over + all primitives belonging to the volume's surfaces, while PrimRef maps each + BVH primitive back to its surface and surface-local primitive. The nested DD + type is the compact, non-owning device-data view used during traversal. +*/ +struct CuBQLVolumeGroup { + struct SurfaceDD { + CuBQLSurfaceMesh::DD mesh; + + bool reverse_sense {false}; + MeshID next_volume {ID_NONE}; + SurfaceBoundaryCondition boundary_condition {UNSET}; + }; + + // Identifies a primitive within the volume group's local surface array. + // TODO - Should this exist outside of volume group as an indpendent struct? + // TODO - can we think of a better name to distinguish between CuBQLSurfaceMesh::primitive_ids and this struct? + struct PrimRef { + std::uint32_t surface_index {0}; + std::uint32_t primitive_index {0}; + }; + + struct DD { + const SurfaceDD* surfaces {nullptr}; + const PrimRef* prim_refs {nullptr}; + cuBQL::bvh3f bvh; + }; + + cuBQL::bvh3f bvh; + SurfaceDD* d_surfaces {nullptr}; + PrimRef* d_prim_refs {nullptr}; + + std::uint32_t num_surfaces {0}; + std::uint32_t num_primitives {0}; + int gpu_id {0}; + + // Return the non-owning view used by device traversal and intersection code. + DD get_device_data() const + { + return {d_surfaces, d_prim_refs, bvh}; + } + + void release(); +}; + +// Future acceleration structure over instances of volume groups. +struct CuBQLInstanceGroup; + +} // namespace xdg + +#endif // include guard diff --git a/include/xdg/device_ray.h b/include/xdg/device_ray.h new file mode 100644 index 00000000..8b113aa9 --- /dev/null +++ b/include/xdg/device_ray.h @@ -0,0 +1,35 @@ +#ifndef _XDG_DEVICE_RAY_H +#define _XDG_DEVICE_RAY_H + +#include +#include + +namespace xdg { + +struct XDGRayHit { + double origin[3]; + double direction[3]; + double t_min; + double t_max; + std::int32_t volume; + std::int32_t last_hit_primitive {-1}; + + double distance; + std::int32_t surface; + std::int32_t primitive; + std::int32_t point_in_volume; + std::int32_t next_volume; + std::int32_t boundary_condition; + double normal[3]; +}; + +// Light wrapper for count and device id associated with pointer +struct XDGRayHitBuffer { + XDGRayHit* data {nullptr}; + std::size_t count {0}; + int device_id {-1}; +}; + +} // 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/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h index e0cf3cb7..d593cbb9 100644 --- a/include/xdg/ray_tracing_interface.h +++ b/include/xdg/ray_tracing_interface.h @@ -1,11 +1,13 @@ #ifndef _XDG_RAY_TRACING_INTERFACE_H #define _XDG_RAY_TRACING_INTERFACE_H +#include #include #include #include #include "xdg/constants.h" +#include "xdg/device_ray.h" #include "xdg/mesh_manager_interface.h" #include "xdg/primitive_ref.h" #include "xdg/geometry_data.h" @@ -85,6 +87,13 @@ class RayTracer { HitOrientation orientation = HitOrientation::EXITING, std::vector* const exclude_primitives = nullptr) = 0; + virtual XDGRayHitBuffer allocate_ray_hits(std::size_t count) const; + + virtual void free_ray_hits(XDGRayHitBuffer& ray_hits) const; + + virtual void ray_fire_batch(const XDGRayHitBuffer& ray_hits, + HitOrientation hit_orientation = HitOrientation::EXITING) const; + /** * @brief Finds the element containing a given point using the global element tree. * @@ -149,4 +158,4 @@ class RayTracer { } // namespace xdg -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/shared_enums.h b/include/xdg/shared_enums.h index f6198f60..aa0f6fc1 100644 --- a/include/xdg/shared_enums.h +++ b/include/xdg/shared_enums.h @@ -14,6 +14,13 @@ namespace xdg { ENTERING = 1, }; + enum SurfaceBoundaryCondition : int { + UNSET = -1, + TRANSMISSION = 0, // Cross into the next volume + VACUUM = 1, // Kill particle on crossing + REFLECTIVE = 2 // Reflect particle using the surface normal + }; + } -#endif // XDG_SHARED_ENUMS_H \ No newline at end of file +#endif // XDG_SHARED_ENUMS_H diff --git a/include/xdg/util/rng.h b/include/xdg/util/rng.h index 539fba9c..e2069327 100644 --- a/include/xdg/util/rng.h +++ b/include/xdg/util/rng.h @@ -1,6 +1,8 @@ #ifndef XDG_UTIL_RNG_H #define XDG_UTIL_RNG_H +#include +#include #include namespace xdg { @@ -20,6 +22,39 @@ inline double rand_double(double min=0.0, double max=1.0) return min + (max - min) * dis(gen); } +#ifdef _OPENMP +#pragma omp declare target +#endif + +inline double lcg_rand01(std::uint32_t& state) +{ + state = state * 1664525u + 1013904223u; + return static_cast(state) * (1.0 / 4294967296.0); +} + +inline void random_unit_direction_lcg(std::uint32_t& state, + double direction[3]) +{ + double x1; + double x2; + double s; + + do { + x1 = lcg_rand01(state) * 2.0 - 1.0; + x2 = lcg_rand01(state) * 2.0 - 1.0; + s = x1 * x1 + x2 * x2; + } while (s <= 0.0 || s >= 1.0); + + const double t = 2.0 * std::sqrt(1.0 - s); + direction[0] = x1 * t; + direction[1] = x2 * t; + direction[2] = 1.0 - 2.0 * s; +} + +#ifdef _OPENMP +#pragma omp end declare target +#endif + } // namespace xdg #endif // XDG_UTIL_RNG_H diff --git a/include/xdg/xdg.h b/include/xdg/xdg.h index 82c1f8f2..bcae8a30 100644 --- a/include/xdg/xdg.h +++ b/include/xdg/xdg.h @@ -1,6 +1,7 @@ #ifndef _XDG_INTERFACE_H #define _XDG_INTERFACE_H +#include #include #include @@ -75,6 +76,23 @@ std::pair ray_fire(MeshID volume, HitOrientation orientation = HitOrientation::EXITING, std::vector* const exclude_primitives = nullptr) const; +//! Allocates a backend-owned device buffer for batch ray-fire records. +//! @param count Number of XDGRayHit records to allocate. +//! @return Device buffer handle containing the pointer, record count, and device id. +//! Release the buffer with free_ray_hits(). +XDGRayHitBuffer allocate_ray_hits(std::size_t count) const; + +//! Releases a device ray-hit buffer allocated by allocate_ray_hits(). +//! @param ray_hits Buffer handle to release. The handle is cleared after release. +void free_ray_hits(XDGRayHitBuffer& ray_hits) const; + +//! Fires all rays stored in a device ray-hit buffer using the selected backend. +//! @param ray_hits Device buffer whose XDGRayHit records have been populated by the caller. +//! Hit result fields are written back into the same records. +//! @param hit_orientation Orientation filter applied to every ray in the batch. +void ray_fire_batch(const XDGRayHitBuffer& ray_hits, + HitOrientation hit_orientation = HitOrientation::EXITING) const; + std::pair closest(MeshID volume, const Position& origin) const; diff --git a/src/cuBQL/intersection.cpp b/src/cuBQL/intersection.cpp new file mode 100644 index 00000000..e469803f --- /dev/null +++ b/src/cuBQL/intersection.cpp @@ -0,0 +1,238 @@ +#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 flattened +// volume group. +#pragma omp declare target +static inline float reject_candidate(const cuBQL::ray3f& traversal_ray) +{ + return traversal_ray.tMax; +} + +static inline void intersect_surface_tree(CuBQLVolumeGroup::DD volume_group, + CuBQLRay intersection_ray, + CuBQLSurfaceHit* hit, + int orientation, + MeshID last_hit_primitive, + 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); + + auto intersect_prim = [=, &traversal_ray] + (std::uint32_t bvh_primitive_index) -> float + { + const auto ref = volume_group.prim_refs[bvh_primitive_index]; + const auto surface = volume_group.surfaces[ref.surface_index]; + const auto mesh = surface.mesh; + const auto local_index = ref.primitive_index; + const MeshID primitive_id = mesh.primitive_ids[local_index]; + + // Reject the previously hit primitive to avoid immediate self-intersection. + if (primitive_id == last_hit_primitive) { + return reject_candidate(traversal_ray); + } + + // Scalar queries may provide an arbitrary primitive exclusion history. + // TODO - Think about how to provide arbitrary history checks for batch queries. + for (int i = 0; i < exclude_count; ++i) { + if (exclude_primitives[i] == primitive_id) { + return reject_candidate(traversal_ray); + } + } + + const cuBQL::vec3i vertex_indices = mesh.indices[local_index]; + + cuBQL::vec3d vertices[3] = { + mesh.vertices[vertex_indices.x], + mesh.vertices[vertex_indices.y], + mesh.vertices[vertex_indices.z] + }; + + cuBQL::vec3d normal = cuBQL::cross(vertices[1] - vertices[0], + vertices[2] - vertices[0]); + + double normal_dot_direction = dot(normal, intersection_ray.direction); + + if (surface.reverse_sense) { + normal_dot_direction = -normal_dot_direction; + } + + if (orientation_cull(normal_dot_direction, static_cast(orientation))) { + return reject_candidate(traversal_ray); + } + + auto intersection = plucker_ray_tri_intersect(vertices, + intersection_ray.origin, + intersection_ray.direction, + hit->distance, + intersection_ray.tMin, + false, + 0); + + // store ray payload if hit found + if (intersection.hit) { + hit->distance = intersection.t; + hit->surface = mesh.surface_id; + hit->primitive = primitive_id; + hit->piv = normal_dot_direction > 0.0 ? INSIDE : OUTSIDE; + hit->next_volume = surface.next_volume; + hit->boundary_condition = surface.boundary_condition; + hit->normal = normal; + 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 reject_candidate(traversal_ray); + }; + + // Single level traversal call for a shrinking ray query against the flattened BVH of the volume group. + cuBQL::shrinkingRayQuery::forEachPrim(intersect_prim, volume_group.bvh, traversal_ray); +} +#pragma omp end declare target + +void +intersect_surface_tree_scalar(const cubql::Context& context, + const CuBQLVolumeGroup& volume_group, + 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_group_dd = volume_group.get_device_data(); + const int orientation = static_cast(hit_orientation); + + #pragma omp target device(gpu_id) \ + is_device_ptr(d_exclude_primitives, d_surface_hit) + { + intersect_surface_tree(volume_group_dd, + ray, + d_surface_hit, + orientation, + ID_NONE, + 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 CuBQLVolumeGroup::DD* d_volume_to_group, + XDGRayHit* d_ray_hits, + std::size_t num_rays, + HitOrientation hit_orientation) +{ + if (num_rays == 0) return; + + if (!d_volume_to_group || !d_ray_hits) { + fatal_error("Invalid cuBQL batch intersection buffers"); + } + + const int gpu_id = context.gpuID; + + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(d_volume_to_group, d_ray_hits) + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + const XDGRayHit ray_hit = d_ray_hits[ray_id]; + + CuBQLSurfaceHit hit; + hit.distance = ray_hit.t_max; + hit.surface = ID_NONE; + hit.primitive = ID_NONE; + hit.piv = OUTSIDE; + hit.next_volume = ID_NONE; + + if (ray_hit.volume != ID_NONE) { + CuBQLRay ray; + ray.origin = cuBQL::vec3d(ray_hit.origin[0], ray_hit.origin[1], ray_hit.origin[2]); + ray.direction = cuBQL::vec3d(ray_hit.direction[0], ray_hit.direction[1], ray_hit.direction[2]); + ray.tMin = ray_hit.t_min; + ray.tMax = ray_hit.t_max; + ray.volume = ray_hit.volume; + + const CuBQLVolumeGroup::DD volume_group = d_volume_to_group[ray.volume]; + + intersect_surface_tree(volume_group, + ray, + &hit, + static_cast(hit_orientation), + ray_hit.last_hit_primitive, + nullptr, + 0); + } + + d_ray_hits[ray_id].distance = hit.distance; + d_ray_hits[ray_id].surface = hit.surface; + d_ray_hits[ray_id].primitive = hit.primitive; + d_ray_hits[ray_id].point_in_volume = static_cast(hit.piv); + d_ray_hits[ray_id].next_volume = hit.next_volume; + d_ray_hits[ray_id].boundary_condition = static_cast(hit.boundary_condition); + d_ray_hits[ray_id].normal[0] = hit.normal.x; + d_ray_hits[ray_id].normal[1] = hit.normal.y; + d_ray_hits[ray_id].normal[2] = hit.normal.z; + } +} + +} // namespace xdg diff --git a/src/cuBQL/ray_tracer.cpp b/src/cuBQL/ray_tracer.cpp new file mode 100644 index 00000000..fc25677a --- /dev/null +++ b/src/cuBQL/ray_tracer.cpp @@ -0,0 +1,462 @@ +#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 +#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_group_) { + omp_target_free(d_volume_to_group_, context_.gpuID); + d_volume_to_group_ = nullptr; + } + + for (auto& [tree, group] : tree_to_volume_group_) { + group.release(); + } + + for (auto& [surface, mesh] : surface_to_mesh_) { + mesh.release(); + } +} + +void CuBQLRayTracer::init() +{ + upload_volume_to_group_table_(); + initialized_ = true; +} + +void CuBQLRayTracer::upload_volume_to_group_table_() +{ + if (d_volume_to_group_) { + omp_target_free(d_volume_to_group_, context_.gpuID); + d_volume_to_group_ = nullptr; + } + + if (volume_to_group_.empty()) { + return; + } + + d_volume_to_group_ = static_cast + (omp_target_alloc(volume_to_group_.size() * sizeof(CuBQLVolumeGroup::DD), + context_.gpuID)); + omp_target_memcpy(d_volume_to_group_, + volume_to_group_.data(), + volume_to_group_.size() * sizeof(CuBQLVolumeGroup::DD), + 0, + 0, + context_.gpuID, + 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}; +} + +CuBQLSurfaceMesh +CuBQLRayTracer::register_surface(const std::shared_ptr& mesh_manager, + MeshID surface_id) +{ + const auto num_faces = mesh_manager->num_surface_faces(surface_id); + auto vertices = mesh_manager->get_surface_vertices(surface_id); + auto indices = mesh_manager->get_surface_connectivity(surface_id); + auto h_primitive_ids = mesh_manager->get_surface_faces(surface_id); + + std::vector h_vertices; + h_vertices.reserve(vertices.size()); + 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]); + } + + // 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_ids = static_cast + (omp_target_alloc(h_primitive_ids.size() * sizeof(MeshID), context_.gpuID)); + omp_target_memcpy(d_primitive_ids, + h_primitive_ids.data(), + h_primitive_ids.size() * sizeof(MeshID), + 0, + 0, + context_.gpuID, + context_.hostID); + + 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_ids = d_primitive_ids; + surface_mesh.num_vertices = static_cast(h_vertices.size()); + surface_mesh.num_primitives = static_cast(num_faces); + surface_mesh.gpu_id = context_.gpuID; + + return surface_mesh; +} + +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); + + if (volume_surfaces.empty()) { + fatal_error("Volume {} has no surfaces; cannot build cuBQL surface tree", volume_id); + } + + const double volume_bump = bounding_box_bump(mesh_manager, volume_id); + + std::uint64_t num_volume_primitives = 0; + for (const MeshID surface : volume_surfaces) { + if (!surface_to_mesh_.count(surface)) { + surface_to_mesh_.emplace(surface, register_surface(mesh_manager, surface)); + } + + auto& surface_mesh = surface_to_mesh_.at(surface); + surface_mesh.max_parent_volume_bump = std::max(surface_mesh.max_parent_volume_bump, volume_bump); + num_volume_primitives += surface_mesh.num_primitives; + } + + std::vector h_surfaces; + std::vector h_prim_refs; + h_surfaces.reserve(volume_surfaces.size()); + h_prim_refs.reserve(static_cast(num_volume_primitives)); + + for (const MeshID surface : volume_surfaces) { + const auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surface); + const CuBQLSurfaceMesh& surface_mesh = surface_to_mesh_.at(surface); + + CuBQLVolumeGroup::SurfaceDD surface_data; + surface_data.mesh = surface_mesh.get_device_data(); + + if (volume_id == forward_parent) { + surface_data.next_volume = reverse_parent; + } else if (volume_id == reverse_parent) { + surface_data.reverse_sense = true; + surface_data.next_volume = forward_parent; + } else { + fatal_error("Volume {} is not a parent of surface {}", volume_id, surface); + } + + const auto property = mesh_manager->get_surface_property(surface, PropertyType::BOUNDARY_CONDITION); + if (property.value == "vacuum") { + surface_data.boundary_condition = VACUUM; + } else if (property.value == "reflecting" || property.value == "reflective") { + surface_data.boundary_condition = REFLECTIVE; + } else if (property.value == "transmission") { + surface_data.boundary_condition = TRANSMISSION; + } else { + fatal_error("Unsupported boundary condition '{}' on surface {}", property.value, surface); + } + + const auto surface_index = static_cast(h_surfaces.size()); + h_surfaces.push_back(surface_data); + + for (std::uint32_t prim_index = 0; prim_index < surface_mesh.num_primitives; ++prim_index) { + h_prim_refs.push_back({surface_index, prim_index}); + } + } + + auto* d_aabbs = static_cast + (omp_target_alloc(h_prim_refs.size() * sizeof(cuBQL::box3f), + context_.gpuID)); + auto* d_surfaces = static_cast + (omp_target_alloc(h_surfaces.size() * sizeof(CuBQLVolumeGroup::SurfaceDD), + context_.gpuID)); + omp_target_memcpy(d_surfaces, + h_surfaces.data(), + h_surfaces.size() * sizeof(CuBQLVolumeGroup::SurfaceDD), + 0, + 0, + context_.gpuID, + context_.hostID); + + auto* d_prim_refs = static_cast + (omp_target_alloc(h_prim_refs.size() * sizeof(CuBQLVolumeGroup::PrimRef), + context_.gpuID)); + omp_target_memcpy(d_prim_refs, + h_prim_refs.data(), + h_prim_refs.size() * sizeof(CuBQLVolumeGroup::PrimRef), + 0, + 0, + context_.gpuID, + context_.hostID); + + const auto num_primitives = static_cast(h_prim_refs.size()); + const int gpu_id = context_.gpuID; + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(d_aabbs, d_surfaces, d_prim_refs) + for (std::uint32_t primID = 0; primID < num_primitives; ++primID) { + const auto primitive = d_prim_refs[primID]; + const auto surface = d_surfaces[primitive.surface_index]; + const auto mesh = surface.mesh; + const auto indices = mesh.indices[primitive.primitive_index]; + + cuBQL::box3d aabb; + aabb.extend(mesh.vertices[indices.x]); + aabb.extend(mesh.vertices[indices.y]); + aabb.extend(mesh.vertices[indices.z]); + + const cuBQL::vec3d bump(mesh.max_parent_volume_bump); + aabb.lower = aabb.lower - bump; + aabb.upper = aabb.upper + bump; + d_aabbs[primID] = cuBQL::box3f(aabb); + } + + CuBQLVolumeGroup volume_group; + volume_group.d_surfaces = d_surfaces; + volume_group.d_prim_refs = d_prim_refs; + volume_group.num_surfaces = static_cast(h_surfaces.size()); + volume_group.num_primitives = num_primitives; + volume_group.gpu_id = context_.gpuID; + + cuBQL::BuildConfig build_params; + cuBQL::build_omp_target(volume_group.bvh, + d_aabbs, + num_primitives, + build_params, + context_.gpuID); + + omp_target_free(d_aabbs, context_.gpuID); + + // Retain owning objects for scalar TreeID lookups and allocation lifetime. + auto result = tree_to_volume_group_.emplace(tree, std::move(volume_group)); + auto it = result.first; + + // Keep a dense host-side MeshID -> group device-data table for batch queries. + const auto volume_index = static_cast(volume_id); + if (volume_index >= volume_to_group_.size()) { + volume_to_group_.resize(volume_index + 1); + } + volume_to_group_[volume_index] = it->second.get_device_data(); + + if (initialized_) { + upload_volume_to_group_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 CuBQLVolumeGroup& volume_group = tree_to_volume_group_.at(tree); + + // Use provided direction or if Direction == nulptr use default direction + Direction directionUsed = (direction != nullptr) ? Direction{direction->x, direction->y, direction->z} + : 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_group, + ray, + surface_hit, + HitOrientation::ANY, + exclude_primitives); + + // if the ray hit nothing the point must be outside the volume + if (surface_hit.primitive == ID_NONE) return false; + + 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 CuBQLVolumeGroup& volume_group = tree_to_volume_group_.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_group, + 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}; +} + +XDGRayHitBuffer CuBQLRayTracer::allocate_ray_hits(std::size_t count) const +{ + if (count == 0) { + warning("Request to allocate 0 cuBQL XDG ray-hit buffer; returning empty buffer"); + return {}; + } + + auto* d_ray_hits = static_cast + (omp_target_alloc(count * sizeof(XDGRayHit), context_.gpuID)); + + if (!d_ray_hits) { + fatal_error("Failed to allocate cuBQL XDG ray-hit buffer"); + } + + return {d_ray_hits, count, context_.gpuID}; +} + +void CuBQLRayTracer::free_ray_hits(XDGRayHitBuffer& ray_hits) const +{ + if (!ray_hits.data) { + warning("Request to free empty cuBQL XDG ray-hit buffer; ignoring"); + return; + } + + omp_target_free(ray_hits.data, ray_hits.device_id); + ray_hits = {}; +} + +void +CuBQLRayTracer::ray_fire_batch(const XDGRayHitBuffer& ray_hits, + HitOrientation hit_orientation) const +{ + if (ray_hits.count == 0) return; + + if (!ray_hits.data) { + fatal_error("Invalid cuBQL XDG ray-hit buffer"); + } + + if (!d_volume_to_group_) { + fatal_error("cuBQL volume-group lookup table has not been uploaded"); + } + + intersect_surface_tree_batch(context_, + d_volume_to_group_, + ray_hits.data, + ray_hits.count, + hit_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..9a3cf5ac --- /dev/null +++ b/src/cuBQL/triangles.cpp @@ -0,0 +1,51 @@ +#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_ids) { + omp_target_free(d_primitive_ids, gpu_id); + d_primitive_ids = nullptr; + } + + num_vertices = 0; + num_primitives = 0; +} + +void CuBQLVolumeGroup::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_surfaces) { + omp_target_free(d_surfaces, gpu_id); + d_surfaces = nullptr; + } + if (d_prim_refs) { + omp_target_free(d_prim_refs, gpu_id); + d_prim_refs = nullptr; + } + + bvh.numNodes = 0; + bvh.numPrims = 0; + num_surfaces = 0; + num_primitives = 0; +} + +} // namespace xdg diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index b2e6776d..b979d6d8 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -1,6 +1,6 @@ #include "xdg/gprt/ray_tracer.h" #include "gprt/gprt.h" -#include "xdg/gprt/vulkan_probe.h" +#include "xdg/available_device_probe.h" namespace xdg { diff --git a/src/ray_tracing_interface.cpp b/src/ray_tracing_interface.cpp index 66e0adb8..fe9fb6eb 100644 --- a/src/ray_tracing_interface.cpp +++ b/src/ray_tracing_interface.cpp @@ -1,5 +1,6 @@ #include #include "xdg/ray_tracing_interface.h" +#include "xdg/error.h" // Any methods which are identical for all RT backends should be defined here @@ -7,6 +8,23 @@ namespace xdg { RayTracer::~RayTracer() {} +XDGRayHitBuffer RayTracer::allocate_ray_hits(std::size_t) const +{ + fatal_error("Selected ray tracer does not support device batch ray fire"); + return {}; +} + +void RayTracer::free_ray_hits(XDGRayHitBuffer&) const +{ + fatal_error("Selected ray tracer does not support device batch ray fire"); +} + +void RayTracer::ray_fire_batch(const XDGRayHitBuffer&, + HitOrientation) const +{ + fatal_error("Selected ray tracer does not support device batch ray fire"); +} + SurfaceTreeID RayTracer::next_surface_tree_id() { return ++next_surface_tree_id_; @@ -23,4 +41,4 @@ const double RayTracer::bounding_box_bump(const std::shared_ptr mes return std::max(volume_bounding_box.dilation(), numerical_precision_); } -} // namespace xdg \ No newline at end of file +} // namespace xdg diff --git a/src/xdg.cpp b/src/xdg.cpp index 75face29..30a84710 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); }; @@ -246,6 +260,22 @@ XDG::ray_fire(MeshID volume, return ray_tracing_interface()->ray_fire(scene, origin, direction, dist_limit, orientation, exclude_primitives); } +XDGRayHitBuffer XDG::allocate_ray_hits(std::size_t count) const +{ + return ray_tracing_interface()->allocate_ray_hits(count); +} + +void XDG::free_ray_hits(XDGRayHitBuffer& ray_hits) const +{ + ray_tracing_interface()->free_ray_hits(ray_hits); +} + +void XDG::ray_fire_batch(const XDGRayHitBuffer& ray_hits, + HitOrientation hit_orientation) const +{ + ray_tracing_interface()->ray_fire_batch(ray_hits, hit_orientation); +} + std::pair XDG::closest(MeshID volume, const Position& origin) const { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 23b95f1d..fc7e0837 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -39,6 +39,10 @@ if (XDG_ENABLE_MOAB) list(APPEND TEST_NAMES test_ray_tracer_cross_check) endif() +if (XDG_ENABLE_CUBQL) + list(APPEND TEST_NAMES test_xdg_batch_api) +endif() + if (XDG_ENABLE_MOAB AND XDG_BUILD_TOOLS) list(APPEND TEST_NAMES test_overlap_check) endif() @@ -55,6 +59,9 @@ foreach(test ${TEST_NAMES}) if (XDG_ENABLE_MOAB) target_link_libraries(${test} PRIVATE MOAB) endif() + if (XDG_ENABLE_CUBQL) + target_link_libraries(${test} PRIVATE $) + endif() set_target_properties(${test} PROPERTIES BUILD_RPATH "$") catch_discover_tests(${test} diff --git a/tests/test_point_in_volume.cpp b/tests/test_point_in_volume.cpp index bb3b1451..b11a25d8 100644 --- a/tests/test_point_in_volume.cpp +++ b/tests/test_point_in_volume.cpp @@ -16,7 +16,8 @@ using namespace xdg::test; TEMPLATE_TEST_CASE("Point-in-volume on MockedTriTetMesh", "[piv][mock]", Embree_Raytracer, - GPRT_Raytracer) + GPRT_Raytracer, + CuBQL_Raytracer) { constexpr auto rt_backend = TestType::value; diff --git a/tests/test_ray_fire.cpp b/tests/test_ray_fire.cpp index 501fa5f6..de6c0106 100644 --- a/tests/test_ray_fire.cpp +++ b/tests/test_ray_fire.cpp @@ -18,7 +18,8 @@ using namespace xdg::test; TEMPLATE_TEST_CASE("Ray Fire on MockedTriTetMesh (per-backend sections)", "[rayfire][mock]", Embree_Raytracer, - GPRT_Raytracer) + GPRT_Raytracer, + CuBQL_Raytracer) { // Generate one test run per enabled backend constexpr auto rt_backend = TestType::value; diff --git a/tests/test_ray_tracer_cross_check.cpp b/tests/test_ray_tracer_cross_check.cpp index c75b0273..d2407edf 100644 --- a/tests/test_ray_tracer_cross_check.cpp +++ b/tests/test_ray_tracer_cross_check.cpp @@ -57,7 +57,7 @@ TEST_CASE("Test Pincell RT libraries Cross-Check ray_fire queries", "[moab][rayf } srand48(12345); // set fixed seed for rng - std::vector directions(1000); + std::vector directions(100); for (auto &dir : directions) { dir = rand_dir(); } @@ -86,4 +86,59 @@ TEST_CASE("Test Pincell RT libraries Cross-Check ray_fire queries", "[moab][rayf } } +TEST_CASE("Test Pincell RT libraries Cross-Check point_in_volume queries", "[moab][piv][cross-check]") +{ + const auto rt_cases = make_rt_cases("pincell.h5m"); + if (rt_cases.size() < 2) { + SKIP("Fewer than two ray tracing backends are available; skipping cross-check."); + } + + const std::array directions { + Direction {1.0, 0.0, 0.0}, // axis-aligned x ray + Direction {0.0, 1.0, 0.0}, // axis-aligned y ray + Direction {0.0, 0.0, 1.0}, // axis-aligned z ray + Direction {0.371390676, 0.557086014, 0.742781353} // non-axis ray + }; + + // Points generated by codex with a h5dump to test various edge cases around the pincell geometry + const std::array points { + Position {0.0, 0.0, 0.0}, // pincell center + Position {8.5, 0.25, 0.0}, // inside inner cylinder, r < 9 + Position {9.5, 0.25, 0.0}, // between cylinder radii, 9 < r < 10 + Position {10.5, 0.25, 0.0}, // outside outer cylinder, r > 10 + Position {20.0, 0.25, 0.0}, // inside square cell away from cylinder + Position {24.5, 0.25, 0.0}, // just inside x = 25 square boundary + Position {25.5, 0.25, 0.0}, // between x = 25 and outer x = 27.5 + Position {27.0, 0.25, 0.0}, // just inside outer x = 27.5 boundary + Position {28.0, 0.25, 0.0}, // outside outer x boundary + Position {0.25, 0.0, 19.5}, // just inside positive z = 20 cap + Position {0.25, 0.0, 20.5}, // just outside positive z = 20 cap + Position {0.25, 0.0, -19.5}, // just inside negative z = -20 cap + Position {0.25, 0.0, -20.5}, // just outside negative z = -20 cap + Position {0.25, 24.5, 0.0}, // just inside y = 25 square boundary + Position {0.25, 28.0, 0.0} // outside outer y boundary + }; + + const auto& reference_case = rt_cases.front(); + + for (const auto& volume : reference_case.xdg->mesh_manager()->volumes()) { + for (const auto& point : points) { + for (const auto& direction : directions) { + const auto reference_result = + reference_case.xdg->point_in_volume(volume, point, &direction); + + for (size_t i = 1; i < rt_cases.size(); ++i) { + const auto& candidate = rt_cases[i]; + const auto candidate_result = + candidate.xdg->point_in_volume(volume, point, &direction); + + CAPTURE(volume, point, direction, reference_case.name, candidate.name, + reference_result, candidate_result); + REQUIRE(candidate_result == reference_result); + } + } + } + } +} + // TODO - Add all of the other queries diff --git a/tests/test_xdg_batch_api.cpp b/tests/test_xdg_batch_api.cpp new file mode 100644 index 00000000..b2322693 --- /dev/null +++ b/tests/test_xdg_batch_api.cpp @@ -0,0 +1,264 @@ +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "mesh_mocks.h" +#include "util.h" +#include "xdg/constants.h" +#include "xdg/device_ray.h" +#include "xdg/util/rng.h" +#include "xdg/xdg.h" + +using namespace xdg; +using namespace xdg::test; + +constexpr double tolerance = 1e-6; + +TEST_CASE("XDG batch ray-hit buffer lifecycle", "[batch][cubql]") +{ + check_ray_tracer_supported(RTLibrary::CUBQL); + + auto mesh_manager = std::make_shared(false); + auto xdg = std::make_shared(mesh_manager, RTLibrary::CUBQL); + + const XDGRayHitBuffer empty_buffer = xdg->allocate_ray_hits(0); + REQUIRE(empty_buffer.data == nullptr); + REQUIRE(empty_buffer.count == 0); + REQUIRE(empty_buffer.device_id == -1); + + XDGRayHitBuffer ray_hits = xdg->allocate_ray_hits(64); + REQUIRE(ray_hits.data != nullptr); + REQUIRE(ray_hits.count == 64); + REQUIRE(ray_hits.device_id >= 0); + xdg->free_ray_hits(ray_hits); + REQUIRE(ray_hits.data == nullptr); + REQUIRE(ray_hits.count == 0); + REQUIRE(ray_hits.device_id == -1); +} + +TEST_CASE("XDG batch ray fire accepts an empty buffer", + "[rayfire][batch][cubql]") +{ + check_ray_tracer_supported(RTLibrary::CUBQL); + + auto mesh_manager = std::make_shared(false); + auto xdg = std::make_shared(mesh_manager, RTLibrary::CUBQL); + + REQUIRE_NOTHROW(xdg->ray_fire_batch({})); +} + +// The batch ray-fire API is currently implemented only by cuBQL. +// This test also populates its input buffer using OpenMP target offload, +// so support for another backend may require a different population path. +TEST_CASE("XDG batch ray fire matches scalar queries on MeshMock", + "[rayfire][batch][cubql]") +{ + check_ray_tracer_supported(RTLibrary::CUBQL); + + auto mesh_manager = std::make_shared(false); + mesh_manager->init(); + + auto xdg = std::make_shared(mesh_manager, RTLibrary::CUBQL); + xdg->prepare_raytracer(); + + const std::size_t batch_sizes[] = {1, 1024}; + for (std::size_t num_rays : batch_sizes) { + DYNAMIC_SECTION("N=" << num_rays) + { + XDGRayHitBuffer ray_hits = xdg->allocate_ray_hits(num_rays); + + XDGRayHit* d_ray_hits = ray_hits.data; + const int device_id = ray_hits.device_id; + const std::uint32_t seed = 12345u; + const MeshID volume = mesh_manager->volumes()[0]; + + #pragma omp target teams distribute parallel for device(device_id) \ + is_device_ptr(d_ray_hits) + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + std::uint32_t state = seed ^ static_cast(ray_id); + + XDGRayHit ray_hit {}; + ray_hit.origin[0] = 0.0; + ray_hit.origin[1] = 0.0; + ray_hit.origin[2] = 0.0; + random_unit_direction_lcg(state, ray_hit.direction); + ray_hit.t_min = 0.0; + ray_hit.t_max = INFTY; + ray_hit.volume = volume; + ray_hit.last_hit_primitive = ID_NONE; + d_ray_hits[ray_id] = ray_hit; + } + + xdg->ray_fire_batch(ray_hits); // Perform device side ray tracing + + std::vector batch_results(num_rays); + const int copy_status = omp_target_memcpy(batch_results.data(), + ray_hits.data, + num_rays * sizeof(XDGRayHit), + 0, + 0, + omp_get_initial_device(), + ray_hits.device_id); + REQUIRE(copy_status == 0); // Ensure omp memcpy was successful before checking results + + xdg->free_ray_hits(ray_hits); + + // Compare cuBQL batch results with cuBQL scalar ray_fire results + // Note that cuBQL scalar ray_fire is cross-checked against embree in test_ray_tracer_cross_check.cpp + const Position origin {0.0, 0.0, 0.0}; + for (std::size_t ray_id = 0; ray_id < num_rays; ++ray_id) { + std::uint32_t state = seed ^ static_cast(ray_id); + double direction_values[3]; + random_unit_direction_lcg(state, direction_values); + const Direction direction {direction_values[0], direction_values[1], direction_values[2]}; + const auto scalar_result = xdg->ray_fire(volume, origin, direction); + + INFO("ray index: " << ray_id); + INFO("direction: " << direction); + INFO("batch result: surface=" << batch_results[ray_id].surface << ", distance=" << std::setprecision(17) << batch_results[ray_id].distance); + INFO("scalar result: surface=" << scalar_result.second << ", distance=" << std::setprecision(17) << scalar_result.first); + REQUIRE(batch_results[ray_id].surface == scalar_result.second); + REQUIRE_THAT(batch_results[ray_id].distance, Catch::Matchers::WithinAbs(scalar_result.first, tolerance)); + } + } + } +} + +TEST_CASE("XDG batch ray fire multi-volume check full rayhit results", + "[rayfire][batch][cubql]") +{ + check_ray_tracer_supported(RTLibrary::CUBQL); + + auto xdg = XDG::create(MeshLibrary::MOAB, RTLibrary::CUBQL); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("pwr_pincell.h5m"); + mesh_manager->init(); + mesh_manager->parse_metadata(); + xdg->prepare_raytracer(); + + struct VolumeCase { + MeshID volume; + Position origin; + }; + + const std::array volume_cases {{ + {1, {0.0, 0.01, 0.1}}, + {2, {0.42, 0.01, 0.1}}, + {3, {0.50, 0.01, 0.1}} + }}; + + for (const auto& volume_case : volume_cases) { + CAPTURE(volume_case.volume, volume_case.origin); + REQUIRE(xdg->point_in_volume(volume_case.volume, volume_case.origin)); + } + + const std::array axis_aligned_directions {{ + {1.0, 0.0, 0.0}, + {-1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, -1.0, 0.0}, + {0.0, 0.0, 1.0}, + {0.0, 0.0, -1.0} + }}; + + // Build host-side ray-hit buffer and copy to device to handle every ray in one batch + std::vector host_ray_hits; + host_ray_hits.reserve(volume_cases.size() * axis_aligned_directions.size()); + + for (const auto& direction : axis_aligned_directions) { + for (const auto& volume_case : volume_cases) { + XDGRayHit ray_hit {}; + ray_hit.origin[0] = volume_case.origin.x; + ray_hit.origin[1] = volume_case.origin.y; + ray_hit.origin[2] = volume_case.origin.z; + ray_hit.direction[0] = direction.x; + ray_hit.direction[1] = direction.y; + ray_hit.direction[2] = direction.z; + ray_hit.t_min = 0.0; + ray_hit.t_max = INFTY; + ray_hit.volume = volume_case.volume; + ray_hit.last_hit_primitive = ID_NONE; + host_ray_hits.push_back(ray_hit); + } + } + + XDGRayHitBuffer ray_hits = xdg->allocate_ray_hits(host_ray_hits.size()); + const std::size_t buffer_size = host_ray_hits.size() * sizeof(XDGRayHit); + REQUIRE(omp_target_memcpy(ray_hits.data, + host_ray_hits.data(), + buffer_size, + 0, + 0, + ray_hits.device_id, + omp_get_initial_device()) == 0); + + xdg->ray_fire_batch(ray_hits); + + REQUIRE(omp_target_memcpy(host_ray_hits.data(), + ray_hits.data, + buffer_size, + 0, + 0, + omp_get_initial_device(), + ray_hits.device_id) == 0); + + xdg->free_ray_hits(ray_hits); + + // Loop over each ray and perform a scalar ray_fire to compare against the batch results + for (std::size_t ray_id = 0; ray_id < host_ray_hits.size(); ++ray_id) { + const auto& batch_ray_hit = host_ray_hits[ray_id]; // recover rayhit for this ray from batch results + const Position origin {batch_ray_hit.origin[0], batch_ray_hit.origin[1], batch_ray_hit.origin[2]}; + const Direction direction {batch_ray_hit.direction[0], batch_ray_hit.direction[1], batch_ray_hit.direction[2]}; + std::vector last_hit_prim; + const auto scalar_hit = xdg->ray_fire(batch_ray_hit.volume, origin, direction, INFTY, HitOrientation::EXITING, &last_hit_prim); + + REQUIRE(scalar_hit.second != ID_NONE); + REQUIRE(last_hit_prim.size() == 1); + + const MeshID expected_surface = scalar_hit.second; + const double expected_distance = scalar_hit.first; + const MeshID expected_primitive = last_hit_prim.back(); + const PointInVolume expected_point_in_volume = xdg->point_in_volume(batch_ray_hit.volume, origin, &direction) ? INSIDE : OUTSIDE; + const MeshID expected_next_volume = mesh_manager->next_volume(batch_ray_hit.volume, expected_surface); + const auto boundary_property = mesh_manager->get_surface_property(expected_surface, PropertyType::BOUNDARY_CONDITION); + SurfaceBoundaryCondition expected_boundary_condition = UNSET; + if (boundary_property.value == "transmission") { + expected_boundary_condition = TRANSMISSION; + } else if (boundary_property.value == "reflecting") { + expected_boundary_condition = REFLECTIVE; + } + const Direction expected_normal = mesh_manager->face_normal(expected_primitive); + + // Outputting full rayhit state on failure + INFO("ray index: " << ray_id); + INFO("input volume: " << batch_ray_hit.volume); + INFO("origin: " << origin); + INFO("direction: " << direction); + INFO("surface: batch=" << batch_ray_hit.surface << ", scalar=" << expected_surface); + INFO("distance: batch=" << std::setprecision(17) << batch_ray_hit.distance << ", scalar=" << expected_distance); + INFO("primitive: batch=" << batch_ray_hit.primitive << ", scalar=" << expected_primitive); + INFO("point in volume: batch=" << batch_ray_hit.point_in_volume << ", scalar=" << expected_point_in_volume); + INFO("next volume: batch=" << batch_ray_hit.next_volume << ", scalar=" << expected_next_volume); + INFO("boundary condition: batch=" << batch_ray_hit.boundary_condition << ", scalar=" << expected_boundary_condition << " (" << boundary_property.value << ")"); + INFO("normal: batch=[" << batch_ray_hit.normal[0] << " " << batch_ray_hit.normal[1] << " " << batch_ray_hit.normal[2] << "], scalar=" << expected_normal); + + // Perform assertions per ray + REQUIRE(expected_boundary_condition != UNSET); + REQUIRE(batch_ray_hit.surface == expected_surface); + REQUIRE_THAT(batch_ray_hit.distance, Catch::Matchers::WithinAbs(expected_distance, tolerance)); + REQUIRE(batch_ray_hit.primitive == expected_primitive); + REQUIRE(batch_ray_hit.point_in_volume == expected_point_in_volume); + REQUIRE(batch_ray_hit.next_volume == expected_next_volume); + REQUIRE(batch_ray_hit.boundary_condition == expected_boundary_condition); + REQUIRE_THAT(batch_ray_hit.normal[0], Catch::Matchers::WithinAbs(expected_normal.x, tolerance)); + REQUIRE_THAT(batch_ray_hit.normal[1], Catch::Matchers::WithinAbs(expected_normal.y, tolerance)); + REQUIRE_THAT(batch_ray_hit.normal[2], Catch::Matchers::WithinAbs(expected_normal.z, tolerance)); + } +} 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..dab27410 100644 --- a/tools/ray_benchmark.cpp +++ b/tools/ray_benchmark.cpp @@ -19,7 +19,6 @@ #include "ray_benchmark.h" - using namespace xdg; int main(int argc, char** argv) @@ -60,7 +59,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 +97,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 +172,113 @@ 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) { + // Generate random rays directly on the target device. + generation_timer.start(); + + XDGRayHitBuffer ray_hits = xdg->allocate_ray_hits(num_rays); + XDGRayHit* d_ray_hits = ray_hits.data; + const std::size_t ray_count = ray_hits.count; + const int gpu_id = ray_hits.device_id; + + const double origin_x = origin.x; + const double origin_y = origin.y; + const double origin_z = origin.z; + + #pragma omp target teams distribute parallel for device(gpu_id) is_device_ptr(d_ray_hits) + for (std::size_t ray_id = 0; ray_id < ray_count; ++ray_id) { + std::uint32_t state = seed ^ static_cast(ray_id); + + auto sample = tools::benchmark::random_spherical_source(origin_x, + origin_y, + origin_z, + state, + source_radius); + + XDGRayHit ray_hit; + ray_hit.origin[0] = sample.position[0]; + ray_hit.origin[1] = sample.position[1]; + ray_hit.origin[2] = sample.position[2]; + ray_hit.direction[0] = sample.direction[0]; + ray_hit.direction[1] = sample.direction[1]; + ray_hit.direction[2] = sample.direction[2]; + ray_hit.t_min = 0.0; + ray_hit.t_max = INFTY; + ray_hit.volume = volume; + ray_hit.last_hit_primitive = ID_NONE; + ray_hit.distance = INFTY; + ray_hit.surface = ID_NONE; + ray_hit.primitive = ID_NONE; + ray_hit.point_in_volume = OUTSIDE; + + d_ray_hits[ray_id] = ray_hit; + } + + generation_timer.stop(); + + // Trace rays and count hits on the target device. + trace_timer.start(); + xdg->ray_fire_batch(ray_hits); + trace_timer.stop(); + + #pragma omp target teams distribute parallel for device(gpu_id) \ + is_device_ptr(d_ray_hits) reduction(+:num_hits) + for (std::size_t ray_id = 0; ray_id < ray_count; ++ray_id) { + if (d_ray_hits[ray_id].surface != ID_NONE) num_hits++; + } + + xdg->free_ray_hits(ray_hits); } - - trace_timer.stop(); const std::size_t num_misses = num_rays - num_hits; const double hit_fraction = num_rays > 0 ? static_cast(num_hits) / static_cast(num_rays) : 0.0; + const double generation_time = generation_timer.elapsed(); const double trace_time = trace_timer.elapsed(); const double end_to_end_time = generation_time + trace_time; diff --git a/tools/ray_benchmark.h b/tools/ray_benchmark.h index 65a48339..c255935d 100644 --- a/tools/ray_benchmark.h +++ b/tools/ray_benchmark.h @@ -4,6 +4,8 @@ #include #include +#include "xdg/util/rng.h" + namespace xdg::tools::benchmark { struct SourceSample { @@ -15,30 +17,6 @@ struct SourceSample { #pragma omp declare target #endif -inline double rand01(std::uint32_t& state) -{ - state = state * 1664525u + 1013904223u; - return static_cast(state) * (1.0 / 4294967296.0); -} - -inline void random_unit_dir_lcg(std::uint32_t& state, double direction[3]) -{ - double x1; - double x2; - double s; - - do { - x1 = rand01(state) * 2.0 - 1.0; - x2 = rand01(state) * 2.0 - 1.0; - s = x1 * x1 + x2 * x2; - } while (s <= 0.0 || s >= 1.0); - - const double t = 2.0 * std::sqrt(1.0 - s); - direction[0] = x1 * t; - direction[1] = x2 * t; - direction[2] = 1.0 - 2.0 * s; -} - inline SourceSample random_spherical_source(double origin_x, double origin_y, double origin_z, @@ -46,14 +24,14 @@ inline SourceSample random_spherical_source(double origin_x, double source_radius) { SourceSample sample; - random_unit_dir_lcg(state, sample.direction); + random_unit_direction_lcg(state, sample.direction); sample.position[0] = origin_x; sample.position[1] = origin_y; sample.position[2] = origin_z; if (source_radius > 0.0) { - const double radius = source_radius * std::cbrt(rand01(state)); + const double radius = source_radius * std::cbrt(lcg_rand01(state)); sample.position[0] += sample.direction[0] * radius; sample.position[1] += sample.direction[1] * radius; sample.position[2] += sample.direction[2] * radius; 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