Retry locating particles that start on a surface coincident with another

src/particle.cpp

@@ -223,7 +223,20 @@ void Particle::event_calculate_xs()
   // initiate a search for the current cell. This generally happens at the
   // beginning of the history and again for any secondary particles
   if (lowest_coord().cell() == C_NONE) {
-    if (!exhaustive_find_cell(*this)) {
+    bool found = exhaustive_find_cell(*this);
+    if (!found && surface() != SURFACE_NONE) {
+      // A particle that starts on a surface, e.g., one split by weight windows
+      // at a surface crossing, is only known to be on that surface. If another
+      // surface coincides with it, roundoff in the particle's position may
+      // place the particle on the wrong side of the other surface so that no
+      // cell contains it. As when crossing a surface, move the particle
+      // forward a tiny bit and search again.
+      surface() = SURFACE_NONE;
+      n_coord() = 1;
+      r() += TINY_BIT * u();
+      found = exhaustive_find_cell(*this);
+    }
+    if (!found) {
       mark_as_lost(
         "Could not find the cell containing particle " + std::to_string(id()));
       return;

tests/unit_tests/test_lost_particles.py

@@ -47,3 +47,50 @@ def test_max_write_lost_particles(model: openmc.Model, run_in_tmpdir):
     n_procs = int(config['mpi_np']) if config['mpi'] else 1
     assert len(lost_particle_files) == model.settings.max_write_lost_particles * n_procs
 
+
+
+def test_split_on_coincident_surfaces(run_in_tmpdir):
+    """Particles split by weight windows on a surface that coincides with
+    another surface must not be lost (GitHub issue #4140)"""
+    s1 = openmc.XPlane(-57452.33336021505)
+    s2 = openmc.XPlane(43403.32187479845)
+    s3 = openmc.XPlane(62125.04976135607)
+    s4 = openmc.XPlane(62125.04976135607)
+    s5 = openmc.XPlane(100000.0, boundary_type='vacuum')
+    inner = openmc.Universe(cells=[
+        openmc.Cell(region=-s1 | -s2 | -s3),
+        openmc.Cell(region=+s3 & -s5),
+    ])
+    root = openmc.Cell(fill=inner, region=-s1 | -s2 | -s3 | +s4)
+
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([root])
+    model.settings.run_mode = 'fixed source'
+    model.settings.particles = 1
+    model.settings.batches = 1
+    model.settings.max_lost_particles = 100
+    model.settings.rel_max_lost_particles = 0.999999
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point((-85626.45049221347, 0.0, 0.0)),
+        angle=openmc.stats.Monodirectional(
+            (0.9182609440159194, 0.39597580569397484, 0.0)),
+        energy=openmc.stats.delta_function(1.0e6),
+    )
+
+    # The particle is split into ten when crossing s3 into the second mesh bin
+    mesh = openmc.RegularMesh()
+    mesh.lower_left = (-1.0e5, -1.0e5, -1.0)
+    mesh.upper_right = (2.2e5, 1.0e5, 1.0)
+    mesh.dimension = (2, 1, 1)
+    model.settings.weight_windows = openmc.WeightWindows(
+        mesh, lower_ww_bounds=[0.5, 0.01], upper_ww_bounds=[5.0, 0.1],
+        energy_bounds=[0.0, 2.0e7])
+    model.settings.weight_window_checkpoints = {
+        'collision': False, 'surface': True}
+
+    sp_path = model.run()
+
+    # All weight should leak out of the vacuum boundary
+    with openmc.StatePoint(sp_path) as sp:
+        leakage = sp.global_tallies[sp.global_tallies['name'] == b'leakage']
+        assert leakage['mean'][0] == pytest.approx(1.0)