Subcritical Multiplication Run Mode

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,100 @@
+.. _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-methods:
+----------------------------------
+Subcritical Multiplication Methods
+----------------------------------
+The most straightforward method for performing a subcritical multiplication calculation is to perform a fixed source
+simulation. For near critical systems, however, excessive secondary particle production can result in memory issues if
+not properly managed. In very near critical systems, performing an eigenvalue simulation may be sufficient,
+as the flux distribution will be close to the fundamental mode, and :math:`k` will be close to :math:`k_{eff}`, just note that
+tally results must be multiplied by :math:`1/(1-k)`. For more deeply subcritical systems, there are several tailor-made methods
+for subcritical multiplication calculations that limit excessive secondary particle production while equivalent results to fixed source simulations.
+The method currently implemented in OpenMC is that of [Forget]_, which uses a modified power iteration to converge a source
+distribution taking into account the external source and all generations of fission neutrons, which is then sampled from in generations
+to accumulate statistics about the system. This method resembles the standard eigenvalue simulation (see :ref:`methods_eigenvalue`) in every way except 
+that tallies have a pre-defined normalization by :math:`1/(1-k)` (rather than normalized to an arbitrary power distribution), the :math:`k`-eigenvalue is the 
+subcritical multiplication factor (rather than the fundamental mode eigenvalue :math:`k_{eff}`), and the results are determined by the prescribed 
+source distribution.
+
+
+.. [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

@@ -33,6 +33,12 @@ be specified:
   Runs a fixed-source calculation with a specified external source, specified in
   the :attr:`Settings.source` attribute.
 
+'subcritical multiplication'
+  Runs a subcritical multiplication calculation (as described in :ref:`methods_subcritical-multiplication`). 
+  In this mode, the :attr:`Settings.source` specifies an external source distribution that is used for all 
+  generations. Note that tallies must be manually normalized by :math:`1/(1-k)` for direct comparison with fixed
+  source results.
+
 'volume'
   Runs a stochastic volume calculation.
 

include/openmc/constants.h

@@ -356,6 +356,7 @@ enum class RunMode {
   UNSET, // default value, OpenMC throws error if left to this
   FIXED_SOURCE,
   EIGENVALUE,
+  SUBCRITICAL_MULTIPLICATION,
   PLOTTING,
   PARTICLE,
   VOLUME

include/openmc/settings.h

@@ -198,6 +198,8 @@ extern "C" int verbosity;          //!< How verbose to make output
 extern double weight_cutoff;       //!< Weight cutoff for Russian roulette
 extern double weight_survive;      //!< Survival weight after Russian roulette
 
+bool eigenvalue_like();
+
 } // namespace settings
 
 //==============================================================================

openmc/settings.py

@@ -22,6 +22,7 @@
 
 class RunMode(Enum):
     EIGENVALUE = 'eigenvalue'
+    SUBCRITICAL_MULTIPLICATION = 'subcritical multiplication'
     FIXED_SOURCE = 'fixed source'
     PLOT = 'plot'
     VOLUME = 'volume'
@@ -238,7 +239,7 @@ class Settings:
         The 'nuclides' list indicates what nuclides the method should be applied
         to. In its absence, the method will be applied to all nuclides with 0 K
         elastic scattering data present.
-    run_mode : {'eigenvalue', 'fixed source', 'plot', 'volume', 'particle restart'}
+    run_mode : {'eigenvalue', 'fixed source', 'subcritical multiplication', 'plot', 'volume', 'particle restart'}
         The type of calculation to perform (default is 'eigenvalue')
     seed : int
         Seed for the linear congruential pseudorandom number generator

openmc/statepoint.py

@@ -212,7 +212,7 @@ def date_and_time(self):
 
     @property
     def entropy(self):
-        if self.run_mode == 'eigenvalue':
+        if 'entropy' in self._f:
             return self._f['entropy'][()]
         else:
             return None
@@ -233,7 +233,7 @@ def filters(self):
 
     @property
     def generations_per_batch(self):
-        if self.run_mode == 'eigenvalue':
+        if 'generations_per_batch' in self._f:
             return self._f['generations_per_batch'][()]
         else:
             return None
@@ -268,14 +268,14 @@ def k_cmfd(self):
 
     @property
     def k_generation(self):
-        if self.run_mode == 'eigenvalue':
+        if 'k_generation' in self._f:
             return self._f['k_generation'][()]
         else:
             return None
 
     @property
     def keff(self):
-        if self.run_mode == 'eigenvalue':
+        if 'k_combined' in self._f:
             return ufloat(*self._f['k_combined'][()])
         else:
             return None
@@ -290,21 +290,21 @@ def k_combined(self):
 
     @property
     def k_col_abs(self):
-        if self.run_mode == 'eigenvalue':
+        if 'k_col_abs' in self._f:
             return self._f['k_col_abs'][()]
         else:
             return None
 
     @property
     def k_col_tra(self):
-        if self.run_mode == 'eigenvalue':
+        if 'k_col_tra' in self._f:
             return self._f['k_col_tra'][()]
         else:
             return None
 
     @property
     def k_abs_tra(self):
-        if self.run_mode == 'eigenvalue':
+        if 'k_abs_tra' in self._f:
             return self._f['k_abs_tra'][()]
         else:
             return None
@@ -329,7 +329,7 @@ def n_batches(self):
 
     @property
     def n_inactive(self):
-        if self.run_mode == 'eigenvalue':
+        if 'n_inactive' in self._f:
             return self._f['n_inactive'][()]
         else:
             return None

src/main.cpp

@@ -42,6 +42,13 @@ int main(int argc, char* argv[])
       break;
     }
     break;
+  case RunMode::SUBCRITICAL_MULTIPLICATION:
+    switch (settings::solver_type) {
+    case SolverType::MONTE_CARLO:
+      err = openmc_run();
+      break;
+    }
+    break;
   case RunMode::PLOTTING:
     err = openmc_plot_geometry();
     break;

src/output.cpp

@@ -444,11 +444,11 @@ void print_runtime()
     show_time("Collisions", time_event_collision.elapsed(), 2);
     show_time("Particle death", time_event_death.elapsed(), 2);
   }
-  if (settings::run_mode == RunMode::EIGENVALUE) {
+  if (settings::eigenvalue_like()) {
     show_time("Time in inactive batches", time_inactive.elapsed(), 1);
   }
   show_time("Time in active batches", time_active.elapsed(), 1);
-  if (settings::run_mode == RunMode::EIGENVALUE) {
+  if (settings::eigenvalue_like()) {
     show_time("Time synchronizing fission bank", time_bank.elapsed(), 1);
     show_time("Sampling source sites", time_bank_sample.elapsed(), 2);
     show_time("SEND/RECV source sites", time_bank_sendrecv.elapsed(), 2);
@@ -535,7 +535,7 @@ void print_results()
   const auto& gt = simulation::global_tallies;
   double mean, stdev;
   if (n > 1) {
-    if (settings::run_mode == RunMode::EIGENVALUE) {
+    if (settings::eigenvalue_like()) {
       std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_COLLISION, 0), n);
       fmt::print(" k-effective (Collision)     = {:.5f} +/- {:.5f}\n", mean,
         t_n1 * stdev);
@@ -560,7 +560,7 @@ void print_results()
       warning("Could not compute uncertainties -- only one "
               "active batch simulated!");
 
-    if (settings::run_mode == RunMode::EIGENVALUE) {
+    if (settings::eigenvalue_like()) {
       fmt::print(" k-effective (Collision)    = {:.5f}\n",
         gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n);
       fmt::print(" k-effective (Track-length) = {:.5f}\n",

src/particle.cpp

@@ -279,7 +279,7 @@ void Particle::event_advance()
   }
 
   // Score track-length estimate of k-eff
-  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
+  if (settings::eigenvalue_like() && type().is_neutron()) {
     keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
   }
 
@@ -341,7 +341,7 @@ void Particle::event_collide()
 {
 
   // Score collision estimate of keff
-  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
+  if (settings::eigenvalue_like() && type().is_neutron()) {
     keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
   }
 
@@ -518,7 +518,7 @@ void Particle::event_death()
 
   // Record the number of progeny created by this particle.
   // This data will be used to efficiently sort the fission bank.
-  if (settings::run_mode == RunMode::EIGENVALUE) {
+  if (settings::eigenvalue_like()) {
     int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
     simulation::progeny_per_particle[offset] = n_progeny();
   }
@@ -840,7 +840,7 @@ void Particle::write_restart() const
     write_dataset(file_id, "type", type().pdg_number());
 
     int64_t i = current_work();
-    if (settings::run_mode == RunMode::EIGENVALUE) {
+    if (settings::eigenvalue_like()) {
       // take source data from primary bank for eigenvalue simulation
       write_dataset(file_id, "weight", simulation::source_bank[i - 1].wgt);
       write_dataset(file_id, "energy", simulation::source_bank[i - 1].E);

src/particle_restart.cpp

@@ -51,6 +51,8 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode)
     previous_run_mode = RunMode::EIGENVALUE;
   } else if (mode == "fixed source") {
     previous_run_mode = RunMode::FIXED_SOURCE;
+  } else if (mode == "subcritical multiplication") {
+    previous_run_mode = RunMode::SUBCRITICAL_MULTIPLICATION;
   }
   read_dataset(file_id, "id", p.id());
   int type;
@@ -110,6 +112,7 @@ void run_particle_restart()
   int64_t particle_seed;
   switch (previous_run_mode) {
   case RunMode::EIGENVALUE:
+  case RunMode::SUBCRITICAL_MULTIPLICATION:
   case RunMode::FIXED_SOURCE:
     particle_seed = (simulation::total_gen + overall_generation() - 1) *
                       settings::n_particles +

src/physics.cpp

@@ -110,7 +110,7 @@ void sample_neutron_reaction(Particle& p)
 
   if (nuc->fissionable_ && p.neutron_xs(i_nuclide).fission > 0.0) {
     auto& rx = sample_fission(i_nuclide, p);
-    if (settings::run_mode == RunMode::EIGENVALUE) {
+    if (settings::eigenvalue_like()) {
       create_fission_sites(p, i_nuclide, rx);
     } else if (settings::run_mode == RunMode::FIXED_SOURCE &&
                settings::create_fission_neutrons) {
@@ -201,7 +201,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
 
   // Determine whether to place fission sites into the shared fission bank
   // or the secondary particle bank.
-  bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);
+  bool use_fission_bank = (settings::eigenvalue_like());
 
   // Counter for the number of fission sites successfully stored to the shared
   // fission bank or the secondary particle bank
@@ -649,7 +649,7 @@ 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::eigenvalue_like()) {
       p.keff_tally_absorption() += wgt_absorb *
                                    p.neutron_xs(i_nuclide).nu_fission /
                                    p.neutron_xs(i_nuclide).absorption;
@@ -659,7 +659,7 @@ 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::eigenvalue_like()) {
         p.keff_tally_absorption() += p.wgt() *
                                      p.neutron_xs(i_nuclide).nu_fission /
                                      p.neutron_xs(i_nuclide).absorption;
@@ -1212,7 +1212,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide)
     // Stedry, "Self-consistent energy normalization for quasistatic reactor
     // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020.
     double wgt = photon_wgt;
-    if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) {
+    if (settings::eigenvalue_like() && !is_fission(rx->mt_)) {
       wgt *= simulation::keff;
     }
 

src/physics_mg.cpp

@@ -49,7 +49,7 @@ void sample_reaction(Particle& p)
   // absorption (including fission)
 
   if (model::materials[p.material()]->fissionable()) {
-    if (settings::run_mode == RunMode::EIGENVALUE ||
+    if (settings::eigenvalue_like() ||
         (settings::run_mode == RunMode::FIXED_SOURCE &&
           settings::create_fission_neutrons)) {
       create_fission_sites(p);
@@ -127,7 +127,7 @@ void create_fission_sites(Particle& p)
 
   // Determine whether to place fission sites into the shared fission bank
   // or the secondary particle bank.
-  bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);
+  bool use_fission_bank = (settings::eigenvalue_like());
 
   // Counter for the number of fission sites successfully stored to the shared
   // fission bank or the secondary particle bank