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
This commit is contained in:
Eden Rochman 2026-06-04 13:34:36 +02:00
parent 111eb77066
commit 2199229b15
3 changed files with 126 additions and 1 deletions

View file

@ -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 {

View file

@ -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;

View file

@ -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"
)