crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Pull out neutron heating scoring to standalone function

Browse source 
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.
This commit is contained in:
Andrew Johnson 2019-09-04 16:00:28 -05:00
parent f2e14dc6be
commit 390953050b
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
1 changed files with 70 additions and 99 deletions
Download patch file Download diff file Expand all files Collapse all files

169
src/tallies/tally_scoring.cpp
View file

@ -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