diff --git a/CMakeLists.txt b/CMakeLists.txt index 461abe508a..8cb03a1a1d 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/docs/source/usersguide/index.rst b/docs/source/usersguide/index.rst index aef9b1a1c1..5f8e0197e7 100644 --- a/docs/source/usersguide/index.rst +++ b/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 diff --git a/docs/source/usersguide/kinetics.rst b/docs/source/usersguide/kinetics.rst new file mode 100644 index 0000000000..bdf26d341b --- /dev/null +++ b/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 `_. 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). diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 02328af58c..e3b4e508bc 100644 --- a/docs/source/usersguide/tallies.rst +++ b/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: diff --git a/include/openmc/bank.h b/include/openmc/bank.h index 95386514d7..fd8fbd73ee 100644 --- a/include/openmc/bank.h +++ b/include/openmc/bank.h @@ -22,6 +22,14 @@ extern SharedArray surf_source_bank; extern SharedArray fission_bank; +extern vector> ifp_source_delayed_group_bank; + +extern vector> ifp_source_lifetime_bank; + +extern vector> ifp_fission_delayed_group_bank; + +extern vector> ifp_fission_lifetime_bank; + extern vector progeny_per_particle; } // namespace simulation diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 2a7ce19904..252194528c 100644 --- a/include/openmc/constants.h +++ b/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, diff --git a/include/openmc/ifp.h b/include/openmc/ifp.h new file mode 100644 index 0000000000..633a262d5f --- /dev/null +++ b/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 +void resize_ifp_data(vector& delayed_groups, vector& 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 +vector _ifp(const T& value, const vector& data) +{ + vector 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& delayed_groups, vector& 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>& delayed_groups, + const vector>& 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& requests, const vector>& delayed_groups, + vector& send_delayed_groups, const vector>& lifetimes, + vector& 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& requests, vector& delayed_groups, + vector& lifetimes, vector& 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>& delayed_groups, + const vector>& 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& deserialization, + const vector& delayed_groups, const vector& 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>& delayed_groups, + const vector>& 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>& delayed_groups, vector>& 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* delayed_groups_ptr, const vector* lifetimes_ptr); + +} // namespace openmc + +#endif // OPENMC_IFP_H diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 15ae57893a..1570b780bb 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -454,6 +454,9 @@ private: 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 @@ public: 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_; } diff --git a/include/openmc/settings.h b/include/openmc/settings.h index f3aff08b0c..9017b2d080 100644 --- a/include/openmc/settings.h +++ b/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_on; //!< 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? @@ -112,6 +121,10 @@ extern array energy_cutoff; //!< Energy cutoff in [eV] for each particle type extern array 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 diff --git a/openmc/lib/tally.py b/openmc/lib/tally.py index d0b34aedc2..c17b16597f 100644 --- a/openmc/lib/tally.py +++ b/openmc/lib/tally.py @@ -104,7 +104,9 @@ _SCORES = { -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' diff --git a/openmc/settings.py b/openmc/settings.py index a2bf51a51e..6fb36f3fc1 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -86,6 +86,9 @@ class Settings: .. versionadded:: 0.12 generations_per_batch : int Number of generations per batch + ifp_n_generation : int + Number of generations to consider for the Iterated Fission Probability + method. max_lost_particles : int Maximum number of lost particles @@ -375,6 +378,9 @@ class Settings: self._output = None + # Iterated Fission Probability + self._ifp_n_generation = None + # Output options self._statepoint = {} self._sourcepoint = {} @@ -826,6 +832,17 @@ class Settings: cv.check_less_than('verbosity', verbosity, 10, True) self._verbosity = verbosity + @property + def ifp_n_generation(self) -> int: + return self._ifp_n_generation + + @ifp_n_generation.setter + def ifp_n_generation(self, ifp_n_generation: int): + if ifp_n_generation is not None: + cv.check_type("number of generations", ifp_n_generation, Integral) + cv.check_greater_than("number of generations", ifp_n_generation, 0) + self._ifp_n_generation = ifp_n_generation + @property def tabular_legendre(self) -> dict: return self._tabular_legendre @@ -1455,6 +1472,11 @@ class Settings: element = ET.SubElement(root, "no_reduce") element.text = str(self._no_reduce).lower() + def _create_ifp_n_generation_subelement(self, root): + if self._ifp_n_generation is not None: + element = ET.SubElement(root, "ifp_n_generation") + element.text = str(self._ifp_n_generation) + def _create_tabular_legendre_subelements(self, root): if self.tabular_legendre: element = ET.SubElement(root, "tabular_legendre") @@ -1888,6 +1910,11 @@ class Settings: if text is not None: self.verbosity = int(text) + def _ifp_n_generation_from_xml_element(self, root): + text = get_text(root, 'ifp_n_generation') + if text is not None: + self.ifp_n_generation = int(text) + def _tabular_legendre_from_xml_element(self, root): elem = root.find('tabular_legendre') if elem is not None: @@ -2116,6 +2143,7 @@ class Settings: self._create_trigger_subelement(element) self._create_no_reduce_subelement(element) self._create_verbosity_subelement(element) + self._create_ifp_n_generation_subelement(element) self._create_tabular_legendre_subelements(element) self._create_temperature_subelements(element) self._create_trace_subelement(element) @@ -2225,6 +2253,7 @@ class Settings: settings._trigger_from_xml_element(elem) settings._no_reduce_from_xml_element(elem) settings._verbosity_from_xml_element(elem) + settings._ifp_n_generation_from_xml_element(elem) settings._tabular_legendre_from_xml_element(elem) settings._temperature_from_xml_element(elem) settings._trace_from_xml_element(elem) diff --git a/src/bank.cpp b/src/bank.cpp index 8d00d54409..9955939f6e 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -1,6 +1,7 @@ #include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/error.h" +#include "openmc/ifp.h" #include "openmc/message_passing.h" #include "openmc/simulation.h" #include "openmc/vector.h" @@ -26,6 +27,14 @@ SharedArray surf_source_bank; // function. SharedArray fission_bank; +vector> ifp_source_delayed_group_bank; + +vector> ifp_source_lifetime_bank; + +vector> ifp_fission_delayed_group_bank; + +vector> ifp_fission_lifetime_bank; + // Each entry in this vector corresponds to the number of progeny produced // this generation for the particle located at that index. This vector is // used to efficiently sort the fission bank after each iteration. @@ -43,6 +52,10 @@ void free_memory_bank() simulation::surf_source_bank.clear(); simulation::fission_bank.clear(); simulation::progeny_per_particle.clear(); + simulation::ifp_source_delayed_group_bank.clear(); + simulation::ifp_source_lifetime_bank.clear(); + simulation::ifp_fission_delayed_group_bank.clear(); + simulation::ifp_fission_lifetime_bank.clear(); } void init_fission_bank(int64_t max) @@ -82,6 +95,8 @@ void sort_fission_bank() // over provisioned, so we can use that as scratch space. SourceSite* sorted_bank; vector sorted_bank_holder; + vector> sorted_ifp_delayed_group_bank; + vector> sorted_ifp_lifetime_bank; // If there is not enough space, allocate a temporary vector and point to it if (simulation::fission_bank.size() > @@ -92,6 +107,11 @@ void sort_fission_bank() sorted_bank = &simulation::fission_bank[simulation::fission_bank.size()]; } + if (settings::ifp_on) { + allocate_temporary_vector_ifp( + sorted_ifp_delayed_group_bank, sorted_ifp_lifetime_bank); + } + // Use parent and progeny indices to sort fission bank for (int64_t i = 0; i < simulation::fission_bank.size(); i++) { const auto& site = simulation::fission_bank[i]; @@ -102,11 +122,19 @@ void sort_fission_bank() "shared fission bank size."); } sorted_bank[idx] = site; + if (settings::ifp_on) { + copy_ifp_data_from_fission_banks( + i, sorted_ifp_delayed_group_bank[idx], sorted_ifp_lifetime_bank[idx]); + } } // Copy sorted bank into the fission bank std::copy(sorted_bank, sorted_bank + simulation::fission_bank.size(), simulation::fission_bank.data()); + if (settings::ifp_on) { + copy_ifp_data_to_fission_banks( + sorted_ifp_delayed_group_bank.data(), sorted_ifp_lifetime_bank.data()); + } } //============================================================================== diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 8669d76f94..2685bbe98a 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -11,6 +11,7 @@ #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" +#include "openmc/ifp.h" #include "openmc/math_functions.h" #include "openmc/mesh.h" #include "openmc/message_passing.h" @@ -153,8 +154,17 @@ void synchronize_bank() // Allocate temporary source bank -- we don't really know how many fission // sites were created, so overallocate by a factor of 3 int64_t index_temp = 0; + vector temp_sites(3 * simulation::work_per_rank); + // Temporary banks for IFP + vector> temp_delayed_groups; + vector> temp_lifetimes; + if (settings::ifp_on) { + resize_ifp_data( + temp_delayed_groups, temp_lifetimes, 3 * simulation::work_per_rank); + } + for (int64_t i = 0; i < simulation::fission_bank.size(); i++) { const auto& site = simulation::fission_bank[i]; @@ -165,6 +175,10 @@ void synchronize_bank() if (total < settings::n_particles) { for (int64_t j = 1; j <= settings::n_particles / total; ++j) { temp_sites[index_temp] = site; + if (settings::ifp_on) { + copy_ifp_data_from_fission_banks( + i, temp_delayed_groups[index_temp], temp_lifetimes[index_temp]); + } ++index_temp; } } @@ -172,6 +186,10 @@ void synchronize_bank() // Randomly sample sites needed if (prn(&seed) < p_sample) { temp_sites[index_temp] = site; + if (settings::ifp_on) { + copy_ifp_data_from_fission_banks( + i, temp_delayed_groups[index_temp], temp_lifetimes[index_temp]); + } ++index_temp; } } @@ -188,6 +206,8 @@ void synchronize_bank() MPI_Exscan(&index_temp, &start, 1, MPI_INT64_T, MPI_SUM, mpi::intracomm); finish = start + index_temp; + // TODO: protect for MPI_Exscan at rank 0 + // Allocate space for bank_position if this hasn't been done yet int64_t bank_position[mpi::n_procs]; MPI_Allgather( @@ -211,9 +231,15 @@ void synchronize_bank() // If we have too few sites, repeat sites from the very end of the // fission bank sites_needed = settings::n_particles - finish; + // TODO: sites_needed > simulation::fission_bank.size() or other test to + // make sure we don't need info from other proc for (int i = 0; i < sites_needed; ++i) { int i_bank = simulation::fission_bank.size() - sites_needed + i; temp_sites[index_temp] = simulation::fission_bank[i_bank]; + if (settings::ifp_on) { + copy_ifp_data_from_fission_banks(i_bank, + temp_delayed_groups[index_temp], temp_lifetimes[index_temp]); + } ++index_temp; } } @@ -229,15 +255,32 @@ void synchronize_bank() // ========================================================================== // SEND BANK SITES TO NEIGHBORS + // IFP number of generation + int ifp_n_generation; + if (settings::ifp_on) { + broadcast_ifp_n_generation( + ifp_n_generation, temp_delayed_groups, temp_lifetimes); + } + int64_t index_local = 0; vector requests; + // IFP send buffers + vector send_delayed_groups; + vector send_lifetimes; + if (start < settings::n_particles) { // Determine the index of the processor which has the first part of the // source_bank for the local processor int neighbor = upper_bound_index( simulation::work_index.begin(), simulation::work_index.end(), start); + // Resize IFP send buffers + if (settings::ifp_on && mpi::n_procs > 1) { + resize_ifp_data(send_delayed_groups, send_lifetimes, + ifp_n_generation * 3 * simulation::work_per_rank); + } + while (start < finish) { // Determine the number of sites to send int64_t n = @@ -250,6 +293,13 @@ void synchronize_bank() MPI_Isend(&temp_sites[index_local], static_cast(n), mpi::source_site, neighbor, mpi::rank, mpi::intracomm, &requests.back()); + + if (settings::ifp_on) { + // Send IFP data + send_ifp_info(index_local, n, ifp_n_generation, neighbor, requests, + temp_delayed_groups, send_delayed_groups, temp_lifetimes, + send_lifetimes); + } } // Increment all indices @@ -271,6 +321,11 @@ void synchronize_bank() start = simulation::work_index[mpi::rank]; index_local = 0; + // IFP receive buffers + vector recv_delayed_groups; + vector recv_lifetimes; + vector deserialization_info; + // Determine what process has the source sites that will need to be stored at // the beginning of this processor's source bank. @@ -282,6 +337,12 @@ void synchronize_bank() upper_bound_index(bank_position, bank_position + mpi::n_procs, start); } + // Resize IFP receive buffers + if (settings::ifp_on && mpi::n_procs > 1) { + resize_ifp_data(recv_delayed_groups, recv_lifetimes, + ifp_n_generation * simulation::work_per_rank); + } + while (start < simulation::work_index[mpi::rank + 1]) { // Determine how many sites need to be received int64_t n; @@ -301,13 +362,24 @@ void synchronize_bank() MPI_Irecv(&simulation::source_bank[index_local], static_cast(n), mpi::source_site, neighbor, neighbor, mpi::intracomm, &requests.back()); + if (settings::ifp_on) { + // Receive IFP data + receive_ifp_data(index_local, n, ifp_n_generation, neighbor, requests, + recv_delayed_groups, recv_lifetimes, deserialization_info); + } + } else { - // If the source sites are on this procesor, we can simply copy them + // If the source sites are on this processor, we can simply copy them // from the temp_sites bank index_temp = start - bank_position[mpi::rank]; std::copy(&temp_sites[index_temp], &temp_sites[index_temp + n], &simulation::source_bank[index_local]); + + if (settings::ifp_on) { + copy_partial_ifp_data_to_source_banks( + index_temp, n, index_local, temp_delayed_groups, temp_lifetimes); + } } // Increment all indices @@ -323,9 +395,17 @@ void synchronize_bank() int n_request = requests.size(); MPI_Waitall(n_request, requests.data(), MPI_STATUSES_IGNORE); + if (settings::ifp_on) { + deserialize_ifp_info(ifp_n_generation, deserialization_info, + recv_delayed_groups, recv_lifetimes); + } + #else std::copy(temp_sites.data(), temp_sites.data() + settings::n_particles, simulation::source_bank.begin()); + if (settings::ifp_on) { + copy_complete_ifp_data_to_source_banks(temp_delayed_groups, temp_lifetimes); + } #endif simulation::time_bank_sendrecv.stop(); diff --git a/src/finalize.cpp b/src/finalize.cpp index ed3d0e7f46..54aa1d1d9e 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -171,8 +171,9 @@ int openmc_finalize() // Free all MPI types #ifdef OPENMC_MPI - if (mpi::source_site != MPI_DATATYPE_NULL) + if (mpi::source_site != MPI_DATATYPE_NULL) { MPI_Type_free(&mpi::source_site); + } #endif openmc_reset_random_ray(); diff --git a/src/ifp.cpp b/src/ifp.cpp new file mode 100644 index 0000000000..1f81f26f6e --- /dev/null +++ b/src/ifp.cpp @@ -0,0 +1,216 @@ +#include "openmc/ifp.h" + +#include "openmc/bank.h" +#include "openmc/message_passing.h" +#include "openmc/particle.h" +#include "openmc/particle_data.h" +#include "openmc/settings.h" +#include "openmc/simulation.h" +#include "openmc/vector.h" + +namespace openmc { + +bool is_beta_effective_or_both() +{ + if (settings::ifp_parameter == IFPParameter::BetaEffective || + settings::ifp_parameter == IFPParameter::Both) { + return true; + } + return false; +} + +bool is_generation_time_or_both() +{ + if (settings::ifp_parameter == IFPParameter::GenerationTime || + settings::ifp_parameter == IFPParameter::Both) { + return true; + } + return false; +} + +void ifp(const Particle& p, const SourceSite& site, int64_t idx) +{ + if (is_beta_effective_or_both()) { + const auto& delayed_groups = + simulation::ifp_source_delayed_group_bank[p.current_work() - 1]; + simulation::ifp_fission_delayed_group_bank[idx] = + _ifp(site.delayed_group, delayed_groups); + } + if (is_generation_time_or_both()) { + const auto& lifetimes = + simulation::ifp_source_lifetime_bank[p.current_work() - 1]; + simulation::ifp_fission_lifetime_bank[idx] = _ifp(p.lifetime(), lifetimes); + } +} + +void resize_simulation_ifp_banks() +{ + resize_ifp_data(simulation::ifp_source_delayed_group_bank, + simulation::ifp_source_lifetime_bank, simulation::work_per_rank); + resize_ifp_data(simulation::ifp_fission_delayed_group_bank, + simulation::ifp_fission_lifetime_bank, 3 * simulation::work_per_rank); +} + +void copy_ifp_data_from_fission_banks( + int i_bank, vector& delayed_groups, vector& lifetimes) +{ + if (is_beta_effective_or_both()) { + delayed_groups = simulation::ifp_fission_delayed_group_bank[i_bank]; + } + if (is_generation_time_or_both()) { + lifetimes = simulation::ifp_fission_lifetime_bank[i_bank]; + } +} + +#ifdef OPENMC_MPI + +void broadcast_ifp_n_generation(int& n_generation, + const vector>& delayed_groups, + const vector>& lifetimes) +{ + if (mpi::rank == 0) { + if (is_beta_effective_or_both()) { + n_generation = static_cast(delayed_groups[0].size()); + } else { + n_generation = static_cast(lifetimes[0].size()); + } + } + MPI_Bcast(&n_generation, 1, MPI_INT, 0, mpi::intracomm); +} + +void send_ifp_info(int64_t idx, int64_t n, int n_generation, int neighbor, + vector& requests, const vector>& delayed_groups, + vector& send_delayed_groups, const vector>& lifetimes, + vector& send_lifetimes) +{ + // Copy data in send buffers + for (int i = idx; i < idx + n; i++) { + if (is_beta_effective_or_both()) { + std::copy(delayed_groups[i].begin(), delayed_groups[i].end(), + send_delayed_groups.begin() + i * n_generation); + } + if (is_generation_time_or_both()) { + std::copy(lifetimes[i].begin(), lifetimes[i].end(), + send_lifetimes.begin() + i * n_generation); + } + } + // Send delayed groups + if (is_beta_effective_or_both()) { + requests.emplace_back(); + MPI_Isend(&send_delayed_groups[n_generation * idx], + n_generation * static_cast(n), MPI_INT, neighbor, mpi::rank, + mpi::intracomm, &requests.back()); + } + // Send lifetimes + if (is_generation_time_or_both()) { + requests.emplace_back(); + MPI_Isend(&send_lifetimes[n_generation * idx], + n_generation * static_cast(n), MPI_DOUBLE, neighbor, mpi::rank, + mpi::intracomm, &requests.back()); + } +} + +void receive_ifp_data(int64_t idx, int64_t n, int n_generation, int neighbor, + vector& requests, vector& delayed_groups, + vector& lifetimes, vector& deserialization) +{ + // Receive delayed groups + if (is_beta_effective_or_both()) { + requests.emplace_back(); + MPI_Irecv(&delayed_groups[n_generation * idx], + n_generation * static_cast(n), MPI_INT, neighbor, neighbor, + mpi::intracomm, &requests.back()); + } + // Receive lifetimes + if (is_generation_time_or_both()) { + requests.emplace_back(); + MPI_Irecv(&lifetimes[n_generation * idx], + n_generation * static_cast(n), MPI_DOUBLE, neighbor, neighbor, + mpi::intracomm, &requests.back()); + } + // Deserialization info to reconstruct data later + DeserializationInfo info = {idx, n}; + deserialization.push_back(info); +} + +void copy_partial_ifp_data_to_source_banks(int64_t idx, int n, int64_t i_bank, + const vector>& delayed_groups, + const vector>& lifetimes) +{ + if (is_beta_effective_or_both()) { + std::copy(&delayed_groups[idx], &delayed_groups[idx + n], + &simulation::ifp_source_delayed_group_bank[i_bank]); + } + if (is_generation_time_or_both()) { + std::copy(&lifetimes[idx], &lifetimes[idx + n], + &simulation::ifp_source_lifetime_bank[i_bank]); + } +} + +void deserialize_ifp_info(int n_generation, + const vector& deserialization, + const vector& delayed_groups, const vector& lifetimes) +{ + for (auto info : deserialization) { + int64_t index_local = info.index_local; + int64_t n = info.n; + + for (int i = index_local; i < index_local + n; i++) { + if (is_beta_effective_or_both()) { + vector delayed_groups_received( + delayed_groups.begin() + n_generation * i, + delayed_groups.begin() + n_generation * (i + 1)); + simulation::ifp_source_delayed_group_bank[i] = delayed_groups_received; + } + if (is_generation_time_or_both()) { + vector lifetimes_received(lifetimes.begin() + n_generation * i, + lifetimes.begin() + n_generation * (i + 1)); + simulation::ifp_source_lifetime_bank[i] = lifetimes_received; + } + } + } +} + +#endif + +void copy_complete_ifp_data_to_source_banks( + const vector>& delayed_groups, + const vector>& lifetimes) +{ + if (is_beta_effective_or_both()) { + std::copy(delayed_groups.data(), + delayed_groups.data() + settings::n_particles, + simulation::ifp_source_delayed_group_bank.begin()); + } + if (is_generation_time_or_both()) { + std::copy(lifetimes.data(), lifetimes.data() + settings::n_particles, + simulation::ifp_source_lifetime_bank.begin()); + } +} + +void allocate_temporary_vector_ifp( + vector>& delayed_groups, vector>& lifetimes) +{ + if (is_beta_effective_or_both()) { + delayed_groups.resize(simulation::fission_bank.size()); + } + if (is_generation_time_or_both()) { + lifetimes.resize(simulation::fission_bank.size()); + } +} + +void copy_ifp_data_to_fission_banks(const vector* const delayed_groups_ptr, + const vector* lifetimes_ptr) +{ + if (is_beta_effective_or_both()) { + std::copy(delayed_groups_ptr, + delayed_groups_ptr + simulation::fission_bank.size(), + simulation::ifp_fission_delayed_group_bank.data()); + } + if (is_generation_time_or_both()) { + std::copy(lifetimes_ptr, lifetimes_ptr + simulation::fission_bank.size(), + simulation::ifp_fission_lifetime_bank.data()); + } +} + +} // namespace openmc diff --git a/src/output.cpp b/src/output.cpp index e20868efbb..baaf682a15 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -593,6 +593,9 @@ const std::unordered_map score_names = { {SCORE_FISS_Q_RECOV, "Recoverable fission power"}, {SCORE_CURRENT, "Current"}, {SCORE_PULSE_HEIGHT, "pulse-height"}, + {SCORE_IFP_TIME_NUM, "IFP lifetime numerator"}, + {SCORE_IFP_BETA_NUM, "IFP delayed fraction numerator"}, + {SCORE_IFP_DENOM, "IFP common denominator"}, }; //! Create an ASCII output file showing all tally results. diff --git a/src/particle.cpp b/src/particle.cpp index c51011d6b7..324276d15f 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -116,6 +116,7 @@ void Particle::from_source(const SourceSite* src) n_collision() = 0; fission() = false; zero_flux_derivs(); + lifetime() = 0.0; // Copy attributes from source bank site type() = src->particle; @@ -235,7 +236,9 @@ void Particle::event_advance() for (int j = 0; j < n_coord(); ++j) { coord(j).r += distance * coord(j).u; } - this->time() += distance / this->speed(); + double dt = distance / this->speed(); + this->time() += dt; + this->lifetime() += dt; // Kill particle if its time exceeds the cutoff bool hit_time_boundary = false; @@ -243,6 +246,7 @@ void Particle::event_advance() if (time() > time_cutoff) { double dt = time() - time_cutoff; time() = time_cutoff; + lifetime() = time_cutoff; double push_back_distance = speed() * dt; this->move_distance(-push_back_distance); diff --git a/src/physics.cpp b/src/physics.cpp index 72c04a5caf..a8e5b9e813 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -8,6 +8,7 @@ #include "openmc/eigenvalue.h" #include "openmc/endf.h" #include "openmc/error.h" +#include "openmc/ifp.h" #include "openmc/material.h" #include "openmc/math_functions.h" #include "openmc/message_passing.h" @@ -233,6 +234,10 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // Break out of loop as no more sites can be added to fission bank break; } + // Iterated Fission Probability (IFP) method + if (settings::ifp_on) { + ifp(p, site, idx); + } } else { p.secondary_bank().push_back(site); } diff --git a/src/reaction.cpp b/src/reaction.cpp index 9714734383..d96790c6d4 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -202,6 +202,9 @@ std::unordered_map REACTION_NAME_MAP { {SCORE_FISS_Q_PROMPT, "fission-q-prompt"}, {SCORE_FISS_Q_RECOV, "fission-q-recoverable"}, {SCORE_PULSE_HEIGHT, "pulse-height"}, + {SCORE_IFP_TIME_NUM, "ifp-time-numerator"}, + {SCORE_IFP_BETA_NUM, "ifp-beta-numerator"}, + {SCORE_IFP_DENOM, "ifp-denominator"}, // Normal ENDF-based reactions {TOTAL_XS, "(n,total)"}, {ELASTIC, "(n,elastic)"}, diff --git a/src/settings.cpp b/src/settings.cpp index 53fcdec38b..c135fa3391 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -52,6 +52,7 @@ bool create_fission_neutrons {true}; bool delayed_photon_scaling {true}; bool entropy_on {false}; bool event_based {false}; +bool ifp_on {false}; bool legendre_to_tabular {true}; bool material_cell_offsets {true}; bool output_summary {true}; @@ -106,6 +107,8 @@ int max_particle_events {1000000}; ElectronTreatment electron_treatment {ElectronTreatment::TTB}; array energy_cutoff {0.0, 1000.0, 0.0, 0.0}; array time_cutoff {INFTY, INFTY, INFTY, INFTY}; +int ifp_n_generation {-1}; +IFPParameter ifp_parameter {IFPParameter::None}; int legendre_to_tabular_points {C_NONE}; int max_order {0}; int n_log_bins {8000}; @@ -1059,6 +1062,20 @@ void read_settings_xml(pugi::xml_node root) temperature_range[1] = range.at(1); } + // Check for user value for the number of generation of the Iterated Fission + // Probability (IFP) method + if (check_for_node(root, "ifp_n_generation")) { + ifp_n_generation = std::stoi(get_node_value(root, "ifp_n_generation")); + if (ifp_n_generation <= 0) { + fatal_error("'ifp_n_generation' must be greater than 0."); + } + // Avoid tallying 0 if IFP logs are not complete when active cycles start + if (ifp_n_generation > n_inactive) { + fatal_error("'ifp_n_generation' must be lower than or equal to the " + "number of inactive cycles."); + } + } + // Check for tabular_legendre options if (check_for_node(root, "tabular_legendre")) { // Get pointer to tabular_legendre node diff --git a/src/simulation.cpp b/src/simulation.cpp index 74a7dbe639..e4521c06dc 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -7,6 +7,7 @@ #include "openmc/error.h" #include "openmc/event.h" #include "openmc/geometry_aux.h" +#include "openmc/ifp.h" #include "openmc/material.h" #include "openmc/mcpl_interface.h" #include "openmc/message_passing.h" @@ -335,6 +336,11 @@ void allocate_banks() // Allocate fission bank init_fission_bank(3 * simulation::work_per_rank); + + // Allocate IFP bank + if (settings::ifp_on) { + resize_simulation_ifp_banks(); + } } if (settings::surf_source_write) { diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 96d684f71a..35805b20d7 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -180,6 +180,64 @@ Tally::Tally(pugi::xml_node node) fatal_error(fmt::format("No scores specified on tally {}.", id_)); } + // Set IFP if needed + if (!settings::ifp_on) { + // Determine if this tally has an IFP score + bool has_ifp_score = false; + for (int score : scores_) { + if (score == SCORE_IFP_TIME_NUM || score == SCORE_IFP_BETA_NUM || + score == SCORE_IFP_DENOM) { + has_ifp_score = true; + break; + } + } + + // Check for errors + if (has_ifp_score) { + if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::ifp_n_generation < 0) { + settings::ifp_n_generation = DEFAULT_IFP_N_GENERATION; + warning(fmt::format( + "{} generations will be used for IFP (default value). It can be " + "changed using the 'ifp_n_generation' settings.", + settings::ifp_n_generation)); + } + if (settings::ifp_n_generation > settings::n_inactive) { + fatal_error("'ifp_n_generation' must be lower than or equal to the " + "number of inactive cycles."); + } + settings::ifp_on = true; + } else { + fatal_error( + "Iterated Fission Probability can only be used in an eigenvalue " + "calculation."); + } + } + } + + // Set IFP parameters if needed + if (settings::ifp_on) { + for (int score : scores_) { + switch (score) { + case SCORE_IFP_TIME_NUM: + if (settings::ifp_parameter == IFPParameter::None) { + settings::ifp_parameter = IFPParameter::GenerationTime; + } else if (settings::ifp_parameter == IFPParameter::BetaEffective) { + settings::ifp_parameter = IFPParameter::Both; + } + break; + case SCORE_IFP_BETA_NUM: + case SCORE_IFP_DENOM: + if (settings::ifp_parameter == IFPParameter::None) { + settings::ifp_parameter = IFPParameter::BetaEffective; + } else if (settings::ifp_parameter == IFPParameter::GenerationTime) { + settings::ifp_parameter = IFPParameter::Both; + } + break; + } + } + } + // Check if tally is compatible with particle type if (!settings::photon_transport) { for (int score : scores_) { @@ -585,7 +643,12 @@ void Tally::set_scores(const vector& scores) } } } + break; + case SCORE_IFP_TIME_NUM: + case SCORE_IFP_BETA_NUM: + case SCORE_IFP_DENOM: + estimator_ = TallyEstimator::COLLISION; break; } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 02cb485671..5c6386fb42 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -4,6 +4,7 @@ #include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/error.h" +#include "openmc/ifp.h" #include "openmc/material.h" #include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" @@ -890,6 +891,56 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, score_fission_q(p, score_bin, tally, flux, i_nuclide, atom_density); break; + case SCORE_IFP_TIME_NUM: + if (settings::ifp_on) { + if ((p.type() == Type::neutron) && (p.fission())) { + if (is_generation_time_or_both()) { + const auto& lifetimes = + simulation::ifp_source_lifetime_bank[p.current_work() - 1]; + if (lifetimes.size() == settings::ifp_n_generation) { + score = lifetimes[0] * p.wgt_last(); + } + } + } + } + break; + + case SCORE_IFP_BETA_NUM: + if (settings::ifp_on) { + if ((p.type() == Type::neutron) && (p.fission())) { + if (is_beta_effective_or_both()) { + const auto& delayed_groups = + simulation::ifp_source_delayed_group_bank[p.current_work() - 1]; + if (delayed_groups.size() == settings::ifp_n_generation) { + if (delayed_groups[0] > 0) { + score = p.wgt_last(); + } + } + } + } + } + break; + + case SCORE_IFP_DENOM: + if (settings::ifp_on) { + if ((p.type() == Type::neutron) && (p.fission())) { + int ifp_data_size; + if (is_beta_effective_or_both()) { + ifp_data_size = static_cast( + simulation::ifp_source_delayed_group_bank[p.current_work() - 1] + .size()); + } else { + ifp_data_size = static_cast( + simulation::ifp_source_lifetime_bank[p.current_work() - 1] + .size()); + } + if (ifp_data_size == settings::ifp_n_generation) { + score = p.wgt_last(); + } + } + } + break; + case N_2N: case N_3N: case N_4N: diff --git a/tests/regression_tests/ifp/__init__.py b/tests/regression_tests/ifp/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/ifp/inputs_true.dat b/tests/regression_tests/ifp/inputs_true.dat new file mode 100644 index 0000000000..a3a3f1d77e --- /dev/null +++ b/tests/regression_tests/ifp/inputs_true.dat @@ -0,0 +1,33 @@ + + + + + + + + + + + + + + eigenvalue + 1000 + 20 + 5 + + + -10.0 -10.0 -10.0 10.0 10.0 10.0 + + + true + + + 5 + + + + ifp-time-numerator ifp-beta-numerator ifp-denominator + + + diff --git a/tests/regression_tests/ifp/results_true.dat b/tests/regression_tests/ifp/results_true.dat new file mode 100644 index 0000000000..a74e2bd78b --- /dev/null +++ b/tests/regression_tests/ifp/results_true.dat @@ -0,0 +1,9 @@ +k-combined: +1.007452E+00 5.705278E-03 +tally 1: +8.996235E-08 +5.461421E-16 +4.800000E-02 +4.680000E-04 +1.512000E+01 +1.526063E+01 diff --git a/tests/regression_tests/ifp/test.py b/tests/regression_tests/ifp/test.py new file mode 100644 index 0000000000..6969a54c49 --- /dev/null +++ b/tests/regression_tests/ifp/test.py @@ -0,0 +1,45 @@ +"""Test the Iterated Fission Probability (IFP) method to compute adjoint-weighted +kinetics parameters using dedicated tallies.""" + +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture() +def ifp_model(): + model = openmc.Model() + + # Material + material = openmc.Material(name="core") + material.add_nuclide("U235", 1.0) + material.set_density('g/cm3', 16.0) + + # Geometry + radius = 10.0 + sphere = openmc.Sphere(r=radius, boundary_type="vacuum") + cell = openmc.Cell(region=-sphere, fill=material) + model.geometry = openmc.Geometry([cell]) + + # Settings + model.settings.particles = 1000 + model.settings.batches = 20 + model.settings.inactive = 5 + model.settings.ifp_n_generation = 5 + + space = openmc.stats.Box(*cell.bounding_box) + model.settings.source = openmc.IndependentSource( + space=space, constraints={'fissionable': True}) + + # Tally IFP scores + tally = openmc.Tally(name="ifp-scores") + tally.scores = ["ifp-time-numerator", "ifp-beta-numerator", "ifp-denominator"] + model.tallies = [tally] + + return model + + +def test_iterated_fission_probability(ifp_model): + harness = PyAPITestHarness("statepoint.20.h5", model=ifp_model) + harness.main() diff --git a/tests/unit_tests/test_ifp.py b/tests/unit_tests/test_ifp.py new file mode 100644 index 0000000000..8d0fd98010 --- /dev/null +++ b/tests/unit_tests/test_ifp.py @@ -0,0 +1,49 @@ +"""Test the Iterated Fission Probability (IFP) method to compute adjoint-weighted +kinetics parameters using dedicated tallies.""" + +import pytest +import openmc + + +def test_xml_serialization(run_in_tmpdir): + """Check that a simple use case can be written and read in XML.""" + parameter = 5 + settings = openmc.Settings() + settings.ifp_n_generation = parameter + settings.export_to_xml() + + read_settings = openmc.Settings.from_xml() + assert read_settings.ifp_n_generation == parameter + + +@pytest.fixture(scope="module") +def geometry(): + openmc.reset_auto_ids() + material = openmc.Material() + material.add_nuclide("U235", 1.0) + sphere = openmc.Sphere(r=1.0, boundary_type="vacuum") + cell = openmc.Cell(region=-sphere, fill=material) + return openmc.Geometry([cell]) + + +@pytest.mark.parametrize( + "options, error", + [ + ({"ifp_n_generation": 0}, ValueError), + ({"ifp_n_generation": -1}, ValueError), + ({"run_mode": "fixed source"}, RuntimeError), + ({"inactive": 5, "ifp_n_generation": 6}, RuntimeError), + ({"inactive": 9}, RuntimeError) + ], +) +def test_exceptions(options, error, run_in_tmpdir, geometry): + """Test settings configuration that should return an error.""" + with pytest.raises(error): + settings = openmc.Settings(**options) + settings.particles = 100 + settings.batches = 15 + tally = openmc.Tally(name="ifp-scores") + tally.scores = ["ifp-time-numerator", "ifp-beta-numerator", "ifp-denominator"] + tallies = openmc.Tallies([tally]) + model = openmc.Model(geometry=geometry, settings=settings, tallies=tallies) + model.run()