diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 43acdbba54..b01238611a 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -242,15 +242,6 @@ enum ReactionType { constexpr array 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}; diff --git a/include/openmc/urr.h b/include/openmc/urr.h index 5e3d467a0f..3978c7b86a 100644 --- a/include/openmc/urr.h +++ b/include/openmc/urr.h @@ -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 energy_; //!< incident energies - xt::xtensor prob_; //!< Actual probability tables + + vector 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 cdf_values_; // Note: must be row major! + xt::xtensor 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 diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 6d34f309b4..ae3085aff2 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -247,7 +247,7 @@ Nuclide::Nuclide(hid_t group, const vector& 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(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.; } diff --git a/src/urr.cpp b/src/urr.cpp index 280fb01efb..02cef22809 100644 --- a/src/urr.cpp +++ b/src/urr.cpp @@ -1,5 +1,6 @@ #include "openmc/urr.h" +#include // any_of #include 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 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 \ No newline at end of file +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