From 2199229b15f53a83a5c34bb76928cae7d625c823 Mon Sep 17 00:00:00 2001 From: Eden Rochman Date: Thu, 4 Jun 2026 13:34:36 +0200 Subject: [PATCH] Preserve collision count when splitting particles When weight windows split a particle, the child particles now inherit the parent's collision count. Previously, children were created with n_collision=0 via from_source(), which caused CollisionFilter to incorrectly classify scattered children as uncollided, producing biased tallies. Closes #3916 --- include/openmc/particle_data.h | 1 + src/particle.cpp | 3 +- tests/unit_tests/test_collision_filter_ww.py | 123 +++++++++++++++++++ 3 files changed, 126 insertions(+), 1 deletion(-) create mode 100644 tests/unit_tests/test_collision_filter_ww.py diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index f72948f6eb..786b3f639e 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -55,6 +55,7 @@ struct SourceSite { double wgt_born {1.0}; double wgt_ww_born {-1.0}; int64_t n_split {0}; + int n_collision {0}; }; struct CollisionTrackSite { diff --git a/src/particle.cpp b/src/particle.cpp index 779ae18da9..fa07ea8284 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -135,6 +135,7 @@ void Particle::split(double wgt) bank.wgt_born = wgt_born(); bank.wgt_ww_born = wgt_ww_born(); bank.n_split = n_split(); + bank.n_collision = n_collision(); bank.parent_id = current_work(); if (settings::use_shared_secondary_bank) { bank.progeny_id = n_progeny()++; @@ -150,7 +151,7 @@ void Particle::from_source(const SourceSite* src) surface() = SURFACE_NONE; cell_born() = C_NONE; material() = C_NONE; - n_collision() = 0; + n_collision() = src->n_collision; fission() = false; zero_flux_derivs(); lifetime() = 0.0; diff --git a/tests/unit_tests/test_collision_filter_ww.py b/tests/unit_tests/test_collision_filter_ww.py new file mode 100644 index 0000000000..e762d6b0c2 --- /dev/null +++ b/tests/unit_tests/test_collision_filter_ww.py @@ -0,0 +1,123 @@ +"""Test that CollisionFilter produces unbiased results with weight windows. + +Regression test for https://github.com/openmc-dev/openmc/issues/3916 +When weight windows split a particle, the child particles must inherit +the parent's collision count. Otherwise CollisionFilter bins are biased. +""" + +import numpy as np +import openmc +from tests import cdtemp + + +def test_collision_filter_weight_windows(): + """Verify that total flux from CollisionFilter bins matches + the total flux tally when weight windows are active. + + The test runs two simulations (analog and weight-windowed) and + checks that the uncollided flux (CollisionFilter bin 0) is + statistically consistent between them. Before the fix, split + particles had their collision count reset to zero, inflating the + bin-0 tally. + """ + + model = openmc.Model() + + # Concrete-like material for good scattering + mat = openmc.Material() + mat.set_density('g/cc', 2.3) + mat.add_nuclide('H1', 0.168) + mat.add_nuclide('O16', 0.562) + mat.add_nuclide('Si28', 0.187) + mat.add_nuclide('Ca40', 0.018) + mat.add_nuclide('Fe56', 0.004) + + sphere_inner = openmc.Sphere(r=100) + sphere_outer = openmc.Sphere(r=120, boundary_type='vacuum') + cell = openmc.Cell(fill=mat, region=-sphere_inner) + void = openmc.Cell(region=+sphere_inner & -sphere_outer) + model.geometry = openmc.Geometry([cell, void]) + + settings = openmc.Settings() + settings.run_mode = 'fixed source' + settings.particles = 2000 + settings.batches = 5 + settings.source = openmc.IndependentSource( + space=openmc.stats.Point((0, 0, 0)), + energy=openmc.stats.Discrete([14e6], [1.0]), + ) + settings.survival_biasing = True + + # Weight windows: decreasing values radially to force splitting + ww_mesh = openmc.RegularMesh() + ww_mesh.lower_left = (-120, -120, -120) + ww_mesh.upper_right = (120, 120, 120) + ww_mesh.dimension = (3, 3, 3) + + n_bins = 3 * 3 * 3 + lower_ww = np.logspace(-2, -5, n_bins) + ww = openmc.WeightWindows( + ww_mesh, lower_ww, None, 5.0, + energy_bounds=[0.0, 20e6], + survival_ratio=3.0, + ) + + cell_filter = openmc.CellFilter([cell]) + collision_filter = openmc.CollisionFilter(list(range(10))) + + # Tally 1: total flux + tally_total = openmc.Tally() + tally_total.filters = [cell_filter] + tally_total.scores = ['flux'] + + # Tally 2: flux by collision number + tally_collision = openmc.Tally() + tally_collision.filters = [cell_filter, collision_filter] + tally_collision.scores = ['flux'] + + model.tallies = openmc.Tallies([tally_total, tally_collision]) + + with cdtemp(): + # Run analog (no weight windows) + settings.weight_windows_on = False + model.settings = settings + sp_analog = model.run() + + import os + os.rename(sp_analog, 'statepoint.analog.h5') + + # Run with weight windows + settings.weight_windows = [ww] + settings.weight_windows_on = True + settings.max_history_splits = 1000 + model.settings = settings + sp_ww = model.run() + + # Compare uncollided flux (bin 0) between analog and WW runs + sp_a = openmc.StatePoint('statepoint.analog.h5') + sp_w = openmc.StatePoint(sp_ww) + + # Get uncollided flux (collision bin 0) from both runs + uncollided_analog = sp_a.tallies[tally_collision.id].get_values( + filters=[openmc.CollisionFilter], filter_bins=[(0,)] + ).sum() + uncollided_ww = sp_w.tallies[tally_collision.id].get_values( + filters=[openmc.CollisionFilter], filter_bins=[(0,)] + ).sum() + + total_analog = sp_a.tallies[tally_total.id].mean.sum() + total_ww = sp_w.tallies[tally_total.id].mean.sum() + + # The uncollided fraction should be similar in both runs. + # Before the fix, the WW run had an inflated uncollided fraction + # because split particles were incorrectly tagged as uncollided. + frac_analog = uncollided_analog / total_analog if total_analog > 0 else 0 + frac_ww = uncollided_ww / total_ww if total_ww > 0 else 0 + + # Allow generous tolerance due to Monte Carlo statistics, but + # the bug caused 2-4x inflation so a 50% tolerance catches it + assert frac_ww < frac_analog * 1.5, ( + f"Uncollided fraction with WW ({frac_ww:.4f}) is much larger " + f"than analog ({frac_analog:.4f}), suggesting collision count " + f"is not preserved during particle splitting" + )