From 113e4e92830bab9bacd0d54fc5935000315ab290 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 2 Oct 2026 23:42:12 +0000 Subject: [PATCH] Retry locating particles that start on a surface coincident with another Particles banked on a surface, such as particles split by weight windows at a surface crossing, are only known to be on that one surface. If another surface coincides with it, roundoff in the banked position can place the particle on the wrong side of the other surface, beyond FP_COINCIDENT for large coordinates, so that no cell contains it and the particle is lost (GitHub issue #4140). When a particle with a surface set cannot be located at the start of its history, move it forward by TINY_BIT with the surface cleared and search again, as Particle::cross_surface already does after crossing a surface. Fixes #4140 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014RN7JroEkrbkLKVHbMDap9 --- src/particle.cpp | 15 +++++++- tests/unit_tests/test_lost_particles.py | 47 +++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 1 deletion(-) diff --git a/src/particle.cpp b/src/particle.cpp index 998b71883a8..a1778cf3520 100644 --- a/src/particle.cpp +++ b/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; diff --git a/tests/unit_tests/test_lost_particles.py b/tests/unit_tests/test_lost_particles.py index 478ddbaf0f7..a9f3ebea094 100644 --- a/tests/unit_tests/test_lost_particles.py +++ b/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)