Delayed fission sites rejected by the time cutoff break analog fission tallies #4179

Bug Description

With a neutron time cutoff (Settings.cutoff['time_neutron']), the analog fission tallies count the fission neutrons of a collision inconsistently, and some of them read past the end of a vector.

create_fission_sites() rejects a delayed fission site whose emission time is past the cutoff. Both copies of the function do this: src/physics.cpp (continuous energy) and src/physics_mg.cpp (multigroup). The continuous-energy version, abridged:

for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) {
  ...
  // Reject site if it exceeds time cutoff
  if (site.delayed_group > 0) {
    double t_cutoff = settings::time_cutoff[site.particle.transport_index()];
    if (site.time > t_cutoff) {
      continue;
    }
  }
  ...
  // Increment the number of neutrons born delayed
  if (site.delayed_group > 0) {
    nu_d[site.delayed_group - 1]++;
  }

  // Write fission particles to nuBank
  NuBank& nu_bank_entry = p.nu_bank().emplace_back();
  ...
}

// If shared fission bank was full, and no fissions could be added,
// set the particle fission flag to false.
if (n_sites_stored == 0) {
  p.fission() = false;
  return;
}
...
nu = n_sites_stored;

// Store the total weight banked for analog fission tallies
p.n_bank() = nu;
p.wgt_bank() = nu / weight;

The continue still runs the loop increment. A rejected site is therefore counted in n_bank and wgt_bank, but it gets no nu_bank entry and no n_delayed_bank count. After a collision with r rejected sites, n_bank == nu_bank.size() + r, and n_sites_stored == 0 never holds when a site was rejected, so fission() stays true even if every site was rejected. The analog scores in src/tallies/tally_scoring.cpp then go wrong in four ways:

1. Out-of-bounds reads. score_fission_eout() is used for analog nu-fission, prompt-nu-fission and delayed-nu-fission with an EnergyoutFilter. It and the analog decay-rate score loop over p.nu_bank(i) for i < p.n_bank(), and nu_bank(i) is unchecked. They read past the end of the vector. Within its capacity they get a stale entry from an earlier collision, and the result is silently wrong. Beyond the capacity they read other memory, and the run can crash.
2. Rejected delayed neutrons are scored as prompt. Analog prompt-nu-fission uses 1 - n_delayed / n_bank.
3. They are missing from delayed-nu-fission. Analog delayed-nu-fission (with or without a DelayedGroupFilter) counts only the delayed neutrons emitted before the cutoff. The track-length and collision estimates of the same score use nu(E) * fission, so they count all of them.
4. fission() stays true with an empty nu_bank. This happens when every site of a collision is rejected, and the loops of item 1 then read from an empty vector. In the runs counted below, nearly every collision that produced fission neutrons produced one, and every collision that lost a site lost all of its sites.

How often it happens was counted on develop with a temporary instrumentation patch, which is not part of any PR. The cutoff was 1e-2 s in every run. Collisions with a rejected site, out of all collisions that produced fission neutrons:

• 300 of 1034 in the first generation of a 2-group MG slab eigenvalue run with delayed fractions 0.1 and 0.2;
• 55 of 9027 in a CE U-235 sphere eigenvalue run;
• 4974 of 4975 in the fixed-source script below.

The rejection came with the delayed neutron emission time (#2898, merged as part of #2992), so it was first released in v0.15.3 and is also in v0.16.0. Before that, a delayed site carried its parent's time and could not be past the cutoff. Without a time cutoff (the default is infinite) no site is rejected, and none of the above happens.

Steps to Reproduce

The script below is a multigroup model that needs no external data. Every fission neutron is delayed (the prompt nu-fission is zero). It runs in fixed-source mode with a time cutoff of 1e-2 s:

import numpy as np
import openmc
from openmc.examples import slab_mg

# Two energy groups and two delayed groups. Every fission neutron is
# delayed: the prompt nu-fission cross section is zero.
groups = openmc.mgxs.EnergyGroups([0.0, 0.625, 20.0e6])
fission = np.array([0.002817, 0.097])
nu_fission = 2.5 * fission
xsdata = openmc.XSdata('mat_1', groups, num_delayed_groups=2)
xsdata.order = 1
xsdata.set_total([0.33588, 0.54628])
xsdata.set_absorption(np.array([0.008708, 0.02518]) + fission)
xsdata.set_scatter_matrix(np.array(
    [[[0.31980, 0.06694], [0.004555, -0.0003972]],
     [[0.00000, 0.00000], [0.424100, 0.05439000]]]))
xsdata.set_fission(fission)
xsdata.set_prompt_nu_fission(np.zeros(2))
xsdata.set_delayed_nu_fission(np.outer([1/3, 2/3], nu_fission))
xsdata.set_chi_prompt([1.0, 0.0])
xsdata.set_chi_delayed([[0.8, 0.2], [0.6, 0.4]])
xsdata.set_decay_rate([0.0127, 0.0317])
xsdata.set_inverse_velocity([1.0e-9, 4.5e-6])
library = openmc.MGXSLibrary(groups, num_delayed_groups=2)
library.add_xsdata(xsdata)
library.export_to_hdf5('2g.h5')

model = slab_mg(mgxslib_name='2g.h5')
model.settings.run_mode = 'fixed source'
model.settings.batches = 5
model.settings.particles = 1000
model.settings.seed = 1
model.settings.cutoff = {'time_neutron': 1.0e-2}

tally = openmc.Tally(name='analog')
tally.scores = ['nu-fission', 'prompt-nu-fission', 'delayed-nu-fission']
tally.estimator = 'analog'
model.tallies = [tally]

sp_path = model.run()
with openmc.StatePoint(sp_path) as sp:
    t = sp.get_tally(name='analog')
    for score in t.scores:
        print(f'{score:20s} {t.get_values(scores=[score]).sum():.4f}')

Output on develop (a5bc348, 2 OpenMP threads):

nu-fission           0.9950
prompt-nu-fission    0.9948
delayed-nu-fission   0.0002

The material has no prompt fission neutrons, so prompt-nu-fission should be 0. Of the 4975 fission neutrons, 4974 are delayed neutrons that would be emitted after the cutoff. They are counted in nu-fission, scored as prompt, and missing from delayed-nu-fission.

To see the out-of-bounds read, add an analog tally with an EnergyoutFilter, i.e. replace model.tallies = [tally] with:

eout = openmc.Tally(name='energyout')
eout.filters = [openmc.EnergyoutFilter([0.0, 0.625, 20.0e6])]
eout.scores = ['nu-fission']
eout.estimator = 'analog'
model.tallies = [tally, eout]

On a RelWithDebInfo build of develop, this run crashed with SIGSEGV in score_fission_eout() (tally_scoring.cpp:405, the first read of the entry). A build configured with -DCMAKE_CXX_FLAGS=-D_GLIBCXX_ASSERTIONS, which checks std::vector::operator[], aborts with:

/usr/include/c++/13/bits/stl_vector.h:1128: std::vector<_Tp, _Alloc>::reference std::vector<_Tp, _Alloc>::operator[](size_type) [with _Tp = openmc::NuBank; _Alloc = std::allocator<openmc::NuBank>; reference = openmc::NuBank&; size_type = long unsigned int]: Assertion '__n < this->size()' failed.

#8  openmc::ParticleData::nu_bank (i=0, ...) at include/openmc/particle_data.h:735
#9  openmc::score_fission_eout (...) at src/tallies/tally_scoring.cpp:402

It is called from score_general_mg(), score_analog_tally_mg() and Particle::event_collide(). Continuous-energy mode aborts the same way, with score_general_ce_analog() as the caller. One such case is the first generation of an eigenvalue run of a bare U-235 sphere with the cutoff and analog fission tallies with an EnergyoutFilter.

A second check needs no assumption about the expected values. In the first generation of an eigenvalue run, no neutron history lasts 1e-2 s, so the transport is the same with and without the cutoff: k and the flux agree to round-off. Every analog fission tally should then agree too. On develop they do not:

• 2-group MG slab (slab_mg with the library above, but with delayed fractions 0.1 and 0.2 and prompt fission), 1000 particles, cutoff versus no cutoff:
  • prompt-nu-fission is 1.034 versus 0.734;
  • delayed-nu-fission is 0 versus 0.3;
  • nu-fission is 1.034 in both runs.
• CE U-235 sphere (16 g/cm3, r = 10 cm), 10000 particles:
  • prompt-nu-fission is 0.9039 versus 0.8984;
  • delayed-nu-fission is 0.0001 versus 0.0056;
  • nu-fission is 0.904 in both runs.

Environment

• OpenMC develop at a5bc348 (openmc --version: 0.16.1-dev57), built from source with CMake (RelWithDebInfo, OPENMC_ENABLE_STRICT_FP=ON); for the assertion, the same with -DCMAKE_CXX_FLAGS=-D_GLIBCXX_ASSERTIONS.
• Ubuntu 24.04.5, GCC 13.3.0, 2 OpenMP threads. History-based and event-based transport are both affected.
• Data: the multigroup library is made by the script. Continuous energy uses the NNDC HDF5 library (nndc_hdf5).
• Reading the code shows the same rejection in v0.15.3 and v0.16.0; nothing was run on those versions.