Increase numerical robustness of cylindrical and spherical distance to boundary checks

[GuySten]

Description

This PR apply the same fixes that were done in #3792 to the spherical and cylindrical surfaces themselves.

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)

[gonuke]

It is not clear to my how this improves robustness. I would appreciate an explanation for what motivates these changes. I looked back at #3792 and didn't find any explanation there. There was a link to a discourse post that it was addressing, but there is a leap between the symptoms discussed in that post and the particular solution implemented there.

The changes here look mostly subtle except for the substantive change on the tolerance for determining that a point is coincident with a cylindrical/spherical surface. I'd like to understand the scaling of that change and why it's the necessary/best approach.

[GuySten]

@pshriwise, wrote a nice summary of this issue:
https://openmc.discourse.group/uploads/short-url/rBPUb4HPRzM7Cj8zgfS4yH26qYz.pdf

My motivation came from noticing that for small spheres, openmc falsely detect overlap for points close to the boundary (because the distance to boundary check is not robust enough). These changes fix that problem.

[lewisgross1296]

Would be good to have a test which fails with the old OpenMC and passes with these changes. For example, a test of the overlap condition in geometry.cpp which fails due to lack of robustness that now passes.

[gonuke]

@pshriwise, wrote a nice summary of this issue:

https://openmc.discourse.group/uploads/short-url/rBPUb4HPRzM7Cj8zgfS4yH26qYz.pdf

My motivation came from noticing that for small spheres, openmc falsely detect overlap for points close to the boundary (because the distance to boundary check is not robust enough). These changes fix that problem.

Thanks for sharing that summary by @pshriwise. Two follow-ups:

• could we not just compare |x|^2 - R^2 to a smaller threshold than FP_COINCIDENT and get most of the same benefit? This is different from the recommendation to compare things differently at different radii, but rather, compare small r^2 to a different threshold than small r
• in addition to changing this to a test of |x| - R instead of |x|^2 - R^2, this also changes the threshold from FP_COINCIDENT to FP_COINCIDENT * (1 + R). Why this change - that is not documented in the @pshriwise summary

[GuySten]

Regarding your second point, the difference is that the new way is using a relative tolerance of FP_COINCIDENT in addition to an absolute tolerance of FP_COINCIDENT.
You can find more info in the documentaion of python isclose function https://docs.python.org/3/library/math.html#math.isclose.

Regarding your first point, I am not sure. I think the way it is done now is more "natural" but your suggestion might also work.

[lewisgross1296]

I'm wondering if this test should use openmc.lib to actually confirm what cell it's in and what cell it thinks it's in, similar to the conditions in geometry.cpp to confirm that there is no overlap detected

openmc/src/geometry.cpp

Lines 44 to 52 in fd1bc26

	for (auto index_cell : univ.cells_) {
	Cell& c = *model::cells[index_cell];
	if (c.contains(p.coord(j).r(), p.coord(j).u(), p.surface())) {
	if (index_cell != p.coord(j).cell()) {
	if (error) {
	fatal_error(
	fmt::format("Overlapping cells detected: {}, {} on universe {}",
	c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_));
	}

I'm not sure if any tests actually test using openmc -g but that could be another way to do the same thing

[GuySten]

I am using geometry_debug=True which is equivalent to openmc -g.

[lewisgross1296]

Missed that, my bad. Looks good!

[lewisgross1296]

Some good news. My previous overlap check failed on batch 57, and I've just now made it to batch 58! The overlap checking is definitely going slower, which is to be expected due to the additional sqrt() operation as noted in Patrick's write up; however, as a debug feature I think this is fine.

I am going to let the overlap checker keep going just to be sure. The eigenvalues seem to be within statistical noise of each other, so I'm guessing this means some point containment decisions resulted in different answers with this new change, but not enough to drastically change results.