Preserve collision count when splitting particles

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 {

openmc/lib/core.py

@@ -34,7 +34,8 @@ class _SourceSite(Structure):
                 ('progeny_id', c_int64),
                 ('wgt_born', c_double),
                 ('wgt_ww_born', c_double),
-                ('n_split', c_int64)]
+                ('n_split', c_int64),
+                ('n_collision', c_int)]
 
 # Define input type for numpy arrays that will be passed into C++ functions
 # Must be an int or double array, with single dimension that is contiguous

src/initialize.cpp

@@ -161,7 +161,7 @@ void initialize_mpi(MPI_Comm intracomm)
 
   // Create bank datatype
   SourceSite b;
-  MPI_Aint disp[14];
+  MPI_Aint disp[15];
   MPI_Get_address(&b.r, &disp[0]);
   MPI_Get_address(&b.u, &disp[1]);
   MPI_Get_address(&b.E, &disp[2]);
@@ -176,12 +176,13 @@ void initialize_mpi(MPI_Comm intracomm)
   MPI_Get_address(&b.wgt_born, &disp[11]);
   MPI_Get_address(&b.wgt_ww_born, &disp[12]);
   MPI_Get_address(&b.n_split, &disp[13]);
-  for (int i = 13; i >= 0; --i) {
+  MPI_Get_address(&b.n_collision, &disp[14]);
+  for (int i = 14; i >= 0; --i) {
     disp[i] -= disp[0];
   }
 
   // Block counts for each field
-  int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
+  int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
 
   // Types for each field
   MPI_Datatype types[] = {
@@ -198,10 +199,11 @@ void initialize_mpi(MPI_Comm intracomm)
     MPI_INT64_T, // progeny_id
     MPI_DOUBLE,  // wgt_born
     MPI_DOUBLE,  // wgt_ww_born
-    MPI_INT64_T  // n_split
+    MPI_INT64_T, // n_split
+    MPI_INT      // n_collision
   };
 
-  MPI_Type_create_struct(14, blocks, disp, types, &mpi::source_site);
+  MPI_Type_create_struct(15, blocks, disp, types, &mpi::source_site);
   MPI_Type_commit(&mpi::source_site);
 
   CollisionTrackSite bc;

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;

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