diff --git a/docs/source/usersguide/kinetics.rst b/docs/source/usersguide/kinetics.rst index 9024ff8227a..85d5035a264 100644 --- a/docs/source/usersguide/kinetics.rst +++ b/docs/source/usersguide/kinetics.rst @@ -84,6 +84,11 @@ total :math:`\beta_{\text{eff}}` specified by providing a 6-group # Add DelayedGroupFilter to enable group-wise tallies beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, 7)))] +A breakdown on isotopes of the :math:`\beta_{\text{eff}}` can be estimated, +e.g., for :math:`^{U235}U` and :math:`^{U238}U`, by defining:: + + beta_tally.nuclides = ['U235', 'U238'] + Here is an example showing how to declare the three available IFP scores in a single tally:: @@ -113,7 +118,7 @@ for ``ifp-denominator``: \beta_{\text{eff}} = \frac{S_{\text{ifp-beta-numerator}}}{S_{\text{ifp-denominator}}} The kinetics parameters can be retrieved directly from a statepoint file using -the :meth:`openmc.StatePoint.ifp_results` method:: +the :meth:`openmc.StatePoint.get_kinetics_parameters` method:: with openmc.StatePoint(output_path) as sp: generation_time, beta_eff = sp.get_kinetics_parameters() diff --git a/include/openmc/bank.h b/include/openmc/bank.h index c4e940bc877..6ce396b2ad1 100644 --- a/include/openmc/bank.h +++ b/include/openmc/bank.h @@ -26,10 +26,14 @@ extern SharedArray fission_bank; extern vector> ifp_source_delayed_group_bank; +extern vector> ifp_source_ancestor_nuclide_bank; + extern vector> ifp_source_lifetime_bank; extern vector> ifp_fission_delayed_group_bank; +extern vector> ifp_fission_ancestor_nuclide_bank; + extern vector> ifp_fission_lifetime_bank; extern vector progeny_per_particle; diff --git a/include/openmc/ifp.h b/include/openmc/ifp.h index 01904d13c94..694a16b475e 100644 --- a/include/openmc/ifp.h +++ b/include/openmc/ifp.h @@ -21,13 +21,16 @@ bool is_generation_time_or_both(); //! Resize IFP vectors //! //! \param[in,out] delayed_groups List of delayed group numbers +//! \param[in,out] ancestors List of event nuclides //! \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) +void resize_ifp_data(vector& delayed_groups, vector& ancestors, + vector& lifetimes, int64_t n) { if (is_beta_effective_or_both()) { delayed_groups.resize(n); + ancestors.resize(n); } if (is_generation_time_or_both()) { lifetimes.resize(n); @@ -84,9 +87,10 @@ void resize_simulation_ifp_banks(); //! //! \param[in] i_bank Index in the fission banks //! \param[in,out] delayed_groups Delayed group numbers +//! \param[in,out] ancestors List of Ancestor nuclides //! \param[in,out] lifetimes Lifetimes lists -void copy_ifp_data_from_fission_banks( - int i_bank, vector& delayed_groups, vector& lifetimes); +void copy_ifp_data_from_fission_banks(int i_bank, vector& delayed_groups, + vector& ancestors, vector& lifetimes); #ifdef OPENMC_MPI @@ -101,9 +105,11 @@ struct DeserializationInfo { //! //! \param[in] n_generation Number of generations //! \param[in] delayed_groups List of delayed group numbers lists +//! \param[in] ancestors List of Ancestor nuclide //! \param[in] lifetimes List of lifetimes lists void broadcast_ifp_n_generation(int& n_generation, const vector>& delayed_groups, + const vector>& ancestors, const vector>& lifetimes); //! Send IFP data using MPI. @@ -115,11 +121,14 @@ void broadcast_ifp_n_generation(int& n_generation, //! \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] ancestors List of event nuclide lists +//! \param[out] send_ancestors Event nuclide 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_delayed_groups, const vector>& ancestors, + vector& send_ancestors, const vector>& lifetimes, vector& send_lifetimes); //! Receive IFP data using MPI. @@ -130,11 +139,13 @@ void send_ifp_info(int64_t idx, int64_t n, int n_generation, int neighbor, //! \param[in] neighbor Index of the neighboring processor //! \param[in] requests MPI requests //! \param[in] delayed_groups List of delayed group numbers +//! \param[in] ancestors List of event nuclide //! \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); + vector& ancestors, vector& lifetimes, + vector& deserialization); //! Copy partial IFP data from local lists to source banks. //! @@ -142,9 +153,11 @@ void receive_ifp_data(int64_t idx, int64_t n, int n_generation, int neighbor, //! \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] ancestors List of event nuclide //! \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>& ancestors, const vector>& lifetimes); //! Deserialize IFP information received using MPI and store it in @@ -153,34 +166,40 @@ void copy_partial_ifp_data_to_source_banks(int64_t idx, int n, int64_t i_bank, //! \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] ancestors List of event nuclide //! \param[in] lifetimes List of lifetimes void deserialize_ifp_info(int n_generation, const vector& deserialization, - const vector& delayed_groups, const vector& lifetimes); + const vector& delayed_groups, const vector& ancestors, + const vector& lifetimes); #endif //! Copy IFP temporary vectors to source banks. //! //! \param[in] delayed_groups List of delayed group numbers lists +//! \param[in] ancestors List of event nuclide lists //! \param[in] lifetimes List of lifetimes lists void copy_complete_ifp_data_to_source_banks( const vector>& delayed_groups, + const vector>& ancestors, const vector>& lifetimes); //! Allocate temporary vectors for IFP data. //! //! \param[in,out] delayed_groups List of delayed group numbers lists +//! \param[in,out] ancestors List of event nuclide lists //! \param[in,out] lifetimes List of delayed group numbers lists -void allocate_temporary_vector_ifp( - vector>& delayed_groups, vector>& lifetimes); +void allocate_temporary_vector_ifp(vector>& delayed_groups, + vector>& ancestors, vector>& lifetimes); //! Copy local IFP data to IFP fission banks. //! //! \param[in] delayed_groups_ptr Pointer to delayed group numbers +//! \param[in] ancestors_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); +void copy_ifp_data_to_fission_banks(const vector* delayed_groups_ptr, + const vector* ancestors_ptr, const vector* lifetimes_ptr); } // namespace openmc diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index fdacfa765b3..8c3047290b7 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -52,6 +52,7 @@ struct SourceSite { // Extra attributes that don't show up in source written to file int parent_nuclide {-1}; + int ancestor_nuclide {-1}; int64_t parent_id; int64_t progeny_id; }; @@ -65,6 +66,7 @@ struct CollisionTrackSite { double wgt {1.0}; int event_mt {0}; int delayed_group {0}; + int ancestor_nuclide {0}; int cell_id {0}; int nuclide_id; int material_id {0}; @@ -517,6 +519,7 @@ class ParticleData : public GeometryState { int event_mt_; int delayed_group_ {0}; int parent_nuclide_ {-1}; + int ancestor_nuclide_ {-1}; int n_bank_ {0}; double bank_second_E_ {0.0}; @@ -652,8 +655,10 @@ class ParticleData : public GeometryState { const int& event_mt() const { return event_mt_; } int& delayed_group() { return delayed_group_; } // delayed group const int& delayed_group() const { return delayed_group_; } - const int& parent_nuclide() const { return parent_nuclide_; } int& parent_nuclide() { return parent_nuclide_; } // Parent nuclide + const int& parent_nuclide() const { return parent_nuclide_; } + int& ancestor_nuclide() { return ancestor_nuclide_; } // Ancestor nuclide + const int& ancestor_nuclide() const { return ancestor_nuclide_; } // Post-collision data double& bank_second_E() diff --git a/include/openmc/settings.h b/include/openmc/settings.h index b369c99fef8..4cbd1dfea35 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -71,6 +71,7 @@ extern "C" bool entropy_on; //!< calculate Shannon entropy? extern "C" bool event_based; //!< use event-based mode (instead of history-based) extern bool ifp_on; //!< Use IFP for kinetics parameters? +extern bool ifp_beta_nuclide; //!< Use IFP for beta calculation by nuclide? 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? diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 02c7784d1cc..3e9414b042e 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -30,6 +30,7 @@ class _SourceSite(Structure): ('surf_id', c_int), ('particle', c_int), ('parent_nuclide', c_int), + ('ancestor_nuclide', c_int), ('parent_id', c_int64), ('progeny_id', c_int64)] diff --git a/openmc/model/model.py b/openmc/model/model.py index a9aaa481d5b..e2060958fab 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -235,13 +235,13 @@ def _get_all_materials(self) -> dict[int, openmc.Material]: return materials - def add_kinetics_parameters_tallies(self, num_groups: int | None = None): + def add_kinetics_parameters_tallies(self, num_groups: int | None = None, nuclides: Sequence[str] | None = None): """Add tallies for calculating kinetics parameters using the IFP method. This method adds tallies to the model for calculating two kinetics parameters, the generation time and the effective delayed neutron fraction (beta effective). After a model is run, these parameters can be - determined through the :meth:`openmc.StatePoint.ifp_results` method. + determined through the :meth:`openmc.StatePoint.get_kinetics_parameters` method. Parameters ---------- @@ -249,6 +249,9 @@ def add_kinetics_parameters_tallies(self, num_groups: int | None = None): Number of precursor groups to filter the delayed neutron fraction. If None, only the total effective delayed neutron fraction is tallied. + nuclides : int, optional + Nuclides to calculate separate kinetic parameters for. + If None, do not separate kinetic parameters per nuclide. """ if not any('ifp-time-numerator' in t.scores for t in self.tallies): @@ -260,6 +263,8 @@ def add_kinetics_parameters_tallies(self, num_groups: int | None = None): beta_tally.scores = ['ifp-beta-numerator'] if num_groups is not None: beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))] + if nuclides: + beta_tally.nuclides = list(nuclides) self.tallies.append(beta_tally) if not any('ifp-denominator' in t.scores for t in self.tallies): denom_tally = openmc.Tally(name='IFP denominator') diff --git a/openmc/statepoint.py b/openmc/statepoint.py index a10ec3a839d..2ad87a01c4d 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -761,21 +761,17 @@ def get_ufloat(tally, score): return uarray(tally.get_values(scores=[score]), tally.get_values(scores=[score], value='std_dev')) - denom_values = get_ufloat(denom_tally, 'ifp-denominator') + denom_value = get_ufloat(denom_tally, 'ifp-denominator').squeeze() if gen_time_tally is None: generation_time = None else: - gen_time_values = get_ufloat(gen_time_tally, 'ifp-time-numerator') - gen_time_values /= denom_values*self.keff - generation_time = gen_time_values.flatten()[0] + generation_time = get_ufloat(gen_time_tally, 'ifp-time-numerator').squeeze() + generation_time /= denom_value*self.keff if beta_tally is None: beta_effective = None else: - beta_values = get_ufloat(beta_tally, 'ifp-beta-numerator') - beta_values /= denom_values - beta_effective = beta_values.flatten() - if beta_effective.size == 1: - beta_effective = beta_effective[0] + beta_effective = get_ufloat(beta_tally, 'ifp-beta-numerator').squeeze() + beta_effective /= denom_value return KineticsParameters(generation_time, beta_effective) diff --git a/src/bank.cpp b/src/bank.cpp index 33790379b85..c495be353eb 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -31,10 +31,14 @@ SharedArray fission_bank; vector> ifp_source_delayed_group_bank; +vector> ifp_source_ancestor_nuclide_bank; + vector> ifp_source_lifetime_bank; vector> ifp_fission_delayed_group_bank; +vector> ifp_fission_ancestor_nuclide_bank; + vector> ifp_fission_lifetime_bank; // Each entry in this vector corresponds to the number of progeny produced @@ -56,8 +60,10 @@ void free_memory_bank() simulation::fission_bank.clear(); simulation::progeny_per_particle.clear(); simulation::ifp_source_delayed_group_bank.clear(); + simulation::ifp_source_ancestor_nuclide_bank.clear(); simulation::ifp_source_lifetime_bank.clear(); simulation::ifp_fission_delayed_group_bank.clear(); + simulation::ifp_fission_ancestor_nuclide_bank.clear(); simulation::ifp_fission_lifetime_bank.clear(); } @@ -92,6 +98,7 @@ void sort_fission_bank() SourceSite* sorted_bank; vector sorted_bank_holder; vector> sorted_ifp_delayed_group_bank; + vector> sorted_ifp_ancestor_nuclide_bank; vector> sorted_ifp_lifetime_bank; // If there is not enough space, allocate a temporary vector and point to it @@ -104,8 +111,8 @@ void sort_fission_bank() } if (settings::ifp_on) { - allocate_temporary_vector_ifp( - sorted_ifp_delayed_group_bank, sorted_ifp_lifetime_bank); + allocate_temporary_vector_ifp(sorted_ifp_delayed_group_bank, + sorted_ifp_ancestor_nuclide_bank, sorted_ifp_lifetime_bank); } // Use parent and progeny indices to sort fission bank @@ -119,8 +126,8 @@ void sort_fission_bank() } 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_ifp_data_from_fission_banks(i, sorted_ifp_delayed_group_bank[idx], + sorted_ifp_ancestor_nuclide_bank[idx], sorted_ifp_lifetime_bank[idx]); } } @@ -128,8 +135,8 @@ void sort_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()); + copy_ifp_data_to_fission_banks(sorted_ifp_delayed_group_bank.data(), + sorted_ifp_ancestor_nuclide_bank.data(), sorted_ifp_lifetime_bank.data()); } } diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 8412cbd3b47..6ecf60bf47d 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -137,10 +137,11 @@ void synchronize_bank() // Temporary banks for IFP vector> temp_delayed_groups; + vector> temp_ancestors; vector> temp_lifetimes; if (settings::ifp_on) { - resize_ifp_data( - temp_delayed_groups, temp_lifetimes, 3 * simulation::work_per_rank); + resize_ifp_data(temp_delayed_groups, temp_ancestors, temp_lifetimes, + 3 * simulation::work_per_rank); } // ========================================================================== @@ -168,8 +169,8 @@ void synchronize_bank() int64_t idx = std::floor(tooth) - start; temp_sites[index_temp] = simulation::fission_bank[idx]; if (settings::ifp_on) { - copy_ifp_data_from_fission_banks( - idx, temp_delayed_groups[index_temp], temp_lifetimes[index_temp]); + copy_ifp_data_from_fission_banks(idx, temp_delayed_groups[index_temp], + temp_ancestors[index_temp], temp_lifetimes[index_temp]); } ++index_temp; @@ -211,7 +212,7 @@ void synchronize_bank() int ifp_n_generation; if (settings::ifp_on) { broadcast_ifp_n_generation( - ifp_n_generation, temp_delayed_groups, temp_lifetimes); + ifp_n_generation, temp_delayed_groups, temp_ancestors, temp_lifetimes); } int64_t index_local = 0; @@ -219,6 +220,7 @@ void synchronize_bank() // IFP send buffers vector send_delayed_groups; + vector send_ancestors; vector send_lifetimes; if (start < settings::n_particles) { @@ -229,7 +231,7 @@ void synchronize_bank() // Resize IFP send buffers if (settings::ifp_on && mpi::n_procs > 1) { - resize_ifp_data(send_delayed_groups, send_lifetimes, + resize_ifp_data(send_delayed_groups, send_ancestors, send_lifetimes, ifp_n_generation * 3 * simulation::work_per_rank); } @@ -249,8 +251,8 @@ void synchronize_bank() 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); + temp_delayed_groups, send_delayed_groups, temp_ancestors, + send_ancestors, temp_lifetimes, send_lifetimes); } } @@ -275,6 +277,7 @@ void synchronize_bank() // IFP receive buffers vector recv_delayed_groups; + vector recv_ancestors; vector recv_lifetimes; vector deserialization_info; @@ -291,7 +294,7 @@ void synchronize_bank() // Resize IFP receive buffers if (settings::ifp_on && mpi::n_procs > 1) { - resize_ifp_data(recv_delayed_groups, recv_lifetimes, + resize_ifp_data(recv_delayed_groups, recv_ancestors, recv_lifetimes, ifp_n_generation * simulation::work_per_rank); } @@ -317,7 +320,8 @@ void synchronize_bank() 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); + recv_delayed_groups, recv_ancestors, recv_lifetimes, + deserialization_info); } } else { @@ -329,8 +333,8 @@ void synchronize_bank() &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); + copy_partial_ifp_data_to_source_banks(index_temp, n, index_local, + temp_delayed_groups, temp_ancestors, temp_lifetimes); } } @@ -349,14 +353,15 @@ void synchronize_bank() if (settings::ifp_on) { deserialize_ifp_info(ifp_n_generation, deserialization_info, - recv_delayed_groups, recv_lifetimes); + recv_delayed_groups, recv_ancestors, 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); + copy_complete_ifp_data_to_source_banks( + temp_delayed_groups, temp_ancestors, temp_lifetimes); } #endif diff --git a/src/ifp.cpp b/src/ifp.cpp index cc4a76538b1..af3a92f116e 100644 --- a/src/ifp.cpp +++ b/src/ifp.cpp @@ -35,6 +35,12 @@ void ifp(const Particle& p, int64_t idx) simulation::ifp_source_delayed_group_bank[p.current_work() - 1]; simulation::ifp_fission_delayed_group_bank[idx] = _ifp(p.delayed_group(), delayed_groups); + if (settings::ifp_beta_nuclide) { + const auto& ancestor_nuclides = + simulation::ifp_source_ancestor_nuclide_bank[p.current_work() - 1]; + simulation::ifp_fission_ancestor_nuclide_bank[idx] = + _ifp(p.event_nuclide(), ancestor_nuclides); + } } if (is_generation_time_or_both()) { const auto& lifetimes = @@ -46,16 +52,19 @@ void ifp(const Particle& p, int64_t idx) void resize_simulation_ifp_banks() { resize_ifp_data(simulation::ifp_source_delayed_group_bank, + simulation::ifp_source_ancestor_nuclide_bank, simulation::ifp_source_lifetime_bank, simulation::work_per_rank); resize_ifp_data(simulation::ifp_fission_delayed_group_bank, + simulation::ifp_fission_ancestor_nuclide_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) +void copy_ifp_data_from_fission_banks(int i_bank, vector& delayed_groups, + vector& ancestors, vector& lifetimes) { if (is_beta_effective_or_both()) { delayed_groups = simulation::ifp_fission_delayed_group_bank[i_bank]; + ancestors = simulation::ifp_fission_ancestor_nuclide_bank[i_bank]; } if (is_generation_time_or_both()) { lifetimes = simulation::ifp_fission_lifetime_bank[i_bank]; @@ -66,7 +75,7 @@ void copy_ifp_data_from_fission_banks( void broadcast_ifp_n_generation(int& n_generation, const vector>& delayed_groups, - const vector>& lifetimes) + const vector>& ancestors, const vector>& lifetimes) { if (mpi::rank == 0) { if (is_beta_effective_or_both()) { @@ -80,7 +89,8 @@ void broadcast_ifp_n_generation(int& n_generation, 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_delayed_groups, const vector>& ancestors, + vector& send_ancestors, const vector>& lifetimes, vector& send_lifetimes) { // Copy data in send buffers @@ -88,6 +98,8 @@ void send_ifp_info(int64_t idx, int64_t n, int n_generation, int neighbor, if (is_beta_effective_or_both()) { std::copy(delayed_groups[i].begin(), delayed_groups[i].end(), send_delayed_groups.begin() + i * n_generation); + std::copy(ancestors[i].begin(), ancestors[i].end(), + send_ancestors.begin() + i * n_generation); } if (is_generation_time_or_both()) { std::copy(lifetimes[i].begin(), lifetimes[i].end(), @@ -100,6 +112,9 @@ void send_ifp_info(int64_t idx, int64_t n, int n_generation, int neighbor, MPI_Isend(&send_delayed_groups[n_generation * idx], n_generation * static_cast(n), MPI_INT, neighbor, mpi::rank, mpi::intracomm, &requests.back()); + MPI_Isend(&send_ancestors[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()) { @@ -112,7 +127,8 @@ void send_ifp_info(int64_t idx, int64_t n, int n_generation, int neighbor, void receive_ifp_data(int64_t idx, int64_t n, int n_generation, int neighbor, vector& requests, vector& delayed_groups, - vector& lifetimes, vector& deserialization) + vector& ancestors, vector& lifetimes, + vector& deserialization) { // Receive delayed groups if (is_beta_effective_or_both()) { @@ -120,6 +136,9 @@ void receive_ifp_data(int64_t idx, int64_t n, int n_generation, int neighbor, MPI_Irecv(&delayed_groups[n_generation * idx], n_generation * static_cast(n), MPI_INT, neighbor, neighbor, mpi::intracomm, &requests.back()); + MPI_Irecv(&ancestors[n_generation * idx], + n_generation * static_cast(n), MPI_INT, neighbor, neighbor, + mpi::intracomm, &requests.back()); } // Receive lifetimes if (is_generation_time_or_both()) { @@ -135,11 +154,13 @@ void receive_ifp_data(int64_t idx, int64_t n, int n_generation, int neighbor, void copy_partial_ifp_data_to_source_banks(int64_t idx, int n, int64_t i_bank, const vector>& delayed_groups, - const vector>& lifetimes) + const vector>& ancestors, 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]); + std::copy(&ancestors[idx], &ancestors[idx + n], + &simulation::ifp_source_ancestor_nuclide_bank[i_bank]); } if (is_generation_time_or_both()) { std::copy(&lifetimes[idx], &lifetimes[idx + n], @@ -149,7 +170,8 @@ void copy_partial_ifp_data_to_source_banks(int64_t idx, int n, int64_t i_bank, void deserialize_ifp_info(int n_generation, const vector& deserialization, - const vector& delayed_groups, const vector& lifetimes) + const vector& delayed_groups, const vector& ancestors, + const vector& lifetimes) { for (auto info : deserialization) { int64_t index_local = info.index_local; @@ -161,6 +183,9 @@ void deserialize_ifp_info(int n_generation, delayed_groups.begin() + n_generation * i, delayed_groups.begin() + n_generation * (i + 1)); simulation::ifp_source_delayed_group_bank[i] = delayed_groups_received; + vector ancestors_received(ancestors.begin() + n_generation * i, + ancestors.begin() + n_generation * (i + 1)); + simulation::ifp_source_ancestor_nuclide_bank[i] = ancestors_received; } if (is_generation_time_or_both()) { vector lifetimes_received(lifetimes.begin() + n_generation * i, @@ -175,12 +200,14 @@ void deserialize_ifp_info(int n_generation, void copy_complete_ifp_data_to_source_banks( const vector>& delayed_groups, - const vector>& lifetimes) + const vector>& ancestors, 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()); + std::copy(ancestors.data(), ancestors.data() + settings::n_particles, + simulation::ifp_source_ancestor_nuclide_bank.begin()); } if (is_generation_time_or_both()) { std::copy(lifetimes.data(), lifetimes.data() + settings::n_particles, @@ -188,11 +215,12 @@ void copy_complete_ifp_data_to_source_banks( } } -void allocate_temporary_vector_ifp( - vector>& delayed_groups, vector>& lifetimes) +void allocate_temporary_vector_ifp(vector>& delayed_groups, + vector>& ancestors, vector>& lifetimes) { if (is_beta_effective_or_both()) { delayed_groups.resize(simulation::fission_bank.size()); + ancestors.resize(simulation::fission_bank.size()); } if (is_generation_time_or_both()) { lifetimes.resize(simulation::fission_bank.size()); @@ -200,12 +228,14 @@ void allocate_temporary_vector_ifp( } void copy_ifp_data_to_fission_banks(const vector* const delayed_groups_ptr, - const vector* lifetimes_ptr) + const vector* const ancestors_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()); + std::copy(ancestors_ptr, ancestors_ptr + simulation::fission_bank.size(), + simulation::ifp_fission_ancestor_nuclide_bank.data()); } if (is_generation_time_or_both()) { std::copy(lifetimes_ptr, lifetimes_ptr + simulation::fission_bank.size(), diff --git a/src/initialize.cpp b/src/initialize.cpp index e2a5b974338..17cabad0adf 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -157,7 +157,7 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype SourceSite b; - MPI_Aint disp[11]; + MPI_Aint disp[12]; MPI_Get_address(&b.r, &disp[0]); MPI_Get_address(&b.u, &disp[1]); MPI_Get_address(&b.E, &disp[2]); @@ -167,16 +167,18 @@ void initialize_mpi(MPI_Comm intracomm) MPI_Get_address(&b.surf_id, &disp[6]); MPI_Get_address(&b.particle, &disp[7]); MPI_Get_address(&b.parent_nuclide, &disp[8]); - MPI_Get_address(&b.parent_id, &disp[9]); - MPI_Get_address(&b.progeny_id, &disp[10]); - for (int i = 10; i >= 0; --i) { + MPI_Get_address(&b.ancestor_nuclide, &disp[9]); + MPI_Get_address(&b.parent_id, &disp[10]); + MPI_Get_address(&b.progeny_id, &disp[11]); + for (int i = 11; i >= 0; --i) { disp[i] -= disp[0]; } - int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, - MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; - MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site); + MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, + MPI_LONG}; + MPI_Type_create_struct(12, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); CollisionTrackSite bc; diff --git a/src/particle.cpp b/src/particle.cpp index 2d70d715e22..7b004d061fb 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -148,6 +148,7 @@ void Particle::from_source(const SourceSite* src) time() = src->time; time_last() = src->time; parent_nuclide() = src->parent_nuclide; + ancestor_nuclide() = src->ancestor_nuclide; delayed_group() = src->delayed_group; // Convert signed surface ID to signed index diff --git a/src/settings.cpp b/src/settings.cpp index d47e9b5e6bd..54059418b5e 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -56,6 +56,7 @@ bool delayed_photon_scaling {true}; bool entropy_on {false}; bool event_based {false}; bool ifp_on {false}; +bool ifp_beta_nuclide {false}; bool legendre_to_tabular {true}; bool material_cell_offsets {true}; bool output_summary {true}; diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 6eef1da9cfd..134a3b4196d 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -235,6 +235,10 @@ Tally::Tally(pugi::xml_node node) } break; case SCORE_IFP_BETA_NUM: + // This bool indicates if the ancestor_nuclide vector has to be filled + if (nuclides_[0] > -1) { + settings::ifp_beta_nuclide = true; + } case SCORE_IFP_DENOM: if (settings::ifp_parameter == IFPParameter::None) { settings::ifp_parameter = IFPParameter::BetaEffective; diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 67e851644a9..7dc237d0c52 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -962,14 +962,23 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, 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) { + int j = -1; + if (delayed_groups[0] > 0 && tally.nuclides_[0] == -1) + j = 0; + if (delayed_groups[1] > 0 && tally.nuclides_[0] > -1) { + const auto& ancestor_event_nuclide = simulation:: + ifp_source_ancestor_nuclide_bank[p.current_work() - 1]; + if (ancestor_event_nuclide[0] == i_nuclide) + j = 1; + } + if (j > -1) { score = p.wgt_last(); if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; const DelayedGroupFilter& filt { *dynamic_cast( model::tally_filters[i_dg_filt].get())}; - score_fission_delayed_dg(i_tally, delayed_groups[0] - 1, + score_fission_delayed_dg(i_tally, delayed_groups[j] - 1, score, score_index, p.filter_matches()); continue; } diff --git a/tests/regression_tests/ifp/nuclidewise/__init__.py b/tests/regression_tests/ifp/nuclidewise/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/ifp/nuclidewise/inputs_true.dat b/tests/regression_tests/ifp/nuclidewise/inputs_true.dat new file mode 100644 index 00000000000..68fc5888589 --- /dev/null +++ b/tests/regression_tests/ifp/nuclidewise/inputs_true.dat @@ -0,0 +1,45 @@ + + + + + + + + + + + + + + + eigenvalue + 100000 + 20 + 5 + + + -10.0 -10.0 -10.0 10.0 10.0 10.0 + + + true + + + 5 + + + + 1 2 3 4 5 6 + + + ifp-time-numerator + + + 1 + U235 Pu239 + ifp-beta-numerator + + + ifp-denominator + + + diff --git a/tests/regression_tests/ifp/nuclidewise/results_true.dat b/tests/regression_tests/ifp/nuclidewise/results_true.dat new file mode 100644 index 00000000000..d3fa8333047 --- /dev/null +++ b/tests/regression_tests/ifp/nuclidewise/results_true.dat @@ -0,0 +1,33 @@ +k-combined: +1.256195E+00 5.261787E-04 +tally 1: +7.942195E-08 +4.207211E-16 +tally 2: +1.710000E-03 +4.295000E-07 +9.400000E-04 +1.794000E-07 +5.970000E-03 +2.607900E-06 +4.550000E-03 +1.625100E-06 +5.070000E-03 +2.047900E-06 +2.600000E-03 +5.692000E-07 +1.236000E-02 +1.062980E-05 +4.020000E-03 +1.379200E-06 +5.690000E-03 +3.302300E-06 +3.180000E-03 +8.316000E-07 +2.100000E-03 +4.486000E-07 +1.300000E-03 +2.008000E-07 +tally 3: +1.495260E+01 +1.490559E+01 diff --git a/tests/regression_tests/ifp/nuclidewise/test.py b/tests/regression_tests/ifp/nuclidewise/test.py new file mode 100644 index 00000000000..05251b3c31c --- /dev/null +++ b/tests/regression_tests/ifp/nuclidewise/test.py @@ -0,0 +1,41 @@ +"""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(): + # Material + material = openmc.Material(name="core") + material.add_nuclide("U235", 1.0) + material.add_nuclide("Pu239", 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) + geometry = openmc.Geometry([cell]) + + # Settings + settings = openmc.Settings() + settings.particles = 100000 + settings.batches = 20 + settings.inactive = 5 + settings.ifp_n_generation = 5 + + model = openmc.Model(settings=settings, geometry=geometry) + + space = openmc.stats.Box(*cell.bounding_box) + model.settings.source = openmc.IndependentSource( + space=space, constraints={'fissionable': True}) + model.add_kinetics_parameters_tallies(num_groups=6, nuclides = ["U235", "Pu239"]) + 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 index e527f162460..5952eaf1efd 100644 --- a/tests/unit_tests/test_ifp.py +++ b/tests/unit_tests/test_ifp.py @@ -21,6 +21,7 @@ def geometry(): openmc.reset_auto_ids() material = openmc.Material() material.add_nuclide("U235", 1.0) + material.add_nuclide("Pu239", 1.0) sphere = openmc.Sphere(r=1.0, boundary_type="vacuum") cell = openmc.Cell(region=-sphere, fill=material) return openmc.Geometry([cell]) @@ -50,15 +51,19 @@ def test_exceptions(options, error, run_in_tmpdir, geometry): @pytest.mark.parametrize( - "num_groups, use_auto_tallies", + "num_groups, nuclides, use_auto_tallies", [ - (None, True), - (None, False), - (6, True), - (6, False), + (None, None, True), + (None, None, False), + (6, None, True), + (6, None, False), + (None, ["U235", "Pu239"], True), + (None, ["U235", "Pu239"], False), + (6, ["U235", "Pu239"], True), + (6, ["U235", "Pu239"], False), ], ) -def test_get_kinetics_parameters(run_in_tmpdir, geometry, num_groups, use_auto_tallies): +def test_get_kinetics_parameters(run_in_tmpdir, geometry, num_groups, nuclides, use_auto_tallies): # Create basic model model = openmc.Model(geometry=geometry) model.settings.particles = 1000 @@ -68,13 +73,16 @@ def test_get_kinetics_parameters(run_in_tmpdir, geometry, num_groups, use_auto_t # Add IFP tallies either via the convenience method or manually if use_auto_tallies: - model.add_kinetics_parameters_tallies(num_groups=num_groups) + model.add_kinetics_parameters_tallies(num_groups=num_groups, nuclides=nuclides) else: for score in ["ifp-time-numerator", "ifp-beta-numerator", "ifp-denominator"]: tally = openmc.Tally() tally.scores = [score] - if score == "ifp-beta-numerator" and num_groups is not None: - tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))] + if score == "ifp-beta-numerator": + if num_groups is not None: + tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))] + if nuclides: + tally.nuclides = nuclides model.tallies.append(tally) # Run and get kinetics parameters @@ -85,4 +93,9 @@ def test_get_kinetics_parameters(run_in_tmpdir, geometry, num_groups, use_auto_t assert params.generation_time is not None assert params.beta_effective is not None if num_groups is not None: - assert len(params.beta_effective) == num_groups + if nuclides: + assert params.beta_effective.shape == (num_groups, len(nuclides)) + else: + assert len(params.beta_effective) == num_groups + elif nuclides: + assert len(params.beta_effective) == len(nuclides)