Fix get_index_in_direction for regular meshes

[JoffreyDorville]

Description

For regular meshes, the function get_index_in_direction is returning 0 instead of 1 when the tested coordinate is coincident with the lower boundary. This results in points being incorrectly seen outside the mesh when they lie on surfaces corresponding to the lower boundaries of the mesh.

For rectilinear meshes, the function get_index_in_direction works as expected.

I added a cpp unit test to illustrate the problem and to show that this PR fixes it.

This fix also corrected the results of the weightwindows/survival_biasing regression tests which have been updated. I cross-checked that using an equivalent rectilinear mesh instead of the regular mesh in tests/regression_tests/weightwindows/survival_biasing/test.py produces the same results as the updated ones with the fix.

Rectilinear mesh used:

mesh = openmc.RectilinearMesh()
mesh.x_grid = np.linspace(0.0, 160.0, 21)
mesh.y_grid = np.linspace(0.0, 160.0, 21)
mesh.z_grid = np.linspace(0.0, 160.0, 2)

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)

[GuySten]

Looks good to me.

[paulromano]

I'm not sure I agree with the change. When the position being tested is on the boundary between two mesh elements, it is ambiguous as to whether it should be considered as being in the mesh element to the left or to the right. To me it seems more consistent to always pick one direction.

We have other ways of handling this ambiguity elsewhere in the code such as 1) pushing the position forward by TINY_BIT to move it off the boundary, or 2) Looking at the dot product of u and the surface normal to decide which one to choose.

[GuySten]

I think using the particle direction is a good solution.
In any case, equivalent regular and rectilinear meshes should have the same behavior.
Currently they behave differently.

[JoffreyDorville]

Thanks for the feedback @paulromano!

I just added the direction check for regular and rectilinear meshes. The solution is even simpler than pushing the particle or using the surface normal.

For cylindrical and spherical meshes, the task is a bit more involved because of specific cases (i.e., non-convex shapes). I have some ideas on how to handle these cases but it can be on a different PR as it is not required for the current state of the temperature field PR.

I have not looked into unstructured meshes for the moment.

I changed the signature of several functions that needed the new direction parameter, which also had the drawback of changing the signature of functions that did not need it. Let me know what you prefer for the Mesh interface.

[paulromano]

I made some comments below regarding the current implementation. However, I think the bigger issue is that this seems like an overly complicated solution for a relatively minor consistency issue. Sorry to go back and forth here, but I would prefer either:

1. Reverting to the original fix, if (r == lower_left_[i]) return 1, or
2. Change the implementation of RectilinearMesh::get_index_in_direction to use std::lower_bound directly (std::lower_bound(grid_[i].begin(), grid_[i].end(), r) - grid_[i].begin()). The reason it currently behaves differently than RegularMesh is because lower_bound_index has a special case for the first value in an array.

I lean toward solution 2 because it would treat all mesh cells equivalently with respect to behavior on the boundary (each mesh cells owns its upper boundary but not its lower).

[paulromano]

This doesn't really match mentally what I was expecting. I would have thought that a particle that is coincident with any mesh boundary would account for the direction, not just the upper boundary. Additionally, the check for the upper boundary is done using floating-point equality, which seems like it may fail even when the particle is exactly on the boundary because of rounding/truncation.

[paulromano]

If r is below the first grid point, we would have that idx == 0, which means that grid_[i][idx - 1] would be out of bounds.

[paulromano]

Same consideration here: if r is above the last grid point, this will be out of bounds.