From 2e815c47765a9085aaee71ae4e4afa48bb15bc3f Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 7 Oct 2026 03:55:33 -0400 Subject: [PATCH 1/4] Fix analog fission tallies for sites rejected by the time cutoff With a neutron time cutoff, create_fission_sites() rejects a delayed fission site whose emission time exceeds the cutoff, since the neutron would be killed at birth. The rejected site was still counted in n_bank and wgt_bank, because the loop counter that becomes n_bank is also incremented when a site is rejected, but it got no nu_bank entry and no n_delayed_bank count. Analog fission tallies then read past the end of nu_bank in score_fission_eout() and in the analog decay-rate scores, scored the rejected delayed neutrons as prompt neutrons, since the prompt fraction is taken from n_bank and n_delayed_bank, and left fission() true with an empty nu_bank when every site of a collision was rejected. The continuous-energy and the multigroup versions of create_fission_sites() both had the problem. A rejected delayed neutron is now recorded like any other fission neutron, in nu_bank and in the delayed group counts, but it is still not banked: it takes no progeny ID and enters neither the fission bank nor the secondary bank. n_bank and wgt_bank are taken from the number of neutrons recorded, so n_bank equals the size of nu_bank whenever fission() is set. Analog nu-fission, prompt-nu-fission, delayed-nu-fission and decay-rate scores then count every fission neutron produced, as the track-length and collision estimators do, whether or not a time cutoff is set. Nothing changes without a time cutoff, and the random number stream is the same. Add tests/unit_tests/test_time_cutoff_fission.py. In the first generation of a multigroup and of a continuous-energy eigenvalue run, every analog fission tally must be the same with and without a time cutoff, and a tally with an EnergyoutFilter must sum to the same tally without the filter. In a multigroup fixed-source problem in which every fission neutron is delayed, no prompt-nu-fission may be scored. --- docs/source/io_formats/settings.rst | 4 +- include/openmc/particle_data.h | 11 +- openmc/settings.py | 5 +- src/physics.cpp | 49 ++-- src/physics_mg.cpp | 49 ++-- tests/unit_tests/test_time_cutoff_fission.py | 224 +++++++++++++++++++ 6 files changed, 292 insertions(+), 50 deletions(-) create mode 100644 tests/unit_tests/test_time_cutoff_fission.py diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 01021e8d22d..18fa60b3e31 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -206,7 +206,9 @@ time. *Default*: 0.0 :time_neutron - The time above which neutrons will be killed. + The time above which neutrons will be killed. Delayed fission neutrons that + would be emitted after this time are not banked, but analog fission tallies + count them as produced, as the track-length and collision estimators do. *Default*: Infinity diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 24b7eb53c0a..28b2cc7c1b8 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -669,16 +669,19 @@ class ParticleData : public GeometryState { } // energy of last reaction secondaries const double& bank_second_E() const { return bank_second_E_; } - int& n_bank() { return n_bank_; } // number of banked fission sites - double& wgt_bank() { return wgt_bank_; } // weight of banked fission sites + // Fission neutrons recorded in nu_bank() for analog fission tallies: those + // stored in the fission or secondary bank and delayed neutrons that are not + // banked because they would be emitted after the time cutoff + int& n_bank() { return n_bank_; } // number of fission neutrons + double& wgt_bank() { return wgt_bank_; } // weight of fission neutrons int* n_delayed_bank() { return n_delayed_bank_; - } // number of delayed fission sites + } // number of delayed fission neutrons int& n_delayed_bank(int i) { return n_delayed_bank_[i]; - } // number of delayed fission sites + } // number of delayed fission neutrons // Index of cell particle is born in int& cell_born() { return cell_born_; } diff --git a/openmc/settings.py b/openmc/settings.py index 5d71df091c5..e152c65fe3b 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -84,7 +84,10 @@ class Settings: assigned to particles that are not killed after Russian roulette. Value of energy should be a float indicating energy in eV below which particle type will be killed. Value of time should be a float in seconds. - Particles will be killed exactly at the specified time. Value for + Particles will be killed exactly at the specified time. Delayed fission + neutrons that would be emitted after the 'time_neutron' cutoff are not + banked, but analog fission tallies count them as produced, as the + track-length and collision estimators do. Value for 'survival_normalization' is a bool indicating whether or not the weight cutoff parameters will be applied relative to the particle's starting weight or to its current weight. diff --git a/src/physics.cpp b/src/physics.cpp index 71dae94b590..e4852a97cfa 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -204,11 +204,21 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // or the secondary particle bank. bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); - // Counter for the number of fission sites successfully stored to the shared - // fission bank or the secondary particle bank - int n_sites_stored; + // Record a fission neutron for analog fission tallies + auto record_fission_neutron = [&](const SourceSite& site) { + // 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(); + nu_bank_entry.wgt = site.wgt; + nu_bank_entry.E = site.E; + nu_bank_entry.delayed_group = site.delayed_group; + }; - for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { + for (int i = 0; i < nu; ++i) { // Initialize fission site object with particle data SourceSite site; site.r = p.r(); @@ -220,10 +230,13 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // Sample delayed group and angle/energy for fission reaction sample_fission_neutron(i_nuclide, rx, &site, p); - // Reject site if it exceeds time cutoff + // Reject site if it exceeds time cutoff. The neutron would be killed at + // birth, so it is not banked, but it was produced by this fission and is + // recorded for analog fission tallies like any other fission neutron. if (site.delayed_group > 0) { double t_cutoff = settings::time_cutoff[site.particle.transport_index()]; if (site.time > t_cutoff) { + record_fission_neutron(site); continue; } } @@ -261,30 +274,22 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) p.n_secondaries()++; } - // 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(); - nu_bank_entry.wgt = site.wgt; - nu_bank_entry.E = site.E; - nu_bank_entry.delayed_group = site.delayed_group; + record_fission_neutron(site); } - // If shared fission bank was full, and no fissions could be added, - // set the particle fission flag to false. - if (n_sites_stored == 0) { + // If shared fission bank was full, and no fission neutrons could be + // recorded, set the particle fission flag to false. + if (p.nu_bank().empty()) { p.fission() = false; return; } - // Set nu to the number of fission sites successfully stored. If the fission - // bank was not found to be full then these values are already equivalent. - nu = n_sites_stored; + // Set nu to the number of fission neutrons recorded, i.e., the sites stored + // in a bank plus those rejected by the time cutoff. If the fission bank was + // not found to be full then these values are already equivalent. + nu = p.nu_bank().size(); - // Store the total weight banked for analog fission tallies + // Store the total weight recorded for analog fission tallies p.n_bank() = nu; p.wgt_bank() = nu / weight; for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 212ba765ca0..2280612f5bd 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -129,11 +129,21 @@ void create_fission_sites(Particle& p) // or the secondary particle bank. bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); - // Counter for the number of fission sites successfully stored to the shared - // fission bank or the secondary particle bank - int n_sites_stored; + // Record a fission neutron for analog fission tallies + auto record_fission_neutron = [&](const SourceSite& site) { + // 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(); + nu_bank_entry.wgt = site.wgt; + nu_bank_entry.E = site.E; + nu_bank_entry.delayed_group = site.delayed_group; + }; - for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { + for (int i = 0; i < nu; ++i) { // Initialize fission site object with particle data SourceSite site; site.r = p.r(); @@ -171,9 +181,12 @@ void create_fission_sites(Particle& p) macro_xs.get_xs(MgxsType::DECAY_RATE, 0, nullptr, nullptr, &dg, 0, 0); site.time -= std::log(prn(p.current_seed())) / decay_rate; - // Reject site if it exceeds time cutoff + // Reject site if it exceeds time cutoff. The neutron would be killed at + // birth, so it is not banked, but it was produced by this fission and + // is recorded for analog fission tallies like any other fission neutron. double t_cutoff = settings::time_cutoff[site.particle.transport_index()]; if (site.time > t_cutoff) { + record_fission_neutron(site); continue; } } @@ -210,30 +223,22 @@ void create_fission_sites(Particle& p) // Set the delayed group on the particle as well p.delayed_group() = dg + 1; - // Increment the number of neutrons born delayed - if (p.delayed_group() > 0) { - nu_d[dg]++; - } - - // Write fission particles to nuBank - NuBank& nu_bank_entry = p.nu_bank().emplace_back(); - nu_bank_entry.wgt = site.wgt; - nu_bank_entry.E = site.E; - nu_bank_entry.delayed_group = site.delayed_group; + record_fission_neutron(site); } - // If shared fission bank was full, and no fissions could be added, - // set the particle fission flag to false. - if (n_sites_stored == 0) { + // If shared fission bank was full, and no fission neutrons could be + // recorded, set the particle fission flag to false. + if (p.nu_bank().empty()) { p.fission() = false; return; } - // Set nu to the number of fission sites successfully stored. If the fission - // bank was not found to be full then these values are already equivalent. - nu = n_sites_stored; + // Set nu to the number of fission neutrons recorded, i.e., the sites stored + // in a bank plus those rejected by the time cutoff. If the fission bank was + // not found to be full then these values are already equivalent. + nu = p.nu_bank().size(); - // Store the total weight banked for analog fission tallies + // Store the total weight recorded for analog fission tallies p.n_bank() = nu; p.wgt_bank() = nu / weight; for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { diff --git a/tests/unit_tests/test_time_cutoff_fission.py b/tests/unit_tests/test_time_cutoff_fission.py new file mode 100644 index 00000000000..1939b51f656 --- /dev/null +++ b/tests/unit_tests/test_time_cutoff_fission.py @@ -0,0 +1,224 @@ +"""Analog fission tallies with a neutron time cutoff. + +A delayed fission neutron that would be emitted after the time cutoff is not +banked, since it would be killed at birth, but it was still produced by the +fission. Analog fission tallies must count it, as the track-length and +collision estimators do, and must keep the record of the fission neutrons of a +collision consistent with the number of neutrons they loop over. + +""" + +from pathlib import Path + +import numpy as np +import openmc +import pytest +from openmc.examples import slab_mg + + +T_CUTOFF = 1.0e-2 + + +def create_library(all_delayed=False): + """Two-group cross sections with two delayed groups. + + The delayed neutron fractions are unphysically large so that a short run + produces many delayed neutrons, nearly all of which are emitted after the + time cutoff. With all_delayed, every fission neutron is delayed. + + """ + groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6]) + fission = np.array([0.002817, 0.097]) + nu_fission = 2.5 * fission + capture = np.array([0.008708, 0.02518]) + scatter = np.array( + [[[0.31980, 0.06694], [0.004555, -0.0003972]], + [[0.00000, 0.00000], [0.424100, 0.05439000]]]) + if all_delayed: + beta = np.array([1/3, 2/3]) + prompt_nu_fission = np.zeros(2) + else: + beta = np.array([0.1, 0.2]) + prompt_nu_fission = (1.0 - beta.sum()) * nu_fission + + xsdata = openmc.XSdata('mat_1', groups, num_delayed_groups=2) + xsdata.order = 1 + xsdata.set_total([0.33588, 0.54628]) + xsdata.set_absorption(capture + fission) + xsdata.set_scatter_matrix(scatter) + xsdata.set_fission(fission) + xsdata.set_prompt_nu_fission(prompt_nu_fission) + xsdata.set_delayed_nu_fission(np.outer(beta, 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') + return Path('2g.h5').resolve() + + +def mg_model(all_delayed=False): + # The library path is absolute because the model may run in another + # directory + model = slab_mg(mgxslib_name=str(create_library(all_delayed))) + model.settings.seed = 1 + return model + + +def ce_model(): + model = openmc.Model() + material = openmc.Material() + material.add_nuclide('U235', 1.0) + material.set_density('g/cm3', 16.0) + sphere = openmc.Sphere(r=10.0, boundary_type='vacuum') + cell = openmc.Cell(region=-sphere, fill=material) + model.geometry = openmc.Geometry([cell]) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box(*cell.bounding_box), + constraints={'fissionable': True}) + model.settings.seed = 1 + return model + + +def add_tallies(model, energyout_edges, delayed_groups): + scores = ['nu-fission', 'prompt-nu-fission', 'delayed-nu-fission'] + for name, filters, tally_scores in [ + ('analog', [], scores), + ('energyout', [openmc.EnergyoutFilter(energyout_edges)], scores), + ('delayed', [openmc.DelayedGroupFilter(delayed_groups)], + ['delayed-nu-fission', 'decay-rate']), + ]: + tally = openmc.Tally(name=name) + tally.filters = filters + tally.scores = tally_scores + tally.estimator = 'analog' + model.tallies.append(tally) + tally = openmc.Tally(name='flux') + tally.scores = ['flux'] + tally.estimator = 'tracklength' + model.tallies.append(tally) + + +def get_mean(sp, name, score): + tally = sp.get_tally(name=name) + return tally.get_values(scores=[score]).ravel() + + +def assert_equal(a, b): + """Equal to round-off in the sums of the tally results.""" + a = np.asarray(a) + b = np.asarray(b) + assert np.abs(a - b).max() <= 1e-10 * np.abs(b).max() + + +@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +@pytest.mark.parametrize('energy_mode', ['mg', 'ce']) +def test_time_cutoff_eigenvalue(run_in_tmpdir, energy_mode, event_based): + """The first generation with and without a time cutoff. + + The source particles are transported the same way whether or not the + cutoff is set, because no history lasts as long as the cutoff, so analog + tallies of the fission neutrons produced must be the same. Only the + fission bank differs: with the cutoff, delayed neutrons emitted after it + are not banked. + + """ + if energy_mode == 'mg': + model = mg_model() + model.settings.particles = 1000 + energyout_edges = [0.0, 0.625, 20.0e6] + delayed_groups = [1, 2] + else: + model = ce_model() + model.settings.particles = 10000 + # Fission neutron energies are resampled below the maximum energy of + # the data, so these bins cover every fission neutron + energyout_edges = [0.0, 1.0e6, 1.0e8] + delayed_groups = [1, 2, 3, 4, 5, 6] + model.settings.batches = 1 + model.settings.inactive = 0 + model.settings.event_based = event_based + add_tallies(model, energyout_edges, delayed_groups) + + results = {} + for run in ('no-cutoff', 'cutoff'): + if run == 'cutoff': + model.settings.cutoff = {'time_neutron': T_CUTOFF} + sp_path = model.run(cwd=run) + with openmc.StatePoint(sp_path) as sp: + res = {'k': sp.k_generation, 'source': sp.source} + for name in ('flux', 'analog', 'energyout', 'delayed'): + for score in sp.get_tally(name=name).scores: + res[name, score] = get_mean(sp, name, score) + results[run] = res + a = results['cutoff'] + b = results['no-cutoff'] + + # The transport is the same in both runs + assert_equal(a['k'], b['k']) + assert_equal(a['flux', 'flux'], b['flux', 'flux']) + + # Delayed neutrons emitted after the cutoff are produced in both runs, but + # banked only without the cutoff + source = b['source'] + assert np.any((source['delayed_group'] > 0) & (source['time'] > T_CUTOFF)) + assert np.all(a['source']['time'] <= T_CUTOFF) + + for key in a: + if key[0] in ('analog', 'energyout', 'delayed'): + assert_equal(a[key], b[key]) + + for res in (a, b): + # The tally with an EnergyoutFilter loops over the fission neutrons of + # each collision, so it must sum to the tally without the filter + for score in ('nu-fission', 'prompt-nu-fission', 'delayed-nu-fission'): + assert_equal(res['energyout', score].sum(), res['analog', score]) + + if energy_mode == 'mg': + # The decay rate loops over the fission neutrons of each collision, + # and the delayed-nu-fission score counts them + decay_rate = np.array([0.0127, 0.0317]) + assert_equal(res['delayed', 'decay-rate'], + decay_rate * res['delayed', 'delayed-nu-fission']) + + if energy_mode == 'mg': + assert np.all(a['delayed', 'delayed-nu-fission'] > 0.0) + else: + assert a['analog', 'delayed-nu-fission'].sum() > 0.0 + + +@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +def test_time_cutoff_all_delayed(run_in_tmpdir, event_based): + """No fission neutron is prompt, so no prompt-nu-fission may be scored.""" + model = mg_model(all_delayed=True) + model.settings.run_mode = 'fixed source' + model.settings.batches = 5 + model.settings.particles = 1000 + model.settings.event_based = event_based + model.settings.cutoff = {'time_neutron': T_CUTOFF} + scores = ['nu-fission', 'prompt-nu-fission', 'delayed-nu-fission'] + for name, filters in [ + ('analog', []), + ('energyout', [openmc.EnergyoutFilter([0.0, 0.625, 20.0e6])]), + ]: + tally = openmc.Tally(name=name) + tally.filters = filters + tally.scores = scores + tally.estimator = 'analog' + model.tallies.append(tally) + + sp_path = model.run() + with openmc.StatePoint(sp_path) as sp: + for name in ('analog', 'energyout'): + nu_fission = get_mean(sp, name, 'nu-fission') + prompt = get_mean(sp, name, 'prompt-nu-fission') + delayed = get_mean(sp, name, 'delayed-nu-fission') + assert np.all(prompt == 0.0) + assert_equal(delayed, nu_fission) + assert np.all(delayed > 0.0) + for score in scores: + assert_equal(get_mean(sp, 'energyout', score).sum(), + get_mean(sp, 'analog', score)) From 09fefd0f0bcbd95d654e6a946e58eab9bde76eb4 Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 7 Oct 2026 05:24:30 -0400 Subject: [PATCH 2/4] Test time cutoff fission tallies with full banks, UFS and fixed source Extend tests/unit_tests/test_time_cutoff_fission.py with checks that do not depend on comparing two runs: - With a multigroup library in which every collision produces exactly one delayed fission neutron, emitted long after the time cutoff, the analog nu-fission and delayed-nu-fission estimates must equal the collision estimates collision by collision. No fission neutron may enter the simulation (ParticleProductionFilter, track-length flux), and collision tracking must show that the colliding particle keeps the delayed group of its source site and takes no progeny ID. - In a continuous-energy fixed-source problem, the flux must be positive and the analog nu-fission, prompt-nu-fission and delayed-nu-fission estimates must agree with the track-length estimates. - When the fission bank fills up in the first generation (multigroup, eight neutrons per collision; continuous energy, Cf-249), the analog tallies must count exactly the sites in the bank, plus, in the multigroup run with a time cutoff, the delayed neutrons it rejected. - With uniform fission site weighting, the analog tallies with an EnergyoutFilter must sum to the tallies without it, in multigroup and continuous-energy problems. Also set inactive = 0 explicitly in the all-delayed fixed-source test. --- tests/unit_tests/test_time_cutoff_fission.py | 275 +++++++++++++++++++ 1 file changed, 275 insertions(+) diff --git a/tests/unit_tests/test_time_cutoff_fission.py b/tests/unit_tests/test_time_cutoff_fission.py index 1939b51f656..90f62013466 100644 --- a/tests/unit_tests/test_time_cutoff_fission.py +++ b/tests/unit_tests/test_time_cutoff_fission.py @@ -196,6 +196,7 @@ def test_time_cutoff_all_delayed(run_in_tmpdir, event_based): model = mg_model(all_delayed=True) model.settings.run_mode = 'fixed source' model.settings.batches = 5 + model.settings.inactive = 0 model.settings.particles = 1000 model.settings.event_based = event_based model.settings.cutoff = {'time_neutron': T_CUTOFF} @@ -222,3 +223,277 @@ def test_time_cutoff_all_delayed(run_in_tmpdir, event_based): for score in scores: assert_equal(get_mean(sp, 'energyout', score).sum(), get_mean(sp, 'analog', score)) + + +# Decay rates so small that every delayed neutron is emitted long after the +# time cutoff +SLOW_DECAY_RATE = np.array([1.0e-6, 2.0e-6]) + + +def create_exact_library(nu, prompt_fraction): + """Two-group cross sections in which a collision produces nu neutrons. + + nu-fission / total equals nu in both groups and every number is a binary + fraction, so a collision of a particle of unit weight samples exactly nu + fission neutrons while keff is 1, i.e., in a fixed-source run and in the + first generation of an eigenvalue run. + + """ + groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6]) + nu_fission = np.full(2, 0.5 * nu) + xsdata = openmc.XSdata('mat_1', groups, num_delayed_groups=2) + xsdata.order = 0 + xsdata.set_total([0.5, 0.5]) + xsdata.set_absorption([0.375, 0.375]) + xsdata.set_scatter_matrix([[[0.0625], [0.0625]], [[0.0], [0.125]]]) + xsdata.set_fission([0.25, 0.25]) + xsdata.set_prompt_nu_fission(prompt_fraction * nu_fission) + xsdata.set_delayed_nu_fission( + np.outer([0.5, 0.5], (1.0 - prompt_fraction) * nu_fission)) + xsdata.set_chi_prompt([1.0, 0.0]) + xsdata.set_chi_delayed([[0.75, 0.25], [0.5, 0.5]]) + xsdata.set_decay_rate(SLOW_DECAY_RATE) + 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('exact.h5') + return Path('exact.h5').resolve() + + +@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +def test_time_cutoff_one_neutron_per_collision(run_in_tmpdir, event_based): + """Every collision produces one delayed neutron, which is never banked. + + Each collision of a fixed-source run produces exactly one fission neutron, + which is delayed and emitted long after the time cutoff. The analog + estimates must then equal the collision estimates collision by collision, + and no fission neutron may enter the simulation or change the particle + that produced it. + + """ + model = slab_mg(mgxslib_name=str(create_exact_library(1, 0.0))) + model.settings.run_mode = 'fixed source' + model.settings.batches = 2 + model.settings.inactive = 0 + model.settings.particles = 1000 + model.settings.seed = 1 + model.settings.event_based = event_based + model.settings.cutoff = {'time_neutron': T_CUTOFF} + model.settings.collision_track = {'max_collisions': 100000} + + scores = ['nu-fission', 'prompt-nu-fission', 'delayed-nu-fission'] + for name, filters, tally_scores, estimator in [ + ('analog', [], scores, 'analog'), + ('collision', [], ['nu-fission', 'delayed-nu-fission'], 'collision'), + ('energyout', [openmc.EnergyoutFilter([0.0, 0.625, 20.0e6])], scores, + 'analog'), + ('delayed', [openmc.DelayedGroupFilter([1, 2])], + ['delayed-nu-fission', 'decay-rate'], 'analog'), + ('flux', [], ['flux'], 'tracklength'), + ('production', [openmc.ParticleProductionFilter('neutron')], + ['events'], 'analog'), + ]: + tally = openmc.Tally(name=name) + tally.filters = filters + tally.scores = tally_scores + tally.estimator = estimator + model.tallies.append(tally) + + sp_path = model.run() + with openmc.StatePoint(sp_path) as sp: + nu_fission = get_mean(sp, 'analog', 'nu-fission') + assert np.all(nu_fission > 0.0) + + # One neutron per collision, all of them delayed + assert_equal(nu_fission, get_mean(sp, 'collision', 'nu-fission')) + assert_equal(get_mean(sp, 'collision', 'delayed-nu-fission'), + nu_fission) + assert_equal(get_mean(sp, 'analog', 'delayed-nu-fission'), nu_fission) + assert np.all(get_mean(sp, 'analog', 'prompt-nu-fission') == 0.0) + assert np.all(get_mean(sp, 'energyout', 'prompt-nu-fission') == 0.0) + for score in scores: + assert_equal(get_mean(sp, 'energyout', score).sum(), + get_mean(sp, 'analog', score)) + delayed = get_mean(sp, 'delayed', 'delayed-nu-fission') + assert_equal(delayed.sum(), nu_fission) + assert_equal(get_mean(sp, 'delayed', 'decay-rate'), + SLOW_DECAY_RATE * delayed) + + # No fission neutron is transported + assert np.all(get_mean(sp, 'production', 'events') == 0.0) + assert np.all(get_mean(sp, 'flux', 'flux') > 0.0) + + # The particle that produced the fission neutrons is unchanged by them: it + # keeps the delayed group of its source site and takes no progeny + collisions = openmc.read_collision_track_hdf5('collision_track.h5') + assert collisions.size > 0 + assert np.all(collisions['delayed_group'] == 0) + assert np.all(collisions['progeny_id'] == 0) + + +@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +def test_time_cutoff_fixed_source_ce(run_in_tmpdir, event_based): + """Analog and track-length estimates agree in a subcritical sphere. + + Only fission neutrons emitted before the cutoff may enter the simulation: + one emitted after it would be moved back to the cutoff and score a + negative track length, so the flux must be positive. The analog estimates + of the fission neutrons produced must agree with the track-length + estimates, which count every neutron produced whatever the cutoff. + + """ + model = openmc.Model() + material = openmc.Material() + material.add_nuclide('U235', 1.0) + material.set_density('g/cm3', 16.0) + sphere = openmc.Sphere(r=8.0, boundary_type='vacuum') + cell = openmc.Cell(region=-sphere, fill=material) + model.geometry = openmc.Geometry([cell]) + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point(), energy=openmc.stats.delta_function(2.0e6)) + model.settings.batches = 10 + model.settings.particles = 1000 + model.settings.seed = 1 + model.settings.event_based = event_based + model.settings.cutoff = {'time_neutron': T_CUTOFF} + + scores = ['nu-fission', 'prompt-nu-fission', 'delayed-nu-fission'] + for name, filters, tally_scores, estimator in [ + ('analog', [], scores, 'analog'), + ('tracklength', [], scores + ['flux'], 'tracklength'), + ('energyout', [openmc.EnergyoutFilter([0.0, 1.0e6, 1.0e8])], scores, + 'analog'), + ]: + tally = openmc.Tally(name=name) + tally.filters = filters + tally.scores = tally_scores + tally.estimator = estimator + model.tallies.append(tally) + + sp_path = model.run() + with openmc.StatePoint(sp_path) as sp: + assert np.all(get_mean(sp, 'tracklength', 'flux') > 0.0) + for score in scores: + analog = get_mean(sp, 'analog', score) + tracklength = get_mean(sp, 'tracklength', score) + # The relative standard deviation of the analog estimate is below + # 6% for delayed neutrons and below 3% otherwise + assert np.all(np.abs(analog / tracklength - 1.0) < 0.5) + assert_equal(get_mean(sp, 'energyout', score).sum(), analog) + + +@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +@pytest.mark.parametrize('energy_mode', ['mg', 'ce']) +def test_time_cutoff_full_fission_bank(run_in_tmpdir, energy_mode, + event_based): + """The fission bank fills up in the first generation. + + The fission bank holds three sites per source particle. The site that + fails to enter the full bank and the later sites of its collision are + neither banked nor counted, so analog tallies count exactly the sites in + the bank, plus the delayed neutrons that the time cutoff rejected before + that site. With the multigroup library and the cutoff, every delayed + neutron is rejected, so the bank holds prompt neutrons only. + + """ + if energy_mode == 'mg': + # Every collision produces eight neutrons and the bank size, 1500, is + # not a multiple of eight, so the bank fills in the middle of a + # collision + model = slab_mg(mgxslib_name=str(create_exact_library(8, 0.5))) + runs = ('no-cutoff', 'cutoff') + else: + # k-infinity is about 4, so the fission sites of the first generation + # do not fit in the bank + model = openmc.Model() + material = openmc.Material() + material.add_nuclide('Cf249', 1.0) + material.set_density('g/cm3', 15.0) + box = openmc.model.RectangularParallelepiped( + -5.0, 5.0, -5.0, 5.0, -5.0, 5.0, boundary_type='reflective') + cell = openmc.Cell(region=-box, fill=material) + model.geometry = openmc.Geometry([cell]) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box(*cell.bounding_box)) + # A delayed neutron of the evaluated data can be emitted before the + # cutoff and banked, so the count is exact only without a cutoff + runs = ('no-cutoff',) + model.settings.batches = 1 + model.settings.inactive = 0 + model.settings.particles = 500 + model.settings.seed = 1 + model.settings.event_based = event_based + scores = ['nu-fission', 'prompt-nu-fission', 'delayed-nu-fission'] + tally = openmc.Tally(name='analog') + tally.scores = scores + tally.estimator = 'analog' + model.tallies.append(tally) + + for run in runs: + if run == 'cutoff': + model.settings.cutoff = {'time_neutron': T_CUTOFF} + sp_path = model.run(cwd=run) + with openmc.StatePoint(sp_path) as sp: + nu_fission, prompt, delayed = ( + get_mean(sp, 'analog', score) for score in scores) + assert_equal(prompt + delayed, nu_fission) + if run == 'no-cutoff': + # The tallies are normalized by the number of source particles + assert_equal(nu_fission, [3.0]) + else: + assert_equal(prompt, [3.0]) + assert np.all(delayed > 0.0) + + +@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +@pytest.mark.parametrize('energy_mode', ['mg', 'ce']) +def test_time_cutoff_ufs(run_in_tmpdir, energy_mode, event_based): + """Fission neutrons whose weight is set by uniform fission site weighting. + + From the second generation on, a fission neutron has the weight 1/w of + its collision's UFS weight w, so the analog tallies that loop over the + fission neutrons of a collision must sum to the tallies that use the + total weight of the collision's fission neutrons. + + """ + mesh = openmc.RegularMesh() + if energy_mode == 'mg': + model = mg_model() + model.settings.particles = 1000 + mesh.lower_left = (0.0, -1.0e6, -1.0e6) + mesh.upper_right = (929.45, 1.0e6, 1.0e6) + mesh.dimension = (4, 1, 1) + energyout_edges = [0.0, 0.625, 20.0e6] + delayed_groups = [1, 2] + else: + model = ce_model() + model.settings.particles = 2000 + mesh.lower_left = (-10.0, -10.0, -10.0) + mesh.upper_right = (10.0, 10.0, 10.0) + mesh.dimension = (3, 3, 3) + energyout_edges = [0.0, 1.0e6, 1.0e8] + delayed_groups = [1, 2, 3, 4, 5, 6] + model.settings.batches = 3 + model.settings.inactive = 1 + model.settings.event_based = event_based + model.settings.cutoff = {'time_neutron': T_CUTOFF} + model.settings.ufs_mesh = mesh + add_tallies(model, energyout_edges, delayed_groups) + + sp_path = model.run() + with openmc.StatePoint(sp_path) as sp: + # The weights are not all one + assert np.any(sp.source['wgt'] != 1.0) + for score in ('nu-fission', 'prompt-nu-fission', 'delayed-nu-fission'): + assert_equal(get_mean(sp, 'energyout', score).sum(), + get_mean(sp, 'analog', score)) + delayed = get_mean(sp, 'delayed', 'delayed-nu-fission') + assert delayed.sum() > 0.0 + assert_equal(delayed.sum(), + get_mean(sp, 'analog', 'delayed-nu-fission')) + if energy_mode == 'mg': + decay_rate = np.array([0.0127, 0.0317]) + assert_equal(get_mean(sp, 'delayed', 'decay-rate'), + decay_rate * delayed) From f92ed975bb4ac29dc4e759f31de9b7b88155c6d0 Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 7 Oct 2026 11:03:38 -0400 Subject: [PATCH 3/4] Wrap long lines in the time cutoff tests Wrap the six @pytest.mark.parametrize('event_based', ...) lines of tests/unit_tests/test_time_cutoff_fission.py that exceed 79 characters. In the multigroup full fission bank test with a time cutoff, also check that the delayed neutrons sampled after the site that failed to enter the bank are not counted: the analog delayed-nu-fission estimate must stay well below the collision estimate, which counts all of them. --- tests/unit_tests/test_time_cutoff_fission.py | 32 ++++++++++++++------ 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/tests/unit_tests/test_time_cutoff_fission.py b/tests/unit_tests/test_time_cutoff_fission.py index 90f62013466..d8b21eee460 100644 --- a/tests/unit_tests/test_time_cutoff_fission.py +++ b/tests/unit_tests/test_time_cutoff_fission.py @@ -114,7 +114,8 @@ def assert_equal(a, b): assert np.abs(a - b).max() <= 1e-10 * np.abs(b).max() -@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +@pytest.mark.parametrize('event_based', [False, True], + ids=['history', 'event']) @pytest.mark.parametrize('energy_mode', ['mg', 'ce']) def test_time_cutoff_eigenvalue(run_in_tmpdir, energy_mode, event_based): """The first generation with and without a time cutoff. @@ -190,7 +191,8 @@ def test_time_cutoff_eigenvalue(run_in_tmpdir, energy_mode, event_based): assert a['analog', 'delayed-nu-fission'].sum() > 0.0 -@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +@pytest.mark.parametrize('event_based', [False, True], + ids=['history', 'event']) def test_time_cutoff_all_delayed(run_in_tmpdir, event_based): """No fission neutron is prompt, so no prompt-nu-fission may be scored.""" model = mg_model(all_delayed=True) @@ -261,7 +263,8 @@ def create_exact_library(nu, prompt_fraction): return Path('exact.h5').resolve() -@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +@pytest.mark.parametrize('event_based', [False, True], + ids=['history', 'event']) def test_time_cutoff_one_neutron_per_collision(run_in_tmpdir, event_based): """Every collision produces one delayed neutron, which is never banked. @@ -332,7 +335,8 @@ def test_time_cutoff_one_neutron_per_collision(run_in_tmpdir, event_based): assert np.all(collisions['progeny_id'] == 0) -@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +@pytest.mark.parametrize('event_based', [False, True], + ids=['history', 'event']) def test_time_cutoff_fixed_source_ce(run_in_tmpdir, event_based): """Analog and track-length estimates agree in a subcritical sphere. @@ -384,7 +388,8 @@ def test_time_cutoff_fixed_source_ce(run_in_tmpdir, event_based): assert_equal(get_mean(sp, 'energyout', score).sum(), analog) -@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +@pytest.mark.parametrize('event_based', [False, True], + ids=['history', 'event']) @pytest.mark.parametrize('energy_mode', ['mg', 'ce']) def test_time_cutoff_full_fission_bank(run_in_tmpdir, energy_mode, event_based): @@ -426,10 +431,11 @@ def test_time_cutoff_full_fission_bank(run_in_tmpdir, energy_mode, model.settings.seed = 1 model.settings.event_based = event_based scores = ['nu-fission', 'prompt-nu-fission', 'delayed-nu-fission'] - tally = openmc.Tally(name='analog') - tally.scores = scores - tally.estimator = 'analog' - model.tallies.append(tally) + for estimator in ('analog', 'collision'): + tally = openmc.Tally(name=estimator) + tally.scores = scores + tally.estimator = estimator + model.tallies.append(tally) for run in runs: if run == 'cutoff': @@ -438,6 +444,7 @@ def test_time_cutoff_full_fission_bank(run_in_tmpdir, energy_mode, with openmc.StatePoint(sp_path) as sp: nu_fission, prompt, delayed = ( get_mean(sp, 'analog', score) for score in scores) + collision_delayed = get_mean(sp, 'collision', 'delayed-nu-fission') assert_equal(prompt + delayed, nu_fission) if run == 'no-cutoff': # The tallies are normalized by the number of source particles @@ -445,9 +452,14 @@ def test_time_cutoff_full_fission_bank(run_in_tmpdir, energy_mode, else: assert_equal(prompt, [3.0]) assert np.all(delayed > 0.0) + # The delayed neutrons sampled after the site that failed to enter + # the bank are not counted either, so the analog estimate stays + # well below the collision estimate, which counts all of them + assert np.all(delayed < 0.85 * collision_delayed) -@pytest.mark.parametrize('event_based', [False, True], ids=['history', 'event']) +@pytest.mark.parametrize('event_based', [False, True], + ids=['history', 'event']) @pytest.mark.parametrize('energy_mode', ['mg', 'ce']) def test_time_cutoff_ufs(run_in_tmpdir, energy_mode, event_based): """Fission neutrons whose weight is set by uniform fission site weighting. From da4eae03e47829ad22a29a283948edb002df6e46 Mon Sep 17 00:00:00 2001 From: Ondrej Chvala Date: Wed, 7 Oct 2026 11:26:05 -0400 Subject: [PATCH 4/4] Reword a comment in the time cutoff full bank test Once a site fails to enter the full fission bank, the loop over the fission sites of the collision ends, so the neutrons of the remaining sites are never sampled. Say so in the comment of the multigroup check. --- tests/unit_tests/test_time_cutoff_fission.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/unit_tests/test_time_cutoff_fission.py b/tests/unit_tests/test_time_cutoff_fission.py index d8b21eee460..bf5c5840379 100644 --- a/tests/unit_tests/test_time_cutoff_fission.py +++ b/tests/unit_tests/test_time_cutoff_fission.py @@ -452,9 +452,10 @@ def test_time_cutoff_full_fission_bank(run_in_tmpdir, energy_mode, else: assert_equal(prompt, [3.0]) assert np.all(delayed > 0.0) - # The delayed neutrons sampled after the site that failed to enter - # the bank are not counted either, so the analog estimate stays - # well below the collision estimate, which counts all of them + # The delayed neutrons that the sites after the one that failed to + # enter the bank would have produced are neither sampled nor + # counted, so the analog estimate stays well below the collision + # estimate, which counts all of them assert np.all(delayed < 0.85 * collision_delayed)