Fix virtual surface crossing

[GuySten]

Description

Currently, when using complex cells openmc sometimes decide that a particle cross surface that is not on the boundary of the complex cell.

In the illustration below a particle in cell 2 cross a surface and still stays in cell 2 (a virtual crossing).

This PR fix that behavior by adding another check in complex cells that after the surface crossing the particle should be outside the cell.

Fixes #3924
Fixes #3186

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[HunterBelanger]

I haven't looked at this in detail at all, but I think this would also fix the problem I described in #3186 if I am not mistaken.

[GuySten]

I think you are right. Are you interested in reviewing this PR?

[HunterBelanger]

Except for the really bad statistics, it looks like the flux suppression in the reflector is now gone, so I think it might have fixed that problem too. I can try to take a look at this late this week.

[HunterBelanger]

Other than my one comment, this generally looks okay to me and seems to solve the problem I observed, and should also fix the other current bug which another user reported. We should see what someone else thinks though, before merging such an important change which has geometry changes.

[GuySten]

@HunterBelanger, could you approve this PR?
Afterwards I am planning to label this PR with a merging soon label and will merge it next week only if no one else will have any comments until then.

[paulromano]

Thanks a lot for proposing this @GuySten, and it's encouraging to see that it does address the behaviors that were documented in the linked issues! However, I am concerned about the performance impact as described below. Also, I tried this out on a few complex models that I have access to and on at least one of them (JET), it seems to result in an infinite loop. I'll have to dig into why this is but that will need to be resolved before this could move forward.

[paulromano]

I'm concerned about the performance impact of this. In the worst case scenario, this means that for a cell with N half-spaces bounding it, you could end up doing N extra cell containment checks. Obviously the average case will be less, but nevertheless it's a cause for concern. I'm going to look into alternatives for how to handle this.

[GuySten]

I've reduced the number of point containments.
I tried to have the virtual surface crossing computation inside the distance_complex function. I hope that makes sense to you.
I think now there should not be a performance degradation.
(You should keep in mind that when there are virtual surface crossing the distance complex function is called more times)

[GuySten]

@paulromano, I've tried to fix the infinite loop you mentioned in the JET problem by improving robustness of the algorithm in general and especially in the case of large geometry.
Could you check if the latest changes fix the problem?

[paulromano]

@GuySten Thanks a lot for your work on this. The original implementation correctly recognized that a geometrical surface intersection is not necessarily a boundary of a complex CSG region, but it did not fully rescan the ray after a virtual crossing. In particular, each search continued to use the original ray position, and each quadratic surface supplied only its nearest positive root. If the first root of a sphere/cylinder was virtual, the implementation could never discover the second root, even when that second intersection was the actual region boundary. For example, if you had a region like the following:

x = openmc.XPlane()
s = openmc.Sphere(r=1)
region = -x | -s

and a particle at (-2, 0, 0) traveling to the right, it would never see the correct boundary because it corresponds to the larger root of the quadratic formed from the intersection with the sphere.

To mitigate this, I've updated the implementation to use an "advance-and-rescan" algorithm for complex regions:

1. Determine whether the starting point is inside the region.
2. Find the nearest intersection with any surface in the expression from the current ray position.
3. Advance the ray to that intersection and accumulate the traveled distance.
4. Determine the side entered using the surface normal and ray direction.
5. Re-evaluate region membership at the intersection.
6. If membership changed, return the accumulated distance and surface. Otherwise, the crossing was virtual, so continue scanning from that intersection.

Because each rescan starts at the newly reached surface, quadratic surfaces can return their next root when appropriate.

[paulromano]

After further testing, I'm still seeing the infinite loop on the JET model but I've narrowed it down and should have a fix soon.

[paulromano]

Ok, it turned out the infinite loop in the JET geometry was due to two surfaces that were effectively redundant, a z-cylinder and a general quadric. When converting from the original MCNP geometry, they weren't de-duplicated by virtue of being different surface types. In the distance_complex method, it kept alternating between those two surfaces. To avoid this, I've added additional logic to handle the case where two surfaces are equidistant, which should also help for other pathological cases (e.g., a plane tangent to a cylinder).

I did some performance testing on a few different problems. For a geometry-heavy problem with only simple cells, the performance hit is only ~4%. For the JET model, I am seeing a 30% drop in performance. For the ITER E-lite model, I see a 10% drop in performance. While I wish this had less impact on performance for those cases, I think this is a necessary step for robustness and consistency. And, I have ideas for further optimization so hopefully we can get back some of that performance loss.