diff --git a/src/state_point.cpp b/src/state_point.cpp index 608918cf428..ff9111a3cb0 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -368,7 +368,11 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) void restart_set_keff() { if (simulation::restart_batch > settings::n_inactive) { - for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) { + // 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); } diff --git a/tests/unit_tests/test_restart_keff.py b/tests/unit_tests/test_restart_keff.py new file mode 100644 index 00000000000..7c5ad8f7d93 --- /dev/null +++ b/tests/unit_tests/test_restart_keff.py @@ -0,0 +1,67 @@ +"""Tests for resuming an eigenvalue simulation from a statepoint.""" + +import re + +import numpy as np +import pytest + +import openmc + + +# One line is printed per generation, e.g. +# 9/1 1.18342 1.18055 +/- 0.00287 +# capturing the batch, the generation and the running average k-effective. +_GENERATION_LINE = re.compile( + r"^\s*(\d+)/(\d+)\s+[\d.]+\s+([\d.]+)\s+\+/-\s+[\d.]+\s*$", re.MULTILINE +) + + +@pytest.mark.parametrize("generations_per_batch", [1, 5]) +def test_restart_keff(run_in_tmpdir, capsys, generations_per_batch): + """A restarted run must resume with the correct accumulated k-effective. + + ``simulation::k_generation`` holds one entry per generation rather than one + per batch, so the accumulation in ``restart_set_keff()`` has to be indexed + by generation. Indexing it by batch sums the wrong entries -- and too few of + them -- leaving the restarted run with a source-normalization k-effective + roughly a factor of ``generations_per_batch`` too low, which over-produces + fission sites and fills the shared fission bank. + """ + model = openmc.examples.pwr_pin_cell() + model.settings.particles = 500 + model.settings.inactive = 3 + model.settings.batches = 8 + model.settings.generations_per_batch = generations_per_batch + + statepoint = model.run() + + with openmc.StatePoint(statepoint) as sp: + assert sp.generations_per_batch == generations_per_batch + n_inactive = sp.n_inactive + n_loaded = len(sp.k_generation) + + # Restart from the statepoint, asking for two more batches + model.settings.batches = 10 + model.export_to_model_xml() + capsys.readouterr() + openmc.run(restart_file=statepoint) + output = capsys.readouterr().out + + # Symptom of a low restart k-effective, asserted separately so that the + # failure mode is obvious + assert "fission bank is full" not in output + + with openmc.StatePoint("statepoint.10.h5") as sp: + k_generation = np.asarray(sp.k_generation) + + average_k = [float(m[2]) for m in _GENERATION_LINE.findall(output)] + assert len(average_k) == 2 * generations_per_batch + + # The running average reported for each generation of the restarted run must + # be the mean over every active generation simulated so far, including the + # ones loaded from the statepoint. The comparison is exact to the precision + # of the printed value. + first_active = generations_per_batch * n_inactive + for i, reported in enumerate(average_k): + expected = k_generation[first_active:n_loaded + i + 1].mean() + assert reported == pytest.approx(expected, abs=1e-5)