Retry locating particles that start on a surface coincident with another

src/particle.cpp

@@ -40,6 +40,36 @@
 
 namespace openmc {
 
+namespace {
+
+//! Find the cell containing a particle that has no cell yet
+//
+//! 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. Move the particle forward a tiny bit and search again. The surface is
+//! kept set, so that its sense is still taken from the side the particle is
+//! on rather than from its position: on a near-grazing direction the move may
+//! be too small to get past roundoff, and the particle must then be lost
+//! rather than placed in the cell on the other side.
+//
+//! \param p Particle to locate
+//! \param verbose Whether to print verbose search information
+//! \return Whether a cell containing the particle was found
+bool find_cell_from_surface(Particle& p, bool verbose = false)
+{
+  bool found = exhaustive_find_cell(p, verbose);
+  if (!found && p.surface() != SURFACE_NONE) {
+    p.n_coord() = 1;
+    p.r() += TINY_BIT * p.u();
+    found = exhaustive_find_cell(p, verbose);
+  }
+  return found;
+}
+
+} // namespace
+
 //==============================================================================
 // Particle implementation
 //==============================================================================
@@ -223,7 +253,7 @@ 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)) {
+    if (!find_cell_from_surface(*this)) {
       mark_as_lost(
         "Could not find the cell containing particle " + std::to_string(id()));
       return;
@@ -520,7 +550,7 @@ void Particle::event_revive_from_secondary(const SourceSite& site)
     // removed from the pulse-height of this cell.
     if (lowest_coord().cell() == C_NONE) {
       bool verbose = settings::verbosity >= 10 || trace();
-      if (!exhaustive_find_cell(*this, verbose)) {
+      if (!find_cell_from_surface(*this, verbose)) {
         mark_as_lost("Could not find the cell containing particle " +
                      std::to_string(id()));
         return;

tests/unit_tests/test_lost_particles.py

@@ -1,3 +1,4 @@
+import math
 from pathlib import Path
 
 import openmc
@@ -47,3 +48,85 @@ 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 _coincident_surfaces_model(particle,
+                               direction=(0.9182609440159194,
+                                          0.39597580569397484, 0.0)):
+    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(
+        particle=particle,
+        space=openmc.stats.Point((-85626.45049221347, 0.0, 0.0)),
+        angle=openmc.stats.Monodirectional(direction),
+        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.0e12, -1.0)
+    mesh.upper_right = (2.2e5, 1.0e12, 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], particle_type=particle)
+    model.settings.weight_window_checkpoints = {
+        'collision': False, 'surface': True}
+    return model
+
+
+@pytest.mark.parametrize('u_x', [0.9182609440159194, 1.0e-4])
+def test_split_on_coincident_surfaces(run_in_tmpdir, u_x):
+    """Particles split by weight windows on a surface that coincides with
+    another surface must not be lost (GitHub issue #4140), including on a
+    near-grazing crossing where the retry's small move is comparable to
+    roundoff in the position"""
+    model = _coincident_surfaces_model(
+        'neutron', direction=(u_x, math.sqrt(1.0 - u_x**2), 0.0))
+    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)
+
+
+def test_split_photons_on_coincident_surfaces_pulse_height(run_in_tmpdir):
+    """Split photons are located before their pulse-height contribution is
+    removed, and that search must not lose them either"""
+    model = _coincident_surfaces_model('photon')
+    model.settings.photon_transport = True
+    cells = [c for c in model.geometry.get_all_cells().values()
+             if not isinstance(c.fill, openmc.UniverseBase)]
+    # Photon data is only loaded for materials in the model, so fill one cell
+    # with a material too thin to interact
+    mat = openmc.Material()
+    mat.add_element('H', 1.0)
+    mat.set_density('g/cm3', 1.0e-12)
+    cells[-1].fill = mat
+    tally = openmc.Tally()
+    tally.filters = [openmc.CellFilter(cells),
+                     openmc.EnergyFilter([0.0, 1.0e3, 2.0e7])]
+    tally.scores = ['pulse-height']
+    model.tallies = [tally]
+    sp_path = model.run()
+
+    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)