Fix analog fission tallies for sites rejected by the time cutoff

[ondrejch]

Description

With a neutron time cutoff, create_fission_sites() rejects a delayed fission site whose emission time is past the cutoff. This happens in the continuous-energy version (src/physics.cpp) and in the multigroup version (src/physics_mg.cpp). The loop counter that becomes n_bank is incremented by the continue as well. The rejected site is therefore counted in n_bank and wgt_bank, but it gets no nu_bank entry and no n_delayed_bank count. As described in #4179, the analog fission tallies then:

• read past the end of nu_bank, in score_fission_eout() and in the analog decay-rate score;
• score the rejected delayed neutrons as prompt;
• see fission() true with an empty nu_bank when every site of a collision was rejected.

The rejection was added with the delayed neutron emission time (#2898, merged as part of #2992), so it was first released in v0.15.3. Before that, a delayed site had its parent's time and was never past the cutoff, so every fission neutron produced was banked and counted. The rejection itself has to stay. Since #3526, Particle::event_advance() limits a flight to (t_cutoff - time) * speed. A neutron banked past the cutoff would therefore be moved a negative distance and score a negative track length before it is killed.

The change

A rejected delayed neutron is now recorded like any other fission neutron, in nu_bank and in the delayed-group counts. It is still not banked:

• it takes no progeny ID;
• it does not call ifp();
• it enters neither the fission bank nor the secondary bank;
• in multigroup mode, it does not set p.delayed_group().

The loop index no longer doubles as the counter. n_bank and wgt_bank are taken from nu_bank.size(), so they agree with nu_bank and n_delayed_bank whenever fission() is set.

The full-bank path is unchanged. The site that cannot be banked and all later ones are dropped, with the existing warning. Neutrons rejected before that point are recorded, so n_bank is the same as on develop.

The rejection consumes no random numbers and nothing is resampled, so the random number stream does not change. tally_scoring.cpp, file formats, MPI types and SourceSite are not touched. In particle_data.h, the comments on n_bank, wgt_bank and n_delayed_bank now say what they count. One sentence under time_neutron in settings.rst, and a matching one in the Settings.cutoff docstring, states the rule.

For a collision with a rejected site:

Analog score	develop	this PR
nu-fission	counts the rejected neutrons	unchanged
prompt-nu-fission	counts the rejected delayed neutrons as prompt	prompt neutrons only
delayed-nu-fission, with or without DelayedGroupFilter	delayed neutrons emitted before the cutoff only	all delayed neutrons produced
any of the three with an EnergyoutFilter	reads past the end of nu_bank	one entry per neutron counted above
decay-rate	reads past the end of nu_bank	all delayed neutrons produced
fission() gate, IFP scores, absorption event MT	collision is a fission	unchanged

Several things are also unchanged:

• the fission and secondary banks;
• the progeny IDs and the sort_bank() check;
• the k estimators and combing;
• every non-analog estimator.

ParticleProductionFilter counts the secondaries that are banked. In fixed-source runs with a cutoff it therefore stays below analog nu-fission by the rejected neutrons, as on develop. The default cutoff is infinite, so without a time_neutron cutoff no site is rejected and nothing changes (see Testing).

Counting the rejected neutrons, and the alternative

The alternative, (a), counts only the kept sites. n_bank, wgt_bank and n_delayed_bank would come from the banked sites, and fission() would be false when all sites were rejected. With the structure of this PR, (a) is the same patch without the record_fission_neutron(site) call before the continue. I chose to count the rejected neutrons, for these reasons:

• The same quantity for every estimator. The track-length and collision estimates of nu-fission, prompt-nu-fission, delayed-nu-fission and decay-rate use nu(E) * fission, which includes delayed neutrons whatever the cutoff. With this PR the analog estimates estimate the same quantities (see Testing). Under (a), the analog delayed scores would be scaled by the probability that the neutron is emitted before the cutoff. The same score would then mean different things for different estimators, and OpenMC switches a tally to the analog estimator when it has an EnergyoutFilter.
• Time bins. Analog fission scores are binned at the collision time. Under (a), a TimeFilter bin well before the cutoff would still depend on where the cutoff is, because its delayed neutrons are emitted later.
• The rejection is a shortcut for "bank the neutron and kill it at birth", which Fix a bug in time cutoff behavior #3526 makes impossible to do literally. A neutron killed by the cutoff after its birth is counted where it was produced. A neutron killed at birth should be too.
• The behaviour before Combing for fission site sampling, and delayed neutron emission time #2992. Every fission neutron produced was counted then. This PR keeps that, and delayed-nu-fission keeps the meaning it had.
• A smaller behaviour change. Several things stay as on develop: analog nu-fission, fission(), the IFP scores, and the absorption event MT (seen by collision tracks, ReactionFilter and the MT scores). Under (a), all of them would change for collisions whose sites were all rejected.
• Physics switches live in nu. Where OpenMC removes delayed neutrons on purpose (create_delayed_neutrons = False), it does so in Nuclide::nu(), and every estimator agrees. The time cutoff is a transport truncation, not a physics switch.

The case for (a):

• It is the smallest change.
• The comments in tally_scoring.cpp describe the analog estimate as "the number of particles that were banked in the fission bank". That matches (a) literally. This PR leaves tally_scoring.cpp unchanged.
• It agrees with ParticleProductionFilter, and with the full-bank path, which drops the sites it cannot bank without counting them.
• It matches how DNP drift (regular mesh only) #3944 counts its new rejections (below).

If reviewers prefer (a), the code change is only the two record_fission_neutron(site) calls in the rejection branches. The tests change more: every case that counts the rejected neutrons. With that variant, 16 of the 18 cases fail; only the two continuous-energy full-bank cases, which run without a cutoff, still pass. For example, the eigenvalue test would have to require equal prompt-nu-fission, and smaller delayed-nu-fission and decay-rate, with the cutoff.

Related PRs and issues

• DNP drift (regular mesh only) #3944 (open, DNP drift) rewrites the same continuous-energy loop:
  • It also separates the loop index from the counter, but it increments n_sites_stored only after a site is banked. That is (a) counting, and on its own it would make the continuous-energy counts consistent that way.
  • It does not touch physics_mg.cpp.
  • It adds two more rejections of delayed sites before the time-cutoff check: a precursor that leaves the model, and a site that is not usable for the next generation.
  • The two PRs conflict textually in three places: the loop header, the rejection block, and the end of the loop body.
  • Whichever lands second has to settle how those two rejections count. Either they are recorded as produced, like the time-cutoff one here (with a record_fission_neutron(site) call before their continue), or they are not counted at all.
• Bundle each IFP data stream with the flag and banks that govern it #4097 (open) changes the IFP call inside the same loop (settings::ifp_on() to ifp_on()). This PR does not touch that line.
• Issue Adding a score for the particles killed by the cutoff #2531 (open) asks for a score of the particles killed by the cutoff, which is a different request.

Testing

The new tests/unit_tests/test_time_cutoff_fission.py has three test configurations. Each runs in history and event mode, so there are six cases:

• test_time_cutoff_eigenvalue[mg|ce] runs the first generation of an eigenvalue problem with and without a cutoff of 1e-2 s.
  • Models: for MG, slab_mg with a 2-group library with delayed fractions 0.1 and 0.2 (1000 particles); for CE, a bare U-235 sphere (10000 particles).
  • No history lasts 1e-2 s, so the transport is the same in both runs, and the test checks k_generation and the flux.
  • Every analog fission tally must then be equal in the two runs: unfiltered, with an EnergyoutFilter, and with a DelayedGroupFilter (including decay-rate).
  • In each run, the EnergyoutFilter tally must sum to the unfiltered one. In MG, decay-rate must equal lambda times delayed-nu-fission.
  • The source bank of the run without the cutoff must hold delayed sites past the cutoff, and the run with the cutoff must hold none.
• test_time_cutoff_all_delayed is an MG fixed-source problem in which every fission neutron is delayed. prompt-nu-fission must be exactly 0, and delayed-nu-fission must equal nu-fission, with and without an EnergyoutFilter.

On develop all six cases fail:

• In history mode, the MG eigenvalue case gives prompt-nu-fission 1.034 with the cutoff versus 0.734 without, and delayed-nu-fission 0 versus 0.3.
• The CE case gives 0.9039 versus 0.8984, and 0.0001 versus 0.0056.
• The event-mode eigenvalue cases (in the run with the cutoff) and both fixed-source cases crash with SIGSEGV in score_fission_eout().

With this PR all six pass, in 0.75 s at 2 threads.

A second commit adds six more configurations to the same module, again in both modes, so twelve more cases. They do not compare two runs:

• test_time_cutoff_one_neutron_per_collision uses an MG fixed-source library in which every collision produces exactly one delayed neutron, emitted long after the cutoff.
  • Analog nu-fission and delayed-nu-fission must equal the collision estimates, and prompt-nu-fission must be 0.
  • No fission neutron may enter the simulation: ParticleProductionFilter scores nothing, and the track-length flux must be positive.
  • Collision tracks must show delayed group 0 and progeny ID 0 for every collision.
• test_time_cutoff_fixed_source_ce uses a CE fixed-source U-235 sphere. The track-length flux must be positive, and each analog estimate must be within 50 % of the track-length one.
• test_time_cutoff_full_fission_bank[mg|ce] fills the fission bank in the first generation: MG with eight neutrons per collision, CE with a reflected Cf-249 cube.
  • Without the cutoff, analog nu-fission must equal the number of banked sites, three per source particle. In every run, prompt plus delayed must equal nu-fission.
  • In MG with the cutoff, analog prompt-nu-fission must equal the banked sites, and delayed-nu-fission, the rejected neutrons, must be positive.
  • A third commit, which also wraps long lines, adds a check to that MG run. Analog delayed-nu-fission must stay below 0.85 times the collision estimate, because the neutrons the collision would have produced after the site that failed to enter the bank are neither sampled nor counted.
• test_time_cutoff_ufs[mg|ce] uses uniform fission site weighting. The EnergyoutFilter tallies must sum to the unfiltered ones. In MG, decay-rate must also equal lambda times delayed-nu-fission.

With all three commits, all 18 cases pass with this PR, in 2.0 s at 2 threads. They also pass with -D_GLIBCXX_ASSERTIONS and with 2 MPI ranks. On develop, 16 fail. The two that pass there are the CE full-bank cases, which run without a cutoff.

Other checks were run on Ubuntu 24.04.5 with GCC 13.3.0, built RelWithDebInfo with OPENMC_ENABLE_STRICT_FP=ON, with 2 OpenMP threads unless noted and the NNDC data:

• Assertions. Built with -D_GLIBCXX_ASSERTIONS, develop aborts at nu_bank(i=0) in score_fission_eout() in all six cases. With this PR all six pass, with no assertion.

• Invariants. A temporary instrumentation patch (not included) checked these after every create_fission_sites() call:

  • n_bank == nu_bank.size();
  • n_delayed_bank of each group equals the number of nu_bank entries of that group;
  • wgt_bank equals the sum of the nu_bank weights;
  • fission() == (n_bank > 0);
  • the progeny count equals the number of sites banked.

  With this PR there was no violation, in history and event mode, in the test models and in the fixed-source models below (seed 1). There was none either in two extra models with several neutrons per collision; in one of them the fission bank fills up and more than 11000 calls find it full. On develop, n_bank == nu_bank.size() fails in every call with a rejected site: 300 of the 1034 calls that produced fission neutrons in the MG eigenvalue case with the cutoff, and 55 of 9027 in CE. In every run whose fission bank does not fill, the numbers of calls and of sampled and rejected sites are identical with and without this PR. With a full bank at 2 threads, they vary from run to run on either build.

• Analog versus track-length, in fixed source with the cutoff:

  • Models: MG (slab_mg with the library with delayed fractions 0.1 and 0.2) and CE (a bare U-235 sphere, r = 8 cm, with a 2 MeV point source at t = 0). Each ran 10 seeds of 10 batches x 1000 particles.
  • Method: for each seed, r = analog / track-length - 1, for nu-fission, prompt-nu-fission and delayed-nu-fission. Then a t-test of mean(r) = 0, with Šidák control of the six tests at a family level of 0.01 (|t| <= 4.42).
  • With this PR, the largest |t| is 0.997.
  • On develop, analog delayed-nu-fission is 99.97 % (MG) and 99.68 % (CE) below the track-length estimate, and prompt-nu-fission is 42.8 % (MG) and 0.61 % (CE) above it. nu-fission agrees with the track-length estimate on both builds.

• Without a time cutoff:

  • pytest --update over tests/regression_tests regenerates every results_true.dat and inputs_true.dat byte for byte the same on develop and with this PR.
  • Rerun at one thread, the other generated files of the thread-order-dependent tests agree too. The exceptions are the version and commit attributes, the depletion wall-clock time, and surface_source_write/case-a01, which forces two threads and whose bank holds the same sites in a different order.
  • Against the committed results, and over tests/unit_tests, both builds give the same passes and failures, apart from the new module. The failures on both come from the test environment: no installed OpenMC for cpp_driver and the dlopen tests, no openpyxl, and no OPENMC_ENDF_DATA.

• Event mode. The ten mg_* regression tests pass with --event.

• MPI (OpenMPI 4.1.6, 2 ranks of 1 thread). The checks of the MG eigenvalue test hold in history and event mode. mg_basic, mg_basic_delayed, mg_tallies and mg_fixed_source_ww_fission_shared_secondary pass with --mpi --mpi-np 2.

• Other suites. ctest passes for the serial and the MPI builds. test_time_filter.py, test_filter_delayedgroup.py, test_settings.py, test_mg_inverse_velocity.py and test_mesh.py pass. clang-format 18.1.8 reports no changes.

Fixes #4179

AI Assistance

Code, tests, analysis and this text:

• Harness: Claude Code
• Model: Claude Opus 5.5 (claude-opus-5-5)
• Reasoning effort: xhigh

Review of the change:

• Harness: OpenCode
• Model: grok-4.7
• Reasoning effort: not recorded

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)