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

Merge pull request #1904 from gridley/urr_improvements

Browse source 
improve URR lookup data layout and performance
This commit is contained in:
Paul Romano 2021-11-11 22:12:12 -06:00 • committed by GitHub
commit c8003e0243
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
4 changed files with 127 additions and 63 deletions
Download patch file Download diff file Expand all files Collapse all files

9
include/openmc/constants.h
View file

@ -242,15 +242,6 @@ enum ReactionType {
constexpr array<int, 6> DEPLETION_RX {N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N};
enum class URRTableParam {
CUM_PROB,
TOTAL,
ELASTIC,
FISSION,
N_GAMMA,
HEATING
};
// Maximum number of partial fission reactions
constexpr int PARTIAL_FISSION_MAX {4};

37
include/openmc/urr.h
View file

@ -7,6 +7,7 @@
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/vector.h"
namespace openmc {
@ -16,16 +17,46 @@ namespace openmc {
class UrrData {
public:
// Since we access all of these at once, we want
// them contiguous in memory.
struct XSSet {
double total;
double elastic;
double fission;
double n_gamma;
double heating;
};
Interpolation interp_; //!< interpolation type
int inelastic_flag_; //!< inelastic competition flag
int absorption_flag_; //!< other absorption flag
bool multiply_smooth_; //!< multiply by smooth cross section?
int n_energy_; //!< number of energy points
xt::xtensor<double, 1> energy_; //!< incident energies
xt::xtensor<double, 3> prob_; //!< Actual probability tables
vector<double> energy_; //!< incident energies
auto n_energy() const { return energy_.size(); }
/* The row indexes correspond to the incident energy table, and column
* indices correspond to values of the CDF at that energy. For the CDF matrix
* below, obviously, values of the CDF are stored. For the xs_values
* variable, the columns line up with the index of cdf_values.
*/
xt::xtensor<double, 2> cdf_values_; // Note: must be row major!
xt::xtensor<XSSet, 2> xs_values_;
// Number of points in the CDF
auto n_cdf() const { return cdf_values_.shape()[1]; }
//! \brief Load the URR data from the provided HDF5 group
explicit UrrData(hid_t group_id);
// Checks if any negative CDF or XS values are present
bool has_negative() const;
// Checks if the passed energy is within the bounds of the URR table
bool energy_in_bounds(double E) const
{
return energy_.front() < E && E < energy_.back();
}
};
} // namespace openmc

82
src/nuclide.cpp
View file

@ -247,7 +247,7 @@ Nuclide::Nuclide(hid_t group, const vector<double>& temperature)
close_group(urr_group);
// Check for negative values
if (xt::any(urr_data_[i].prob_ < 0.) && mpi::master) {
if (urr_data_[i].has_negative() && mpi::master) {
warning("Negative value(s) found on probability table for nuclide " +
name_ + " at " + temp_str);
}
@ -775,11 +775,8 @@ void Nuclide::calculate_xs(
// If the particle is in the unresolved resonance range and there are
// probability tables, we need to determine cross sections from the table
if (settings::urr_ptables_on && urr_present_ && !use_mp) {
int n = urr_data_[micro.index_temp].n_energy_;
if ((p.E() > urr_data_[micro.index_temp].energy_(0)) &&
(p.E() < urr_data_[micro.index_temp].energy_(n - 1))) {
if (urr_data_[micro.index_temp].energy_in_bounds(p.E()))
this->calculate_urr_xs(micro.index_temp, p);
}
}
micro.last_E = p.E();
@ -825,10 +822,8 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
const auto& urr = urr_data_[i_temp];
// Determine the energy table
int i_energy = 0;
while (p.E() >= urr.energy_(i_energy + 1)) {
++i_energy;
};
int i_energy =
lower_bound_index(urr.energy_.begin(), urr.energy_.end(), p.E());
// Sample the probability table using the cumulative distribution
@ -840,15 +835,13 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
double r = future_prn(static_cast<int64_t>(index_), *p.current_seed());
p.stream() = STREAM_TRACKING;
int i_low = 0;
while (urr.prob_(i_energy, URRTableParam::CUM_PROB, i_low) <= r) {
++i_low;
};
int i_up = 0;
while (urr.prob_(i_energy + 1, URRTableParam::CUM_PROB, i_up) <= r) {
++i_up;
};
// Warning: this assumes row-major order of cdf_values_
int i_low = upper_bound_index(&urr.cdf_values_(i_energy, 0),
&urr.cdf_values_(i_energy, 0) + urr.n_cdf(), r) +
1;
int i_up = upper_bound_index(&urr.cdf_values_(i_energy + 1, 0),
&urr.cdf_values_(i_energy + 1, 0) + urr.n_cdf(), r) +
1;
// Determine elastic, fission, and capture cross sections from the
// probability table
@ -858,49 +851,46 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
double f;
if (urr.interp_ == Interpolation::lin_lin) {
// Determine the interpolation factor on the table
f = (p.E() - urr.energy_(i_energy)) /
(urr.energy_(i_energy + 1) - urr.energy_(i_energy));
f = (p.E() - urr.energy_[i_energy]) /
(urr.energy_[i_energy + 1] - urr.energy_[i_energy]);
elastic = (1. - f) * urr.prob_(i_energy, URRTableParam::ELASTIC, i_low) +
f * urr.prob_(i_energy + 1, URRTableParam::ELASTIC, i_up);
fission = (1. - f) * urr.prob_(i_energy, URRTableParam::FISSION, i_low) +
f * urr.prob_(i_energy + 1, URRTableParam::FISSION, i_up);
capture = (1. - f) * urr.prob_(i_energy, URRTableParam::N_GAMMA, i_low) +
f * urr.prob_(i_energy + 1, URRTableParam::N_GAMMA, i_up);
elastic = (1. - f) * urr.xs_values_(i_energy, i_low).elastic +
f * urr.xs_values_(i_energy + 1, i_up).elastic;
fission = (1. - f) * urr.xs_values_(i_energy, i_low).fission +
f * urr.xs_values_(i_energy + 1, i_up).fission;
capture = (1. - f) * urr.xs_values_(i_energy, i_low).n_gamma +
f * urr.xs_values_(i_energy + 1, i_up).n_gamma;
} else if (urr.interp_ == Interpolation::log_log) {
// Determine interpolation factor on the table
f = std::log(p.E() / urr.energy_(i_energy)) /
std::log(urr.energy_(i_energy + 1) / urr.energy_(i_energy));
f = std::log(p.E() / urr.energy_[i_energy]) /
std::log(urr.energy_[i_energy + 1] / urr.energy_[i_energy]);
// Calculate the elastic cross section/factor
if ((urr.prob_(i_energy, URRTableParam::ELASTIC, i_low) > 0.) &&
(urr.prob_(i_energy + 1, URRTableParam::ELASTIC, i_up) > 0.)) {
elastic = std::exp(
(1. - f) *
std::log(urr.prob_(i_energy, URRTableParam::ELASTIC, i_low)) +
f * std::log(urr.prob_(i_energy + 1, URRTableParam::ELASTIC, i_up)));
if ((urr.xs_values_(i_energy, i_low).elastic > 0.) &&
(urr.xs_values_(i_energy + 1, i_up).elastic > 0.)) {
elastic =
std::exp((1. - f) * std::log(urr.xs_values_(i_energy, i_low).elastic) +
f * std::log(urr.xs_values_(i_energy + 1, i_up).elastic));
} else {
elastic = 0.;
}
// Calculate the fission cross section/factor
if ((urr.prob_(i_energy, URRTableParam::FISSION, i_low) > 0.) &&
(urr.prob_(i_energy + 1, URRTableParam::FISSION, i_up) > 0.)) {
fission = std::exp(
(1. - f) *
std::log(urr.prob_(i_energy, URRTableParam::FISSION, i_low)) +
f * std::log(urr.prob_(i_energy + 1, URRTableParam::FISSION, i_up)));
if ((urr.xs_values_(i_energy, i_low).fission > 0.) &&
(urr.xs_values_(i_energy + 1, i_up).fission > 0.)) {
fission =
std::exp((1. - f) * std::log(urr.xs_values_(i_energy, i_low).fission) +
f * std::log(urr.xs_values_(i_energy + 1, i_up).fission));
} else {
fission = 0.;
}
// Calculate the capture cross section/factor
if ((urr.prob_(i_energy, URRTableParam::N_GAMMA, i_low) > 0.) &&
(urr.prob_(i_energy + 1, URRTableParam::N_GAMMA, i_up) > 0.)) {
capture = std::exp(
(1. - f) *
std::log(urr.prob_(i_energy, URRTableParam::N_GAMMA, i_low)) +
f * std::log(urr.prob_(i_energy + 1, URRTableParam::N_GAMMA, i_up)));
if ((urr.xs_values_(i_energy, i_low).n_gamma > 0.) &&
(urr.xs_values_(i_energy + 1, i_up).n_gamma > 0.)) {
capture =
std::exp((1. - f) * std::log(urr.xs_values_(i_energy, i_low).n_gamma) +
f * std::log(urr.xs_values_(i_energy + 1, i_up).n_gamma));
} else {
capture = 0.;
}

62
src/urr.cpp
View file

@ -1,5 +1,6 @@
#include "openmc/urr.h"
#include <algorithm> // any_of
#include <iostream>
namespace openmc {
@ -21,11 +22,62 @@ UrrData::UrrData(hid_t group_id)
// read the energies at which tables exist
read_dataset(group_id, "energy", energy_);
// Set n_energy_
n_energy_ = energy_.shape()[0];
// Read URR tables. The HDF5 format is a little
// different from how we want it laid out in memory.
// This array used to be called "prob_".
xt::xtensor<double, 3> tmp_prob;
read_dataset(group_id, "table", tmp_prob);
auto shape = tmp_prob.shape();
// Read URR tables
read_dataset(group_id, "table", prob_);
// We separate out into two matrices (one with CDF values,
// the other with cross section sets) in order to improve
// contiguity of memory accesses.
const auto n_energy = shape[0];
const auto n_cdf_values = shape[2];
cdf_values_.resize({n_energy, n_cdf_values});
xs_values_.resize({n_energy, n_cdf_values});
// Now fill in the values. Using manual loops here since we might
// not have fancy xtensor slicing code written for GPU tensors.
// The below enum gives how URR tables are laid out in our HDF5 tables.
enum class URRTableParam {
CUM_PROB,
TOTAL,
ELASTIC,
FISSION,
N_GAMMA,
HEATING
};
for (int i_energy = 0; i_energy < n_energy; ++i_energy) {
for (int i_cdf = 0; i_cdf < n_cdf_values; ++i_cdf) {
cdf_values_(i_energy, i_cdf) =
tmp_prob(i_energy, URRTableParam::CUM_PROB, i_cdf);
xs_values_(i_energy, i_cdf).total =
tmp_prob(i_energy, URRTableParam::TOTAL, i_cdf);
xs_values_(i_energy, i_cdf).elastic =
tmp_prob(i_energy, URRTableParam::ELASTIC, i_cdf);
xs_values_(i_energy, i_cdf).fission =
tmp_prob(i_energy, URRTableParam::FISSION, i_cdf);
xs_values_(i_energy, i_cdf).n_gamma =
tmp_prob(i_energy, URRTableParam::N_GAMMA, i_cdf);
xs_values_(i_energy, i_cdf).heating =
tmp_prob(i_energy, URRTableParam::HEATING, i_cdf);
}
}
}
} // namespace openmc
bool UrrData::has_negative() const
{
// Lambda checks if any value in XSSset is negative
auto xs_set_negative = [](const XSSet& xs) {
return xs.total < 0.0 || xs.elastic < 0.0 || xs.fission < 0.0 ||
xs.n_gamma < 0.0 || xs.heating < 0.0;
};
return std::any_of(cdf_values_.begin(), cdf_values_.end(), [](double x) {
return x < 0.0;
}) || std::any_of(xs_values_.begin(), xs_values_.end(), xs_set_negative);
}
} // namespace openmc