Kinetics parameters using Iterated Fission Probability

CMakeLists.txt

@@ -343,6 +343,7 @@ list(APPEND libopenmc_SOURCES
   src/geometry.cpp
   src/geometry_aux.cpp
   src/hdf5_interface.cpp
+  src/ifp.cpp
   src/initialize.cpp
   src/lattice.cpp
   src/material.cpp

docs/source/usersguide/index.rst

@@ -22,6 +22,7 @@ essential aspects of using OpenMC to perform simulations.
     plots
     depletion
     decay_sources
+    kinetics
     scripts
     processing
     parallel

docs/source/usersguide/kinetics.rst

@@ -0,0 +1,110 @@
+.. _kinetics:
+
+===================
+Kinetics parameters
+===================
+
+OpenMC has the capability to estimate the following adjoint-weighted effective
+generation time :math:`\Lambda_{\text{eff}}` and the effective delayed neutron
+fraction :math:`\beta_{\text{eff}}`. These parameters are calculated using the
+iterated fission probability (IFP) method [Hurwitz_1964]_ based on a similar
+approach as in `Serpent 2 <https://doi.org/10.1016/j.anucene.2013.10.032>`_. The
+implementation in OpenMC is limited to eigenvalue calculations and is described
+in more details in [Dorville_2025]_.
+
+----------------------------------
+Iterated Fission Probability (IFP)
+----------------------------------
+
+With IFP, additional information needs to be recorded during the simulation
+compared to a typical eigenvalue calculation. OpenMC stores an additional
+set of values (neutron lifetime or delayed neutron group number for
+:math:`\Lambda_{\text{eff}}` or :math:`\beta_{\text{eff}}`, respectively)
+for every fission neutron simulated. Each set of values corresponds to
+the values that are associated to the :math:`N_{\text{gen}}` direct ancestors
+of any given fission neutron.
+
+:math:`N_{\text{gen}}` is referred to as the number of generations in the
+IFP method and corresponds to the number of generations between the birth of
+a fission neutron and the time its score is added to the IFP tally. By default,
+OpenMC considers 10 generations but this value can be modified by the user via
+the ``ifp_n_generation`` settings in the Python API::
+
+    settings.ifp_n_generation = 5
+
+``ifp_n_generation`` should be greater than 0, but should also be lower than
+or equal to the number of inactive batches declared for the calculation.
+The respect of these constraints is verified by OpenMC before any calculation.
+
+OpenMC will automatically detect the type of data that needs to be stored based
+on the tally scores selected by the user. This guarantees that only information
+of interest are stored during a simulation and avoids using extra memory when
+only one parameter is needed. The following table shows the tally scores that
+are needed to compute kinetics parameters in OpenMC:
+
+.. table:: **OpenMC tally scores needed to calculate adjoint-weighted kinetics parameters**
+    :align: center
+
+    =============================== ============================ ========================== ========
+    OpenMC tally score \\ Parameter :math:`\Lambda_{\text{eff}}` :math:`\beta_{\text{eff}}` Both
+    =============================== ============================ ========================== ========
+    ``ifp-time-numerator``          X                                                       X
+    ``ifp-beta-numerator``                                       X                          X
+    ``ifp-denominator``             X                            X                          X
+    =============================== ============================ ========================== ========
+
+|
+
+.. note:: Because the memory footprint of additional data is generally non-negligible
+    with IFP, it is recommended to choose the value for ``ifp_n_generation`` carefully.
+    For example, using one generation for both kinetics parameters corresponds to store
+    one additional integer (for the delayed neutron group number used with
+    :math:`\beta_{\text{eff}}`) and one floating point value (for the neutron lifetime
+    used with :math:`\Lambda_{\text{eff}}`) for every fission neutron simulated once the
+    asymptotic regime is reached.
+
+-----------------------------
+Obtaining kinetics parameters
+-----------------------------
+
+Here is an example showing how to declare the three available IFP scores in a
+single tally::
+
+    tally = openmc.Tally(name="ifp-scores")
+    tally.scores = [
+        "ifp-time-numerator",
+        "ifp-beta-numerator",
+        "ifp-denominator"
+    ]
+
+The effective generation time :math:`\Lambda_{\text{eff}}` is calculated
+by dividing the result of the ``ifp-time-numerator`` score by the one obtained
+for ``ifp-denominator`` and by the :math:`k_{\text{eff}}` of the simulation:
+
+.. math::
+    :label: lambda_eff
+
+    \Lambda_{\text{eff}} = \frac{S_{\text{ifp-time-numerator}}}{S_{\text{ifp-denominator}} \times k_{\text{eff}}}
+
+The effective delayed neutron fraction :math:`\beta_{\text{eff}}` is calculated
+by dividing the result of the ``ifp-beta-numerator`` score by the one obtained
+for ``ifp-denominator``:
+
+.. math::
+    :label: beta_eff
+
+    \beta_{\text{eff}} = \frac{S_{\text{ifp-beta-numerator}}}{S_{\text{ifp-denominator}}}
+
+.. only:: html
+
+   .. rubric:: References
+
+.. [Hurwitz_1964] H. Hurwitz Jr., "Naval Reactors Physics Handbook", volume 1, p. 864.
+    Radkowsky, A. (Ed.), Naval Reactors, Division of Reactor Development, U.S.
+    Atomic Energy Commission (1964).
+
+.. [Dorville_2025] J. Dorville, L. Labrie-Cleary, and P. K. Romano, "Implementation
+    of the Iterated Fission Probability Method in OpenMC to Compute Adjoint-Weighted
+    Kinetics Parameters", International Conference on Mathematics and Computational
+    Methods Applied to Nuclear Science and Engineering (M&C 2025), Denver, April 27-30,
+    2025 (to be presented).

docs/source/usersguide/tallies.rst

@@ -322,6 +322,24 @@ The following tables show all valid scores:
     |                      |particle. Note that this score can only be combined|
     |                      |with a cell filter and an energy filter.           |
     +----------------------+---------------------------------------------------+
+    |ifp-time-numerator    |Adjoint-weighted lifetime of neutron produced by   |
+    |                      |fission in units of seconds per source particle.   |
+    |                      |This score is used to compute kinetics parameters  |
+    |                      |using the iterated fission probability (IFP)       |
+    |                      |method.                                            |
+    +----------------------+---------------------------------------------------+
+    |ifp-beta-numerator    |Adjoint-weighted number of delayed fission events  |
+    |                      |in units of number of delayed fission event per    |
+    |                      |source particle. This score is used to compute     |
+    |                      |kinetics parameters using the iterated fission     |
+    |                      |probability (IFP) method.                          |
+    +----------------------+---------------------------------------------------+
+    |ifp-denominator       |Weights corresponding to the number of fission     |
+    |                      |events in units of number of fission event per     |
+    |                      |source particle. This score is used to compute     |
+    |                      |kinetics parameters using the iterated fission     |
+    |                      |probability (IFP) method.                          |
+    +----------------------+---------------------------------------------------+
 
 .. _usersguide_tally_normalization:
 

include/openmc/bank.h

@@ -22,6 +22,14 @@ extern SharedArray<SourceSite> surf_source_bank;
 
 extern SharedArray<SourceSite> fission_bank;
 
+extern vector<vector<int>> ifp_source_delayed_group_bank;
+
+extern vector<vector<double>> ifp_source_lifetime_bank;
+
+extern vector<vector<int>> ifp_fission_delayed_group_bank;
+
+extern vector<vector<double>> ifp_fission_lifetime_bank;
+
 extern vector<int64_t> progeny_per_particle;
 
 } // namespace simulation

include/openmc/constants.h

@@ -312,7 +312,10 @@ enum TallyScore {
   SCORE_FISS_Q_PROMPT = -14,      // prompt fission Q-value
   SCORE_FISS_Q_RECOV = -15,       // recoverable fission Q-value
   SCORE_DECAY_RATE = -16,         // delayed neutron precursor decay rate
-  SCORE_PULSE_HEIGHT = -17        // pulse-height
+  SCORE_PULSE_HEIGHT = -17,       // pulse-height
+  SCORE_IFP_TIME_NUM = -18,       // IFP lifetime numerator
+  SCORE_IFP_BETA_NUM = -19,       // IFP delayed fraction numerator
+  SCORE_IFP_DENOM = -20           // IFP common denominator
 };
 
 // Global tally parameters
@@ -322,6 +325,9 @@ enum class GlobalTally { K_COLLISION, K_ABSORPTION, K_TRACKLENGTH, LEAKAGE };
 // Miscellaneous
 constexpr int C_NONE {-1};
 
+// Default value of generation for IFP
+constexpr int DEFAULT_IFP_N_GENERATION {10};
+
 // Interpolation rules
 enum class Interpolation {
   histogram = 1,

include/openmc/ifp.h

@@ -0,0 +1,188 @@
+#ifndef OPENMC_IFP_H
+#define OPENMC_IFP_H
+
+#include "openmc/message_passing.h"
+#include "openmc/particle.h"
+#include "openmc/particle_data.h"
+#include "openmc/settings.h"
+
+namespace openmc {
+
+//! Check the value of the IFP parameter for beta effective or both.
+//!
+//! \return true if "BetaEffective" or "Both", false otherwise.
+bool is_beta_effective_or_both();
+
+//! Check the value of the IFP parameter for generation time or both.
+//!
+//! \return true if "GenerationTime" or "Both", false otherwise.
+bool is_generation_time_or_both();
+
+//! Resize IFP vectors
+//!
+//! \param[in,out] delayed_groups List of delayed group numbers
+//! \param[in,out] lifetimes List of lifetimes
+//! \param[in] n  Dimension to resize vectors
+template<typename T, typename U>
+void resize_ifp_data(vector<T>& delayed_groups, vector<U>& lifetimes, int64_t n)
+{
+  if (is_beta_effective_or_both()) {
+    delayed_groups.resize(n);
+  }
+  if (is_generation_time_or_both()) {
+    lifetimes.resize(n);
+  }
+}
+
+//! Update a list of values by adding a new value if the size
+//! of the list can accomodate the new value or by shifting all
+//! values to the left (removing the first value of the list
+//! and adding the new value at the end of the list).
+//!
+//! \param[in] value Value to add to the list
+//! \param[in] data Initial version of the list
+//! \return Updated list
+template<typename T>
+vector<T> _ifp(const T& value, const vector<T>& data)
+{
+  vector<T> updated;
+  size_t source_idx = data.size();
+
+  if (source_idx < settings::ifp_n_generation) {
+    updated.resize(source_idx + 1);
+    for (size_t i = 0; i < source_idx; i++) {
+      updated[i] = data[i];
+    }
+    updated[source_idx] = value;
+  } else if (source_idx == settings::ifp_n_generation) {
+    updated.resize(source_idx);
+    for (size_t i = 0; i < source_idx - 1; i++) {
+      updated[i] = data[i + 1];
+    }
+    updated[source_idx - 1] = value;
+  }
+  return updated;
+}
+
+//! \brief Iterated Fission Probability (IFP) method.
+//!
+//! Add the IFP information in the IFP banks using the same index
+//! as the one used to append the fission site to the fission bank.
+//! Multithreading protection is guaranteed by the index returned by the
+//! thread_safe_append call in physics.cpp.
+//!
+//! Needs to be done after the delayed group is found.
+//!
+//! \param[in] p Particle
+//! \param[in] site Fission site
+//! \param[in] idx Bank index from the thread_safe_append call in physics.cpp
+void ifp(const Particle& p, const SourceSite& site, int64_t idx);
+
+//! Resize the IFP banks used in the simulation
+void resize_simulation_ifp_banks();
+
+//! Retrieve IFP data from the IFP fission banks.
+//!
+//! \param[in] i_bank Index in the fission banks
+//! \param[in,out] delayed_groups Delayed group numbers
+//! \param[in,out] lifetimes Lifetimes lists
+void copy_ifp_data_from_fission_banks(
+  int i_bank, vector<int>& delayed_groups, vector<double>& lifetimes);
+
+#ifdef OPENMC_MPI
+
+//! Deserialization information for transfer of IFP data using MPI
+struct DeserializationInfo {
+  int64_t index_local; //!< local index
+  int64_t n;           //!< number of sites sent
+};
+
+//! Broadcast the number of generation determined by the size of the first
+//! element on the first processor.
+//!
+//! \param[in] n_generation Number of generations
+//! \param[in] delayed_groups List of delayed group numbers lists
+//! \param[in] lifetimes List of lifetimes lists
+void broadcast_ifp_n_generation(int& n_generation,
+  const vector<vector<int>>& delayed_groups,
+  const vector<vector<double>>& lifetimes);
+
+//! Send IFP data using MPI.
+//!
+//! \param[in] idx Index of the first site
+//! \param[in] n Number of sites to send
+//! \param[in] n_generation Number of generations
+//! \param[in] neighbor Index of the neighboring processor
+//! \param[in] requests MPI requests
+//! \param[in] delayed_groups List of delayed group numbers lists
+//! \param[out] send_delayed_groups Delayed group numbers buffer
+//! \param[in] lifetimes List of lifetimes lists
+//! \param[out] send_lifetimes Lifetimes buffer
+void send_ifp_info(int64_t idx, int64_t n, int n_generation, int neighbor,
+  vector<MPI_Request>& requests, const vector<vector<int>>& delayed_groups,
+  vector<int>& send_delayed_groups, const vector<vector<double>>& lifetimes,
+  vector<double>& send_lifetimes);
+
+//! Receive IFP data using MPI.
+//!
+//! \param[in] idx Index of the first site
+//! \param[in] n Number of sites to receive
+//! \param[in] n_generation Number of generations
+//! \param[in] neighbor Index of the neighboring processor
+//! \param[in] requests MPI requests
+//! \param[in] delayed_groups List of delayed group numbers
+//! \param[in] lifetimes List of lifetimes
+//! \param[out] deserialization Information to deserialize the received data
+void receive_ifp_data(int64_t idx, int64_t n, int n_generation, int neighbor,
+  vector<MPI_Request>& requests, vector<int>& delayed_groups,
+  vector<double>& lifetimes, vector<DeserializationInfo>& deserialization);
+
+//! Copy partial IFP data from local lists to source banks.
+//!
+//! \param[in] idx Index of the first site
+//! \param[in] n Number of sites to copy
+//! \param[in] i_bank Index in the IFP source banks
+//! \param[in] delayed_groups List of delayed group numbers lists
+//! \param[in] lifetimes List of lifetimes lists
+void copy_partial_ifp_data_to_source_banks(int64_t idx, int n, int64_t i_bank,
+  const vector<vector<int>>& delayed_groups,
+  const vector<vector<double>>& lifetimes);
+
+//! Deserialize IFP information received using MPI and store it in
+//! the IFP source banks.
+//!
+//! \param[in] n_generation Number of generations
+//! \param[out] deserialization Information to deserialize the received data
+//! \param[in] delayed_groups List of delayed group numbers
+//! \param[in] lifetimes List of lifetimes
+void deserialize_ifp_info(int n_generation,
+  const vector<DeserializationInfo>& deserialization,
+  const vector<int>& delayed_groups, const vector<double>& lifetimes);
+
+#endif
+
+//! Copy IFP temporary vectors to source banks.
+//!
+//! \param[in] delayed_groups List of delayed group numbers lists
+//! \param[in] lifetimes List of lifetimes lists
+void copy_complete_ifp_data_to_source_banks(
+  const vector<vector<int>>& delayed_groups,
+  const vector<vector<double>>& lifetimes);
+
+//! Allocate temporary vectors for IFP data.
+//!
+//! \param[in,out] delayed_groups List of delayed group numbers lists
+//! \param[in,out] lifetimes List of delayed group numbers lists
+void allocate_temporary_vector_ifp(
+  vector<vector<int>>& delayed_groups, vector<vector<double>>& lifetimes);
+
+//! Copy local IFP data to IFP fission banks.
+//!
+//! \param[in] delayed_groups_ptr Pointer to delayed group numbers
+//! \param[in] lifetimes_ptr Pointer to lifetimes
+void copy_ifp_data_to_fission_banks(
+  const vector<int>* delayed_groups_ptr, const vector<double>* lifetimes_ptr);
+
+} // namespace openmc
+
+#endif // OPENMC_IFP_H

include/openmc/particle_data.h

@@ -454,6 +454,9 @@ class ParticleData : public GeometryState {
 
   int cell_born_ {-1};
 
+  // Iterated Fission Probability
+  double lifetime_ {0.0}; //!< neutron lifetime [s]
+
   int n_collision_ {0};
 
   bool write_track_ {false};
@@ -560,6 +563,10 @@ class ParticleData : public GeometryState {
   double& time_last() { return time_last_; }
   const double& time_last() const { return time_last_; }
 
+  // Particle lifetime
+  double& lifetime() { return lifetime_; }
+  const double& lifetime() const { return lifetime_; }
+
   // What event took place, described in greater detail below
   TallyEvent& event() { return event_; }
   const TallyEvent& event() const { return event_; }