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

Merge pull request #1978 from paulromano/reaction-xs-method

Browse source 
Add Reaction::xs methods
This commit is contained in:
John Tramm 2022-02-22 16:07:38 -06:00 • committed by GitHub
commit 82e0a5116e
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
4 changed files with 56 additions and 57 deletions
Download patch file Download diff file Expand all files Collapse all files

13
include/openmc/reaction.h
View file

@ -9,6 +9,7 @@
#include "hdf5.h"
#include <gsl/gsl-lite.hpp>
#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<int>& 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

45
src/physics.cpp
View file

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

17
src/reaction.cpp
View file

@ -65,6 +65,23 @@ Reaction::Reaction(hid_t group, const vector<int>& 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<const double> energy, gsl::span<const double> flux,
const vector<double>& grid) const

38
src/tallies/tally_scoring.cpp
View file

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