From 50195517b49bbbe3a648522e191997c74f281531 Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 26 Jan 2026 23:59:37 +0200 Subject: [PATCH 01/24] resolve conflicts with weightwindows and survival biasing --- include/openmc/physics_common.h | 4 ++++ src/physics.cpp | 15 +++------------ src/physics_common.cpp | 16 ++++++++++++++++ src/physics_mg.cpp | 15 +++------------ src/weight_windows.cpp | 8 ++++++-- 5 files changed, 32 insertions(+), 26 deletions(-) diff --git a/include/openmc/physics_common.h b/include/openmc/physics_common.h index e38a3c7f883..9822b367fcd 100644 --- a/include/openmc/physics_common.h +++ b/include/openmc/physics_common.h @@ -13,5 +13,9 @@ namespace openmc { //! \param[in] weight_survive Weight assigned to particles that survive void russian_roulette(Particle& p, double weight_survive); +//! \brief Performs the global russian roulette operation for a neutron +//! \param[in,out] p Particle object +void apply_neutron_russian_roulette(Particle& p); + } // namespace openmc #endif // OPENMC_PHYSICS_COMMON_H diff --git a/src/physics.cpp b/src/physics.cpp index 41509af97be..bc8723d99d9 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -153,18 +153,9 @@ void sample_neutron_reaction(Particle& p) advance_prn_seed(data::nuclides.size(), &p.seeds(STREAM_URR_PTABLE)); } - // Play russian roulette if survival biasing is turned on - if (settings::survival_biasing) { - // if survival normalization is on, use normalized weight cutoff and - // normalized weight survive - if (settings::survival_normalization) { - if (p.wgt() < settings::weight_cutoff * p.wgt_born()) { - russian_roulette(p, settings::weight_survive * p.wgt_born()); - } - } else if (p.wgt() < settings::weight_cutoff) { - russian_roulette(p, settings::weight_survive); - } - } + // Play russian roulette if neutrons have no weight windows + if (!settings::weight_window_checkpoint_collision) + apply_neutron_russian_roulette(p); } void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) diff --git a/src/physics_common.cpp b/src/physics_common.cpp index e25ae6b97ab..7001a03c7c8 100644 --- a/src/physics_common.cpp +++ b/src/physics_common.cpp @@ -18,4 +18,20 @@ void russian_roulette(Particle& p, double weight_survive) } } +void apply_neutron_russian_roulette(Particle& p) +{ + // Play russian roulette if survival biasing is turned on + if (settings::survival_biasing) { + // if survival normalization is on, use normalized weight cutoff and + // normalized weight survive + if (settings::survival_normalization) { + if (p.wgt() < settings::weight_cutoff * p.wgt_born()) { + russian_roulette(p, settings::weight_survive * p.wgt_born()); + } + } else if (p.wgt() < settings::weight_cutoff) { + russian_roulette(p, settings::weight_survive); + } + } +} + } // namespace openmc diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 4c28cb1795a..9e2aceb60c6 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -66,18 +66,9 @@ void sample_reaction(Particle& p) // Sample a scattering event to determine the energy of the exiting neutron scatter(p); - // Play Russian roulette if survival biasing is turned on - if (settings::survival_biasing) { - // if survival normalization is applicable, use normalized weight cutoff and - // normalized weight survive - if (settings::survival_normalization) { - if (p.wgt() < settings::weight_cutoff * p.wgt_born()) { - russian_roulette(p, settings::weight_survive * p.wgt_born()); - } - } else if (p.wgt() < settings::weight_cutoff) { - russian_roulette(p, settings::weight_survive); - } - } + // Play russian roulette if neutrons have no weight windows + if (!settings::weight_window_checkpoint_collision) + apply_neutron_russian_roulette(p); } void scatter(Particle& p) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 4838e459152..e203fdd414e 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -86,9 +86,13 @@ void apply_weight_windows(Particle& p) } } - // particle is not in any of the ww domains, do nothing - if (!weight_window.is_valid()) + // If particle is not in any of the ww domains + if (!weight_window.is_valid()) { + // If particle is a neutron do russian roulette. + if (p.type() == ParticleType::neutron) + apply_neutron_russian_roulette(p); return; + } // Normalize weight windows based on particle's starting weight // and the value of the weight window the particle was born in. From eb851da6d582bb216b7336e68cdd3fc94cd1c3b7 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 27 Jan 2026 00:05:59 +0200 Subject: [PATCH 02/24] fix comments --- src/physics.cpp | 2 +- src/physics_mg.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index bc8723d99d9..f37f9ea29f3 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -153,7 +153,7 @@ void sample_neutron_reaction(Particle& p) advance_prn_seed(data::nuclides.size(), &p.seeds(STREAM_URR_PTABLE)); } - // Play russian roulette if neutrons have no weight windows + // Play russian roulette if there are no weight windows if (!settings::weight_window_checkpoint_collision) apply_neutron_russian_roulette(p); } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 9e2aceb60c6..864c8419b98 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -66,7 +66,7 @@ void sample_reaction(Particle& p) // Sample a scattering event to determine the energy of the exiting neutron scatter(p); - // Play russian roulette if neutrons have no weight windows + // Play russian roulette if there are no weight windows if (!settings::weight_window_checkpoint_collision) apply_neutron_russian_roulette(p); } From 29801645cce4ccb2694ec17e022eba4d0651bb9d Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 27 Jan 2026 10:11:57 +0200 Subject: [PATCH 03/24] fix ww bug --- src/physics.cpp | 3 ++- src/physics_mg.cpp | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index f37f9ea29f3..a355a4d9899 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -154,7 +154,8 @@ void sample_neutron_reaction(Particle& p) } // Play russian roulette if there are no weight windows - if (!settings::weight_window_checkpoint_collision) + if (!settings::weight_windows_on || + !settings::weight_window_checkpoint_collision) apply_neutron_russian_roulette(p); } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 864c8419b98..8b3263cfe30 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -67,7 +67,8 @@ void sample_reaction(Particle& p) scatter(p); // Play russian roulette if there are no weight windows - if (!settings::weight_window_checkpoint_collision) + if (!settings::weight_windows_on || + !settings::weight_window_checkpoint_collision) apply_neutron_russian_roulette(p); } From a3ba86f01773a5663e69a55aa312727e04b70cba Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 27 Jan 2026 23:17:22 +0200 Subject: [PATCH 04/24] apply PR suggestions --- include/openmc/physics_common.h | 4 ++-- src/physics.cpp | 2 +- src/physics_common.cpp | 24 ++++++++++++------------ src/physics_mg.cpp | 2 +- src/weight_windows.cpp | 2 +- 5 files changed, 17 insertions(+), 17 deletions(-) diff --git a/include/openmc/physics_common.h b/include/openmc/physics_common.h index 9822b367fcd..b0e395c1ea4 100644 --- a/include/openmc/physics_common.h +++ b/include/openmc/physics_common.h @@ -13,9 +13,9 @@ namespace openmc { //! \param[in] weight_survive Weight assigned to particles that survive void russian_roulette(Particle& p, double weight_survive); -//! \brief Performs the global russian roulette operation for a neutron +//! \brief Performs the global russian roulette operation on a particle //! \param[in,out] p Particle object -void apply_neutron_russian_roulette(Particle& p); +void apply_russian_roulette(Particle& p); } // namespace openmc #endif // OPENMC_PHYSICS_COMMON_H diff --git a/src/physics.cpp b/src/physics.cpp index a355a4d9899..ece201dab47 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -156,7 +156,7 @@ void sample_neutron_reaction(Particle& p) // Play russian roulette if there are no weight windows if (!settings::weight_windows_on || !settings::weight_window_checkpoint_collision) - apply_neutron_russian_roulette(p); + apply_russian_roulette(p); } void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) diff --git a/src/physics_common.cpp b/src/physics_common.cpp index 7001a03c7c8..a40c6f20ce9 100644 --- a/src/physics_common.cpp +++ b/src/physics_common.cpp @@ -18,20 +18,20 @@ void russian_roulette(Particle& p, double weight_survive) } } -void apply_neutron_russian_roulette(Particle& p) +void apply_russian_roulette(Particle& p) { - // Play russian roulette if survival biasing is turned on - if (settings::survival_biasing) { - // if survival normalization is on, use normalized weight cutoff and - // normalized weight survive - if (settings::survival_normalization) { - if (p.wgt() < settings::weight_cutoff * p.wgt_born()) { - russian_roulette(p, settings::weight_survive * p.wgt_born()); - } - } else if (p.wgt() < settings::weight_cutoff) { - russian_roulette(p, settings::weight_survive); + // Exit if survival biasing is turned off + if (settings::survival_biasing) + return; + + // if survival normalization is on, use normalized weight cutoff and + // normalized weight survive + if (settings::survival_normalization) { + if (p.wgt() < settings::weight_cutoff * p.wgt_born()) { + russian_roulette(p, settings::weight_survive * p.wgt_born()); } + } else if (p.wgt() < settings::weight_cutoff) { + russian_roulette(p, settings::weight_survive); } } - } // namespace openmc diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 8b3263cfe30..9ff725a3214 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -69,7 +69,7 @@ void sample_reaction(Particle& p) // Play russian roulette if there are no weight windows if (!settings::weight_windows_on || !settings::weight_window_checkpoint_collision) - apply_neutron_russian_roulette(p); + apply_russian_roulette(p); } void scatter(Particle& p) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index e203fdd414e..c299f1a4eb1 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -90,7 +90,7 @@ void apply_weight_windows(Particle& p) if (!weight_window.is_valid()) { // If particle is a neutron do russian roulette. if (p.type() == ParticleType::neutron) - apply_neutron_russian_roulette(p); + apply_russian_roulette(p); return; } From 928f9472316d2c932b959a9e73f270ada19584e5 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 00:43:03 +0200 Subject: [PATCH 05/24] restructured code --- include/openmc/weight_windows.h | 18 +++++++++ src/physics.cpp | 14 +++++-- src/physics_mg.cpp | 14 +++++-- src/settings.cpp | 7 ++++ src/weight_windows.cpp | 70 +++++++++++++++++---------------- 5 files changed, 82 insertions(+), 41 deletions(-) diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 7638155228e..6b79bc4ce80 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -33,9 +33,18 @@ constexpr double DEFAULT_WEIGHT_CUTOFF {1.0e-38}; // default low weight cutoff //! \param[in] p Particle to apply weight windows to void apply_weight_windows(Particle& p); +//! Apply weight window to a particle +//! \param[in] p Particle to apply weight window to +//! \param[in] ww WeightWindow to apply +void apply_weight_window(Particle& p, WeightWindow& ww); + //! Free memory associated with weight windows void free_memory_weight_windows(); +//! Search weight window that apply to a particle +//! \param[in] p Particle to search weight window for +WeightWindow* search_weight_window(Particle& p); + //============================================================================== // Global variables //============================================================================== @@ -135,10 +144,19 @@ class WeightWindows { //! \param[in] group HDF5 group to write to void to_hdf5(hid_t group) const; + //! Get particle bin inside weight window mesh + //! \param[in] p Particle to get weight window bin for + const int get_mesh_bin(const Particle& p) const; + //! Retrieve the weight window for a particle //! \param[in] p Particle to get weight window for WeightWindow get_weight_window(const Particle& p) const; + //! Retrieve the weight window for a mesh bin + //! \param[in] double E particle energy + //! \param[in] int mesh_bin mesh bin + WeightWindow get_weight_window(double E, const int mesh_bin) const; + std::array bounds_size() const; const vector& energy_bounds() const { return energy_bounds_; } diff --git a/src/physics.cpp b/src/physics.cpp index ece201dab47..8a508c042d3 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -61,8 +61,15 @@ void collision(Particle& p) break; } - if (settings::weight_window_checkpoint_collision) - apply_weight_windows(p); + if (settings::weight_windows_on) { + WeightWindow* ww = search_weight_window(p); + if (ww != nullptr) { + if (settings::weight_window_checkpoint_collision) + apply_weight_window(p, *ww); + } else if (p.type() == ParticleType::neutron) { + apply_russian_roullete(p); + } + } // Kill particle if energy falls below cutoff int type = static_cast(p.type()); @@ -154,8 +161,7 @@ void sample_neutron_reaction(Particle& p) } // Play russian roulette if there are no weight windows - if (!settings::weight_windows_on || - !settings::weight_window_checkpoint_collision) + if (!settings::weight_windows_on) apply_russian_roulette(p); } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 9ff725a3214..5936fc75ed4 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -31,8 +31,15 @@ void collision_mg(Particle& p) // Sample the reaction type sample_reaction(p); - if (settings::weight_window_checkpoint_collision) - apply_weight_windows(p); + if (settings::weight_windows_on) { + WeightWindow* ww = search_weight_window(p); + if (ww != nullptr) { + if (settings::weight_window_checkpoint_collision) + apply_weight_window(p, *ww); + } else if (p.type() == ParticleType::neutron) { + apply_russian_roullete(p); + } + } // Display information about collision if ((settings::verbosity >= 10) || p.trace()) { @@ -67,8 +74,7 @@ void sample_reaction(Particle& p) scatter(p); // Play russian roulette if there are no weight windows - if (!settings::weight_windows_on || - !settings::weight_window_checkpoint_collision) + if (!settings::weight_windows_on) apply_russian_roulette(p); } diff --git a/src/settings.cpp b/src/settings.cpp index 5b472468fc4..65790277d80 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -1264,6 +1264,13 @@ void read_settings_xml(pugi::xml_node root) } } + if (weight_windows_on) { + if (!weight_window_checkpoint_surface && + !weight_window_checkpoint_collision) + fatal_error( + "Weight Windows are enabled but there are no valid checkpoints."); + } + if (check_for_node(root, "use_decay_photons")) { settings::use_decay_photons = get_node_value_bool(root, "use_decay_photons"); diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index c299f1a4eb1..54ed73989ff 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -53,28 +53,32 @@ openmc::vector> weight_windows_generators; // Non-member functions //============================================================================== +WeightWindow* search_weight_window(Particle& p) +{ + // TODO: this is a linear search - should do something more clever + WeightWindow weight_window; + for (const auto& ww : variance_reduction::weight_windows) { + int mesh_bin = ww->get_mesh_bin(p); + if (mesh_bin > 0) + return ww->get_weight_window(p.E(), mesh_bin); + } +} + void apply_weight_windows(Particle& p) { if (!settings::weight_windows_on) return; - // WW on photon and neutron only if (p.type() != ParticleType::neutron && p.type() != ParticleType::photon) return; - // skip dead or no energy - if (p.E() <= 0 || !p.alive()) - return; - - bool in_domain = false; - // TODO: this is a linear search - should do something more clever - WeightWindow weight_window; - for (const auto& ww : variance_reduction::weight_windows) { - weight_window = ww->get_weight_window(p); - if (weight_window.is_valid()) - break; - } + WeightWindow* ww = search_weight_window(p); + if (ww != nullptr) + apply_weight_window(p, *ww); +} +void apply_weight_window(Particle& p, WeightWindow& ww) +{ // If particle has not yet had its birth weight window value set, set it to // the current weight window (or 1.0 if not born in a weight window). if (p.wgt_ww_born() == -1.0) { @@ -86,14 +90,6 @@ void apply_weight_windows(Particle& p) } } - // If particle is not in any of the ww domains - if (!weight_window.is_valid()) { - // If particle is a neutron do russian roulette. - if (p.type() == ParticleType::neutron) - apply_russian_roulette(p); - return; - } - // Normalize weight windows based on particle's starting weight // and the value of the weight window the particle was born in. weight_window.scale(p.wgt_born() / p.wgt_ww_born()); @@ -376,26 +372,34 @@ void WeightWindows::set_mesh(const Mesh* mesh) set_mesh(model::mesh_map[mesh->id_]); } -WeightWindow WeightWindows::get_weight_window(const Particle& p) const +const int WeightWindows::get_mesh_bin(const Particle& p) const { // check for particle type - if (particle_type_ != p.type()) { - return {}; - } - - // Get mesh index for particle's position - const auto& mesh = this->mesh(); - int mesh_bin = mesh->get_bin(p.r()); - - // particle is outside the weight window mesh - if (mesh_bin < 0) - return {}; + if (particle_type_ != p.type()) + return C_NONE; // particle energy double E = p.E(); // check to make sure energy is in range, expects sorted energy values if (E < energy_bounds_.front() || E > energy_bounds_.back()) + return C_NONE; + + // Get mesh index for particle's position + const auto& mesh = this->mesh(); + return mesh->get_bin(p.r()); +} + +WeightWindow WeightWindows::get_weight_window(const Particle& p) const +{ + return get_weight_window(p.E(), get_mesh_bin(p)); +} + +WeightWindow WeightWindows::get_weight_window( + const double E, const int mesh_bin) const +{ + // particle is outside the weight window mesh + if (mesh_bin < 0) return {}; // get the mesh bin in energy group From 8cdca055262231000872c6382ee2067df497b71b Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 00:49:38 +0200 Subject: [PATCH 06/24] rfix typos --- include/openmc/weight_windows.h | 40 ++++++++++++++++----------------- src/physics.cpp | 2 +- src/physics_mg.cpp | 2 +- 3 files changed, 22 insertions(+), 22 deletions(-) diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 6b79bc4ce80..d83b01a5c88 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -25,26 +25,6 @@ enum class WeightWindowUpdateMethod { MAGIC, FW_CADIS }; constexpr double DEFAULT_WEIGHT_CUTOFF {1.0e-38}; // default low weight cutoff -//============================================================================== -// Non-member functions -//============================================================================== - -//! Apply weight windows to a particle -//! \param[in] p Particle to apply weight windows to -void apply_weight_windows(Particle& p); - -//! Apply weight window to a particle -//! \param[in] p Particle to apply weight window to -//! \param[in] ww WeightWindow to apply -void apply_weight_window(Particle& p, WeightWindow& ww); - -//! Free memory associated with weight windows -void free_memory_weight_windows(); - -//! Search weight window that apply to a particle -//! \param[in] p Particle to search weight window for -WeightWindow* search_weight_window(Particle& p); - //============================================================================== // Global variables //============================================================================== @@ -255,6 +235,26 @@ class WeightWindowsGenerator { double ratio_ {5.0}; // Date: Wed, 28 Jan 2026 00:51:43 +0200 Subject: [PATCH 07/24] fix poblem --- src/weight_windows.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 54ed73989ff..3d9121ace37 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -68,10 +68,15 @@ void apply_weight_windows(Particle& p) { if (!settings::weight_windows_on) return; + // WW on photon and neutron only if (p.type() != ParticleType::neutron && p.type() != ParticleType::photon) return; + // skip dead or no energy + if (p.E() <= 0 || !p.alive()) + return; + WeightWindow* ww = search_weight_window(p); if (ww != nullptr) apply_weight_window(p, *ww); From e89f1b798cc47c95c0d92d154e2e6f47233d9d52 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 00:59:03 +0200 Subject: [PATCH 08/24] fix another problem --- include/openmc/weight_windows.h | 4 ++-- src/weight_windows.cpp | 5 ++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index d83b01a5c88..783c1854a09 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -245,8 +245,8 @@ void apply_weight_windows(Particle& p); //! Apply weight window to a particle //! \param[in] p Particle to apply weight window to -//! \param[in] ww WeightWindow to apply -void apply_weight_window(Particle& p, WeightWindow& ww); +//! \param[in] weight_window WeightWindow to apply +void apply_weight_window(Particle& p, WeightWindow& weight_window); //! Free memory associated with weight windows void free_memory_weight_windows(); diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 3d9121ace37..42c09e8c586 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -56,11 +56,10 @@ openmc::vector> weight_windows_generators; WeightWindow* search_weight_window(Particle& p) { // TODO: this is a linear search - should do something more clever - WeightWindow weight_window; for (const auto& ww : variance_reduction::weight_windows) { int mesh_bin = ww->get_mesh_bin(p); if (mesh_bin > 0) - return ww->get_weight_window(p.E(), mesh_bin); + return &(ww->get_weight_window(p.E(), mesh_bin)); } } @@ -82,7 +81,7 @@ void apply_weight_windows(Particle& p) apply_weight_window(p, *ww); } -void apply_weight_window(Particle& p, WeightWindow& ww) +void apply_weight_window(Particle& p, WeightWindow& weight_window) { // If particle has not yet had its birth weight window value set, set it to // the current weight window (or 1.0 if not born in a weight window). From 11b3f41f916fa253d827ef8c4d41ee344804769d Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 01:22:45 +0200 Subject: [PATCH 09/24] fix another problem --- include/openmc/weight_windows.h | 7 +- src/physics.cpp | 5 +- src/physics_mg.cpp | 5 +- src/weight_windows.cpp | 210 ++++++++++++++++---------------- 4 files changed, 114 insertions(+), 113 deletions(-) diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 783c1854a09..bf49e5ce33f 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -246,14 +246,15 @@ void apply_weight_windows(Particle& p); //! Apply weight window to a particle //! \param[in] p Particle to apply weight window to //! \param[in] weight_window WeightWindow to apply -void apply_weight_window(Particle& p, WeightWindow& weight_window); +void apply_weight_window(Particle& p, WeightWindow weight_window); //! Free memory associated with weight windows void free_memory_weight_windows(); //! Search weight window that apply to a particle -//! \param[in] p Particle to search weight window for -WeightWindow* search_weight_window(Particle& p); +//! \param[in] p Particle to search weight window for +//! \param[out] mesh_bin int mesh bin inside WeightWindows +const WeightWindows* search_weight_window(Particle& p, int& mesh_bin); //! Finalize variance reduction objects after all inputs have been read void finalize_variance_reduction(); diff --git a/src/physics.cpp b/src/physics.cpp index f40f96a87b7..16a1b5706b0 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -62,10 +62,11 @@ void collision(Particle& p) } if (settings::weight_windows_on) { - WeightWindow* ww = search_weight_window(p); + int mesh_bin; + const auto& ww = search_weight_window(p, mesh_bin); if (ww != nullptr) { if (settings::weight_window_checkpoint_collision) - apply_weight_window(p, *ww); + apply_weight_window(p, ww->get_weight_window(p.E(), mesh_bin)); } else if (p.type() == ParticleType::neutron) { apply_russian_roulette(p); } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index c1034a95181..4549dcffcd6 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -32,10 +32,11 @@ void collision_mg(Particle& p) sample_reaction(p); if (settings::weight_windows_on) { - WeightWindow* ww = search_weight_window(p); + int mesh_bin; + const auto& ww = search_weight_window(p, mesh_bin); if (ww != nullptr) { if (settings::weight_window_checkpoint_collision) - apply_weight_window(p, *ww); + apply_weight_window(p, ww->get_weight_window(p.E(), mesh_bin)); } else if (p.type() == ParticleType::neutron) { apply_russian_roulette(p); } diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 42c09e8c586..c3db4710458 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -49,112 +49,6 @@ openmc::vector> weight_windows_generators; } // namespace variance_reduction -//============================================================================== -// Non-member functions -//============================================================================== - -WeightWindow* search_weight_window(Particle& p) -{ - // TODO: this is a linear search - should do something more clever - for (const auto& ww : variance_reduction::weight_windows) { - int mesh_bin = ww->get_mesh_bin(p); - if (mesh_bin > 0) - return &(ww->get_weight_window(p.E(), mesh_bin)); - } -} - -void apply_weight_windows(Particle& p) -{ - if (!settings::weight_windows_on) - return; - - // WW on photon and neutron only - if (p.type() != ParticleType::neutron && p.type() != ParticleType::photon) - return; - - // skip dead or no energy - if (p.E() <= 0 || !p.alive()) - return; - - WeightWindow* ww = search_weight_window(p); - if (ww != nullptr) - apply_weight_window(p, *ww); -} - -void apply_weight_window(Particle& p, WeightWindow& weight_window) -{ - // If particle has not yet had its birth weight window value set, set it to - // the current weight window (or 1.0 if not born in a weight window). - if (p.wgt_ww_born() == -1.0) { - if (weight_window.is_valid()) { - p.wgt_ww_born() = - (weight_window.lower_weight + weight_window.upper_weight) / 2; - } else { - p.wgt_ww_born() = 1.0; - } - } - - // Normalize weight windows based on particle's starting weight - // and the value of the weight window the particle was born in. - weight_window.scale(p.wgt_born() / p.wgt_ww_born()); - - // get the paramters - double weight = p.wgt(); - - // first check to see if particle should be killed for weight cutoff - if (p.wgt() < weight_window.weight_cutoff) { - p.wgt() = 0.0; - return; - } - - // check if particle is far above current weight window - // only do this if the factor is not already set on the particle and a - // maximum lower bound ratio is specified - if (p.ww_factor() == 0.0 && weight_window.max_lb_ratio > 1.0 && - p.wgt() > weight_window.lower_weight * weight_window.max_lb_ratio) { - p.ww_factor() = - p.wgt() / (weight_window.lower_weight * weight_window.max_lb_ratio); - } - - // move weight window closer to the particle weight if needed - if (p.ww_factor() > 1.0) - weight_window.scale(p.ww_factor()); - - // if particle's weight is above the weight window split until they are within - // the window - if (weight > weight_window.upper_weight) { - // do not further split the particle if above the limit - if (p.n_split() >= settings::max_history_splits) - return; - - double n_split = std::ceil(weight / weight_window.upper_weight); - double max_split = weight_window.max_split; - n_split = std::min(n_split, max_split); - - p.n_split() += n_split; - - // Create secondaries and divide weight among all particles - int i_split = std::round(n_split); - for (int l = 0; l < i_split - 1; l++) { - p.split(weight / n_split); - } - // remaining weight is applied to current particle - p.wgt() = weight / n_split; - - } else if (weight <= weight_window.lower_weight) { - // if the particle weight is below the window, play Russian roulette - double weight_survive = - std::min(weight * weight_window.max_split, weight_window.survival_weight); - russian_roulette(p, weight_survive); - } // else particle is in the window, continue as normal -} - -void free_memory_weight_windows() -{ - variance_reduction::ww_map.clear(); - variance_reduction::weight_windows.clear(); -} - //============================================================================== // WeightWindowSettings implementation //============================================================================== @@ -1009,6 +903,110 @@ void WeightWindowsGenerator::update() const // Non-member functions //============================================================================== +const WeightWindows* search_weight_window(const Particle& p, int& mesh_bin) +{ + // TODO: this is a linear search - should do something more clever + for (const auto& ww : variance_reduction::weight_windows) { + mesh_bin = ww->get_mesh_bin(p); + if (mesh_bin > 0) + return ww.get(); + } + return nullptr; +} + +void apply_weight_windows(Particle& p) +{ + if (!settings::weight_windows_on) + return; + + // WW on photon and neutron only + if (p.type() != ParticleType::neutron && p.type() != ParticleType::photon) + return; + + // skip dead or no energy + if (p.E() <= 0 || !p.alive()) + return; + + int mesh_bin; + const auto& ww = search_weight_window(p, mesh_bin); + if (ww != nullptr) + apply_weight_window(p, ww->get_weight_window(p.E(), mesh_bin)); +} + +void apply_weight_window(Particle& p, WeightWindow weight_window) +{ + // If particle has not yet had its birth weight window value set, set it to + // the current weight window (or 1.0 if not born in a weight window). + if (p.wgt_ww_born() == -1.0) { + if (weight_window.is_valid()) { + p.wgt_ww_born() = + (weight_window.lower_weight + weight_window.upper_weight) / 2; + } else { + p.wgt_ww_born() = 1.0; + } + } + + // Normalize weight windows based on particle's starting weight + // and the value of the weight window the particle was born in. + weight_window.scale(p.wgt_born() / p.wgt_ww_born()); + + // get the paramters + double weight = p.wgt(); + + // first check to see if particle should be killed for weight cutoff + if (p.wgt() < weight_window.weight_cutoff) { + p.wgt() = 0.0; + return; + } + + // check if particle is far above current weight window + // only do this if the factor is not already set on the particle and a + // maximum lower bound ratio is specified + if (p.ww_factor() == 0.0 && weight_window.max_lb_ratio > 1.0 && + p.wgt() > weight_window.lower_weight * weight_window.max_lb_ratio) { + p.ww_factor() = + p.wgt() / (weight_window.lower_weight * weight_window.max_lb_ratio); + } + + // move weight window closer to the particle weight if needed + if (p.ww_factor() > 1.0) + weight_window.scale(p.ww_factor()); + + // if particle's weight is above the weight window split until they are within + // the window + if (weight > weight_window.upper_weight) { + // do not further split the particle if above the limit + if (p.n_split() >= settings::max_history_splits) + return; + + double n_split = std::ceil(weight / weight_window.upper_weight); + double max_split = weight_window.max_split; + n_split = std::min(n_split, max_split); + + p.n_split() += n_split; + + // Create secondaries and divide weight among all particles + int i_split = std::round(n_split); + for (int l = 0; l < i_split - 1; l++) { + p.split(weight / n_split); + } + // remaining weight is applied to current particle + p.wgt() = weight / n_split; + + } else if (weight <= weight_window.lower_weight) { + // if the particle weight is below the window, play Russian roulette + double weight_survive = + std::min(weight * weight_window.max_split, weight_window.survival_weight); + russian_roulette(p, weight_survive); + } // else particle is in the window, continue as normal +} + +void free_memory_weight_windows() +{ + variance_reduction::ww_map.clear(); + variance_reduction::weight_windows.clear(); +} + void finalize_variance_reduction() { for (const auto& wwg : variance_reduction::weight_windows_generators) { From bd66b8ffd094507d25f0165a5fec541f30079ff5 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 01:27:26 +0200 Subject: [PATCH 10/24] fix typo --- src/weight_windows.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index c3db4710458..47fa85f2377 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -294,7 +294,7 @@ WeightWindow WeightWindows::get_weight_window(const Particle& p) const } WeightWindow WeightWindows::get_weight_window( - const double E, const int mesh_bin) const + double E, const int mesh_bin) const { // particle is outside the weight window mesh if (mesh_bin < 0) From 3aaad3fcf7e3986ddd8d2db4fafe9a43fba7be43 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 01:30:53 +0200 Subject: [PATCH 11/24] fix typo --- include/openmc/weight_windows.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index bf49e5ce33f..d9fb1dc735a 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -254,7 +254,7 @@ void free_memory_weight_windows(); //! Search weight window that apply to a particle //! \param[in] p Particle to search weight window for //! \param[out] mesh_bin int mesh bin inside WeightWindows -const WeightWindows* search_weight_window(Particle& p, int& mesh_bin); +const WeightWindows* search_weight_window(const Particle& p, int& mesh_bin); //! Finalize variance reduction objects after all inputs have been read void finalize_variance_reduction(); From f3c63ec80847dc46f3d1c3c29576d88205f51d43 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 06:59:13 +0200 Subject: [PATCH 12/24] apply PR suggestions --- include/openmc/weight_windows.h | 3 +-- src/physics.cpp | 7 +++---- src/physics_mg.cpp | 7 +++---- src/weight_windows.cpp | 12 ++++++------ 4 files changed, 13 insertions(+), 16 deletions(-) diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index d9fb1dc735a..819a7ae821b 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -253,8 +253,7 @@ void free_memory_weight_windows(); //! Search weight window that apply to a particle //! \param[in] p Particle to search weight window for -//! \param[out] mesh_bin int mesh bin inside WeightWindows -const WeightWindows* search_weight_window(const Particle& p, int& mesh_bin); +WeightWindow search_weight_window(const Particle& p); //! Finalize variance reduction objects after all inputs have been read void finalize_variance_reduction(); diff --git a/src/physics.cpp b/src/physics.cpp index 16a1b5706b0..e68f6b96958 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -62,11 +62,10 @@ void collision(Particle& p) } if (settings::weight_windows_on) { - int mesh_bin; - const auto& ww = search_weight_window(p, mesh_bin); - if (ww != nullptr) { + auto ww = search_weight_window(p, mesh_bin); + if (ww.is_valid()) { if (settings::weight_window_checkpoint_collision) - apply_weight_window(p, ww->get_weight_window(p.E(), mesh_bin)); + apply_weight_window(p, ww); } else if (p.type() == ParticleType::neutron) { apply_russian_roulette(p); } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 4549dcffcd6..967e0809868 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -32,11 +32,10 @@ void collision_mg(Particle& p) sample_reaction(p); if (settings::weight_windows_on) { - int mesh_bin; - const auto& ww = search_weight_window(p, mesh_bin); - if (ww != nullptr) { + auto ww = search_weight_window(p, mesh_bin); + if (ww.is_valid()) { if (settings::weight_window_checkpoint_collision) - apply_weight_window(p, ww->get_weight_window(p.E(), mesh_bin)); + apply_weight_window(p, ww); } else if (p.type() == ParticleType::neutron) { apply_russian_roulette(p); } diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 47fa85f2377..86c6e53563f 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -903,15 +903,15 @@ void WeightWindowsGenerator::update() const // Non-member functions //============================================================================== -const WeightWindows* search_weight_window(const Particle& p, int& mesh_bin) +WeightWindow search_weight_window(const Particle& p) { // TODO: this is a linear search - should do something more clever for (const auto& ww : variance_reduction::weight_windows) { mesh_bin = ww->get_mesh_bin(p); if (mesh_bin > 0) - return ww.get(); + return ww->get_weight_window(p.E(), mesh_bin); } - return nullptr; + return {}; } void apply_weight_windows(Particle& p) @@ -928,9 +928,9 @@ void apply_weight_windows(Particle& p) return; int mesh_bin; - const auto& ww = search_weight_window(p, mesh_bin); - if (ww != nullptr) - apply_weight_window(p, ww->get_weight_window(p.E(), mesh_bin)); + auto ww = search_weight_window(p, mesh_bin); + if (ww.is_valid()) + apply_weight_window(p, ww); } void apply_weight_window(Particle& p, WeightWindow weight_window) From d55bfa42008c9eb0e4c70a9bb0c9facd3ae333d2 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 07:02:01 +0200 Subject: [PATCH 13/24] fix typos --- src/physics.cpp | 2 +- src/physics_mg.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index e68f6b96958..57b6aa9f0f1 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -62,7 +62,7 @@ void collision(Particle& p) } if (settings::weight_windows_on) { - auto ww = search_weight_window(p, mesh_bin); + auto ww = search_weight_window(p); if (ww.is_valid()) { if (settings::weight_window_checkpoint_collision) apply_weight_window(p, ww); diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 967e0809868..9b0e423dff0 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -32,7 +32,7 @@ void collision_mg(Particle& p) sample_reaction(p); if (settings::weight_windows_on) { - auto ww = search_weight_window(p, mesh_bin); + auto ww = search_weight_window(p); if (ww.is_valid()) { if (settings::weight_window_checkpoint_collision) apply_weight_window(p, ww); From 1ac07f58db1da6b4c7fc5a6383ca59280f0333b4 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 07:07:42 +0200 Subject: [PATCH 14/24] more fixes --- src/weight_windows.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 86c6e53563f..61f7f4ca154 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -906,6 +906,7 @@ void WeightWindowsGenerator::update() const WeightWindow search_weight_window(const Particle& p) { // TODO: this is a linear search - should do something more clever + int mesh_bin; for (const auto& ww : variance_reduction::weight_windows) { mesh_bin = ww->get_mesh_bin(p); if (mesh_bin > 0) @@ -927,8 +928,7 @@ void apply_weight_windows(Particle& p) if (p.E() <= 0 || !p.alive()) return; - int mesh_bin; - auto ww = search_weight_window(p, mesh_bin); + auto ww = search_weight_window(p); if (ww.is_valid()) apply_weight_window(p, ww); } From 24bf0f1584ed3772255fedadb387d9843fd31eb5 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 17:34:00 +0200 Subject: [PATCH 15/24] fix typo --- src/physics_common.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/physics_common.cpp b/src/physics_common.cpp index a40c6f20ce9..10760ce2da1 100644 --- a/src/physics_common.cpp +++ b/src/physics_common.cpp @@ -21,7 +21,7 @@ void russian_roulette(Particle& p, double weight_survive) void apply_russian_roulette(Particle& p) { // Exit if survival biasing is turned off - if (settings::survival_biasing) + if (!settings::survival_biasing) return; // if survival normalization is on, use normalized weight cutoff and From b17aa1a24dd894d56893384b078e05130916eaaf Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 20:44:15 +0200 Subject: [PATCH 16/24] fix another typo --- src/weight_windows.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 61f7f4ca154..beddf0e49ce 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -909,7 +909,7 @@ WeightWindow search_weight_window(const Particle& p) int mesh_bin; for (const auto& ww : variance_reduction::weight_windows) { mesh_bin = ww->get_mesh_bin(p); - if (mesh_bin > 0) + if (mesh_bin >= 0) return ww->get_weight_window(p.E(), mesh_bin); } return {}; From 2daffd1e76a5549df246371c8204b0994fd7b456 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 20:51:12 +0200 Subject: [PATCH 17/24] another shortcut --- src/weight_windows.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index beddf0e49ce..0dd33306e2d 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -935,6 +935,10 @@ void apply_weight_windows(Particle& p) void apply_weight_window(Particle& p, WeightWindow weight_window) { + // skip dead or no energy + if (p.E() <= 0 || !p.alive()) + return; + // If particle has not yet had its birth weight window value set, set it to // the current weight window (or 1.0 if not born in a weight window). if (p.wgt_ww_born() == -1.0) { From 5fe66d08a31dbf155baa702bdaf2e49fb0430a1d Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 21:51:36 +0200 Subject: [PATCH 18/24] improve code structure --- src/physics.cpp | 1961 ++++++++++++++++++++++---------------------- src/physics_mg.cpp | 395 ++++----- 2 files changed, 1192 insertions(+), 1164 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index 57b6aa9f0f1..ee94d068c78 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -61,1165 +61,1182 @@ void collision(Particle& p) break; } - if (settings::weight_windows_on) { + // If collision checkpoint is enabled, apply weight window + // if valid and apply global russian roulette if not. + if (settings::weight_window_checkpoint_collision) { auto ww = search_weight_window(p); if (ww.is_valid()) { - if (settings::weight_window_checkpoint_collision) - apply_weight_window(p, ww); - } else if (p.type() == ParticleType::neutron) { - apply_russian_roulette(p); + apply_weight_window(p, ww); + else if (p.type() == ParticleType::neutron) + { + apply_russian_roulette(p) + } + // If collision checkpoint is disabled, apply russian roulette if + // the particle is a neutron and it is outside the weight window domain + } else if (settings::weight_windows_on && + p.type() == ParticleType::neutron) { + auto ww = search_weight_window(p); + if (!ww.is_valid()) + apply_russian_roulette(p); } - } - - // Kill particle if energy falls below cutoff - int type = static_cast(p.type()); - if (p.E() < settings::energy_cutoff[type]) { - p.wgt() = 0.0; - } - // Display information about collision - if (settings::verbosity >= 10 || p.trace()) { - std::string msg; - if (p.event() == TallyEvent::KILL) { - msg = fmt::format(" Killed. Energy = {} eV.", p.E()); - } else if (p.type() == ParticleType::neutron) { - msg = fmt::format(" {} with {}. Energy = {} eV.", - reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_, - p.E()); - } else if (p.type() == ParticleType::photon) { - msg = fmt::format(" {} with {}. Energy = {} eV.", - reaction_name(p.event_mt()), - to_element(data::nuclides[p.event_nuclide()]->name_), p.E()); - } else { - msg = fmt::format(" Disappeared. Energy = {} eV.", p.E()); + // Kill particle if energy falls below cutoff + int type = static_cast(p.type()); + if (p.E() < settings::energy_cutoff[type]) { + p.wgt() = 0.0; } - write_message(msg, 1); - } -} -void sample_neutron_reaction(Particle& p) -{ - // Sample a nuclide within the material - int i_nuclide = sample_nuclide(p); - - // Save which nuclide particle had collision with - p.event_nuclide() = i_nuclide; - - // Create fission bank sites. Note that while a fission reaction is sampled, - // it never actually "happens", i.e. the weight of the particle does not - // change when sampling fission sites. The following block handles all - // absorption (including fission) - - const auto& nuc {data::nuclides[i_nuclide]}; - - if (nuc->fissionable_ && p.neutron_xs(i_nuclide).fission > 0.0) { - auto& rx = sample_fission(i_nuclide, p); - if (settings::run_mode == RunMode::EIGENVALUE) { - create_fission_sites(p, i_nuclide, rx); - } else if (settings::run_mode == RunMode::FIXED_SOURCE && - settings::create_fission_neutrons) { - create_fission_sites(p, i_nuclide, rx); - - // Make sure particle population doesn't grow out of control for - // subcritical multiplication problems. - if (p.secondary_bank().size() >= settings::max_secondaries) { - fatal_error( - "The secondary particle bank appears to be growing without " - "bound. You are likely running a subcritical multiplication problem " - "with k-effective close to or greater than one."); + // Display information about collision + if (settings::verbosity >= 10 || p.trace()) { + std::string msg; + if (p.event() == TallyEvent::KILL) { + msg = fmt::format(" Killed. Energy = {} eV.", p.E()); + } else if (p.type() == ParticleType::neutron) { + msg = fmt::format(" {} with {}. Energy = {} eV.", + reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_, + p.E()); + } else if (p.type() == ParticleType::photon) { + msg = fmt::format(" {} with {}. Energy = {} eV.", + reaction_name(p.event_mt()), + to_element(data::nuclides[p.event_nuclide()]->name_), p.E()); + } else { + msg = fmt::format(" Disappeared. Energy = {} eV.", p.E()); } + write_message(msg, 1); } - p.event_mt() = rx.mt_; } - // Create secondary photons - if (settings::photon_transport) { - sample_secondary_photons(p, i_nuclide); - } + void sample_neutron_reaction(Particle & p) + { + // Sample a nuclide within the material + int i_nuclide = sample_nuclide(p); - // If survival biasing is being used, the following subroutine adjusts the - // weight of the particle. Otherwise, it checks to see if absorption occurs + // Save which nuclide particle had collision with + p.event_nuclide() = i_nuclide; - if (p.neutron_xs(i_nuclide).absorption > 0.0) { - absorption(p, i_nuclide); - } - if (!p.alive()) - return; - - // Sample a scattering reaction and determine the secondary energy of the - // exiting neutron - const auto& ncrystal_mat = model::materials[p.material()]->ncrystal_mat(); - if (ncrystal_mat && p.E() < NCRYSTAL_MAX_ENERGY) { - ncrystal_mat.scatter(p); - } else { - scatter(p, i_nuclide); - } - - // Advance URR seed stream 'N' times after energy changes - if (p.E() != p.E_last()) { - advance_prn_seed(data::nuclides.size(), &p.seeds(STREAM_URR_PTABLE)); - } + // Create fission bank sites. Note that while a fission reaction is sampled, + // it never actually "happens", i.e. the weight of the particle does not + // change when sampling fission sites. The following block handles all + // absorption (including fission) - // Play russian roulette if there are no weight windows - if (!settings::weight_windows_on) - apply_russian_roulette(p); -} + const auto& nuc {data::nuclides[i_nuclide]}; -void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) -{ - // If uniform fission source weighting is turned on, we increase or decrease - // the expected number of fission sites produced - double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; - - // Determine the expected number of neutrons produced - double nu_t = p.wgt() / simulation::keff * weight * - p.neutron_xs(i_nuclide).nu_fission / - p.neutron_xs(i_nuclide).total; - - // Sample the number of neutrons produced - int nu = static_cast(nu_t); - if (prn(p.current_seed()) <= (nu_t - nu)) - ++nu; - - // If no neutrons were produced then don't continue - if (nu == 0) - return; - - // Initialize the counter of delayed neutrons encountered for each delayed - // group. - double nu_d[MAX_DELAYED_GROUPS] = {0.}; - - // Clear out particle's nu fission bank - p.nu_bank().clear(); - - p.fission() = true; - - // Determine whether to place fission sites into the shared fission bank - // or the secondary particle bank. - bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); - - // Counter for the number of fission sites successfully stored to the shared - // fission bank or the secondary particle bank - int n_sites_stored; - - for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { - // Initialize fission site object with particle data - SourceSite site; - site.r = p.r(); - site.particle = ParticleType::neutron; - site.time = p.time(); - site.wgt = 1. / weight; - site.surf_id = 0; - - // Sample delayed group and angle/energy for fission reaction - sample_fission_neutron(i_nuclide, rx, &site, p); - - // Reject site if it exceeds time cutoff - if (site.delayed_group > 0) { - double t_cutoff = settings::time_cutoff[static_cast(site.particle)]; - if (site.time > t_cutoff) { - continue; + if (nuc->fissionable_ && p.neutron_xs(i_nuclide).fission > 0.0) { + auto& rx = sample_fission(i_nuclide, p); + if (settings::run_mode == RunMode::EIGENVALUE) { + create_fission_sites(p, i_nuclide, rx); + } else if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::create_fission_neutrons) { + create_fission_sites(p, i_nuclide, rx); + + // Make sure particle population doesn't grow out of control for + // subcritical multiplication problems. + if (p.secondary_bank().size() >= settings::max_secondaries) { + fatal_error( + "The secondary particle bank appears to be growing without " + "bound. You are likely running a subcritical multiplication " + "problem " + "with k-effective close to or greater than one."); + } } + p.event_mt() = rx.mt_; } - // Set parent and progeny IDs - site.parent_id = p.id(); - site.progeny_id = p.n_progeny()++; - - // Store fission site in bank - if (use_fission_bank) { - int64_t idx = simulation::fission_bank.thread_safe_append(site); - if (idx == -1) { - warning( - "The shared fission bank is full. Additional fission sites created " - "in this generation will not be banked. Results may be " - "non-deterministic."); - - // Decrement number of particle progeny as storage was unsuccessful. - // This step is needed so that the sum of all progeny is equal to the - // size of the shared fission bank. - p.n_progeny()--; - - // 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, idx); - } + // Create secondary photons + if (settings::photon_transport) { + sample_secondary_photons(p, i_nuclide); + } + + // If survival biasing is being used, the following subroutine adjusts the + // weight of the particle. Otherwise, it checks to see if absorption occurs + + if (p.neutron_xs(i_nuclide).absorption > 0.0) { + absorption(p, i_nuclide); + } + if (!p.alive()) + return; + + // Sample a scattering reaction and determine the secondary energy of the + // exiting neutron + const auto& ncrystal_mat = model::materials[p.material()]->ncrystal_mat(); + if (ncrystal_mat && p.E() < NCRYSTAL_MAX_ENERGY) { + ncrystal_mat.scatter(p); } else { - p.secondary_bank().push_back(site); + scatter(p, i_nuclide); } - // Increment the number of neutrons born delayed - if (site.delayed_group > 0) { - nu_d[site.delayed_group - 1]++; + // Advance URR seed stream 'N' times after energy changes + if (p.E() != p.E_last()) { + advance_prn_seed(data::nuclides.size(), &p.seeds(STREAM_URR_PTABLE)); } - // Write fission particles to nuBank - NuBank& nu_bank_entry = p.nu_bank().emplace_back(); - nu_bank_entry.wgt = site.wgt; - nu_bank_entry.E = site.E; - nu_bank_entry.delayed_group = site.delayed_group; + // Play russian roulette if there are no weight windows + if (!settings::weight_windows_on) + apply_russian_roulette(p); } - // If shared fission bank was full, and no fissions could be added, - // set the particle fission flag to false. - if (n_sites_stored == 0) { - p.fission() = false; - return; - } + void create_fission_sites(Particle & p, int i_nuclide, const Reaction& rx) + { + // If uniform fission source weighting is turned on, we increase or decrease + // the expected number of fission sites produced + double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; - // Set nu to the number of fission sites successfully stored. If the fission - // bank was not found to be full then these values are already equivalent. - nu = n_sites_stored; + // Determine the expected number of neutrons produced + double nu_t = p.wgt() / simulation::keff * weight * + p.neutron_xs(i_nuclide).nu_fission / + p.neutron_xs(i_nuclide).total; - // Store the total weight banked for analog fission tallies - p.n_bank() = nu; - p.wgt_bank() = nu / weight; - for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { - p.n_delayed_bank(d) = nu_d[d]; - } -} + // Sample the number of neutrons produced + int nu = static_cast(nu_t); + if (prn(p.current_seed()) <= (nu_t - nu)) + ++nu; -void sample_photon_reaction(Particle& p) -{ - // Kill photon if below energy cutoff -- an extra check is made here because - // photons with energy below the cutoff may have been produced by neutrons - // reactions or atomic relaxation - int photon = static_cast(ParticleType::photon); - if (p.E() < settings::energy_cutoff[photon]) { - p.E() = 0.0; - p.wgt() = 0.0; - return; - } + // If no neutrons were produced then don't continue + if (nu == 0) + return; - // Sample element within material - int i_element = sample_element(p); - const auto& micro {p.photon_xs(i_element)}; - const auto& element {*data::elements[i_element]}; + // Initialize the counter of delayed neutrons encountered for each delayed + // group. + double nu_d[MAX_DELAYED_GROUPS] = {0.}; + + // Clear out particle's nu fission bank + p.nu_bank().clear(); + + p.fission() = true; + + // Determine whether to place fission sites into the shared fission bank + // or the secondary particle bank. + bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); + + // Counter for the number of fission sites successfully stored to the shared + // fission bank or the secondary particle bank + int n_sites_stored; + + for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { + // Initialize fission site object with particle data + SourceSite site; + site.r = p.r(); + site.particle = ParticleType::neutron; + site.time = p.time(); + site.wgt = 1. / weight; + site.surf_id = 0; + + // Sample delayed group and angle/energy for fission reaction + sample_fission_neutron(i_nuclide, rx, &site, p); + + // Reject site if it exceeds time cutoff + if (site.delayed_group > 0) { + double t_cutoff = + settings::time_cutoff[static_cast(site.particle)]; + if (site.time > t_cutoff) { + continue; + } + } - // Calculate photon energy over electron rest mass equivalent - double alpha = p.E() / MASS_ELECTRON_EV; + // Set parent and progeny IDs + site.parent_id = p.id(); + site.progeny_id = p.n_progeny()++; + + // Store fission site in bank + if (use_fission_bank) { + int64_t idx = simulation::fission_bank.thread_safe_append(site); + if (idx == -1) { + warning( + "The shared fission bank is full. Additional fission sites created " + "in this generation will not be banked. Results may be " + "non-deterministic."); + + // Decrement number of particle progeny as storage was unsuccessful. + // This step is needed so that the sum of all progeny is equal to the + // size of the shared fission bank. + p.n_progeny()--; + + // 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, idx); + } + } else { + p.secondary_bank().push_back(site); + } - // For tallying purposes, this routine might be called directly. In that - // case, we need to sample a reaction via the cutoff variable - double prob = 0.0; - double cutoff = prn(p.current_seed()) * micro.total; + // Increment the number of neutrons born delayed + if (site.delayed_group > 0) { + nu_d[site.delayed_group - 1]++; + } - // Coherent (Rayleigh) scattering - prob += micro.coherent; - if (prob > cutoff) { - p.mu() = element.rayleigh_scatter(alpha, p.current_seed()); - p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); - p.event() = TallyEvent::SCATTER; - p.event_mt() = COHERENT; - return; - } + // Write fission particles to nuBank + NuBank& nu_bank_entry = p.nu_bank().emplace_back(); + nu_bank_entry.wgt = site.wgt; + nu_bank_entry.E = site.E; + nu_bank_entry.delayed_group = site.delayed_group; + } - // Incoherent (Compton) scattering - prob += micro.incoherent; - if (prob > cutoff) { - double alpha_out; - int i_shell; - element.compton_scatter( - alpha, true, &alpha_out, &p.mu(), &i_shell, p.current_seed()); - - // Determine binding energy of shell. The binding energy is 0.0 if - // doppler broadening is not used. - double e_b; - if (i_shell == -1) { - e_b = 0.0; - } else { - e_b = element.binding_energy_[i_shell]; + // If shared fission bank was full, and no fissions could be added, + // set the particle fission flag to false. + if (n_sites_stored == 0) { + p.fission() = false; + return; } - // Create Compton electron - double phi = uniform_distribution(0., 2.0 * PI, p.current_seed()); - double E_electron = (alpha - alpha_out) * MASS_ELECTRON_EV - e_b; - int electron = static_cast(ParticleType::electron); - if (E_electron >= settings::energy_cutoff[electron]) { - double mu_electron = (alpha - alpha_out * p.mu()) / - std::sqrt(alpha * alpha + alpha_out * alpha_out - - 2.0 * alpha * alpha_out * p.mu()); - Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed()); - p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); + // Set nu to the number of fission sites successfully stored. If the fission + // bank was not found to be full then these values are already equivalent. + nu = n_sites_stored; + + // Store the total weight banked for analog fission tallies + p.n_bank() = nu; + p.wgt_bank() = nu / weight; + for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { + p.n_delayed_bank(d) = nu_d[d]; } + } - // Allow electrons to fill orbital and produce Auger electrons and - // fluorescent photons. Since Compton subshell data does not match atomic - // relaxation data, use the mapping between the data to find the subshell - if (i_shell >= 0 && element.subshell_map_[i_shell] >= 0) { - element.atomic_relaxation(element.subshell_map_[i_shell], p); + void sample_photon_reaction(Particle & p) + { + // Kill photon if below energy cutoff -- an extra check is made here because + // photons with energy below the cutoff may have been produced by neutrons + // reactions or atomic relaxation + int photon = static_cast(ParticleType::photon); + if (p.E() < settings::energy_cutoff[photon]) { + p.E() = 0.0; + p.wgt() = 0.0; + return; } - phi += PI; - p.E() = alpha_out * MASS_ELECTRON_EV; - p.u() = rotate_angle(p.u(), p.mu(), &phi, p.current_seed()); - p.event() = TallyEvent::SCATTER; - p.event_mt() = INCOHERENT; - return; - } + // Sample element within material + int i_element = sample_element(p); + const auto& micro {p.photon_xs(i_element)}; + const auto& element {*data::elements[i_element]}; - // Photoelectric effect - double prob_after = prob + micro.photoelectric; + // Calculate photon energy over electron rest mass equivalent + double alpha = p.E() / MASS_ELECTRON_EV; - if (prob_after > cutoff) { - // Get grid index, interpolation factor, and bounding subshell - // cross sections - int i_grid = micro.index_grid; - double f = micro.interp_factor; - const auto& xs_lower = xt::row(element.cross_sections_, i_grid); - const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1); + // For tallying purposes, this routine might be called directly. In that + // case, we need to sample a reaction via the cutoff variable + double prob = 0.0; + double cutoff = prn(p.current_seed()) * micro.total; - for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) { - const auto& shell {element.shells_[i_shell]}; + // Coherent (Rayleigh) scattering + prob += micro.coherent; + if (prob > cutoff) { + p.mu() = element.rayleigh_scatter(alpha, p.current_seed()); + p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); + p.event() = TallyEvent::SCATTER; + p.event_mt() = COHERENT; + return; + } - // Check threshold of reaction - if (xs_lower(i_shell) == 0) - continue; + // Incoherent (Compton) scattering + prob += micro.incoherent; + if (prob > cutoff) { + double alpha_out; + int i_shell; + element.compton_scatter( + alpha, true, &alpha_out, &p.mu(), &i_shell, p.current_seed()); + + // Determine binding energy of shell. The binding energy is 0.0 if + // doppler broadening is not used. + double e_b; + if (i_shell == -1) { + e_b = 0.0; + } else { + e_b = element.binding_energy_[i_shell]; + } - // Evaluation subshell photoionization cross section - prob += std::exp( - xs_lower(i_shell) + f * (xs_upper(i_shell) - xs_lower(i_shell))); + // Create Compton electron + double phi = uniform_distribution(0., 2.0 * PI, p.current_seed()); + double E_electron = (alpha - alpha_out) * MASS_ELECTRON_EV - e_b; + int electron = static_cast(ParticleType::electron); + if (E_electron >= settings::energy_cutoff[electron]) { + double mu_electron = (alpha - alpha_out * p.mu()) / + std::sqrt(alpha * alpha + alpha_out * alpha_out - + 2.0 * alpha * alpha_out * p.mu()); + Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed()); + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); + } - if (prob > cutoff) { - // Determine binding energy based on whether atomic relaxation data is - // present (if not, use value from Compton profile data) - double binding_energy = element.has_atomic_relaxation_ - ? shell.binding_energy - : element.binding_energy_[i_shell]; - - // Determine energy of secondary electron - double E_electron = p.E() - binding_energy; - - // Sample mu using non-relativistic Sauter distribution. - // See Eqns 3.19 and 3.20 in "Implementing a photon physics - // model in Serpent 2" by Toni Kaltiaisenaho - double mu; - while (true) { - double r = prn(p.current_seed()); - if (4.0 * (1.0 - r) * r >= prn(p.current_seed())) { - double rel_vel = - std::sqrt(E_electron * (E_electron + 2.0 * MASS_ELECTRON_EV)) / - (E_electron + MASS_ELECTRON_EV); - mu = - (2.0 * r + rel_vel - 1.0) / (2.0 * rel_vel * r - rel_vel + 1.0); - break; + // Allow electrons to fill orbital and produce Auger electrons and + // fluorescent photons. Since Compton subshell data does not match atomic + // relaxation data, use the mapping between the data to find the subshell + if (i_shell >= 0 && element.subshell_map_[i_shell] >= 0) { + element.atomic_relaxation(element.subshell_map_[i_shell], p); + } + + phi += PI; + p.E() = alpha_out * MASS_ELECTRON_EV; + p.u() = rotate_angle(p.u(), p.mu(), &phi, p.current_seed()); + p.event() = TallyEvent::SCATTER; + p.event_mt() = INCOHERENT; + return; + } + + // Photoelectric effect + double prob_after = prob + micro.photoelectric; + + if (prob_after > cutoff) { + // Get grid index, interpolation factor, and bounding subshell + // cross sections + int i_grid = micro.index_grid; + double f = micro.interp_factor; + const auto& xs_lower = xt::row(element.cross_sections_, i_grid); + const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1); + + for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) { + const auto& shell {element.shells_[i_shell]}; + + // Check threshold of reaction + if (xs_lower(i_shell) == 0) + continue; + + // Evaluation subshell photoionization cross section + prob += std::exp( + xs_lower(i_shell) + f * (xs_upper(i_shell) - xs_lower(i_shell))); + + if (prob > cutoff) { + // Determine binding energy based on whether atomic relaxation data is + // present (if not, use value from Compton profile data) + double binding_energy = element.has_atomic_relaxation_ + ? shell.binding_energy + : element.binding_energy_[i_shell]; + + // Determine energy of secondary electron + double E_electron = p.E() - binding_energy; + + // Sample mu using non-relativistic Sauter distribution. + // See Eqns 3.19 and 3.20 in "Implementing a photon physics + // model in Serpent 2" by Toni Kaltiaisenaho + double mu; + while (true) { + double r = prn(p.current_seed()); + if (4.0 * (1.0 - r) * r >= prn(p.current_seed())) { + double rel_vel = + std::sqrt(E_electron * (E_electron + 2.0 * MASS_ELECTRON_EV)) / + (E_electron + MASS_ELECTRON_EV); + mu = + (2.0 * r + rel_vel - 1.0) / (2.0 * rel_vel * r - rel_vel + 1.0); + break; + } } + + double phi = uniform_distribution(0., 2.0 * PI, p.current_seed()); + Direction u; + u.x = mu; + u.y = std::sqrt(1.0 - mu * mu) * std::cos(phi); + u.z = std::sqrt(1.0 - mu * mu) * std::sin(phi); + + // Create secondary electron + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); + + // Allow electrons to fill orbital and produce auger electrons + // and fluorescent photons + element.atomic_relaxation(i_shell, p); + p.event() = TallyEvent::ABSORB; + p.event_mt() = 533 + shell.index_subshell; + p.wgt() = 0.0; + p.E() = 0.0; + return; } + } + } + prob = prob_after; - double phi = uniform_distribution(0., 2.0 * PI, p.current_seed()); - Direction u; - u.x = mu; - u.y = std::sqrt(1.0 - mu * mu) * std::cos(phi); - u.z = std::sqrt(1.0 - mu * mu) * std::sin(phi); + // Pair production + prob += micro.pair_production; + if (prob > cutoff) { + double E_electron, E_positron; + double mu_electron, mu_positron; + element.pair_production(alpha, &E_electron, &E_positron, &mu_electron, + &mu_positron, p.current_seed()); - // Create secondary electron - p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); + // Create secondary electron + Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed()); + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); - // Allow electrons to fill orbital and produce auger electrons - // and fluorescent photons - element.atomic_relaxation(i_shell, p); - p.event() = TallyEvent::ABSORB; - p.event_mt() = 533 + shell.index_subshell; - p.wgt() = 0.0; - p.E() = 0.0; - return; - } + // Create secondary positron + u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed()); + p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron); + + p.event() = TallyEvent::ABSORB; + p.event_mt() = PAIR_PROD; + p.wgt() = 0.0; + p.E() = 0.0; } } - prob = prob_after; - - // Pair production - prob += micro.pair_production; - if (prob > cutoff) { - double E_electron, E_positron; - double mu_electron, mu_positron; - element.pair_production(alpha, &E_electron, &E_positron, &mu_electron, - &mu_positron, p.current_seed()); - // Create secondary electron - Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed()); - p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); + void sample_electron_reaction(Particle & p) + { + // TODO: create reaction types - // Create secondary positron - u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed()); - p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron); + if (settings::electron_treatment == ElectronTreatment::TTB) { + double E_lost; + thick_target_bremsstrahlung(p, &E_lost); + } - p.event() = TallyEvent::ABSORB; - p.event_mt() = PAIR_PROD; - p.wgt() = 0.0; p.E() = 0.0; + p.wgt() = 0.0; + p.event() = TallyEvent::ABSORB; } -} -void sample_electron_reaction(Particle& p) -{ - // TODO: create reaction types + void sample_positron_reaction(Particle & p) + { + // TODO: create reaction types - if (settings::electron_treatment == ElectronTreatment::TTB) { - double E_lost; - thick_target_bremsstrahlung(p, &E_lost); - } + if (settings::electron_treatment == ElectronTreatment::TTB) { + double E_lost; + thick_target_bremsstrahlung(p, &E_lost); + } - p.E() = 0.0; - p.wgt() = 0.0; - p.event() = TallyEvent::ABSORB; -} + // Sample angle isotropically + Direction u = isotropic_direction(p.current_seed()); -void sample_positron_reaction(Particle& p) -{ - // TODO: create reaction types + // Create annihilation photon pair traveling in opposite directions + p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon); + p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon); - if (settings::electron_treatment == ElectronTreatment::TTB) { - double E_lost; - thick_target_bremsstrahlung(p, &E_lost); + p.E() = 0.0; + p.wgt() = 0.0; + p.event() = TallyEvent::ABSORB; } - // Sample angle isotropically - Direction u = isotropic_direction(p.current_seed()); + int sample_nuclide(Particle & p) + { + // Sample cumulative distribution function + double cutoff = prn(p.current_seed()) * p.macro_xs().total; - // Create annihilation photon pair traveling in opposite directions - p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon); - p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon); + // Get pointers to nuclide/density arrays + const auto& mat {model::materials[p.material()]}; + int n = mat->nuclide_.size(); - p.E() = 0.0; - p.wgt() = 0.0; - p.event() = TallyEvent::ABSORB; -} + double prob = 0.0; + for (int i = 0; i < n; ++i) { + // Get atom density + int i_nuclide = mat->nuclide_[i]; + double atom_density = mat->atom_density(i, p.density_mult()); + + // Increment probability to compare to cutoff + prob += atom_density * p.neutron_xs(i_nuclide).total; + if (prob >= cutoff) + return i_nuclide; + } -int sample_nuclide(Particle& p) -{ - // Sample cumulative distribution function - double cutoff = prn(p.current_seed()) * p.macro_xs().total; - - // Get pointers to nuclide/density arrays - const auto& mat {model::materials[p.material()]}; - int n = mat->nuclide_.size(); - - double prob = 0.0; - for (int i = 0; i < n; ++i) { - // Get atom density - int i_nuclide = mat->nuclide_[i]; - double atom_density = mat->atom_density(i, p.density_mult()); - - // Increment probability to compare to cutoff - prob += atom_density * p.neutron_xs(i_nuclide).total; - if (prob >= cutoff) - return i_nuclide; + // If we reach here, no nuclide was sampled + p.write_restart(); + throw std::runtime_error {"Did not sample any nuclide during collision."}; } - // If we reach here, no nuclide was sampled - p.write_restart(); - throw std::runtime_error {"Did not sample any nuclide during collision."}; -} + int sample_element(Particle & p) + { + // Sample cumulative distribution function + double cutoff = prn(p.current_seed()) * p.macro_xs().total; -int sample_element(Particle& p) -{ - // Sample cumulative distribution function - double cutoff = prn(p.current_seed()) * p.macro_xs().total; + // Get pointers to elements, densities + const auto& mat {model::materials[p.material()]}; - // Get pointers to elements, densities - const auto& mat {model::materials[p.material()]}; - - double prob = 0.0; - for (int i = 0; i < mat->element_.size(); ++i) { - // Find atom density - int i_element = mat->element_[i]; - double atom_density = mat->atom_density(i, p.density_mult()); + double prob = 0.0; + for (int i = 0; i < mat->element_.size(); ++i) { + // Find atom density + int i_element = mat->element_[i]; + double atom_density = mat->atom_density(i, p.density_mult()); - // Determine microscopic cross section - double sigma = atom_density * p.photon_xs(i_element).total; + // Determine microscopic cross section + double sigma = atom_density * p.photon_xs(i_element).total; - // Increment probability to compare to cutoff - prob += sigma; - if (prob > cutoff) { - // Save which nuclide particle had collision with for tally purpose - p.event_nuclide() = mat->nuclide_[i]; + // Increment probability to compare to cutoff + prob += sigma; + if (prob > cutoff) { + // Save which nuclide particle had collision with for tally purpose + p.event_nuclide() = mat->nuclide_[i]; - return i_element; + return i_element; + } } - } - - // If we made it here, no element was sampled - p.write_restart(); - fatal_error("Did not sample any element during collision."); -} -Reaction& sample_fission(int i_nuclide, Particle& p) -{ - // Get pointer to nuclide - const auto& nuc {data::nuclides[i_nuclide]}; - - // If we're in the URR, by default use the first fission reaction. We also - // default to the first reaction if we know that there are no partial fission - // reactions - if (p.neutron_xs(i_nuclide).use_ptable || !nuc->has_partial_fission_) { - return *nuc->fission_rx_[0]; + // If we made it here, no element was sampled + p.write_restart(); + fatal_error("Did not sample any element during collision."); } - // Check to see if we are in a windowed multipole range. WMP only supports - // the first fission reaction. - if (nuc->multipole_) { - if (p.E() >= nuc->multipole_->E_min_ && p.E() <= nuc->multipole_->E_max_) { + Reaction& sample_fission(int i_nuclide, Particle& p) + { + // Get pointer to nuclide + const auto& nuc {data::nuclides[i_nuclide]}; + + // If we're in the URR, by default use the first fission reaction. We also + // default to the first reaction if we know that there are no partial + // fission reactions + if (p.neutron_xs(i_nuclide).use_ptable || !nuc->has_partial_fission_) { return *nuc->fission_rx_[0]; } - } - // Get grid index and interpolation factor and sample fission cdf - const auto& micro = p.neutron_xs(i_nuclide); - double cutoff = prn(p.current_seed()) * p.neutron_xs(i_nuclide).fission; - double prob = 0.0; + // Check to see if we are in a windowed multipole range. WMP only supports + // the first fission reaction. + if (nuc->multipole_) { + if (p.E() >= nuc->multipole_->E_min_ && + p.E() <= nuc->multipole_->E_max_) { + return *nuc->fission_rx_[0]; + } + } - // Loop through each partial fission reaction type - for (auto& rx : nuc->fission_rx_) { - // add to cumulative probability - prob += rx->xs(micro); + // Get grid index and interpolation factor and sample fission cdf + const auto& micro = p.neutron_xs(i_nuclide); + double cutoff = prn(p.current_seed()) * p.neutron_xs(i_nuclide).fission; + double prob = 0.0; - // Create fission bank sites if fission occurs - if (prob > cutoff) - return *rx; + // Loop through each partial fission reaction type + for (auto& rx : nuc->fission_rx_) { + // add to cumulative probability + prob += rx->xs(micro); + + // Create fission bank sites if fission occurs + if (prob > cutoff) + return *rx; + } + + // If we reached here, no reaction was sampled + throw std::runtime_error { + "No fission reaction was sampled for " + nuc->name_}; } - // If we reached here, no reaction was sampled - throw std::runtime_error { - "No fission reaction was sampled for " + nuc->name_}; -} + void sample_photon_product( + int i_nuclide, Particle& p, int* i_rx, int* i_product) + { + // Get grid index and interpolation factor and sample photon production cdf + const auto& micro = p.neutron_xs(i_nuclide); + double cutoff = prn(p.current_seed()) * micro.photon_prod; + double prob = 0.0; + + // Loop through each reaction type + const auto& nuc {data::nuclides[i_nuclide]}; + for (int i = 0; i < nuc->reactions_.size(); ++i) { + // Evaluate neutron cross section + const auto& rx = nuc->reactions_[i]; + double xs = rx->xs(micro); -void sample_photon_product( - int i_nuclide, Particle& p, int* i_rx, int* i_product) -{ - // Get grid index and interpolation factor and sample photon production cdf - const auto& micro = p.neutron_xs(i_nuclide); - double cutoff = prn(p.current_seed()) * micro.photon_prod; - double prob = 0.0; - - // Loop through each reaction type - const auto& nuc {data::nuclides[i_nuclide]}; - for (int i = 0; i < nuc->reactions_.size(); ++i) { - // Evaluate neutron cross section - const auto& rx = nuc->reactions_[i]; - double xs = rx->xs(micro); - - // if cross section is zero for this reaction, skip it - if (xs == 0.0) - continue; - - for (int j = 0; j < rx->products_.size(); ++j) { - if (rx->products_[j].particle_ == ParticleType::photon) { - // For fission, artificially increase the photon yield to account - // for delayed photons - double f = 1.0; - if (settings::delayed_photon_scaling) { - if (is_fission(rx->mt_)) { - if (nuc->prompt_photons_ && nuc->delayed_photons_) { - double energy_prompt = (*nuc->prompt_photons_)(p.E()); - double energy_delayed = (*nuc->delayed_photons_)(p.E()); - f = (energy_prompt + energy_delayed) / (energy_prompt); + // if cross section is zero for this reaction, skip it + if (xs == 0.0) + continue; + + for (int j = 0; j < rx->products_.size(); ++j) { + if (rx->products_[j].particle_ == ParticleType::photon) { + // For fission, artificially increase the photon yield to account + // for delayed photons + double f = 1.0; + if (settings::delayed_photon_scaling) { + if (is_fission(rx->mt_)) { + if (nuc->prompt_photons_ && nuc->delayed_photons_) { + double energy_prompt = (*nuc->prompt_photons_)(p.E()); + double energy_delayed = (*nuc->delayed_photons_)(p.E()); + f = (energy_prompt + energy_delayed) / (energy_prompt); + } } } - } - // add to cumulative probability - prob += f * (*rx->products_[j].yield_)(p.E()) * xs; + // add to cumulative probability + prob += f * (*rx->products_[j].yield_)(p.E()) * xs; - *i_rx = i; - *i_product = j; - if (prob > cutoff) - return; + *i_rx = i; + *i_product = j; + if (prob > cutoff) + return; + } } } } -} -void absorption(Particle& p, int i_nuclide) -{ - if (settings::survival_biasing) { - // Determine weight absorbed in survival biasing - const double wgt_absorb = p.wgt() * p.neutron_xs(i_nuclide).absorption / - p.neutron_xs(i_nuclide).total; - - // Adjust weight of particle by probability of absorption - p.wgt() -= wgt_absorb; - - // Score implicit absorption estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE) { - p.keff_tally_absorption() += wgt_absorb * - p.neutron_xs(i_nuclide).nu_fission / - p.neutron_xs(i_nuclide).absorption; - } - } else { - // See if disappearance reaction happens - if (p.neutron_xs(i_nuclide).absorption > - prn(p.current_seed()) * p.neutron_xs(i_nuclide).total) { - // Score absorption estimate of keff + void absorption(Particle & p, int i_nuclide) + { + if (settings::survival_biasing) { + // Determine weight absorbed in survival biasing + const double wgt_absorb = p.wgt() * p.neutron_xs(i_nuclide).absorption / + p.neutron_xs(i_nuclide).total; + + // Adjust weight of particle by probability of absorption + p.wgt() -= wgt_absorb; + + // Score implicit absorption estimate of keff if (settings::run_mode == RunMode::EIGENVALUE) { - p.keff_tally_absorption() += p.wgt() * + p.keff_tally_absorption() += wgt_absorb * p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).absorption; } + } else { + // See if disappearance reaction happens + if (p.neutron_xs(i_nuclide).absorption > + prn(p.current_seed()) * p.neutron_xs(i_nuclide).total) { + // Score absorption estimate of keff + if (settings::run_mode == RunMode::EIGENVALUE) { + p.keff_tally_absorption() += p.wgt() * + p.neutron_xs(i_nuclide).nu_fission / + p.neutron_xs(i_nuclide).absorption; + } - p.wgt() = 0.0; - p.event() = TallyEvent::ABSORB; - if (!p.fission()) { - p.event_mt() = N_DISAPPEAR; + p.wgt() = 0.0; + p.event() = TallyEvent::ABSORB; + if (!p.fission()) { + p.event_mt() = N_DISAPPEAR; + } } } } -} -void scatter(Particle& p, int i_nuclide) -{ - // copy incoming direction - Direction u_old {p.u()}; - - // Get pointer to nuclide and grid index/interpolation factor - const auto& nuc {data::nuclides[i_nuclide]}; - const auto& micro {p.neutron_xs(i_nuclide)}; - int i_temp = micro.index_temp; - - // For tallying purposes, this routine might be called directly. In that - // case, we need to sample a reaction via the cutoff variable - double cutoff = prn(p.current_seed()) * (micro.total - micro.absorption); - bool sampled = false; - - // Calculate elastic cross section if it wasn't precalculated - if (micro.elastic == CACHE_INVALID) { - nuc->calculate_elastic_xs(p); - } + void scatter(Particle & p, int i_nuclide) + { + // copy incoming direction + Direction u_old {p.u()}; - double prob = micro.elastic - micro.thermal; - if (prob > cutoff) { - // ======================================================================= - // NON-S(A,B) ELASTIC SCATTERING + // Get pointer to nuclide and grid index/interpolation factor + const auto& nuc {data::nuclides[i_nuclide]}; + const auto& micro {p.neutron_xs(i_nuclide)}; + int i_temp = micro.index_temp; - // Determine temperature - double kT = nuc->multipole_ ? p.sqrtkT() * p.sqrtkT() : nuc->kTs_[i_temp]; + // For tallying purposes, this routine might be called directly. In that + // case, we need to sample a reaction via the cutoff variable + double cutoff = prn(p.current_seed()) * (micro.total - micro.absorption); + bool sampled = false; - // Perform collision physics for elastic scattering - elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p); + // Calculate elastic cross section if it wasn't precalculated + if (micro.elastic == CACHE_INVALID) { + nuc->calculate_elastic_xs(p); + } - p.event_mt() = ELASTIC; - sampled = true; - } + double prob = micro.elastic - micro.thermal; + if (prob > cutoff) { + // ======================================================================= + // NON-S(A,B) ELASTIC SCATTERING - prob = micro.elastic; - if (prob > cutoff && !sampled) { - // ======================================================================= - // S(A,B) SCATTERING + // Determine temperature + double kT = nuc->multipole_ ? p.sqrtkT() * p.sqrtkT() : nuc->kTs_[i_temp]; - sab_scatter(i_nuclide, micro.index_sab, p); + // Perform collision physics for elastic scattering + elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p); - p.event_mt() = ELASTIC; - sampled = true; - } + p.event_mt() = ELASTIC; + sampled = true; + } - if (!sampled) { - // ======================================================================= - // INELASTIC SCATTERING + prob = micro.elastic; + if (prob > cutoff && !sampled) { + // ======================================================================= + // S(A,B) SCATTERING - int n = nuc->index_inelastic_scatter_.size(); - int i = 0; - for (int j = 0; j < n && prob < cutoff; ++j) { - i = nuc->index_inelastic_scatter_[j]; + sab_scatter(i_nuclide, micro.index_sab, p); - // add to cumulative probability - prob += nuc->reactions_[i]->xs(micro); + p.event_mt() = ELASTIC; + sampled = true; } - // Perform collision physics for inelastic scattering - const auto& rx {nuc->reactions_[i]}; - inelastic_scatter(*nuc, *rx, p); - p.event_mt() = rx->mt_; - } + if (!sampled) { + // ======================================================================= + // INELASTIC SCATTERING - // Set event component - p.event() = TallyEvent::SCATTER; - - // Sample new outgoing angle for isotropic-in-lab scattering - const auto& mat {model::materials[p.material()]}; - if (!mat->p0_.empty()) { - int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide]; - if (mat->p0_[i_nuc_mat]) { - // Sample isotropic-in-lab outgoing direction - p.u() = isotropic_direction(p.current_seed()); - p.mu() = u_old.dot(p.u()); - } - } -} + int n = nuc->index_inelastic_scatter_.size(); + int i = 0; + for (int j = 0; j < n && prob < cutoff; ++j) { + i = nuc->index_inelastic_scatter_[j]; -void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p) -{ - // get pointer to nuclide - const auto& nuc {data::nuclides[i_nuclide]}; + // add to cumulative probability + prob += nuc->reactions_[i]->xs(micro); + } - double vel = std::sqrt(p.E()); - double awr = nuc->awr_; + // Perform collision physics for inelastic scattering + const auto& rx {nuc->reactions_[i]}; + inelastic_scatter(*nuc, *rx, p); + p.event_mt() = rx->mt_; + } - // Neutron velocity in LAB - Direction v_n = vel * p.u(); + // Set event component + p.event() = TallyEvent::SCATTER; - // Sample velocity of target nucleus - Direction v_t {}; - if (!p.neutron_xs(i_nuclide).use_ptable) { - v_t = sample_target_velocity(*nuc, p.E(), p.u(), v_n, - p.neutron_xs(i_nuclide).elastic, kT, p.current_seed()); + // Sample new outgoing angle for isotropic-in-lab scattering + const auto& mat {model::materials[p.material()]}; + if (!mat->p0_.empty()) { + int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide]; + if (mat->p0_[i_nuc_mat]) { + // Sample isotropic-in-lab outgoing direction + p.u() = isotropic_direction(p.current_seed()); + p.mu() = u_old.dot(p.u()); + } + } } - // Velocity of center-of-mass - Direction v_cm = (v_n + awr * v_t) / (awr + 1.0); + void elastic_scatter( + int i_nuclide, const Reaction& rx, double kT, Particle& p) + { + // get pointer to nuclide + const auto& nuc {data::nuclides[i_nuclide]}; - // Transform to CM frame - v_n -= v_cm; + double vel = std::sqrt(p.E()); + double awr = nuc->awr_; - // Find speed of neutron in CM - vel = v_n.norm(); + // Neutron velocity in LAB + Direction v_n = vel * p.u(); - // Sample scattering angle, checking if angle distribution is present (assume - // isotropic otherwise) - double mu_cm; - auto& d = rx.products_[0].distribution_[0]; - auto d_ = dynamic_cast(d.get()); - if (!d_->angle().empty()) { - mu_cm = d_->angle().sample(p.E(), p.current_seed()); - } else { - mu_cm = uniform_distribution(-1., 1., p.current_seed()); - } + // Sample velocity of target nucleus + Direction v_t {}; + if (!p.neutron_xs(i_nuclide).use_ptable) { + v_t = sample_target_velocity(*nuc, p.E(), p.u(), v_n, + p.neutron_xs(i_nuclide).elastic, kT, p.current_seed()); + } - // Determine direction cosines in CM - Direction u_cm = v_n / vel; + // Velocity of center-of-mass + Direction v_cm = (v_n + awr * v_t) / (awr + 1.0); - // Rotate neutron velocity vector to new angle -- note that the speed of the - // neutron in CM does not change in elastic scattering. However, the speed - // will change when we convert back to LAB - v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p.current_seed()); + // Transform to CM frame + v_n -= v_cm; - // Transform back to LAB frame - v_n += v_cm; + // Find speed of neutron in CM + vel = v_n.norm(); - p.E() = v_n.dot(v_n); - vel = std::sqrt(p.E()); + // Sample scattering angle, checking if angle distribution is present + // (assume isotropic otherwise) + double mu_cm; + auto& d = rx.products_[0].distribution_[0]; + auto d_ = dynamic_cast(d.get()); + if (!d_->angle().empty()) { + mu_cm = d_->angle().sample(p.E(), p.current_seed()); + } else { + mu_cm = uniform_distribution(-1., 1., p.current_seed()); + } - // compute cosine of scattering angle in LAB frame by taking dot product of - // neutron's pre- and post-collision angle - p.mu() = p.u().dot(v_n) / vel; + // Determine direction cosines in CM + Direction u_cm = v_n / vel; - // Set energy and direction of particle in LAB frame - p.u() = v_n / vel; + // Rotate neutron velocity vector to new angle -- note that the speed of the + // neutron in CM does not change in elastic scattering. However, the speed + // will change when we convert back to LAB + v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p.current_seed()); - // Because of floating-point roundoff, it may be possible for mu_lab to be - // outside of the range [-1,1). In these cases, we just set mu_lab to exactly - // -1 or 1 - if (std::abs(p.mu()) > 1.0) - p.mu() = std::copysign(1.0, p.mu()); -} + // Transform back to LAB frame + v_n += v_cm; -void sab_scatter(int i_nuclide, int i_sab, Particle& p) -{ - // Determine temperature index - const auto& micro {p.neutron_xs(i_nuclide)}; - int i_temp = micro.index_temp_sab; - - // Sample energy and angle - double E_out; - data::thermal_scatt[i_sab]->data_[i_temp].sample( - micro, p.E(), &E_out, &p.mu(), p.current_seed()); - - // Set energy to outgoing, change direction of particle - p.E() = E_out; - p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); -} - -Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, - Direction v_neut, double xs_eff, double kT, uint64_t* seed) -{ - // check if nuclide is a resonant scatterer - ResScatMethod sampling_method; - if (nuc.resonant_) { + p.E() = v_n.dot(v_n); + vel = std::sqrt(p.E()); - // sampling method to use - sampling_method = settings::res_scat_method; + // compute cosine of scattering angle in LAB frame by taking dot product of + // neutron's pre- and post-collision angle + p.mu() = p.u().dot(v_n) / vel; - // upper resonance scattering energy bound (target is at rest above this E) - if (E > settings::res_scat_energy_max) { - return {}; + // Set energy and direction of particle in LAB frame + p.u() = v_n / vel; - // lower resonance scattering energy bound (should be no resonances below) - } else if (E < settings::res_scat_energy_min) { - sampling_method = ResScatMethod::cxs; - } + // Because of floating-point roundoff, it may be possible for mu_lab to be + // outside of the range [-1,1). In these cases, we just set mu_lab to + // exactly -1 or 1 + if (std::abs(p.mu()) > 1.0) + p.mu() = std::copysign(1.0, p.mu()); + } - // otherwise, use free gas model - } else { - if (E >= settings::free_gas_threshold * kT && nuc.awr_ > 1.0) { - return {}; - } else { - sampling_method = ResScatMethod::cxs; - } + void sab_scatter(int i_nuclide, int i_sab, Particle& p) + { + // Determine temperature index + const auto& micro {p.neutron_xs(i_nuclide)}; + int i_temp = micro.index_temp_sab; + + // Sample energy and angle + double E_out; + data::thermal_scatt[i_sab]->data_[i_temp].sample( + micro, p.E(), &E_out, &p.mu(), p.current_seed()); + + // Set energy to outgoing, change direction of particle + p.E() = E_out; + p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); } - // use appropriate target velocity sampling method - switch (sampling_method) { - case ResScatMethod::cxs: - - // sample target velocity with the constant cross section (cxs) approx. - return sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); - - case ResScatMethod::dbrc: - case ResScatMethod::rvs: { - double E_red = std::sqrt(nuc.awr_ * E / kT); - double E_low = std::pow(std::max(0.0, E_red - 4.0), 2) * kT / nuc.awr_; - double E_up = (E_red + 4.0) * (E_red + 4.0) * kT / nuc.awr_; - - // find lower and upper energy bound indices - // lower index - int i_E_low; - if (E_low < nuc.energy_0K_.front()) { - i_E_low = 0; - } else if (E_low > nuc.energy_0K_.back()) { - i_E_low = nuc.energy_0K_.size() - 2; - } else { - i_E_low = - lower_bound_index(nuc.energy_0K_.begin(), nuc.energy_0K_.end(), E_low); - } + Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, + Direction v_neut, double xs_eff, double kT, uint64_t* seed) + { + // check if nuclide is a resonant scatterer + ResScatMethod sampling_method; + if (nuc.resonant_) { + + // sampling method to use + sampling_method = settings::res_scat_method; + + // upper resonance scattering energy bound (target is at rest above this + // E) + if (E > settings::res_scat_energy_max) { + return {}; + + // lower resonance scattering energy bound (should be no resonances + // below) + } else if (E < settings::res_scat_energy_min) { + sampling_method = ResScatMethod::cxs; + } - // upper index - int i_E_up; - if (E_up < nuc.energy_0K_.front()) { - i_E_up = 0; - } else if (E_up > nuc.energy_0K_.back()) { - i_E_up = nuc.energy_0K_.size() - 2; + // otherwise, use free gas model } else { - i_E_up = - lower_bound_index(nuc.energy_0K_.begin(), nuc.energy_0K_.end(), E_up); + if (E >= settings::free_gas_threshold * kT && nuc.awr_ > 1.0) { + return {}; + } else { + sampling_method = ResScatMethod::cxs; + } } - if (i_E_up == i_E_low) { - // Handle degenerate case -- if the upper/lower bounds occur for the same - // index, then using cxs is probably a good approximation + // use appropriate target velocity sampling method + switch (sampling_method) { + case ResScatMethod::cxs: + + // sample target velocity with the constant cross section (cxs) approx. return sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); - } - if (sampling_method == ResScatMethod::dbrc) { - // interpolate xs since we're not exactly at the energy indices - double xs_low = nuc.elastic_0K_[i_E_low]; - double m = (nuc.elastic_0K_[i_E_low + 1] - xs_low) / - (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]); - xs_low += m * (E_low - nuc.energy_0K_[i_E_low]); - double xs_up = nuc.elastic_0K_[i_E_up]; - m = (nuc.elastic_0K_[i_E_up + 1] - xs_up) / - (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]); - xs_up += m * (E_up - nuc.energy_0K_[i_E_up]); - - // get max 0K xs value over range of practical relative energies - double xs_max = *std::max_element( - &nuc.elastic_0K_[i_E_low + 1], &nuc.elastic_0K_[i_E_up + 1]); - xs_max = std::max({xs_low, xs_max, xs_up}); - - while (true) { - double E_rel; - Direction v_target; - while (true) { - // sample target velocity with the constant cross section (cxs) - // approx. - v_target = sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); - Direction v_rel = v_neut - v_target; - E_rel = v_rel.dot(v_rel); - if (E_rel < E_up) - break; - } + case ResScatMethod::dbrc: + case ResScatMethod::rvs: { + double E_red = std::sqrt(nuc.awr_ * E / kT); + double E_low = std::pow(std::max(0.0, E_red - 4.0), 2) * kT / nuc.awr_; + double E_up = (E_red + 4.0) * (E_red + 4.0) * kT / nuc.awr_; + + // find lower and upper energy bound indices + // lower index + int i_E_low; + if (E_low < nuc.energy_0K_.front()) { + i_E_low = 0; + } else if (E_low > nuc.energy_0K_.back()) { + i_E_low = nuc.energy_0K_.size() - 2; + } else { + i_E_low = lower_bound_index( + nuc.energy_0K_.begin(), nuc.energy_0K_.end(), E_low); + } - // perform Doppler broadening rejection correction (dbrc) - double xs_0K = nuc.elastic_xs_0K(E_rel); - double R = xs_0K / xs_max; - if (prn(seed) < R) - return v_target; + // upper index + int i_E_up; + if (E_up < nuc.energy_0K_.front()) { + i_E_up = 0; + } else if (E_up > nuc.energy_0K_.back()) { + i_E_up = nuc.energy_0K_.size() - 2; + } else { + i_E_up = + lower_bound_index(nuc.energy_0K_.begin(), nuc.energy_0K_.end(), E_up); } - } else if (sampling_method == ResScatMethod::rvs) { - // interpolate xs CDF since we're not exactly at the energy indices - // cdf value at lower bound attainable energy - double cdf_low = 0.0; - if (E_low > nuc.energy_0K_.front()) { - double m = (nuc.xs_cdf_[i_E_low + 1] - nuc.xs_cdf_[i_E_low]) / - (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]); - cdf_low = nuc.xs_cdf_[i_E_low] + m * (E_low - nuc.energy_0K_[i_E_low]); + if (i_E_up == i_E_low) { + // Handle degenerate case -- if the upper/lower bounds occur for the + // same index, then using cxs is probably a good approximation + return sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); } - // cdf value at upper bound attainable energy - double m = (nuc.xs_cdf_[i_E_up + 1] - nuc.xs_cdf_[i_E_up]) / - (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]); - double cdf_up = nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]); - - while (true) { - // directly sample Maxwellian - double E_t = -kT * std::log(prn(seed)); - - // sample a relative energy using the xs cdf - double cdf_rel = cdf_low + prn(seed) * (cdf_up - cdf_low); - int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low, - nuc.xs_cdf_.begin() + i_E_up + 2, cdf_rel); - double E_rel = nuc.energy_0K_[i_E_low + i_E_rel]; - double m = (nuc.xs_cdf_[i_E_low + i_E_rel + 1] - - nuc.xs_cdf_[i_E_low + i_E_rel]) / - (nuc.energy_0K_[i_E_low + i_E_rel + 1] - - nuc.energy_0K_[i_E_low + i_E_rel]); - E_rel += (cdf_rel - nuc.xs_cdf_[i_E_low + i_E_rel]) / m; - - // perform rejection sampling on cosine between - // neutron and target velocities - double mu = (E_t + nuc.awr_ * (E - E_rel)) / - (2.0 * std::sqrt(nuc.awr_ * E * E_t)); - - if (std::abs(mu) < 1.0) { - // set and accept target velocity - E_t /= nuc.awr_; - return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, seed); + if (sampling_method == ResScatMethod::dbrc) { + // interpolate xs since we're not exactly at the energy indices + double xs_low = nuc.elastic_0K_[i_E_low]; + double m = (nuc.elastic_0K_[i_E_low + 1] - xs_low) / + (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]); + xs_low += m * (E_low - nuc.energy_0K_[i_E_low]); + double xs_up = nuc.elastic_0K_[i_E_up]; + m = (nuc.elastic_0K_[i_E_up + 1] - xs_up) / + (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]); + xs_up += m * (E_up - nuc.energy_0K_[i_E_up]); + + // get max 0K xs value over range of practical relative energies + double xs_max = *std::max_element( + &nuc.elastic_0K_[i_E_low + 1], &nuc.elastic_0K_[i_E_up + 1]); + xs_max = std::max({xs_low, xs_max, xs_up}); + + while (true) { + double E_rel; + Direction v_target; + while (true) { + // sample target velocity with the constant cross section (cxs) + // approx. + v_target = sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); + Direction v_rel = v_neut - v_target; + E_rel = v_rel.dot(v_rel); + if (E_rel < E_up) + break; + } + + // perform Doppler broadening rejection correction (dbrc) + double xs_0K = nuc.elastic_xs_0K(E_rel); + double R = xs_0K / xs_max; + if (prn(seed) < R) + return v_target; } - } - } - } // case RVS, DBRC - } // switch (sampling_method) - UNREACHABLE(); -} + } else if (sampling_method == ResScatMethod::rvs) { + // interpolate xs CDF since we're not exactly at the energy indices + // cdf value at lower bound attainable energy + double cdf_low = 0.0; + if (E_low > nuc.energy_0K_.front()) { + double m = (nuc.xs_cdf_[i_E_low + 1] - nuc.xs_cdf_[i_E_low]) / + (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]); + cdf_low = + nuc.xs_cdf_[i_E_low] + m * (E_low - nuc.energy_0K_[i_E_low]); + } -Direction sample_cxs_target_velocity( - double awr, double E, Direction u, double kT, uint64_t* seed) -{ - double beta_vn = std::sqrt(awr * E / kT); - double alpha = 1.0 / (1.0 + std::sqrt(PI) * beta_vn / 2.0); + // cdf value at upper bound attainable energy + double m = (nuc.xs_cdf_[i_E_up + 1] - nuc.xs_cdf_[i_E_up]) / + (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]); + double cdf_up = + nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]); - double beta_vt_sq; - double mu; - while (true) { - // Sample two random numbers - double r1 = prn(seed); - double r2 = prn(seed); + while (true) { + // directly sample Maxwellian + double E_t = -kT * std::log(prn(seed)); + + // sample a relative energy using the xs cdf + double cdf_rel = cdf_low + prn(seed) * (cdf_up - cdf_low); + int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low, + nuc.xs_cdf_.begin() + i_E_up + 2, cdf_rel); + double E_rel = nuc.energy_0K_[i_E_low + i_E_rel]; + double m = (nuc.xs_cdf_[i_E_low + i_E_rel + 1] - + nuc.xs_cdf_[i_E_low + i_E_rel]) / + (nuc.energy_0K_[i_E_low + i_E_rel + 1] - + nuc.energy_0K_[i_E_low + i_E_rel]); + E_rel += (cdf_rel - nuc.xs_cdf_[i_E_low + i_E_rel]) / m; + + // perform rejection sampling on cosine between + // neutron and target velocities + double mu = (E_t + nuc.awr_ * (E - E_rel)) / + (2.0 * std::sqrt(nuc.awr_ * E * E_t)); + + if (std::abs(mu) < 1.0) { + // set and accept target velocity + E_t /= nuc.awr_; + return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, seed); + } + } + } + } // case RVS, DBRC + } // switch (sampling_method) - if (prn(seed) < alpha) { - // With probability alpha, we sample the distribution p(y) = - // y*e^(-y). This can be done with sampling scheme C45 from the Monte - // Carlo sampler + UNREACHABLE(); + } - beta_vt_sq = -std::log(r1 * r2); + Direction sample_cxs_target_velocity( + double awr, double E, Direction u, double kT, uint64_t* seed) + { + double beta_vn = std::sqrt(awr * E / kT); + double alpha = 1.0 / (1.0 + std::sqrt(PI) * beta_vn / 2.0); - } else { - // With probability 1-alpha, we sample the distribution p(y) = y^2 * - // e^(-y^2). This can be done with sampling scheme C61 from the Monte - // Carlo sampler + double beta_vt_sq; + double mu; + while (true) { + // Sample two random numbers + double r1 = prn(seed); + double r2 = prn(seed); - double c = std::cos(PI / 2.0 * prn(seed)); - beta_vt_sq = -std::log(r1) - std::log(r2) * c * c; - } + if (prn(seed) < alpha) { + // With probability alpha, we sample the distribution p(y) = + // y*e^(-y). This can be done with sampling scheme C45 from the Monte + // Carlo sampler - // Determine beta * vt - double beta_vt = std::sqrt(beta_vt_sq); + beta_vt_sq = -std::log(r1 * r2); - // Sample cosine of angle between neutron and target velocity - mu = uniform_distribution(-1., 1., seed); + } else { + // With probability 1-alpha, we sample the distribution p(y) = y^2 * + // e^(-y^2). This can be done with sampling scheme C61 from the Monte + // Carlo sampler - // Determine rejection probability - double accept_prob = - std::sqrt(beta_vn * beta_vn + beta_vt_sq - 2 * beta_vn * beta_vt * mu) / - (beta_vn + beta_vt); + double c = std::cos(PI / 2.0 * prn(seed)); + beta_vt_sq = -std::log(r1) - std::log(r2) * c * c; + } - // Perform rejection sampling on vt and mu - if (prn(seed) < accept_prob) - break; - } + // Determine beta * vt + double beta_vt = std::sqrt(beta_vt_sq); - // Determine speed of target nucleus - double vt = std::sqrt(beta_vt_sq * kT / awr); + // Sample cosine of angle between neutron and target velocity + mu = uniform_distribution(-1., 1., seed); - // Determine velocity vector of target nucleus based on neutron's velocity - // and the sampled angle between them - return vt * rotate_angle(u, mu, nullptr, seed); -} + // Determine rejection probability + double accept_prob = + std::sqrt(beta_vn * beta_vn + beta_vt_sq - 2 * beta_vn * beta_vt * mu) / + (beta_vn + beta_vt); -void sample_fission_neutron( - int i_nuclide, const Reaction& rx, SourceSite* site, Particle& p) -{ - // Get attributes of particle - double E_in = p.E(); - uint64_t* seed = p.current_seed(); - - // Determine total nu, delayed nu, and delayed neutron fraction - const auto& nuc {data::nuclides[i_nuclide]}; - double nu_t = nuc->nu(E_in, Nuclide::EmissionMode::total); - double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed); - double beta = nu_d / nu_t; - - if (prn(seed) < beta) { - // ==================================================================== - // DELAYED NEUTRON SAMPLED - - // sampled delayed precursor group - double xi = prn(seed) * nu_d; - double prob = 0.0; - int group; - for (group = 1; group < nuc->n_precursor_; ++group) { - // determine delayed neutron precursor yield for group j - double yield = (*rx.products_[group].yield_)(E_in); - - // Check if this group is sampled - prob += yield; - if (xi < prob) + // Perform rejection sampling on vt and mu + if (prn(seed) < accept_prob) break; } - // if the sum of the probabilities is slightly less than one and the - // random number is greater, j will be greater than nuc % - // n_precursor -- check for this condition - group = std::min(group, nuc->n_precursor_); + // Determine speed of target nucleus + double vt = std::sqrt(beta_vt_sq * kT / awr); + + // Determine velocity vector of target nucleus based on neutron's velocity + // and the sampled angle between them + return vt * rotate_angle(u, mu, nullptr, seed); + } + + void sample_fission_neutron( + int i_nuclide, const Reaction& rx, SourceSite* site, Particle& p) + { + // Get attributes of particle + double E_in = p.E(); + uint64_t* seed = p.current_seed(); + + // Determine total nu, delayed nu, and delayed neutron fraction + const auto& nuc {data::nuclides[i_nuclide]}; + double nu_t = nuc->nu(E_in, Nuclide::EmissionMode::total); + double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed); + double beta = nu_d / nu_t; + + if (prn(seed) < beta) { + // ==================================================================== + // DELAYED NEUTRON SAMPLED + + // sampled delayed precursor group + double xi = prn(seed) * nu_d; + double prob = 0.0; + int group; + for (group = 1; group < nuc->n_precursor_; ++group) { + // determine delayed neutron precursor yield for group j + double yield = (*rx.products_[group].yield_)(E_in); + + // Check if this group is sampled + prob += yield; + if (xi < prob) + break; + } - // set the delayed group for the particle born from fission - site->delayed_group = group; + // if the sum of the probabilities is slightly less than one and the + // random number is greater, j will be greater than nuc % + // n_precursor -- check for this condition + group = std::min(group, nuc->n_precursor_); - // Sample time of emission based on decay constant of precursor - double decay_rate = rx.products_[site->delayed_group].decay_rate_; - site->time -= std::log(prn(p.current_seed())) / decay_rate; + // set the delayed group for the particle born from fission + site->delayed_group = group; - } else { - // ==================================================================== - // PROMPT NEUTRON SAMPLED + // Sample time of emission based on decay constant of precursor + double decay_rate = rx.products_[site->delayed_group].decay_rate_; + site->time -= std::log(prn(p.current_seed())) / decay_rate; - // set the delayed group for the particle born from fission to 0 - site->delayed_group = 0; - } + } else { + // ==================================================================== + // PROMPT NEUTRON SAMPLED - // sample from prompt neutron energy distribution - int n_sample = 0; - double mu; - while (true) { - rx.products_[site->delayed_group].sample(E_in, site->E, mu, seed); - - // resample if energy is greater than maximum neutron energy - constexpr int neutron = static_cast(ParticleType::neutron); - if (site->E < data::energy_max[neutron]) - break; - - // check for large number of resamples - ++n_sample; - if (n_sample == MAX_SAMPLE) { - // particle_write_restart(p) - fatal_error("Resampled energy distribution maximum number of times " - "for nuclide " + - nuc->name_); + // set the delayed group for the particle born from fission to 0 + site->delayed_group = 0; } - } - // Sample azimuthal angle uniformly in [0, 2*pi) and assign angle - site->u = rotate_angle(p.u(), mu, nullptr, seed); -} + // sample from prompt neutron energy distribution + int n_sample = 0; + double mu; + while (true) { + rx.products_[site->delayed_group].sample(E_in, site->E, mu, seed); -void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) -{ - // copy energy of neutron - double E_in = p.E(); - - // sample outgoing energy and scattering cosine - double E; - double mu; - rx.products_[0].sample(E_in, E, mu, p.current_seed()); - - // if scattering system is in center-of-mass, transfer cosine of scattering - // angle and outgoing energy from CM to LAB - if (rx.scatter_in_cm_) { - double E_cm = E; - - // determine outgoing energy in lab - double A = nuc.awr_; - E = E_cm + (E_in + 2.0 * mu * (A + 1.0) * std::sqrt(E_in * E_cm)) / - ((A + 1.0) * (A + 1.0)); - - // determine outgoing angle in lab - mu = mu * std::sqrt(E_cm / E) + 1.0 / (A + 1.0) * std::sqrt(E_in / E); - } + // resample if energy is greater than maximum neutron energy + constexpr int neutron = static_cast(ParticleType::neutron); + if (site->E < data::energy_max[neutron]) + break; - // Because of floating-point roundoff, it may be possible for mu to be - // outside of the range [-1,1). In these cases, we just set mu to exactly -1 - // or 1 - if (std::abs(mu) > 1.0) - mu = std::copysign(1.0, mu); - - // Set outgoing energy and scattering angle - p.E() = E; - p.mu() = mu; - - // change direction of particle - p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed()); - - // evaluate yield - double yield = (*rx.products_[0].yield_)(E_in); - if (std::floor(yield) == yield && yield > 0) { - // If yield is integral, create exactly that many secondary particles - for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { - p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron); + // check for large number of resamples + ++n_sample; + if (n_sample == MAX_SAMPLE) { + // particle_write_restart(p) + fatal_error("Resampled energy distribution maximum number of times " + "for nuclide " + + nuc->name_); + } } - } else { - // Otherwise, change weight of particle based on yield - p.wgt() *= yield; - } -} -void sample_secondary_photons(Particle& p, int i_nuclide) -{ - // Sample the number of photons produced - double y_t = - p.neutron_xs(i_nuclide).photon_prod / p.neutron_xs(i_nuclide).total; - double photon_wgt = p.wgt(); - int y = 1; - - if (settings::use_decay_photons) { - // For decay photons, sample a single photon and modify the weight - if (y_t <= 0.0) - return; - photon_wgt *= y_t; - } else { - // For prompt photons, sample an integral number of photons with weight - // equal to the neutron's weight - y = static_cast(y_t); - if (prn(p.current_seed()) <= y_t - y) - ++y; + // Sample azimuthal angle uniformly in [0, 2*pi) and assign angle + site->u = rotate_angle(p.u(), mu, nullptr, seed); } - // Sample each secondary photon - for (int i = 0; i < y; ++i) { - // Sample the reaction and product - int i_rx; - int i_product; - sample_photon_product(i_nuclide, p, &i_rx, &i_product); + void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) + { + // copy energy of neutron + double E_in = p.E(); - // Sample the outgoing energy and angle - auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx]; + // sample outgoing energy and scattering cosine double E; double mu; - rx->products_[i_product].sample(p.E(), E, mu, p.current_seed()); - - // Sample the new direction - Direction u = rotate_angle(p.u(), mu, nullptr, p.current_seed()); - - // In a k-eigenvalue simulation, it's necessary to provide higher weight to - // secondary photons from non-fission reactions to properly balance energy - // release and deposition. See D. P. Griesheimer, S. J. Douglass, and M. H. - // Stedry, "Self-consistent energy normalization for quasistatic reactor - // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. - double wgt = photon_wgt; - if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) { - wgt *= simulation::keff; + rx.products_[0].sample(E_in, E, mu, p.current_seed()); + + // if scattering system is in center-of-mass, transfer cosine of scattering + // angle and outgoing energy from CM to LAB + if (rx.scatter_in_cm_) { + double E_cm = E; + + // determine outgoing energy in lab + double A = nuc.awr_; + E = E_cm + (E_in + 2.0 * mu * (A + 1.0) * std::sqrt(E_in * E_cm)) / + ((A + 1.0) * (A + 1.0)); + + // determine outgoing angle in lab + mu = mu * std::sqrt(E_cm / E) + 1.0 / (A + 1.0) * std::sqrt(E_in / E); + } + + // Because of floating-point roundoff, it may be possible for mu to be + // outside of the range [-1,1). In these cases, we just set mu to exactly -1 + // or 1 + if (std::abs(mu) > 1.0) + mu = std::copysign(1.0, mu); + + // Set outgoing energy and scattering angle + p.E() = E; + p.mu() = mu; + + // change direction of particle + p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed()); + + // evaluate yield + double yield = (*rx.products_[0].yield_)(E_in); + if (std::floor(yield) == yield && yield > 0) { + // If yield is integral, create exactly that many secondary particles + for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { + p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron); + } + } else { + // Otherwise, change weight of particle based on yield + p.wgt() *= yield; + } + } + + void sample_secondary_photons(Particle & p, int i_nuclide) + { + // Sample the number of photons produced + double y_t = + p.neutron_xs(i_nuclide).photon_prod / p.neutron_xs(i_nuclide).total; + double photon_wgt = p.wgt(); + int y = 1; + + if (settings::use_decay_photons) { + // For decay photons, sample a single photon and modify the weight + if (y_t <= 0.0) + return; + photon_wgt *= y_t; + } else { + // For prompt photons, sample an integral number of photons with weight + // equal to the neutron's weight + y = static_cast(y_t); + if (prn(p.current_seed()) <= y_t - y) + ++y; } - // Create the secondary photon - bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon); + // Sample each secondary photon + for (int i = 0; i < y; ++i) { + // Sample the reaction and product + int i_rx; + int i_product; + sample_photon_product(i_nuclide, p, &i_rx, &i_product); + + // Sample the outgoing energy and angle + auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx]; + double E; + double mu; + rx->products_[i_product].sample(p.E(), E, mu, p.current_seed()); + + // Sample the new direction + Direction u = rotate_angle(p.u(), mu, nullptr, p.current_seed()); + + // In a k-eigenvalue simulation, it's necessary to provide higher weight + // to secondary photons from non-fission reactions to properly balance + // energy release and deposition. See D. P. Griesheimer, S. J. Douglass, + // and M. H. Stedry, "Self-consistent energy normalization for quasistatic + // reactor calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. + double wgt = photon_wgt; + if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) { + wgt *= simulation::keff; + } + + // Create the secondary photon + bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon); - // Tag secondary particle with parent nuclide - if (created_photon && settings::use_decay_photons) { - p.secondary_bank().back().parent_nuclide = - rx->products_[i_product].parent_nuclide_; + // Tag secondary particle with parent nuclide + if (created_photon && settings::use_decay_photons) { + p.secondary_bank().back().parent_nuclide = + rx->products_[i_product].parent_nuclide_; + } } } -} } // namespace openmc diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 9b0e423dff0..e2e3f687dd4 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -31,229 +31,240 @@ void collision_mg(Particle& p) // Sample the reaction type sample_reaction(p); - if (settings::weight_windows_on) { + // If collision checkpoint is enabled, apply weight window + // if valid and apply global russian roulette if not. + if (settings::weight_window_checkpoint_collision) { auto ww = search_weight_window(p); if (ww.is_valid()) { - if (settings::weight_window_checkpoint_collision) - apply_weight_window(p, ww); - } else if (p.type() == ParticleType::neutron) { - apply_russian_roulette(p); + apply_weight_window(p, ww); + else if (p.type() == ParticleType::neutron) + { + apply_russian_roulette(p) + } + // If collision checkpoint is disabled, apply russian roulette if + // the particle is outside the weight window domain + } else if (settings::weight_windows_on) { + auto ww = search_weight_window(p); + if (!ww.is_valid()) + apply_russian_roulette(p); } - } - // Display information about collision - if ((settings::verbosity >= 10) || p.trace()) { - write_message(fmt::format(" Energy Group = {}", p.g()), 1); + // Display information about collision + if ((settings::verbosity >= 10) || p.trace()) { + write_message(fmt::format(" Energy Group = {}", p.g()), 1); + } } -} -void sample_reaction(Particle& p) -{ - // Create fission bank sites. Note that while a fission reaction is sampled, - // it never actually "happens", i.e. the weight of the particle does not - // change when sampling fission sites. The following block handles all - // absorption (including fission) - - if (model::materials[p.material()]->fissionable()) { - if (settings::run_mode == RunMode::EIGENVALUE || - (settings::run_mode == RunMode::FIXED_SOURCE && - settings::create_fission_neutrons)) { - create_fission_sites(p); + void sample_reaction(Particle & p) + { + // Create fission bank sites. Note that while a fission reaction is sampled, + // it never actually "happens", i.e. the weight of the particle does not + // change when sampling fission sites. The following block handles all + // absorption (including fission) + + if (model::materials[p.material()]->fissionable()) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::create_fission_neutrons)) { + create_fission_sites(p); + } } - } - // If survival biasing is being used, the following subroutine adjusts the - // weight of the particle. Otherwise, it checks to see if absorption occurs. - if (p.macro_xs().absorption > 0.) { - absorption(p); + // If survival biasing is being used, the following subroutine adjusts the + // weight of the particle. Otherwise, it checks to see if absorption occurs. + if (p.macro_xs().absorption > 0.) { + absorption(p); + } + if (!p.alive()) + return; + + // Sample a scattering event to determine the energy of the exiting neutron + scatter(p); + + // Play russian roulette if there are no weight windows + if (!settings::weight_windows_on) + apply_russian_roulette(p); } - if (!p.alive()) - return; - // Sample a scattering event to determine the energy of the exiting neutron - scatter(p); + void scatter(Particle & p) + { + data::mg.macro_xs_[p.material()].sample_scatter(p.g_last(), p.g(), p.mu(), + p.wgt(), p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a); - // Play russian roulette if there are no weight windows - if (!settings::weight_windows_on) - apply_russian_roulette(p); -} + // Rotate the angle + p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); -void scatter(Particle& p) -{ - data::mg.macro_xs_[p.material()].sample_scatter(p.g_last(), p.g(), p.mu(), - p.wgt(), p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a); + // Update energy value for downstream compatability (in tallying) + p.E() = data::mg.energy_bin_avg_[p.g()]; - // Rotate the angle - p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); + // Set event component + p.event() = TallyEvent::SCATTER; + } - // Update energy value for downstream compatability (in tallying) - p.E() = data::mg.energy_bin_avg_[p.g()]; + void create_fission_sites(Particle & p) + { + // If uniform fission source weighting is turned on, we increase or decrease + // the expected number of fission sites produced + double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; - // Set event component - p.event() = TallyEvent::SCATTER; -} + // Determine the expected number of neutrons produced + double nu_t = p.wgt() / simulation::keff * weight * + p.macro_xs().nu_fission / p.macro_xs().total; -void create_fission_sites(Particle& p) -{ - // If uniform fission source weighting is turned on, we increase or decrease - // the expected number of fission sites produced - double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; - - // Determine the expected number of neutrons produced - double nu_t = p.wgt() / simulation::keff * weight * p.macro_xs().nu_fission / - p.macro_xs().total; - - // Sample the number of neutrons produced - int nu = static_cast(nu_t); - if (prn(p.current_seed()) <= (nu_t - int(nu_t))) { - nu++; - } + // Sample the number of neutrons produced + int nu = static_cast(nu_t); + if (prn(p.current_seed()) <= (nu_t - int(nu_t))) { + nu++; + } - // If no neutrons were produced then don't continue - if (nu == 0) - return; - - // Initialize the counter of delayed neutrons encountered for each delayed - // group. - double nu_d[MAX_DELAYED_GROUPS] = {0.}; - - // Clear out particle's nu fission bank - p.nu_bank().clear(); - - p.fission() = true; - - // Determine whether to place fission sites into the shared fission bank - // or the secondary particle bank. - bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); - - // Counter for the number of fission sites successfully stored to the shared - // fission bank or the secondary particle bank - int n_sites_stored; - - for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { - // Initialize fission site object with particle data - SourceSite site; - site.r = p.r(); - site.particle = ParticleType::neutron; - site.time = p.time(); - site.wgt = 1. / weight; - - // Sample the cosine of the angle, assuming fission neutrons are emitted - // isotropically - double mu = 2. * prn(p.current_seed()) - 1.; - - // Sample the azimuthal angle uniformly in [0, 2.pi) - double phi = 2. * PI * prn(p.current_seed()); - site.u.x = mu; - site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi); - site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi); - - // Sample secondary energy distribution for the fission reaction - int dg; - int gout; - data::mg.macro_xs_[p.material()].sample_fission_energy( - p.g(), dg, gout, p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a); - - // Store the energy and delayed groups on the fission bank - site.E = gout; - - // We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest - // of the code, 0 is prompt. - site.delayed_group = dg + 1; - - // If delayed product production, sample time of emission - if (dg != -1) { - auto& macro_xs = data::mg.macro_xs_[p.material()]; - double decay_rate = - macro_xs.get_xs(MgxsType::DECAY_RATE, 0, nullptr, nullptr, &dg, 0, 0); - site.time -= std::log(prn(p.current_seed())) / decay_rate; - - // Reject site if it exceeds time cutoff - double t_cutoff = settings::time_cutoff[static_cast(site.particle)]; - if (site.time > t_cutoff) { - continue; + // If no neutrons were produced then don't continue + if (nu == 0) + return; + + // Initialize the counter of delayed neutrons encountered for each delayed + // group. + double nu_d[MAX_DELAYED_GROUPS] = {0.}; + + // Clear out particle's nu fission bank + p.nu_bank().clear(); + + p.fission() = true; + + // Determine whether to place fission sites into the shared fission bank + // or the secondary particle bank. + bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); + + // Counter for the number of fission sites successfully stored to the shared + // fission bank or the secondary particle bank + int n_sites_stored; + + for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { + // Initialize fission site object with particle data + SourceSite site; + site.r = p.r(); + site.particle = ParticleType::neutron; + site.time = p.time(); + site.wgt = 1. / weight; + + // Sample the cosine of the angle, assuming fission neutrons are emitted + // isotropically + double mu = 2. * prn(p.current_seed()) - 1.; + + // Sample the azimuthal angle uniformly in [0, 2.pi) + double phi = 2. * PI * prn(p.current_seed()); + site.u.x = mu; + site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi); + site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi); + + // Sample secondary energy distribution for the fission reaction + int dg; + int gout; + data::mg.macro_xs_[p.material()].sample_fission_energy(p.g(), dg, gout, + p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a); + + // Store the energy and delayed groups on the fission bank + site.E = gout; + + // We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest + // of the code, 0 is prompt. + site.delayed_group = dg + 1; + + // If delayed product production, sample time of emission + if (dg != -1) { + auto& macro_xs = data::mg.macro_xs_[p.material()]; + double decay_rate = + macro_xs.get_xs(MgxsType::DECAY_RATE, 0, nullptr, nullptr, &dg, 0, 0); + site.time -= std::log(prn(p.current_seed())) / decay_rate; + + // Reject site if it exceeds time cutoff + double t_cutoff = + settings::time_cutoff[static_cast(site.particle)]; + if (site.time > t_cutoff) { + continue; + } } - } - // Set parent and progeny ID - site.parent_id = p.id(); - site.progeny_id = p.n_progeny()++; - - // Store fission site in bank - if (use_fission_bank) { - int64_t idx = simulation::fission_bank.thread_safe_append(site); - if (idx == -1) { - warning( - "The shared fission bank is full. Additional fission sites created " - "in this generation will not be banked. Results may be " - "non-deterministic."); - - // Decrement number of particle progeny as storage was unsuccessful. - // This step is needed so that the sum of all progeny is equal to the - // size of the shared fission bank. - p.n_progeny()--; - - // Break out of loop as no more sites can be added to fission bank - break; + // Set parent and progeny ID + site.parent_id = p.id(); + site.progeny_id = p.n_progeny()++; + + // Store fission site in bank + if (use_fission_bank) { + int64_t idx = simulation::fission_bank.thread_safe_append(site); + if (idx == -1) { + warning( + "The shared fission bank is full. Additional fission sites created " + "in this generation will not be banked. Results may be " + "non-deterministic."); + + // Decrement number of particle progeny as storage was unsuccessful. + // This step is needed so that the sum of all progeny is equal to the + // size of the shared fission bank. + p.n_progeny()--; + + // Break out of loop as no more sites can be added to fission bank + break; + } + } else { + p.secondary_bank().push_back(site); } - } else { - p.secondary_bank().push_back(site); - } - // Set the delayed group on the particle as well - p.delayed_group() = dg + 1; + // Set the delayed group on the particle as well + p.delayed_group() = dg + 1; - // Increment the number of neutrons born delayed - if (p.delayed_group() > 0) { - nu_d[dg]++; + // Increment the number of neutrons born delayed + if (p.delayed_group() > 0) { + nu_d[dg]++; + } + + // Write fission particles to nuBank + NuBank& nu_bank_entry = p.nu_bank().emplace_back(); + nu_bank_entry.wgt = site.wgt; + nu_bank_entry.E = site.E; + nu_bank_entry.delayed_group = site.delayed_group; } - // Write fission particles to nuBank - NuBank& nu_bank_entry = p.nu_bank().emplace_back(); - nu_bank_entry.wgt = site.wgt; - nu_bank_entry.E = site.E; - nu_bank_entry.delayed_group = site.delayed_group; - } + // If shared fission bank was full, and no fissions could be added, + // set the particle fission flag to false. + if (n_sites_stored == 0) { + p.fission() = false; + return; + } + + // Set nu to the number of fission sites successfully stored. If the fission + // bank was not found to be full then these values are already equivalent. + nu = n_sites_stored; - // If shared fission bank was full, and no fissions could be added, - // set the particle fission flag to false. - if (n_sites_stored == 0) { - p.fission() = false; - return; + // Store the total weight banked for analog fission tallies + p.n_bank() = nu; + p.wgt_bank() = nu / weight; + for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { + p.n_delayed_bank(d) = nu_d[d]; + } } - // Set nu to the number of fission sites successfully stored. If the fission - // bank was not found to be full then these values are already equivalent. - nu = n_sites_stored; + void absorption(Particle & p) + { + if (settings::survival_biasing) { + // Determine weight absorbed in survival biasing + double wgt_absorb = + p.wgt() * p.macro_xs().absorption / p.macro_xs().total; - // Store the total weight banked for analog fission tallies - p.n_bank() = nu; - p.wgt_bank() = nu / weight; - for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { - p.n_delayed_bank(d) = nu_d[d]; - } -} + // Adjust weight of particle by the probability of absorption + p.wgt() -= wgt_absorb; -void absorption(Particle& p) -{ - if (settings::survival_biasing) { - // Determine weight absorbed in survival biasing - double wgt_absorb = p.wgt() * p.macro_xs().absorption / p.macro_xs().total; - - // Adjust weight of particle by the probability of absorption - p.wgt() -= wgt_absorb; - - // Score implicit absorpion estimate of keff - p.keff_tally_absorption() += - wgt_absorb * p.macro_xs().nu_fission / p.macro_xs().absorption; - } else { - if (p.macro_xs().absorption > prn(p.current_seed()) * p.macro_xs().total) { + // Score implicit absorpion estimate of keff p.keff_tally_absorption() += - p.wgt() * p.macro_xs().nu_fission / p.macro_xs().absorption; - p.wgt() = 0.0; - p.event() = TallyEvent::ABSORB; + wgt_absorb * p.macro_xs().nu_fission / p.macro_xs().absorption; + } else { + if (p.macro_xs().absorption > + prn(p.current_seed()) * p.macro_xs().total) { + p.keff_tally_absorption() += + p.wgt() * p.macro_xs().nu_fission / p.macro_xs().absorption; + p.wgt() = 0.0; + p.event() = TallyEvent::ABSORB; + } } } -} } // namespace openmc From 11ce0f64f4d1e67a966bb29c66c54f01eb44663c Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 21:54:24 +0200 Subject: [PATCH 19/24] fix typos --- src/physics.cpp | 1967 ++++++++++++++++++++++---------------------- src/physics_mg.cpp | 394 +++++---- 2 files changed, 1175 insertions(+), 1186 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index ee94d068c78..299aba473a4 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -67,1176 +67,1169 @@ void collision(Particle& p) auto ww = search_weight_window(p); if (ww.is_valid()) { apply_weight_window(p, ww); - else if (p.type() == ParticleType::neutron) - { - apply_russian_roulette(p) - } - // If collision checkpoint is disabled, apply russian roulette if - // the particle is a neutron and it is outside the weight window domain - } else if (settings::weight_windows_on && - p.type() == ParticleType::neutron) { - auto ww = search_weight_window(p); - if (!ww.is_valid()) - apply_russian_roulette(p); + } else if (p.type() == ParticleType::neutron) { + apply_russian_roulette(p); } + // If collision checkpoint is disabled, apply russian roulette if + // the particle is a neutron and it is outside the weight window domain + } else if (settings::weight_windows_on && p.type() == ParticleType::neutron) { + auto ww = search_weight_window(p); + if (!ww.is_valid()) + apply_russian_roulette(p); + } - // Kill particle if energy falls below cutoff - int type = static_cast(p.type()); - if (p.E() < settings::energy_cutoff[type]) { - p.wgt() = 0.0; + // Kill particle if energy falls below cutoff + int type = static_cast(p.type()); + if (p.E() < settings::energy_cutoff[type]) { + p.wgt() = 0.0; + } + + // Display information about collision + if (settings::verbosity >= 10 || p.trace()) { + std::string msg; + if (p.event() == TallyEvent::KILL) { + msg = fmt::format(" Killed. Energy = {} eV.", p.E()); + } else if (p.type() == ParticleType::neutron) { + msg = fmt::format(" {} with {}. Energy = {} eV.", + reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_, + p.E()); + } else if (p.type() == ParticleType::photon) { + msg = fmt::format(" {} with {}. Energy = {} eV.", + reaction_name(p.event_mt()), + to_element(data::nuclides[p.event_nuclide()]->name_), p.E()); + } else { + msg = fmt::format(" Disappeared. Energy = {} eV.", p.E()); } + write_message(msg, 1); + } +} - // Display information about collision - if (settings::verbosity >= 10 || p.trace()) { - std::string msg; - if (p.event() == TallyEvent::KILL) { - msg = fmt::format(" Killed. Energy = {} eV.", p.E()); - } else if (p.type() == ParticleType::neutron) { - msg = fmt::format(" {} with {}. Energy = {} eV.", - reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_, - p.E()); - } else if (p.type() == ParticleType::photon) { - msg = fmt::format(" {} with {}. Energy = {} eV.", - reaction_name(p.event_mt()), - to_element(data::nuclides[p.event_nuclide()]->name_), p.E()); - } else { - msg = fmt::format(" Disappeared. Energy = {} eV.", p.E()); +void sample_neutron_reaction(Particle& p) +{ + // Sample a nuclide within the material + int i_nuclide = sample_nuclide(p); + + // Save which nuclide particle had collision with + p.event_nuclide() = i_nuclide; + + // Create fission bank sites. Note that while a fission reaction is sampled, + // it never actually "happens", i.e. the weight of the particle does not + // change when sampling fission sites. The following block handles all + // absorption (including fission) + + const auto& nuc {data::nuclides[i_nuclide]}; + + if (nuc->fissionable_ && p.neutron_xs(i_nuclide).fission > 0.0) { + auto& rx = sample_fission(i_nuclide, p); + if (settings::run_mode == RunMode::EIGENVALUE) { + create_fission_sites(p, i_nuclide, rx); + } else if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::create_fission_neutrons) { + create_fission_sites(p, i_nuclide, rx); + + // Make sure particle population doesn't grow out of control for + // subcritical multiplication problems. + if (p.secondary_bank().size() >= settings::max_secondaries) { + fatal_error( + "The secondary particle bank appears to be growing without " + "bound. You are likely running a subcritical multiplication " + "problem " + "with k-effective close to or greater than one."); } - write_message(msg, 1); } + p.event_mt() = rx.mt_; } - void sample_neutron_reaction(Particle & p) - { - // Sample a nuclide within the material - int i_nuclide = sample_nuclide(p); - - // Save which nuclide particle had collision with - p.event_nuclide() = i_nuclide; - - // Create fission bank sites. Note that while a fission reaction is sampled, - // it never actually "happens", i.e. the weight of the particle does not - // change when sampling fission sites. The following block handles all - // absorption (including fission) + // Create secondary photons + if (settings::photon_transport) { + sample_secondary_photons(p, i_nuclide); + } - const auto& nuc {data::nuclides[i_nuclide]}; + // If survival biasing is being used, the following subroutine adjusts the + // weight of the particle. Otherwise, it checks to see if absorption occurs - if (nuc->fissionable_ && p.neutron_xs(i_nuclide).fission > 0.0) { - auto& rx = sample_fission(i_nuclide, p); - if (settings::run_mode == RunMode::EIGENVALUE) { - create_fission_sites(p, i_nuclide, rx); - } else if (settings::run_mode == RunMode::FIXED_SOURCE && - settings::create_fission_neutrons) { - create_fission_sites(p, i_nuclide, rx); - - // Make sure particle population doesn't grow out of control for - // subcritical multiplication problems. - if (p.secondary_bank().size() >= settings::max_secondaries) { - fatal_error( - "The secondary particle bank appears to be growing without " - "bound. You are likely running a subcritical multiplication " - "problem " - "with k-effective close to or greater than one."); - } - } - p.event_mt() = rx.mt_; - } + if (p.neutron_xs(i_nuclide).absorption > 0.0) { + absorption(p, i_nuclide); + } + if (!p.alive()) + return; + + // Sample a scattering reaction and determine the secondary energy of the + // exiting neutron + const auto& ncrystal_mat = model::materials[p.material()]->ncrystal_mat(); + if (ncrystal_mat && p.E() < NCRYSTAL_MAX_ENERGY) { + ncrystal_mat.scatter(p); + } else { + scatter(p, i_nuclide); + } - // Create secondary photons - if (settings::photon_transport) { - sample_secondary_photons(p, i_nuclide); - } + // Advance URR seed stream 'N' times after energy changes + if (p.E() != p.E_last()) { + advance_prn_seed(data::nuclides.size(), &p.seeds(STREAM_URR_PTABLE)); + } - // If survival biasing is being used, the following subroutine adjusts the - // weight of the particle. Otherwise, it checks to see if absorption occurs + // Play russian roulette if there are no weight windows + if (!settings::weight_windows_on) + apply_russian_roulette(p); +} - if (p.neutron_xs(i_nuclide).absorption > 0.0) { - absorption(p, i_nuclide); +void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) +{ + // If uniform fission source weighting is turned on, we increase or decrease + // the expected number of fission sites produced + double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; + + // Determine the expected number of neutrons produced + double nu_t = p.wgt() / simulation::keff * weight * + p.neutron_xs(i_nuclide).nu_fission / + p.neutron_xs(i_nuclide).total; + + // Sample the number of neutrons produced + int nu = static_cast(nu_t); + if (prn(p.current_seed()) <= (nu_t - nu)) + ++nu; + + // If no neutrons were produced then don't continue + if (nu == 0) + return; + + // Initialize the counter of delayed neutrons encountered for each delayed + // group. + double nu_d[MAX_DELAYED_GROUPS] = {0.}; + + // Clear out particle's nu fission bank + p.nu_bank().clear(); + + p.fission() = true; + + // Determine whether to place fission sites into the shared fission bank + // or the secondary particle bank. + bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); + + // Counter for the number of fission sites successfully stored to the shared + // fission bank or the secondary particle bank + int n_sites_stored; + + for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { + // Initialize fission site object with particle data + SourceSite site; + site.r = p.r(); + site.particle = ParticleType::neutron; + site.time = p.time(); + site.wgt = 1. / weight; + site.surf_id = 0; + + // Sample delayed group and angle/energy for fission reaction + sample_fission_neutron(i_nuclide, rx, &site, p); + + // Reject site if it exceeds time cutoff + if (site.delayed_group > 0) { + double t_cutoff = settings::time_cutoff[static_cast(site.particle)]; + if (site.time > t_cutoff) { + continue; + } } - if (!p.alive()) - return; - // Sample a scattering reaction and determine the secondary energy of the - // exiting neutron - const auto& ncrystal_mat = model::materials[p.material()]->ncrystal_mat(); - if (ncrystal_mat && p.E() < NCRYSTAL_MAX_ENERGY) { - ncrystal_mat.scatter(p); + // Set parent and progeny IDs + site.parent_id = p.id(); + site.progeny_id = p.n_progeny()++; + + // Store fission site in bank + if (use_fission_bank) { + int64_t idx = simulation::fission_bank.thread_safe_append(site); + if (idx == -1) { + warning( + "The shared fission bank is full. Additional fission sites created " + "in this generation will not be banked. Results may be " + "non-deterministic."); + + // Decrement number of particle progeny as storage was unsuccessful. + // This step is needed so that the sum of all progeny is equal to the + // size of the shared fission bank. + p.n_progeny()--; + + // 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, idx); + } } else { - scatter(p, i_nuclide); + p.secondary_bank().push_back(site); } - // Advance URR seed stream 'N' times after energy changes - if (p.E() != p.E_last()) { - advance_prn_seed(data::nuclides.size(), &p.seeds(STREAM_URR_PTABLE)); + // Increment the number of neutrons born delayed + if (site.delayed_group > 0) { + nu_d[site.delayed_group - 1]++; } - // Play russian roulette if there are no weight windows - if (!settings::weight_windows_on) - apply_russian_roulette(p); + // Write fission particles to nuBank + NuBank& nu_bank_entry = p.nu_bank().emplace_back(); + nu_bank_entry.wgt = site.wgt; + nu_bank_entry.E = site.E; + nu_bank_entry.delayed_group = site.delayed_group; } - void create_fission_sites(Particle & p, int i_nuclide, const Reaction& rx) - { - // If uniform fission source weighting is turned on, we increase or decrease - // the expected number of fission sites produced - double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; + // If shared fission bank was full, and no fissions could be added, + // set the particle fission flag to false. + if (n_sites_stored == 0) { + p.fission() = false; + return; + } - // Determine the expected number of neutrons produced - double nu_t = p.wgt() / simulation::keff * weight * - p.neutron_xs(i_nuclide).nu_fission / - p.neutron_xs(i_nuclide).total; + // Set nu to the number of fission sites successfully stored. If the fission + // bank was not found to be full then these values are already equivalent. + nu = n_sites_stored; - // Sample the number of neutrons produced - int nu = static_cast(nu_t); - if (prn(p.current_seed()) <= (nu_t - nu)) - ++nu; + // Store the total weight banked for analog fission tallies + p.n_bank() = nu; + p.wgt_bank() = nu / weight; + for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { + p.n_delayed_bank(d) = nu_d[d]; + } +} - // If no neutrons were produced then don't continue - if (nu == 0) - return; +void sample_photon_reaction(Particle& p) +{ + // Kill photon if below energy cutoff -- an extra check is made here because + // photons with energy below the cutoff may have been produced by neutrons + // reactions or atomic relaxation + int photon = static_cast(ParticleType::photon); + if (p.E() < settings::energy_cutoff[photon]) { + p.E() = 0.0; + p.wgt() = 0.0; + return; + } - // Initialize the counter of delayed neutrons encountered for each delayed - // group. - double nu_d[MAX_DELAYED_GROUPS] = {0.}; - - // Clear out particle's nu fission bank - p.nu_bank().clear(); - - p.fission() = true; - - // Determine whether to place fission sites into the shared fission bank - // or the secondary particle bank. - bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); - - // Counter for the number of fission sites successfully stored to the shared - // fission bank or the secondary particle bank - int n_sites_stored; - - for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { - // Initialize fission site object with particle data - SourceSite site; - site.r = p.r(); - site.particle = ParticleType::neutron; - site.time = p.time(); - site.wgt = 1. / weight; - site.surf_id = 0; - - // Sample delayed group and angle/energy for fission reaction - sample_fission_neutron(i_nuclide, rx, &site, p); - - // Reject site if it exceeds time cutoff - if (site.delayed_group > 0) { - double t_cutoff = - settings::time_cutoff[static_cast(site.particle)]; - if (site.time > t_cutoff) { - continue; - } - } + // Sample element within material + int i_element = sample_element(p); + const auto& micro {p.photon_xs(i_element)}; + const auto& element {*data::elements[i_element]}; - // Set parent and progeny IDs - site.parent_id = p.id(); - site.progeny_id = p.n_progeny()++; - - // Store fission site in bank - if (use_fission_bank) { - int64_t idx = simulation::fission_bank.thread_safe_append(site); - if (idx == -1) { - warning( - "The shared fission bank is full. Additional fission sites created " - "in this generation will not be banked. Results may be " - "non-deterministic."); - - // Decrement number of particle progeny as storage was unsuccessful. - // This step is needed so that the sum of all progeny is equal to the - // size of the shared fission bank. - p.n_progeny()--; - - // 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, idx); - } - } else { - p.secondary_bank().push_back(site); - } + // Calculate photon energy over electron rest mass equivalent + double alpha = p.E() / MASS_ELECTRON_EV; - // Increment the number of neutrons born delayed - if (site.delayed_group > 0) { - nu_d[site.delayed_group - 1]++; - } + // For tallying purposes, this routine might be called directly. In that + // case, we need to sample a reaction via the cutoff variable + double prob = 0.0; + double cutoff = prn(p.current_seed()) * micro.total; - // Write fission particles to nuBank - NuBank& nu_bank_entry = p.nu_bank().emplace_back(); - nu_bank_entry.wgt = site.wgt; - nu_bank_entry.E = site.E; - nu_bank_entry.delayed_group = site.delayed_group; - } + // Coherent (Rayleigh) scattering + prob += micro.coherent; + if (prob > cutoff) { + p.mu() = element.rayleigh_scatter(alpha, p.current_seed()); + p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); + p.event() = TallyEvent::SCATTER; + p.event_mt() = COHERENT; + return; + } - // If shared fission bank was full, and no fissions could be added, - // set the particle fission flag to false. - if (n_sites_stored == 0) { - p.fission() = false; - return; + // Incoherent (Compton) scattering + prob += micro.incoherent; + if (prob > cutoff) { + double alpha_out; + int i_shell; + element.compton_scatter( + alpha, true, &alpha_out, &p.mu(), &i_shell, p.current_seed()); + + // Determine binding energy of shell. The binding energy is 0.0 if + // doppler broadening is not used. + double e_b; + if (i_shell == -1) { + e_b = 0.0; + } else { + e_b = element.binding_energy_[i_shell]; } - // Set nu to the number of fission sites successfully stored. If the fission - // bank was not found to be full then these values are already equivalent. - nu = n_sites_stored; - - // Store the total weight banked for analog fission tallies - p.n_bank() = nu; - p.wgt_bank() = nu / weight; - for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { - p.n_delayed_bank(d) = nu_d[d]; + // Create Compton electron + double phi = uniform_distribution(0., 2.0 * PI, p.current_seed()); + double E_electron = (alpha - alpha_out) * MASS_ELECTRON_EV - e_b; + int electron = static_cast(ParticleType::electron); + if (E_electron >= settings::energy_cutoff[electron]) { + double mu_electron = (alpha - alpha_out * p.mu()) / + std::sqrt(alpha * alpha + alpha_out * alpha_out - + 2.0 * alpha * alpha_out * p.mu()); + Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed()); + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); } - } - void sample_photon_reaction(Particle & p) - { - // Kill photon if below energy cutoff -- an extra check is made here because - // photons with energy below the cutoff may have been produced by neutrons - // reactions or atomic relaxation - int photon = static_cast(ParticleType::photon); - if (p.E() < settings::energy_cutoff[photon]) { - p.E() = 0.0; - p.wgt() = 0.0; - return; + // Allow electrons to fill orbital and produce Auger electrons and + // fluorescent photons. Since Compton subshell data does not match atomic + // relaxation data, use the mapping between the data to find the subshell + if (i_shell >= 0 && element.subshell_map_[i_shell] >= 0) { + element.atomic_relaxation(element.subshell_map_[i_shell], p); } - // Sample element within material - int i_element = sample_element(p); - const auto& micro {p.photon_xs(i_element)}; - const auto& element {*data::elements[i_element]}; - - // Calculate photon energy over electron rest mass equivalent - double alpha = p.E() / MASS_ELECTRON_EV; - - // For tallying purposes, this routine might be called directly. In that - // case, we need to sample a reaction via the cutoff variable - double prob = 0.0; - double cutoff = prn(p.current_seed()) * micro.total; + phi += PI; + p.E() = alpha_out * MASS_ELECTRON_EV; + p.u() = rotate_angle(p.u(), p.mu(), &phi, p.current_seed()); + p.event() = TallyEvent::SCATTER; + p.event_mt() = INCOHERENT; + return; + } - // Coherent (Rayleigh) scattering - prob += micro.coherent; - if (prob > cutoff) { - p.mu() = element.rayleigh_scatter(alpha, p.current_seed()); - p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); - p.event() = TallyEvent::SCATTER; - p.event_mt() = COHERENT; - return; - } + // Photoelectric effect + double prob_after = prob + micro.photoelectric; - // Incoherent (Compton) scattering - prob += micro.incoherent; - if (prob > cutoff) { - double alpha_out; - int i_shell; - element.compton_scatter( - alpha, true, &alpha_out, &p.mu(), &i_shell, p.current_seed()); - - // Determine binding energy of shell. The binding energy is 0.0 if - // doppler broadening is not used. - double e_b; - if (i_shell == -1) { - e_b = 0.0; - } else { - e_b = element.binding_energy_[i_shell]; - } + if (prob_after > cutoff) { + // Get grid index, interpolation factor, and bounding subshell + // cross sections + int i_grid = micro.index_grid; + double f = micro.interp_factor; + const auto& xs_lower = xt::row(element.cross_sections_, i_grid); + const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1); - // Create Compton electron - double phi = uniform_distribution(0., 2.0 * PI, p.current_seed()); - double E_electron = (alpha - alpha_out) * MASS_ELECTRON_EV - e_b; - int electron = static_cast(ParticleType::electron); - if (E_electron >= settings::energy_cutoff[electron]) { - double mu_electron = (alpha - alpha_out * p.mu()) / - std::sqrt(alpha * alpha + alpha_out * alpha_out - - 2.0 * alpha * alpha_out * p.mu()); - Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed()); - p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); - } + for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) { + const auto& shell {element.shells_[i_shell]}; - // Allow electrons to fill orbital and produce Auger electrons and - // fluorescent photons. Since Compton subshell data does not match atomic - // relaxation data, use the mapping between the data to find the subshell - if (i_shell >= 0 && element.subshell_map_[i_shell] >= 0) { - element.atomic_relaxation(element.subshell_map_[i_shell], p); - } + // Check threshold of reaction + if (xs_lower(i_shell) == 0) + continue; - phi += PI; - p.E() = alpha_out * MASS_ELECTRON_EV; - p.u() = rotate_angle(p.u(), p.mu(), &phi, p.current_seed()); - p.event() = TallyEvent::SCATTER; - p.event_mt() = INCOHERENT; - return; - } + // Evaluation subshell photoionization cross section + prob += std::exp( + xs_lower(i_shell) + f * (xs_upper(i_shell) - xs_lower(i_shell))); - // Photoelectric effect - double prob_after = prob + micro.photoelectric; - - if (prob_after > cutoff) { - // Get grid index, interpolation factor, and bounding subshell - // cross sections - int i_grid = micro.index_grid; - double f = micro.interp_factor; - const auto& xs_lower = xt::row(element.cross_sections_, i_grid); - const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1); - - for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) { - const auto& shell {element.shells_[i_shell]}; - - // Check threshold of reaction - if (xs_lower(i_shell) == 0) - continue; - - // Evaluation subshell photoionization cross section - prob += std::exp( - xs_lower(i_shell) + f * (xs_upper(i_shell) - xs_lower(i_shell))); - - if (prob > cutoff) { - // Determine binding energy based on whether atomic relaxation data is - // present (if not, use value from Compton profile data) - double binding_energy = element.has_atomic_relaxation_ - ? shell.binding_energy - : element.binding_energy_[i_shell]; - - // Determine energy of secondary electron - double E_electron = p.E() - binding_energy; - - // Sample mu using non-relativistic Sauter distribution. - // See Eqns 3.19 and 3.20 in "Implementing a photon physics - // model in Serpent 2" by Toni Kaltiaisenaho - double mu; - while (true) { - double r = prn(p.current_seed()); - if (4.0 * (1.0 - r) * r >= prn(p.current_seed())) { - double rel_vel = - std::sqrt(E_electron * (E_electron + 2.0 * MASS_ELECTRON_EV)) / - (E_electron + MASS_ELECTRON_EV); - mu = - (2.0 * r + rel_vel - 1.0) / (2.0 * rel_vel * r - rel_vel + 1.0); - break; - } + if (prob > cutoff) { + // Determine binding energy based on whether atomic relaxation data is + // present (if not, use value from Compton profile data) + double binding_energy = element.has_atomic_relaxation_ + ? shell.binding_energy + : element.binding_energy_[i_shell]; + + // Determine energy of secondary electron + double E_electron = p.E() - binding_energy; + + // Sample mu using non-relativistic Sauter distribution. + // See Eqns 3.19 and 3.20 in "Implementing a photon physics + // model in Serpent 2" by Toni Kaltiaisenaho + double mu; + while (true) { + double r = prn(p.current_seed()); + if (4.0 * (1.0 - r) * r >= prn(p.current_seed())) { + double rel_vel = + std::sqrt(E_electron * (E_electron + 2.0 * MASS_ELECTRON_EV)) / + (E_electron + MASS_ELECTRON_EV); + mu = + (2.0 * r + rel_vel - 1.0) / (2.0 * rel_vel * r - rel_vel + 1.0); + break; } - - double phi = uniform_distribution(0., 2.0 * PI, p.current_seed()); - Direction u; - u.x = mu; - u.y = std::sqrt(1.0 - mu * mu) * std::cos(phi); - u.z = std::sqrt(1.0 - mu * mu) * std::sin(phi); - - // Create secondary electron - p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); - - // Allow electrons to fill orbital and produce auger electrons - // and fluorescent photons - element.atomic_relaxation(i_shell, p); - p.event() = TallyEvent::ABSORB; - p.event_mt() = 533 + shell.index_subshell; - p.wgt() = 0.0; - p.E() = 0.0; - return; } - } - } - prob = prob_after; - // Pair production - prob += micro.pair_production; - if (prob > cutoff) { - double E_electron, E_positron; - double mu_electron, mu_positron; - element.pair_production(alpha, &E_electron, &E_positron, &mu_electron, - &mu_positron, p.current_seed()); - - // Create secondary electron - Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed()); - p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); + double phi = uniform_distribution(0., 2.0 * PI, p.current_seed()); + Direction u; + u.x = mu; + u.y = std::sqrt(1.0 - mu * mu) * std::cos(phi); + u.z = std::sqrt(1.0 - mu * mu) * std::sin(phi); - // Create secondary positron - u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed()); - p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron); + // Create secondary electron + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); - p.event() = TallyEvent::ABSORB; - p.event_mt() = PAIR_PROD; - p.wgt() = 0.0; - p.E() = 0.0; + // Allow electrons to fill orbital and produce auger electrons + // and fluorescent photons + element.atomic_relaxation(i_shell, p); + p.event() = TallyEvent::ABSORB; + p.event_mt() = 533 + shell.index_subshell; + p.wgt() = 0.0; + p.E() = 0.0; + return; + } } } + prob = prob_after; - void sample_electron_reaction(Particle & p) - { - // TODO: create reaction types + // Pair production + prob += micro.pair_production; + if (prob > cutoff) { + double E_electron, E_positron; + double mu_electron, mu_positron; + element.pair_production(alpha, &E_electron, &E_positron, &mu_electron, + &mu_positron, p.current_seed()); - if (settings::electron_treatment == ElectronTreatment::TTB) { - double E_lost; - thick_target_bremsstrahlung(p, &E_lost); - } + // Create secondary electron + Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed()); + p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron); + + // Create secondary positron + u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed()); + p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron); - p.E() = 0.0; - p.wgt() = 0.0; p.event() = TallyEvent::ABSORB; + p.event_mt() = PAIR_PROD; + p.wgt() = 0.0; + p.E() = 0.0; } +} - void sample_positron_reaction(Particle & p) - { - // TODO: create reaction types +void sample_electron_reaction(Particle& p) +{ + // TODO: create reaction types - if (settings::electron_treatment == ElectronTreatment::TTB) { - double E_lost; - thick_target_bremsstrahlung(p, &E_lost); - } + if (settings::electron_treatment == ElectronTreatment::TTB) { + double E_lost; + thick_target_bremsstrahlung(p, &E_lost); + } - // Sample angle isotropically - Direction u = isotropic_direction(p.current_seed()); + p.E() = 0.0; + p.wgt() = 0.0; + p.event() = TallyEvent::ABSORB; +} - // Create annihilation photon pair traveling in opposite directions - p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon); - p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon); +void sample_positron_reaction(Particle& p) +{ + // TODO: create reaction types - p.E() = 0.0; - p.wgt() = 0.0; - p.event() = TallyEvent::ABSORB; + if (settings::electron_treatment == ElectronTreatment::TTB) { + double E_lost; + thick_target_bremsstrahlung(p, &E_lost); } - int sample_nuclide(Particle & p) - { - // Sample cumulative distribution function - double cutoff = prn(p.current_seed()) * p.macro_xs().total; + // Sample angle isotropically + Direction u = isotropic_direction(p.current_seed()); - // Get pointers to nuclide/density arrays - const auto& mat {model::materials[p.material()]}; - int n = mat->nuclide_.size(); + // Create annihilation photon pair traveling in opposite directions + p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon); + p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon); - double prob = 0.0; - for (int i = 0; i < n; ++i) { - // Get atom density - int i_nuclide = mat->nuclide_[i]; - double atom_density = mat->atom_density(i, p.density_mult()); - - // Increment probability to compare to cutoff - prob += atom_density * p.neutron_xs(i_nuclide).total; - if (prob >= cutoff) - return i_nuclide; - } + p.E() = 0.0; + p.wgt() = 0.0; + p.event() = TallyEvent::ABSORB; +} - // If we reach here, no nuclide was sampled - p.write_restart(); - throw std::runtime_error {"Did not sample any nuclide during collision."}; +int sample_nuclide(Particle& p) +{ + // Sample cumulative distribution function + double cutoff = prn(p.current_seed()) * p.macro_xs().total; + + // Get pointers to nuclide/density arrays + const auto& mat {model::materials[p.material()]}; + int n = mat->nuclide_.size(); + + double prob = 0.0; + for (int i = 0; i < n; ++i) { + // Get atom density + int i_nuclide = mat->nuclide_[i]; + double atom_density = mat->atom_density(i, p.density_mult()); + + // Increment probability to compare to cutoff + prob += atom_density * p.neutron_xs(i_nuclide).total; + if (prob >= cutoff) + return i_nuclide; } - int sample_element(Particle & p) - { - // Sample cumulative distribution function - double cutoff = prn(p.current_seed()) * p.macro_xs().total; + // If we reach here, no nuclide was sampled + p.write_restart(); + throw std::runtime_error {"Did not sample any nuclide during collision."}; +} - // Get pointers to elements, densities - const auto& mat {model::materials[p.material()]}; +int sample_element(Particle& p) +{ + // Sample cumulative distribution function + double cutoff = prn(p.current_seed()) * p.macro_xs().total; - double prob = 0.0; - for (int i = 0; i < mat->element_.size(); ++i) { - // Find atom density - int i_element = mat->element_[i]; - double atom_density = mat->atom_density(i, p.density_mult()); + // Get pointers to elements, densities + const auto& mat {model::materials[p.material()]}; - // Determine microscopic cross section - double sigma = atom_density * p.photon_xs(i_element).total; + double prob = 0.0; + for (int i = 0; i < mat->element_.size(); ++i) { + // Find atom density + int i_element = mat->element_[i]; + double atom_density = mat->atom_density(i, p.density_mult()); - // Increment probability to compare to cutoff - prob += sigma; - if (prob > cutoff) { - // Save which nuclide particle had collision with for tally purpose - p.event_nuclide() = mat->nuclide_[i]; + // Determine microscopic cross section + double sigma = atom_density * p.photon_xs(i_element).total; - return i_element; - } - } + // Increment probability to compare to cutoff + prob += sigma; + if (prob > cutoff) { + // Save which nuclide particle had collision with for tally purpose + p.event_nuclide() = mat->nuclide_[i]; - // If we made it here, no element was sampled - p.write_restart(); - fatal_error("Did not sample any element during collision."); + return i_element; + } } - Reaction& sample_fission(int i_nuclide, Particle& p) - { - // Get pointer to nuclide - const auto& nuc {data::nuclides[i_nuclide]}; + // If we made it here, no element was sampled + p.write_restart(); + fatal_error("Did not sample any element during collision."); +} - // If we're in the URR, by default use the first fission reaction. We also - // default to the first reaction if we know that there are no partial - // fission reactions - if (p.neutron_xs(i_nuclide).use_ptable || !nuc->has_partial_fission_) { - return *nuc->fission_rx_[0]; - } +Reaction& sample_fission(int i_nuclide, Particle& p) +{ + // Get pointer to nuclide + const auto& nuc {data::nuclides[i_nuclide]}; + + // If we're in the URR, by default use the first fission reaction. We also + // default to the first reaction if we know that there are no partial + // fission reactions + if (p.neutron_xs(i_nuclide).use_ptable || !nuc->has_partial_fission_) { + return *nuc->fission_rx_[0]; + } - // Check to see if we are in a windowed multipole range. WMP only supports - // the first fission reaction. - if (nuc->multipole_) { - if (p.E() >= nuc->multipole_->E_min_ && - p.E() <= nuc->multipole_->E_max_) { - return *nuc->fission_rx_[0]; - } + // Check to see if we are in a windowed multipole range. WMP only supports + // the first fission reaction. + if (nuc->multipole_) { + if (p.E() >= nuc->multipole_->E_min_ && p.E() <= nuc->multipole_->E_max_) { + return *nuc->fission_rx_[0]; } + } - // Get grid index and interpolation factor and sample fission cdf - const auto& micro = p.neutron_xs(i_nuclide); - double cutoff = prn(p.current_seed()) * p.neutron_xs(i_nuclide).fission; - double prob = 0.0; - - // Loop through each partial fission reaction type - for (auto& rx : nuc->fission_rx_) { - // add to cumulative probability - prob += rx->xs(micro); + // Get grid index and interpolation factor and sample fission cdf + const auto& micro = p.neutron_xs(i_nuclide); + double cutoff = prn(p.current_seed()) * p.neutron_xs(i_nuclide).fission; + double prob = 0.0; - // Create fission bank sites if fission occurs - if (prob > cutoff) - return *rx; - } + // Loop through each partial fission reaction type + for (auto& rx : nuc->fission_rx_) { + // add to cumulative probability + prob += rx->xs(micro); - // If we reached here, no reaction was sampled - throw std::runtime_error { - "No fission reaction was sampled for " + nuc->name_}; + // Create fission bank sites if fission occurs + if (prob > cutoff) + return *rx; } - void sample_photon_product( - int i_nuclide, Particle& p, int* i_rx, int* i_product) - { - // Get grid index and interpolation factor and sample photon production cdf - const auto& micro = p.neutron_xs(i_nuclide); - double cutoff = prn(p.current_seed()) * micro.photon_prod; - double prob = 0.0; - - // Loop through each reaction type - const auto& nuc {data::nuclides[i_nuclide]}; - for (int i = 0; i < nuc->reactions_.size(); ++i) { - // Evaluate neutron cross section - const auto& rx = nuc->reactions_[i]; - double xs = rx->xs(micro); - - // if cross section is zero for this reaction, skip it - if (xs == 0.0) - continue; + // If we reached here, no reaction was sampled + throw std::runtime_error { + "No fission reaction was sampled for " + nuc->name_}; +} - for (int j = 0; j < rx->products_.size(); ++j) { - if (rx->products_[j].particle_ == ParticleType::photon) { - // For fission, artificially increase the photon yield to account - // for delayed photons - double f = 1.0; - if (settings::delayed_photon_scaling) { - if (is_fission(rx->mt_)) { - if (nuc->prompt_photons_ && nuc->delayed_photons_) { - double energy_prompt = (*nuc->prompt_photons_)(p.E()); - double energy_delayed = (*nuc->delayed_photons_)(p.E()); - f = (energy_prompt + energy_delayed) / (energy_prompt); - } +void sample_photon_product( + int i_nuclide, Particle& p, int* i_rx, int* i_product) +{ + // Get grid index and interpolation factor and sample photon production cdf + const auto& micro = p.neutron_xs(i_nuclide); + double cutoff = prn(p.current_seed()) * micro.photon_prod; + double prob = 0.0; + + // Loop through each reaction type + const auto& nuc {data::nuclides[i_nuclide]}; + for (int i = 0; i < nuc->reactions_.size(); ++i) { + // Evaluate neutron cross section + const auto& rx = nuc->reactions_[i]; + double xs = rx->xs(micro); + + // if cross section is zero for this reaction, skip it + if (xs == 0.0) + continue; + + for (int j = 0; j < rx->products_.size(); ++j) { + if (rx->products_[j].particle_ == ParticleType::photon) { + // For fission, artificially increase the photon yield to account + // for delayed photons + double f = 1.0; + if (settings::delayed_photon_scaling) { + if (is_fission(rx->mt_)) { + if (nuc->prompt_photons_ && nuc->delayed_photons_) { + double energy_prompt = (*nuc->prompt_photons_)(p.E()); + double energy_delayed = (*nuc->delayed_photons_)(p.E()); + f = (energy_prompt + energy_delayed) / (energy_prompt); } } + } - // add to cumulative probability - prob += f * (*rx->products_[j].yield_)(p.E()) * xs; + // add to cumulative probability + prob += f * (*rx->products_[j].yield_)(p.E()) * xs; - *i_rx = i; - *i_product = j; - if (prob > cutoff) - return; - } + *i_rx = i; + *i_product = j; + if (prob > cutoff) + return; } } } +} - void absorption(Particle & p, int i_nuclide) - { - if (settings::survival_biasing) { - // Determine weight absorbed in survival biasing - const double wgt_absorb = p.wgt() * p.neutron_xs(i_nuclide).absorption / - p.neutron_xs(i_nuclide).total; - - // Adjust weight of particle by probability of absorption - p.wgt() -= wgt_absorb; - - // Score implicit absorption estimate of keff +void absorption(Particle& p, int i_nuclide) +{ + if (settings::survival_biasing) { + // Determine weight absorbed in survival biasing + const double wgt_absorb = p.wgt() * p.neutron_xs(i_nuclide).absorption / + p.neutron_xs(i_nuclide).total; + + // Adjust weight of particle by probability of absorption + p.wgt() -= wgt_absorb; + + // Score implicit absorption estimate of keff + if (settings::run_mode == RunMode::EIGENVALUE) { + p.keff_tally_absorption() += wgt_absorb * + p.neutron_xs(i_nuclide).nu_fission / + p.neutron_xs(i_nuclide).absorption; + } + } else { + // See if disappearance reaction happens + if (p.neutron_xs(i_nuclide).absorption > + prn(p.current_seed()) * p.neutron_xs(i_nuclide).total) { + // Score absorption estimate of keff if (settings::run_mode == RunMode::EIGENVALUE) { - p.keff_tally_absorption() += wgt_absorb * + p.keff_tally_absorption() += p.wgt() * p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).absorption; } - } else { - // See if disappearance reaction happens - if (p.neutron_xs(i_nuclide).absorption > - prn(p.current_seed()) * p.neutron_xs(i_nuclide).total) { - // Score absorption estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE) { - p.keff_tally_absorption() += p.wgt() * - p.neutron_xs(i_nuclide).nu_fission / - p.neutron_xs(i_nuclide).absorption; - } - p.wgt() = 0.0; - p.event() = TallyEvent::ABSORB; - if (!p.fission()) { - p.event_mt() = N_DISAPPEAR; - } + p.wgt() = 0.0; + p.event() = TallyEvent::ABSORB; + if (!p.fission()) { + p.event_mt() = N_DISAPPEAR; } } } +} - void scatter(Particle & p, int i_nuclide) - { - // copy incoming direction - Direction u_old {p.u()}; - - // Get pointer to nuclide and grid index/interpolation factor - const auto& nuc {data::nuclides[i_nuclide]}; - const auto& micro {p.neutron_xs(i_nuclide)}; - int i_temp = micro.index_temp; - - // For tallying purposes, this routine might be called directly. In that - // case, we need to sample a reaction via the cutoff variable - double cutoff = prn(p.current_seed()) * (micro.total - micro.absorption); - bool sampled = false; - - // Calculate elastic cross section if it wasn't precalculated - if (micro.elastic == CACHE_INVALID) { - nuc->calculate_elastic_xs(p); - } - - double prob = micro.elastic - micro.thermal; - if (prob > cutoff) { - // ======================================================================= - // NON-S(A,B) ELASTIC SCATTERING +void scatter(Particle& p, int i_nuclide) +{ + // copy incoming direction + Direction u_old {p.u()}; + + // Get pointer to nuclide and grid index/interpolation factor + const auto& nuc {data::nuclides[i_nuclide]}; + const auto& micro {p.neutron_xs(i_nuclide)}; + int i_temp = micro.index_temp; + + // For tallying purposes, this routine might be called directly. In that + // case, we need to sample a reaction via the cutoff variable + double cutoff = prn(p.current_seed()) * (micro.total - micro.absorption); + bool sampled = false; + + // Calculate elastic cross section if it wasn't precalculated + if (micro.elastic == CACHE_INVALID) { + nuc->calculate_elastic_xs(p); + } - // Determine temperature - double kT = nuc->multipole_ ? p.sqrtkT() * p.sqrtkT() : nuc->kTs_[i_temp]; + double prob = micro.elastic - micro.thermal; + if (prob > cutoff) { + // ======================================================================= + // NON-S(A,B) ELASTIC SCATTERING - // Perform collision physics for elastic scattering - elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p); + // Determine temperature + double kT = nuc->multipole_ ? p.sqrtkT() * p.sqrtkT() : nuc->kTs_[i_temp]; - p.event_mt() = ELASTIC; - sampled = true; - } + // Perform collision physics for elastic scattering + elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p); - prob = micro.elastic; - if (prob > cutoff && !sampled) { - // ======================================================================= - // S(A,B) SCATTERING + p.event_mt() = ELASTIC; + sampled = true; + } - sab_scatter(i_nuclide, micro.index_sab, p); + prob = micro.elastic; + if (prob > cutoff && !sampled) { + // ======================================================================= + // S(A,B) SCATTERING - p.event_mt() = ELASTIC; - sampled = true; - } + sab_scatter(i_nuclide, micro.index_sab, p); - if (!sampled) { - // ======================================================================= - // INELASTIC SCATTERING + p.event_mt() = ELASTIC; + sampled = true; + } - int n = nuc->index_inelastic_scatter_.size(); - int i = 0; - for (int j = 0; j < n && prob < cutoff; ++j) { - i = nuc->index_inelastic_scatter_[j]; + if (!sampled) { + // ======================================================================= + // INELASTIC SCATTERING - // add to cumulative probability - prob += nuc->reactions_[i]->xs(micro); - } + int n = nuc->index_inelastic_scatter_.size(); + int i = 0; + for (int j = 0; j < n && prob < cutoff; ++j) { + i = nuc->index_inelastic_scatter_[j]; - // Perform collision physics for inelastic scattering - const auto& rx {nuc->reactions_[i]}; - inelastic_scatter(*nuc, *rx, p); - p.event_mt() = rx->mt_; + // add to cumulative probability + prob += nuc->reactions_[i]->xs(micro); } - // Set event component - p.event() = TallyEvent::SCATTER; + // Perform collision physics for inelastic scattering + const auto& rx {nuc->reactions_[i]}; + inelastic_scatter(*nuc, *rx, p); + p.event_mt() = rx->mt_; + } - // Sample new outgoing angle for isotropic-in-lab scattering - const auto& mat {model::materials[p.material()]}; - if (!mat->p0_.empty()) { - int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide]; - if (mat->p0_[i_nuc_mat]) { - // Sample isotropic-in-lab outgoing direction - p.u() = isotropic_direction(p.current_seed()); - p.mu() = u_old.dot(p.u()); - } + // Set event component + p.event() = TallyEvent::SCATTER; + + // Sample new outgoing angle for isotropic-in-lab scattering + const auto& mat {model::materials[p.material()]}; + if (!mat->p0_.empty()) { + int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide]; + if (mat->p0_[i_nuc_mat]) { + // Sample isotropic-in-lab outgoing direction + p.u() = isotropic_direction(p.current_seed()); + p.mu() = u_old.dot(p.u()); } } +} - void elastic_scatter( - int i_nuclide, const Reaction& rx, double kT, Particle& p) - { - // get pointer to nuclide - const auto& nuc {data::nuclides[i_nuclide]}; - - double vel = std::sqrt(p.E()); - double awr = nuc->awr_; +void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p) +{ + // get pointer to nuclide + const auto& nuc {data::nuclides[i_nuclide]}; - // Neutron velocity in LAB - Direction v_n = vel * p.u(); + double vel = std::sqrt(p.E()); + double awr = nuc->awr_; - // Sample velocity of target nucleus - Direction v_t {}; - if (!p.neutron_xs(i_nuclide).use_ptable) { - v_t = sample_target_velocity(*nuc, p.E(), p.u(), v_n, - p.neutron_xs(i_nuclide).elastic, kT, p.current_seed()); - } + // Neutron velocity in LAB + Direction v_n = vel * p.u(); - // Velocity of center-of-mass - Direction v_cm = (v_n + awr * v_t) / (awr + 1.0); + // Sample velocity of target nucleus + Direction v_t {}; + if (!p.neutron_xs(i_nuclide).use_ptable) { + v_t = sample_target_velocity(*nuc, p.E(), p.u(), v_n, + p.neutron_xs(i_nuclide).elastic, kT, p.current_seed()); + } - // Transform to CM frame - v_n -= v_cm; + // Velocity of center-of-mass + Direction v_cm = (v_n + awr * v_t) / (awr + 1.0); - // Find speed of neutron in CM - vel = v_n.norm(); + // Transform to CM frame + v_n -= v_cm; - // Sample scattering angle, checking if angle distribution is present - // (assume isotropic otherwise) - double mu_cm; - auto& d = rx.products_[0].distribution_[0]; - auto d_ = dynamic_cast(d.get()); - if (!d_->angle().empty()) { - mu_cm = d_->angle().sample(p.E(), p.current_seed()); - } else { - mu_cm = uniform_distribution(-1., 1., p.current_seed()); - } + // Find speed of neutron in CM + vel = v_n.norm(); - // Determine direction cosines in CM - Direction u_cm = v_n / vel; + // Sample scattering angle, checking if angle distribution is present + // (assume isotropic otherwise) + double mu_cm; + auto& d = rx.products_[0].distribution_[0]; + auto d_ = dynamic_cast(d.get()); + if (!d_->angle().empty()) { + mu_cm = d_->angle().sample(p.E(), p.current_seed()); + } else { + mu_cm = uniform_distribution(-1., 1., p.current_seed()); + } - // Rotate neutron velocity vector to new angle -- note that the speed of the - // neutron in CM does not change in elastic scattering. However, the speed - // will change when we convert back to LAB - v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p.current_seed()); + // Determine direction cosines in CM + Direction u_cm = v_n / vel; - // Transform back to LAB frame - v_n += v_cm; + // Rotate neutron velocity vector to new angle -- note that the speed of the + // neutron in CM does not change in elastic scattering. However, the speed + // will change when we convert back to LAB + v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p.current_seed()); - p.E() = v_n.dot(v_n); - vel = std::sqrt(p.E()); + // Transform back to LAB frame + v_n += v_cm; - // compute cosine of scattering angle in LAB frame by taking dot product of - // neutron's pre- and post-collision angle - p.mu() = p.u().dot(v_n) / vel; + p.E() = v_n.dot(v_n); + vel = std::sqrt(p.E()); - // Set energy and direction of particle in LAB frame - p.u() = v_n / vel; + // compute cosine of scattering angle in LAB frame by taking dot product of + // neutron's pre- and post-collision angle + p.mu() = p.u().dot(v_n) / vel; - // Because of floating-point roundoff, it may be possible for mu_lab to be - // outside of the range [-1,1). In these cases, we just set mu_lab to - // exactly -1 or 1 - if (std::abs(p.mu()) > 1.0) - p.mu() = std::copysign(1.0, p.mu()); - } + // Set energy and direction of particle in LAB frame + p.u() = v_n / vel; - void sab_scatter(int i_nuclide, int i_sab, Particle& p) - { - // Determine temperature index - const auto& micro {p.neutron_xs(i_nuclide)}; - int i_temp = micro.index_temp_sab; + // Because of floating-point roundoff, it may be possible for mu_lab to be + // outside of the range [-1,1). In these cases, we just set mu_lab to + // exactly -1 or 1 + if (std::abs(p.mu()) > 1.0) + p.mu() = std::copysign(1.0, p.mu()); +} - // Sample energy and angle - double E_out; - data::thermal_scatt[i_sab]->data_[i_temp].sample( - micro, p.E(), &E_out, &p.mu(), p.current_seed()); +void sab_scatter(int i_nuclide, int i_sab, Particle& p) +{ + // Determine temperature index + const auto& micro {p.neutron_xs(i_nuclide)}; + int i_temp = micro.index_temp_sab; + + // Sample energy and angle + double E_out; + data::thermal_scatt[i_sab]->data_[i_temp].sample( + micro, p.E(), &E_out, &p.mu(), p.current_seed()); + + // Set energy to outgoing, change direction of particle + p.E() = E_out; + p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); +} + +Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, + Direction v_neut, double xs_eff, double kT, uint64_t* seed) +{ + // check if nuclide is a resonant scatterer + ResScatMethod sampling_method; + if (nuc.resonant_) { + + // sampling method to use + sampling_method = settings::res_scat_method; + + // upper resonance scattering energy bound (target is at rest above this + // E) + if (E > settings::res_scat_energy_max) { + return {}; + + // lower resonance scattering energy bound (should be no resonances + // below) + } else if (E < settings::res_scat_energy_min) { + sampling_method = ResScatMethod::cxs; + } - // Set energy to outgoing, change direction of particle - p.E() = E_out; - p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); + // otherwise, use free gas model + } else { + if (E >= settings::free_gas_threshold * kT && nuc.awr_ > 1.0) { + return {}; + } else { + sampling_method = ResScatMethod::cxs; + } } - Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, - Direction v_neut, double xs_eff, double kT, uint64_t* seed) - { - // check if nuclide is a resonant scatterer - ResScatMethod sampling_method; - if (nuc.resonant_) { - - // sampling method to use - sampling_method = settings::res_scat_method; - - // upper resonance scattering energy bound (target is at rest above this - // E) - if (E > settings::res_scat_energy_max) { - return {}; - - // lower resonance scattering energy bound (should be no resonances - // below) - } else if (E < settings::res_scat_energy_min) { - sampling_method = ResScatMethod::cxs; - } - - // otherwise, use free gas model + // use appropriate target velocity sampling method + switch (sampling_method) { + case ResScatMethod::cxs: + + // sample target velocity with the constant cross section (cxs) approx. + return sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); + + case ResScatMethod::dbrc: + case ResScatMethod::rvs: { + double E_red = std::sqrt(nuc.awr_ * E / kT); + double E_low = std::pow(std::max(0.0, E_red - 4.0), 2) * kT / nuc.awr_; + double E_up = (E_red + 4.0) * (E_red + 4.0) * kT / nuc.awr_; + + // find lower and upper energy bound indices + // lower index + int i_E_low; + if (E_low < nuc.energy_0K_.front()) { + i_E_low = 0; + } else if (E_low > nuc.energy_0K_.back()) { + i_E_low = nuc.energy_0K_.size() - 2; } else { - if (E >= settings::free_gas_threshold * kT && nuc.awr_ > 1.0) { - return {}; - } else { - sampling_method = ResScatMethod::cxs; - } + i_E_low = + lower_bound_index(nuc.energy_0K_.begin(), nuc.energy_0K_.end(), E_low); } - // use appropriate target velocity sampling method - switch (sampling_method) { - case ResScatMethod::cxs: + // upper index + int i_E_up; + if (E_up < nuc.energy_0K_.front()) { + i_E_up = 0; + } else if (E_up > nuc.energy_0K_.back()) { + i_E_up = nuc.energy_0K_.size() - 2; + } else { + i_E_up = + lower_bound_index(nuc.energy_0K_.begin(), nuc.energy_0K_.end(), E_up); + } - // sample target velocity with the constant cross section (cxs) approx. + if (i_E_up == i_E_low) { + // Handle degenerate case -- if the upper/lower bounds occur for the + // same index, then using cxs is probably a good approximation return sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); + } - case ResScatMethod::dbrc: - case ResScatMethod::rvs: { - double E_red = std::sqrt(nuc.awr_ * E / kT); - double E_low = std::pow(std::max(0.0, E_red - 4.0), 2) * kT / nuc.awr_; - double E_up = (E_red + 4.0) * (E_red + 4.0) * kT / nuc.awr_; - - // find lower and upper energy bound indices - // lower index - int i_E_low; - if (E_low < nuc.energy_0K_.front()) { - i_E_low = 0; - } else if (E_low > nuc.energy_0K_.back()) { - i_E_low = nuc.energy_0K_.size() - 2; - } else { - i_E_low = lower_bound_index( - nuc.energy_0K_.begin(), nuc.energy_0K_.end(), E_low); - } - - // upper index - int i_E_up; - if (E_up < nuc.energy_0K_.front()) { - i_E_up = 0; - } else if (E_up > nuc.energy_0K_.back()) { - i_E_up = nuc.energy_0K_.size() - 2; - } else { - i_E_up = - lower_bound_index(nuc.energy_0K_.begin(), nuc.energy_0K_.end(), E_up); - } + if (sampling_method == ResScatMethod::dbrc) { + // interpolate xs since we're not exactly at the energy indices + double xs_low = nuc.elastic_0K_[i_E_low]; + double m = (nuc.elastic_0K_[i_E_low + 1] - xs_low) / + (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]); + xs_low += m * (E_low - nuc.energy_0K_[i_E_low]); + double xs_up = nuc.elastic_0K_[i_E_up]; + m = (nuc.elastic_0K_[i_E_up + 1] - xs_up) / + (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]); + xs_up += m * (E_up - nuc.energy_0K_[i_E_up]); + + // get max 0K xs value over range of practical relative energies + double xs_max = *std::max_element( + &nuc.elastic_0K_[i_E_low + 1], &nuc.elastic_0K_[i_E_up + 1]); + xs_max = std::max({xs_low, xs_max, xs_up}); + + while (true) { + double E_rel; + Direction v_target; + while (true) { + // sample target velocity with the constant cross section (cxs) + // approx. + v_target = sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); + Direction v_rel = v_neut - v_target; + E_rel = v_rel.dot(v_rel); + if (E_rel < E_up) + break; + } - if (i_E_up == i_E_low) { - // Handle degenerate case -- if the upper/lower bounds occur for the - // same index, then using cxs is probably a good approximation - return sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); + // perform Doppler broadening rejection correction (dbrc) + double xs_0K = nuc.elastic_xs_0K(E_rel); + double R = xs_0K / xs_max; + if (prn(seed) < R) + return v_target; } - if (sampling_method == ResScatMethod::dbrc) { - // interpolate xs since we're not exactly at the energy indices - double xs_low = nuc.elastic_0K_[i_E_low]; - double m = (nuc.elastic_0K_[i_E_low + 1] - xs_low) / + } else if (sampling_method == ResScatMethod::rvs) { + // interpolate xs CDF since we're not exactly at the energy indices + // cdf value at lower bound attainable energy + double cdf_low = 0.0; + if (E_low > nuc.energy_0K_.front()) { + double m = (nuc.xs_cdf_[i_E_low + 1] - nuc.xs_cdf_[i_E_low]) / (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]); - xs_low += m * (E_low - nuc.energy_0K_[i_E_low]); - double xs_up = nuc.elastic_0K_[i_E_up]; - m = (nuc.elastic_0K_[i_E_up + 1] - xs_up) / - (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]); - xs_up += m * (E_up - nuc.energy_0K_[i_E_up]); - - // get max 0K xs value over range of practical relative energies - double xs_max = *std::max_element( - &nuc.elastic_0K_[i_E_low + 1], &nuc.elastic_0K_[i_E_up + 1]); - xs_max = std::max({xs_low, xs_max, xs_up}); - - while (true) { - double E_rel; - Direction v_target; - while (true) { - // sample target velocity with the constant cross section (cxs) - // approx. - v_target = sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); - Direction v_rel = v_neut - v_target; - E_rel = v_rel.dot(v_rel); - if (E_rel < E_up) - break; - } + cdf_low = nuc.xs_cdf_[i_E_low] + m * (E_low - nuc.energy_0K_[i_E_low]); + } - // perform Doppler broadening rejection correction (dbrc) - double xs_0K = nuc.elastic_xs_0K(E_rel); - double R = xs_0K / xs_max; - if (prn(seed) < R) - return v_target; + // cdf value at upper bound attainable energy + double m = (nuc.xs_cdf_[i_E_up + 1] - nuc.xs_cdf_[i_E_up]) / + (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]); + double cdf_up = nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]); + + while (true) { + // directly sample Maxwellian + double E_t = -kT * std::log(prn(seed)); + + // sample a relative energy using the xs cdf + double cdf_rel = cdf_low + prn(seed) * (cdf_up - cdf_low); + int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low, + nuc.xs_cdf_.begin() + i_E_up + 2, cdf_rel); + double E_rel = nuc.energy_0K_[i_E_low + i_E_rel]; + double m = (nuc.xs_cdf_[i_E_low + i_E_rel + 1] - + nuc.xs_cdf_[i_E_low + i_E_rel]) / + (nuc.energy_0K_[i_E_low + i_E_rel + 1] - + nuc.energy_0K_[i_E_low + i_E_rel]); + E_rel += (cdf_rel - nuc.xs_cdf_[i_E_low + i_E_rel]) / m; + + // perform rejection sampling on cosine between + // neutron and target velocities + double mu = (E_t + nuc.awr_ * (E - E_rel)) / + (2.0 * std::sqrt(nuc.awr_ * E * E_t)); + + if (std::abs(mu) < 1.0) { + // set and accept target velocity + E_t /= nuc.awr_; + return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, seed); } + } + } + } // case RVS, DBRC + } // switch (sampling_method) - } else if (sampling_method == ResScatMethod::rvs) { - // interpolate xs CDF since we're not exactly at the energy indices - // cdf value at lower bound attainable energy - double cdf_low = 0.0; - if (E_low > nuc.energy_0K_.front()) { - double m = (nuc.xs_cdf_[i_E_low + 1] - nuc.xs_cdf_[i_E_low]) / - (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]); - cdf_low = - nuc.xs_cdf_[i_E_low] + m * (E_low - nuc.energy_0K_[i_E_low]); - } + UNREACHABLE(); +} - // cdf value at upper bound attainable energy - double m = (nuc.xs_cdf_[i_E_up + 1] - nuc.xs_cdf_[i_E_up]) / - (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]); - double cdf_up = - nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]); +Direction sample_cxs_target_velocity( + double awr, double E, Direction u, double kT, uint64_t* seed) +{ + double beta_vn = std::sqrt(awr * E / kT); + double alpha = 1.0 / (1.0 + std::sqrt(PI) * beta_vn / 2.0); - while (true) { - // directly sample Maxwellian - double E_t = -kT * std::log(prn(seed)); - - // sample a relative energy using the xs cdf - double cdf_rel = cdf_low + prn(seed) * (cdf_up - cdf_low); - int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low, - nuc.xs_cdf_.begin() + i_E_up + 2, cdf_rel); - double E_rel = nuc.energy_0K_[i_E_low + i_E_rel]; - double m = (nuc.xs_cdf_[i_E_low + i_E_rel + 1] - - nuc.xs_cdf_[i_E_low + i_E_rel]) / - (nuc.energy_0K_[i_E_low + i_E_rel + 1] - - nuc.energy_0K_[i_E_low + i_E_rel]); - E_rel += (cdf_rel - nuc.xs_cdf_[i_E_low + i_E_rel]) / m; - - // perform rejection sampling on cosine between - // neutron and target velocities - double mu = (E_t + nuc.awr_ * (E - E_rel)) / - (2.0 * std::sqrt(nuc.awr_ * E * E_t)); - - if (std::abs(mu) < 1.0) { - // set and accept target velocity - E_t /= nuc.awr_; - return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, seed); - } - } - } - } // case RVS, DBRC - } // switch (sampling_method) + double beta_vt_sq; + double mu; + while (true) { + // Sample two random numbers + double r1 = prn(seed); + double r2 = prn(seed); - UNREACHABLE(); - } + if (prn(seed) < alpha) { + // With probability alpha, we sample the distribution p(y) = + // y*e^(-y). This can be done with sampling scheme C45 from the Monte + // Carlo sampler - Direction sample_cxs_target_velocity( - double awr, double E, Direction u, double kT, uint64_t* seed) - { - double beta_vn = std::sqrt(awr * E / kT); - double alpha = 1.0 / (1.0 + std::sqrt(PI) * beta_vn / 2.0); + beta_vt_sq = -std::log(r1 * r2); - double beta_vt_sq; - double mu; - while (true) { - // Sample two random numbers - double r1 = prn(seed); - double r2 = prn(seed); + } else { + // With probability 1-alpha, we sample the distribution p(y) = y^2 * + // e^(-y^2). This can be done with sampling scheme C61 from the Monte + // Carlo sampler - if (prn(seed) < alpha) { - // With probability alpha, we sample the distribution p(y) = - // y*e^(-y). This can be done with sampling scheme C45 from the Monte - // Carlo sampler + double c = std::cos(PI / 2.0 * prn(seed)); + beta_vt_sq = -std::log(r1) - std::log(r2) * c * c; + } - beta_vt_sq = -std::log(r1 * r2); + // Determine beta * vt + double beta_vt = std::sqrt(beta_vt_sq); - } else { - // With probability 1-alpha, we sample the distribution p(y) = y^2 * - // e^(-y^2). This can be done with sampling scheme C61 from the Monte - // Carlo sampler + // Sample cosine of angle between neutron and target velocity + mu = uniform_distribution(-1., 1., seed); - double c = std::cos(PI / 2.0 * prn(seed)); - beta_vt_sq = -std::log(r1) - std::log(r2) * c * c; - } + // Determine rejection probability + double accept_prob = + std::sqrt(beta_vn * beta_vn + beta_vt_sq - 2 * beta_vn * beta_vt * mu) / + (beta_vn + beta_vt); - // Determine beta * vt - double beta_vt = std::sqrt(beta_vt_sq); + // Perform rejection sampling on vt and mu + if (prn(seed) < accept_prob) + break; + } - // Sample cosine of angle between neutron and target velocity - mu = uniform_distribution(-1., 1., seed); + // Determine speed of target nucleus + double vt = std::sqrt(beta_vt_sq * kT / awr); - // Determine rejection probability - double accept_prob = - std::sqrt(beta_vn * beta_vn + beta_vt_sq - 2 * beta_vn * beta_vt * mu) / - (beta_vn + beta_vt); + // Determine velocity vector of target nucleus based on neutron's velocity + // and the sampled angle between them + return vt * rotate_angle(u, mu, nullptr, seed); +} - // Perform rejection sampling on vt and mu - if (prn(seed) < accept_prob) +void sample_fission_neutron( + int i_nuclide, const Reaction& rx, SourceSite* site, Particle& p) +{ + // Get attributes of particle + double E_in = p.E(); + uint64_t* seed = p.current_seed(); + + // Determine total nu, delayed nu, and delayed neutron fraction + const auto& nuc {data::nuclides[i_nuclide]}; + double nu_t = nuc->nu(E_in, Nuclide::EmissionMode::total); + double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed); + double beta = nu_d / nu_t; + + if (prn(seed) < beta) { + // ==================================================================== + // DELAYED NEUTRON SAMPLED + + // sampled delayed precursor group + double xi = prn(seed) * nu_d; + double prob = 0.0; + int group; + for (group = 1; group < nuc->n_precursor_; ++group) { + // determine delayed neutron precursor yield for group j + double yield = (*rx.products_[group].yield_)(E_in); + + // Check if this group is sampled + prob += yield; + if (xi < prob) break; } - // Determine speed of target nucleus - double vt = std::sqrt(beta_vt_sq * kT / awr); + // if the sum of the probabilities is slightly less than one and the + // random number is greater, j will be greater than nuc % + // n_precursor -- check for this condition + group = std::min(group, nuc->n_precursor_); - // Determine velocity vector of target nucleus based on neutron's velocity - // and the sampled angle between them - return vt * rotate_angle(u, mu, nullptr, seed); - } + // set the delayed group for the particle born from fission + site->delayed_group = group; - void sample_fission_neutron( - int i_nuclide, const Reaction& rx, SourceSite* site, Particle& p) - { - // Get attributes of particle - double E_in = p.E(); - uint64_t* seed = p.current_seed(); - - // Determine total nu, delayed nu, and delayed neutron fraction - const auto& nuc {data::nuclides[i_nuclide]}; - double nu_t = nuc->nu(E_in, Nuclide::EmissionMode::total); - double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed); - double beta = nu_d / nu_t; - - if (prn(seed) < beta) { - // ==================================================================== - // DELAYED NEUTRON SAMPLED - - // sampled delayed precursor group - double xi = prn(seed) * nu_d; - double prob = 0.0; - int group; - for (group = 1; group < nuc->n_precursor_; ++group) { - // determine delayed neutron precursor yield for group j - double yield = (*rx.products_[group].yield_)(E_in); - - // Check if this group is sampled - prob += yield; - if (xi < prob) - break; - } - - // if the sum of the probabilities is slightly less than one and the - // random number is greater, j will be greater than nuc % - // n_precursor -- check for this condition - group = std::min(group, nuc->n_precursor_); + // Sample time of emission based on decay constant of precursor + double decay_rate = rx.products_[site->delayed_group].decay_rate_; + site->time -= std::log(prn(p.current_seed())) / decay_rate; - // set the delayed group for the particle born from fission - site->delayed_group = group; + } else { + // ==================================================================== + // PROMPT NEUTRON SAMPLED - // Sample time of emission based on decay constant of precursor - double decay_rate = rx.products_[site->delayed_group].decay_rate_; - site->time -= std::log(prn(p.current_seed())) / decay_rate; - - } else { - // ==================================================================== - // PROMPT NEUTRON SAMPLED + // set the delayed group for the particle born from fission to 0 + site->delayed_group = 0; + } - // set the delayed group for the particle born from fission to 0 - site->delayed_group = 0; + // sample from prompt neutron energy distribution + int n_sample = 0; + double mu; + while (true) { + rx.products_[site->delayed_group].sample(E_in, site->E, mu, seed); + + // resample if energy is greater than maximum neutron energy + constexpr int neutron = static_cast(ParticleType::neutron); + if (site->E < data::energy_max[neutron]) + break; + + // check for large number of resamples + ++n_sample; + if (n_sample == MAX_SAMPLE) { + // particle_write_restart(p) + fatal_error("Resampled energy distribution maximum number of times " + "for nuclide " + + nuc->name_); } + } - // sample from prompt neutron energy distribution - int n_sample = 0; - double mu; - while (true) { - rx.products_[site->delayed_group].sample(E_in, site->E, mu, seed); + // Sample azimuthal angle uniformly in [0, 2*pi) and assign angle + site->u = rotate_angle(p.u(), mu, nullptr, seed); +} - // resample if energy is greater than maximum neutron energy - constexpr int neutron = static_cast(ParticleType::neutron); - if (site->E < data::energy_max[neutron]) - break; +void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) +{ + // copy energy of neutron + double E_in = p.E(); + + // sample outgoing energy and scattering cosine + double E; + double mu; + rx.products_[0].sample(E_in, E, mu, p.current_seed()); + + // if scattering system is in center-of-mass, transfer cosine of scattering + // angle and outgoing energy from CM to LAB + if (rx.scatter_in_cm_) { + double E_cm = E; + + // determine outgoing energy in lab + double A = nuc.awr_; + E = E_cm + (E_in + 2.0 * mu * (A + 1.0) * std::sqrt(E_in * E_cm)) / + ((A + 1.0) * (A + 1.0)); + + // determine outgoing angle in lab + mu = mu * std::sqrt(E_cm / E) + 1.0 / (A + 1.0) * std::sqrt(E_in / E); + } - // check for large number of resamples - ++n_sample; - if (n_sample == MAX_SAMPLE) { - // particle_write_restart(p) - fatal_error("Resampled energy distribution maximum number of times " - "for nuclide " + - nuc->name_); - } + // Because of floating-point roundoff, it may be possible for mu to be + // outside of the range [-1,1). In these cases, we just set mu to exactly -1 + // or 1 + if (std::abs(mu) > 1.0) + mu = std::copysign(1.0, mu); + + // Set outgoing energy and scattering angle + p.E() = E; + p.mu() = mu; + + // change direction of particle + p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed()); + + // evaluate yield + double yield = (*rx.products_[0].yield_)(E_in); + if (std::floor(yield) == yield && yield > 0) { + // If yield is integral, create exactly that many secondary particles + for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { + p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron); } + } else { + // Otherwise, change weight of particle based on yield + p.wgt() *= yield; + } +} - // Sample azimuthal angle uniformly in [0, 2*pi) and assign angle - site->u = rotate_angle(p.u(), mu, nullptr, seed); +void sample_secondary_photons(Particle& p, int i_nuclide) +{ + // Sample the number of photons produced + double y_t = + p.neutron_xs(i_nuclide).photon_prod / p.neutron_xs(i_nuclide).total; + double photon_wgt = p.wgt(); + int y = 1; + + if (settings::use_decay_photons) { + // For decay photons, sample a single photon and modify the weight + if (y_t <= 0.0) + return; + photon_wgt *= y_t; + } else { + // For prompt photons, sample an integral number of photons with weight + // equal to the neutron's weight + y = static_cast(y_t); + if (prn(p.current_seed()) <= y_t - y) + ++y; } - void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) - { - // copy energy of neutron - double E_in = p.E(); + // Sample each secondary photon + for (int i = 0; i < y; ++i) { + // Sample the reaction and product + int i_rx; + int i_product; + sample_photon_product(i_nuclide, p, &i_rx, &i_product); - // sample outgoing energy and scattering cosine + // Sample the outgoing energy and angle + auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx]; double E; double mu; - rx.products_[0].sample(E_in, E, mu, p.current_seed()); - - // if scattering system is in center-of-mass, transfer cosine of scattering - // angle and outgoing energy from CM to LAB - if (rx.scatter_in_cm_) { - double E_cm = E; - - // determine outgoing energy in lab - double A = nuc.awr_; - E = E_cm + (E_in + 2.0 * mu * (A + 1.0) * std::sqrt(E_in * E_cm)) / - ((A + 1.0) * (A + 1.0)); - - // determine outgoing angle in lab - mu = mu * std::sqrt(E_cm / E) + 1.0 / (A + 1.0) * std::sqrt(E_in / E); - } - - // Because of floating-point roundoff, it may be possible for mu to be - // outside of the range [-1,1). In these cases, we just set mu to exactly -1 - // or 1 - if (std::abs(mu) > 1.0) - mu = std::copysign(1.0, mu); - - // Set outgoing energy and scattering angle - p.E() = E; - p.mu() = mu; - - // change direction of particle - p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed()); - - // evaluate yield - double yield = (*rx.products_[0].yield_)(E_in); - if (std::floor(yield) == yield && yield > 0) { - // If yield is integral, create exactly that many secondary particles - for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { - p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron); - } - } else { - // Otherwise, change weight of particle based on yield - p.wgt() *= yield; - } - } - - void sample_secondary_photons(Particle & p, int i_nuclide) - { - // Sample the number of photons produced - double y_t = - p.neutron_xs(i_nuclide).photon_prod / p.neutron_xs(i_nuclide).total; - double photon_wgt = p.wgt(); - int y = 1; - - if (settings::use_decay_photons) { - // For decay photons, sample a single photon and modify the weight - if (y_t <= 0.0) - return; - photon_wgt *= y_t; - } else { - // For prompt photons, sample an integral number of photons with weight - // equal to the neutron's weight - y = static_cast(y_t); - if (prn(p.current_seed()) <= y_t - y) - ++y; + rx->products_[i_product].sample(p.E(), E, mu, p.current_seed()); + + // Sample the new direction + Direction u = rotate_angle(p.u(), mu, nullptr, p.current_seed()); + + // In a k-eigenvalue simulation, it's necessary to provide higher weight + // to secondary photons from non-fission reactions to properly balance + // energy release and deposition. See D. P. Griesheimer, S. J. Douglass, + // and M. H. Stedry, "Self-consistent energy normalization for quasistatic + // reactor calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. + double wgt = photon_wgt; + if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) { + wgt *= simulation::keff; } - // Sample each secondary photon - for (int i = 0; i < y; ++i) { - // Sample the reaction and product - int i_rx; - int i_product; - sample_photon_product(i_nuclide, p, &i_rx, &i_product); - - // Sample the outgoing energy and angle - auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx]; - double E; - double mu; - rx->products_[i_product].sample(p.E(), E, mu, p.current_seed()); - - // Sample the new direction - Direction u = rotate_angle(p.u(), mu, nullptr, p.current_seed()); - - // In a k-eigenvalue simulation, it's necessary to provide higher weight - // to secondary photons from non-fission reactions to properly balance - // energy release and deposition. See D. P. Griesheimer, S. J. Douglass, - // and M. H. Stedry, "Self-consistent energy normalization for quasistatic - // reactor calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. - double wgt = photon_wgt; - if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) { - wgt *= simulation::keff; - } - - // Create the secondary photon - bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon); + // Create the secondary photon + bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon); - // Tag secondary particle with parent nuclide - if (created_photon && settings::use_decay_photons) { - p.secondary_bank().back().parent_nuclide = - rx->products_[i_product].parent_nuclide_; - } + // Tag secondary particle with parent nuclide + if (created_photon && settings::use_decay_photons) { + p.secondary_bank().back().parent_nuclide = + rx->products_[i_product].parent_nuclide_; } } +} } // namespace openmc diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index e2e3f687dd4..e2bc0529ae1 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -37,234 +37,230 @@ void collision_mg(Particle& p) auto ww = search_weight_window(p); if (ww.is_valid()) { apply_weight_window(p, ww); - else if (p.type() == ParticleType::neutron) - { - apply_russian_roulette(p) - } - // If collision checkpoint is disabled, apply russian roulette if - // the particle is outside the weight window domain - } else if (settings::weight_windows_on) { - auto ww = search_weight_window(p); - if (!ww.is_valid()) - apply_russian_roulette(p); - } - - // Display information about collision - if ((settings::verbosity >= 10) || p.trace()) { - write_message(fmt::format(" Energy Group = {}", p.g()), 1); + } else if (p.type() == ParticleType::neutron) { + apply_russian_roulette(p) } + // If collision checkpoint is disabled, apply russian roulette if + // the particle is outside the weight window domain + } else if (settings::weight_windows_on) { + auto ww = search_weight_window(p); + if (!ww.is_valid()) + apply_russian_roulette(p); } - void sample_reaction(Particle & p) - { - // Create fission bank sites. Note that while a fission reaction is sampled, - // it never actually "happens", i.e. the weight of the particle does not - // change when sampling fission sites. The following block handles all - // absorption (including fission) - - if (model::materials[p.material()]->fissionable()) { - if (settings::run_mode == RunMode::EIGENVALUE || - (settings::run_mode == RunMode::FIXED_SOURCE && - settings::create_fission_neutrons)) { - create_fission_sites(p); - } - } + // Display information about collision + if ((settings::verbosity >= 10) || p.trace()) { + write_message(fmt::format(" Energy Group = {}", p.g()), 1); + } +} - // If survival biasing is being used, the following subroutine adjusts the - // weight of the particle. Otherwise, it checks to see if absorption occurs. - if (p.macro_xs().absorption > 0.) { - absorption(p); +void sample_reaction(Particle& p) +{ + // Create fission bank sites. Note that while a fission reaction is sampled, + // it never actually "happens", i.e. the weight of the particle does not + // change when sampling fission sites. The following block handles all + // absorption (including fission) + + if (model::materials[p.material()]->fissionable()) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::create_fission_neutrons)) { + create_fission_sites(p); } - if (!p.alive()) - return; - - // Sample a scattering event to determine the energy of the exiting neutron - scatter(p); - - // Play russian roulette if there are no weight windows - if (!settings::weight_windows_on) - apply_russian_roulette(p); } - void scatter(Particle & p) - { - data::mg.macro_xs_[p.material()].sample_scatter(p.g_last(), p.g(), p.mu(), - p.wgt(), p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a); + // If survival biasing is being used, the following subroutine adjusts the + // weight of the particle. Otherwise, it checks to see if absorption occurs. + if (p.macro_xs().absorption > 0.) { + absorption(p); + } + if (!p.alive()) + return; - // Rotate the angle - p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); + // Sample a scattering event to determine the energy of the exiting neutron + scatter(p); - // Update energy value for downstream compatability (in tallying) - p.E() = data::mg.energy_bin_avg_[p.g()]; + // Play russian roulette if there are no weight windows + if (!settings::weight_windows_on) + apply_russian_roulette(p); +} - // Set event component - p.event() = TallyEvent::SCATTER; - } +void scatter(Particle& p) +{ + data::mg.macro_xs_[p.material()].sample_scatter(p.g_last(), p.g(), p.mu(), + p.wgt(), p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a); - void create_fission_sites(Particle & p) - { - // If uniform fission source weighting is turned on, we increase or decrease - // the expected number of fission sites produced - double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; + // Rotate the angle + p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); - // Determine the expected number of neutrons produced - double nu_t = p.wgt() / simulation::keff * weight * - p.macro_xs().nu_fission / p.macro_xs().total; + // Update energy value for downstream compatability (in tallying) + p.E() = data::mg.energy_bin_avg_[p.g()]; - // Sample the number of neutrons produced - int nu = static_cast(nu_t); - if (prn(p.current_seed()) <= (nu_t - int(nu_t))) { - nu++; - } + // Set event component + p.event() = TallyEvent::SCATTER; +} - // If no neutrons were produced then don't continue - if (nu == 0) - return; - - // Initialize the counter of delayed neutrons encountered for each delayed - // group. - double nu_d[MAX_DELAYED_GROUPS] = {0.}; - - // Clear out particle's nu fission bank - p.nu_bank().clear(); - - p.fission() = true; - - // Determine whether to place fission sites into the shared fission bank - // or the secondary particle bank. - bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); - - // Counter for the number of fission sites successfully stored to the shared - // fission bank or the secondary particle bank - int n_sites_stored; - - for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { - // Initialize fission site object with particle data - SourceSite site; - site.r = p.r(); - site.particle = ParticleType::neutron; - site.time = p.time(); - site.wgt = 1. / weight; - - // Sample the cosine of the angle, assuming fission neutrons are emitted - // isotropically - double mu = 2. * prn(p.current_seed()) - 1.; - - // Sample the azimuthal angle uniformly in [0, 2.pi) - double phi = 2. * PI * prn(p.current_seed()); - site.u.x = mu; - site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi); - site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi); - - // Sample secondary energy distribution for the fission reaction - int dg; - int gout; - data::mg.macro_xs_[p.material()].sample_fission_energy(p.g(), dg, gout, - p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a); - - // Store the energy and delayed groups on the fission bank - site.E = gout; - - // We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest - // of the code, 0 is prompt. - site.delayed_group = dg + 1; - - // If delayed product production, sample time of emission - if (dg != -1) { - auto& macro_xs = data::mg.macro_xs_[p.material()]; - double decay_rate = - macro_xs.get_xs(MgxsType::DECAY_RATE, 0, nullptr, nullptr, &dg, 0, 0); - site.time -= std::log(prn(p.current_seed())) / decay_rate; - - // Reject site if it exceeds time cutoff - double t_cutoff = - settings::time_cutoff[static_cast(site.particle)]; - if (site.time > t_cutoff) { - continue; - } - } +void create_fission_sites(Particle& p) +{ + // If uniform fission source weighting is turned on, we increase or decrease + // the expected number of fission sites produced + double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; + + // Determine the expected number of neutrons produced + double nu_t = p.wgt() / simulation::keff * weight * p.macro_xs().nu_fission / + p.macro_xs().total; + + // Sample the number of neutrons produced + int nu = static_cast(nu_t); + if (prn(p.current_seed()) <= (nu_t - int(nu_t))) { + nu++; + } - // Set parent and progeny ID - site.parent_id = p.id(); - site.progeny_id = p.n_progeny()++; - - // Store fission site in bank - if (use_fission_bank) { - int64_t idx = simulation::fission_bank.thread_safe_append(site); - if (idx == -1) { - warning( - "The shared fission bank is full. Additional fission sites created " - "in this generation will not be banked. Results may be " - "non-deterministic."); - - // Decrement number of particle progeny as storage was unsuccessful. - // This step is needed so that the sum of all progeny is equal to the - // size of the shared fission bank. - p.n_progeny()--; - - // Break out of loop as no more sites can be added to fission bank - break; - } - } else { - p.secondary_bank().push_back(site); + // If no neutrons were produced then don't continue + if (nu == 0) + return; + + // Initialize the counter of delayed neutrons encountered for each delayed + // group. + double nu_d[MAX_DELAYED_GROUPS] = {0.}; + + // Clear out particle's nu fission bank + p.nu_bank().clear(); + + p.fission() = true; + + // Determine whether to place fission sites into the shared fission bank + // or the secondary particle bank. + bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); + + // Counter for the number of fission sites successfully stored to the shared + // fission bank or the secondary particle bank + int n_sites_stored; + + for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) { + // Initialize fission site object with particle data + SourceSite site; + site.r = p.r(); + site.particle = ParticleType::neutron; + site.time = p.time(); + site.wgt = 1. / weight; + + // Sample the cosine of the angle, assuming fission neutrons are emitted + // isotropically + double mu = 2. * prn(p.current_seed()) - 1.; + + // Sample the azimuthal angle uniformly in [0, 2.pi) + double phi = 2. * PI * prn(p.current_seed()); + site.u.x = mu; + site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi); + site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi); + + // Sample secondary energy distribution for the fission reaction + int dg; + int gout; + data::mg.macro_xs_[p.material()].sample_fission_energy( + p.g(), dg, gout, p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a); + + // Store the energy and delayed groups on the fission bank + site.E = gout; + + // We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest + // of the code, 0 is prompt. + site.delayed_group = dg + 1; + + // If delayed product production, sample time of emission + if (dg != -1) { + auto& macro_xs = data::mg.macro_xs_[p.material()]; + double decay_rate = + macro_xs.get_xs(MgxsType::DECAY_RATE, 0, nullptr, nullptr, &dg, 0, 0); + site.time -= std::log(prn(p.current_seed())) / decay_rate; + + // Reject site if it exceeds time cutoff + double t_cutoff = settings::time_cutoff[static_cast(site.particle)]; + if (site.time > t_cutoff) { + continue; } + } - // Set the delayed group on the particle as well - p.delayed_group() = dg + 1; - - // Increment the number of neutrons born delayed - if (p.delayed_group() > 0) { - nu_d[dg]++; + // Set parent and progeny ID + site.parent_id = p.id(); + site.progeny_id = p.n_progeny()++; + + // Store fission site in bank + if (use_fission_bank) { + int64_t idx = simulation::fission_bank.thread_safe_append(site); + if (idx == -1) { + warning( + "The shared fission bank is full. Additional fission sites created " + "in this generation will not be banked. Results may be " + "non-deterministic."); + + // Decrement number of particle progeny as storage was unsuccessful. + // This step is needed so that the sum of all progeny is equal to the + // size of the shared fission bank. + p.n_progeny()--; + + // Break out of loop as no more sites can be added to fission bank + break; } - - // Write fission particles to nuBank - NuBank& nu_bank_entry = p.nu_bank().emplace_back(); - nu_bank_entry.wgt = site.wgt; - nu_bank_entry.E = site.E; - nu_bank_entry.delayed_group = site.delayed_group; + } else { + p.secondary_bank().push_back(site); } - // If shared fission bank was full, and no fissions could be added, - // set the particle fission flag to false. - if (n_sites_stored == 0) { - p.fission() = false; - return; + // Set the delayed group on the particle as well + p.delayed_group() = dg + 1; + + // Increment the number of neutrons born delayed + if (p.delayed_group() > 0) { + nu_d[dg]++; } - // Set nu to the number of fission sites successfully stored. If the fission - // bank was not found to be full then these values are already equivalent. - nu = n_sites_stored; + // Write fission particles to nuBank + NuBank& nu_bank_entry = p.nu_bank().emplace_back(); + nu_bank_entry.wgt = site.wgt; + nu_bank_entry.E = site.E; + nu_bank_entry.delayed_group = site.delayed_group; + } - // Store the total weight banked for analog fission tallies - p.n_bank() = nu; - p.wgt_bank() = nu / weight; - for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { - p.n_delayed_bank(d) = nu_d[d]; - } + // If shared fission bank was full, and no fissions could be added, + // set the particle fission flag to false. + if (n_sites_stored == 0) { + p.fission() = false; + return; } - void absorption(Particle & p) - { - if (settings::survival_biasing) { - // Determine weight absorbed in survival biasing - double wgt_absorb = - p.wgt() * p.macro_xs().absorption / p.macro_xs().total; + // Set nu to the number of fission sites successfully stored. If the fission + // bank was not found to be full then these values are already equivalent. + nu = n_sites_stored; - // Adjust weight of particle by the probability of absorption - p.wgt() -= wgt_absorb; + // Store the total weight banked for analog fission tallies + p.n_bank() = nu; + p.wgt_bank() = nu / weight; + for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { + p.n_delayed_bank(d) = nu_d[d]; + } +} - // Score implicit absorpion estimate of keff +void absorption(Particle& p) +{ + if (settings::survival_biasing) { + // Determine weight absorbed in survival biasing + double wgt_absorb = p.wgt() * p.macro_xs().absorption / p.macro_xs().total; + + // Adjust weight of particle by the probability of absorption + p.wgt() -= wgt_absorb; + + // Score implicit absorpion estimate of keff + p.keff_tally_absorption() += + wgt_absorb * p.macro_xs().nu_fission / p.macro_xs().absorption; + } else { + if (p.macro_xs().absorption > prn(p.current_seed()) * p.macro_xs().total) { p.keff_tally_absorption() += - wgt_absorb * p.macro_xs().nu_fission / p.macro_xs().absorption; - } else { - if (p.macro_xs().absorption > - prn(p.current_seed()) * p.macro_xs().total) { - p.keff_tally_absorption() += - p.wgt() * p.macro_xs().nu_fission / p.macro_xs().absorption; - p.wgt() = 0.0; - p.event() = TallyEvent::ABSORB; - } + p.wgt() * p.macro_xs().nu_fission / p.macro_xs().absorption; + p.wgt() = 0.0; + p.event() = TallyEvent::ABSORB; } } +} } // namespace openmc From 50ce042de7412a254e675b1752b47d64840cdc38 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 21:55:53 +0200 Subject: [PATCH 20/24] fix typo --- src/physics_mg.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index e2bc0529ae1..4d928fe26d9 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -38,10 +38,10 @@ void collision_mg(Particle& p) if (ww.is_valid()) { apply_weight_window(p, ww); } else if (p.type() == ParticleType::neutron) { - apply_russian_roulette(p) + apply_russian_roulette(p); } - // If collision checkpoint is disabled, apply russian roulette if - // the particle is outside the weight window domain + // If collision checkpoint is disabled, apply russian roulette if + // the particle is outside the weight window domain } else if (settings::weight_windows_on) { auto ww = search_weight_window(p); if (!ww.is_valid()) From 5b06d2839a5d3ce7749b21d4f8deaa360fbea14b Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 21:59:20 +0200 Subject: [PATCH 21/24] more fixes --- src/physics.cpp | 40 +++++++++++++++++++--------------------- src/physics_mg.cpp | 5 +++-- 2 files changed, 22 insertions(+), 23 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index 299aba473a4..e5fca9d8826 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -70,8 +70,9 @@ void collision(Particle& p) } else if (p.type() == ParticleType::neutron) { apply_russian_roulette(p); } - // If collision checkpoint is disabled, apply russian roulette if - // the particle is a neutron and it is outside the weight window domain + // If collision checkpoint is disabled and weight windows are enabled, + // apply russian roulette if the particle is a neutron and it is outside + // the weight window domain } else if (settings::weight_windows_on && p.type() == ParticleType::neutron) { auto ww = search_weight_window(p); if (!ww.is_valid()) @@ -132,8 +133,7 @@ void sample_neutron_reaction(Particle& p) if (p.secondary_bank().size() >= settings::max_secondaries) { fatal_error( "The secondary particle bank appears to be growing without " - "bound. You are likely running a subcritical multiplication " - "problem " + "bound. You are likely running a subcritical multiplication problem " "with k-effective close to or greater than one."); } } @@ -561,8 +561,8 @@ Reaction& sample_fission(int i_nuclide, Particle& p) const auto& nuc {data::nuclides[i_nuclide]}; // If we're in the URR, by default use the first fission reaction. We also - // default to the first reaction if we know that there are no partial - // fission reactions + // default to the first reaction if we know that there are no partial fission + // reactions if (p.neutron_xs(i_nuclide).use_ptable || !nuc->has_partial_fission_) { return *nuc->fission_rx_[0]; } @@ -784,8 +784,8 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p) // Find speed of neutron in CM vel = v_n.norm(); - // Sample scattering angle, checking if angle distribution is present - // (assume isotropic otherwise) + // Sample scattering angle, checking if angle distribution is present (assume + // isotropic otherwise) double mu_cm; auto& d = rx.products_[0].distribution_[0]; auto d_ = dynamic_cast(d.get()); @@ -817,8 +817,8 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p) p.u() = v_n / vel; // Because of floating-point roundoff, it may be possible for mu_lab to be - // outside of the range [-1,1). In these cases, we just set mu_lab to - // exactly -1 or 1 + // outside of the range [-1,1). In these cases, we just set mu_lab to exactly + // -1 or 1 if (std::abs(p.mu()) > 1.0) p.mu() = std::copysign(1.0, p.mu()); } @@ -849,13 +849,11 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, // sampling method to use sampling_method = settings::res_scat_method; - // upper resonance scattering energy bound (target is at rest above this - // E) + // upper resonance scattering energy bound (target is at rest above this E) if (E > settings::res_scat_energy_max) { return {}; - // lower resonance scattering energy bound (should be no resonances - // below) + // lower resonance scattering energy bound (should be no resonances below) } else if (E < settings::res_scat_energy_min) { sampling_method = ResScatMethod::cxs; } @@ -906,8 +904,8 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, } if (i_E_up == i_E_low) { - // Handle degenerate case -- if the upper/lower bounds occur for the - // same index, then using cxs is probably a good approximation + // Handle degenerate case -- if the upper/lower bounds occur for the same + // index, then using cxs is probably a good approximation return sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed); } @@ -1211,11 +1209,11 @@ void sample_secondary_photons(Particle& p, int i_nuclide) // Sample the new direction Direction u = rotate_angle(p.u(), mu, nullptr, p.current_seed()); - // In a k-eigenvalue simulation, it's necessary to provide higher weight - // to secondary photons from non-fission reactions to properly balance - // energy release and deposition. See D. P. Griesheimer, S. J. Douglass, - // and M. H. Stedry, "Self-consistent energy normalization for quasistatic - // reactor calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. + // In a k-eigenvalue simulation, it's necessary to provide higher weight to + // secondary photons from non-fission reactions to properly balance energy + // release and deposition. See D. P. Griesheimer, S. J. Douglass, and M. H. + // Stedry, "Self-consistent energy normalization for quasistatic reactor + // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. double wgt = photon_wgt; if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) { wgt *= simulation::keff; diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 4d928fe26d9..6acb09d9002 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -40,8 +40,9 @@ void collision_mg(Particle& p) } else if (p.type() == ParticleType::neutron) { apply_russian_roulette(p); } - // If collision checkpoint is disabled, apply russian roulette if - // the particle is outside the weight window domain + // If collision checkpoint is disabled and weight windows are enabled, + // apply russian roulette if the particle is outside the weight window + // domain } else if (settings::weight_windows_on) { auto ww = search_weight_window(p); if (!ww.is_valid()) From 9a0565ab81ccf26eba35f98c537aa7f06896a18e Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 28 Jan 2026 22:54:42 +0200 Subject: [PATCH 22/24] another fix --- src/weight_windows.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 0dd33306e2d..6ad0a5e340f 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -907,10 +907,14 @@ WeightWindow search_weight_window(const Particle& p) { // TODO: this is a linear search - should do something more clever int mesh_bin; + WeightWindow weight_window; for (const auto& ww : variance_reduction::weight_windows) { mesh_bin = ww->get_mesh_bin(p); - if (mesh_bin >= 0) - return ww->get_weight_window(p.E(), mesh_bin); + if (mesh_bin < 0) + continue; + weight_window = ww->get_weight_window(p.E(), mesh_bin); + if (weight_window.is_valid()) + return weight_window; } return {}; } From 5f6c40d246400997165623be452a4f901110a954 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 00:34:25 +0200 Subject: [PATCH 23/24] simplify code and fix bug --- src/physics.cpp | 30 ++++++++++++++++-------------- src/physics_mg.cpp | 16 ++++------------ 2 files changed, 20 insertions(+), 26 deletions(-) diff --git a/src/physics.cpp b/src/physics.cpp index e5fca9d8826..aa49a6e54c5 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -61,22 +61,24 @@ void collision(Particle& p) break; } - // If collision checkpoint is enabled, apply weight window - // if valid and apply global russian roulette if not. - if (settings::weight_window_checkpoint_collision) { - auto ww = search_weight_window(p); - if (ww.is_valid()) { - apply_weight_window(p, ww); + if (settings::weight_windows_on) { + if (settings::weight_window_checkpoint_collision) { + // If collision checkpoint is enabled, apply weight window + // if valid and apply global russian roulette if not. + auto ww = search_weight_window(p); + if (ww.is_valid()) { + apply_weight_window(p, ww); + } else if (p.type() == ParticleType::neutron) { + apply_russian_roulette(p); + } } else if (p.type() == ParticleType::neutron) { - apply_russian_roulette(p); + // If collision checkpoint is disabled and weight windows are enabled, + // apply russian roulette if the particle is a neutron and it is outside + // the weight window domain + auto ww = search_weight_window(p); + if (!ww.is_valid()) + apply_russian_roulette(p); } - // If collision checkpoint is disabled and weight windows are enabled, - // apply russian roulette if the particle is a neutron and it is outside - // the weight window domain - } else if (settings::weight_windows_on && p.type() == ParticleType::neutron) { - auto ww = search_weight_window(p); - if (!ww.is_valid()) - apply_russian_roulette(p); } // Kill particle if energy falls below cutoff diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 6acb09d9002..c8bd65eff1f 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -31,22 +31,14 @@ void collision_mg(Particle& p) // Sample the reaction type sample_reaction(p); - // If collision checkpoint is enabled, apply weight window - // if valid and apply global russian roulette if not. - if (settings::weight_window_checkpoint_collision) { + if (settings::weight_windows_on) { auto ww = search_weight_window(p); if (ww.is_valid()) { - apply_weight_window(p, ww); - } else if (p.type() == ParticleType::neutron) { + if (settings::weight_window_checkpoint_collision) + apply_weight_window(p, ww); + } else { apply_russian_roulette(p); } - // If collision checkpoint is disabled and weight windows are enabled, - // apply russian roulette if the particle is outside the weight window - // domain - } else if (settings::weight_windows_on) { - auto ww = search_weight_window(p); - if (!ww.is_valid()) - apply_russian_roulette(p); } // Display information about collision From a14249281469b7acac59953960ae4ad310847cb8 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 29 Jan 2026 00:40:04 +0200 Subject: [PATCH 24/24] simplify code further --- src/weight_windows.cpp | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 6ad0a5e340f..becf3fc36f2 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -933,8 +933,12 @@ void apply_weight_windows(Particle& p) return; auto ww = search_weight_window(p); - if (ww.is_valid()) + if (ww.is_valid()) { apply_weight_window(p, ww); + } else { + if (p.wgt_ww_born() == -1.0) + p.wgt_ww_born() = 1.0; + } } void apply_weight_window(Particle& p, WeightWindow weight_window) @@ -944,15 +948,10 @@ void apply_weight_window(Particle& p, WeightWindow weight_window) return; // If particle has not yet had its birth weight window value set, set it to - // the current weight window (or 1.0 if not born in a weight window). - if (p.wgt_ww_born() == -1.0) { - if (weight_window.is_valid()) { - p.wgt_ww_born() = - (weight_window.lower_weight + weight_window.upper_weight) / 2; - } else { - p.wgt_ww_born() = 1.0; - } - } + // the current weight window. + if (p.wgt_ww_born() == -1.0) + p.wgt_ww_born() = + (weight_window.lower_weight + weight_window.upper_weight) / 2; // Normalize weight windows based on particle's starting weight // and the value of the weight window the particle was born in.