Multi-group capability for kinetics parameter calculations with Iterated Fission Probability

CMakeLists.txt

@@ -353,6 +353,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

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

@@ -308,7 +308,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
@@ -318,6 +321,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,145 @@
+#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 {
+
+// ---------------------------------------------------------
+// Common helpers
+// ---------------------------------------------------------
+
+bool is_beta_effective_or_both();
+
+bool is_generation_time_or_both();
+
+// ---------------------------------------------------------
+// physics.cpp functions
+// ---------------------------------------------------------
+
+template<typename T>
+vector<T> _ifp(const T& value, const vector<T>& data)
+{
+  vector<T> updated_data;
+  size_t ifp_idx = data.size();
+  if (ifp_idx < settings::ifp_n_generation) {
+    updated_data.resize(ifp_idx + 1);
+    for (size_t i = 0; i < ifp_idx; i++) {
+      updated_data[i] = data[i];
+    }
+    updated_data[ifp_idx] = value;
+  } else if (ifp_idx == settings::ifp_n_generation) {
+    updated_data.resize(ifp_idx);
+    for (size_t i = 0; i < ifp_idx - 1; i++) {
+      updated_data[i] = data[i + 1];
+    }
+    updated_data[ifp_idx - 1] = value;
+  }
+  return updated_data;
+}
+
+//! Iterated Fission Probability (IFP) method
+//!
+//! Needs to be done after the delayed group is found.
+//!
+//! 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.
+void ifp(const Particle& p, const SourceSite& site, int64_t idx);
+
+// ---------------------------------------------------------
+// simulation.cpp functions
+// ---------------------------------------------------------
+
+//! Resize the IFP banks used in the simulation
+void resize_simulation_ifp_banks();
+
+// ---------------------------------------------------------
+// eigenvalue.cpp functions
+// ---------------------------------------------------------
+
+//! Resize IFP vectors
+template<typename T, typename U>
+void resize_ifp_data(
+  vector<T>& delayed_groups, vector<U>& lifetimes, int64_t n_dg, int64_t n_l)
+{
+  if (is_beta_effective_or_both()) {
+    delayed_groups.resize(n_dg);
+  }
+  if (is_generation_time_or_both()) {
+    lifetimes.resize(n_l);
+  }
+}
+
+void initialize_ifp_pointers(int64_t i_bank,
+  const vector<int>*& delayed_groups_ptr, const vector<double>*& lifetimes_ptr);
+
+void add_ifp_data(int64_t idx, vector<vector<int>>& delayed_groups,
+  const vector<int>* delayed_groups_ptr, vector<vector<double>>& lifetimes,
+  const vector<double>* lifetimes_ptr);
+
+void retrieve_ifp_data_from_fission_banks(int64_t idx, int i_bank,
+  vector<vector<int>>& delayed_groups, vector<vector<double>>& lifetimes);
+
+#ifdef OPENMC_MPI
+
+//! Deserialization info for IFP transfer via MPI
+struct DeserializationInfo {
+  int64_t index_local;
+  int64_t n;
+};
+
+//! Broadcast the number of generation from the size of the first element
+void broadcast_ifp_n_generation(int& n_generation,
+  const vector<vector<int>>& delayed_groups,
+  const vector<vector<double>>& lifetimes);
+
+//! Send IFP info via MPI
+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 info through MPI
+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);
+
+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
+void deserialize_ifp_info(int n_generation,
+  const vector<DeserializationInfo>& deserialization,
+  const vector<int>& delayed_groups, const vector<double>& lifetimes);
+
+#endif
+
+//! Copy IFP temporary vector to source banks
+void copy_complete_ifp_data_to_source_banks(
+  const vector<vector<int>>& delayed_groups,
+  const vector<vector<double>>& lifetimes);
+
+// ---------------------------------------------------------
+// bank.cpp functions
+// ---------------------------------------------------------
+
+void allocate_temporary_vector_ifp(
+  vector<vector<int>>& delayed_group_bank_holder,
+  vector<int>*& delayed_group_bank,
+  vector<vector<double>>& lifetime_bank_holder, vector<double>*& lifetime_bank);
+
+void sort_ifp_data_from_fission_banks(int64_t i_bank, int64_t idx,
+  vector<int>*& delayed_group_bank, vector<double>*& lifetime_bank);
+
+void copy_ifp_data_to_fission_banks(
+  vector<int>* const& delayed_group_bank, vector<double>* const& lifetime_bank);
+
+} // namespace openmc
+
+#endif // OPENMC_IFP_H

include/openmc/particle_data.h

@@ -424,6 +424,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};
@@ -522,6 +525,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_; }

include/openmc/settings.h

@@ -24,6 +24,14 @@ enum class SSWCellType {
   To,
 };
 
+// Type of IFP parameters
+enum class IFPParameter {
+  None,
+  Both,
+  BetaEffective,
+  GenerationTime,
+};
+
 //==============================================================================
 // Global variable declarations
 //==============================================================================
@@ -42,7 +50,8 @@ extern bool
   delayed_photon_scaling;   //!< Scale fission photon yield to include delayed
 extern "C" bool entropy_on; //!< calculate Shannon entropy?
 extern "C" bool
-  event_based; //!< use event-based mode (instead of history-based)
+  event_based;   //!< use event-based mode (instead of history-based)
+extern bool ifp; //!< Use IFP for kinetics parameters?
 extern bool legendre_to_tabular; //!< convert Legendre distributions to tabular?
 extern bool material_cell_offsets; //!< create material cells offsets?
 extern "C" bool output_summary;    //!< write summary.h5?
@@ -109,6 +118,10 @@ extern array<double, 4>
   energy_cutoff; //!< Energy cutoff in [eV] for each particle type
 extern array<double, 4>
   time_cutoff; //!< Time cutoff in [s] for each particle type
+extern int
+  ifp_n_generation; //!< Number of generation for Iterated Fission Probability
+extern IFPParameter
+  ifp_parameter; //!< Parameter to calculate for Iterated Fission Probability
 extern int
   legendre_to_tabular_points; //!< number of points to convert Legendres
 extern int max_order;         //!< Maximum Legendre order for multigroup data

openmc/lib/tally.py

@@ -104,7 +104,9 @@
     -5: 'absorption', -6: 'fission', -7: 'nu-fission', -8: 'kappa-fission',
     -9: 'current', -10: 'events', -11: 'delayed-nu-fission',
     -12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt',
-    -15: 'fission-q-recoverable', -16: 'decay-rate', -17: 'pulse-height'
+    -15: 'fission-q-recoverable', -16: 'decay-rate', -17: 'pulse-height',
+    -18: 'ifp-time-numerator', -19: 'ifp-beta-numerator',
+    -20: 'ifp-denominator',
 }
 _ESTIMATORS = {
     0: 'analog', 1: 'tracklength', 2: 'collision'

openmc/model/model.py

@@ -1,6 +1,6 @@
 from __future__ import annotations
 from collections.abc import Iterable
-from functools import lru_cache
+from functools import cache
 import os
 from pathlib import Path
 from numbers import Integral
@@ -81,6 +81,9 @@ def __init__(self, geometry=None, materials=None, settings=None,
         if plots is not None:
             self.plots = plots
 
+        if settings.iterated_fission_probability:
+            self._init_ifp()
+
     @property
     def geometry(self) -> openmc.Geometry | None:
         return self._geometry
@@ -150,7 +153,7 @@ def is_initialized(self) -> bool:
             return False
 
     @property
-    @lru_cache(maxsize=None)
+    @cache
     def _materials_by_id(self) -> dict:
         """Dictionary mapping material ID --> material"""
         if self.materials:
@@ -160,14 +163,14 @@ def _materials_by_id(self) -> dict:
         return {mat.id: mat for mat in mats}
 
     @property
-    @lru_cache(maxsize=None)
+    @cache
     def _cells_by_id(self) -> dict:
         """Dictionary mapping cell ID --> cell"""
         cells = self.geometry.get_all_cells()
         return {cell.id: cell for cell in cells.values()}
 
     @property
-    @lru_cache(maxsize=None)
+    @cache
     def _cells_by_name(self) -> dict[int, openmc.Cell]:
         # Get the names maps, but since names are not unique, store a set for
         # each name key. In this way when the user requests a change by a name,
@@ -180,7 +183,7 @@ def _cells_by_name(self) -> dict[int, openmc.Cell]:
         return result
 
     @property
-    @lru_cache(maxsize=None)
+    @cache
     def _materials_by_name(self) -> dict[int, openmc.Material]:
         if self.materials is None:
             mats = self.geometry.get_all_materials().values()
@@ -193,6 +196,31 @@ def _materials_by_name(self) -> dict[int, openmc.Material]:
             result[mat.name].add(mat)
         return result
 
+    def _init_ifp(self) -> None:
+        """Automate tally creation for calculating Iterated Fission Probability kinetics parameters"""
+        param = self.settings.iterated_fission_probability['parameter']
+        if param == 'both' or param == 'generation_time':
+            gen_time_tally = openmc.Tally()
+            gen_time_tally.scores = ['ifp-time-numerator']
+            self._tallies.append(gen_time_tally)
+
+        if param == 'both' or param == 'beta_effective':
+            beta_tally = openmc.Tally()
+            beta_tally.scores = ['ifp-beta-numerator']
+
+            has_dg_filter = False
+            for tally in self._tallies:
+                if tally.contains_filter(openmc.DelayedGroupFilter):
+                    has_dg_filter = True
+                    beta_tally.filters = [tally.find_filter(openmc.DelayedGroupFilter)]
+            if not has_dg_filter:
+                beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1,7)))]
+            self._tallies.append(beta_tally)
+
+        denom_tally = openmc.Tally()
+        denom_tally.scores = ['ifp-denominator']
+        self._tallies.append(denom_tally)
+
     @classmethod
     def from_xml(cls, geometry='geometry.xml', materials='materials.xml',
                  settings='settings.xml', tallies='tallies.xml',