diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index 46705acf5e..93987d9d7e 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -9,6 +9,7 @@ #include "hdf5.h" #include +#include "openmc/particle_data.h" #include "openmc/reaction_product.h" #include "openmc/vector.h" @@ -27,6 +28,18 @@ public: //! \param[in] temperatures Desired temperatures for cross sections explicit Reaction(hid_t group, const vector& temperatures); + //! Calculate cross section given temperautre/grid index, interpolation factor + // + //! \param[in] i_temp Temperature index + //! \param[in] i_grid Energy grid index + //! \param[in] interp_factor Interpolation factor between grid points + double xs(gsl::index i_temp, gsl::index i_grid, double interp_factor) const; + + //! Calculate cross section + // + //! \param[in] micro Microscopic cross section cache + double xs(const NuclideMicroXS& micro) const; + //! \brief Calculate reaction rate based on group-wise flux distribution // //! \param[in] i_temp Temperature index diff --git a/src/physics.cpp b/src/physics.cpp index 63d764b677..edc581c427 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -541,23 +541,15 @@ Reaction& sample_fission(int i_nuclide, Particle& p) } } - // Get grid index and interpolatoin factor and sample fission cdf - int i_temp = p.neutron_xs(i_nuclide).index_temp; - int i_grid = p.neutron_xs(i_nuclide).index_grid; - double f = p.neutron_xs(i_nuclide).interp_factor; + // 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_) { - // if energy is below threshold for this reaction, skip it - int threshold = rx->xs_[i_temp].threshold; - if (i_grid < threshold) - continue; - // add to cumulative probability - prob += (1.0 - f) * rx->xs_[i_temp].value[i_grid - threshold] + - f * rx->xs_[i_temp].value[i_grid - threshold + 1]; + prob += rx->xs(micro); // Create fission bank sites if fission occurs if (prob > cutoff) @@ -573,25 +565,20 @@ 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 - int i_temp = p.neutron_xs(i_nuclide).index_temp; - int i_grid = p.neutron_xs(i_nuclide).index_grid; - double f = p.neutron_xs(i_nuclide).interp_factor; - double cutoff = prn(p.current_seed()) * p.neutron_xs(i_nuclide).photon_prod; + 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) { - const auto& rx = nuc->reactions_[i]; - int threshold = rx->xs_[i_temp].threshold; - - // if energy is below threshold for this reaction, skip it - if (i_grid < threshold) - continue; - // Evaluate neutron cross section - double xs = ((1.0 - f) * rx->xs_[i_temp].value[i_grid - threshold] + - f * (rx->xs_[i_temp].value[i_grid - threshold + 1])); + 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) { @@ -663,8 +650,6 @@ void scatter(Particle& p, int i_nuclide) const auto& nuc {data::nuclides[i_nuclide]}; const auto& micro {p.neutron_xs(i_nuclide)}; int i_temp = micro.index_temp; - int i_grid = micro.index_grid; - double f = micro.interp_factor; // For tallying purposes, this routine might be called directly. In that // case, we need to sample a reaction via the cutoff variable @@ -718,14 +703,8 @@ void scatter(Particle& p, int i_nuclide) fatal_error("Did not sample any reaction for nuclide " + nuc->name_); } - // if energy is below threshold for this reaction, skip it - const auto& xs {nuc->reactions_[i]->xs_[i_temp]}; - if (i_grid < xs.threshold) - continue; - // add to cumulative probability - prob += (1.0 - f) * xs.value[i_grid - xs.threshold] + - f * xs.value[i_grid - xs.threshold + 1]; + prob += nuc->reactions_[i]->xs(micro); } // Perform collision physics for inelastic scattering diff --git a/src/reaction.cpp b/src/reaction.cpp index c929677f2c..ec9d32f9a1 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -65,6 +65,23 @@ Reaction::Reaction(hid_t group, const vector& temperatures) } } +double Reaction::xs( + gsl::index i_temp, gsl::index i_grid, double interp_factor) const +{ + // If energy is below threshold, return 0. Otherwise interpolate between + // nearest grid points + const auto& x = xs_[i_temp]; + return (i_grid < x.threshold) + ? 0.0 + : (1.0 - interp_factor) * x.value[i_grid - x.threshold] + + interp_factor * x.value[i_grid - x.threshold + 1]; +} + +double Reaction::xs(const NuclideMicroXS& micro) const +{ + return this->xs(micro.index_temp, micro.index_grid, micro.interp_factor); +} + double Reaction::collapse_rate(gsl::index i_temp, gsl::span energy, gsl::span flux, const vector& grid) const diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index bd516ca18a..b06570b992 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -251,20 +251,16 @@ double get_nuclide_neutron_heating( if (mt == C_NONE) return 0.0; - auto i_temp = p.neutron_xs(i_nuclide).index_temp; + const auto& micro = p.neutron_xs(i_nuclide); + auto i_temp = micro.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; - // Determine total kerma - auto f = p.neutron_xs(i_nuclide).interp_factor; - double kerma = (1.0 - f) * xs.value[i_grid - xs.threshold] + - f * xs.value[i_grid - xs.threshold + 1]; + const auto& rx {*nuc.reactions_[mt]}; + double kerma = rx.xs(micro); + if (kerma == 0.0) + return 0.0; if (settings::run_mode == RunMode::EIGENVALUE) { // Determine kerma for fission as (EFR + EB)*sigma_f @@ -477,7 +473,7 @@ double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin) auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) return 0.0; - const auto& rxn {*nuc.reactions_[m]}; + const auto& rx {*nuc.reactions_[m]}; const auto& micro {p.neutron_xs(i_nuclide)}; // In the URR, the (n,gamma) cross section is sampled randomly from @@ -494,14 +490,7 @@ double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin) auto f = micro.interp_factor; // Calculate interpolated cross section - const auto& xs {rxn.xs_[i_temp]}; - double value; - if (i_grid >= xs.threshold) { - value = ((1.0 - f) * xs.value[i_grid - xs.threshold] + - f * xs.value[i_grid - xs.threshold + 1]); - } else { - value = 0.0; - } + double xs = rx.xs(micro); if (settings::run_mode == RunMode::EIGENVALUE && score_bin == HEATING_LOCAL) { @@ -515,18 +504,18 @@ double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin) : 0.0; // Determine non-fission kerma as difference - double kerma_non_fission = value - kerma_fission; + double kerma_non_fission = xs - kerma_fission; // Re-weight non-fission kerma by keff 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. - value = simulation::keff * kerma_non_fission + kerma_fission; + xs = simulation::keff * kerma_non_fission + kerma_fission; } - return value; + return xs; } else { // For multipole, calculate (n,gamma) from other reactions - return rxn.mt_ == N_GAMMA ? micro.absorption - micro.fission : 0.0; + return rx.mt_ == N_GAMMA ? micro.absorption - micro.fission : 0.0; } return 0.0; } @@ -1497,7 +1486,8 @@ void score_general_mg(Particle& p, int i_tally, int start_index, } // For shorthand, assign pointers to the material and nuclide xs set - auto& nuc_xs = (i_nuclide >= 0) ? data::mg.nuclides_[i_nuclide] : data::mg.macro_xs_[p.material()]; + auto& nuc_xs = (i_nuclide >= 0) ? data::mg.nuclides_[i_nuclide] + : data::mg.macro_xs_[p.material()]; auto& macro_xs = data::mg.macro_xs_[p.material()]; // Find the temperature and angle indices of interest