Skip to content

Gprt vol tracing - #205

Open
Waqar-ukaea wants to merge 17 commits into
xdg-org:mainfrom
Waqar-ukaea:GPRT-vol-tracing
Open

Gprt vol tracing#205
Waqar-ukaea wants to merge 17 commits into
xdg-org:mainfrom
Waqar-ukaea:GPRT-vol-tracing

Conversation

@Waqar-ukaea

@Waqar-ukaea Waqar-ukaea commented Feb 27, 2026

Copy link
Copy Markdown
Collaborator

This PR adds all of the necessary setup to allow GPRT to perform ray tracing against volumetric (tet) meshes.

Changes:

  • Added all required GPRT setup/shader plumbing to handle construction of BVHs on tet meshes as well as shaders to trace against those BVHs
  • Added GPRTRayTracer::find_element() and updated test_find_element to test across both Embree and GPRT backends. Note - important to stress that this is a single ray dispatch version and an batch based dispatch will come in a later PR. Single dispatch versions do remain very useful for testing parity with Embree however.
  • For device side AABB population I have implemented a slang generic (equivalent to C++ templating) which will be nice for future PRs that add support for elements with differing numbers of nodes.
  • Restructured file structure for shaders to split triangle and tet shaders
  • Made plucker_tet_containment_test cross compatible between slang and C++. In doing so I have opted to make use of Cramer's rule rather than matrix inversion since I was having trouble with matrices between host and device. GPRT's double3x3 required implementing my own matrix inversion function and had a different row-column ordering compared to linalg's implementation.
  • Some minor QOL changes across GPRTRayTracer

It looks like a much larger PR than it really is since this PR also makes some structural changes to the file layout across shaders. Namely it splits shaders into three .slang files rt_common.slangh, triangle_rt_shaders.slang and tetrahedron_shaders.slang. So around 🟩 +180 additions and 🟥 −270 deletions are solely from moving the existing shaders for surface ray tracing into the new files.

Extra notes:

  • This is not an implementation of the full volumetric tracking algorithm we have in place that allows us to walk elements
  • But rather just all of the RT side BVH setup and shader setup to do element point containment checks against the BVH in find_element()
  • Since the PR for batch based array queries isn't in place this is just an initial attempt with each raygen shader launching a single ray. 

I guess in theory should we ever need to actually perform ray tracing operations against tets it also allows for that (would just need to define a new custom intersection shader for double precision) but as far as i know there is no use case for this required by OpenMC?

@pshriwise pshriwise left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few thoughts and comments here.

One larger question: Where is the box expansion occurring in the GPRTRayTracer? I see the expansion using FLT_EPSILON in populate_aabb, but for larger volumes this expansion may not be enough to guarantee we don't run into false negatives when traversing the tree.

Comment thread include/xdg/geometry/dp_math.h Outdated
const dp::vec3 v2,
const dp::vec3 v3) {

// TODO - I've decided to use Cramer's rule here instead of matrix inversion to make it easier to implement a cross-compilable version of this function

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In general this looks fine to me. In my branch for quads/hexes I've moved to a Plucker coordinate check against face and I might recommend we change to that broadly as it naturally extends to linear elements with an arbitrary number of faces, but I'm curious, what was incompatible in the previous version?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's been a while, but I think the main issue was related to row/column convention differences between the matrix API used by GPRT/Slang and the linalg API used on the C++ side. The old method (still in main) was doing this matrix multiplication to solve for our barycentric coords:

    // Solve T * [λ1, λ2, λ3] = rhs
    double3 lambda123 = mul(inverse(T),{rhs.x, rhs.y, rhs.z});

I initially tried to implement the same logic on the Slang side, but I got different barycentric coordinates on CPU vs GPU. My understanding is that the two math APIs were interpreting the matrix row/column ordering differently which was leading to differences when I would call mul(...).

I switched to Cramer's rule to avoid the matrix inversion and multiplication path entirely. That keeps the shared Slang/C++ implementation to just dot and cross products, which is easier to keep consistent across both C++ and slang.

I might recommend we change to that broadly as it naturally extends to linear elements with an arbitrary number of faces, but I'm curious, what was incompatible in the previous version?

Agreed, seems reasonable to move towards a more arbitrary method.

Comment thread include/xdg/gprt/ray_tracer.h Outdated

// Ray Generation parameters
uint32_t numRayTypes_ = 1; // <! Number of ray types. Allows multiple shaders to be set to the same geometery
uint32_t numRayTypes_ = 2; // <! Number of ray types. 0=surface, 1=volume

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's create a new enum for this property if there isn't one already.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I see there are some constants in shared_structs.h. Let's apply RT_NUM_RAY_TYPES here.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I see there are some constants in shared_structs.h. Let's apply RT_NUM_RAY_TYPES here.

Yeah, not sure what the best place for these would be since they're constants but will only ever be used inside of GPRTRayTracer and the slang shaders so it didn't feel appropriate to put them in constants.h. But open to suggestions on that front.

Comment thread include/xdg/gprt/rt_common.slangh Outdated
Comment thread src/gprt/ray_tracer.cpp Outdated
Comment thread src/gprt/tetrahedron_rt_shaders.slang Outdated
Comment thread src/gprt/tetrahedron_rt_shaders.slang Outdated
Comment thread src/gprt/triangle_rt_shaders.slang Outdated
uint rayID = DispatchRaysIndex().x;

// There is some logic for handling next volumes inside the h5m-reader which I could make use of too
// TODO : Should the dblHit struct return the next volume ID for the ray back to the host

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Responding to this comment randomly: We could do but I don't see a reason to necessarily. As long as we know what volume we're in and what surface was hit it's a really cheap operation to find the next volume and in the case of particle transport we may not use it if a collision will occur before the particle reaches that surface.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense. I'll remove this comment then.

Comment thread src/gprt/triangle_rt_shaders.slang Outdated
Comment thread src/gprt/triangle_rt_shaders.slang Outdated

// ------------------------------------------------ CUSTOM INTERSECTION SHADERS ------------------------------------------------

/* 1D ray generation intersection with a double precision triangle using the Plucker intersection algorithm*/

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for the stray: The term "1D" threw me here for a sec.

Suggested change
/* 1D ray generation intersection with a double precision triangle using the Plucker intersection algorithm*/
/* Single ray generation intersection with a double precision triangle using the Plucker intersection algorithm*/

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed it is a little misleading. I've changed the comment to this in the latest commit:

// Custom FP64 Plucker ray-triangle intersection algorithm for each ray-triangle pair
[shader("intersection")]
void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record)
{

@Waqar-ukaea

Copy link
Copy Markdown
Collaborator Author

@pshriwise bb61bf6 addresses some of the comments you left. Most of the changes are self-explanatory, but i'll highlight one here:

I added a new shared_constants.h header so the shader code can use xdg::ID_NONE instead of hard-coded -1 values as you suggested. constants.h currently pulls in C++-only dependencies so it cant be included directly in the Slang shader path so I made a new header shared_constants.h which doesn't contain any C++ specific stuff.

I also moved the enums from shared_enums.h into that header and removed the old shared_enums.h.

A couple of small syntax changes were needed to keep the header Slang-compatible: constexpr was replaced with static const, and brace initialization like MeshID ID_NONE {-1} was replaced with MeshID ID_NONE = -1.

@Waqar-ukaea

Copy link
Copy Markdown
Collaborator Author

One larger question: Where is the box expansion occurring in the GPRTRayTracer? I see the expansion using FLT_EPSILON in populate_aabb, but for larger volumes this expansion may not be enough to guarantee we don't run into false negatives when traversing the tree.

I see your point, this is missing from GPRT atm and we only do the FLT_EPSILON expansion you saw. I'll think about how best to handle this and get a separate PR up for it once I figure it out

Comment thread include/xdg/shared_constants.h Outdated
#ifndef XDG_SHARED_CONSTANTS_H
#define XDG_SHARED_CONSTANTS_H

namespace xdg {

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As discussed live this week, probably best to avoid changing the infrastructure too drastically to accommodate slang. Sorry for the mislead here.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No problem! I'll revert back to what I had before

@pshriwise pshriwise left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A few more thoughts as I look over this today.

// Geometry Type and Instances
std::vector<gprt::Instance> globalBlasInstances_; //<! List of every BLAS instance stored in this ray tracer
GPRTGeomTypeOf<DPTriangleGeomData> trianglesGeomType_; //<! Geometry type for triangles
GPRTGeomTypeOf<DPTetrahedronGeomData> tetrahedraGeomType_; //<! Geometry type for tetrahedra (not currently supported)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not currently supported?


// Ray Generation parameters
uint32_t numRayTypes_ = 1; // <! Number of ray types. Allows multiple shaders to be set to the same geometery
uint32_t numRayTypes_ = RT_NUM_RAY_TYPES; // <! Number of ray types. 0=surface, 1=volume

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's remove the values in the comments here lest they change in the code itself and we forget to update the values here.

Comment on lines +184 to +186
const double barycentric_min = -dp::PLUCKER_ZERO_TOL;
const double barycentric_max = 1.0 + dp::PLUCKER_ZERO_TOL;

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems like these could be compile-time constants, no?

int volume_tree; // TreeID of the volume being queried
SurfaceAccelerationStructure volume_accel; // The volume accel
SurfaceAccelerationStructure volume_accel_surf; // The volume accel for surface acceleration structure
SolidAccelerationStructure volume_accel_solid; // The volume accel for solid acceleration structure

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The use of solid here feels a little off. What about "ElementAccelerationStrucutre"?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do agree but both SurfaceAccelerationStructure and SolidAccelerationStructure are types defined in GPRT here - https://github.com/gprt-org/GPRT/blob/master/gprt/gprt_builtins.slang

So I would need to put in a PR to make the name change upstream to make that change

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or I suppose I could typedef them here locally in XDG to have that extra consistency with the naming scheme we have.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, I missed that. Thanks for clarifying. typedef's sound reasonable to me. I think it would help with readability and allow us to use consistent terminology. Let's try it on.

@Waqar-ukaea Waqar-ukaea added the GPRT GPRT specific changes which have no effect on the existing Embree interface label Apr 14, 2026
@Waqar-ukaea

Waqar-ukaea commented Apr 21, 2026

Copy link
Copy Markdown
Collaborator Author

One larger question: Where is the box expansion occurring in the GPRTRayTracer? I see the expansion using FLT_EPSILON in populate_aabb, but for larger volumes this expansion may not be enough to guarantee we don't run into false negatives when traversing the tree.

This has now been addressed in #222 and I recommend that #222 PR is merged before this PR is so that I can rebase those changes properly into this branch first.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

GPRT GPRT specific changes which have no effect on the existing Embree interface

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants