diff --git a/include/xdg/gprt/ray.h b/include/xdg/gprt/ray.h new file mode 100644 index 00000000..d82eb8d5 --- /dev/null +++ b/include/xdg/gprt/ray.h @@ -0,0 +1,40 @@ +#ifndef _XDG_GPRT_RAY_H +#define _XDG_GPRT_RAY_H + +#include "gprt.h" +#include "../shared_enums.h" + +/* + * Double-precision ray and hit structures used by the GPRT backend. + * + * These types are not inherently GPRT-specific, but we keep them here for now + * since GPRT is the only GPU backend. If another GPU backend is added, these + * can be reused. Unifying them with the CPU/Embree types is possible, but may + * not be worth the added complexity at this stage. + */ + +namespace xdg { + +struct dblRay +{ + double3 origin; + double3 direction; + int volume_mesh_id; // MeshID of the volume this ray will be traced against + uint enabled; // Flag to indicate if the ray is active + int32_t* exclude_primitives; // Optional for excluding primitives + int32_t exclude_count; // Number of excluded primitives +}; + + +struct dblHit +{ + double distance; + int surf_id; + int primitive_id; + PointInVolume piv; // Point in volume check result (0 for outside, 1 for inside) +}; + +} + + +#endif diff --git a/include/xdg/gprt/ray_tracer.h b/include/xdg/gprt/ray_tracer.h index 8d24d107..03580e42 100644 --- a/include/xdg/gprt/ray_tracer.h +++ b/include/xdg/gprt/ray_tracer.h @@ -7,12 +7,9 @@ #include "xdg/constants.h" #include "xdg/mesh_manager_interface.h" -#include "xdg/primitive_ref.h" -#include "xdg/geometry_data.h" #include "xdg/ray_tracing_interface.h" -#include "xdg/ray.h" #include "xdg/error.h" -#include "gprt/gprt.h" + #include "shared_structs.h" extern GPRTProgram dbl_deviceCode; @@ -26,17 +23,16 @@ enum class RayGenType { }; struct gprtRayHit { - size_t capacity = 1; // Max number of rays allocated - size_t size = 0; // Current number of active rays + DeviceRayHitBuffers view; // external facing POD for rayhit buffers + size_t size = 0; // Current number of active rays GPRTBufferOf ray = nullptr; GPRTBufferOf hit = nullptr; - dblRay* devRayAddr = nullptr; - dblHit* devHitAddr = nullptr; - bool is_valid() const { return capacity > 0 && ray && hit && devRayAddr && devHitAddr; } + bool is_valid() const { + return view.capacity > 0 && ray && hit && view.rayDevPtr && view.hitDevPtr; + } }; - class GPRTRayTracer : public RayTracer { public: GPRTRayTracer(); @@ -83,15 +79,40 @@ class GPRTRayTracer : public RayTracer { const Direction* direction = nullptr, const std::vector* exclude_primitives = nullptr) const override; + void point_in_volume(TreeID tree, + const Position* points, + const size_t num_points, + uint8_t* results, + const Direction* directions = nullptr, + std::vector* exclude_primitives = nullptr) override; + std::pair ray_fire(TreeID scene, 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(TreeID tree, + const Position* origins, + const Direction* directions, + const size_t num_rays, + double* hitDistances, + MeshID* surfaceIDs, + const double dist_limit = INFTY, + HitOrientation orientation = HitOrientation::EXITING, + std::vector* const exclude_primitives = nullptr) override; + + void ray_fire_prepared(const size_t num_rays, + const double dist_limit = INFTY, + HitOrientation orientation = HitOrientation::EXITING) override; + + void point_in_volume_prepared(const size_t num_rays) override; std::pair closest(TreeID scene, - const Position& origin) override {}; + const Position& origin) override { + fatal_error("Closest queries are not currently supported with GPRT ray tracer"); + return {INFTY, ID_NONE}; + }; bool occluded(TreeID scene, const Position& origin, @@ -100,9 +121,41 @@ class GPRTRayTracer : public RayTracer { fatal_error("Occlusion queries are not currently supported with GPRT ray tracer"); return false; } - + + // Check to see if buffers large enough and resize if not + void check_rayhit_buffer_capacity(const size_t N) override; + + // Method to expose device ray and hit buffers for external population + DeviceRayHitBuffers get_device_rayhit_buffers(const size_t N) override; + + /** + * @brief Allocate device buffers and invoke a callback to populate them + * + * This method enables downstream applications to populate ray buffers using + * any compute API (GPRT, CUDA, HIP, etc.) without XDG needing to know the details. + */ + void populate_rays_external(size_t numRays, + const RayPopulationCallback& callback) override; + + void transfer_hits_buffer_to_host(const size_t num_rays, + std::vector& hits) override; + + GPRTContext context() + { + return context_; + } + + SurfaceAccelerationStructure* tlas_handle_device_ptr() const + { + return gprtBufferGetDevicePointer(tlas_handle_buffer_); + } + + size_t tlas_handle_count() const + { + return tlas_handles_.size(); + } + private: - void check_ray_buffer_capacity(size_t N); // GPRT objects GPRTContext context_; @@ -133,7 +186,31 @@ class GPRTRayTracer : public RayTracer { // Internal GPRT Mappings std::unordered_map surface_volume_tree_to_accel_map; // Map from XDG::TreeID to GPRTAccel for volume TLAS - std::vector blas_handles_; // Store BLAS handles so that they can be explicitly referenced in destructor + std::unordered_map surface_tree_to_volume_map_; + std::vector tlas_handles_; // Host side storage of TLAS device addresses + GPRTBufferOf tlas_handle_buffer_; // Device buffer for TLAS addresses + std::vector meshid_to_sense_; // Host-side MeshID -> sense map + GPRTBufferOf meshid_to_sense_buffer_ {nullptr}; // Device buffer for MeshID -> sense map + bool initialized_ {false}; // flag to indicate if init() has been called + + void update_tlas_table_(); + void update_meshid_to_sense_(); + + template + void upload_device_buffer_(GPRTBufferOf& buf, const std::vector& host_data) + { + if (host_data.empty()) return; + + if (!buf) { + buf = gprtDeviceBufferCreate(context_, host_data.size(), host_data.data()); + return; + } + + gprtBufferResize(context_, buf, host_data.size(), false); + gprtBufferMap(buf); + std::copy(host_data.begin(), host_data.end(), gprtBufferGetHostPointer(buf)); + gprtBufferUnmap(buf); + } // Global Tree IDs GPRTAccel global_surface_accel_ {nullptr}; @@ -142,5 +219,4 @@ class GPRTRayTracer : public RayTracer { }; } // namespace xdg - -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/gprt/shared_structs.h b/include/xdg/gprt/shared_structs.h index 3855775a..915898f5 100644 --- a/include/xdg/gprt/shared_structs.h +++ b/include/xdg/gprt/shared_structs.h @@ -1,32 +1,15 @@ +#ifndef XDG_GPRT_SHARED_STRUCTS_H +#define XDG_GPRT_SHARED_STRUCTS_H + #include "gprt.h" #include "../shared_enums.h" +#include "ray.h" struct GPRTPrimitiveRef { int id; // ID of the primitive - int sense; -}; - -struct dblRay -{ - double3 origin; - double3 direction; - double tMin; // Minimum distance for ray intersection - double tMax; // Maximum distance for ray intersection - int32_t* exclude_primitives; // Optional for excluding primitives - int32_t exclude_count; // Number of excluded primitives - xdg::HitOrientation hitOrientation; - int volume_tree; // TreeID of the volume being queried - SurfaceAccelerationStructure volume_accel; // The volume accel }; -struct dblHit -{ - double distance; - int surf_id; - int primitive_id; - xdg::PointInVolume piv; // Point in volume check result (0 for outside, 1 for inside) -}; /* variables for double precision triangle mesh geometry */ struct DPTriangleGeomData { @@ -35,20 +18,17 @@ struct DPTriangleGeomData { uint3 *index; // index buffer double3 *normals; // normals buffer int surf_id; - int2 vols; - int forward_vol; - int reverse_vol; - dblRay *ray; // double precision rays + int* meshid_to_sense; // MeshID -> sense (+1 forward, -1 reverse) + xdg::dblRay *ray; // double precision rays xdg::HitOrientation hitOrientation; - int forward_tree; // TreeID of the forward volume - int reverse_tree; // TreeID of the reverse volume GPRTPrimitiveRef* primitive_refs; int num_faces; // Number of faces in the geometry }; struct dblRayGenData { - dblRay *ray; - dblHit *hit; + xdg::dblRay *ray; + xdg::dblHit *hit; + SurfaceAccelerationStructure* meshid_to_accel_address; // MeshID->TLAS address table to recover volume to trace against }; /* A small structure of constants that can change every frame without rebuilding the @@ -57,4 +37,7 @@ struct dblRayGenData { struct dblRayFirePushConstants { double tMax; double tMin; + xdg::HitOrientation hitOrientation; }; + +#endif diff --git a/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h index 9d938978..37266e4e 100644 --- a/include/xdg/ray_tracing_interface.h +++ b/include/xdg/ray_tracing_interface.h @@ -4,16 +4,61 @@ #include #include #include +#include +#include "xdg/error.h" #include "xdg/constants.h" #include "xdg/embree_interface.h" #include "xdg/mesh_manager_interface.h" #include "xdg/primitive_ref.h" #include "xdg/geometry_data.h" + namespace xdg { +struct dblHit; // forward declaration for dblHit + +/** + * @brief Device ray/hit buffer descriptor + * + * This structure provides access to device-allocated ray and hit buffers + * in a backend-agnostic way. The buffers contain XDG's standard ray and hit + * data structures (dblRay and dblHit), regardless of which compute backend + * is being used. + * + * Key design principle: + * - Device pointers are opaque (void*) + * - The data layout is always the XDG types dblRay and dblHit + * - Downstream code can write to these buffers (hopefully) using any compute API + * + * For type-safe access in downstream code: + * - Cast rayDevPtr to (dblRay*) when using C++ or kernels + * - Cast hitDevPtr to (dblHit*) when reading hit results + */ +struct DeviceRayHitBuffers { + void* rayDevPtr; + void* hitDevPtr; + size_t capacity; // Number of rays the buffer can hold + size_t rayStride; // Bytes between ray elements - sizeof(dblRay) + size_t hitStride; // Bytes between hit elements - sizeof(dblHit) +}; + +/** + * @brief Callback alias for external ray population + * + * Allows downstream applications to populate ray buffers using their own compute backend + * (GPRT, CUDA, OpenMP) without XDG needing to know the specifics. + * + * The callback receives opaque device pointers and should interpret them according to + * the buffer metadata (stride information). Alternatively, users can rely on the standard + * dblRay/dblHit layouts if they don't need custom padding/alignment. + * + * @param buffer Device ray buffer descriptor with opaque pointers and metadata + * @param numRays Number of rays to generate/populate + */ +using RayPopulationCallback = std::function; + class RayTracer { public: // Constructors/Destructors @@ -73,12 +118,41 @@ class RayTracer { */ virtual void create_global_element_tree() = 0; - // Query Methods + /** + * @brief Check whether a point lies in a specified volume + * + * This method performs a check to see whether a given point is inside a volume provided. + * It computes this by firing a ray from the point and checking whether or not the ray is Entering or Exiting + * the volume boundary. If no direction is provided, a default direction will be used. + * Note - zero length direction vectors are not explicitly checked for internally and should be avoided to avoid causing undefined behavior. + * + * @param[in] tree The TreeID of the volume we are querying against + * @param[in] point The point to be queried + * @param[in] direction (optional) direction to launch a ray in a specified direction - must be non-zero length + * @param[in] exclude_primitives (optional) vector of surface element MeshIDs to exclude from intersection tests + * @return Boolean result of point in volume check + */ virtual bool point_in_volume(TreeID tree, const Position& point, const Direction* direction = nullptr, const std::vector* exclude_primitives = nullptr) const = 0; - + + /** + * @brief Fire a ray against a given volume and return the first hit + * + * This method fires a ray from a given origin in a specified direction against the surfaces of a volume. + * It returns the distance to the closest hit and the MeshID of the surface hit. The user can specify + * a distance limit and whether Entering/Exiting hits should be rejected. + * Note - zero length direction vectors are not explicitly checked for internally and should be avoided to avoid causing undefined behavior. + * + * @param[in] tree The TreeID of the volume we are querying against + * @param[in] origin An array of Position objects representing the starting points of the rays + * @param[in] direction (optional) Direction object to launch a ray in a specified direction + * @param[in] dist_limit (optional) maximum distance to consider for intersections + * @param[in] orientation (optional) flag to consider whether Entering/Exiting hits should be rejected. Defaults to EXITING + * @param[in] exclude_primitives (optional) vector of surface element MeshIDs to exclude from intersection tests + * @return A pair containing the distance to the closest hit and the MeshID of the surface hit + */ virtual std::pair ray_fire(TreeID tree, const Position& origin, const Direction& direction, @@ -123,6 +197,89 @@ class RayTracer { int num_registered_surface_trees() const { return surface_trees_.size(); }; int num_registered_element_trees() const { return element_trees_.size(); }; + + // GPU Ray Tracing Support + + virtual void point_in_volume(TreeID tree, + const Position* points, + const size_t num_points, + uint8_t* results, + const Direction* directions = nullptr, + std::vector* exclude_primitives = nullptr) + { + fatal_error("GPU ray tracing not supported with this RayTracer backend"); + } + + virtual void ray_fire(TreeID tree, + const Position* origins, + const Direction* directions, + const size_t num_rays, + double* hitDistances, + MeshID* surfaceIDs, + const double dist_limit = INFTY, + HitOrientation orientation = HitOrientation::EXITING, + std::vector* const exclude_primitives = nullptr) + { + fatal_error("GPU ray tracing not supported with this RayTracer backend"); + } + + virtual void ray_fire_prepared(const size_t num_rays, + const double dist_limit = INFTY, + HitOrientation orientation = HitOrientation::EXITING) + { + fatal_error("GPU ray tracing not supported with this RayTracer backend"); + } + + virtual void point_in_volume_prepared(const size_t num_points) + { + fatal_error("GPU ray tracing not supported with this RayTracer backend"); + } + + /** + * @brief Check whether the current ray buffer capacity is sufficient for the number of rays requested + * @param[in] num_rays The number of rays to be processed + */ + virtual void check_rayhit_buffer_capacity(const size_t num_rays) { + fatal_error("GPU ray tracing not supported with this RayTracer backend"); + } + + /** + * @brief return device pointers to ray and hit buffers for GPU ray tracing + * @return DeviceRayHitBuffers struct containing device pointers to ray and hit buffers + */ + virtual DeviceRayHitBuffers get_device_rayhit_buffers(const size_t num_rays) { + fatal_error("GPU ray tracing not supported with this RayTracer backend"); + return {}; + } + + /** + * @brief Allocate device ray buffers and populate them via a user-provided callback + * + * This method allows downstream applications to populate ray buffers using any compute + * backend (GPRT, CUDA, HIP, OpenCL, etc.) without coupling them to XDG's internals. + * + * The workflow: + * 1. XDG allocates device memory for rays (if not already large enough) + * 2. XDG passes device pointers to the callback + * 3. User's callback populates the buffers using their preferred compute kernel/shader + * 4. User's callback returns (XDG assumes buffers are now populated) + * + * This avoids unnecessary host-device transfers by allowing users to write directly + * to XDG's device buffers without any host-side transfers. + * + * @param numRays Number of rays to allocate space for + * @param callback Function that will populate the ray buffer. Receives the allocated buffer and ray count. + */ + virtual void populate_rays_external(size_t numRays, + const RayPopulationCallback& callback) { + fatal_error("GPU ray tracing not supported with this RayTracer backend"); + } + + virtual void transfer_hits_buffer_to_host(const size_t num_rays, + std::vector& hits) { + fatal_error("GPU ray tracing not supported with this RayTracer backend"); + } + protected: // Common functions across RayTracers const double bounding_box_bump(const std::shared_ptr mesh_manager, MeshID volume_id); // return a bump value based on the size of a bounding box (minimum 1e-3). Should this be a part of mesh_manager? @@ -150,4 +307,4 @@ class RayTracer { } // namespace xdg -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/xdg.h b/include/xdg/xdg.h index 82c1f8f2..9dbc8ad2 100644 --- a/include/xdg/xdg.h +++ b/include/xdg/xdg.h @@ -10,6 +10,8 @@ namespace xdg { +struct DeviceRayHitBuffers; // forward declaration +struct dblHit; // forward declaration class XDG { public: @@ -63,11 +65,40 @@ next_element(MeshID current_element, const Position& r, const Direction& u) const; +/** + * @brief Check whether a point lies in a specified volume + * + * This method performs a check to see whether a given point is inside a volume provided. + * It computes this by firing a ray from the point and checking whether or not the ray is Entering or Exiting + * the volume boundary. If no direction is provided, a default direction will be used. + * Note - zero length direction vectors are not explicitly checked for internally and should be avoided to avoid causing undefined behavior. + * + * @param[in] tree The TreeID of the volume we are querying against + * @param[in] point The point to be queried + * @param[in] direction (optional) direction to launch a ray in a specified direction - must be non-zero length + * @param[in] exclude_primitives (optional) vector of surface element MeshIDs to exclude from intersection tests + * @return Boolean result of point in volume check + */ bool point_in_volume(MeshID volume, const Position point, const Direction* direction = nullptr, const std::vector* exclude_primitives = nullptr) const; +/** + * @brief Fire a ray against a given volume and return the first hit + * + * This method fires a ray from a given origin in a specified direction against the surfaces of a volume. + * It returns the distance to the closest hit and the MeshID of the surface hit. The user can specify + * a distance limit and whether Entering/Exiting hits should be rejected. + * + * @param[in] volume The MeshID of the volume we are querying against + * @param[in] origin An array of Position objects representing the starting points of the rays + * @param[in] direction (optional) Direction object to launch a ray in a specified direction + * @param[in] dist_limit (optional) maximum distance to consider for intersections + * @param[in] orientation (optional) flag to consider whether Entering/Exiting hits should be rejected. Defaults to EXITING + * @param[in] exclude_primitives (optional) vector of surface element MeshIDs to exclude from intersection tests + * @return A pair containing the distance to the closest hit and the MeshID of the surface hit + */ std::pair ray_fire(MeshID volume, const Position& origin, const Direction& direction, @@ -75,6 +106,82 @@ std::pair ray_fire(MeshID volume, HitOrientation orientation = HitOrientation::EXITING, std::vector* const exclude_primitives = nullptr) const; +/** + * @brief Array based version of point_in_volume query + * + * This method performs a set of point_in_volume queries on a batch of rays defined by their origins and directions. + * It computes whether or not a point lies in a given volume for each point in the batch. With GPRT ray tracing + * this launches the RT pipeline with the number of rays provided. + * + * @param[in] tree The TreeID of the volume we are querying against + * @param[in] points An array of points to query + * @param[in] num_points The number of points to be processed in the batch + * @param[out] results An output array to store the computed results for each point (1 if inside volume, 0 if outside) + * @param[in] directions (optional) array of directions to launch rays in explicit directions per point - these must be non-zero length + * @param[in] exclude_primitives (optional) vector of surface element MeshIDs to exclude from intersection tests + * @return Void. Outputs stored in results array + */ +void point_in_volume(MeshID volume, + const Position* points, + const size_t num_points, + uint8_t* results, + const Direction* directions = nullptr, + std::vector* exclude_primitives = nullptr) const; + +/** + * @brief Array based version of ray_fire query + * + * This method performs a set of ray fire queries on a batch of rays defined by their origins and directions. + * It computes the intersection distances and surface IDs for each ray in the batch. With GPRT ray tracing + * this launches the RT pipeline with the number of rays provided. + * + * @param[in] tree The TreeID of the volume we are querying against + * @param[in] origins An array of Position objects representing the starting points of the rays + * @param[in] directions An array of Direction objects representing the directions of the rays + * @param[in] num_rays The number of rays to be processed in the batch + * @param[out] hitDistances An output array to store the computed intersection distances for each ray + * @param[out] surfaceIDs An output array to store the MeshIDs of the surfaces hit by each ray + * @param[in] dist_limit (optional) maximum distance to consider for intersections + * @param[in] orientation (optional) flag to consider whether Entering/Exiting hits should be rejected. Defaults to EXITING + * @param[in] exclude_primitives (optional) vector of surface element MeshIDs to exclude from intersection tests + * @return Void. Outputs stored in hitDistances and surfaceIDs arrays + */ +void ray_fire(MeshID volume, + const Position* origins, + const Direction* directions, + const size_t num_rays, + double* hitDistances, + MeshID* surfaceIDs, + const double dist_limit = INFTY, + HitOrientation orientation = HitOrientation::EXITING, + std::vector* const exclude_primitives = nullptr); + +/** + * @brief Call ray fire on pre-populated ray buffers + * + * This method performs a set of ray fire queries on a set of rays that have already been populated on device + * via the external ray population callback method. With GPRT ray tracing this launches the RT pipeline with the number of rays provided. + * + * @param[in] num_rays The number of rays to be processed in the batch + * @param[in] dist_limit (optional) maximum distance to consider for intersections + * @param[in] orientation (optional) flag to consider whether Entering/Exiting hits should be rejected. Defaults to EXITING + * @return Void. Outputs stored in dblHit buffer on device. And can be recovered on host via transfer_hits_buffer_to_host method + */ +void ray_fire_prepared(const size_t num_rays, + const double dist_limit = INFTY, + HitOrientation orientation = HitOrientation::EXITING); + +/** + * @brief Call point_in_volume on pre-populated ray buffers + * + * This method performs a set of point_in_volume queries on a set of points that have already been populated on device + * via the external ray population callback method. With GPRT ray tracing this launches the RT pipeline with the number of points provided. + * + * @param[in] num_points The number of points to be processed in the batch + * @return Void. Outputs stored in dblHit buffer on device. And can be recovered on host via transfer_hits_buffer_to_host method + */ +void point_in_volume_prepared(const size_t num_points); + std::pair closest(MeshID volume, const Position& origin) const; @@ -105,6 +212,23 @@ Direction surface_normal(MeshID surface, ray_tracing_interface_ = ray_tracing_interface; } + DeviceRayHitBuffers get_device_rayhit_buffers(const size_t requiredCapacity) + { + return ray_tracing_interface()->get_device_rayhit_buffers(requiredCapacity); + } + + void populate_rays_external(size_t numRays, + const RayPopulationCallback& callback) + { + return ray_tracing_interface()->populate_rays_external(numRays, callback); + } + +// Device to host transfer of hit buffers (GPRT only for now) +#ifdef XDG_ENABLE_GPRT + void transfer_hits_buffer_to_host(const size_t num_rays, + std::vector& hits); +#endif + // Accessors const std::shared_ptr& ray_tracing_interface() const { return ray_tracing_interface_; @@ -113,6 +237,7 @@ Direction surface_normal(MeshID surface, const std::shared_ptr& mesh_manager() const { return mesh_manager_; } + // Private methods private: double _triangle_volume_contribution(const PrimitiveRef& triangle) const; diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index ff082a45..0ba8a490 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -53,22 +53,27 @@ void ray_fire_miss(inout RayFirePayload payload) { void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { RayFirePayload payload; uint rayID = DispatchRaysIndex().x; + xdg::dblRay ray = record.ray[rayID]; // recover ray from buffer // Trace the ray into the scene RayDesc rayDesc; - rayDesc.Origin = float3(record.ray[rayID].origin); - rayDesc.Direction = normalize(float3(record.ray[rayID].direction)); - rayDesc.TMin = float(record.ray[rayID].tMin); - rayDesc.TMax = float(record.ray[rayID].tMax); + rayDesc.Origin = float3(ray.origin); + rayDesc.Direction = normalize(float3(ray.direction)); + rayDesc.TMin = float(PC.tMin); + rayDesc.TMax = float(PC.tMax); - SurfaceAccelerationStructure world = record.ray[rayID].volume_accel; + int mesh_id = ray.volume_mesh_id; + // Recover the TLAS we are tracing against for this ray + SurfaceAccelerationStructure world = record.meshid_to_accel_address[mesh_id]; - // Pass the ray's origin and direction to the payload + // Set payload default values payload.distance = -1.0f; payload.surf_id = -1; payload.tlas = world; - TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); + if (ray.enabled == 1u) { + TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); + } // Store the distance to the hit point and the surface ID in buffers for CPU record.hit[rayID].distance = payload.distance; @@ -80,22 +85,26 @@ void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { RayFirePayload payload; uint rayID = DispatchRaysIndex().x; + xdg::dblRay ray = record.ray[rayID]; // recover ray from buffer // Trace the ray into the scene RayDesc rayDesc; - rayDesc.Origin = float3(record.ray[rayID].origin); - rayDesc.Direction = float3(normalize(record.ray[rayID].direction)); - rayDesc.TMin = float(record.ray[rayID].tMin); - rayDesc.TMax = float(record.ray[rayID].tMax); - - SurfaceAccelerationStructure world = record.ray[rayID].volume_accel; - + rayDesc.Origin = float3(ray.origin); + rayDesc.Direction = float3(normalize(ray.direction)); + rayDesc.TMin = float(PC.tMin); + rayDesc.TMax = float(PC.tMax); + + int mesh_id = ray.volume_mesh_id; + // Recover the TLAS we are tracing against for this ray + SurfaceAccelerationStructure world = record.meshid_to_accel_address[mesh_id]; // Pass the ray's origin and direction to the payload payload.surf_id = -1; payload.tlas = world; payload.piv = xdg::PointInVolume::OUTSIDE; // Initialize point in volume check result to outside (0) - TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); + if (ray.enabled == 1u) { // skip RT pipeline for rays that are marked as disabled + TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); + } record.hit[rayID].surf_id = payload.surf_id; record.hit[rayID].piv = payload.piv; // Point in volume check result @@ -104,10 +113,14 @@ void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData me // ------------------------------------------------- Compute Shaders ------------------------------------------------- /* A shader to compute and store AABB min/maxes in single precision using double precision coords*/ [shader("compute")] -[numthreads(1, 1, 1)] -void -populate_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGeomData record) { +[numthreads(64, 1, 1)] +void populate_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGeomData record) { int primID = DispatchThreadID.x; + + // guard against more threads than primitives + if (primID >= record.num_faces) + return; + int3 indices = record.index[primID]; double3 A = record.vertex[indices[0]]; double3 B = record.vertex[indices[1]]; @@ -123,7 +136,6 @@ populate_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGe // ------------------------------------------------ CUSTOM INTERSECTION SHADERS ------------------------------------------------ - /* 1D ray generation intersection with a double precision triangle using the Plucker intersection algorithm*/ [shader("intersection")] void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) @@ -134,6 +146,7 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) uint rayID = DispatchRaysIndex().x; uint nRays = DispatchRaysDimensions().x; uint flags = RayFlags(); + xdg::dblRay ray = record.ray[rayID]; // recover ray from buffer if (rayID >= nRays) { return; @@ -155,11 +168,11 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) double3 v1 = record.vertex[indices[1]]; double3 v2 = record.vertex[indices[2]]; - double3 origin = record.ray[rayID].origin; - double3 direction = record.ray[rayID].direction; + double3 origin = ray.origin; + double3 direction = ray.direction; - double tMin = record.ray[rayID].tMin; - double tMax = record.ray[rayID].tMax; + double tMin = PC.tMin; + double tMax = PC.tMax; const double3 raya = direction; const double3 rayb = cross(direction, origin); @@ -228,8 +241,9 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) double3 norm = record.normals[primID]; // recover double precision normal. TODO - Should we calculate from vertices instead? - // sense adjustment of normal - if (record.ray[rayID].volume_tree == record.reverse_tree) + // sense adjustment of normal (per MeshID) + int mesh_sense = record.meshid_to_sense[ray.volume_mesh_id]; // +1 forward, -1 reverse + if (mesh_sense < 0) { norm = -norm; } @@ -238,16 +252,16 @@ void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) uint hit_kind = norm_dot_dir < 0 ? HIT_KIND_TRIANGLE_FRONT_FACE : HIT_KIND_TRIANGLE_BACK_FACE; - xdg::HitOrientation hitOrientation = record.ray[rayID].hitOrientation; + xdg::HitOrientation hitOrientation = PC.hitOrientation; if (orientation_cull(direction, norm, hitOrientation)) { return; } - for (int i = 0; i < record.ray[rayID].exclude_count; ++i) + for (int i = 0; i < ray.exclude_count; ++i) { - if (record.ray[rayID].exclude_primitives[i] == global_prim_id) { + if (ray.exclude_primitives[i] == global_prim_id) { return; } } @@ -312,4 +326,4 @@ float next_after(float a) { a_++; } return asfloat(a_); -} \ No newline at end of file +} diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index bb3ece22..1d76283e 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -1,22 +1,27 @@ #include "xdg/gprt/ray_tracer.h" #include "gprt/gprt.h" - namespace xdg { GPRTRayTracer::GPRTRayTracer() { + gprtRequestRayTypeCount(numRayTypes_); // Set the number of shaders which can be set to the same geometry context_ = gprtContextCreate(); module_ = gprtModuleCreate(context_, dbl_deviceCode); - rayHitBuffers_.capacity = 1; // Preallocate space for 1 ray - rayHitBuffers_.ray = gprtDeviceBufferCreate(context_, rayHitBuffers_.capacity); - rayHitBuffers_.hit = gprtDeviceBufferCreate(context_, rayHitBuffers_.capacity); - rayHitBuffers_.devRayAddr = gprtBufferGetDevicePointer(rayHitBuffers_.ray); - rayHitBuffers_.devHitAddr = gprtBufferGetDevicePointer(rayHitBuffers_.hit); + // Buffer setup + rayHitBuffers_.view.capacity = 1e6; // Preallocate space for 1m rays + rayHitBuffers_.ray = gprtDeviceBufferCreate(context_, rayHitBuffers_.view.capacity); + rayHitBuffers_.hit = gprtDeviceBufferCreate(context_, rayHitBuffers_.view.capacity); + rayHitBuffers_.view.rayDevPtr = gprtBufferGetDevicePointer(rayHitBuffers_.ray); + rayHitBuffers_.view.hitDevPtr = gprtBufferGetDevicePointer(rayHitBuffers_.hit); + rayHitBuffers_.view.rayStride = sizeof(dblRay); + rayHitBuffers_.view.hitStride = sizeof(dblHit); excludePrimitivesBuffer_ = gprtDeviceBufferCreate(context_); // initialise buffer of size 1 + tlas_handle_buffer_ = gprtDeviceBufferCreate(context_); // initialise buffer of size 1 + setup_shaders(); @@ -32,6 +37,8 @@ GPRTRayTracer::GPRTRayTracer() // Set up build parameters for acceleration structures buildParams_.buildMode = GPRT_BUILD_MODE_FAST_BUILD_NO_UPDATE; + + } GPRTRayTracer::~GPRTRayTracer() @@ -46,11 +53,6 @@ GPRTRayTracer::~GPRTRayTracer() gprtAccelDestroy(accel); } - // Destroy BLAS structures - for (const auto& blas : blas_handles_) { - gprtAccelDestroy(blas); - } - // Destroy Geoms and Types for (const auto& [surf, geom] : surface_to_geometry_map_) { gprtGeomDestroy(geom); @@ -61,6 +63,7 @@ GPRTRayTracer::~GPRTRayTracer() gprtBufferDestroy(rayHitBuffers_.ray); gprtBufferDestroy(rayHitBuffers_.hit); gprtBufferDestroy(excludePrimitivesBuffer_); + gprtBufferDestroy(tlas_handle_buffer_); // Destroy module and context gprtModuleDestroy(module_); @@ -85,9 +88,13 @@ void GPRTRayTracer::setup_shaders() void GPRTRayTracer::init() { + update_tlas_table_(); + // Build the shader binding table (SBT) after all shader programs and acceleration structures are set up gprtBuildShaderBindingTable(context_, GPRT_SBT_ALL); // Note that should we need to update any shaders or acceleration structures, we must rebuild the SBT + + initialized_ = true; } std::pair @@ -164,8 +171,12 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana geom_data->normals = gprtBufferGetDevicePointer(normal_buffer); geom_data->primitive_refs = gprtBufferGetDevicePointer(primitive_refs_buffer); geom_data->num_faces = num_faces; - - gprtComputeLaunch(aabbPopulationProgram_, {num_faces, 1, 1}, {1, 1, 1}, *geom_data); + // meshid_to_sense pointer is set after meshid_to_sense_buffer_ is created + + constexpr uint32_t threadsPerGroup = 64; // must match [numthreads(64,1,1)] + uint32_t numGroupsX = (num_faces + threadsPerGroup - 1) / threadsPerGroup; + + gprtComputeLaunch(aabbPopulationProgram_, {numGroupsX, 1, 1}, {threadsPerGroup, 1, 1}, *geom_data); GPRTAccel blas = gprtAABBAccelCreate(context_, triangleGeom, buildParams_.buildMode); @@ -184,14 +195,14 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana surfaceBlasInstances.push_back(instance); globalBlasInstances_.push_back(instance); - // Always update per-volume info + // Always update per-volume info and MeshID -> sparse sense mapping auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); if (volume_id == forward_parent) { - geom_data->forward_vol = forward_parent; - geom_data->forward_tree = tree; + meshid_to_sense_.resize(static_cast(forward_parent) + 1, 1); + meshid_to_sense_[forward_parent] = 1; } else if (volume_id == reverse_parent) { - geom_data->reverse_vol = reverse_parent; - geom_data->reverse_tree = tree; + meshid_to_sense_.resize(static_cast(reverse_parent) + 1, 1); + meshid_to_sense_[reverse_parent] = -1; } else { fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); } @@ -202,7 +213,19 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana GPRTAccel volume_tlas = gprtInstanceAccelCreate(context_, surfaceBlasInstances.size(), instanceBuffer); gprtAccelBuild(context_, volume_tlas, buildParams_); surface_volume_tree_to_accel_map[tree] = volume_tlas; - + surface_tree_to_volume_map_[tree] = volume_id; + if (volume_id >= tlas_handles_.size()) { + tlas_handles_.resize(volume_id + 1, SurfaceAccelerationStructure{}); + } + tlas_handles_[volume_id] = gprtAccelGetDeviceAddress(volume_tlas); + + update_meshid_to_sense_(); + + if (initialized_) { + update_tlas_table_(); + gprtBuildShaderBindingTable(context_, GPRT_SBT_ALL); + } + return tree; } @@ -222,19 +245,24 @@ bool GPRTRayTracer::point_in_volume(SurfaceTreeID tree, auto rayGen = rayGenPrograms_.at(RayGenType::POINT_IN_VOLUME); dblRayGenData* rayGenPIVData = gprtRayGenGetParameters(rayGen); + const Direction defaultDir = Direction{1. / std::sqrt(2.0), 1. / std::sqrt(2.0), 0.0}; + // 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}; + : defaultDir; + + // Catch directions with zero length + const double l2 = directionUsed.x*directionUsed.x + + directionUsed.y*directionUsed.y + + directionUsed.z*directionUsed.z; + if (l2 == 0.0) directionUsed = defaultDir; gprtBufferMap(rayHitBuffers_.ray); // Update the ray input buffer dblRay* ray = gprtBufferGetHostPointer(rayHitBuffers_.ray); - ray[0].volume_accel = gprtAccelGetDeviceAddress(volume); ray[0].origin = {point.x, point.y, point.z}; ray[0].direction = {directionUsed.x, directionUsed.y, directionUsed.z}; - ray[0].tMax = INFTY; // Set a large distance limit - ray[0].tMin = 0.0; - ray[0].volume_tree = tree; // Set the TreeID of the volume being queried - ray[0].hitOrientation = HitOrientation::ANY; // No orientation culling for point-in-volume check + ray[0].volume_mesh_id = surface_tree_to_volume_map_.at(tree); + ray[0].enabled = 1; // Ensure the ray is enabled if (exclude_primitives) { if (!exclude_primitives->empty()) gprtBufferResize(context_, excludePrimitivesBuffer_, exclude_primitives->size(), false); @@ -252,7 +280,12 @@ bool GPRTRayTracer::point_in_volume(SurfaceTreeID tree, } gprtBufferUnmap(rayHitBuffers_.ray); // required to sync buffer back on GPU? - gprtRayGenLaunch1D(context_, rayGen, 1); // Launch raygen shader (entry point to RT pipeline) + dblRayFirePushConstants pushConstants; + pushConstants.hitOrientation = HitOrientation::ANY; + pushConstants.tMax = INFTY; + pushConstants.tMin = 0.0; + + gprtRayGenLaunch1D(context_, rayGen, 1, pushConstants); // Launch raygen shader (entry point to RT pipeline) gprtGraphicsSynchronize(context_); // Ensure all GPU operations are complete before returning control flow to CPU // Retrieve the hit from the dblHit buffer @@ -282,16 +315,13 @@ std::pair GPRTRayTracer::ray_fire(SurfaceTreeID tree, GPRTAccel volume = surface_volume_tree_to_accel_map.at(tree); auto rayGen = rayGenPrograms_.at(RayGenType::RAY_FIRE); dblRayGenData* rayGenData = gprtRayGenGetParameters(rayGen); - + gprtBufferMap(rayHitBuffers_.ray); // Update the ray input buffer dblRay* ray = gprtBufferGetHostPointer(rayHitBuffers_.ray); - ray[0].volume_accel = gprtAccelGetDeviceAddress(volume); ray[0].origin = {origin.x, origin.y, origin.z}; ray[0].direction = {direction.x, direction.y, direction.z}; - ray[0].tMax = dist_limit; - ray[0].tMin = 0.0; - ray[0].hitOrientation = orientation; // Set orientation for the ray - ray[0].volume_tree = tree; // Set the TreeID of the volume being queried + ray[0].volume_mesh_id = surface_tree_to_volume_map_.at(tree); + ray[0].enabled = 1; // Ensure the ray is enabled if (exclude_primitives) { if (!exclude_primitives->empty()) gprtBufferResize(context_, excludePrimitivesBuffer_, exclude_primitives->size(), false); @@ -309,7 +339,13 @@ std::pair GPRTRayTracer::ray_fire(SurfaceTreeID tree, } gprtBufferUnmap(rayHitBuffers_.ray); // required to sync buffer back on GPU? - gprtRayGenLaunch1D(context_, rayGen, 1); // Launch raygen shader (entry point to RT pipeline) + // Set push constants (same for every ray) + dblRayFirePushConstants pushConstants; + pushConstants.hitOrientation = orientation; + pushConstants.tMax = dist_limit; + pushConstants.tMin = 0.0; + + gprtRayGenLaunch1D(context_, rayGen, 1, pushConstants); // Launch raygen shader (entry point to RT pipeline) gprtGraphicsSynchronize(context_); // Ensure all GPU operations are complete before returning control flow to CPU // Retrieve the hit from the dblHit buffer @@ -326,7 +362,179 @@ std::pair GPRTRayTracer::ray_fire(SurfaceTreeID tree, if (exclude_primitives) exclude_primitives->push_back(primitive_id); return {distance, surface}; } - + +void GPRTRayTracer::point_in_volume(TreeID tree, + const Position* points, + const size_t num_points, + uint8_t* results, + const Direction* directions, + std::vector* exclude_primitives) +{ + if (num_points == 0) { + warning("Warning number of points passed to point_in_volume is 0. No work to be done."); + return; + } + + GPRTAccel volume = surface_volume_tree_to_accel_map.at(tree); + auto rayGen = rayGenPrograms_.at(RayGenType::POINT_IN_VOLUME); + dblRayGenData* rayGenData = gprtRayGenGetParameters(rayGen); + check_rayhit_buffer_capacity(num_points); + + // TODO - handle exclude_primitives for batch version + + // Set a default direction to be used if no direction is provided + const Direction defaultDir = Direction{1. / std::sqrt(2.0), 1. / std::sqrt(2.0), 0.0}; + + // Map the region start + gprtBufferMap(rayHitBuffers_.ray); + dblRay* ray = gprtBufferGetHostPointer(rayHitBuffers_.ray); + + // Common ray params + for (size_t i = 0; i < num_points; ++i) { + ray[i].origin = {points[i].x, points[i].y, points[i].z}; + ray[i].exclude_primitives = nullptr; + ray[i].volume_mesh_id = surface_tree_to_volume_map_.at(tree); + ray[i].enabled = 1; // Ensure the ray is enabled + } + + // Directions + if (!directions) { + for (size_t i = 0; i < num_points; ++i) + ray[i].direction = double3{ defaultDir.x, defaultDir.y, defaultDir.z }; + } else { + for (size_t i = 0; i < num_points; ++i) + ray[i].direction = double3{ directions[i].x, directions[i].y, directions[i].z }; + } + + gprtBufferUnmap(rayHitBuffers_.ray); // required to sync buffer back on GPU? + + // Set push constants (same for every ray) + dblRayFirePushConstants pushConstants; + pushConstants.hitOrientation = HitOrientation::ANY; + pushConstants.tMax = INFTY; + pushConstants.tMin = 0.0; + + gprtRayGenLaunch1D(context_, rayGen, num_points, pushConstants); + gprtGraphicsSynchronize(context_); + + // Retrieve the output from the ray output buffer + gprtBufferMap(rayHitBuffers_.hit); + dblHit* hit = gprtBufferGetHostPointer(rayHitBuffers_.hit); + for (size_t i = 0; i < num_points; ++i) { + auto piv = hit[i].piv; // Point in volume check result + results[i] = static_cast(piv); + } + gprtBufferUnmap(rayHitBuffers_.hit); // required to sync buffer back on GPU? Maybe this second unmap isn't actually needed since we dont need to resyncrhonize after retrieving the data from device + + return; +} + +// Array version of ray_fire +void GPRTRayTracer::ray_fire(TreeID tree, + const Position* origins, + const Direction* directions, + const size_t num_rays, + double* hitDistances, + MeshID* surfaceIDs, + const double dist_limit, + HitOrientation orientation, + std::vector* const exclude_primitives) +{ + if (num_rays == 0) { + warning("Warning number of rays passed to ray_fire is 0. No work to be done."); + return; + } + + GPRTAccel volume = surface_volume_tree_to_accel_map.at(tree); + auto rayGen = rayGenPrograms_.at(RayGenType::RAY_FIRE); + dblRayGenData* rayGenData = gprtRayGenGetParameters(rayGen); + check_rayhit_buffer_capacity(num_rays); + + gprtBufferMap(rayHitBuffers_.ray); + dblRay* ray = gprtBufferGetHostPointer(rayHitBuffers_.ray); + // Set per ray values + for (size_t i = 0; i < num_rays; ++i) { + const auto& origin = origins[i]; + const auto& direction = directions[i]; + + ray[i].origin = {origin.x, origin.y, origin.z}; + ray[i].direction = {direction.x, direction.y, direction.z}; + ray[i].exclude_primitives = nullptr; // Not currently supported in batch version + ray[i].volume_mesh_id = surface_tree_to_volume_map_.at(tree); + ray[i].enabled = 1; // Ensure the ray is enabled + } + + gprtBufferUnmap(rayHitBuffers_.ray); // required to sync buffer back on GPU? + + // Set push constants (same for every ray) + dblRayFirePushConstants pushConstants; + pushConstants.hitOrientation = orientation; + pushConstants.tMax = dist_limit; + pushConstants.tMin = 0.0; + + // Launch the ray generation shader with push constants and buffer bindings + gprtRayGenLaunch1D(context_, rayGen, num_rays, pushConstants); + gprtGraphicsSynchronize(context_); + + // Retrieve the output from the ray output buffer + gprtBufferMap(rayHitBuffers_.hit); + dblHit* hit = gprtBufferGetHostPointer(rayHitBuffers_.hit); + // populate the result arrays + for (size_t i = 0; i < num_rays; ++i) { + const MeshID surfaceHit = hit[i].surf_id; + if (surfaceHit == ID_NONE) { + hitDistances[i] = INFTY; + surfaceIDs[i] = ID_NONE; + } + else { + hitDistances[i] = hit[i].distance; + surfaceIDs[i] = surfaceHit; + // TODO - handle exclude_primitives for batch version + } + } + gprtBufferUnmap(rayHitBuffers_.hit); // required to sync buffer back on GPU? Maybe this second unmap isn't actually needed since we dont need to resyncrhonize after retrieving the data from device + return; +} + +void +GPRTRayTracer::ray_fire_prepared(const size_t num_rays, + const double dist_limit, + HitOrientation orientation) +{ + if (num_rays == 0) return; // no work to do. Early exit + + check_rayhit_buffer_capacity(num_rays); + + auto rayGen = rayGenPrograms_.at(RayGenType::RAY_FIRE); + + dblRayFirePushConstants pushConstants; + pushConstants.tMax = dist_limit; + pushConstants.tMin = 0.0; + pushConstants.hitOrientation = orientation; // Set orientation for the ray + + gprtRayGenLaunch1D(context_, rayGen, num_rays, pushConstants); + gprtGraphicsSynchronize(context_); + return; +} + +void +GPRTRayTracer::point_in_volume_prepared(const size_t num_points) +{ + if (num_points == 0) return; // no work to do. Early exit + + check_rayhit_buffer_capacity(num_points); + auto rayGen = rayGenPrograms_.at(RayGenType::POINT_IN_VOLUME); + + dblRayFirePushConstants pushConstants; + pushConstants.tMax = INFTY; + pushConstants.tMin = 0.0; + pushConstants.hitOrientation = HitOrientation::ANY; // Set orientation for the ray + + gprtRayGenLaunch1D(context_, rayGen, num_points, pushConstants); + gprtGraphicsSynchronize(context_); + return; +} + void GPRTRayTracer::create_global_surface_tree() { // Create a TLAS (Top-Level Acceleration Structure) for all the volumes @@ -341,16 +549,22 @@ void GPRTRayTracer::create_global_surface_tree() global_surface_accel_ = global_accel; } -void GPRTRayTracer::check_ray_buffer_capacity(size_t N) +void GPRTRayTracer::check_rayhit_buffer_capacity(const size_t N) { - if (N <= rayHitBuffers_.capacity) return; // current capacity is sufficient + if (N <= rayHitBuffers_.view.capacity) return; // current capacity is sufficient - // Resize buffers to accommodate N rays - size_t newCapacity = std::max(N, rayHitBuffers_.capacity * 2); // double the capacity or set to N, whichever is larger + // Resize buffers to accommodate N rays - double the capacity or set to N, whichever is larger + size_t newCapacity = std::max(N, rayHitBuffers_.view.capacity * 2); gprtBufferResize(context_, rayHitBuffers_.ray, newCapacity, false); gprtBufferResize(context_, rayHitBuffers_.hit, newCapacity, false); - rayHitBuffers_.capacity = newCapacity; + rayHitBuffers_.view.capacity = newCapacity; + + // Get fresh device pointers after resize + rayHitBuffers_.view.rayDevPtr = gprtBufferGetDevicePointer(rayHitBuffers_.ray); + rayHitBuffers_.view.hitDevPtr = gprtBufferGetDevicePointer(rayHitBuffers_.hit); + rayHitBuffers_.view.rayStride = sizeof(dblRay); + rayHitBuffers_.view.hitStride = sizeof(dblHit); // Since we have resized the ray buffers, we need to update the geom_data->rayIn pointers in all geometries too for (auto const& [surf, geom] : surface_to_geometry_map_) { @@ -363,9 +577,75 @@ void GPRTRayTracer::check_ray_buffer_capacity(size_t N) dblRayGenData* rayGenData = gprtRayGenGetParameters(rayGen); rayGenData->ray = gprtBufferGetDevicePointer(rayHitBuffers_.ray); rayGenData->hit = gprtBufferGetDevicePointer(rayHitBuffers_.hit); + rayGenData->meshid_to_accel_address = gprtBufferGetDevicePointer(tlas_handle_buffer_); } gprtBuildShaderBindingTable(context_, static_cast(GPRT_SBT_GEOM | GPRT_SBT_RAYGEN)); } +// Update the TLAS table (MeshID -> SurfaceAccelerationStructure) buffer on the device +void GPRTRayTracer::update_tlas_table_() +{ + upload_device_buffer_(tlas_handle_buffer_, tlas_handles_); + + for (auto type : {RayGenType::RAY_FIRE, RayGenType::POINT_IN_VOLUME}) { + auto* raygendata = gprtRayGenGetParameters(rayGenPrograms_.at(type)); + raygendata->meshid_to_accel_address = gprtBufferGetDevicePointer(tlas_handle_buffer_); + } +} + +void GPRTRayTracer::update_meshid_to_sense_() +{ + upload_device_buffer_(meshid_to_sense_buffer_, meshid_to_sense_); + + for (auto const& [surf, geom] : surface_to_geometry_map_) { + DPTriangleGeomData* geom_data = gprtGeomGetParameters(geom); + geom_data->meshid_to_sense = gprtBufferGetDevicePointer(meshid_to_sense_buffer_); + } +} + + +DeviceRayHitBuffers GPRTRayTracer::get_device_rayhit_buffers(const size_t N) +{ + check_rayhit_buffer_capacity(N); + return rayHitBuffers_.view; +} + +void GPRTRayTracer::populate_rays_external(size_t numRays, + const RayPopulationCallback& callback) +{ + if (numRays == 0) { + warning("Warning number of rays passed to populate_rays_external is 0. No work to be done."); + return; + } + + // Ensure device buffers are large enough + check_rayhit_buffer_capacity(numRays); + + // Use the user callback to populate the rays directly on the device + callback(rayHitBuffers_.view, numRays); + + // After callback returns, we assume the ray buffer is populated and ready to trace + // Note: The callback is responsible for synchronization if using an async API +} + +void GPRTRayTracer::transfer_hits_buffer_to_host(const size_t num_rays, + std::vector& hits) +{ + hits.clear(); // Ensure hits vector is empty before populating + if (num_rays == 0) { + warning("Warning number of rays passed to transfer_hits_buffer_to_host is 0. No work to be done."); + return; + } + if (num_rays > rayHitBuffers_.view.capacity) { + fatal_error("Requested {} hits, but hit buffer capacity is {}", num_rays, rayHitBuffers_.view.capacity); + } + + hits.resize(num_rays); + gprtBufferMap(rayHitBuffers_.hit); + dblHit* hit = gprtBufferGetHostPointer(rayHitBuffers_.hit); + std::copy(hit, hit + num_rays, hits.begin()); + gprtBufferUnmap(rayHitBuffers_.hit); +} + } // namespace xdg diff --git a/src/tetrahedron_contain.cpp b/src/tetrahedron_contain.cpp index 743ea1d0..98757768 100644 --- a/src/tetrahedron_contain.cpp +++ b/src/tetrahedron_contain.cpp @@ -13,7 +13,9 @@ bool plucker_tet_containment_test(const Position& point, const Position& v1, const Position& v2, const Position& v3) { - using namespace linalg::aliases; + using linalg::aliases::double3x3; + using linalg::aliases::double3; + using linalg::aliases::double4; // Create matrix T = [v1 - v0, v2 - v0, v3 - v0] Vec3da e0 = v1 - v0; Vec3da e1 = v2 - v0; diff --git a/src/xdg.cpp b/src/xdg.cpp index 371a9a1c..8bc8a2ec 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -52,6 +52,14 @@ void XDG::prepare_volume_for_raytracing(MeshID volume) { volume_to_point_location_tree_map_[volume] = volume_tree; } +#ifdef XDG_ENABLE_GPRT +void XDG::transfer_hits_buffer_to_host(const size_t num_rays, + std::vector& hits) +{ + ray_tracing_interface()->transfer_hits_buffer_to_host(num_rays, hits); +} +#endif + std::shared_ptr XDG::create(MeshLibrary mesh_lib, RTLibrary ray_tracing_lib) { std::shared_ptr xdg = std::make_shared(); @@ -110,6 +118,17 @@ bool XDG::point_in_volume(MeshID volume, return ray_tracing_interface()->point_in_volume(tree, point, direction, exclude_primitives); } +void XDG::point_in_volume(MeshID volume, + const Position* points, + const size_t num_points, + uint8_t* results, + const Direction* directions, + std::vector* exclude_primitives) const +{ + TreeID tree = volume_to_surface_tree_map_.at(volume); + ray_tracing_interface()->point_in_volume(tree, points, num_points, results, directions, exclude_primitives); +} + MeshID XDG::find_volume(const Position& point, const Direction& direction) const { @@ -235,8 +254,38 @@ XDG::ray_fire(MeshID volume, HitOrientation orientation, std::vector* const exclude_primitives) const { - TreeID scene = volume_to_surface_tree_map_.at(volume); - return ray_tracing_interface()->ray_fire(scene, origin, direction, dist_limit, orientation, exclude_primitives); + TreeID tree = volume_to_surface_tree_map_.at(volume); + return ray_tracing_interface()->ray_fire(tree, origin, direction, dist_limit, orientation, exclude_primitives); +} + +// Array version of ray_fire +void +XDG::ray_fire(MeshID volume, + const Position* origins, + const Direction* directions, + const size_t num_rays, + double* hitDistances, + MeshID* surfaceIDs, + const double dist_limit, + HitOrientation orientation, + std::vector* const exclude_primitives) +{ + TreeID tree = volume_to_surface_tree_map_.at(volume); + return ray_tracing_interface()->ray_fire(tree, origins, directions, num_rays, hitDistances, surfaceIDs, dist_limit, orientation, exclude_primitives); +} + +void +XDG::ray_fire_prepared(const size_t num_rays, + const double dist_limit, + HitOrientation orientation) +{ + return ray_tracing_interface()->ray_fire_prepared(num_rays, dist_limit, orientation); +} + +void +XDG::point_in_volume_prepared(const size_t num_points) +{ + return ray_tracing_interface()->point_in_volume_prepared(num_points); } std::pair XDG::closest(MeshID volume, @@ -325,4 +374,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 d7286555..8166eee9 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -37,6 +37,11 @@ if (XDG_ENABLE_MOAB) list(APPEND TEST_NAMES test_overlap_check) endif() +# This test really should be appended when any GPU library is enabled +if (XDG_ENABLE_GPRT) + list(APPEND TEST_NAMES test_direct_ray_buffer_access) +endif() + foreach(test ${TEST_NAMES}) add_executable(${test} ${test}.cpp) target_link_libraries(${test} xdg Catch2::Catch2WithMain) @@ -48,6 +53,18 @@ foreach(test ${TEST_NAMES}) TEST_PREFIX "${test}::") endforeach() +if (XDG_ENABLE_GPRT) + embed_devicecode( + OUTPUT_TARGET + test_direct_ray_buffer_access_deviceCode + HEADERS + ${CMAKE_CURRENT_SOURCE_DIR}/test_direct_ray_buffer_access_shared.h + SOURCES + ${CMAKE_CURRENT_SOURCE_DIR}/test_direct_ray_buffer_access_deviceCode.slang + ) + target_link_libraries(test_direct_ray_buffer_access test_direct_ray_buffer_access_deviceCode) +endif() + set( TEST_FILES diff --git a/tests/test_direct_ray_buffer_access.cpp b/tests/test_direct_ray_buffer_access.cpp new file mode 100644 index 00000000..97a96a11 --- /dev/null +++ b/tests/test_direct_ray_buffer_access.cpp @@ -0,0 +1,201 @@ +// for testing +#include +#include +#include +#include + +// xdg includes +#include "xdg/constants.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/gprt/ray_tracer.h" +#include "xdg/xdg.h" +#include "mesh_mock.h" +#include "test_direct_ray_buffer_access_shared.h" +#include "util.h" +#include "gprt.h" + +#include + +using namespace xdg; +using namespace xdg::test; + +extern GPRTProgram test_direct_ray_buffer_access_deviceCode; + +static RayPopulationCallback make_populate_callback(const std::vector& origins, + const std::vector& directions, + const std::vector& volume_ids, + GPRTContext context, + GPRTComputeOf packRays) { + return [&origins, &directions, volume_ids, context, packRays] + (const DeviceRayHitBuffers& buffer, size_t num_rays) { + + // When passing arrays to the callback, ensure they are the correct size + assert(origins.size() == num_rays); + assert(directions.size() == num_rays); + + // Convert to double3 for use on GPU + std::vector origins_device(num_rays); + std::vector directions_device(num_rays); + for (size_t i = 0; i < num_rays; ++i) { + origins_device[i] = {origins[i].x, origins[i].y, origins[i].z}; + directions_device[i] = {directions[i].x, directions[i].y, directions[i].z}; + } + + auto origins_buffer = gprtDeviceBufferCreate(context, num_rays, origins_device.data()); + auto directions_buffer = gprtDeviceBufferCreate(context, num_rays, directions_device.data()); + auto volume_ids_buffer = gprtDeviceBufferCreate(context, num_rays, volume_ids.data()); + + constexpr uint32_t threads_per_group = 256; + const uint32_t groups = static_cast((num_rays + threads_per_group - 1) / threads_per_group); + + ExternalRayParams params = {}; + params.xdgRays = static_cast(buffer.rayDevPtr); + params.origins = gprtBufferGetDevicePointer(origins_buffer); + params.directions = gprtBufferGetDevicePointer(directions_buffer); + params.num_rays = static_cast(num_rays); + params.total_threads = groups * threads_per_group; + params.volume_mesh_ids = gprtBufferGetDevicePointer(volume_ids_buffer); // Pass array of volume IDs to compute shader + params.enabled = 1u; + + gprtComputeLaunch(packRays, + { groups, 1, 1 }, + { threads_per_group, 1, 1 }, + params); + gprtComputeSynchronize(context); + + gprtBufferDestroy(origins_buffer); + gprtBufferDestroy(directions_buffer); + if (volume_ids_buffer) { + gprtBufferDestroy(volume_ids_buffer); + } + }; +} + +// This is a GPU only test - skip if no GPU ray tracing backends are enabled +TEMPLATE_TEST_CASE("Ray Fire with external populated rays", "[rayfire][mock]", + GPRT_Raytracer) +{ + constexpr auto rt_backend = TestType::value; + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB, rt_backend); + REQUIRE(xdg->ray_tracing_interface()->library() == rt_backend); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MOAB); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("jezebel.h5m"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + std::vector origins; + std::vector directions; + size_t N = 64; + make_rays(N, origins, directions); + + auto gprt_rt = std::dynamic_pointer_cast(xdg->ray_tracing_interface()); + REQUIRE(gprt_rt); + + std::vector volumes = mesh_manager->volumes(); + REQUIRE(volumes.size() >= 2); + const MeshID volume_id = volumes[0]; + const MeshID volume_id_alt = volumes[1]; + GPRTContext context = gprt_rt->context(); + GPRTModule module = gprtModuleCreate(context, test_direct_ray_buffer_access_deviceCode); + auto packRays = gprtComputeCreate(context, module, "pack_external_rays"); + + std::vector expected_distances(N, INFTY); + std::vector expected_surfaces(N, ID_NONE); + std::vector volume_ids(N, volume_id); + for (size_t i = 0; i < N; ++i) { + volume_ids[i] = (i % 2 == 0) ? volume_id : volume_id_alt; // Volume IDs alternating between two volumes + const auto [dist, surf] = xdg->ray_fire(volume_ids[i], origins[i], directions[i]); + expected_distances[i] = dist; + expected_surfaces[i] = surf; + } + + // Create callback to populate rays on device + RayPopulationCallback populate_callback = make_populate_callback(origins, + directions, + volume_ids, + context, + packRays); + + // Populate rays via external API + xdg->populate_rays_external(N, populate_callback); + + xdg->ray_fire_prepared(N); + std::vector hits; + xdg->transfer_hits_buffer_to_host(N, hits); + + REQUIRE(hits.size() == N); + for (size_t i = 0; i < N; ++i) { + REQUIRE(hits[i].surf_id == expected_surfaces[i]); + if (expected_surfaces[i] != ID_NONE) { + REQUIRE_THAT(hits[i].distance, Catch::Matchers::WithinAbs(expected_distances[i], 1e-6)); + } + } + + gprtComputeDestroy(packRays); + gprtModuleDestroy(module); + } +} + +TEMPLATE_TEST_CASE("Point-in-volume with external populated rays", "[piv][mock]", + GPRT_Raytracer) +{ + constexpr auto rt_backend = TestType::value; + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB, rt_backend); + REQUIRE(xdg->ray_tracing_interface()->library() == rt_backend); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MOAB); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("jezebel.h5m"); + mesh_manager->init(); + xdg->prepare_raytracer(); + + std::vector points; + std::vector directions; + size_t N = 64; + make_points(N, points, directions); + + auto gprt_rt = std::dynamic_pointer_cast(xdg->ray_tracing_interface()); + REQUIRE(gprt_rt); + + std::vector volumes = mesh_manager->volumes(); + REQUIRE(volumes.size() >= 2); + const MeshID volume_id = volumes[0]; + const MeshID volume_id_alt = volumes[1]; + GPRTContext context = gprt_rt->context(); + GPRTModule module = gprtModuleCreate(context, test_direct_ray_buffer_access_deviceCode); + auto packRays = gprtComputeCreate(context, module, "pack_external_rays"); + + std::vector expected_piv(N, 0); + std::vector volume_ids(N, volume_id); + for (size_t i = 0; i < N; ++i) { + volume_ids[i] = (i % 2 == 0) ? volume_id : volume_id_alt; + expected_piv[i] = static_cast(xdg->point_in_volume(volume_ids[i], points[i], &directions[i])); + } + + RayPopulationCallback populate_callback = make_populate_callback(points, + directions, + volume_ids, + context, + packRays); + xdg->populate_rays_external(N, populate_callback); + + xdg->point_in_volume_prepared(N); + std::vector hits; + xdg->transfer_hits_buffer_to_host(N, hits); + + REQUIRE(hits.size() == N); + for (size_t i = 0; i < N; ++i) { + const auto expected = expected_piv[i] ? xdg::PointInVolume::INSIDE : xdg::PointInVolume::OUTSIDE; // convert back to enum + REQUIRE(hits[i].piv == expected); + } + + gprtComputeDestroy(packRays); + gprtModuleDestroy(module); + } +} diff --git a/tests/test_direct_ray_buffer_access_deviceCode.slang b/tests/test_direct_ray_buffer_access_deviceCode.slang new file mode 100644 index 00000000..b61e6126 --- /dev/null +++ b/tests/test_direct_ray_buffer_access_deviceCode.slang @@ -0,0 +1,24 @@ +#include "test_direct_ray_buffer_access_shared.h" + +[shader("compute")] +[numthreads(256, 1, 1)] +void pack_external_rays(uint3 DispatchThreadID: SV_DispatchThreadID, + uniform ExternalRayParams extParams) +{ + uint globalThreadID = DispatchThreadID.x; + uint stride = extParams.total_threads; // Groups * 256 + + // Grid-stride loop: each thread handles ray idx, idx+stride, idx+2*stride, ... + for (uint idx = globalThreadID; idx < extParams.num_rays; idx += stride) + { + xdg::dblRay r; + r.origin = extParams.origins[idx]; + r.direction = extParams.directions[idx]; + r.exclude_primitives = nullptr; + r.exclude_count = 0; + r.volume_mesh_id = extParams.volume_mesh_ids[idx]; // Set volume mesh ID per ray + r.enabled = extParams.enabled; + + extParams.xdgRays[idx] = r; // Write to device ray buffer + } +} diff --git a/tests/test_direct_ray_buffer_access_shared.h b/tests/test_direct_ray_buffer_access_shared.h new file mode 100644 index 00000000..e8d5ceb3 --- /dev/null +++ b/tests/test_direct_ray_buffer_access_shared.h @@ -0,0 +1,13 @@ +#include "gprt.h" + +#include "../include/xdg/gprt/ray.h" + +struct ExternalRayParams { + xdg::dblRay* xdgRays; + double3* origins; + double3* directions; + uint num_rays; + uint total_threads; + int32_t* volume_mesh_ids; + uint enabled; +}; diff --git a/tests/test_files b/tests/test_files index a3caf0af..ca579198 160000 --- a/tests/test_files +++ b/tests/test_files @@ -1 +1 @@ -Subproject commit a3caf0af3f128944c4d6eac93b481df6e4efd97c +Subproject commit ca57919851224047ef86fab177a0bfe9fa920127 diff --git a/tests/test_point_in_volume.cpp b/tests/test_point_in_volume.cpp index ae34e823..acb62a41 100644 --- a/tests/test_point_in_volume.cpp +++ b/tests/test_point_in_volume.cpp @@ -24,7 +24,6 @@ TEMPLATE_TEST_CASE("Point-in-volume on MeshMock", "[piv][mock]", check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time auto rti = create_raytracer(rt_backend); REQUIRE(rti); - rti->init(); // Keep MeshMock usage consistent across backends auto mm = std::make_shared(false); @@ -78,3 +77,71 @@ TEMPLATE_TEST_CASE("Point-in-volume on MeshMock", "[piv][mock]", REQUIRE(result == false); } } + +TEMPLATE_TEST_CASE("Batch API Point-in-volume on MeshMock", "[piv][mock][batch]", + Embree_Raytracer, + GPRT_Raytracer) +{ + constexpr auto rt_backend = TestType::value; + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); + if (rt_backend == RTLibrary::EMBREE) { + SKIP("Skipping PIV batch for Embree: batch API not implemented yet"); + } + + auto rti = create_raytracer(rt_backend); + REQUIRE(rti); + + auto mm = std::make_shared(false); + mm->init(); + REQUIRE(mm->mesh_library() == MeshLibrary::MOCK); + + auto [volume_tree, element_tree] = rti->register_volume(mm, mm->volumes()[0]); + REQUIRE(volume_tree != ID_NONE); + REQUIRE(element_tree == ID_NONE); + + rti->init(); + + std::vector points; + std::vector directions; + std::vector has_dir; + size_t N; + + SECTION("N=0 no-op") { + rti->point_in_volume(volume_tree, + nullptr, /*points*/ + 0, /*num_points*/ + nullptr /*results*/); + SUCCEED("N=0 completed without error"); + } + + SECTION("N=1") { + N = 1; + make_points(N, points, directions); + + auto scalar_result = static_cast(rti->point_in_volume(volume_tree, points[0], &directions[0])); + + std::vector batch_result(N, 0xFF); + rti->point_in_volume(volume_tree, points.data(), N, batch_result.data(), directions.data()); + REQUIRE(batch_result[0] == scalar_result); + } + + SECTION("N=64") { + N = 64; + make_points(N, points, directions); + + // Store results of scalar point_in_volume calls to verify batch against scalar + std::vector scalar_results(N, 0); + for (size_t i = 0; i < N; ++i) { + scalar_results[i] = static_cast(rti->point_in_volume(volume_tree, points[i], &directions[i])); + } + + std::vector batch_results(N, 0xFF); + rti->point_in_volume(volume_tree, points.data(), N, batch_results.data(), directions.data()); + for (size_t i = 0; i < points.size(); ++i) { + REQUIRE(batch_results[i] == scalar_results[i]); + } + } + } +} diff --git a/tests/test_ray_fire.cpp b/tests/test_ray_fire.cpp index be816c36..552a4ed6 100644 --- a/tests/test_ray_fire.cpp +++ b/tests/test_ray_fire.cpp @@ -4,13 +4,14 @@ #include #include - // xdg includes #include "xdg/constants.h" #include "xdg/mesh_manager_interface.h" #include "mesh_mock.h" #include "util.h" +#include + using namespace xdg; using namespace xdg::test; @@ -109,4 +110,81 @@ TEMPLATE_TEST_CASE("Ray Fire on MeshMock (per-backend sections)", "[rayfire][moc intersection = rti->ray_fire(volume_tree, origin, direction, INFTY, HitOrientation::EXITING, &exclude_primitives); REQUIRE(intersection.second == ID_NONE); } -} \ No newline at end of file +} + +TEMPLATE_TEST_CASE("Batch API Ray Fire on MeshMock", "[rayfire][mock][batch]", + Embree_Raytracer, + GPRT_Raytracer) +{ + constexpr auto rt_backend = TestType::value; + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); + if (rt_backend == RTLibrary::EMBREE) { + SKIP("Skipping batch query mechanics test for Embree: batch API not implemented."); + } + + auto rti = create_raytracer(rt_backend); + REQUIRE(rti); + + auto mm = std::make_shared(false); + mm->init(); + REQUIRE(mm->mesh_library() == MeshLibrary::MOCK); + + auto [volume_tree, element_tree] = rti->register_volume(mm, mm->volumes()[0]); + REQUIRE(volume_tree != ID_NONE); + REQUIRE(element_tree == ID_NONE); + + rti->init(); + + std::vector origins; + std::vector directions; + size_t N; + + // ---- N = 0 ---- + SECTION("N=0 no-op") { + rti->ray_fire(volume_tree, nullptr, nullptr, 0, nullptr, nullptr, + INFTY, HitOrientation::EXITING, nullptr); + SUCCEED("N=0 completed without error"); + } + + // ---- N = 1 ---- + SECTION("N=1 equals scalar") { + N = 1; + make_rays(N, origins, directions); + + auto [dist_scalar, id_scalar] = rti->ray_fire(volume_tree, origins[0], directions[0], INFTY, HitOrientation::EXITING); + + double dist_batch = -1.0; + MeshID id_batch = ID_NONE; + rti->ray_fire(volume_tree, origins.data(), directions.data(), 1, + &dist_batch, &id_batch, INFTY, HitOrientation::EXITING, nullptr); + + REQUIRE(id_batch == id_scalar); + REQUIRE_THAT(dist_batch, Catch::Matchers::WithinAbs(dist_scalar, 1e-6)); + } + + // ---- N = 64 ---- + SECTION("N=64 matches scalar for all") { + N = 64; + make_rays(N, origins, directions); + + std::vector dist_scalar(64, INFTY); + std::vector id_scalar(64, ID_NONE); + for (size_t i = 0; i < 64; ++i) { + auto [d, id] = rti->ray_fire(volume_tree, origins[i], directions[i], INFTY, HitOrientation::EXITING); + dist_scalar[i] = d; id_scalar[i] = id; + } + + std::vector dist_batch(64, -1.0); + std::vector id_batch(64, ID_NONE); + rti->ray_fire(volume_tree, origins.data(), directions.data(), origins.size(), + dist_batch.data(), id_batch.data(), INFTY, HitOrientation::EXITING, nullptr); + + for (size_t i = 0; i < 64; ++i) { + REQUIRE(id_batch[i] == id_scalar[i]); + REQUIRE_THAT(dist_batch[i], Catch::Matchers::WithinAbs(dist_scalar[i], 1e-6)); + } + } + } +} diff --git a/tests/util.h b/tests/util.h index 579ff652..58d81d8f 100644 --- a/tests/util.h +++ b/tests/util.h @@ -96,3 +96,32 @@ create_raytracer(xdg::RTLibrary rt) { return nullptr; } + +inline void make_rays(size_t N, std::vector& origins, std::vector& directions) +{ + origins.clear(); + directions.clear(); + origins.reserve(N); + directions.reserve(N); + for (size_t i = 0; i < N; ++i) { + int axis = static_cast(i % 3); + double s = (i % 2) ? 1.0 : -1.0; + origins.push_back({0.0, 0.0, 0.0}); + if (axis == 0) directions.push_back({s, 0.0, 0.0}); + else if (axis == 1) directions.push_back({0.0, s, 0.0}); + else directions.push_back({0.0, 0.0, s}); + } +} + +inline void make_points(size_t N, std::vector& points, std::vector& directions) +{ + points.resize(N); + directions.resize(N); + for (size_t i = 0; i < N; ++i) { + // evens inside (origin), odds just outside +X + points[i] = (i % 2 == 0) ? xdg::Position{0.0, 0.0, 0.0} : xdg::Position{5.1, 0.0, 0.0}; + // alternate ±X directions + directions[i] = (i % 2 == 0) ? xdg::Direction{1.0, 0.0, 0.0} + : xdg::Direction{-1.0, 0.0, 0.0}; + } +} diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index f23f871e..00d7ecd3 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -1,8 +1,10 @@ set(TOOL_NAMES particle_sim ray_fire +batch_ray_fire find_volume point_in_volume +batch_point_in_volume overlap_check walk_elements tally_segments @@ -30,4 +32,6 @@ foreach(tool ${TOOL_NAMES}) target_compile_definitions(${tool_exec} PUBLIC XDG_OPENMP) endif() install(TARGETS ${tool_exec} DESTINATION ${CMAKE_INSTALL_BINDIR}/tools) -endforeach() \ No newline at end of file +endforeach() + +add_subdirectory(ray_benchmark) diff --git a/tools/batch_point_in_volume.cpp b/tools/batch_point_in_volume.cpp new file mode 100644 index 00000000..56ab46c4 --- /dev/null +++ b/tools/batch_point_in_volume.cpp @@ -0,0 +1,194 @@ +#include +#include +#include +#include +#include +#include + +#include "xdg/error.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/vec3da.h" +#include "xdg/xdg.h" + +#include "argparse/argparse.hpp" + +using namespace xdg; + +int main(int argc, char** argv) { + + argparse::ArgumentParser args("XDG Batch Point In Volume Tool", "1.0", argparse::default_arguments::help); + + args.add_argument("filename") + .help("Path to the input file"); + + args.add_argument("volume") + .help("Volume ID to query").scan<'i', int>(); + + args.add_argument("-l", "--list") + .default_value(false) + .implicit_value(true) + .help("List all volumes in the file and exit"); + + args.add_argument("-o", "-p", "--origin", "--position") + .default_value(std::vector{0.0, 0.0, 0.0}) + .help("Ray origin/position. Repeat to supply multiple origins.") + .scan<'g', double>().nargs(3).append(); + + args.add_argument("-d", "--direction") + .default_value(std::vector{0.0, 0.0, 1.0}) + .help("Ray direction. Repeat to supply multiple directions.") + .scan<'g', double>().nargs(3).append(); + + + args.add_argument("-m", "--mesh-library") + .help("Mesh library to use. One of (MOAB, LIBMESH)") + .default_value("MOAB"); + + args.add_argument("-r", "--rt-library") + .help("Ray tracing library to use. One of (EMBREE, GPRT)") + .default_value("GPRT"); + + // High-level rules in the description + args.add_description( + "Directions are completely optional for this tool but the number provided will effect how the program runs: \n\n" + " Only points (mask all, device default dir used)\n" + " --origin 0 0 0 --origin 5.1 0 0 --origin 0 0 0\n\n" + " One direction (broadcast to all)\n" + " --origin 0 0 0 --origin 5.1 0 0 --direction 1 0 0\n\n" + " Several directions. Match to points and mask remainder\n" + " --origin 0 0 0 --origin 5.1 0 0 --origin 4.999999 0 0 \\\n" + " --direction 1 0 0 --direction -1 0 0\n" + ); + + try { + args.parse_args(argc, argv); + } + catch (const std::runtime_error& err) { + std::cout << err.what() << std::endl; + std::cout << args; + exit(0); + } + + std::string mesh_str = args.get("--mesh-library"); + std::string rt_str = args.get("--rt-library"); + + MeshLibrary mesh_lib; + if (mesh_str == "MOAB") mesh_lib = MeshLibrary::MOAB; + else if (mesh_str == "LIBMESH") fatal_error("LibMesh is not currently supported with GPRT"); + else fatal_error("Invalid mesh library '{}' specified", mesh_str); + + RTLibrary rt_lib; + if (rt_str == "EMBREE") rt_lib = RTLibrary::EMBREE; + else if (rt_str == "GPRT") rt_lib = RTLibrary::GPRT; + else fatal_error("Invalid ray tracing library '{}' specified", rt_str); + + // create a mesh manager + std::shared_ptr xdg = XDG::create(mesh_lib, rt_lib); + const auto& mm = xdg->mesh_manager(); + mm->load_file(args.get("filename")); + mm->init(); + mm->parse_metadata(); + + auto rti = xdg->ray_tracing_interface(); + + if (args.get("--list")) { + std::cout << "Volumes: " << std::endl; + for (auto volume : mm->volumes()) { + std::cout << volume << std::endl; + } + exit(0); + } + + MeshID volume = args.get("volume"); + xdg->prepare_volume_for_raytracing(volume); + + // Gather our inputs and determine which mode of operation the tool will be working in + auto flat_origins = args.get>("--origin"); + auto flat_directions = args.get>("--direction"); + + if (flat_origins.empty()) { + fatal_error("You must supply at least one --origin x y z"); + } + if (flat_origins.size() % 3 != 0) { + fatal_error("Origins must be supplied in groups of 3 numbers."); + } + + // group every 3 into Position / Direction + std::vector> args_origins; + for (size_t i = 0; i < flat_origins.size(); i += 3) { + args_origins.push_back({flat_origins[i], flat_origins[i+1], flat_origins[i+2]}); + } + + std::vector> args_directions; + for (size_t i = 0; i < flat_directions.size(); i += 3) { + args_directions.push_back({flat_directions[i], flat_directions[i+1], flat_directions[i+2]}); + } + + // helper lambdas to convert std::vector to xdg::Position and xdg::Direction types + auto vec_to_pos = [](const std::vector& v) { return Position{v[0], v[1], v[2]}; }; + auto vec_to_dir = [](const std::vector& v) { + Direction dir{v[0], v[1], v[2]}; + dir.normalize(); + return dir; + }; + + const size_t N = args_origins.size(); + size_t num_dirs = args_directions.size(); + + std::vector origins; + origins.reserve(N); + for (const auto& o : args_origins) origins.push_back(vec_to_pos(o)); + + std::vector directions; + std::vector has_dir; // mask to indicate which rays have directions + const Direction* directions_ptr = nullptr; + const uint8_t* has_dir_ptr = nullptr; + + if (num_dirs == 0) { + // No directions let batch API set default direction per point + directions_ptr = nullptr; + has_dir_ptr = nullptr; + } else if (num_dirs == 1) { + // Broadcast one direction to all points (no mask needed) + directions.assign(N, vec_to_dir(args_directions[0])); + directions_ptr = directions.data(); + has_dir_ptr = nullptr; + } else if (num_dirs < N) { + // First k get explicit directions; rest fall back to default via mask + const size_t k = num_dirs; + directions.resize(N); + has_dir.assign(N, 0); + for (size_t i = 0; i < k; ++i) { + directions[i] = vec_to_dir(args_directions[i]); + has_dir[i] = 1; + } + directions_ptr = directions.data(); + has_dir_ptr = has_dir.data(); + } else { + // ≥ N directions → use first N pairwise (no mask needed) + directions.reserve(N); + for (size_t i = 0; i < N; ++i) directions.push_back(vec_to_dir(args_directions[i])); + directions_ptr = directions.data(); + has_dir_ptr = nullptr; + } + + std::vector results(N, 0xFF); + + xdg->point_in_volume(volume, + origins.data(), + N, + results.data(), + directions.data()); + + std::cout << std::endl << "Printing Batch point in volume results..." << std::endl; + + std::cout << "\nPrinting Batch point-in-volume results...\n"; + for (size_t i = 0; i < N; ++i) { + const auto& p = origins[i]; + std::cout << "Point (" << p.x << ", " << p.y << ", " << p.z << ") " + << (results[i] ? "is in " : "is NOT in ") + << "Volume " << volume << "\n"; + } + + return 0; +} diff --git a/tools/batch_ray_fire.cpp b/tools/batch_ray_fire.cpp new file mode 100644 index 00000000..a208ed74 --- /dev/null +++ b/tools/batch_ray_fire.cpp @@ -0,0 +1,226 @@ +#include +#include +#include +#include +#include + +#include "xdg/error.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/vec3da.h" +#include "xdg/xdg.h" + +#include "argparse/argparse.hpp" + +enum class BatchMode { + ORIGIN_BROADCAST, // 1 origin, many directions + DIRECTION_BROADCAST, // many origins, 1 direction + PAIRWISE // equal numbers of origins and directions +}; + +inline const char* to_string(BatchMode mode) { + switch (mode) { + case BatchMode::ORIGIN_BROADCAST: return "ORIGIN_BROADCAST"; + case BatchMode::DIRECTION_BROADCAST: return "DIRECTION_BROADCAST"; + case BatchMode::PAIRWISE: return "PAIRWISE"; + default: return "UNKNOWN"; + } +} + +inline BatchMode deduce_batch_mode(size_t num_origins, size_t num_directions) { + if (num_origins == 0 || num_directions == 0) { + throw std::runtime_error("At least one origin and one direction must be provided."); + } + + if (num_origins == 1 && num_directions > 1) { + return BatchMode::ORIGIN_BROADCAST; + } + else if (num_directions == 1 && num_origins > 1) { + return BatchMode::DIRECTION_BROADCAST; + } + else if (num_origins == num_directions) { + return BatchMode::PAIRWISE; + } + else { + throw std::runtime_error( + "Invalid combination: number of origins (" + std::to_string(num_origins) + + ") does not match number of directions (" + std::to_string(num_directions) + + ") for broadcast or pairwise mode." + ); + } +} + +using namespace xdg; + +int main(int argc, char** argv) { + + argparse::ArgumentParser args("XDG Batch Ray Fire Tool", "1.0", argparse::default_arguments::help); + + args.add_argument("filename") + .help("Path to the input file"); + + args.add_argument("volume") + .help("Volume ID to query").scan<'i', int>(); + + args.add_argument("-l", "--list") + .default_value(false) + .implicit_value(true) + .help("List all volumes in the file and exit"); + + args.add_argument("-o", "-p", "--origin", "--position") + .default_value(std::vector{0.0, 0.0, 0.0}) + .help("Ray origin/position. Repeat to supply multiple origins.") + .scan<'g', double>().nargs(3).append(); + + args.add_argument("-d", "--direction") + .default_value(std::vector{0.0, 0.0, 1.0}) + .help("Ray direction. Repeat to supply multiple directions.") + .scan<'g', double>().nargs(3).append(); + + + args.add_argument("-m", "--mesh-library") + .help("Mesh library to use. One of (MOAB, LIBMESH)") + .default_value("MOAB"); + + args.add_argument("-r", "--rt-library") + .help("Ray tracing library to use. One of (EMBREE, GPRT)") + .default_value("GPRT"); + + // High-level rules in the description + args.add_description( + "This tool supports two modes of operation for batch ray firing: 'Broadcast' and 'Pairwise'\n\n" + "To use 'Broadcast' mode, provide one origin and many directions, or one direction and many origins:\n" + " --origin x y z --direction u1 v1 w1 --direction u2 v2 w2 ...\n" + " --direction u v w --origin x1 y1 z1 --origin x2 y2 z2 ...\n\n" + "To use 'Pairwise' mode, each origin is paired with a corresponding direction in order:\n" + " --origin x1 y1 z1 --direction u1 v1 w1 --origin x2 y2 z2 --direction u2 v2 w2 ...\n" + ); + + try { + args.parse_args(argc, argv); + } + catch (const std::runtime_error& err) { + std::cout << err.what() << std::endl; + std::cout << args; + exit(0); + } + + std::string mesh_str = args.get("--mesh-library"); + std::string rt_str = args.get("--rt-library"); + + MeshLibrary mesh_lib; + if (mesh_str == "MOAB") mesh_lib = MeshLibrary::MOAB; + else if (mesh_str == "LIBMESH") fatal_error("LibMesh is not currently supported with GPRT"); + else fatal_error("Invalid mesh library '{}' specified", mesh_str); + + RTLibrary rt_lib; + if (rt_str == "EMBREE") rt_lib = RTLibrary::EMBREE; + else if (rt_str == "GPRT") rt_lib = RTLibrary::GPRT; + else fatal_error("Invalid ray tracing library '{}' specified", rt_str); + + // create a mesh manager + std::shared_ptr xdg = XDG::create(mesh_lib, rt_lib); + const auto& mm = xdg->mesh_manager(); + mm->load_file(args.get("filename")); + mm->init(); + mm->parse_metadata(); + + auto rti = xdg->ray_tracing_interface(); + + if (args.get("--list")) { + std::cout << "Volumes: " << std::endl; + for (auto volume : mm->volumes()) { + std::cout << volume << std::endl; + } + exit(0); + } + + MeshID volume = args.get("volume"); + xdg->prepare_volume_for_raytracing(volume); + + // Gather our inputs and determine which mode of operation the tool will be working in + auto flat_origins = args.get>("--origin"); + auto flat_directions = args.get>("--direction"); + + if (flat_origins.empty()) { + fatal_error("You must supply at least one --origin x y z"); + } + if (flat_origins.size() % 3 != 0) { + fatal_error("Origins must be supplied in groups of 3 numbers."); + } + if (flat_directions.size() % 3 != 0) { + fatal_error("Directions must be supplied in groups of 3 numbers."); + } + + // group every 3 into Position / Direction + std::vector> args_origins; + for (size_t i = 0; i < flat_origins.size(); i += 3) { + args_origins.push_back({flat_origins[i], flat_origins[i+1], flat_origins[i+2]}); + } + + std::vector> args_directions; + for (size_t i = 0; i < flat_directions.size(); i += 3) { + args_directions.push_back({flat_directions[i], flat_directions[i+1], flat_directions[i+2]}); + } + + // helper lambdas to convert std::vector to xdg::Position and xdg::Direction types + auto vec_to_pos = [](const std::vector& v) { return Position{v[0], v[1], v[2]}; }; + auto vec_to_dir = [](const std::vector& v) { + Direction dir{v[0], v[1], v[2]}; + dir.normalize(); + return dir; + }; + + size_t num_orig = args_origins.size(); + size_t num_dirs = args_directions.size(); + + auto mode = deduce_batch_mode(num_orig, num_dirs); + std::cout << "Running XDG Batch Ray Fire in " << to_string(mode) << " mode" << std::endl; + std::vector origins; + std::vector directions; + + switch (mode) + { + case BatchMode::ORIGIN_BROADCAST: + origins.assign(num_dirs, vec_to_pos(args_origins[0])); + directions.reserve(num_dirs); + for (const auto& dir : args_directions) directions.push_back(vec_to_dir(dir)); + break; + case BatchMode::DIRECTION_BROADCAST: + directions.assign(num_orig, vec_to_dir(args_directions[0])); + origins.reserve(num_orig); + for (const auto& origin : args_origins) origins.push_back(vec_to_pos(origin)); + break; + case BatchMode::PAIRWISE: + origins.reserve(num_orig); + directions.reserve(num_dirs); + for (size_t i = 0; i < num_orig; ++i) + { + origins.push_back(vec_to_pos(args_origins[i])); + directions.push_back(vec_to_dir(args_directions[i])); + } + break; + + default: + fatal_error("You must provide either a single origin and many directions. " + "A single direction and many origins. Or an equal number of origins and directions."); + } + + size_t num_rays = origins.size(); // get number of rays to fire from now aligned arrays + + std::vector hitDistances(num_rays); + std::vector surfacesHit(num_rays); + + xdg->ray_fire(volume, origins.data(), directions.data(), num_rays, hitDistances.data(), surfacesHit.data()); + + std::cout << std::endl << "Printing Batch Ray results..." << std::endl; + + for (size_t i = 0; i < num_rays; ++i) { + std::cout << "Ray[" << i << "] " + << "Origin=(" << origins[i].x << ", " << origins[i].y << ", " << origins[i].z << ") " + << "Dir=(" << directions[i].x << ", " << directions[i].y << ", " << directions[i].z << ") " + << "Distance=" << std::setprecision(17) << hitDistances[i] << " " + << "| Surface=" << surfacesHit[i] << "\n"; + } + + return 0; +} diff --git a/tools/ray_benchmark/CMakeLists.txt b/tools/ray_benchmark/CMakeLists.txt new file mode 100644 index 00000000..9566a620 --- /dev/null +++ b/tools/ray_benchmark/CMakeLists.txt @@ -0,0 +1,33 @@ +#=============================================================================== +# ray-benchmark (special case - requires linking directly to GPRT) +#=============================================================================== +if (XDG_ENABLE_GPRT) + # Embed and compile the device code + embed_devicecode( + OUTPUT_TARGET + ray_benchmark_deviceCode + HEADERS + ${CMAKE_CURRENT_SOURCE_DIR}/ray_benchmark_shared.h + SOURCES + ${CMAKE_CURRENT_SOURCE_DIR}/ray_benchmark_deviceCode.slang + ) + + # Create the ray-benchmark executable + add_executable(ray-benchmark ray_benchmark.cpp) + target_link_libraries(ray-benchmark xdg argparse ray_benchmark_deviceCode) + # Keep the runtime output alongside other tools for single- and multi-config generators. + get_filename_component(TOOLS_BIN_DIR "${CMAKE_CURRENT_BINARY_DIR}" DIRECTORY) + set_target_properties(ray-benchmark PROPERTIES + RUNTIME_OUTPUT_DIRECTORY "${TOOLS_BIN_DIR}" + ) + foreach(config DEBUG RELEASE RELWITHDEBINFO MINSIZEREL) + set_target_properties(ray-benchmark PROPERTIES + RUNTIME_OUTPUT_DIRECTORY_${config} "${TOOLS_BIN_DIR}" + ) + endforeach() + if (OpenMP_CXX_FOUND) + target_link_libraries(ray-benchmark OpenMP::OpenMP_CXX) + target_compile_definitions(ray-benchmark PUBLIC XDG_OPENMP) + endif() + install(TARGETS ray-benchmark DESTINATION ${CMAKE_INSTALL_BINDIR}/tools) +endif() diff --git a/tools/ray_benchmark/ray_benchmark.cpp b/tools/ray_benchmark/ray_benchmark.cpp new file mode 100644 index 00000000..afc975c1 --- /dev/null +++ b/tools/ray_benchmark/ray_benchmark.cpp @@ -0,0 +1,231 @@ +#include +#include +#include +#include +#include + +#include "xdg/error.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/moab/mesh_manager.h" +#include "xdg/vec3da.h" +#include "xdg/xdg.h" +#include "xdg/ray_tracers.h" +#include "xdg/timer.h" + +#include "argparse/argparse.hpp" + +#include "ray_benchmark.h" + +#include + +using namespace xdg; + +int main(int argc, char** argv) { + + argparse::ArgumentParser args("XDG Ray Tracing throughput benchmarking tool", "1.0", argparse::default_arguments::help); + + args.add_argument("filename") + .help("Path to the input file"); + + args.add_argument("volume") + .help("Volume ID to query") + .scan<'i', int>(); + + args.add_argument("-n", "--num-rays") + .default_value(10'000'000) + .help("Number of rays to be cast for the benchmark (default - 10 million)") + .scan<'u', uint32_t>(); + + args.add_argument("-s", "--seed") + .default_value(12345) + .help("Seed for random number generator (default - 12345)") + .scan<'u', uint32_t>(); + + args.add_argument("-o", "-p", "--origin", "--position") + .default_value(std::vector{0.0, 0.0, 0.0}) + .help("Ray origin/position (default - {0.0, 0.0, 0.0} )") + .scan<'g', double>().nargs(3); + + args.add_argument("-m", "--mesh-library") + .help("Mesh library to use. One of (MOAB, LIBMESH)") + .default_value("MOAB"); + + args.add_argument("-rt", "--rt-library") + .help("Ray tracing library to use. One of (EMBREE, GPRT)") + .default_value("EMBREE"); + + args.add_argument("-l", "--list") + .default_value(false) + .implicit_value(true) + .help("List all volumes in the file and exit"); + + args.add_argument("-sr", "--source-radius") + .default_value(0.0) + .help("Radius of a scattered source blob around the origin (0.0 = point source)") + .scan<'g', double>(); + + args.add_description( + "This tool supports can be used to benchmark XDG ray tracing throughput on a given mesh against" + "a given volume \n." + "A single origin/seed point is provided and ray directions are randomly generated in 360 degrees from that position" + ); + + try { + args.parse_args(argc, argv); + } + catch (const std::runtime_error& err) { + std::cout << err.what() << std::endl; + std::cout << args; + return 1; + } + + std::string mesh_str = args.get("--mesh-library"); + std::string rt_str = args.get("--rt-library"); + + RTLibrary rt_lib; + if (rt_str == "EMBREE") + rt_lib = RTLibrary::EMBREE; + else if (rt_str == "GPRT") + rt_lib = RTLibrary::GPRT; + else + fatal_error("Invalid ray tracing library '{}' specified", rt_str); + + MeshLibrary mesh_lib; + if (mesh_str == "MOAB") { + mesh_lib = MeshLibrary::MOAB; + } else if (mesh_str == "LIBMESH") { + mesh_lib = MeshLibrary::LIBMESH; + if (rt_lib == RTLibrary::GPRT) + fatal_error("LibMesh is not currently supported with GPRT"); + } else { + fatal_error("Invalid mesh library '{}' specified", mesh_str); + } + + // Full wall-clock timer (post-argparse) + Timer wall_timer; + wall_timer.start(); + + // Separate timers for setup, generation, and tracing + Timer setup_timer; + Timer gen_timer; + Timer trace_timer; + + // -------------------------- + // XDG setup timing + // -------------------------- + setup_timer.start(); + + std::shared_ptr xdg = XDG::create(mesh_lib, rt_lib); + const auto& mm = xdg->mesh_manager(); + mm->load_file(args.get("filename")); + mm->init(); + + MeshID volume = args.get("volume"); + xdg->prepare_raytracer(); + xdg->prepare_volume_for_raytracing(volume); + auto rti = xdg->ray_tracing_interface(); + + setup_timer.stop(); + + std::size_t N = args.get("--num-rays"); + uint32_t seed = args.get("--seed"); + Position origin = args.get>("--origin"); + double source_radius = args.get("--source-radius"); + + std::cout << "Volume ID: " << volume << " with: " + << mm->num_volume_faces(volume) << " faces" << std::endl; + + + if (rt_lib == RTLibrary::EMBREE) { + int num_threads = omp_get_max_threads(); + rt_str += " (" + std::to_string(num_threads) + " CPU threads)"; + } + std::cout << "Starting ray fire benchmark with " << N << " rays" + << " using " << rt_str << ": \n" << std::endl; + + std::cout << "XDG initalisation Time = " << setup_timer.elapsed() << "s" << std::endl; + + std::shared_ptr gprt_rt; + if (rt_lib == RTLibrary::GPRT) { + // ---- Random ray generation on device via callback method ---- + gen_timer.start(); + + gprt_rt = std::dynamic_pointer_cast(xdg->ray_tracing_interface()); + auto generateRaysCallback = + tools::benchmark::make_generate_rays_callback(gprt_rt->context(), origin, source_radius, seed, volume); + + // Let XDG internally allocate buffers and invoke the callback to populate them + xdg->populate_rays_external(N, generateRaysCallback); + + gen_timer.stop(); + std::cout << "Random ray generation (via external compute shader) Time = " + << gen_timer.elapsed() << "s" << std::endl; + + // ---- Ray tracing on device ---- + trace_timer.start(); + xdg->ray_fire_prepared(N); // ray_fire against pre-packed rays on device + trace_timer.stop(); + + } + else { // EMBREE / CPU backend + + // ---- Random ray generation on host ---- + gen_timer.start(); + std::vector directions(N); + std::vector origins(N); + + #pragma omp parallel for schedule(static) + for (uint32_t i = 0; i < N; ++i) { + uint32_t state = seed ^ i; + auto [pos,dir] = tools::benchmark::random_spherical_source(origin, state, source_radius); + origins[i] = pos; + directions[i] = dir; + } + gen_timer.stop(); + + std::cout << "Random ray generation Time = " + << gen_timer.elapsed() << "s" << std::endl; + + // ---- Ray tracing on host ---- + trace_timer.start(); + #pragma omp parallel for schedule(static) + for (std::size_t i = 0; i < N; ++i) { + auto result = xdg->ray_fire(volume, origins[i], directions[i]); + } + trace_timer.stop(); + } + + // -------------------------- + // Final reporting + // -------------------------- + double setup_time = setup_timer.elapsed(); + double gen_time = gen_timer.elapsed(); + double trace_time = trace_timer.elapsed(); + + double trace_only_rps = (trace_time > 0.0) + ? static_cast(N) / trace_time + : 0.0; + + double end_to_end_time = gen_time + trace_time; + double end_to_end_rps = (end_to_end_time > 0.0) + ? static_cast(N) / end_to_end_time + : 0.0; + + wall_timer.stop(); + double wall_time = wall_timer.elapsed(); + + std::cout << "Generation + tracing time = " << end_to_end_time + << "s" << std::endl; + std::cout << "End-to-end throughput = " << end_to_end_rps + << " rays/s" << std::endl; + std::cout << "Full wall-clock time = " << wall_time + << "s (post-argparse)" << std::endl; + + std::cout << "----------------------------------------" << std::endl; + std::cout << "Ray Tracing Time (trace-only) = " << trace_time + << "s for " << N << " rays" << std::endl; + std::cout << "Trace-only throughput = " << trace_only_rps + << " rays/s" << std::endl; + std::cout << "---------------------------------------- \n" << std::endl; + return 0; +} diff --git a/tools/ray_benchmark/ray_benchmark.h b/tools/ray_benchmark/ray_benchmark.h new file mode 100644 index 00000000..52b29d62 --- /dev/null +++ b/tools/ray_benchmark/ray_benchmark.h @@ -0,0 +1,97 @@ +#ifndef _XDG_RAY_BENCHMARK_H +#define _XDG_RAY_BENCHMARK_H + +#include +#include +#include + +#include "gprt/gprt.h" +#include "xdg/gprt/ray.h" +#include "xdg/gprt/ray_tracer.h" +#include "xdg/vec3da.h" +#include "xdg/xdg.h" + +#include "ray_benchmark_shared.h" + +extern GPRTProgram ray_benchmark_deviceCode; + +namespace xdg::tools::benchmark { + +inline double rand01(uint32_t &state) +{ + state = state * 1664525u + 1013904223u; + return double(state) * (1.0 / 4294967296.0); +} + +inline Direction random_unit_dir_lcg(uint32_t &state) +{ + double x1, x2, 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); + + double t = 2.0 * std::sqrt(1.0 - s); + return { x1 * t, x2 * t, 1.0 - 2.0 * s }; +} + +// Generates a random point cloud with radius (--source-radius) +inline std::pair random_spherical_source(const Position& origin, + std::uint32_t state, + double source_radius) +{ + // Always generate random direction + Direction dir = random_unit_dir_lcg(state); + Position pos = origin; + if (source_radius > 0.0) { + // random origins (spherical source) + double r = source_radius * std::cbrt(rand01(state)); // uniform in ball + pos += dir * r; + } + return {pos, dir}; +} + +// - User creates their own GPU compute API method to populate rays and passes that to XDG +// - In this miniapp we are using GPRT as a demonstration +// - This callback runs inside populate_rays_external and receives XDG's device buffers +inline RayPopulationCallback make_generate_rays_callback(GPRTContext gprt_context, + Position origin, + double source_radius, + uint32_t seed, + MeshID volume) +{ + return [gprt_context, origin, source_radius, seed, volume](const DeviceRayHitBuffers& buffer, size_t numRays) { + GPRTContext context = gprt_context; + GPRTModule module = gprtModuleCreate(context, ray_benchmark_deviceCode); + auto genRandomRays = gprtComputeCreate( + context, module, "generate_random_rays"); + + constexpr int threadsPerGroup = 64; + const int neededGroups = static_cast((numRays + threadsPerGroup - 1) / threadsPerGroup); + const int groups = std::min(neededGroups, WORKGROUP_LIMIT); + + GenerateRandomRayParams randomRayParams = {}; + randomRayParams.rays = static_cast(buffer.rayDevPtr); // Cast opaque pointer to typed dblRay* + randomRayParams.numRays = static_cast(numRays); + randomRayParams.source_radius = source_radius; + randomRayParams.origin = { origin.x, origin.y, origin.z }; + randomRayParams.seed = seed; + randomRayParams.total_threads = static_cast(groups * threadsPerGroup); + randomRayParams.volume_mesh_id = volume; + randomRayParams.enabled = 1u; + + gprtComputeLaunch(genRandomRays, + { static_cast(groups), 1, 1 }, + { static_cast(threadsPerGroup), 1, 1 }, + randomRayParams); + gprtComputeSynchronize(context); + + gprtComputeDestroy(genRandomRays); + gprtModuleDestroy(module); + }; +} + +} // namespace xdg::tools::benchmark + +#endif // _XDG_RAY_BENCHMARK_H diff --git a/tools/ray_benchmark/ray_benchmark_deviceCode.slang b/tools/ray_benchmark/ray_benchmark_deviceCode.slang new file mode 100644 index 00000000..e9a6aefb --- /dev/null +++ b/tools/ray_benchmark/ray_benchmark_deviceCode.slang @@ -0,0 +1,57 @@ +#include "ray_benchmark_shared.h" + +/* +For this simple benchmark case we are mocking what a downstream application would do in terms of populating +ray buffers. The idea is that the downstream application generates rays (origins + directions). +*/ +[shader("compute")] +[numthreads(64, 1, 1)] +void generate_random_rays(uint3 DispatchThreadID: SV_DispatchThreadID, + uniform GenerateRandomRayParams params) +{ + uint globalThreadID = DispatchThreadID.x; + uint stride = params.total_threads; + uint nRays = params.numRays; + + for (uint idx = globalThreadID; idx < nRays; idx += stride) + { + uint state = params.seed ^ idx; + + double3 dir = random_unit_dir_lcg(state); + + double3 pos = params.origin; + if (params.source_radius > 0.0) { + double u = float(rand01(state)); + float r = float(params.source_radius) * pow(float(u), 1.0f / 3.0f); // cbrt(u) + pos += dir * double(r); + } + + params.rays[idx].origin = pos; + params.rays[idx].direction = dir; + params.rays[idx].exclude_primitives = nullptr; + params.rays[idx].exclude_count = 0; + params.rays[idx].enabled = params.enabled; + params.rays[idx].volume_mesh_id = params.volume_mesh_id; + } +} + +// Simple LCG random number generator +double rand01(inout uint state) +{ + state = state * 1664525u + 1013904223u; + return double(state) * double(1.0 / 4294967296.0); +} + +// return random unit dir +double3 random_unit_dir_lcg(inout uint state) +{ + double x1, x2, 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); + + double t = 2.0 * sqrt(1.0 - s); + return double3(x1 * t, x2 * t, 1.0 - 2.0 * s); +} \ No newline at end of file diff --git a/tools/ray_benchmark/ray_benchmark_driver.py b/tools/ray_benchmark/ray_benchmark_driver.py new file mode 100644 index 00000000..1839418f --- /dev/null +++ b/tools/ray_benchmark/ray_benchmark_driver.py @@ -0,0 +1,297 @@ +#!/usr/bin/env python3 +import subprocess +import statistics +import sys +import csv +import os + +# --- CONFIG --- + +BENCHMARK = "./tools/ray-benchmark" +MESH_PATH = "../dagmc_xdg_test.h5m" +VOLUME_ID = "2" +NUM_RAYS = "80000000" +ORIGIN = ["-o", "180", "250", "-27"] # x y z as strings + +# --- PARSING HELPERS --- + +def parse_float_before_s(s: str) -> float: + """ + Given a string like 'XDG initalisation Time = 1.25017s', + pull out 1.25017 as float. + """ + try: + after_eq = s.split('=', 1)[1] + number_str = after_eq.split('s', 1)[0].strip() + return float(number_str) + except Exception as e: + raise ValueError(f"Failed to parse float from line: {s!r}") from e + +def parse_throughput_line(s: str) -> float: + """ + Given a string like 'Trace-only throughput = 2.64065e+09 rays/s', + pull out 2.64065e+09 as float. + """ + try: + after_eq = s.split('=', 1)[1] + number_str = after_eq.split('rays', 1)[0].strip() + return float(number_str) + except Exception as e: + raise ValueError(f"Failed to parse throughput from line: {s!r}") from e + +def parse_benchmark_output(output: str): + """ + Parse the benchmark stdout text and return a dict of metrics. + Expected keys: + - xdg_init + - gen + - gen_trace + - end_to_end + - wall_clock + - trace_only + - trace_only_throughput + """ + metrics = {} + + for line in output.splitlines(): + line = line.strip() + + if line.startswith("XDG initalisation Time"): + metrics["xdg_init"] = parse_float_before_s(line) + + elif line.startswith("Random ray generation"): + metrics["gen"] = parse_float_before_s(line) + + elif line.startswith("Generation + tracing time"): + metrics["gen_trace"] = parse_float_before_s(line) + + elif line.startswith("End-to-end throughput"): + metrics["end_to_end"] = parse_throughput_line(line) + + elif line.startswith("Full wall-clock time"): + metrics["wall_clock"] = parse_float_before_s(line) + + elif line.startswith("Ray Tracing Time (trace-only)"): + metrics["trace_only"] = parse_float_before_s(line) + + elif line.startswith("Trace-only throughput"): + metrics["trace_only_throughput"] = parse_throughput_line(line) + + required = [ + "xdg_init", "gen", "gen_trace", "end_to_end", + "wall_clock", "trace_only", "trace_only_throughput" + ] + missing = [k for k in required if k not in metrics] + if missing: + raise RuntimeError(f"Missing metrics in output: {missing}") + + return metrics + +# --- MAIN DRIVER --- + +def main(): + # Ask for backend + backend_in = input("Choose backend (embree/gprt): ").strip().lower() + if backend_in not in ("embree", "gprt"): + print("Invalid backend, please choose 'embree' or 'gprt'.") + sys.exit(1) + + base_backend = backend_in.upper() # what we pass to -r: EMBREE or GPRT + + # If GPRT, ask for which variant + if backend_in == "gprt": + mode_in = input( + "GPRT mode: [1] GPRT (FP64), [2] GPRT (FP32) + RT cores [1]: " + ).strip() + if mode_in == "2": + variant = "fp32_rt" + label = "GPRT (FP32) + RT cores" + else: + variant = "fp64" + label = "GPRT (FP64)" + else: + # Embree is effectively FP64 for your purposes + variant = "fp64" + label = "Embree" + + runs_str = input("How many runs? ").strip() + try: + num_runs = int(runs_str) + if num_runs <= 0: + raise ValueError + except ValueError: + print("Number of runs must be a positive integer.") + sys.exit(1) + + # Ask for CSV filename + csv_filename = input("CSV output file [benchmarks.csv]: ").strip() + if not csv_filename: + csv_filename = "benchmarks.csv" + + mesh_name = os.path.basename(MESH_PATH) + + all_metrics = { + "xdg_init": [], + "gen": [], + "gen_trace": [], + "end_to_end": [], + "wall_clock": [], + "trace_only": [], + "trace_only_throughput": [], + } + + # CSV header: machine-friendly backend/variant, plus pretty label + header = [ + "backend", # EMBREE / GPRT + "variant", # fp64 / fp32_rt + "label", # Embree / GPRT (FP64) / GPRT (FP32) + RT cores + "mesh_name", + "volume_id", + "num_rays", + "run_index", + "xdg_init", + "gen", + "gen_trace", + "end_to_end", + "wall_clock", + "trace_only", + "trace_only_throughput", + ] + + # Decide whether to append or overwrite + file_exists = os.path.exists(csv_filename) + write_header = False + file_mode = "w" + append_mode = False + + if file_exists: + choice = input( + f"File '{csv_filename}' already exists. " + "[o]verwrite, [a]ppend, or e[x]it? [a]: " + ).strip().lower() + + if choice in ("x", "q"): + print("Aborting, no benchmarks run.") + sys.exit(0) + elif choice in ("", "a"): + file_mode = "a" + write_header = False # assume header already there + append_mode = True + elif choice == "o": + file_mode = "w" + write_header = True + append_mode = False + else: + print("Unrecognized choice, aborting.") + sys.exit(1) + else: + # new file: write header + file_mode = "w" + write_header = True + append_mode = False + + csv_file = open(csv_filename, file_mode, newline="") + + # If appending, add a separation comment line so it's obvious this is a new batch + if append_mode: + csv_file.write( + f"\n# --- New benchmark batch: " + f"label={label}, backend={base_backend}, variant={variant}, " + f"mesh={mesh_name}, volume={VOLUME_ID}, " + f"rays={NUM_RAYS}, runs={num_runs} ---\n" + ) + + writer = csv.writer(csv_file) + + if write_header: + writer.writerow(header) + + try: + for i in range(1, num_runs + 1): + print(f"\n=== Run {i}/{num_runs} ({label}) ===") + + cmd = [ + BENCHMARK, + MESH_PATH, + VOLUME_ID, + "-r", base_backend, # EMBREE or GPRT + "-n", NUM_RAYS, + *ORIGIN, + ] + + print("Running:", " ".join(cmd)) + + try: + result = subprocess.run( + cmd, + check=True, + text=True, + capture_output=True, + ) + except subprocess.CalledProcessError as e: + print("Benchmark command failed!") + print("STDOUT:\n", e.stdout) + print("STDERR:\n", e.stderr) + sys.exit(1) + + try: + metrics = parse_benchmark_output(result.stdout) + except Exception as e: + print("Failed to parse benchmark output:", e) + print("Raw output:\n", result.stdout) + sys.exit(1) + + # store for averages + for k in all_metrics.keys(): + all_metrics[k].append(metrics[k]) + + # write CSV row + writer.writerow([ + base_backend, # backend + variant, # variant + label, # label + mesh_name, + VOLUME_ID, + NUM_RAYS, + i, # run_index + metrics["xdg_init"], + metrics["gen"], + metrics["gen_trace"], + metrics["end_to_end"], + metrics["wall_clock"], + metrics["trace_only"], + metrics["trace_only_throughput"], + ]) + + # per-run summary + print(f"XDG init : {metrics['xdg_init']:.6f} s") + print(f"Generation : {metrics['gen']:.6f} s") + print(f"Gen + trace : {metrics['gen_trace']:.6f} s") + print(f"End-to-end : {metrics['end_to_end']:.3e} rays/s") + print(f"Wall-clock : {metrics['wall_clock']:.6f} s") + print(f"Trace-only : {metrics['trace_only']:.6f} s") + print(f"Trace-only thrpt : {metrics['trace_only_throughput']:.3e} rays/s") + + finally: + csv_file.close() + + # Averages + print( + "\n=== Averages over", + num_runs, + f"runs (label: {label}) ===" + ) + + def avg(key): return statistics.mean(all_metrics[key]) + + print(f"Avg XDG init : {avg('xdg_init'):.6f} s") + print(f"Avg Generation : {avg('gen'):.6f} s") + print(f"Avg Gen + trace : {avg('gen_trace'):.6f} s") + print(f"Avg End-to-end : {avg('end_to_end'):.3e} rays/s") + print(f"Avg Wall-clock : {avg('wall_clock'):.6f} s") + print(f"Avg Trace-only : {avg('trace_only'):.6f} s") + print(f"Avg Trace-only thrpt : {avg('trace_only_throughput'):.3e} rays/s") + print(f"\nResults written to: {csv_filename}") + +if __name__ == "__main__": + main() diff --git a/tools/ray_benchmark/ray_benchmark_shared.h b/tools/ray_benchmark/ray_benchmark_shared.h new file mode 100644 index 00000000..a3fffcf0 --- /dev/null +++ b/tools/ray_benchmark/ray_benchmark_shared.h @@ -0,0 +1,14 @@ +#include "gprt.h" + +#include "../../include/xdg/gprt/ray.h" + +struct GenerateRandomRayParams { + xdg::dblRay* rays; + uint numRays; + double3 origin; + uint seed; + uint total_threads; + double source_radius; + int volume_mesh_id; + uint enabled; +}; diff --git a/vendor/GPRT b/vendor/GPRT index f1e95e41..405d9ee9 160000 --- a/vendor/GPRT +++ b/vendor/GPRT @@ -1 +1 @@ -Subproject commit f1e95e4188cde591547d6b4a33a70bf2afaeec59 +Subproject commit 405d9ee9f5ee8e1a0455f776f9e2c3adffb64160