Point in volume determination for location on triangle #986
Description
Describe the Bug
I'm working with a simple model of nested cubes in OpenMC to resolve a problem with sourcing secondary particles resulting from weight window splitting (openmc-dev/openmc#3601). The original problem was that the particle history was not being reset when the split particles were picked up from the secondary bank under certain conditions, but I'll leave that to the PR.
A new issue arose with this contrived model afterward that I believe is related to DAGMC. I found that the model ran smoothly with weight windows applied using DAGMC + double-down, but that I ran into a problem with DAGMC (no double-down). DAGMC using MOAB's GeomQueryTool fails with the following error:
Simulating batch 1
WARNING: Could not find the cell containing particle 251 at location
(-59.75000000000001, -20.311014690968396, -39.0653767943094)
(direction (-0.7137604953357174, -0.6851083376566464,
-0.14550780381635203))
It took a lot of digging to figure out what was going on here and, in the end, I think it's related to the box check we perform for the point_in_volume query.
Here's an image of the model containing three cubes with widths 20 cm, 115 cm and 120 cm. The last two surfaces are close together to severely expose the secondary particle sourcing issue in OpenMC.
So, here's our current setup:
Note that the particle is just outside the surface (exaggerated in the image). Really, this is in ambiguous/coincident territory so declaring the particle in either volume is valid IMO. Not declaring it in either is incorrect.
Before diving in further it's important to note that the nearest intersection of this ray with either volume 4 or 5 and no orientation culling (which is how our point in volume checks occur) is an intersection distance of -0 with the surface between volumes 4 and 5.
What happens in double-down
We loop over volumes calling point_in_volume. For volume 4 we obtain an exiting intersection at distance -0 with volume 4 and the search is over.
What happens in GeomQueryTool
The GeomQueryTool::point_in_volume call contains a check against the volume's bounding box.
When calling point_in_volume for volume 4, we declare the origin of the ray to be outside the bounding box of volume 4 and return False, the point is not in volume 4.
We then move on to a point_in_volume call for volume 5. It passes the initial bounding box check -- the point is well within the bounding box of volume 5. At this point we use the nearest hit along the ray direction (again, no orientation filtering here) to make a determination about the point's containment
The zero-distance intersection with the surface between 4 and 5 is returned but due to the orientation of the normals this is declared an entering hit and the particle is declared outside the volume.
Here's where things get fun. In the GeomQueryTool::add_intersection method this line is present:
which means that if we were to skip the initial bounding box check when querying volume 4, the GeomQueryTool would also accept the zero-distance hit with that surface and declare the particle to be inside volume 4. I've tried this by disabling the bounding box check and can confirm that it works.
Summary
So the crux of the matter to me is this: The bounding box check is an optimization to avoid a ray query. However, based on this example it can exclude operations for particles on the volume boundary that make DAGMC robust to these ambiguous situations.
FWIW, I think this is somewhat rare and I've stumbled upon a pathology related to situations where the bounding box fits the volume exactly. Normally our
Proposed fix
IMO if our ray kernel/triangle intersector/what have you is structured to accept ever-so-slightly negative distance intersections in these cases, we should not keep them from occurring due to this bounding box check.
The definition of tol in the MOAB code above defines the tolerance for which we accept hits within any distance of the ray origin in the MOAB code above comes from DAGMC's numerical precision attribute1. What I'd propose is that we build some tolerance into the point_in_volume bounding box check to account for this as well to ensure that we aren't excluding these kind of intersections with this optimization.
Link to slides with images in case these turned out too blurry here
Footnotes
-
this should be checking against the absolute value of the distance, not the distance itself if the comment above that line is described the correct/desired behavior ↩