Geometry debug function

[viktormai]

This PR adds a 3D geometry debugging utility to Model for identifying overlap and undefined regions within an OpenMC geometry by sampling a user-defined bounding box.

Main changes

Added Model.geometry_debug():

• Samples a 3D region using a stack of 2D slices obtained through openmc.lib.slice_data.
• Detects both overlap regions (using the overlap sentinel) and internal undefined regions, returning a summary containing the total number of overlap and undefined samples and example coordinates for each type of geometry issue.
• Optionally prints a human-readable summary.

Added _classify_undefined_regions() helper:

• Separates undefined pixels into boundary-connected ("outside") undefined regions and internal undefined regions enclosed within the geometry.
• Uses a breadth-first search starting from the image boundary to identify boundary-connected undefined pixels.
• Issues a warning when undefined regions exist but none are connected to the sampled slice boundary, indicating a potentially enclosed undefined region or insufficient sampling resolution.

This pair of functions provides a simple way to locate common geometry construction problems—particularly overlaps and enclosed undefined regions—by returning representative coordinates that users can inspect directly, rather than relying solely on transport failures or visual inspection. These sample points can also be leveraged by AI agents to easily fix geometry issues in the model. Reusable undefined-region classification logic that can be leveraged by future plotting and geometry diagnostics for plotting undefined regions is also established.

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)

[paulromano]

Thanks for the initial shot at this! Here are my initial thoughts:

[paulromano]

You should also use the new slice_data_overlap_info to get information about the overlaps. Rather than storing an explicit list of every coordinate, I would store a mapping of (univ, cell1, cell2) to a bounding box that covers all points found to overlap for that combination of cells.

[paulromano]

I asked GPT-5.6 Sol whether there was a better way to implement this functionality and here is what it told me:

SciPy’s ndimage module is a much better fit ... and SciPy is already an OpenMC dependency.

The simplest replacement is scipy.ndimage.binary_fill_holes:

from scipy import ndimage

undefined = cell_ids == _NOT_FOUND

# Holes in the defined-pixel mask are internal undefined regions.
internal = ndimage.binary_fill_holes(~undefined) & undefined
outside = undefined & ~internal

Its default connectivity is the same four-neighbor connectivity used by the current BFS. This eliminates the deque, boundary initialization, and Python-level pixel traversal.

In a local benchmark with a representative undefined mask:

Grid	Current BFS	binary_fill_holes
100 × 100	3.8 ms	0.15 ms
500 × 500	95 ms	2.8 ms
1000 × 1000	393 ms	11 ms

That’s roughly 25–35× faster, although slice_data() may still dominate the overall geometry_debug() runtime.

Another semantically direct option is binary_propagation, seeded with undefined boundary pixels. It performed about the same but requires more setup. ndimage.label was fastest in my benchmark—around 3.5 ms at 1000²—but requires additional logic to identify labels touching the boundary and allocates an integer label array.

My recommendation for the PR would therefore be binary_fill_holes: it is considerably shorter, preserves the existing connectivity semantics, adds no dependency, and moves the expensive traversal into compiled SciPy code.

[paulromano]

The return type hint doesn't match all possibilities (namely, possibility of returning tuple of None)

[paulromano]

Elsewhere in the API, when the user provides a single number of samples/points, it gets distributed over the three dimensions rather than being used for each dimension.

[viktormai]

Thanks for checking it out @paulromano. I have added some updates for the bounding box idea for both overlapping and undefined regions, including a warning if undefined regions are under-resolved. Let me know what you think.

[paulromano]

Thanks for the updates @viktormai! Here is another round of comments:

[paulromano]

Under what circumstances would cell_ids be None? Should we not raise an exception in this case?

[paulromano]

There should be a "Returns" section in the docstring here describing what is returned

[paulromano]

Personally, I wouldn't attach this to the Model class as a @staticmethod but just have it be a separate helper function at the top-level of this or another module. Given that we plan on using it in the plotter, it also shouldn't be hidden (leading underscore) but rather exposed as a public API function.

[paulromano]

The return type hint doesn't match all possibilities (namely, possibility of returning tuple of None)

[paulromano]

This gives nx + ny + nz ≈ n_samples, which is not what we want. We need nx × ny × nz ≈ n_samples, which can be achieved with something like:

width = np.asarray(upper_right) - np.asarray(lower_left)
scale = np.cbrt(n_samples / np.prod(width))
nx, ny, nz = np.maximum(1, np.rint(scale * width).astype(int))

[paulromano]

If you store the bounding boxes using the BoundingBox class, you can take advantage of the | operator, which does exactly this logic (take the maximum extent of two boxes).

[paulromano]

I don't understand this warning -- what's wrong with having only internal undefined pixels?