CollisionFilter produces biased results when weight windows are active

[Edvree-dot]

Bug Description

CollisionFilter([0]) produces biased (inflated) flux tallies when weight windows are active. The bias is deterministic, it does not decrease with more particles, and converges to the wrong answer. Total flux tallies without CollisionFilter are unbiased in all configurations.

The inflation scales with the aggressiveness of the weight windows:

Configuration	Uncollided/Free-field	Analytical	Deviation
Analog (no VR)	0.1438	0.1602	-10.2% (statistical)
Survival biasing only	0.1438	0.1602	-10.2% (statistical)
Weight windows only	0.3389	0.1602	+111.6%
WW + survival biasing	0.6034	0.1602	+276.7%

Total flux (without CollisionFilter) is correct (~0.65 ratio to free-field) in all four configurations.

Suspected cause: When weight windows split a particle, child particles are created via create_secondary in weight_windows.cpp with collision_count = 0, regardless of the parent's actual collision count. If a parent particle has scattered multiple times and is then split, the children are tagged as having zero collisions. CollisionFilter([0]) then incorrectly scores these children as uncollided particles.

Impact:

• Any user combining CollisionFilter with weight windows will get biased results for individual collision-count bins
• Uncollided flux measurements are unreliable when weight windows are active
• Total flux tallies (without CollisionFilter) are NOT affected

Steps to Reproduce

The script below runs five simulations: one free-field reference (no slab) and four configurations with a 10 cm carbon slab, varying weight windows and survival biasing. It compares the CollisionFilter([0]) tally against the free-field flux and the known analytical transmission (~0.16).

#!/usr/bin/env python3
"""
Minimal Working Example: CollisionFilter bias with weight windows
Setup:

10 cm carbon slab (1.7 g/cm3) at z=50..60 cm
Detector sphere (r=30 cm) at z=200 cm
Watt fission spectrum source, cosine angular distribution
Five runs: free-field (no slab), then four slab configurations

varying weight windows (WW) and survival biasing (SB)

Expected: Uncollided/free-field ratio ~ 0.16 for all configurations.

Observed: WW-only gives ~0.34 (+112%), WW+SB gives ~0.60 (+277%).

"""

import openmc

import numpy as np

import os

from pathlib import Path
DENSITY = 1.7

SLAB_THICK = 10.0

SLAB_Z = 50.0

DET_Z = 200.0

DET_R = 30.0

SOURCE_R = 100.0

N_PARTICLES = 500_000

N_BATCHES = 10

A_WATT = 0.988e6

B_WATT = 2.249e-6
Spectrum-averaged transmission through 10 cm carbon at 1.7 g/cm3
(computed from ENDF/B-VIII.0 total cross sections sampled over Watt spectrum)
T_ANALYTICAL = 0.1602
def build_and_run(with_slab, name, use_ww=False, use_sb=False):

carbon = openmc.Material(name='carbon')

carbon.add_element('C', 1.0, 'ao')

carbon.set_density('g/cm3', DENSITY)

materials = openmc.Materials([carbon])
z_lo = openmc.ZPlane(z0=SLAB_Z)
z_hi = openmc.ZPlane(z0=SLAB_Z + SLAB_THICK)
slab_region = +z_lo & -z_hi
det_sph = openmc.Sphere(z0=DET_Z, r=DET_R)
det_cell = openmc.Cell(name='detector', fill=None, region=-det_sph)
outer = openmc.Sphere(r=DET_Z + 200, boundary_type='vacuum')

if with_slab:
    slab_cell = openmc.Cell(name='slab', fill=carbon, region=slab_region & -outer)
    void_cell = openmc.Cell(name='void', fill=None,
                            region=-outer & ~slab_region & +det_sph)
    cells = [slab_cell, det_cell, void_cell]
else:
    void_cell = openmc.Cell(name='void', fill=None, region=-outer & +det_sph)
    cells = [det_cell, void_cell]

geometry = openmc.Geometry(cells)

source = openmc.IndependentSource()
source.space = openmc.stats.CylindricalIndependent(
    r=openmc.stats.PowerLaw(0, SOURCE_R, 1),
    phi=openmc.stats.Uniform(0, 2 * np.pi),
    z=openmc.stats.Discrete([0.0], [1.0]),
    origin=(0.0, 0.0, 0.0))
source.energy = openmc.stats.Watt(a=A_WATT, b=B_WATT)
mu = np.linspace(0, 1, 500)
source.angle = openmc.stats.PolarAzimuthal(
    mu=openmc.stats.Tabular(mu, 2 * mu, interpolation='linear-linear'),
    phi=openmc.stats.Uniform(0, 2 * np.pi),
    reference_uvw=(0.0, 0.0, 1.0))

settings = openmc.Settings()
settings.source = source
settings.run_mode = 'fixed source'
settings.batches = N_BATCHES
settings.particles = N_PARTICLES
settings.survival_biasing = use_sb

if use_ww and with_slab:
    z_edges = list(np.linspace(0, SLAB_Z, 5))
    z_edges += list(np.linspace(SLAB_Z, SLAB_Z + SLAB_THICK, 10))
    z_edges += list(np.linspace(SLAB_Z + SLAB_THICK, DET_Z + 50, 5))
    z_edges = sorted(set(z_edges))
    n_z = len(z_edges) - 1
    lb = np.ones(n_z)
    for i in range(n_z):
        zm = (z_edges[i] + z_edges[i + 1]) / 2
        d = max(0, min(zm, SLAB_Z + SLAB_THICK) - SLAB_Z) if zm > SLAB_Z else 0
        geom = (50.0 / max(zm, 10)) ** 2
        lb[i] = geom * np.exp(-0.15 * d)
    lb /= lb[0]
    lb = np.maximum(lb, 1e-8)
    mesh = openmc.RectilinearMesh()
    mesh.x_grid = [-600, 600]
    mesh.y_grid = [-600, 600]
    mesh.z_grid = z_edges
    ww = openmc.WeightWindows(
        mesh=mesh,
        lower_ww_bounds=lb.reshape(1, 1, n_z, 1),
        upper_ww_bounds=(lb * 5).reshape(1, 1, n_z, 1),
        energy_bounds=np.array([0.0, 20e6]),
        particle_type='neutron', survival_ratio=3.0)
    settings.weight_windows = ww

tallies = openmc.Tallies()
t_uncoll = openmc.Tally(name='uncollided')
t_uncoll.filters = [openmc.CellFilter([det_cell]), openmc.CollisionFilter([0])]
t_uncoll.scores = ['flux']
tallies.append(t_uncoll)
t_total = openmc.Tally(name='total')
t_total.filters = [openmc.CellFilter([det_cell])]
t_total.scores = ['flux']
tallies.append(t_total)

model = openmc.Model(geometry, materials, settings, tallies)
out_dir = Path(f"output/collisionfilter_mwe_{name}")
out_dir.mkdir(parents=True, exist_ok=True)
orig = os.getcwd()
os.chdir(out_dir)
model.export_to_xml()
sp_file = model.run()
sp = openmc.StatePoint(sp_file)
u = sp.get_tally(name='uncollided').mean.flatten()[0]
t = sp.get_tally(name='total').mean.flatten()[0]
sp.close()
os.chdir(orig)
return u, t

def main():

print("CollisionFilter + Weight Windows Bias Test")

print(f"  {SLAB_THICK} cm carbon slab, {N_PARTICLES * N_BATCHES:,} total particles")

print(f"  Analytical uncollided transmission: {T_ANALYTICAL:.4f}\n")
_, t_free = build_and_run(with_slab=False, name="free")
print(f"  Free-field total flux: {t_free:.4e}\n")

configs = [
    ("analog",     False, False, "No variance reduction"),
    ("sb_only",    False, True,  "Survival biasing only"),
    ("ww_only",    True,  False, "Weight windows only"),
    ("ww_plus_sb", True,  True,  "WW + survival biasing"),
]

print(f"{'Config':<26} {'Uncoll/Free':>14} {'Total/Free':>14} {'Expected':>14} {'Deviation':>12}")
print("-" * 84)
for label, use_ww, use_sb, desc in configs:
    u, t = build_and_run(with_slab=True, name=label, use_ww=use_ww, use_sb=use_sb)
    ru = u / t_free if t_free > 0 else float('nan')
    rt = t / t_free if t_free > 0 else float('nan')
    dev = (ru - T_ANALYTICAL) / T_ANALYTICAL * 100
    print(f"{desc:<26} {ru:14.4f} {rt:14.4f} {T_ANALYTICAL:14.4f} {dev:+11.1f}%")

print("\nExpected: all Uncoll/Free values ~ 0.16")
print("Observed: WW configurations are inflated by 2-4x")

if name == "main":

main()

Runtime is approximately 3-5 minutes (5 million total particles across 5 runs).

Environment

• OpenMC 0.15.3 (conda-forge)
• ENDF/B-VIII.0 cross sections
• Python 3.12
• Ubuntu 24.04 (WSL2)

[Edvree-dot]

collisionfilter_ww_mwe.py

The .py file for reference.