From 113e4e92830bab9bacd0d54fc5935000315ab290 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 2 Oct 2026 23:42:12 +0000 Subject: [PATCH 1/3] 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) From 80fd7a2fde2848b8b09abbe4c3b6c381fa423ae3 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 12:54:07 +0000 Subject: [PATCH 2/3] Retry the cell search for split photons with pulse-height tallies event_revive_from_secondary locates photons before removing their energy from the pulse-height tally, and that search lost particles on coincident surfaces just like the one in event_calculate_xs. Share the retry between both searches and add a pulse-height regression test. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BK6aD4aznm1HR84bTxGn5v --- src/particle.cpp | 45 ++++++++++++++++--------- tests/unit_tests/test_lost_particles.py | 38 ++++++++++++++++++--- 2 files changed, 64 insertions(+), 19 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index a1778cf3520..5dae9801ef0 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -40,6 +40,34 @@ 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. As when crossing a surface, move the particle forward a tiny bit and +//! search again. +// +//! \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.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,20 +251,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) { - 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) { + if (!find_cell_from_surface(*this)) { mark_as_lost( "Could not find the cell containing particle " + std::to_string(id())); return; @@ -533,7 +548,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; diff --git a/tests/unit_tests/test_lost_particles.py b/tests/unit_tests/test_lost_particles.py index a9f3ebea094..39d3c3a791c 100644 --- a/tests/unit_tests/test_lost_particles.py +++ b/tests/unit_tests/test_lost_particles.py @@ -49,9 +49,7 @@ def test_max_write_lost_particles(model: openmc.Model, run_in_tmpdir): -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)""" +def _coincident_surfaces_model(particle): s1 = openmc.XPlane(-57452.33336021505) s2 = openmc.XPlane(43403.32187479845) s3 = openmc.XPlane(62125.04976135607) @@ -71,6 +69,7 @@ def test_split_on_coincident_surfaces(run_in_tmpdir): 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( (0.9182609440159194, 0.39597580569397484, 0.0)), @@ -84,13 +83,44 @@ def test_split_on_coincident_surfaces(run_in_tmpdir): 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]) + energy_bounds=[0.0, 2.0e7], particle_type=particle) model.settings.weight_window_checkpoints = { 'collision': False, 'surface': True} + return model + +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)""" + model = _coincident_surfaces_model('neutron') 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) From 6b4f8ec326951a9d13eb3f082a2c3ce9310f01e4 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 14:18:25 +0000 Subject: [PATCH 3/3] Keep the surface sense during the coincident-surface cell retry Keep the particle's surface set while retrying the cell search, so the side of that surface still comes from the particle's surface attribute rather than its position. On a near-grazing direction the small move may not get past roundoff, and the particle must then be lost rather than placed in the cell on the other side. Add a near-grazing case to the split test. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BK6aD4aznm1HR84bTxGn5v --- src/particle.cpp | 8 +++++--- tests/unit_tests/test_lost_particles.py | 22 ++++++++++++++-------- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 5dae9801ef0..55bbc302043 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -48,8 +48,11 @@ namespace { //! 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. +//! 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 @@ -58,7 +61,6 @@ bool find_cell_from_surface(Particle& p, bool verbose = false) { bool found = exhaustive_find_cell(p, verbose); if (!found && p.surface() != SURFACE_NONE) { - p.surface() = SURFACE_NONE; p.n_coord() = 1; p.r() += TINY_BIT * p.u(); found = exhaustive_find_cell(p, verbose); diff --git a/tests/unit_tests/test_lost_particles.py b/tests/unit_tests/test_lost_particles.py index 39d3c3a791c..f5fb08d9f4f 100644 --- a/tests/unit_tests/test_lost_particles.py +++ b/tests/unit_tests/test_lost_particles.py @@ -1,3 +1,4 @@ +import math from pathlib import Path import openmc @@ -49,7 +50,9 @@ def test_max_write_lost_particles(model: openmc.Model, run_in_tmpdir): -def _coincident_surfaces_model(particle): +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) @@ -71,15 +74,14 @@ def _coincident_surfaces_model(particle): model.settings.source = openmc.IndependentSource( particle=particle, space=openmc.stats.Point((-85626.45049221347, 0.0, 0.0)), - angle=openmc.stats.Monodirectional( - (0.9182609440159194, 0.39597580569397484, 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.0e5, -1.0) - mesh.upper_right = (2.2e5, 1.0e5, 1.0) + 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], @@ -89,10 +91,14 @@ def _coincident_surfaces_model(particle): return model -def test_split_on_coincident_surfaces(run_in_tmpdir): +@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)""" - model = _coincident_surfaces_model('neutron') + 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