Fixed-Source Subcritical Multiplication Factors

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

docs/source/methods/index.rst

@@ -17,6 +17,7 @@ Theory and Methodology
     charged_particles_physics
     tallies
     eigenvalue
+    subcritical_multiplication
     depletion
     energy_deposition
     parallelization

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.

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.

include/openmc/eigenvalue.h

@@ -1,3 +1,3 @@
 //! \file eigenvalue.h
 //! \brief Data/functions related to k-eigenvalue calculations
 
@@ -22,7 +22,10 @@
 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<double, 2> k_sum; //!< Used to reduce sum and sum_sq
+extern array<double, 2> kq_sum;
 extern vector<double> entropy; //!< Shannon entropy at each generation
 extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
 

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

include/openmc/simulation.h

@@ -1,3 +1,3 @@
 //! \file simulation.h
 //! \brief Variables/functions related to a running simulation
 
@@ -28,11 +28,21 @@
 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* ufs_mesh;
 
 extern vector<double> k_generation;
+extern vector<double> kq_generation;
+extern vector<double> ks_generation;
 extern vector<int64_t> work_index;
 
 } // namespace simulation

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 <utility>
+
+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

include/openmc/tallies/tally.h

@@ -222,6 +222,8 @@ namespace simulation {
 //! Global tallies (such as k-effective estimators)
 extern xt::xtensor_fixed<double, xt::xshape<N_GLOBAL_TALLIES, 3>>
   global_tallies;
+extern xt::xtensor_fixed<double, xt::xshape<N_GLOBAL_TALLIES, 3>>
+  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
 //==============================================================================

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)