From 390953050b88e848a828e3063400854a59664de6 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 4 Sep 2019 16:00:28 -0500 Subject: [PATCH] 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