From abb3e92c2f6954976fd652296321d46d7102603c Mon Sep 17 00:00:00 2001 From: Matthew Louis Date: Wed, 17 Dec 2025 17:03:08 -0500 Subject: [PATCH 1/8] Added new input option for calculating subcritical k in fixed source mode --- include/openmc/settings.h | 2 + openmc/settings.py | 29 +++++++++ src/finalize.cpp | 1 + src/settings.cpp | 11 ++++ .../fixed_source_k/inputs_true.dat | 35 +++++++++++ .../fixed_source_k/results_true.dat | 0 tests/regression_tests/fixed_source_k/test.py | 63 +++++++++++++++++++ 7 files changed, 141 insertions(+) create mode 100644 tests/regression_tests/fixed_source_k/inputs_true.dat create mode 100644 tests/regression_tests/fixed_source_k/results_true.dat create mode 100644 tests/regression_tests/fixed_source_k/test.py diff --git a/include/openmc/settings.h b/include/openmc/settings.h index b369c99fef8..1981aa837c9 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -169,6 +169,8 @@ extern double source_rejection_fraction; //!< Minimum fraction of source sites //!< that must be accepted extern double free_gas_threshold; //!< Threshold multiplier for free gas //!< scattering treatment +extern bool calculate_subcritical_k; //!< Calculate subcritical k in fixed + //!< source mode extern int max_history_splits; //!< maximum number of particle splits for weight windows diff --git a/openmc/settings.py b/openmc/settings.py index 76e191c5e02..848a1d8e2ac 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -365,6 +365,11 @@ class Settings: .. versionadded::0.14.0 write_initial_source : bool Indicate whether to write the initial source distribution to file + + .. versionadded::0.15.4 + calculate_subcritical_k : bool + Indicate whether to calculate and report the subcritical multiplication + factor k in fixed source simulations """ def __init__(self, **kwargs): @@ -377,6 +382,7 @@ def __init__(self, **kwargs): self._max_write_lost_particles = None self._particles = None self._keff_trigger = None + self._calculate_subcritical_k = False # Energy mode subelement self._energy_mode = None @@ -1385,6 +1391,18 @@ def free_gas_threshold(self, free_gas_threshold: float | None): cv.check_greater_than('free gas threshold', free_gas_threshold, 0.0) self._free_gas_threshold = free_gas_threshold + @property + def calculate_subcritical_k(self) -> bool: + return self._calculate_subcritical_k + + @calculate_subcritical_k.setter + def calculate_subcritical_k(self, calculate_subcritical_k: bool): + cv.check_type('calculate subcritical k', calculate_subcritical_k, bool) + if not self._run_mode == RunMode.FIXED_SOURCE: + raise ValueError("calculate_subcritical_k can only be set when " + "run_mode is 'fixed source'") + self._calculate_subcritical_k = calculate_subcritical_k + def _create_run_mode_subelement(self, root): elem = ET.SubElement(root, "run_mode") elem.text = self._run_mode.value @@ -1913,6 +1931,11 @@ def _create_free_gas_threshold_subelement(self, root): element = ET.SubElement(root, "free_gas_threshold") element.text = str(self._free_gas_threshold) + def _create_calculate_subcritical_k_subelement(self, root): + if self._calculate_subcritical_k: + elem = ET.SubElement(root, "calculate_subcritical_k") + elem.text = str(self._calculate_subcritical_k).lower() + def _eigenvalue_from_xml_element(self, root): elem = root.find('eigenvalue') if elem is not None: @@ -2376,6 +2399,11 @@ def _free_gas_threshold_from_xml_element(self, root): if text is not None: self.free_gas_threshold = float(text) + def _calculate_subcritical_k_from_xml_element(self, root): + text = get_text(root, 'calculate_subcritical_k') + if text is not None: + self.calculate_subcritical_k = text in ('true', '1') + def to_xml_element(self, mesh_memo=None): """Create a 'settings' element to be written to an XML file. @@ -2448,6 +2476,7 @@ def to_xml_element(self, mesh_memo=None): self._create_use_decay_photons_subelement(element) self._create_source_rejection_fraction_subelement(element) self._create_free_gas_threshold_subelement(element) + self._create_calculate_subcritical_k_subelement(element) # Clean the indentation in the file to be user-readable clean_indentation(element) diff --git a/src/finalize.cpp b/src/finalize.cpp index 344eaa1a0a7..02ca5a8d5c9 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -88,6 +88,7 @@ int openmc_finalize() settings::entropy_on = false; settings::event_based = false; settings::free_gas_threshold = 400.0; + settings::calculate_subcritical_k = false; settings::gen_per_batch = 1; settings::legendre_to_tabular = true; settings::legendre_to_tabular_points = -1; diff --git a/src/settings.cpp b/src/settings.cpp index 5b472468fc4..3336e2acb2b 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -130,6 +130,7 @@ std::unordered_set sourcepoint_batch; std::unordered_set statepoint_batch; double source_rejection_fraction {0.05}; double free_gas_threshold {400.0}; +bool calculate_subcritical_k {false}; std::unordered_set source_write_surf_id; CollisionTrackConfig collision_track_config {}; int64_t ssw_max_particles; @@ -676,6 +677,16 @@ void read_settings_xml(pugi::xml_node root) free_gas_threshold = std::stod(get_node_value(root, "free_gas_threshold")); } + if (check_for_node(root, "calculate_subcritical_k")) { + if (run_mode == RunMode::FIXED_SOURCE) { + calculate_subcritical_k = + get_node_value_bool(root, "calculate_subcritical_k"); + } else { + fatal_error("The 'calculate_subcritical_k' setting is only valid in " + "fixed source mode."); + } + } + // Survival biasing if (check_for_node(root, "survival_biasing")) { survival_biasing = get_node_value_bool(root, "survival_biasing"); diff --git a/tests/regression_tests/fixed_source_k/inputs_true.dat b/tests/regression_tests/fixed_source_k/inputs_true.dat new file mode 100644 index 00000000000..af662d47749 --- /dev/null +++ b/tests/regression_tests/fixed_source_k/inputs_true.dat @@ -0,0 +1,35 @@ + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + 5 + + + 0 0 0 + + + + + + diff --git a/tests/regression_tests/fixed_source_k/results_true.dat b/tests/regression_tests/fixed_source_k/results_true.dat new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/fixed_source_k/test.py b/tests/regression_tests/fixed_source_k/test.py new file mode 100644 index 00000000000..370ca21297d --- /dev/null +++ b/tests/regression_tests/fixed_source_k/test.py @@ -0,0 +1,63 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + +def create_universe(): + # Define materials + heu = openmc.Material(name='HEU') + heu.add_element('U', 1, enrichment=93.0, enrichment_type='wo') + heu.set_density('g/cm3', 19.1) + + dep_uranium = openmc.Material(name='Depleted Uranium') + dep_uranium.add_nuclide('U238', 1.0) + dep_uranium.set_density('g/cm3', 19.1) + + mats = openmc.Materials([heu, dep_uranium]) + mats.export_to_xml() + + # Geometry + fuel_radius = 4.8 + dep_uranium_radius = 12.2 + + fuel_sphere = openmc.Sphere(r=fuel_radius) + dep_uranium_sphere = openmc.Sphere(r=dep_uranium_radius, boundary_type='vacuum') + + fuel_region = -fuel_sphere + dep_uranium_region = +fuel_sphere & -dep_uranium_sphere + + fuel_cell = openmc.Cell(name='Fuel', region=fuel_region) + fuel_cell.fill = heu + + dep_uranium_cell = openmc.Cell(name='Depleted Uranium', region=dep_uranium_region) + dep_uranium_cell.fill = dep_uranium + + universe = openmc.Universe(cells=[fuel_cell, dep_uranium_cell]) + return universe + +def test_source(): + source_space = openmc.stats.Point((0,0,0)) + source_angle = openmc.stats.Isotropic() + source_energy = openmc.stats.Maxwell(293.6) + source = openmc.IndependentSource(space=source_space, angle=source_angle, energy=source_energy) + return source + +@pytest.fixture +def model(): + model = openmc.Model() + + universe = create_universe() + model.geometry = openmc.Geometry(universe) + model.settings.run_mode = 'fixed source' + model.settings.source = test_source() + + model.settings.batches = 10 + model.settings.inactive = 5 + model.settings.particles = 1000 + + return model + +def test_fixed_source_run(model): + """Test that fixed source run with subcritical k calculation works.""" + harness = PyAPITestHarness("statepoint.10.h5", model, inputs_true='inputs_true.dat') + harness.main() \ No newline at end of file From 3cd6b2753e73378a767669edb98c70de6e828bf3 Mon Sep 17 00:00:00 2001 From: Matthew Louis Date: Thu, 18 Dec 2025 09:32:05 -0500 Subject: [PATCH 2/8] Implemented k_sq calculation and printing. Composite and combined estimators in output are wrong --- src/output.cpp | 35 +++++--- src/particle.cpp | 8 +- src/physics.cpp | 8 +- src/settings.cpp | 4 + src/simulation.cpp | 86 +++++++++++++++++-- .../fixed_source_k/inputs_true.dat | 1 + tests/regression_tests/fixed_source_k/test.py | 1 + 7 files changed, 119 insertions(+), 24 deletions(-) diff --git a/src/output.cpp b/src/output.cpp index 80e2b10ab89..06b4e42aca0 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -396,7 +396,7 @@ void print_generation() fmt::print(" {:>9} {:8.5f}", batch_and_gen, simulation::k_generation[idx]); // write out entropy info - if (settings::entropy_on) { + if (settings::entropy_on && !settings::calculate_subcritical_k) { fmt::print(" {:8.5f}", simulation::entropy[idx]); } @@ -534,21 +534,28 @@ void print_results() // write global tallies const auto& gt = simulation::global_tallies; double mean, stdev; + std::string eigenvalue_name = " k-effective"; + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + eigenvalue_name = " Subcritical k"; + } if (n > 1) { - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { std::tie(mean, stdev) = mean_stdev(>(GlobalTally::K_COLLISION, 0), n); - fmt::print(" k-effective (Collision) = {:.5f} +/- {:.5f}\n", mean, - t_n1 * stdev); + fmt::print(" {} (Collision) = {:.5f} +/- {:.5f}\n", eigenvalue_name, + mean, t_n1 * stdev); std::tie(mean, stdev) = mean_stdev(>(GlobalTally::K_TRACKLENGTH, 0), n); - fmt::print(" k-effective (Track-length) = {:.5f} +/- {:.5f}\n", mean, - t_n1 * stdev); + fmt::print(" {} (Track-length) = {:.5f} +/- {:.5f}\n", eigenvalue_name, + mean, t_n1 * stdev); std::tie(mean, stdev) = mean_stdev(>(GlobalTally::K_ABSORPTION, 0), n); - fmt::print(" k-effective (Absorption) = {:.5f} +/- {:.5f}\n", mean, - t_n1 * stdev); + fmt::print(" {} (Absorption) = {:.5f} +/- {:.5f}\n", eigenvalue_name, + mean, t_n1 * stdev); if (n > 3) { double k_combined[2]; openmc_get_keff(k_combined); - fmt::print(" Combined k-effective = {:.5f} +/- {:.5f}\n", + fmt::print(" Combined {} = {:.5f} +/- {:.5f}\n", eigenvalue_name, k_combined[0], k_combined[1]); } } @@ -560,12 +567,14 @@ void print_results() warning("Could not compute uncertainties -- only one " "active batch simulated!"); - if (settings::run_mode == RunMode::EIGENVALUE) { - fmt::print(" k-effective (Collision) = {:.5f}\n", + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { + fmt::print(" {} (Collision) = {:.5f}\n", eigenvalue_name, gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n); - fmt::print(" k-effective (Track-length) = {:.5f}\n", + fmt::print(" {} (Track-length) = {:.5f}\n", eigenvalue_name, gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n); - fmt::print(" k-effective (Absorption) = {:.5f}\n", + fmt::print(" {} (Absorption) = {:.5f}\n", eigenvalue_name, gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n); } fmt::print(" Leakage Fraction = {:.5f}\n", diff --git a/src/particle.cpp b/src/particle.cpp index 2d70d715e22..38a338b0283 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -269,7 +269,9 @@ void Particle::event_advance() } // Score track-length estimate of k-eff - if (settings::run_mode == RunMode::EIGENVALUE && + if ((settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) && type() == ParticleType::neutron) { keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission; } @@ -331,7 +333,9 @@ void Particle::event_cross_surface() void Particle::event_collide() { // Score collision estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE && + if ((settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) && type() == ParticleType::neutron) { keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total; } diff --git a/src/physics.cpp b/src/physics.cpp index 41509af97be..3e3bd53cd09 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -646,7 +646,9 @@ void absorption(Particle& p, int i_nuclide) p.wgt() -= wgt_absorb; // Score implicit absorption estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { p.keff_tally_absorption() += wgt_absorb * p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).absorption; @@ -656,7 +658,9 @@ void absorption(Particle& p, int i_nuclide) if (p.neutron_xs(i_nuclide).absorption > prn(p.current_seed()) * p.neutron_xs(i_nuclide).total) { // Score absorption estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { p.keff_tally_absorption() += p.wgt() * p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).absorption; diff --git a/src/settings.cpp b/src/settings.cpp index 3336e2acb2b..aec3bf5cab9 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -679,6 +679,10 @@ void read_settings_xml(pugi::xml_node root) if (check_for_node(root, "calculate_subcritical_k")) { if (run_mode == RunMode::FIXED_SOURCE) { + if (solver_type != SolverType::MONTE_CARLO) { + fatal_error("The 'calculate_subcritical_k' setting is only valid in " + "fixed source mode with the Monte Carlo solver."); + } calculate_subcritical_k = get_node_value_bool(root, "calculate_subcritical_k"); } else { diff --git a/src/simulation.cpp b/src/simulation.cpp index b536ae5881f..4c97084644b 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -142,6 +142,9 @@ int openmc_simulation_init() if (settings::run_mode == RunMode::FIXED_SOURCE) { if (settings::solver_type == SolverType::MONTE_CARLO) { header("FIXED SOURCE TRANSPORT SIMULATION", 3); + if (settings::verbosity >= 7 && settings::calculate_subcritical_k) { + print_columns(); + } } else if (settings::solver_type == SolverType::RANDOM_RAY) { header("FIXED SOURCE TRANSPORT SIMULATION (RANDOM RAY SOLVER)", 3); } @@ -366,7 +369,9 @@ void initialize_batch() write_message( 6, "Simulating batch {:<4} (inactive)", simulation::current_batch); } else { - write_message(6, "Simulating batch {}", simulation::current_batch); + if (!settings::calculate_subcritical_k) { + write_message(6, "Simulating batch {}", simulation::current_batch); + } } } @@ -505,7 +510,9 @@ void finalize_batch() void initialize_generation() { - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { // Clear out the fission bank simulation::fission_bank.resize(0); @@ -523,9 +530,9 @@ void finalize_generation() { auto& gt = simulation::global_tallies; - // Update global tallies with the accumulation variables - if (settings::run_mode == RunMode::EIGENVALUE) { - gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision; + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) += global_tally_absorption; gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += @@ -533,8 +540,73 @@ void finalize_generation() } gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage; - // reset tallies - if (settings::run_mode == RunMode::EIGENVALUE) { + // Perform MPI Reduction + // We must reduce global_tallies and total_weight across all ranks + // so that calculate_generation_keff sees the total sum. +#ifdef OPENMC_MPI + // Reduce global tallies + xt::xtensor buffer = gt; + MPI_Reduce(buffer.data(), gt.data(), gt.size(), MPI_DOUBLE, MPI_SUM, 0, + mpi::intracomm); + + // Reduce total weight (Source S) + double total_weight_reduced = 0.0; + MPI_Reduce(&simulation::total_weight, &total_weight_reduced, 1, MPI_DOUBLE, + MPI_SUM, 0, mpi::intracomm); + + // Update local variable on master with the reduced value + if (mpi::master) + simulation::total_weight = total_weight_reduced; +#endif + + // Calculate Combined Estimator and k_sq + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + // We need the Global Sum of Production (P) and Source Weight (S). + // The variables global_tally_tracklength hold the LOCAL production for this + // generation. + double local_P = global_tally_tracklength; + double local_S = simulation::total_weight; + + double global_P = 0.0; + double global_S = 0.0; + +#ifdef OPENMC_MPI + // Reduce Production and Source across all ranks + MPI_Reduce(&local_P, &global_P, 1, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); + MPI_Reduce(&local_S, &global_S, 1, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); +#else + global_P = local_P; + global_S = local_S; +#endif + + if (mpi::master) { + // Calculate k_sq = P / (P + S) + // Note: We use unnormalized totals for both P and S, so the ratio is + // correct. + double k_sq = 0.0; + if (global_P + global_S > 0.0) { + k_sq = global_P / (global_P + global_S); + } + + // Update the history vector + simulation::k_generation.push_back(k_sq); + + // Calculate the running Mean/StdDev of the k_sq history + calculate_average_keff(); + + // Print the standard generation table + if (settings::verbosity >= 7) { + print_generation(); + } + } + } + + // Reset tallies for the next generation + if (settings::run_mode == RunMode::EIGENVALUE || + settings::run_mode == RunMode::SUBCRITICAL_MULTIPLICATION || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { global_tally_collision = 0.0; global_tally_absorption = 0.0; global_tally_tracklength = 0.0; diff --git a/tests/regression_tests/fixed_source_k/inputs_true.dat b/tests/regression_tests/fixed_source_k/inputs_true.dat index af662d47749..cb621dae631 100644 --- a/tests/regression_tests/fixed_source_k/inputs_true.dat +++ b/tests/regression_tests/fixed_source_k/inputs_true.dat @@ -31,5 +31,6 @@ + true diff --git a/tests/regression_tests/fixed_source_k/test.py b/tests/regression_tests/fixed_source_k/test.py index 370ca21297d..ffac7e44a3a 100644 --- a/tests/regression_tests/fixed_source_k/test.py +++ b/tests/regression_tests/fixed_source_k/test.py @@ -50,6 +50,7 @@ def model(): model.geometry = openmc.Geometry(universe) model.settings.run_mode = 'fixed source' model.settings.source = test_source() + model.settings.calculate_subcritical_k = True model.settings.batches = 10 model.settings.inactive = 5 From 004f1eeea911e4bdf6bbec42e0189d058237e15c Mon Sep 17 00:00:00 2001 From: Matthew Louis Date: Mon, 22 Dec 2025 16:12:45 -0500 Subject: [PATCH 3/8] Initial implementation of subcritical k in fixed source mode --- CMakeLists.txt | 1 + include/openmc/eigenvalue.h | 7 +- include/openmc/simulation.h | 6 ++ include/openmc/subcritical.h | 16 ++++ include/openmc/tallies/tally.h | 6 ++ src/eigenvalue.cpp | 38 ++++----- src/output.cpp | 9 +++ src/physics.cpp | 12 ++- src/simulation.cpp | 142 +++++++++++++++++---------------- src/subcritical.cpp | 22 +++++ src/tallies/tally.cpp | 8 +- 11 files changed, 175 insertions(+), 92 deletions(-) create mode 100644 include/openmc/subcritical.h create mode 100644 src/subcritical.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 87b8789d101..9c7c2935e65 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -400,6 +400,7 @@ list(APPEND libopenmc_SOURCES src/source.cpp src/state_point.cpp src/string_utils.cpp + src/subcritical.cpp src/summary.cpp src/surface.cpp src/tallies/derivative.cpp diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index b456fee21eb..2e52f561f2a 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -22,7 +22,10 @@ namespace openmc { namespace simulation { extern double keff_generation; //!< Single-generation k on each processor +extern double kq_generation_val; //!< Single-generation kq on each processor +extern double ks_generation_val; //!< Single-generation ks on each processor extern array k_sum; //!< Used to reduce sum and sum_sq +extern array kq_sum; extern vector entropy; //!< Shannon entropy at each generation extern xt::xtensor source_frac; //!< Source fraction for UFS @@ -33,14 +36,14 @@ extern xt::xtensor source_frac; //!< Source fraction for UFS //============================================================================== //! Collect/normalize the tracklength keff from each process -void calculate_generation_keff(); +void calculate_generation_keff(xt::xtensor_fixed> gt, double& keff_generation, vector& k_generation); //! Calculate mean/standard deviation of keff during active generations //! //! This function sets the global variables keff and keff_std which represent //! the mean and standard deviation of the mean of k-effective over active //! generations. It also broadcasts the value from the master process. -void calculate_average_keff(); +void calculate_average_keff(double& keff, double& keff_std, const vector& k_generation, std::array& k_sum); //! Calculates a minimum variance estimate of k-effective //! diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 9a6cf1b2131..5d744d647cb 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -28,6 +28,10 @@ extern "C" int current_gen; //!< current fission generation extern "C" bool initialized; //!< has simulation been initialized? extern "C" double keff; //!< average k over batches extern "C" double keff_std; //!< standard deviation of average k +extern "C" double kq; //!< average kq over batches +extern "C" double kq_std; //!< standard deviation of average kq +extern "C" double ks; //!< average ks over batches +extern "C" double ks_std; //!< standard deviation of average ks extern "C" double k_col_abs; //!< sum over batches of k_collision * k_absorption extern "C" double k_col_tra; //!< sum over batches of k_collision * k_tracklength @@ -47,6 +51,8 @@ extern const RegularMesh* entropy_mesh; extern const RegularMesh* ufs_mesh; extern vector k_generation; +extern vector kq_generation; +extern vector ks_generation; extern vector work_index; } // namespace simulation diff --git a/include/openmc/subcritical.h b/include/openmc/subcritical.h new file mode 100644 index 00000000000..8383c38fb7d --- /dev/null +++ b/include/openmc/subcritical.h @@ -0,0 +1,16 @@ +//! \file subcritical.h +//! \brief Data/functions related to subcritical (fixed source) multiplication calculations + +#ifndef OPENMC_SUBCRITICAL_H +#define OPENMC_SUBCRITICAL_H + +#include + +namespace openmc { + +void convert_to_subcritical_k(double& k, double& k_std); + +double convert_to_subcritical_k(double k); + +} // namespace openmc +#endif // OPENMC_SUBCRITICAL_H \ No newline at end of file diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 374daff92a0..2ef2302aa09 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -222,6 +222,8 @@ namespace simulation { //! Global tallies (such as k-effective estimators) extern xt::xtensor_fixed> global_tallies; +extern xt::xtensor_fixed> + global_tallies_first_gen; //! Number of realizations for global tallies extern "C" int32_t n_realizations; @@ -232,6 +234,10 @@ extern double global_tally_collision; extern double global_tally_tracklength; extern double global_tally_leakage; +extern double global_tally_absorption_first_gen; +extern double global_tally_collision_first_gen; +extern double global_tally_tracklength_first_gen; + //============================================================================== // Non-member functions //============================================================================== diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 8412cbd3b47..c136afb75f4 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -37,7 +37,10 @@ namespace openmc { namespace simulation { double keff_generation; +double kq_generation_val; +double ks_generation_val; array k_sum; +array kq_sum; vector entropy; xt::xtensor source_frac; @@ -47,20 +50,18 @@ xt::xtensor source_frac; // Non-member functions //============================================================================== -void calculate_generation_keff() +void calculate_generation_keff(xt::xtensor_fixed> gt, double& keff_generation, vector& k_generation) { - const auto& gt = simulation::global_tallies; - // Get keff for this generation by subtracting off the starting value - simulation::keff_generation = + keff_generation = gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) - - simulation::keff_generation; + keff_generation; double keff_reduced; #ifdef OPENMC_MPI if (settings::solver_type != SolverType::RANDOM_RAY) { // Combine values across all processors - MPI_Allreduce(&simulation::keff_generation, &keff_reduced, 1, MPI_DOUBLE, + MPI_Allreduce(&keff_generation, &keff_reduced, 1, MPI_DOUBLE, MPI_SUM, mpi::intracomm); } else { // If using random ray, MPI parallelism is provided by domain replication. @@ -68,10 +69,10 @@ void calculate_generation_keff() // such that all ranks have identical scalar flux vectors, and will all // independently compute the same value of k. Thus, there is no need to // perform any additional MPI reduction here. - keff_reduced = simulation::keff_generation; + keff_reduced = keff_generation; } #else - keff_reduced = simulation::keff_generation; + keff_reduced = keff_generation; #endif // Normalize single batch estimate of k @@ -80,7 +81,7 @@ void calculate_generation_keff() keff_reduced /= settings::n_particles; } - simulation::k_generation.push_back(keff_reduced); + k_generation.push_back(keff_reduced); } void synchronize_bank() @@ -364,7 +365,7 @@ void synchronize_bank() simulation::time_bank.stop(); } -void calculate_average_keff() +void calculate_average_keff(double& keff, double& keff_std, const vector& k_generation, std::array& k_sum) { // Determine overall generation and number of active generations int i = overall_generation() - 1; @@ -379,15 +380,14 @@ void calculate_average_keff() if (n <= 0) { // For inactive generations, use current generation k as estimate for next // generation - simulation::keff = simulation::k_generation[i]; + keff = k_generation[i]; } else { // Sample mean of keff - simulation::k_sum[0] += simulation::k_generation[i]; - simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2); + k_sum[0] += k_generation[i]; + k_sum[1] += std::pow(k_generation[i], 2); // Determine mean - simulation::keff = simulation::k_sum[0] / n; - + keff = k_sum[0] / n; if (n > 1) { double t_value; if (settings::confidence_intervals) { @@ -399,10 +399,10 @@ void calculate_average_keff() } // Standard deviation of the sample mean of k - simulation::keff_std = + keff_std = t_value * std::sqrt( - (simulation::k_sum[1] / n - std::pow(simulation::keff, 2)) / (n - 1)); + (k_sum[1] / n - std::pow(keff, 2)) / (n - 1)); // In some cases (such as an infinite medium problem), random ray // may estimate k exactly and in an unvarying manner between iterations. @@ -410,8 +410,8 @@ void calculate_average_keff() // power operations may cause an extremely small negative value to occur // inside the sqrt operation, leading to NaN. If this occurs, we check for // it and set the std dev to zero. - if (!std::isfinite(simulation::keff_std)) { - simulation::keff_std = 0.0; + if (!std::isfinite(keff_std)) { + keff_std = 0.0; } } } diff --git a/src/output.cpp b/src/output.cpp index 06b4e42aca0..511a9d1b61b 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -373,6 +373,9 @@ void print_columns() if (settings::entropy_on) { fmt::print(" Bat./Gen. k Entropy Average k \n" " ========= ======== ======== ====================\n"); + } else if (settings::calculate_subcritical_k) { + fmt::print(" Bat./Gen. k kq ks Average k Average kq Average ks \n" + " ========= ======== ======== ======== ==================== ==================== ====================\n"); } else { fmt::print(" Bat./Gen. k Average k\n" " ========= ======== ====================\n"); @@ -394,6 +397,9 @@ void print_generation() auto batch_and_gen = std::to_string(simulation::current_batch) + "/" + std::to_string(simulation::current_gen); fmt::print(" {:>9} {:8.5f}", batch_and_gen, simulation::k_generation[idx]); + if (settings::calculate_subcritical_k) { + fmt::print(" {:8.5f} {:8.5f}", simulation::kq_generation[idx], simulation::ks_generation[idx]); + } // write out entropy info if (settings::entropy_on && !settings::calculate_subcritical_k) { @@ -402,6 +408,9 @@ void print_generation() if (n > 1) { fmt::print(" {:8.5f} +/-{:8.5f}", simulation::keff, simulation::keff_std); + if (settings::calculate_subcritical_k) { + fmt::print(" {:8.5f} +/-{:8.5f} {:8.5f} +/-{:8.5f}", simulation::kq, simulation::kq_std, simulation::ks, simulation::ks_std); + } } fmt::print("\n"); std::fflush(stdout); diff --git a/src/physics.cpp b/src/physics.cpp index 3e3bd53cd09..d34d2354d2d 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -174,9 +174,19 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p.wgt() / simulation::keff * weight * + double nu_t; + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + // In fixed source mode, no need to scale by keff + nu_t = p.wgt() * weight * p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).total; + } else { + nu_t = p.wgt() / simulation::keff * weight * + p.neutron_xs(i_nuclide).nu_fission / + p.neutron_xs(i_nuclide).total; + } + // Sample the number of neutrons produced int nu = static_cast(nu_t); diff --git a/src/simulation.cpp b/src/simulation.cpp index 4c97084644b..b29bf6b5238 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -26,6 +26,7 @@ #include "openmc/timer.h" #include "openmc/track_output.h" #include "openmc/weight_windows.h" +#include "openmc/subcritical.h" #ifdef _OPENMP #include @@ -310,6 +311,10 @@ double keff_std; double k_col_abs {0.0}; double k_col_tra {0.0}; double k_abs_tra {0.0}; +double kq; +double kq_std; +double ks; +double ks_std; double log_spacing; int n_lost_particles {0}; bool need_depletion_rx {false}; @@ -324,6 +329,8 @@ const RegularMesh* entropy_mesh {nullptr}; const RegularMesh* ufs_mesh {nullptr}; vector k_generation; +vector kq_generation; +vector ks_generation; vector work_index; } // namespace simulation @@ -419,6 +426,7 @@ void finalize_batch() // Reset global tally results if (simulation::current_batch <= settings::n_inactive) { xt::view(simulation::global_tallies, xt::all()) = 0.0; + xt::view(simulation::global_tallies_first_gen, xt::all()) = 0.0; simulation::n_realizations = 0; } @@ -523,94 +531,61 @@ void initialize_generation() // Store current value of tracklength k simulation::keff_generation = simulation::global_tallies( GlobalTally::K_TRACKLENGTH, TallyResult::VALUE); + + if (settings::calculate_subcritical_k) { + simulation::kq_generation_val = simulation::global_tallies_first_gen( + GlobalTally::K_TRACKLENGTH, TallyResult::VALUE); + } } } void finalize_generation() { auto& gt = simulation::global_tallies; + auto& gt_first_gen = simulation::global_tallies_first_gen; - if (settings::run_mode == RunMode::EIGENVALUE || - (settings::run_mode == RunMode::FIXED_SOURCE && - settings::calculate_subcritical_k)) { + if (settings::run_mode == RunMode::EIGENVALUE) { + gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += + global_tally_collision; gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) += global_tally_absorption; gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += global_tally_tracklength; } - gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage; - - // Perform MPI Reduction - // We must reduce global_tallies and total_weight across all ranks - // so that calculate_generation_keff sees the total sum. -#ifdef OPENMC_MPI - // Reduce global tallies - xt::xtensor buffer = gt; - MPI_Reduce(buffer.data(), gt.data(), gt.size(), MPI_DOUBLE, MPI_SUM, 0, - mpi::intracomm); - - // Reduce total weight (Source S) - double total_weight_reduced = 0.0; - MPI_Reduce(&simulation::total_weight, &total_weight_reduced, 1, MPI_DOUBLE, - MPI_SUM, 0, mpi::intracomm); - - // Update local variable on master with the reduced value - if (mpi::master) - simulation::total_weight = total_weight_reduced; -#endif - - // Calculate Combined Estimator and k_sq if (settings::run_mode == RunMode::FIXED_SOURCE && - settings::calculate_subcritical_k) { - // We need the Global Sum of Production (P) and Source Weight (S). - // The variables global_tally_tracklength hold the LOCAL production for this - // generation. - double local_P = global_tally_tracklength; - double local_S = simulation::total_weight; - - double global_P = 0.0; - double global_S = 0.0; - -#ifdef OPENMC_MPI - // Reduce Production and Source across all ranks - MPI_Reduce(&local_P, &global_P, 1, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); - MPI_Reduce(&local_S, &global_S, 1, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); -#else - global_P = local_P; - global_S = local_S; -#endif - - if (mpi::master) { - // Calculate k_sq = P / (P + S) - // Note: We use unnormalized totals for both P and S, so the ratio is - // correct. - double k_sq = 0.0; - if (global_P + global_S > 0.0) { - k_sq = global_P / (global_P + global_S); - } - - // Update the history vector - simulation::k_generation.push_back(k_sq); - - // Calculate the running Mean/StdDev of the k_sq history - calculate_average_keff(); + settings::calculate_subcritical_k) { + // Update global tallies + gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += + convert_to_subcritical_k(global_tally_collision); + gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) += + convert_to_subcritical_k(global_tally_absorption); + gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += + convert_to_subcritical_k(global_tally_tracklength); - // Print the standard generation table - if (settings::verbosity >= 7) { - print_generation(); - } - } + // Update first generation tallies + gt_first_gen(GlobalTally::K_ABSORPTION, TallyResult::VALUE) += + global_tally_absorption_first_gen; + gt_first_gen(GlobalTally::K_COLLISION, TallyResult::VALUE) += + global_tally_collision_first_gen; + gt_first_gen(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += + global_tally_tracklength_first_gen; } + gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage; - // Reset tallies for the next generation + // Reset tallies if (settings::run_mode == RunMode::EIGENVALUE || - settings::run_mode == RunMode::SUBCRITICAL_MULTIPLICATION || (settings::run_mode == RunMode::FIXED_SOURCE && settings::calculate_subcritical_k)) { global_tally_collision = 0.0; global_tally_absorption = 0.0; global_tally_tracklength = 0.0; } + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + global_tally_absorption_first_gen = 0.0; + global_tally_collision_first_gen = 0.0; + global_tally_tracklength_first_gen = 0.0; + } global_tally_leakage = 0.0; if (settings::run_mode == RunMode::EIGENVALUE && @@ -624,16 +599,32 @@ void finalize_generation() synchronize_bank(); } - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { // Calculate shannon entropy if (settings::entropy_on && - settings::solver_type == SolverType::MONTE_CARLO) + settings::solver_type == SolverType::MONTE_CARLO && + !settings::calculate_subcritical_k) { shannon_entropy(); + } // Collect results and statistics - calculate_generation_keff(); - calculate_average_keff(); + calculate_generation_keff(simulation::global_tallies, simulation::keff_generation, simulation::k_generation); + calculate_average_keff(simulation::keff, simulation::keff_std, simulation::k_generation, simulation::k_sum); + + if (settings::calculate_subcritical_k) { + // Compute kq and average + calculate_generation_keff(simulation::global_tallies_first_gen, simulation::kq_generation_val, simulation::kq_generation); + calculate_average_keff(simulation::kq, simulation::kq_std, simulation::kq_generation, simulation::kq_sum); + + // Calculate ks from k = kq/(1 - ks + kq) -> ks = 1 + kq*(k - 1)/k + simulation::ks_generation_val = 1 + simulation::kq_generation.back() * (simulation::k_generation.back() - 1) / simulation::k_generation.back(); + simulation::ks_generation.push_back(simulation::ks_generation_val); + simulation::ks = 1 + simulation::kq * (simulation::keff - 1) / simulation::keff; + simulation::ks_std = sqrt( pow(simulation::kq_std, 2) + pow(simulation::kq / simulation::keff, 2) * ( pow(simulation::kq_std / simulation::kq,2) + pow(simulation::keff_std / simulation::keff,2) ) ); + } // Write generation output if (mpi::master && settings::verbosity >= 7) { @@ -854,6 +845,8 @@ void broadcast_results() // Also broadcast global tally results auto& gt = simulation::global_tallies; MPI_Bcast(gt.data(), gt.size(), MPI_DOUBLE, 0, mpi::intracomm); + auto& gt_first_gen = simulation::global_tallies_first_gen; + MPI_Bcast(gt_first_gen.data(), gt_first_gen.size(), MPI_DOUBLE, 0, mpi::intracomm); // These guys are needed so that non-master processes can calculate the // combined estimate of k-effective @@ -875,6 +868,7 @@ void free_memory_simulation() void transport_history_based_single_particle(Particle& p) { + bool tally_first_generation = true; while (p.alive()) { p.event_calculate_xs(); if (p.alive()) { @@ -887,6 +881,16 @@ void transport_history_based_single_particle(Particle& p) p.event_collide(); } } + // Check for first generation completion + if ((!p.alive() || p.n_progeny() > 0) && tally_first_generation) { + if (settings::calculate_subcritical_k) { + global_tally_absorption_first_gen += p.keff_tally_absorption(); + global_tally_collision_first_gen += p.keff_tally_collision(); + global_tally_tracklength_first_gen += p.keff_tally_tracklength(); + } + tally_first_generation = false; + } + p.event_revive_from_secondary(); } p.event_death(); diff --git a/src/subcritical.cpp b/src/subcritical.cpp new file mode 100644 index 00000000000..766bcfb9152 --- /dev/null +++ b/src/subcritical.cpp @@ -0,0 +1,22 @@ +#include "openmc/subcritical.h" +#include "openmc/simulation.h" +#include +#include + +namespace openmc{ + +void convert_to_subcritical_k(double& k, double& k_std) { + double k_sub = k / (k + 1); + double k_sub_std = k_sub * sqrt( pow(k_std / k, 2) + pow(k_std / (k + 1), 2) ); + k = k_sub; + k_std = k_sub_std; +} + +double convert_to_subcritical_k(double k) { + // Used to convert keff estimators in fixed source mode, which are multiplicities, specifically = M-1 + // into the corresponding estimators for subcritical k + double k_sub = (k/simulation::total_weight) / (k/simulation::total_weight + 1) * simulation::total_weight; + return k_sub; +} + +} // namespace openmc \ No newline at end of file diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 6eef1da9cfd..c757815aa01 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -70,6 +70,7 @@ vector time_grid; namespace simulation { xt::xtensor_fixed> global_tallies; +xt::xtensor_fixed> global_tallies_first_gen; int32_t n_realizations {0}; } // namespace simulation @@ -78,6 +79,10 @@ double global_tally_collision; double global_tally_tracklength; double global_tally_leakage; +double global_tally_absorption_first_gen; +double global_tally_collision_first_gen; +double global_tally_tracklength_first_gen; + //============================================================================== // Tally object implementation //============================================================================== @@ -1072,7 +1077,8 @@ void accumulate_tallies() if (mpi::master || !settings::reduce_tallies) { auto& gt = simulation::global_tallies; - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && settings::calculate_subcritical_k)) { if (simulation::current_batch > settings::n_inactive) { // Accumulate products of different estimators of k double k_col = gt(GlobalTally::K_COLLISION, TallyResult::VALUE) / From ebf449e2dabbeeabcd99d0cb867dc30fe621eedf Mon Sep 17 00:00:00 2001 From: Matthew Louis Date: Tue, 23 Dec 2025 07:55:12 -0500 Subject: [PATCH 4/8] Implemented combined, generation, and average calculations of subcritical values, and statepoint writing --- include/openmc/simulation.h | 6 + include/openmc/subcritical.h | 9 ++ openmc/statepoint.py | 95 +++++++++++++--- src/finalize.cpp | 3 + src/simulation.cpp | 22 +++- src/state_point.cpp | 4 + src/subcritical.cpp | 206 +++++++++++++++++++++++++++++++++++ src/tallies/tally.cpp | 22 ++++ 8 files changed, 346 insertions(+), 21 deletions(-) diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 5d744d647cb..60cfe4d746c 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -37,6 +37,12 @@ extern "C" double k_col_tra; //!< sum over batches of k_collision * k_tracklength extern "C" double k_abs_tra; //!< sum over batches of k_absorption * k_tracklength +extern "C" double + kq_col_abs; //!< sum over batches of kq_collision * kq_absorption +extern "C" double + kq_col_tra; //!< sum over batches of kq_collision * kq_tracklength +extern "C" double + kq_abs_tra; //!< sum over batches of kq_absorption * kq_tracklength extern double log_spacing; //!< lethargy spacing for energy grid searches extern "C" int n_lost_particles; //!< cumulative number of lost particles extern "C" bool need_depletion_rx; //!< need to calculate depletion rx? diff --git a/include/openmc/subcritical.h b/include/openmc/subcritical.h index 8383c38fb7d..a6df6715d97 100644 --- a/include/openmc/subcritical.h +++ b/include/openmc/subcritical.h @@ -4,6 +4,7 @@ #ifndef OPENMC_SUBCRITICAL_H #define OPENMC_SUBCRITICAL_H +#include "hdf5.h" #include namespace openmc { @@ -12,5 +13,13 @@ void convert_to_subcritical_k(double& k, double& k_std); double convert_to_subcritical_k(double k); +double calculate_ks(double k, double kq); + +double calculate_sigma_ks(double k, double k_std, double kq, double kq_std); + +extern "C" int openmc_get_subcritical_kq(double* k_combined); + +void write_subcritical_hdf5(hid_t group); + } // namespace openmc #endif // OPENMC_SUBCRITICAL_H \ No newline at end of file diff --git a/openmc/statepoint.py b/openmc/statepoint.py index a10ec3a839d..5bedcc61dd1 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -77,10 +77,24 @@ class StatePoint: Cross-product of collision and tracklength estimates of k-effective k_abs_tra : float Cross-product of absorption and tracklength estimates of k-effective + kq_col_abs : float + Cross-product of collision and absorption estimates of kq + kq_col_tra : float + Cross-product of collision and tracklength estimates of kq + kq_abs_tra : float + Cross-product of absorption and tracklength estimates of kq k_generation : numpy.ndarray Estimate of k-effective for each batch/generation + kq_generation : numpy.ndarray + Estimate of kq for each batch/generation + ks_generation : numpy.ndarray + Estimate of ks for each batch/generation keff : uncertainties.UFloat Combined estimator for k-effective + kq : uncertainties.UFloat + Combined estimator for kq + ks : uncertainties.UFloat + Combined estimator for ks .. versionadded:: 0.13.1 meshes : dict @@ -212,9 +226,9 @@ def date_and_time(self): @property def entropy(self): - if self.run_mode == 'eigenvalue': + try: return self._f['entropy'][()] - else: + except KeyError: return None @property @@ -233,9 +247,9 @@ def filters(self): @property def generations_per_batch(self): - if self.run_mode == 'eigenvalue': + try: return self._f['generations_per_batch'][()] - else: + except KeyError: return None @property @@ -268,16 +282,30 @@ def k_cmfd(self): @property def k_generation(self): - if self.run_mode == 'eigenvalue': + try: return self._f['k_generation'][()] - else: + except KeyError: + return None + + @property + def kq_generation(self): + try: + return self._f['kq_generation'][()] + except KeyError: + return None + + @property + def ks_generation(self): + try: + return self._f['ks_generation'][()] + except KeyError: return None @property def keff(self): - if self.run_mode == 'eigenvalue': + try: return ufloat(*self._f['k_combined'][()]) - else: + except KeyError: return None @property @@ -287,26 +315,61 @@ def k_combined(self): "removed in a future version of OpenMC.", FutureWarning ) return self.keff + + @property + def kq(self): + try: + return ufloat(*self._f['kq_combined'][()]) + except KeyError: + return None + + @property + def ks(self): + try: + return ufloat(*self._f['ks_combined'][()]) + except KeyError: + return None @property def k_col_abs(self): - if self.run_mode == 'eigenvalue': + try: return self._f['k_col_abs'][()] - else: + except KeyError: return None @property def k_col_tra(self): - if self.run_mode == 'eigenvalue': + try: return self._f['k_col_tra'][()] - else: + except KeyError: return None @property def k_abs_tra(self): - if self.run_mode == 'eigenvalue': + try: return self._f['k_abs_tra'][()] - else: + except KeyError: + return None + + @property + def kq_col_abs(self): + try: + return self._f['kq_col_abs'][()] + except KeyError: + return None + + @property + def kq_col_tra(self): + try: + return self._f['kq_col_tra'][()] + except KeyError: + return None + + @property + def kq_abs_tra(self): + try: + return self._f['kq_abs_tra'][()] + except KeyError: return None @property @@ -329,9 +392,9 @@ def n_batches(self): @property def n_inactive(self): - if self.run_mode == 'eigenvalue': + try: return self._f['n_inactive'][()] - else: + except KeyError: return None @property diff --git a/src/finalize.cpp b/src/finalize.cpp index 02ca5a8d5c9..6073dd65f16 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -206,6 +206,9 @@ int openmc_reset() simulation::k_col_abs = 0.0; simulation::k_col_tra = 0.0; simulation::k_abs_tra = 0.0; + simulation::kq_col_abs = 0.0; + simulation::kq_col_tra = 0.0; + simulation::kq_abs_tra = 0.0; simulation::k_sum = {0.0, 0.0}; simulation::satisfy_triggers = false; diff --git a/src/simulation.cpp b/src/simulation.cpp index b29bf6b5238..a4512334b88 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -311,6 +311,9 @@ double keff_std; double k_col_abs {0.0}; double k_col_tra {0.0}; double k_abs_tra {0.0}; +double kq_col_abs {0.0}; +double kq_col_tra {0.0}; +double kq_abs_tra {0.0}; double kq; double kq_std; double ks; @@ -618,12 +621,14 @@ void finalize_generation() // Compute kq and average calculate_generation_keff(simulation::global_tallies_first_gen, simulation::kq_generation_val, simulation::kq_generation); calculate_average_keff(simulation::kq, simulation::kq_std, simulation::kq_generation, simulation::kq_sum); - - // Calculate ks from k = kq/(1 - ks + kq) -> ks = 1 + kq*(k - 1)/k - simulation::ks_generation_val = 1 + simulation::kq_generation.back() * (simulation::k_generation.back() - 1) / simulation::k_generation.back(); + + // Calculate ks from k, kq + simulation::ks_generation_val = calculate_ks( + simulation::k_generation.back(), simulation::kq_generation.back()); simulation::ks_generation.push_back(simulation::ks_generation_val); - simulation::ks = 1 + simulation::kq * (simulation::keff - 1) / simulation::keff; - simulation::ks_std = sqrt( pow(simulation::kq_std, 2) + pow(simulation::kq / simulation::keff, 2) * ( pow(simulation::kq_std / simulation::kq,2) + pow(simulation::keff_std / simulation::keff,2) ) ); + simulation::ks = calculate_ks(simulation::keff, simulation::kq); + simulation::ks_std = calculate_sigma_ks(simulation::keff, + simulation::keff_std, simulation::kq, simulation::kq_std); } // Write generation output @@ -856,6 +861,13 @@ void broadcast_results() simulation::k_col_abs = temp[0]; simulation::k_col_tra = temp[1]; simulation::k_abs_tra = temp[2]; + + double temp_kq[] { + simulation::kq_col_abs, simulation::kq_col_tra, simulation::kq_abs_tra}; + MPI_Bcast(temp_kq, 3, MPI_DOUBLE, 0, mpi::intracomm); + simulation::kq_col_abs = temp_kq[0]; + simulation::kq_col_tra = temp_kq[1]; + simulation::kq_abs_tra = temp_kq[2]; } #endif diff --git a/src/state_point.cpp b/src/state_point.cpp index 8ccebeb05a7..1cc1476a0ad 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -24,6 +24,7 @@ #include "openmc/output.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/subcritical.h" #include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" #include "openmc/tallies/filter_mesh.h" @@ -120,6 +121,9 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) if (settings::run_mode == RunMode::EIGENVALUE) write_eigenvalue_hdf5(file_id); + if (settings::calculate_subcritical_k) + write_subcritical_hdf5(file_id); + hid_t tallies_group = create_group(file_id, "tallies"); // Write meshes diff --git a/src/subcritical.cpp b/src/subcritical.cpp index 766bcfb9152..c7627a7bc59 100644 --- a/src/subcritical.cpp +++ b/src/subcritical.cpp @@ -1,5 +1,8 @@ #include "openmc/subcritical.h" +#include "hdf5.h" +#include "openmc/eigenvalue.h" #include "openmc/simulation.h" +#include "openmc/tallies/tally.h" #include #include @@ -19,4 +22,207 @@ double convert_to_subcritical_k(double k) { return k_sub; } +int openmc_get_subcritical_kq(double* k_combined) +{ + k_combined[0] = 0.0; + k_combined[1] = 0.0; + + // Special case for n <=3. Notice that at the end, + // there is a N-3 term in a denominator. + if (simulation::n_realizations <= 3 || + settings::solver_type == SolverType::RANDOM_RAY) { + k_combined[0] = simulation::kq; + k_combined[1] = simulation::kq_std; + if (simulation::n_realizations <= 1) { + k_combined[1] = std::numeric_limits::infinity(); + } + return 0; + } + + // Initialize variables + int64_t n = simulation::n_realizations; + + // Copy estimates of k-effective and its variance (not variance of the mean) + const auto& gt = simulation::global_tallies_first_gen; + + array kv {}; + xt::xtensor cov = xt::zeros({3, 3}); + kv[0] = gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n; + kv[1] = gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n; + kv[2] = gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n; + cov(0, 0) = + (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n * kv[0] * kv[0]) / + (n - 1); + cov(1, 1) = + (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - n * kv[1] * kv[1]) / + (n - 1); + cov(2, 2) = + (gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM_SQ) - n * kv[2] * kv[2]) / + (n - 1); + + // Calculate covariances based on sums with Bessel's correction + cov(0, 1) = (simulation::kq_col_abs - n * kv[0] * kv[1]) / (n - 1); + cov(0, 2) = (simulation::kq_col_tra - n * kv[0] * kv[2]) / (n - 1); + cov(1, 2) = (simulation::kq_abs_tra - n * kv[1] * kv[2]) / (n - 1); + cov(1, 0) = cov(0, 1); + cov(2, 0) = cov(0, 2); + cov(2, 1) = cov(1, 2); + + // Check to see if two estimators are the same; this is guaranteed to happen + // in MG-mode with survival biasing when the collision and absorption + // estimators are the same, but can theoretically happen at anytime. + // If it does, the standard estimators will produce floating-point + // exceptions and an expression specifically derived for the combination of + // two estimators (vice three) should be used instead. + + // First we will identify if there are any matching estimators + int i, j; + bool use_three = false; + if ((std::abs(kv[0] - kv[1]) / kv[0] < FP_REL_PRECISION) && + (std::abs(cov(0, 0) - cov(1, 1)) / cov(0, 0) < FP_REL_PRECISION)) { + // 0 and 1 match, so only use 0 and 2 in our comparisons + i = 0; + j = 2; + + } else if ((std::abs(kv[0] - kv[2]) / kv[0] < FP_REL_PRECISION) && + (std::abs(cov(0, 0) - cov(2, 2)) / cov(0, 0) < FP_REL_PRECISION)) { + // 0 and 2 match, so only use 0 and 1 in our comparisons + i = 0; + j = 1; + + } else if ((std::abs(kv[1] - kv[2]) / kv[1] < FP_REL_PRECISION) && + (std::abs(cov(1, 1) - cov(2, 2)) / cov(1, 1) < FP_REL_PRECISION)) { + // 1 and 2 match, so only use 0 and 1 in our comparisons + i = 0; + j = 1; + + } else { + // No two estimators match, so set boolean to use all three estimators. + use_three = true; + } + + if (use_three) { + // Use three estimators as derived in the paper by Urbatsch + + // Initialize variables + double g = 0.0; + array S {}; + + for (int l = 0; l < 3; ++l) { + // Permutations of estimates + int k; + switch (l) { + case 0: + // i = collision, j = absorption, k = tracklength + i = 0; + j = 1; + k = 2; + break; + case 1: + // i = absortion, j = tracklength, k = collision + i = 1; + j = 2; + k = 0; + break; + case 2: + // i = tracklength, j = collision, k = absorption + i = 2; + j = 0; + k = 1; + break; + } + + // Calculate weighting + double f = cov(j, j) * (cov(k, k) - cov(i, k)) - cov(k, k) * cov(i, j) + + cov(j, k) * (cov(i, j) + cov(i, k) - cov(j, k)); + + // Add to S sums for variance of combined estimate + S[0] += f * cov(0, l); + S[1] += (cov(j, j) + cov(k, k) - 2.0 * cov(j, k)) * kv[l] * kv[l]; + S[2] += (cov(k, k) + cov(i, j) - cov(j, k) - cov(i, k)) * kv[l] * kv[j]; + + // Add to sum for combined k-effective + k_combined[0] += f * kv[l]; + g += f; + } + + // Complete calculations of S sums + for (auto& S_i : S) { + S_i *= (n - 1); + } + S[0] *= (n - 1) * (n - 1); + + // Calculate combined estimate of k-effective + k_combined[0] /= g; + + // Calculate standard deviation of combined estimate + g *= (n - 1) * (n - 1); + k_combined[1] = + std::sqrt(S[0] / (g * n * (n - 3)) * (1 + n * ((S[1] - 2 * S[2]) / g))); + + } else { + // Use only two estimators + // These equations are derived analogously to that done in the paper by + // Urbatsch, but are simpler than for the three estimators case since the + // block matrices of the three estimator equations reduces to scalars here + + // Store the commonly used term + double f = kv[i] - kv[j]; + double g = cov(i, i) + cov(j, j) - 2.0 * cov(i, j); + + // Calculate combined estimate of k-effective + k_combined[0] = kv[i] - (cov(i, i) - cov(i, j)) / g * f; + + // Calculate standard deviation of combined estimate + k_combined[1] = (cov(i, i) * cov(j, j) - cov(i, j) * cov(i, j)) * + (g + n * f * f) / (n * (n - 2) * g * g); + k_combined[1] = std::sqrt(k_combined[1]); + } + return 0; +} + +double calculate_ks(double k, double kq) +{ + // Calculate ks from k = kq/(1 - ks + kq) -> ks = 1 + kq*(k - 1)/k + double ks = 1 + kq * (k - 1) / k; + return ks; +} + +double calculate_sigma_ks(double k, double k_std, double kq, double kq_std) +{ + double sigma_ks = + sqrt(pow(kq_std, 2) + + pow(kq / k, 2) * (pow(kq_std / kq, 2) + pow(k_std / k, 2))); + return sigma_ks; +} + +void write_subcritical_hdf5(hid_t group) +{ + write_dataset(group, "n_inactive", settings::n_inactive); + write_dataset(group, "generations_per_batch", settings::gen_per_batch); + write_dataset(group, "k_generation", simulation::k_generation); + if (settings::entropy_on) { + write_dataset(group, "entropy", simulation::entropy); + } + write_dataset(group, "k_col_abs", simulation::k_col_abs); + write_dataset(group, "k_col_tra", simulation::k_col_tra); + write_dataset(group, "k_abs_tra", simulation::k_abs_tra); + write_dataset(group, "kq_col_abs", simulation::kq_col_abs); + write_dataset(group, "kq_col_tra", simulation::kq_col_tra); + write_dataset(group, "kq_abs_tra", simulation::kq_abs_tra); + array k_combined; + openmc_get_keff(k_combined.data()); + array kq_combined; + openmc_get_subcritical_kq(kq_combined.data()); + write_dataset(group, "k_combined", k_combined); + write_dataset(group, "kq_combined", kq_combined); + array ks_combined; + ks_combined[0] = calculate_ks(k_combined[0], kq_combined[0]); + ks_combined[1] = calculate_sigma_ks( + k_combined[0], k_combined[1], kq_combined[0], kq_combined[1]); + write_dataset(group, "ks_combined", ks_combined); + write_dataset(group, "kq_generation", simulation::kq_generation); + write_dataset(group, "ks_generation", simulation::ks_generation); +} + } // namespace openmc \ No newline at end of file diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index c757815aa01..888e9bf72d5 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -1076,6 +1076,7 @@ void accumulate_tallies() // Accumulate on master only unless run is not reduced then do it on all if (mpi::master || !settings::reduce_tallies) { auto& gt = simulation::global_tallies; + auto& gt_first_gen = simulation::global_tallies_first_gen; if (settings::run_mode == RunMode::EIGENVALUE || (settings::run_mode == RunMode::FIXED_SOURCE && settings::calculate_subcritical_k)) { @@ -1091,6 +1092,21 @@ void accumulate_tallies() simulation::k_col_tra += k_col * k_tra; simulation::k_abs_tra += k_abs * k_tra; } + if (settings::calculate_subcritical_k) { + // Accumulate products of different estimators of k + double k_col = + gt_first_gen(GlobalTally::K_COLLISION, TallyResult::VALUE) / + simulation::total_weight; + double k_abs = + gt_first_gen(GlobalTally::K_ABSORPTION, TallyResult::VALUE) / + simulation::total_weight; + double k_tra = + gt_first_gen(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) / + simulation::total_weight; + simulation::kq_col_abs += k_col * k_abs; + simulation::kq_col_tra += k_col * k_tra; + simulation::kq_abs_tra += k_abs * k_tra; + } } // Accumulate results for global tallies @@ -1099,6 +1115,12 @@ void accumulate_tallies() gt(i, TallyResult::VALUE) = 0.0; gt(i, TallyResult::SUM) += val; gt(i, TallyResult::SUM_SQ) += val * val; + + double val_first_gen = + gt_first_gen(i, TallyResult::VALUE) / simulation::total_weight; + gt_first_gen(i, TallyResult::VALUE) = 0.0; + gt_first_gen(i, TallyResult::SUM) += val_first_gen; + gt_first_gen(i, TallyResult::SUM_SQ) += val_first_gen * val_first_gen; } } From ad0ec0eb7a75cb7d97d09ccc7be7c3c9aa2b8cc0 Mon Sep 17 00:00:00 2001 From: Matthew Louis Date: Tue, 23 Dec 2025 08:12:01 -0500 Subject: [PATCH 5/8] Fixed tally update causing race condition and non determinism in simulation.cpp, updated testing, and fixed redundancies from previous commit --- src/simulation.cpp | 14 ++++--- .../{results_true.dat => __init__.py} | 0 .../{inputs_true.dat => inputs_true_1.dat} | 11 +++-- .../fixed_source_k/inputs_true_2.dat | 40 ++++++++++++++++++ .../fixed_source_k/inputs_true_3.dat | 40 ++++++++++++++++++ .../fixed_source_k/results_true_1.dat | 3 ++ .../fixed_source_k/results_true_2.dat | 9 ++++ tests/regression_tests/fixed_source_k/test.py | 42 ++++++++++++++++--- tests/testing_harness.py | 35 +++++++++++----- 9 files changed, 170 insertions(+), 24 deletions(-) rename tests/regression_tests/fixed_source_k/{results_true.dat => __init__.py} (100%) rename tests/regression_tests/fixed_source_k/{inputs_true.dat => inputs_true_1.dat} (89%) create mode 100644 tests/regression_tests/fixed_source_k/inputs_true_2.dat create mode 100644 tests/regression_tests/fixed_source_k/inputs_true_3.dat create mode 100644 tests/regression_tests/fixed_source_k/results_true_1.dat create mode 100644 tests/regression_tests/fixed_source_k/results_true_2.dat diff --git a/src/simulation.cpp b/src/simulation.cpp index a4512334b88..2027b60650d 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -896,11 +896,15 @@ void transport_history_based_single_particle(Particle& p) // Check for first generation completion if ((!p.alive() || p.n_progeny() > 0) && tally_first_generation) { if (settings::calculate_subcritical_k) { - global_tally_absorption_first_gen += p.keff_tally_absorption(); - global_tally_collision_first_gen += p.keff_tally_collision(); - global_tally_tracklength_first_gen += p.keff_tally_tracklength(); - } - tally_first_generation = false; +// Protect global updates with atomic to prevent data races +#pragma omp atomic + global_tally_absorption_first_gen += p.keff_tally_absorption(); +#pragma omp atomic + global_tally_collision_first_gen += p.keff_tally_collision(); +#pragma omp atomic + global_tally_tracklength_first_gen += p.keff_tally_tracklength(); + } + tally_first_generation = false; } p.event_revive_from_secondary(); diff --git a/tests/regression_tests/fixed_source_k/results_true.dat b/tests/regression_tests/fixed_source_k/__init__.py similarity index 100% rename from tests/regression_tests/fixed_source_k/results_true.dat rename to tests/regression_tests/fixed_source_k/__init__.py diff --git a/tests/regression_tests/fixed_source_k/inputs_true.dat b/tests/regression_tests/fixed_source_k/inputs_true_1.dat similarity index 89% rename from tests/regression_tests/fixed_source_k/inputs_true.dat rename to tests/regression_tests/fixed_source_k/inputs_true_1.dat index cb621dae631..b2d62492020 100644 --- a/tests/regression_tests/fixed_source_k/inputs_true.dat +++ b/tests/regression_tests/fixed_source_k/inputs_true_1.dat @@ -21,9 +21,8 @@ fixed source - 1000 + 100 10 - 5 0 0 0 @@ -31,6 +30,10 @@ - true - + + + nu-fission + + + diff --git a/tests/regression_tests/fixed_source_k/inputs_true_2.dat b/tests/regression_tests/fixed_source_k/inputs_true_2.dat new file mode 100644 index 00000000000..f9fa805299d --- /dev/null +++ b/tests/regression_tests/fixed_source_k/inputs_true_2.dat @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + fixed source + 100 + 10 + + + 0 0 0 + + + + + true + + + + nu-fission + + + diff --git a/tests/regression_tests/fixed_source_k/inputs_true_3.dat b/tests/regression_tests/fixed_source_k/inputs_true_3.dat new file mode 100644 index 00000000000..4c512f95caa --- /dev/null +++ b/tests/regression_tests/fixed_source_k/inputs_true_3.dat @@ -0,0 +1,40 @@ + + + + + + + + + + + + + + + + + + + + + + + fixed source + 100 + 10 + + + 0 0 0 + + + + + true + + + + nu-fission + + + diff --git a/tests/regression_tests/fixed_source_k/results_true_1.dat b/tests/regression_tests/fixed_source_k/results_true_1.dat new file mode 100644 index 00000000000..2fd4547f120 --- /dev/null +++ b/tests/regression_tests/fixed_source_k/results_true_1.dat @@ -0,0 +1,3 @@ +tally 1: +1.069325E+02 +1.167481E+03 diff --git a/tests/regression_tests/fixed_source_k/results_true_2.dat b/tests/regression_tests/fixed_source_k/results_true_2.dat new file mode 100644 index 00000000000..3e1dc66dcad --- /dev/null +++ b/tests/regression_tests/fixed_source_k/results_true_2.dat @@ -0,0 +1,9 @@ +k-combined: +9.119416E-01 4.753186E-03 +kq: +1.060436E+00 1.137258E-02 +ks: +8.976028E-01 1.793289E-02 +tally 1: +1.069325E+02 +1.167481E+03 diff --git a/tests/regression_tests/fixed_source_k/test.py b/tests/regression_tests/fixed_source_k/test.py index ffac7e44a3a..81f7d073c0f 100644 --- a/tests/regression_tests/fixed_source_k/test.py +++ b/tests/regression_tests/fixed_source_k/test.py @@ -46,6 +46,25 @@ def test_source(): def model(): model = openmc.Model() + universe = create_universe() + model.geometry = openmc.Geometry(universe) + model.settings.run_mode = 'fixed source' + model.settings.source = test_source() + + model.settings.batches = 10 + model.settings.particles = 100 + + model.tallies = openmc.Tallies() + tally = openmc.Tally(name='νΣf tally') + tally.scores = ['nu-fission'] + model.tallies.append(tally) + + return model + +@pytest.fixture +def model_k(): + model = openmc.Model() + universe = create_universe() model.geometry = openmc.Geometry(universe) model.settings.run_mode = 'fixed source' @@ -53,12 +72,25 @@ def model(): model.settings.calculate_subcritical_k = True model.settings.batches = 10 - model.settings.inactive = 5 - model.settings.particles = 1000 + model.settings.particles = 100 + + model.tallies = openmc.Tallies() + tally = openmc.Tally(name='νΣf tally') + tally.scores = ['nu-fission'] + model.tallies.append(tally) return model -def test_fixed_source_run(model): - """Test that fixed source run with subcritical k calculation works.""" - harness = PyAPITestHarness("statepoint.10.h5", model, inputs_true='inputs_true.dat') +def test_tally_results(model, model_k): + """Test that fixed source run with subcritical k calculation gives the same tally results.""" + harness = PyAPITestHarness("statepoint.10.h5", model, inputs_true='inputs_true_1.dat', results_true='results_true_1.dat') + harness.main() + + harness_k = PyAPITestHarness("statepoint.10.h5", model_k, inputs_true='inputs_true_2.dat', results_true='results_true_1.dat') + harness_k.main() + +def test_subcritical_k(model_k): + """Test that subcritical k calculation gives a known result.""" + harness = PyAPITestHarness("statepoint.10.h5", model_k, inputs_true='inputs_true_3.dat', results_true='results_true_2.dat', + subcritical_k_results=True) harness.main() \ No newline at end of file diff --git a/tests/testing_harness.py b/tests/testing_harness.py index 1ad91b7a89f..ac04a997a64 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -32,8 +32,10 @@ def colorize(diff): class TestHarness: """General class for running OpenMC regression tests.""" - def __init__(self, statepoint_name): + def __init__(self, statepoint_name, subcritical_k_results=False, results_true='results_true.dat'): self._sp_name = statepoint_name + self.subcritical_k_results = subcritical_k_results + self.results_true = results_true def main(self): """Accept commandline arguments and either run or update tests.""" @@ -94,6 +96,19 @@ def _get_results(self, hash_output=False): outstr += 'k-combined:\n' form = '{0:12.6E} {1:12.6E}\n' outstr += form.format(sp.keff.n, sp.keff.s) + if self.subcritical_k_results: + # Write out subcritical k-combined. + outstr += 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.keff.n, sp.keff.s) + outstr += 'kq:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.kq.n, + sp.kq.s) + outstr += 'ks:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.ks.n, + sp.ks.s) # Write out tally data. for i, tally_ind in enumerate(sp.tallies): @@ -125,15 +140,15 @@ def _write_results(self, results_string): def _overwrite_results(self): """Overwrite the results_true with the results_test.""" - shutil.copyfile('results_test.dat', 'results_true.dat') + shutil.copyfile('results_test.dat', self.results_true) def _compare_results(self): """Make sure the current results agree with the reference.""" - compare = filecmp.cmp('results_test.dat', 'results_true.dat') + compare = filecmp.cmp('results_test.dat', self.results_true) if not compare: - expected = open('results_true.dat').readlines() + expected = open(self.results_true).readlines() actual = open('results_test.dat').readlines() - diff = unified_diff(expected, actual, 'results_true.dat', + diff = unified_diff(expected, actual, self.results_true, 'results_test.dat') print('Result differences:') print(''.join(colorize(diff))) @@ -286,8 +301,8 @@ def _cleanup(self): class PyAPITestHarness(TestHarness): - def __init__(self, statepoint_name, model=None, inputs_true=None): - super().__init__(statepoint_name) + def __init__(self, statepoint_name, model=None, inputs_true=None, results_true='results_true.dat', subcritical_k_results=False): + super().__init__(statepoint_name, subcritical_k_results=subcritical_k_results, results_true=results_true) if model is None: self._model = pwr_core() else: @@ -441,11 +456,11 @@ def compare_tokens(token1, token2): def _compare_results(self): """Make sure the current results agree with the reference.""" compare = self._are_files_equal( - 'results_test.dat', 'results_true.dat', 1e-6) + 'results_test.dat', self.results_true, 1e-6) if not compare: - expected = open('results_true.dat').readlines() + expected = open(self.results_true).readlines() actual = open('results_test.dat').readlines() - diff = unified_diff(expected, actual, 'results_true.dat', + diff = unified_diff(expected, actual, self.results_true, 'results_test.dat') print('Result differences:') print(''.join(colorize(diff))) From f72582cc275f1367f23b0f257df2937733adf9f3 Mon Sep 17 00:00:00 2001 From: Matthew Louis Date: Tue, 23 Dec 2025 09:09:56 -0500 Subject: [PATCH 6/8] Added MPI test and updated logic --- .../fixed_source_k/inputs_true_2.dat | 14 +++++----- .../fixed_source_k/inputs_true_3.dat | 14 +++++----- tests/regression_tests/fixed_source_k/test.py | 28 ++++++++++++++++++- 3 files changed, 41 insertions(+), 15 deletions(-) diff --git a/tests/regression_tests/fixed_source_k/inputs_true_2.dat b/tests/regression_tests/fixed_source_k/inputs_true_2.dat index f9fa805299d..cfd74f2ef29 100644 --- a/tests/regression_tests/fixed_source_k/inputs_true_2.dat +++ b/tests/regression_tests/fixed_source_k/inputs_true_2.dat @@ -1,23 +1,23 @@ - + - + - - - - + + + + fixed source @@ -33,7 +33,7 @@ true - + nu-fission diff --git a/tests/regression_tests/fixed_source_k/inputs_true_3.dat b/tests/regression_tests/fixed_source_k/inputs_true_3.dat index 4c512f95caa..cfd74f2ef29 100644 --- a/tests/regression_tests/fixed_source_k/inputs_true_3.dat +++ b/tests/regression_tests/fixed_source_k/inputs_true_3.dat @@ -1,23 +1,23 @@ - + - + - - - - + + + + fixed source @@ -33,7 +33,7 @@ true - + nu-fission diff --git a/tests/regression_tests/fixed_source_k/test.py b/tests/regression_tests/fixed_source_k/test.py index 81f7d073c0f..639e6f8cd6d 100644 --- a/tests/regression_tests/fixed_source_k/test.py +++ b/tests/regression_tests/fixed_source_k/test.py @@ -1,7 +1,11 @@ +from glob import glob import openmc import pytest from tests.testing_harness import PyAPITestHarness +from tests.regression_tests import config +import os +import glob def create_universe(): # Define materials @@ -63,6 +67,7 @@ def model(): @pytest.fixture def model_k(): + openmc.reset_auto_ids() model = openmc.Model() universe = create_universe() @@ -81,6 +86,10 @@ def model_k(): return model +def remove_xml_files(): + for f in glob.glob("*.xml"): + os.remove(f) + def test_tally_results(model, model_k): """Test that fixed source run with subcritical k calculation gives the same tally results.""" harness = PyAPITestHarness("statepoint.10.h5", model, inputs_true='inputs_true_1.dat', results_true='results_true_1.dat') @@ -88,9 +97,26 @@ def test_tally_results(model, model_k): harness_k = PyAPITestHarness("statepoint.10.h5", model_k, inputs_true='inputs_true_2.dat', results_true='results_true_1.dat') harness_k.main() + remove_xml_files() def test_subcritical_k(model_k): """Test that subcritical k calculation gives a known result.""" harness = PyAPITestHarness("statepoint.10.h5", model_k, inputs_true='inputs_true_3.dat', results_true='results_true_2.dat', subcritical_k_results=True) - harness.main() \ No newline at end of file + harness.main() + remove_xml_files() + +def test_subcritical_k_mpi_specific(model_k): + if not config['mpi']: + remove_xml_files() + pytest.skip("This test is intended for MPI execution only") + + harness = PyAPITestHarness( + "statepoint.10.h5", + model_k, + inputs_true='inputs_true_3.dat', + results_true='results_true_2.dat', + subcritical_k_results=True + ) + harness.main() + remove_xml_files() \ No newline at end of file From e30c2da524e161acbc071261d4fcdf8711af9430 Mon Sep 17 00:00:00 2001 From: Matthew Louis Date: Tue, 23 Dec 2025 09:26:10 -0500 Subject: [PATCH 7/8] Refactored generation k and average calculation to minimize interface at cost of code duplication --- include/openmc/eigenvalue.h | 4 +- include/openmc/subcritical.h | 4 ++ src/eigenvalue.cpp | 35 +++++++------- src/simulation.cpp | 8 ++-- src/subcritical.cpp | 92 +++++++++++++++++++++++++++++++++++- 5 files changed, 120 insertions(+), 23 deletions(-) diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index 2e52f561f2a..1b1c2432bc0 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -36,14 +36,14 @@ extern xt::xtensor source_frac; //!< Source fraction for UFS //============================================================================== //! Collect/normalize the tracklength keff from each process -void calculate_generation_keff(xt::xtensor_fixed> gt, double& keff_generation, vector& k_generation); +void calculate_generation_keff(); //! Calculate mean/standard deviation of keff during active generations //! //! This function sets the global variables keff and keff_std which represent //! the mean and standard deviation of the mean of k-effective over active //! generations. It also broadcasts the value from the master process. -void calculate_average_keff(double& keff, double& keff_std, const vector& k_generation, std::array& k_sum); +void calculate_average_keff(); //! Calculates a minimum variance estimate of k-effective //! diff --git a/include/openmc/subcritical.h b/include/openmc/subcritical.h index a6df6715d97..1479e82ba46 100644 --- a/include/openmc/subcritical.h +++ b/include/openmc/subcritical.h @@ -13,6 +13,10 @@ void convert_to_subcritical_k(double& k, double& k_std); double convert_to_subcritical_k(double k); +void calculate_generation_kq(); + +void calculate_average_kq(); + double calculate_ks(double k, double kq); double calculate_sigma_ks(double k, double k_std, double kq, double kq_std); diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index c136afb75f4..9782bca473c 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -50,18 +50,20 @@ xt::xtensor source_frac; // Non-member functions //============================================================================== -void calculate_generation_keff(xt::xtensor_fixed> gt, double& keff_generation, vector& k_generation) +void calculate_generation_keff() { + const auto& gt = simulation::global_tallies; + // Get keff for this generation by subtracting off the starting value - keff_generation = + simulation::keff_generation = gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) - - keff_generation; + simulation::keff_generation; double keff_reduced; #ifdef OPENMC_MPI if (settings::solver_type != SolverType::RANDOM_RAY) { // Combine values across all processors - MPI_Allreduce(&keff_generation, &keff_reduced, 1, MPI_DOUBLE, + MPI_Allreduce(&simulation::keff_generation, &keff_reduced, 1, MPI_DOUBLE, MPI_SUM, mpi::intracomm); } else { // If using random ray, MPI parallelism is provided by domain replication. @@ -69,10 +71,10 @@ void calculate_generation_keff(xt::xtensor_fixed& k_generation, std::array& k_sum) +void calculate_average_keff() { // Determine overall generation and number of active generations int i = overall_generation() - 1; @@ -380,14 +382,15 @@ void calculate_average_keff(double& keff, double& keff_std, const vector if (n <= 0) { // For inactive generations, use current generation k as estimate for next // generation - keff = k_generation[i]; + simulation::keff = simulation::k_generation[i]; } else { // Sample mean of keff - k_sum[0] += k_generation[i]; - k_sum[1] += std::pow(k_generation[i], 2); + simulation::k_sum[0] += simulation::k_generation[i]; + simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2); // Determine mean - keff = k_sum[0] / n; + simulation::keff = simulation::k_sum[0] / n; + if (n > 1) { double t_value; if (settings::confidence_intervals) { @@ -399,10 +402,10 @@ void calculate_average_keff(double& keff, double& keff_std, const vector } // Standard deviation of the sample mean of k - keff_std = + simulation::keff_std = t_value * std::sqrt( - (k_sum[1] / n - std::pow(keff, 2)) / (n - 1)); + (simulation::k_sum[1] / n - std::pow(simulation::keff, 2)) / (n - 1)); // In some cases (such as an infinite medium problem), random ray // may estimate k exactly and in an unvarying manner between iterations. @@ -410,8 +413,8 @@ void calculate_average_keff(double& keff, double& keff_std, const vector // power operations may cause an extremely small negative value to occur // inside the sqrt operation, leading to NaN. If this occurs, we check for // it and set the std dev to zero. - if (!std::isfinite(keff_std)) { - keff_std = 0.0; + if (!std::isfinite(simulation::keff_std)) { + simulation::keff_std = 0.0; } } } diff --git a/src/simulation.cpp b/src/simulation.cpp index 2027b60650d..57940f10e8e 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -614,13 +614,13 @@ void finalize_generation() } // Collect results and statistics - calculate_generation_keff(simulation::global_tallies, simulation::keff_generation, simulation::k_generation); - calculate_average_keff(simulation::keff, simulation::keff_std, simulation::k_generation, simulation::k_sum); + calculate_generation_keff(); + calculate_average_keff(); if (settings::calculate_subcritical_k) { // Compute kq and average - calculate_generation_keff(simulation::global_tallies_first_gen, simulation::kq_generation_val, simulation::kq_generation); - calculate_average_keff(simulation::kq, simulation::kq_std, simulation::kq_generation, simulation::kq_sum); + calculate_generation_kq(); + calculate_average_kq(); // Calculate ks from k, kq simulation::ks_generation_val = calculate_ks( diff --git a/src/subcritical.cpp b/src/subcritical.cpp index c7627a7bc59..43c84f09c2b 100644 --- a/src/subcritical.cpp +++ b/src/subcritical.cpp @@ -1,8 +1,10 @@ #include "openmc/subcritical.h" -#include "hdf5.h" #include "openmc/eigenvalue.h" +#include "openmc/math_functions.h" #include "openmc/simulation.h" #include "openmc/tallies/tally.h" + +#include "hdf5.h" #include #include @@ -22,6 +24,94 @@ double convert_to_subcritical_k(double k) { return k_sub; } +void calculate_generation_kq() +{ + const auto& gt = simulation::global_tallies_first_gen; + + // Get keff for this generation by subtracting off the starting value + simulation::kq_generation_val = + gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) - + simulation::kq_generation_val; + + double kq_reduced; +#ifdef OPENMC_MPI + if (settings::solver_type != SolverType::RANDOM_RAY) { + // Combine values across all processors + MPI_Allreduce(&simulation::kq_generation_val, &kq_reduced, 1, MPI_DOUBLE, + MPI_SUM, mpi::intracomm); + } else { + // If using random ray, MPI parallelism is provided by domain replication. + // As such, all fluxes will be reduced at the end of each transport sweep, + // such that all ranks have identical scalar flux vectors, and will all + // independently compute the same value of k. Thus, there is no need to + // perform any additional MPI reduction here. + kq_reduced = simulation::kq_generation_val; + } +#else + kq_reduced = simulation::kq_generation_val; +#endif + + // Normalize single batch estimate of k + // TODO: This should be normalized by total_weight, not by n_particles + if (settings::solver_type != SolverType::RANDOM_RAY) { + kq_reduced /= settings::n_particles; + } + + simulation::kq_generation.push_back(kq_reduced); +} + +void calculate_average_kq() +{ + // Determine overall generation and number of active generations + int i = overall_generation() - 1; + int n; + if (simulation::current_batch > settings::n_inactive) { + n = settings::gen_per_batch * simulation::n_realizations + + simulation::current_gen; + } else { + n = 0; + } + + if (n <= 0) { + // For inactive generations, use current generation k as estimate for next + // generation + simulation::kq = simulation::kq_generation[i]; + } else { + // Sample mean of keff + simulation::kq_sum[0] += simulation::kq_generation[i]; + simulation::kq_sum[1] += std::pow(simulation::kq_generation[i], 2); + // Determine mean + simulation::kq = simulation::kq_sum[0] / n; + + if (n > 1) { + double t_value; + if (settings::confidence_intervals) { + // Calculate t-value for confidence intervals + double alpha = 1.0 - CONFIDENCE_LEVEL; + t_value = t_percentile(1.0 - alpha / 2.0, n - 1); + } else { + t_value = 1.0; + } + + // Standard deviation of the sample mean of k + simulation::keff_std = + t_value * + std::sqrt( + (simulation::kq_sum[1] / n - std::pow(simulation::kq, 2)) / (n - 1)); + + // In some cases (such as an infinite medium problem), random ray + // may estimate k exactly and in an unvarying manner between iterations. + // In this case, the floating point roundoff between the division and the + // power operations may cause an extremely small negative value to occur + // inside the sqrt operation, leading to NaN. If this occurs, we check for + // it and set the std dev to zero. + if (!std::isfinite(simulation::kq_std)) { + simulation::kq_std = 0.0; + } + } + } +} + int openmc_get_subcritical_kq(double* k_combined) { k_combined[0] = 0.0; From 1212e508b747a15a8750edc6185ce11967747532 Mon Sep 17 00:00:00 2001 From: Matthew Louis Date: Tue, 23 Dec 2025 12:31:48 -0500 Subject: [PATCH 8/8] Subcritical k calculation documentation --- docs/source/methods/index.rst | 1 + .../methods/subcritical_multiplication.rst | 103 ++++++++++++++++++ docs/source/usersguide/settings.rst | 17 +++ 3 files changed, 121 insertions(+) create mode 100644 docs/source/methods/subcritical_multiplication.rst diff --git a/docs/source/methods/index.rst b/docs/source/methods/index.rst index 121d04b1ded..c761e461308 100644 --- a/docs/source/methods/index.rst +++ b/docs/source/methods/index.rst @@ -17,6 +17,7 @@ Theory and Methodology charged_particles_physics tallies eigenvalue + subcritical_multiplication depletion energy_deposition parallelization diff --git a/docs/source/methods/subcritical_multiplication.rst b/docs/source/methods/subcritical_multiplication.rst new file mode 100644 index 00000000000..7c0eb50cec3 --- /dev/null +++ b/docs/source/methods/subcritical_multiplication.rst @@ -0,0 +1,103 @@ +.. _methods_subcritical-multiplication: + +======================================= +Subcritical Multiplication Calculations +======================================= + +Subcritical multiplication problems are fixed source simulations with a large +amount of multiplication that are still subcritical with respect to +:math:`k_{eff}`. These problems are common in the design of accelerator driven +systems (ADS) which use an accelerator source to drive a subcritical +multiplication reaction in, e.g., spent fuel to transmute nuclear waste and +even generate power [Bowman]_. An ADS is inherently safe and allows for a much +more flexible fuel composition, hence their popularity for waste transmutation. + +For ADS's, the production of fission neutrons is central as these produced +neutrons amplify the external source. For the case of a proton spallation ADS, +source neutrons are produced from spallation reactions in a heavy metal target +bombarded by high-energy protons. These source neutrons then induce fission in +a subcritical core, producing additional neutrons. + + +.. _methods_subcritical-multiplication-factors: + +---------------------------------- +Subcritical Multiplication Factors +---------------------------------- +In a fixed source simulation, +the total neutron production (per source particle) is used to define + +.. math:: + :label: integral_multiplicity + + \begin{align*} + M - 1 &= \frac{1}{N_{source}}\int d\boldsymbol{r} \int d\boldsymbol{\Omega} \int dE \nu \Sigma_f(\boldsymbol{r}, E) \psi(\boldsymbol{r}, \boldsymbol{\Omega}, E)\\ + &\coloneqq k + k^2 + k^3 + ... \\ + &= \frac{k}{1-k}\\ + \implies M &= \frac{1}{1-k} + \end{align*} + +Where :math:`M` is the subcritical multiplicity, and :math:`k` the subcritical +multiplication factor. The identification on the second line comes from the +picture of a single source neutron inducing several generations of fission +neutrons, producing on average :math:`k` neutrons in the first generation, +which in turn produce :math:`k^2` neutrons in the second generation, and so on. +However, the above picture cannot be taken literally, because the neutrons +born from the external source will have a different importance to neutron +production than will fission neutrons, and we have the following alternative +picture for :math:`M` [Kobayashi]_: + +.. math:: + :label: subcritical_k_factors + + \begin{align*} + M-1 &= k_q + k_q k_s + k_q k_s^2 + ... \\ + &= \frac{k_q}{1-k_s} \\ + \implies \frac{k}{1-k} &= \frac{k_q}{1-k_s}\\ + k &= \frac{k_q}{1 - k_s + k_q} + \end{align*} + +Where :math:`k_q` is the multiplication factor of source neutrons, and :math:`k_s` +is the multiplication factor of fission neutrons, which together define an overall +subcritical multiplication factor :math:`k`. From the above it is clear that +:math:`k < 1 \iff k_s < 1`, and :math:`k_s <1` for :math:`k_{eff}<1`, so a +subcritical system will correctly have :math:`k < 1`. It is however not the case +that :math:`k_s = k_{eff}`, because the angular flux of fission neutrons is not +necessarily the fundamental mode calculated in eigenvalue mode, nor is it +true that :math:`k = k_{eff}`. In fact, for deeply subcritical systems, +:math:`k_{eff}` generally underestimates :math:`k` [Forget]_. In addition, we may +have :math:`k_q>1` despite :math:`k_s<1`, and in the limit, we may have :math:`k` +arbitrarily close to 1: an arbitrarily multiplying system that is still +subcritical with respect to fission neutrons. In fact, this is a primary design +consideration of ADS's, where :math:`k_s` is fixed :math:`<1` to ensure +subcritciality, while :math:`k_q` is maximized to achieve optimal multiplication. +It is therefore necessary to perform fixed source simulations to accurately determine +subcritical multiplication and the flux distribution in ADS's. + +.. _methods_subcritical-multiplication-estimating: + +----------------------------------------------------------------------- +Estimating :math:`k`, :math:`k_q`, and :math:`k_s` in Fixed Source Mode +----------------------------------------------------------------------- +The total multiplication factor :math:`k` can be estimated through :eq:`integral_multiplicity`. +The total fission production can be tallied and estiamted using standard collision, absorption, +and track-length estimators over a neutron history, giving :math:`M-1`, which can be used to +compute :math:`k`. + +To estimate :math:`k_q`, we may use its interpretation interpretation as the multiplication +factor of source neutrons. For a given source history, we may tally the neutron production +estimators, and simply stop before simulating any of the secondary fission neutrons. This gives +an estimate of the neutron production due to source neutrons alone, which can be used to compute +:math:`k_q`. :math:`k_s` can then be computed from :math:`k` and :math:`k_q` using +:eq:`subcritical_k_factors`. + +.. [Bowman] Bowman, Charles D. "Accelerator-driven systems for nuclear waste + transmutation." *Annual Review of Nuclear and Particle Science* 48.1 (1998): + 505-556. + +.. [Kobayashi] Kobayashi, Keisuke, and Kenji Nishihara. "Definition of + subcriticality using the importance function for the production of fission + neutrons." *Nuclear science and engineering* 136.2 (2000): 272-281. + +.. [Forget] Forget, Benoit. "An Efficient Subcritical Multiplication Mode for + Monte Carlo Solvers." *Nuclear Science and Engineering* (2025): 1-11. \ No newline at end of file diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index f973f146552..eff8d125e17 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -845,3 +845,20 @@ the inputs were not modified before the restart run), no particles will be transported and OpenMC will exit immediately. .. note:: A statepoint file must match the input model to be successfully used in a restart simulation. + +---------------------------------------------- +Calculating Subcritical Multiplication Factors +---------------------------------------------- +In addition to the standard effective multiplication factor, :math:`k_{eff}`, calculated using eigenvalue +mode, OpenMC can also calculate subcritical multiplication factors: :math:`k`, :math:`k_q`, and :math:`k_s` +in fixed source mode. Please see :ref:`methods_subcritical-multiplication` for theoretical details on +subcritical multiplication factors. To enable the calculation of subcritical multiplication factors in +fixed source mode, set + +.. code-block:: python + + settings.calculate_subcritical_k = True + +this will enable printing of :math:`k`, :math:`k_q`, and :math:`k_s` both generation-wise and cumulative +averages throughout the simulation analogous to eigenvalue mode. The generation statistics as well as +combined estimates for :math:`k`, :math:`k_q`, and :math:`k_s` will also be stored in the statepoint file. \ No newline at end of file