openmc-dev · mlouis9 · Jul 10, 2025 · Jul 10, 2025 · Jul 10, 2025 · Jul 10, 2025
diff --git 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
@@ -0,0 +1,121 @@
+.. _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.
+
+Source Effectiveness
+--------------------
+In addition, one may quantify the relative effectiveness of a source by the amount of
+multiplication it induces relative to the fundamental mode of the system. This is given
+by the source effectiveness factor :math:`g^*` (or :math:`\phi^*`), defined as [Ye]_:
+
+.. math::
+   M = \frac{g^*}{1 - k_{eff}}
+
+which corrects for the error in the multiplicity one would get from using the
+fundamental mode multiplication factor. As :math:`k_{eff} \to 1`, :math:`g^* \to 1`, 
+since as the system approaches criticality, the importance of the source diminishes
+relative to subsequent fisison neutrons.
+
+.. _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 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.
+
+.. [Ye] Ye, Bin, Chao-Wen Yang, and Chun Zheng. "Measurement of k eff by delayed 
+   neutron multiplication in subcritical systems." Nuclear Science and 
+   Techniques 29.2 (2018): 29.
diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst
@@ -851,3 +851,39 @@ 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 or eigenvalue mode, but note these values are only strictly well-defined for fixed source simulations. 
+Please see :ref:`methods_subcritical-multiplication` for theoretical details on
+subcritical multiplication factors. 
+
+Fixed Source Mode
+-----------------
+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.
+To print out all of the subcritical multiplication factors during the transport calculation, set
+
+.. code-block:: python
+
+    settings.print_all_k_factors = True
+
+Eigenvalue Mode
+---------------
+These quantities may also be computed using the eigensolver if `settings.source` is set to the source
+distribution for your corresponding fixed source problem. These quantities are computed by default, and there's not need
+to set `settings.calculate_subcritical_k` to True. The generation statistics as well as combined estimates for :math:`k`,
+:math:`k_q`, and :math:`k_s` will be stored in the statepoint file. In addition, the source effectiveness factor will also
+be computed, since :math:`k_{eff}` is available. Please interpret these values with caution, as they are only strictly 
+well-defined when :math:`k_{eff} < 1`.
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
@@ -330,8 +330,15 @@ enum TallyScore {
 };
 
 // Global tally parameters
-constexpr int N_GLOBAL_TALLIES {4};
-enum class GlobalTally { K_COLLISION, K_ABSORPTION, K_TRACKLENGTH, LEAKAGE };
+enum class GlobalTally {
+  K_COLLISION,
+  K_ABSORPTION,
+  K_TRACKLENGTH,
+  LEAKAGE,
+  K_TRACKLENGTH_SQ, // for calculation of stddev for generational k
+  SIZE
+};
+constexpr int N_GLOBAL_TALLIES {static_cast<int>(GlobalTally::SIZE)};
 
 // Miscellaneous
 constexpr int C_NONE {-1};
@@ -361,6 +368,9 @@ enum class RunMode {
   VOLUME
 };
 
+// Eigenvalue calculation parameters
+enum class KeffType { k, kq, ks };
+
 enum class SolverType { MONTE_CARLO, RANDOM_RAY };
 
 enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID };

diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h
@@ -21,8 +21,12 @@ namespace openmc {
 
 namespace simulation {
 
-extern double keff_generation; //!<  Single-generation k on each processor
+extern array<double, 2>
+  keff_generation; //!<  Single-generation k on each processor
+extern array<double, 2>
+  kq_generation_val;           //!<  Single-generation kq 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
 
@@ -32,15 +36,25 @@ extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
 // Non-member functions
 //==============================================================================
 
-//! Collect/normalize the tracklength keff from each process
+//! Collect/normalize the tracklength keff from each process (keff)
 void calculate_generation_keff();
 
+//! Collect/normalize the tracklength keff from each process
+//!
+//! Depending on KeffType, this function will calculate either the kq or ks
+//! generation keff
+void calculate_generation_keff(KeffType type);
+
+std::pair<double, double> convert_m_to_k(double m, double m_std);
+std::pair<double, double> convert_k_to_m(double k, double k_std);
+
 //! 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(KeffType type);
 
 //! Calculates a minimum variance estimate of k-effective
 //!
@@ -55,6 +69,11 @@ void calculate_average_keff();
 //! \param[out] k_combined Estimate of k-effective and its standard deviation
 //! \return Error status
 extern "C" int openmc_get_keff(double* k_combined);
+extern "C" int openmc_get_kq(double* kq_combined);
+int get_combined_k_from_tallies(double* k_combined,
+  std::array<double, 3>& k_combined_weights, double keff, double keff_std,
+  xt::xtensor_fixed<double, xt::xshape<N_GLOBAL_TALLIES, 3>> global_tallies,
+  double k_col_abs, double k_col_tra, double k_abs_tra);
 
 //! Sample/redistribute source sites from accumulated fission sites
 void synchronize_bank();

diff --git a/include/openmc/output.h b/include/openmc/output.h
@@ -5,6 +5,7 @@
 #define OPENMC_OUTPUT_H
 
 #include <string>
+#include <xtensor/xsort.hpp>
 
 #include "openmc/particle.h"
 
@@ -48,6 +49,8 @@ void print_columns();
 
 //! Display information about a generation of neutrons
 void print_generation();
+void print_generation_values(int idx, int n, std::string batch_and_gen,
+  array<double, 2> k_gen, double k, double k_std);
 
 //! Display time elapsed for various stages of a run
 void print_runtime();

diff --git a/include/openmc/settings.h b/include/openmc/settings.h
@@ -169,7 +169,10 @@ 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 bool
+  print_all_k_factors; //!< Print k-effective for all batches and generations
 extern int
   max_history_splits; //!< maximum number of particle splits for weight windows
 extern int max_secondaries;       //!< maximum number of secondaries in the bank

diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h
@@ -4,11 +4,13 @@
 #ifndef OPENMC_SIMULATION_H
 #define OPENMC_SIMULATION_H
 
+#include "openmc/array.h"
 #include "openmc/mesh.h"
 #include "openmc/particle.h"
 #include "openmc/vector.h"
 
 #include <cstdint>
+#include <sys/types.h>
 
 namespace openmc {
 
@@ -28,13 +30,30 @@ 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
+  k; //!< average k over batches, used for subcritical multiplication problems
+extern "C" double k_std; //!< standard deviation of average k, used for
+                         //!< subcritical multiplication problems
+extern "C" double
+  kq; //!< average kq over batches, used for subcritical multiplication problems
+extern "C" double kq_std; //!< standard deviation of average kq, used for
+                          //!< subcritical multiplication problems
+extern "C" double
+  ks; //!< average ks over batches, used for subcritical multiplication problems
+extern "C" double ks_std;    //!< standard deviation of average ks, used for
+                             //!< subcritical multiplication problems
 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 double log_spacing; //!< lethargy spacing for energy grid searches
-extern "C" int n_lost_particles;   //!< cumulative number of lost particles
+  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;
+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?
 extern "C" int restart_batch;      //!< batch at which a restart job resumed
 extern "C" bool satisfy_triggers;  //!< have tally triggers been satisfied?
@@ -46,9 +65,25 @@ extern int64_t work_per_rank;      //!< number of particles per MPI rank
 extern const RegularMesh* entropy_mesh;
 extern const RegularMesh* ufs_mesh;
 
-extern vector<double> k_generation;
+extern vector<array<double, 2>> k_generation;
+extern vector<array<double, 2>> kq_generation;
+extern vector<array<double, 2>> ks_generation;
 extern vector<int64_t> work_index;
 
+enum class KEstimator : int {
+  K_COLLISION = 0,
+  K_ABSORPTION = 1,
+  K_TRACKLENGTH = 2,
+  N_ESTIMATORS
+};
+
+constexpr int N_K_EST = static_cast<int>(KEstimator::N_ESTIMATORS);
+extern std::array<std::array<double, N_K_EST>, N_K_EST> k_kq_products;
+extern std::array<std::array<double, N_K_EST>, N_K_EST> k_kq_product;
+
+extern std::array<double, 3> k_combined_weights;
+extern std::array<double, 3> kq_combined_weights;
+
 } // namespace simulation
 
 //==============================================================================

diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h
@@ -49,6 +49,8 @@ void write_source_bank(hid_t group_id, span<SourceSite> source_bank,
 
 void read_source_bank(
   hid_t group_id, vector<SourceSite>& sites, bool distribute);
+void write_global_tallies(hid_t file_id);
+void read_global_tallies(hid_t file_id);
 void write_tally_results_nr(hid_t file_id);
 void restart_set_keff();
 void write_unstructured_mesh_results();

diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
@@ -222,6 +222,9 @@ namespace simulation {
 //! Global tallies (such as k-effective estimators)
 extern xt::xtensor_fixed<double, xt::xshape<N_GLOBAL_TALLIES, 3>>
   global_tallies;
+//! Global tallies for first generation (for subcritical multiplication)
+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;
@@ -230,8 +233,14 @@ extern "C" int32_t n_realizations;
 extern double global_tally_absorption;
 extern double global_tally_collision;
 extern double global_tally_tracklength;
+extern double global_tally_tracklength_sq;
 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;
+extern double global_tally_tracklength_sq_first_gen;
+
 //==============================================================================
 // Non-member functions
 //==============================================================================