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/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 diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index b456fee21eb..1b1c2432bc0 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 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/include/openmc/simulation.h b/include/openmc/simulation.h index 9a6cf1b2131..60cfe4d746c 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -28,11 +28,21 @@ 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 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? @@ -47,6 +57,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..1479e82ba46 --- /dev/null +++ b/include/openmc/subcritical.h @@ -0,0 +1,29 @@ +//! \file subcritical.h +//! \brief Data/functions related to subcritical (fixed source) multiplication calculations + +#ifndef OPENMC_SUBCRITICAL_H +#define OPENMC_SUBCRITICAL_H + +#include "hdf5.h" +#include + +namespace openmc { + +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); + +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/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/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/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/eigenvalue.cpp b/src/eigenvalue.cpp index 8412cbd3b47..9782bca473c 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; diff --git a/src/finalize.cpp b/src/finalize.cpp index 344eaa1a0a7..6073dd65f16 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; @@ -205,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/output.cpp b/src/output.cpp index 80e2b10ab89..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,14 +397,20 @@ 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) { + if (settings::entropy_on && !settings::calculate_subcritical_k) { fmt::print(" {:8.5f}", simulation::entropy[idx]); } 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); @@ -534,21 +543,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 +576,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..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); @@ -646,7 +656,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 +668,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 5b472468fc4..aec3bf5cab9 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,20 @@ 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) { + 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 { + 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/src/simulation.cpp b/src/simulation.cpp index b536ae5881f..57940f10e8e 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 @@ -142,6 +143,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); } @@ -307,6 +311,13 @@ 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; +double ks_std; double log_spacing; int n_lost_particles {0}; bool need_depletion_rx {false}; @@ -321,6 +332,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 @@ -366,7 +379,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); + } } } @@ -414,6 +429,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; } @@ -505,7 +521,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); @@ -516,29 +534,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; - // Update global tallies with the accumulation variables if (settings::run_mode == RunMode::EIGENVALUE) { - gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision; + 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; } + if (settings::run_mode == RunMode::FIXED_SOURCE && + 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); + + // 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 - if (settings::run_mode == RunMode::EIGENVALUE) { + // Reset tallies + if (settings::run_mode == RunMode::EIGENVALUE || + (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 && @@ -552,17 +602,35 @@ 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(); + if (settings::calculate_subcritical_k) { + // Compute kq and average + calculate_generation_kq(); + calculate_average_kq(); + + // 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 = 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 if (mpi::master && settings::verbosity >= 7) { print_generation(); @@ -782,6 +850,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 @@ -791,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 @@ -803,6 +880,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()) { @@ -815,6 +893,20 @@ 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) { +// 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(); } p.event_death(); 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 new file mode 100644 index 00000000000..43c84f09c2b --- /dev/null +++ b/src/subcritical.cpp @@ -0,0 +1,318 @@ +#include "openmc/subcritical.h" +#include "openmc/eigenvalue.h" +#include "openmc/math_functions.h" +#include "openmc/simulation.h" +#include "openmc/tallies/tally.h" + +#include "hdf5.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; +} + +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; + 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 6eef1da9cfd..888e9bf72d5 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 //============================================================================== @@ -1071,8 +1076,10 @@ 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) { + 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) / @@ -1085,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 @@ -1093,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; } } diff --git a/tests/regression_tests/fixed_source_k/__init__.py b/tests/regression_tests/fixed_source_k/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/fixed_source_k/inputs_true_1.dat b/tests/regression_tests/fixed_source_k/inputs_true_1.dat new file mode 100644 index 00000000000..b2d62492020 --- /dev/null +++ b/tests/regression_tests/fixed_source_k/inputs_true_1.dat @@ -0,0 +1,39 @@ + + + + + + + + + + + + + + + + + + + + + + + fixed source + 100 + 10 + + + 0 0 0 + + + + + + + + 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..cfd74f2ef29 --- /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..cfd74f2ef29 --- /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 new file mode 100644 index 00000000000..639e6f8cd6d --- /dev/null +++ b/tests/regression_tests/fixed_source_k/test.py @@ -0,0 +1,122 @@ +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 + 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.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(): + openmc.reset_auto_ids() + model = openmc.Model() + + universe = create_universe() + 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.particles = 100 + + model.tallies = openmc.Tallies() + tally = openmc.Tally(name='νΣf tally') + tally.scores = ['nu-fission'] + model.tallies.append(tally) + + 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') + 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() + 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() + 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 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)))