diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 313b841e2..0c940342c 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -271,7 +271,6 @@ void read_dataset(hid_t obj_id, const char* name, xt::xarray& arr, bool indep close_dataset(dset); } - template void read_dataset_as_shape(hid_t obj_id, const char* name, xt::xtensor& arr, bool indep=false) diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 67444a106..455dab5fc 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -107,7 +107,7 @@ public: //! \brief Determines cross sections in the unresolved resonance range //! from probability tables. - void calculate_urr_xs(const int i_temp, const double E); + void calculate_urr_xs(int i_temp, double E); // Data members std::string name_; //! Name of nuclide, e.g. "U235" @@ -175,8 +175,8 @@ extern "C" void set_micro_xs(); extern "C" bool nuclide_wmp_present(int i_nuclide); extern "C" double nuclide_wmp_emin(int i_nuclide); extern "C" double nuclide_wmp_emax(int i_nuclide); -extern "C" void nuclide_calculate_urr_xs(const bool use_mp, const int i_nuclide, - const int i_temp, const double E); +extern "C" void nuclide_calculate_urr_xs(bool use_mp, int i_nuclide, + int i_temp, double E); } // namespace openmc diff --git a/include/openmc/urr.h b/include/openmc/urr.h index 36ea07332..676425442 100644 --- a/include/openmc/urr.h +++ b/include/openmc/urr.h @@ -16,17 +16,16 @@ namespace openmc { class UrrData{ public: - 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 + 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 //! \brief Load the URR data from the provided HDF5 group - void - from_hdf5(hid_t group_id); + explicit UrrData(hid_t group_id); }; } // namespace openmc diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 727968a4d..f7ef87b56 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -196,7 +196,7 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide) // Read unresolved resonance probability tables if present if (object_exists(group, "urr")) { urr_present_ = true; - urr_data_.resize(temps_to_read.size()); + urr_data_.reserve(temps_to_read.size()); for (int i = 0; i < temps_to_read.size(); i++) { // Get temperature as a string @@ -204,7 +204,7 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide) // Read probability tables for i-th temperature hid_t urr_group = open_group(group, ("urr/" + temp_str).c_str()); - urr_data_[i].from_hdf5(urr_group); + urr_data_.emplace_back(urr_group); close_group(urr_group); // Check for negative values @@ -404,25 +404,22 @@ double Nuclide::elastic_xs_0K(double E) const return (1.0 - f)*elastic_0K_[i_grid] + f*elastic_0K_[i_grid + 1]; } -void Nuclide::calculate_urr_xs(const int i_temp, const double E) +void Nuclide::calculate_urr_xs(int i_temp, double E) { auto& micro = simulation::micro_xs[i_nuclide_]; micro.use_ptable = true; // Create a shorthand for the URR data - UrrData* urr = &(urr_data_[i_temp]); + const auto& urr = urr_data_[i_temp]; // Determine the energy table int i_energy = 0; - while (true) { - if (E < urr->energy_(i_energy + 1)) {break;} - i_energy++; - } + while(E >= urr.energy_(i_energy + 1)) {++i_energy;}; // Sample the probability table using the cumulative distribution // Random nmbers for the xs calculation are sampled from a separate stream. - // This guarantees the rnadomness and, at the same time, makes sure we + // This guarantees the randomness and, at the same time, makes sure we // reuse random numbers for the same nuclide at different temperatures, // therefore preserving correlation of temperature in probability tables. prn_set_stream(STREAM_URR_PTABLE); @@ -432,15 +429,10 @@ void Nuclide::calculate_urr_xs(const int i_temp, const double E) prn_set_stream(STREAM_TRACKING); int i_low = 0; - while (true) { - if (urr->prob_(i_energy, URR_CUM_PROB, i_low) > r) {break;} - i_low++; - } + while (urr.prob_(i_energy, URR_CUM_PROB, i_low) <= r) {++i_low;}; + int i_up = 0; - while (true) { - if (urr->prob_(i_energy + 1, URR_CUM_PROB, i_up) > r) {break;} - i_up++; - } + while (urr.prob_(i_energy + 1, URR_CUM_PROB, i_up) <= r) {++i_up;}; // Determine elastic, fission, and capture cross sections from the // probability table @@ -448,51 +440,51 @@ void Nuclide::calculate_urr_xs(const int i_temp, const double E) double fission = 0.; double capture = 0.; double f; - if (urr->interp_ == Interpolation::lin_lin) { + if (urr.interp_ == Interpolation::lin_lin) { // Determine the interpolation factor on the table - f = (E - urr->energy_(i_energy)) / - (urr->energy_(i_energy + 1) - urr->energy_(i_energy)); + f = (E - urr.energy_(i_energy)) / + (urr.energy_(i_energy + 1) - urr.energy_(i_energy)); - elastic = (1. - f) * urr->prob_(i_energy, URR_ELASTIC, i_low) + - f * urr->prob_(i_energy + 1, URR_ELASTIC, i_up); - fission = (1. - f) * urr->prob_(i_energy, URR_FISSION, i_low) + - f * urr->prob_(i_energy + 1, URR_FISSION, i_up); - capture = (1. - f) * urr->prob_(i_energy, URR_N_GAMMA, i_low) + - f * urr->prob_(i_energy + 1, URR_N_GAMMA, i_up); - } else if (urr->interp_ == Interpolation::log_log) { + elastic = (1. - f) * urr.prob_(i_energy, URR_ELASTIC, i_low) + + f * urr.prob_(i_energy + 1, URR_ELASTIC, i_up); + fission = (1. - f) * urr.prob_(i_energy, URR_FISSION, i_low) + + f * urr.prob_(i_energy + 1, URR_FISSION, i_up); + capture = (1. - f) * urr.prob_(i_energy, URR_N_GAMMA, i_low) + + f * urr.prob_(i_energy + 1, URR_N_GAMMA, i_up); + } else if (urr.interp_ == Interpolation::log_log) { // Determine interpolation factor on the table - f = std::log(E / urr->energy_(i_energy)) / - std::log(urr->energy_(i_energy + 1) / urr->energy_(i_energy)); + f = std::log(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, URR_ELASTIC, i_low) > 0.) && - (urr->prob_(i_energy + 1, URR_ELASTIC, i_up) > 0.)) { + if ((urr.prob_(i_energy, URR_ELASTIC, i_low) > 0.) && + (urr.prob_(i_energy + 1, URR_ELASTIC, i_up) > 0.)) { elastic = std::exp((1. - f) * - std::log(urr->prob_(i_energy, URR_ELASTIC, i_low)) + - f * std::log(urr->prob_(i_energy + 1, URR_ELASTIC, i_up))); + std::log(urr.prob_(i_energy, URR_ELASTIC, i_low)) + + f * std::log(urr.prob_(i_energy + 1, URR_ELASTIC, i_up))); } else { elastic = 0.; } // Calculate the fission cross section/factor - if ((urr->prob_(i_energy, URR_FISSION, i_low) > 0.) && - (urr->prob_(i_energy + 1, URR_FISSION, i_up) > 0.)) { + if ((urr.prob_(i_energy, URR_FISSION, i_low) > 0.) && + (urr.prob_(i_energy + 1, URR_FISSION, i_up) > 0.)) { fission = std::exp((1. - f) * - std::log(urr->prob_(i_energy, URR_FISSION, i_low)) + - f * std::log(urr->prob_(i_energy + 1, URR_FISSION, i_up))); + std::log(urr.prob_(i_energy, URR_FISSION, i_low)) + + f * std::log(urr.prob_(i_energy + 1, URR_FISSION, i_up))); } else { fission = 0.; } // Calculate the capture cross section/factor - if ((urr->prob_(i_energy, URR_N_GAMMA, i_low) > 0.) && - (urr->prob_(i_energy + 1, URR_N_GAMMA, i_up) > 0.)) { + if ((urr.prob_(i_energy, URR_N_GAMMA, i_low) > 0.) && + (urr.prob_(i_energy + 1, URR_N_GAMMA, i_up) > 0.)) { capture = std::exp((1. - f) * - std::log(urr->prob_(i_energy, URR_N_GAMMA, i_low)) + - f * std::log(urr->prob_(i_energy + 1, URR_N_GAMMA, i_up))); + std::log(urr.prob_(i_energy, URR_N_GAMMA, i_low)) + + f * std::log(urr.prob_(i_energy + 1, URR_N_GAMMA, i_up))); } else { capture = 0.; } @@ -500,7 +492,7 @@ void Nuclide::calculate_urr_xs(const int i_temp, const double E) // Determine the treatment of inelastic scattering double inelastic = 0.; - if (urr->inelastic_flag_ != C_NONE) { + if (urr.inelastic_flag_ != C_NONE) { // get interpolation factor f = micro.interp_factor; @@ -514,7 +506,7 @@ void Nuclide::calculate_urr_xs(const int i_temp, const double E) } // Multiply by smooth cross-section if needed - if (urr->multiply_smooth_) { + if (urr.multiply_smooth_) { calculate_elastic_xs(); elastic *= micro.elastic; capture *= (micro.absorption - micro.fission); @@ -576,8 +568,7 @@ void set_micro_xs() } extern "C" void -nuclide_calculate_urr_xs(const bool use_mp, const int i_nuclide, - const int i_temp, const double E) +nuclide_calculate_urr_xs(bool use_mp, int i_nuclide, int i_temp, double E) { Nuclide* nuc = data::nuclides[i_nuclide - 1].get(); if (settings::urr_ptables_on && (nuc->urr_present_ && !use_mp)) { diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index b556df665..296a4a6a2 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -193,14 +193,13 @@ module nuclide_header type(C_PTR) :: path end function - subroutine nuclide_calculate_urr_xs_c(use_mp, i_nuclide, i_temp, E) & - bind(C, name='nuclide_calculate_urr_xs') + subroutine nuclide_calculate_urr_xs(use_mp, i_nuclide, i_temp, E) bind(C) import C_BOOL, C_INT, C_DOUBLE logical(C_BOOL), value, intent(in) :: use_mp integer(C_INT), value, intent(in) :: i_nuclide integer(C_INT), value, intent(in) :: i_temp real(C_DOUBLE), value, intent(in) :: E - end subroutine nuclide_calculate_urr_xs_c + end subroutine nuclide_calculate_urr_xs end interface contains @@ -939,7 +938,7 @@ contains ! If the particle is in the unresolved resonance range and there are ! probability tables, we need to determine cross sections from the table - call nuclide_calculate_urr_xs_c(use_mp, this % i_nuclide, i_temp, E) + call nuclide_calculate_urr_xs(use_mp, this % i_nuclide, i_temp, E) micro_xs % last_E = E micro_xs % last_sqrtkT = sqrtkT diff --git a/src/urr.cpp b/src/urr.cpp index 1ed7bca39..4a6643bfc 100644 --- a/src/urr.cpp +++ b/src/urr.cpp @@ -4,28 +4,12 @@ namespace openmc { -void -UrrData::from_hdf5(hid_t group_id) +UrrData::UrrData(hid_t group_id) { // Read interpolation and other flags int interp_temp; read_attribute(group_id, "interpolation", interp_temp); - switch (interp_temp) { - case static_cast(Interpolation::histogram): - interp_ = Interpolation::histogram; - break; - case static_cast(Interpolation::lin_lin): - interp_ = Interpolation::lin_lin; - break; - case static_cast(Interpolation::lin_log): - interp_ = Interpolation::lin_log; - break; - case static_cast(Interpolation::log_lin): - interp_ = Interpolation::log_lin; - break; - case static_cast(Interpolation::log_log): - interp_ = Interpolation::log_log; - } + interp_ = static_cast(interp_temp); // read the metadata read_attribute(group_id, "inelastic", inelastic_flag_);