From 99bff4f23166c94ebe585ef4abfca9f48684dc3f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 30 Aug 2019 14:22:59 -0500 Subject: [PATCH 1/3] Use helper function when scoring fission_q_* Add a brief helper function get_nuc_fission_q in src/tallies/tally_scoring.cpp that is responsible for returning the correct fission q value for a given score. The two allowed scores are SCORE_FISS_Q_PROMPT and SCORE_FISS_Q_RECOV. The helper function examines the correct dataset on the nuclide object and evaluates the function given the incident particle energy or returns zero. The default value of zero is returned if the fission_q_prompt_ or fission_q_recov_ datasets are empty, or the score is not one of the expected types. --- src/tallies/tally_scoring.cpp | 68 ++++++++++++++--------------------- 1 file changed, 26 insertions(+), 42 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 6bf9ca6e9f..3016799026 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -171,6 +171,22 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index) dg_match.bins_[i_bin] = original_bin; } +//! Helper function to retrieve fission q value from a nuclide + +double get_nuc_fission_q(const Nuclide& nuc, const Particle* p, int score_bin) +{ + if (score_bin == SCORE_FISS_Q_PROMPT) { + if (nuc.fission_q_prompt_) { + return (*nuc.fission_q_prompt_)(p->E_last_); + } + } else if (score_bin == SCORE_FISS_Q_RECOV) { + if (nuc.fission_q_recov_) { + return (*nuc.fission_q_recov_)(p->E_last_); + } + } + return 0.0; +} + //! Helper function for nu-fission tallies with energyout filters. // //! In this case, we may need to score to multiple bins if there were multiple @@ -1040,17 +1056,9 @@ score_general_ce(Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - const auto& nuc {*data::nuclides[p->event_nuclide_]}; + const Nuclide& nuc {*data::nuclides[p->event_nuclide_]}; if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { - double q_value = 0.; - if (score_bin == SCORE_FISS_Q_PROMPT) { - if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->E_last_); - } else if (score_bin == SCORE_FISS_Q_RECOV) { - if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->E_last_); - } - score = p->wgt_absorb_ * q_value + score = p->wgt_absorb_ * get_nuc_fission_q(nuc, p, score_bin) * p->neutron_xs_[p->event_nuclide_].fission / p->neutron_xs_[p->event_nuclide_].absorption * flux; } @@ -1060,51 +1068,27 @@ score_general_ce(Particle* p, int i_tally, int start_index, // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - const auto& nuc {*data::nuclides[p->event_nuclide_]}; + const Nuclide& nuc {*data::nuclides[p->event_nuclide_]}; if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { - double q_value = 0.; - if (score_bin == SCORE_FISS_Q_PROMPT) { - if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->E_last_); - } else if (score_bin == SCORE_FISS_Q_RECOV) { - if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->E_last_); - } - score = p->wgt_last_ * q_value + score = p->wgt_last_ * get_nuc_fission_q(nuc, p, score_bin) * p->neutron_xs_[p->event_nuclide_].fission / p->neutron_xs_[p->event_nuclide_].absorption * flux; } } } else { if (i_nuclide >= 0) { - const auto& nuc {*data::nuclides[i_nuclide]}; - double q_value = 0.; - if (score_bin == SCORE_FISS_Q_PROMPT) { - if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->E_last_); - } else if (score_bin == SCORE_FISS_Q_RECOV) { - if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->E_last_); - } - score = q_value * p->neutron_xs_[i_nuclide].fission - * atom_density * flux; + const Nuclide& nuc {*data::nuclides[i_nuclide]}; + score = get_nuc_fission_q(nuc, p, score_bin) + * p->neutron_xs_[i_nuclide].fission * atom_density * flux; } else { if (p->material_ != MATERIAL_VOID) { const Material& material {*model::materials[p->material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - const auto& nuc {*data::nuclides[j_nuclide]}; - double q_value = 0.; - if (score_bin == SCORE_FISS_Q_PROMPT) { - if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->E_last_); - } else if (score_bin == SCORE_FISS_Q_RECOV) { - if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->E_last_); - } - score += q_value * p->neutron_xs_[j_nuclide].fission - * atom_density * flux; + const Nuclide& nuc {*data::nuclides[j_nuclide]}; + score += get_nuc_fission_q(nuc, p, score_bin) + * p->neutron_xs_[j_nuclide].fission * atom_density * flux; } } } From f2e14dc6bef6156c56fe5fa7e2b775183aaf625f Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 30 Aug 2019 15:29:00 -0500 Subject: [PATCH 2/3] Separate fission-q scoring into separate function In order to reduce the amount of duplicated code in building the energy deposition tally, the logic for scoring SCORE_FISS_Q_PROMPT and SCORE_FISS_Q_RECOV has been pulled into a new function, score_fission_q. The code is almost entirely copied over, with some minor cleanup operations. --- src/tallies/tally_scoring.cpp | 97 +++++++++++++++++++---------------- 1 file changed, 52 insertions(+), 45 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 3016799026..bbbd3fefb9 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -187,6 +187,56 @@ double get_nuc_fission_q(const Nuclide& nuc, const Particle* p, int score_bin) return 0.0; } +double score_fission_q(const Particle* p, int score_bin, const Tally& tally, + double flux, int i_nuclide, double atom_density) +{ + if (tally.estimator_ == ESTIMATOR_ANALOG) { + const Nuclide& nuc {*data::nuclides[p->event_nuclide_]}; + if (settings::survival_biasing) { + // No fission events occur if survival biasing is on -- need to + // calculate fraction of absorptions that would have resulted in + // fission scaled by the Q-value + if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { + return p->wgt_absorb_ * get_nuc_fission_q(nuc, p, score_bin) + * p->neutron_xs_[p->event_nuclide_].fission * flux + / p->neutron_xs_[p->event_nuclide_].absorption; + } + } else { + // Skip any non-absorption events + if (p->event_ == EVENT_SCATTER) return 0.0; + // All fission events will contribute, so again we can use particle's + // weight entering the collision as the estimate for the fission + // reaction rate + if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { + return p->wgt_last_ * get_nuc_fission_q(nuc, p, score_bin) + * p->neutron_xs_[p->event_nuclide_].fission * flux + / p->neutron_xs_[p->event_nuclide_].absorption; + } + } + } else { + if (i_nuclide >= 0) { + const Nuclide& nuc {*data::nuclides[i_nuclide]}; + return get_nuc_fission_q(nuc, p, score_bin) * atom_density * flux + * p->neutron_xs_[i_nuclide].fission; + } else { + if (p->material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p->material_]}; + double score {0.0}; + for (auto i = 0; i < material.nuclide_.size(); ++i) { + auto j_nuclide = material.nuclide_[i]; + auto atom_density = material.atom_density_(i); + const Nuclide& nuc {*data::nuclides[j_nuclide]}; + score += get_nuc_fission_q(nuc, p, score_bin) * atom_density + * p->neutron_xs_[j_nuclide].fission; + } + return score * flux; + } + } + } + return 0.0; +} + + //! Helper function for nu-fission tallies with energyout filters. // //! In this case, we may need to score to multiple bins if there were multiple @@ -339,7 +389,7 @@ void score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, int i_nuclide, double atom_density, double flux) { - auto& tally {*model::tallies[i_tally]}; + Tally& tally {*model::tallies[i_tally]}; // Get the pre-collision energy of the particle. auto E = p->E_last_; @@ -1048,51 +1098,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, case SCORE_FISS_Q_PROMPT: case SCORE_FISS_Q_RECOV: - //continue; if (p->macro_xs_.absorption == 0.) continue; - score = 0.; - if (tally.estimator_ == ESTIMATOR_ANALOG) { - if (settings::survival_biasing) { - // No fission events occur if survival biasing is on -- need to - // calculate fraction of absorptions that would have resulted in - // fission scaled by the Q-value - const Nuclide& nuc {*data::nuclides[p->event_nuclide_]}; - if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { - score = p->wgt_absorb_ * get_nuc_fission_q(nuc, p, score_bin) - * p->neutron_xs_[p->event_nuclide_].fission - / p->neutron_xs_[p->event_nuclide_].absorption * flux; - } - } else { - // Skip any non-absorption events - if (p->event_ == EVENT_SCATTER) continue; - // All fission events will contribute, so again we can use particle's - // weight entering the collision as the estimate for the fission - // reaction rate - const Nuclide& nuc {*data::nuclides[p->event_nuclide_]}; - if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { - score = p->wgt_last_ * get_nuc_fission_q(nuc, p, score_bin) - * p->neutron_xs_[p->event_nuclide_].fission - / p->neutron_xs_[p->event_nuclide_].absorption * flux; - } - } - } else { - if (i_nuclide >= 0) { - const Nuclide& nuc {*data::nuclides[i_nuclide]}; - score = get_nuc_fission_q(nuc, p, score_bin) - * p->neutron_xs_[i_nuclide].fission * atom_density * flux; - } else { - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; - for (auto i = 0; i < material.nuclide_.size(); ++i) { - auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); - const Nuclide& nuc {*data::nuclides[j_nuclide]}; - score += get_nuc_fission_q(nuc, p, score_bin) - * p->neutron_xs_[j_nuclide].fission * atom_density * flux; - } - } - } - } + score = score_fission_q(p, score_bin, tally, flux, i_nuclide, atom_density); break; From 390953050b88e848a828e3063400854a59664de6 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Sep 2019 16:00:28 -0500 Subject: [PATCH 3/3] Pull out neutron heating scoring to standalone function Neutron heating logic has been pulled out from the main score_general_ce function into a single helper function score_neutron_heating. This function returns a double that is the neutron heating, using analog, collision, or track-length estimators, incorporating survival biasing, and the possiblity of scoring across all nuclides in a given material. The functionality is identical to the original implementation, just condensed into two functions to reduce repeated code. The second helper function digs through neutron data to obtain and interpolate kerma coefficients for a given nuclide and reaction. --- src/tallies/tally_scoring.cpp | 169 ++++++++++++++-------------------- 1 file changed, 70 insertions(+), 99 deletions(-) diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index bbbd3fefb9..1f6a8010a5 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -187,6 +187,11 @@ double get_nuc_fission_q(const Nuclide& nuc, const Particle* p, int score_bin) return 0.0; } +//! Helper function to score fission energy +// +//! Pulled out to support both the fission_q scores and energy deposition +//! score + double score_fission_q(const Particle* p, int score_bin, const Tally& tally, double flux, int i_nuclide, double atom_density) { @@ -236,6 +241,69 @@ double score_fission_q(const Particle* p, int score_bin, const Tally& tally, return 0.0; } +//! Helper function to obtain the kerma coefficient for a given nuclide + +double get_nuclide_neutron_heating(const Particle* p, const Nuclide& nuc, + int rxn_index, int i_nuclide) +{ + size_t mt = nuc.reaction_index_[rxn_index]; + if (mt == C_NONE) return 0.0; + auto i_temp = p->neutron_xs_[i_nuclide].index_temp; + if (i_temp < 0) return 0.0; // Can be true due to multipole + const auto& rxn {*nuc.reactions_[mt]}; + const auto& xs {rxn.xs_[i_temp]}; + auto i_grid = p->neutron_xs_[i_nuclide].index_grid; + if (i_grid < xs.threshold) return 0.0; + auto f = p->neutron_xs_[i_nuclide].interp_factor; + return (1.0 - f) * xs.value[i_grid-xs.threshold] + + f * xs.value[i_grid-xs.threshold+1]; +} + +//! Helper function to obtain neutron heating [eV] + +double score_neutron_heating(const Particle* p, const Tally& tally, double flux, + int rxn_bin, int i_nuclide, double atom_density) +{ + double score; + // Get heating macroscopic "cross section" + double heating_xs; + if (i_nuclide >= 0) { + const Nuclide& nuc {*data::nuclides[i_nuclide]}; + heating_xs = get_nuclide_neutron_heating(p, nuc, rxn_bin, i_nuclide); + if (tally.estimator_ == ESTIMATOR_ANALOG) { + heating_xs /= p->neutron_xs_[i_nuclide].total; + } else { + heating_xs *= atom_density; + } + } else { + if (p->material_ != MATERIAL_VOID) { + heating_xs = 0.0; + const Material& material {*model::materials[p->material_]}; + for (auto i = 0; i< material.nuclide_.size(); ++i) { + int j_nuclide = material.nuclide_[i]; + double atom_density {material.atom_density_(i)}; + const Nuclide& nuc {*data::nuclides[j_nuclide]}; + heating_xs += atom_density * get_nuclide_neutron_heating(p, nuc, rxn_bin, j_nuclide); + } + if (tally.estimator_ == ESTIMATOR_ANALOG) { + heating_xs /= p->macro_xs_.total; + } + } + } + score = heating_xs * flux; + if (tally.estimator_ == ESTIMATOR_ANALOG) { + // All events score to a heating tally bin. We actually use a + // collision estimator in place of an analog one since there is no + // reaction-wise heating cross section + if (settings::survival_biasing) { + // Account for the fact that some weight has been absorbed + score *= p->wgt_last_ + p->wgt_absorb_; + } else { + score *= p->wgt_last_; + } + } + return score; +} //! Helper function for nu-fission tallies with energyout filters. // @@ -1145,105 +1213,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, case SCORE_HEATING: score = 0.; if (p->type_ == Particle::Type::neutron) { - if (tally.estimator_ == ESTIMATOR_ANALOG) { - // All events score to a heating tally bin. We actually use a - // collision estimator in place of an analog one since there is no - // reaction-wise heating cross section - if (settings::survival_biasing) { - // We need to account for the fact that some weight was already - // absorbed - score = p->wgt_last_ + p->wgt_absorb_; - } else { - score = p->wgt_last_; - } - if (i_nuclide >= 0) { - // Calculate nuclide heating cross section - double macro_heating = 0.; - const auto& nuc {*data::nuclides[i_nuclide]}; - auto m = nuc.reaction_index_[NEUTRON_HEATING]; - if (m == C_NONE) continue; - const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p->neutron_xs_[i_nuclide].index_temp; - if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p->neutron_xs_[i_nuclide].index_grid; - auto f = p->neutron_xs_[i_nuclide].interp_factor; - const auto& xs {rxn.xs_[i_temp]}; - if (i_grid >= xs.threshold) { - macro_heating = ((1.0 - f) * xs.value[i_grid-xs.threshold] - + f * xs.value[i_grid-xs.threshold+1]); - } - } - score *= macro_heating * flux / p->neutron_xs_[i_nuclide].total; - } else { - if (p->material_ != MATERIAL_VOID) { - // Calculate material heating cross section - double macro_heating = 0.; - const Material& material {*model::materials[p->material_]}; - for (auto i = 0; i < material.nuclide_.size(); ++i) { - auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); - const auto& nuc {*data::nuclides[j_nuclide]}; - auto m = nuc.reaction_index_[NEUTRON_HEATING]; - if (m == C_NONE) continue; - const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p->neutron_xs_[j_nuclide].index_temp; - if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p->neutron_xs_[j_nuclide].index_grid; - auto f = p->neutron_xs_[j_nuclide].interp_factor; - const auto& xs {rxn.xs_[i_temp]}; - if (i_grid >= xs.threshold) { - macro_heating += ((1.0 - f) * xs.value[i_grid-xs.threshold] - + f * xs.value[i_grid-xs.threshold+1]) * atom_density; - } - } - } - score *= macro_heating * flux / p->macro_xs_.total; - } else { - score = 0.; - } - } - } else { - // Calculate neutron heating cross section on-the-fly - if (i_nuclide >= 0) { - const auto& nuc {*data::nuclides[i_nuclide]}; - auto m = nuc.reaction_index_[NEUTRON_HEATING]; - if (m == C_NONE) continue; - const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p->neutron_xs_[i_nuclide].index_temp; - if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p->neutron_xs_[i_nuclide].index_grid; - auto f = p->neutron_xs_[i_nuclide].interp_factor; - const auto& xs {rxn.xs_[i_temp]}; - if (i_grid >= xs.threshold) { - score = ((1.0 - f) * xs.value[i_grid-xs.threshold] - + f * xs.value[i_grid-xs.threshold+1]) * atom_density * flux; - } - } - } else { - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; - for (auto i = 0; i < material.nuclide_.size(); ++i) { - auto j_nuclide = material.nuclide_[i]; - auto atom_density = material.atom_density_(i); - const auto& nuc {*data::nuclides[j_nuclide]}; - auto m = nuc.reaction_index_[NEUTRON_HEATING]; - if (m == C_NONE) continue; - const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p->neutron_xs_[j_nuclide].index_temp; - if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p->neutron_xs_[j_nuclide].index_grid; - auto f = p->neutron_xs_[j_nuclide].interp_factor; - const auto& xs {rxn.xs_[i_temp]}; - if (i_grid >= xs.threshold) { - score += ((1.0 - f) * xs.value[i_grid-xs.threshold] - + f * xs.value[i_grid-xs.threshold+1]) * atom_density - * flux; - } - } - } - } - } - } + score = score_neutron_heating(p, tally, flux, NEUTRON_HEATING, + i_nuclide, atom_density); } else if (p->type_ == Particle::Type::photon) { if (tally.estimator_ == ESTIMATOR_ANALOG) { // Score direct energy deposition in the collision