Use direction of travel when looking up weight window mesh bin

src/weight_windows.cpp

@@ -8,6 +8,7 @@
 
 #include "openmc/tensor.h"
 
+#include "openmc/constants.h"
 #include "openmc/error.h"
 #include "openmc/file_utils.h"
 #include "openmc/hdf5_interface.h"
@@ -289,9 +290,13 @@ std::pair<bool, WeightWindow> WeightWindows::get_weight_window(
   if (E < energy_bounds_.front() || E > energy_bounds_.back())
     return {false, {}};
 
-  // Get mesh index for particle's position
+  // Get mesh index for particle's position. The position is nudged along the
+  // direction of travel so that a particle sitting exactly on a mesh boundary
+  // (e.g., at a surface checkpoint where a geometry surface coincides with a
+  // mesh plane) is assigned to the element it is entering rather than always
+  // to the element on the lower-coordinate side.
   const auto& mesh = this->mesh();
-  int mesh_bin = mesh->get_bin(p.r());
+  int mesh_bin = mesh->get_bin(p.r() + TINY_BIT * p.u());
 
   // particle is outside the weight window mesh
   if (mesh_bin < 0)

tests/regression_tests/random_ray_auto_convert_bootstrap/results_true.dat

@@ -1,9 +1,9 @@
 tally 1:
-8.094106E+01
-6.974487E+02
-3.181906E-04
-1.455701E-08
-9.840422E+02
-1.110377E+05
-8.515372E+01
-8.078508E+02
+8.087496E+01
+7.339588E+02
+2.446860E-04
+8.405655E-09
+9.761237E+02
+1.018073E+05
+9.815520E+01
+1.065938E+03

tests/regression_tests/weightwindows/survival_biasing/local/results_true.dat

@@ -1 +1 @@
-0d7b17d4e364bda2be21feb2e5b2aae4be69fc8ed634c4f45062e8efc61b7bd24dcf448837692fd603afe3756501fe8821d134f2d6289179618f85178ddb8793
+7a9d895b430f0c7696fc36f73dc68179d484664a29bd168df3b28a2b0549552d6216e4e9b32de4449d8134a42a84883d7e7e7364f0b9e4cc75174e95ecaf6b1c

tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat

@@ -1 +1 @@
-2f1c8efbf6ab0b0c7319fadd7ee8f27b545c1ebab94ebe15cf1133b54227487542966b6d798970c8b37bc19e7041e687687d5285803241bb83f4c8dc40a4d276
+5ee91c9d9285829528fb2e53a1c5bf069318259e0d4e8768852309a921bc6dd77852de251bd6e91245f68bc33b0fdc2321cfc3b315432d822603976f74079396

tests/unit_tests/weightwindows/test_ww_surface_checkpoint.py

@@ -0,0 +1,91 @@
+import openmc
+import pytest
+
+
+@pytest.mark.parametrize("direction", [1.0, -1.0])
+def test_ww_surface_on_mesh_boundary(run_in_tmpdir, direction):
+    """Weight windows applied at a surface that coincides with a weight window
+    mesh boundary must use the mesh element the particle is entering,
+    regardless of the direction of travel.
+
+    At the surface checkpoint the particle sits exactly on the plane x = 0.
+    Structured mesh index lookups resolve a position exactly on an element
+    boundary to the element on the lower-coordinate side, so when the window
+    was looked up at the particle's position alone, a particle moving in +x
+    was given the window of the element it was leaving. It was then not split
+    and the +x case failed, while the -x case passed only because there the
+    lower-coordinate element happens to be the one being entered. The lookup
+    now nudges the position along the direction of travel.
+    """
+
+    # Void-like one-group material so the only weight window checks happen at
+    # surface crossings
+    groups = openmc.mgxs.EnergyGroups([0.0, 20.0e6])
+    xsdata = openmc.XSdata('void', groups)
+    xsdata.order = 0
+    xsdata.set_total([0.0])
+    xsdata.set_absorption([0.0])
+    xsdata.set_scatter_matrix([[[0.0]]])
+    mg_library = openmc.MGXSLibrary(groups)
+    mg_library.add_xsdata(xsdata)
+    mg_library.export_to_hdf5('mgxs.h5')
+
+    mat = openmc.Material()
+    mat.add_macroscopic('void')
+    materials = openmc.Materials([mat])
+    materials.cross_sections = 'mgxs.h5'
+
+    # Two cells separated by the plane x = 0, which is also a mesh boundary
+    x_min = openmc.XPlane(-1.0, boundary_type='vacuum')
+    x_mid = openmc.XPlane(0.0)
+    x_max = openmc.XPlane(1.0, boundary_type='vacuum')
+    yz = openmc.model.RectangularPrism(2.0, 2.0, axis='x',
+                                       boundary_type='vacuum')
+    left = openmc.Cell(fill=mat, region=+x_min & -x_mid & -yz)
+    right = openmc.Cell(fill=mat, region=+x_mid & -x_max & -yz)
+    geometry = openmc.Geometry([left, right])
+
+    # Particle is born in the window of its birth element and must be split
+    # five ways when it enters the element on the other side of x = 0
+    birth_cell, target_cell = (left, right) if direction > 0 else (right, left)
+    mesh = openmc.RectilinearMesh()
+    mesh.x_grid = [-1.0, 0.0, 1.0]
+    mesh.y_grid = [-1.0, 1.0]
+    mesh.z_grid = [-1.0, 1.0]
+    birth_bounds = (0.5, 1.5)
+    target_bounds = (0.1, 0.2)
+    if direction > 0:
+        lower = [birth_bounds[0], target_bounds[0]]
+        upper = [birth_bounds[1], target_bounds[1]]
+    else:
+        lower = [target_bounds[0], birth_bounds[0]]
+        upper = [target_bounds[1], birth_bounds[1]]
+    ww = openmc.WeightWindows(mesh, lower_ww_bounds=lower,
+                              upper_ww_bounds=upper)
+
+    settings = openmc.Settings()
+    settings.energy_mode = 'multi-group'
+    settings.run_mode = 'fixed source'
+    settings.particles = 10
+    settings.batches = 1
+    settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point((-0.5*direction, 0.0, 0.0)),
+        angle=openmc.stats.Monodirectional((direction, 0.0, 0.0)),
+    )
+    settings.weight_windows = ww
+    settings.weight_window_checkpoints = {'collision': False, 'surface': True}
+
+    tally = openmc.Tally()
+    tally.filters = [
+        openmc.CellFilter([target_cell]),
+        openmc.WeightFilter([0.0, 0.5, 2.0])
+    ]
+    tally.scores = ['flux']
+
+    model = openmc.Model(geometry, materials, settings, openmc.Tallies([tally]))
+    model.run(apply_tally_results=True)
+
+    # All flux in the target cell should be carried by split particles with
+    # weight 0.2; total flux is conserved (1 cm path length per source particle)
+    flux = tally.mean.squeeze()
+    assert flux == pytest.approx([1.0, 0.0])