Fix restart k-effective accumulation when generations_per_batch > 1

[GuySten]

Fix restart k-effective accumulation when generations_per_batch > 1

Description

restart_set_keff() in src/state_point.cpp mixes batch indices and generation
indices. simulation::k_generation contains one entry per generation, but the
accumulation loop runs over batch indices:

for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) {
  simulation::k_sum[0] += simulation::k_generation[i];
  simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
}
int n = settings::gen_per_batch * simulation::n_realizations;
simulation::keff = simulation::k_sum[0] / n;

The normalization n is already a generation count, so the numerator and the
denominator only agree when gen_per_batch == 1 — which is why this has gone
unnoticed. With gen_per_batch > 1 the loop sums restart_batch - n_inactive
entries taken from the inactive part of k_generation and then divides by
gen_per_batch * n_realizations.

For the reported case (batches = 300, inactive = 100,
generations_per_batch = 5), the loop sums k_generation[100:300] — 200
generations that all fall inside the inactive region — and divides by 1000,
giving a restart simulation::keff of about 0.26 against a true active-batch
k-effective of about 1.3.

Two consequences:

1. simulation::keff normalizes fission site production, so a k-effective that
   is ~5x too low produces ~5x too many fission sites on the first restart
   generations. The shared fission bank overflows and OpenMC repeatedly emits
   WARNING: The shared fission bank is full. Additional fission sites created in this generation will not be banked. Results may be non-deterministic.
2. Less visibly, the bad k_sum[0] / k_sum[1] values are carried into
   calculate_average_keff(), so the mean and standard deviation reported for
   the continued active batches are computed from generations that were never
   supposed to be part of the active tally.

Fix

Index k_generation by generation, consistent with calculate_average_keff(),
which accumulates one entry per active generation and normalizes by
gen_per_batch * n_realizations + current_gen:

if (simulation::restart_batch > settings::n_inactive) {
  // k_generation has one entry per generation rather than per batch, so the
  // active portion of the array must be indexed by generation
  int i_start = settings::gen_per_batch * settings::n_inactive;
  int i_end = settings::gen_per_batch * simulation::restart_batch;
  for (int i = i_start; i < i_end; ++i) {
    simulation::k_sum[0] += simulation::k_generation[i];
    simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2);
  }
  int n = settings::gen_per_batch * simulation::n_realizations;
  simulation::keff = simulation::k_sum[0] / n;
}

The loop now contributes exactly gen_per_batch * (restart_batch - n_inactive)
terms, matching the gen_per_batch * n_realizations denominator. The
gen_per_batch == 1 path is bit-for-bit unchanged, so existing restart results
are unaffected.

Reproducing

Run an eigenvalue model with generations_per_batch = 5, inactive = 100,
batches = 300, then restart from statepoint.300.h5 with batches = 500.
Before the fix, the restart begins with a source-normalization k-effective far
from the converged value and the fission bank fills immediately; after the fix,
the restarted run picks up at the statepoint's active-batch k-effective and the
warnings disappear.

Test

tests/unit_tests/test_restart_keff.py runs a short eigenvalue calculation to
batch 8, restarts from that statepoint for two more batches, and checks the
running average k-effective reported for every restarted generation against the
mean over all active generations simulated so far, taken from the restarted
run's own statepoint. The comparison is exact to the precision of the printed
value. It also asserts that no fission-bank-full warning is emitted, so the
reported symptom has its own assertion.

The test is parametrized over generations_per_batch of 1 and 5. The second
case fails on develop -- the reported average is low by about a factor of
five and the bank warning fires on every generation -- while the first guards
the unchanged single-generation path.

The obvious stronger test, asserting that a restarted run reproduces the
uninterrupted run bit-for-bit, is not currently possible:
openmc_statepoint_load() restores the seed and stride but not
simulation::total_gen, which enters the per-particle seed alongside
overall_generation(). A restarted run therefore draws a different random
stream than the original run would have at the same generation. Making restart
bit-reproducible would change the results of every restart run and seems worth
a separate discussion.

Credit

Reported and diagnosed by @A_Bohemian on the OpenMC forum:
https://openmc.discourse.group/t/possible-restart-bug-with-generations-per-batch-1-fission-bank-fills-after-loading-statepoint/6456

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)