From 12540cbc5b9a7e76715447c12dc7a40a546e70f0 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Wed, 5 Dec 2018 13:49:04 -0500 Subject: [PATCH 1/6] initial implementation of Urr in C++ --- CMakeLists.txt | 1 + include/openmc/nuclide.h | 6 ++++++ include/openmc/urr.h | 34 +++++++++++++++++++++++++++++++++ src/nuclide.cpp | 41 ++++++++++++++++++++++++++++++++++++++++ src/urr.cpp | 40 +++++++++++++++++++++++++++++++++++++++ 5 files changed, 122 insertions(+) create mode 100644 include/openmc/urr.h create mode 100644 src/urr.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 55f4d7f4ee..da55400c43 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -417,6 +417,7 @@ add_library(libopenmc SHARED src/position.cpp src/progress_bar.cpp src/pugixml/pugixml_c.cpp + src/urr.cpp src/random_lcg.cpp src/reaction.cpp src/reaction_product.cpp diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 9a99f71e82..72e15a02f1 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -14,6 +14,7 @@ #include "openmc/endf.h" #include "openmc/reaction.h" #include "openmc/reaction_product.h" +#include "urr.h" namespace openmc { @@ -68,6 +69,11 @@ public: std::vector elastic_0K_; std::vector xs_cdf_; + // Unresolved resonance range information + bool urr_present_ {false}; + int urr_inelastic_ {C_NONE}; + std::vector urr_data_; + std::vector> reactions_; //! Reactions std::vector index_inelastic_scatter_; diff --git a/include/openmc/urr.h b/include/openmc/urr.h new file mode 100644 index 0000000000..224c79a358 --- /dev/null +++ b/include/openmc/urr.h @@ -0,0 +1,34 @@ +//! \brief UrrData information for the unresolved resonance treatment + +#ifndef OPENMC_URR_H +#define OPENMC_URR_H + +#include "xtensor/xtensor.hpp" + +#include "openmc/hdf5_interface.h" + +namespace openmc { + +//============================================================================== +//! UrrData contains probability tables for the unresolved resonance range. +//============================================================================== + +class UrrData{ +public: + long unsigned int n_energy_; // # of incident energies + long unsigned int n_prob_; // # of probabilities + int interp_; // interpolation type (2=lin-lin, 5=log-log) + int inelastic_flag_; // inelastic competition flag + int absorption_flag_; // other absorption flag + bool multiply_smooth_; // multiply by smooth cross section? + xt::xtensor energy_; // incident energies + xt::xtensor prob_; // Actual probability tables + + //! Load the URR data from the provided HDF5 group + void + from_hdf5(hid_t group_id); +}; + +} // namespace openmc + +#endif // OPENMC_URR_H diff --git a/src/nuclide.cpp b/src/nuclide.cpp index ea1e213d57..2f00e612f3 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -191,6 +191,47 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n) } close_group(rxs_group); + // Read unresolved resonance probability tables if present + if (object_exists(group, "urr")) { + urr_present_ = true; + urr_data_.resize(temps_to_read.size()); + + for (int i = 0; i < temps_to_read.size(); i++) { + // Get temperature as a string + std::string temp_str {std::to_string(temps_to_read[i]) + "K"}; + + // 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); + close_group(urr_group); + + // Check for negative values + if (xt::any(urr_data_[i].prob_ < 0.) && mpi::master) { + warning("Negative value(s) found on probability table for nuclide " + + name_ + " at " + temp_str); + } + } + + // If the inelastic competition flag indicates that the inelastic cross + // section should be determined from a normal reaction cross section, we + // need to get the index of the reaction. + if (temps_to_read.size() > 0) { + if (urr_data_[0].inelastic_flag_ > 0) { + for (int i = 0; i < reactions_.size(); i++) { + if (reactions_[i]->mt_ == urr_data_[0].inelastic_flag_) { + urr_inelastic_ = i; + } + } + + // Abort if no corresponding inelastic reaction was found + if (urr_inelastic_ == C_NONE) { + fatal_error("Could no find inelastic reaction specified on " + "unresolved resonance probability table."); + } + } + } + } + // Check for nu-total if (object_exists(group, "total_nu")) { // Read total nu data diff --git a/src/urr.cpp b/src/urr.cpp new file mode 100644 index 0000000000..7a33a20e8f --- /dev/null +++ b/src/urr.cpp @@ -0,0 +1,40 @@ +#include "openmc/urr.h" + +namespace openmc { + +void +UrrData::from_hdf5(hid_t group_id) +{ + // Read interpolation and other flags + read_attribute(group_id, "interpolation", interp_); + read_attribute(group_id, "inelastic", inelastic_flag_); + read_attribute(group_id, "absorption", absorption_flag_); + int i; + read_attribute(group_id, "multiply_smooth", i); + multiply_smooth_ = (i == 1); + + // read the enrgies at which tables exist + hid_t dset = open_dataset(group_id, "energy"); + hsize_t dims[1]; + get_shape(dset, dims); + close_dataset(dset); + n_energy_ = static_cast(dims[0]); + energy_ = xt::xtensor({n_energy_}, 0.); + read_dataset_as_shape(group_id, "energy", energy_); + + // Read URR tables + dset = open_dataset(group_id, "table"); + hsize_t dims3[3]; + get_shape(dset, dims3); + close_dataset(dset); + n_prob_ = static_cast(dims3[0]); + xt::xarray temp({n_energy_, 6, n_prob_}); + read_dataset(group_id, "table", temp); + + prob_ = xt::xtensor({n_energy_, 6, n_prob_}, 0.); + prob_ = temp; + //TODO: swap 1st and last indices? + +} + +} \ No newline at end of file From 5eab6dbed0d2d3d0f564781c51f5c41a6382fc1e Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Wed, 5 Dec 2018 21:02:59 -0500 Subject: [PATCH 2/6] Got URR working, now just need a few simplifications of the Fortran code --- include/openmc/constants.h | 14 ++-- include/openmc/nuclide.h | 121 ++++++++++++++-------------- include/openmc/urr.h | 7 +- src/nuclide.cpp | 153 +++++++++++++++++++++++++++++++++++- src/nuclide_header.F90 | 156 +++---------------------------------- src/physics.cpp | 2 +- src/urr.cpp | 43 ++++++---- 7 files changed, 262 insertions(+), 234 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 2064252037..dbb5ed3147 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -257,12 +257,12 @@ constexpr int LIBRARY_PHOTON {3}; constexpr int LIBRARY_MULTIGROUP {4}; // Probability table parameters -constexpr int URR_CUM_PROB {1}; -constexpr int URR_TOTAL {2}; -constexpr int URR_ELASTIC {3}; -constexpr int URR_FISSION {4}; -constexpr int URR_N_GAMMA {5}; -constexpr int URR_HEATING {6}; +constexpr int URR_CUM_PROB {0}; +constexpr int URR_TOTAL {1}; +constexpr int URR_ELASTIC {2}; +constexpr int URR_FISSION {3}; +constexpr int URR_N_GAMMA {4}; +constexpr int URR_HEATING {5}; // Maximum number of partial fission reactions constexpr int PARTIAL_FISSION_MAX {4}; @@ -424,7 +424,7 @@ constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE // Interpolation rules enum class Interpolation { - histogram, lin_lin, lin_log, log_lin, log_log + histogram = 1, lin_lin = 2, lin_log = 3, log_lin = 4, log_log = 5 }; // Run modes diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 72e15a02f1..662a6fa7c2 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -24,63 +24,6 @@ namespace openmc { constexpr double CACHE_INVALID {-1.0}; -//=============================================================================== -// Data for a nuclide -//=============================================================================== - -class Nuclide { -public: - // Types, aliases - using EmissionMode = ReactionProduct::EmissionMode; - struct EnergyGrid { - std::vector grid_index; - std::vector energy; - }; - - // Constructors - Nuclide(hid_t group, const double* temperature, int n); - - // Methods - double nu(double E, EmissionMode mode, int group=0) const; - void calculate_elastic_xs(int i_nuclide) const; - - //! Determines the microscopic 0K elastic cross section at a trial relative - //! energy used in resonance scattering - double elastic_xs_0K(double E) const; - - // Data members - std::string name_; //! Name of nuclide, e.g. "U235" - int Z_; //! Atomic number - int A_; //! Mass number - int metastable_; //! Metastable state - double awr_; //! Atomic weight ratio - std::vector kTs_; //! temperatures in eV (k*T) - std::vector grid_; //! Energy grid at each temperature - - bool fissionable_ {false}; //! Whether nuclide is fissionable - bool has_partial_fission_ {false}; //! has partial fission reactions? - std::vector fission_rx_; //! Fission reactions - int n_precursor_ {0}; //! Number of delayed neutron precursors - std::unique_ptr total_nu_; //! Total neutron yield - - // Resonance scattering information - bool resonant_ {false}; - std::vector energy_0K_; - std::vector elastic_0K_; - std::vector xs_cdf_; - - // Unresolved resonance range information - bool urr_present_ {false}; - int urr_inelastic_ {C_NONE}; - std::vector urr_data_; - - std::vector> reactions_; //! Reactions - std::vector index_inelastic_scatter_; - -private: - void create_derived(); -}; - //=============================================================================== //! Cached microscopic cross sections for a particular nuclide at the current //! energy @@ -138,6 +81,68 @@ struct MaterialMacroXS { double pair_production; //!< macroscopic pair production xs }; +//=============================================================================== +// Data for a nuclide +//=============================================================================== + +class Nuclide { +public: + // Types, aliases + using EmissionMode = ReactionProduct::EmissionMode; + struct EnergyGrid { + std::vector grid_index; + std::vector energy; + }; + + // Constructors + Nuclide(hid_t group, const double* temperature, int n, int i_nuclide); + + // Methods + double nu(double E, EmissionMode mode, int group=0) const; + void calculate_elastic_xs() const; + + //! Determines the microscopic 0K elastic cross section at a trial relative + //! energy used in resonance scattering + double elastic_xs_0K(double E) const; + + //! \brief Determines cross sections in the unresolved resonance range + //! from probability tables. + void calculate_urr_xs(const int i_temp, const double E); + + // Data members + std::string name_; //! Name of nuclide, e.g. "U235" + int Z_; //! Atomic number + int A_; //! Mass number + int metastable_; //! Metastable state + double awr_; //! Atomic weight ratio + std::vector kTs_; //! temperatures in eV (k*T) + std::vector grid_; //! Energy grid at each temperature + int i_nuclide_; //! Index in the nuclides array + + bool fissionable_ {false}; //! Whether nuclide is fissionable + bool has_partial_fission_ {false}; //! has partial fission reactions? + std::vector fission_rx_; //! Fission reactions + int n_precursor_ {0}; //! Number of delayed neutron precursors + std::unique_ptr total_nu_; //! Total neutron yield + + // Resonance scattering information + bool resonant_ {false}; + std::vector energy_0K_; + std::vector elastic_0K_; + std::vector xs_cdf_; + + // Unresolved resonance range information + bool urr_present_ {false}; + int urr_inelastic_ {C_NONE}; + std::vector urr_data_; + + std::vector> reactions_; //! Reactions + std::vector index_inelastic_scatter_; + +private: + void create_derived(); +}; + //============================================================================== // Global variables //============================================================================== @@ -170,6 +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 int i_nuclide, + const int i_temp, const double E); } // namespace openmc diff --git a/include/openmc/urr.h b/include/openmc/urr.h index 224c79a358..acfe4d0590 100644 --- a/include/openmc/urr.h +++ b/include/openmc/urr.h @@ -5,6 +5,7 @@ #include "xtensor/xtensor.hpp" +#include "openmc/constants.h" #include "openmc/hdf5_interface.h" namespace openmc { @@ -15,16 +16,14 @@ namespace openmc { class UrrData{ public: - long unsigned int n_energy_; // # of incident energies - long unsigned int n_prob_; // # of probabilities - int interp_; // interpolation type (2=lin-lin, 5=log-log) + Interpolation interp_; // interpolation type int inelastic_flag_; // inelastic competition flag int absorption_flag_; // other absorption flag bool multiply_smooth_; // multiply by smooth cross section? xt::xtensor energy_; // incident energies xt::xtensor prob_; // Actual probability tables - //! Load the URR data from the provided HDF5 group + //! \brief Load the URR data from the provided HDF5 group void from_hdf5(hid_t group_id); }; diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 2f00e612f3..64c983fd39 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -5,6 +5,7 @@ #include "openmc/error.h" #include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" +#include "openmc/random_lcg.h" #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/string_utils.h" @@ -33,7 +34,7 @@ MaterialMacroXS material_xs; // Nuclide implementation //============================================================================== -Nuclide::Nuclide(hid_t group, const double* temperature, int n) +Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide) { // Get name of nuclide from group, removing leading '/' name_ = object_name(group).substr(1); @@ -42,6 +43,7 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n) read_attribute(group, "A", A_); read_attribute(group, "metastable", metastable_); read_attribute(group, "atomic_weight_ratio", awr_); + i_nuclide_ = i_nuclide; // Determine temperatures available hid_t kT_group = open_group(group, "kTs"); @@ -365,10 +367,10 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const } } -void Nuclide::calculate_elastic_xs(int i_nuclide) const +void Nuclide::calculate_elastic_xs() const { // Get temperature index, grid index, and interpolation factor - auto& micro = simulation::micro_xs[i_nuclide]; + auto& micro = simulation::micro_xs[i_nuclide_]; int i_temp = micro.index_temp - 1; int i_grid = micro.index_grid - 1; double f = micro.interp_factor; @@ -402,6 +404,142 @@ 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) +{ + auto& micro = simulation::micro_xs[i_nuclide_]; + micro.use_ptable = true; + + // Create a shorthand for the URR data + UrrData* 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++; + } + + // 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 + // reuse random numbers for the same nuclide at different temperatures, + // therefore preserving correlation of temperature in probability tables. + prn_set_stream(STREAM_URR_PTABLE); + //TODO: to maintain the same random number stream as the Fortran code this + //replaces, the seed is set with i_nuclide_ + 1 instead of i_nuclide_ + double r = future_prn(static_cast(i_nuclide_ + 1)); + 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++; + } + int i_up = 0; + while (true) { + if (urr->prob_(i_energy + 1, URR_CUM_PROB, i_up) > r) {break;} + i_up++; + } + + // Determine elastic, fission, and capture cross sections from the + // probability table + double elastic = 0.; + double fission = 0.; + double capture = 0.; + double f; + 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)); + + 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)); + + // 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.)) { + 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))); + } 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.)) { + 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))); + } 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.)) { + 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))); + } else { + capture = 0.; + } + } + + // Determine the treatment of inelastic scattering + double inelastic = 0.; + if (urr->inelastic_flag_ != C_NONE) { + // get interpolation factor + f = micro.interp_factor; + + // Determine inelastic scattering cross section + Reaction* rx = reactions_[urr_inelastic_].get(); + int xs_index = micro.index_grid - rx->xs_[i_temp].threshold; + if (xs_index >= 0) { + inelastic = (1. - f) * rx->xs_[i_temp].value[xs_index] + + f * rx->xs_[i_temp].value[xs_index + 1]; + } + } + + // Multiply by smooth cross-section if needed + if (urr->multiply_smooth_) { + calculate_elastic_xs(); + elastic *= micro.elastic; + capture *= (micro.absorption - micro.fission); + fission *= micro.fission; + } + + // Check for negative values + if (elastic < 0.) {elastic = 0.;} + if (fission < 0.) {fission = 0.;} + if (capture < 0.) {capture = 0.;} + + // Set elastic, absorption, fission, and total x/s. Note that the total x/s + // is calculated as a sum of partials instead of the table-provided value + micro.elastic = elastic; + micro.absorption = capture + fission; + micro.fission = fission; + micro.total = elastic + inelastic + capture + fission; + + // Determine nu-fission cross-section + if (fissionable_) { + micro.nu_fission = nu(E, EmissionMode::total) * micro.fission; + } + +} + //============================================================================== // Fortran compatibility functions //============================================================================== @@ -415,7 +553,8 @@ set_particle_energy_bounds(int particle, double E_min, double E_max) extern "C" Nuclide* nuclide_from_hdf5_c(hid_t group, const double* temperature, int n) { - data::nuclides.push_back(std::make_unique(group, temperature, n)); + data::nuclides.push_back(std::make_unique(group, temperature, n, + data::nuclides.size())); return data::nuclides.back().get(); } @@ -436,4 +575,10 @@ void set_micro_xs() } } +extern "C" void +nuclide_calculate_urr_xs(const int i_nuclide, const int i_temp, const double E) +{ + data::nuclides[i_nuclide - 1]->calculate_urr_xs(i_temp - 1, E); +} + } // namespace openmc diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index becb865240..e05a703130 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -198,6 +198,14 @@ module nuclide_header character(kind=C_CHAR), intent(in) :: name(*) type(C_PTR) :: path end function + + subroutine nuclide_calculate_urr_xs_c(i_nuclide, i_temp, E) & + bind(C, name='nuclide_calculate_urr_xs') + import C_INT, C_DOUBLE + 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 interface contains @@ -980,7 +988,7 @@ contains if (urr_ptables_on .and. this % urr_present .and. .not. use_mp) then if (E > this % urr_data(i_temp) % energy(1) .and. E < this % & urr_data(i_temp) % energy(this % urr_data(i_temp) % n_energy)) then - call calculate_urr_xs(this, i_temp, E, micro_xs) + call nuclide_calculate_urr_xs_c(this % i_nuclide, i_temp, E) end if end if @@ -1246,152 +1254,6 @@ contains end if end subroutine multipole_deriv_eval -!=============================================================================== -! CALCULATE_URR_XS determines cross sections in the unresolved resonance range -! from probability tables -!=============================================================================== - - subroutine calculate_urr_xs(this, i_temp, E, micro_xs) - class(Nuclide), intent(in) :: this ! Nuclide object - integer, intent(in) :: i_temp ! temperature index - real(8), intent(in) :: E ! energy - type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache - - integer :: i_energy ! index for energy - integer :: i_low ! band index at lower bounding energy - integer :: i_up ! band index at upper bounding energy - integer :: threshold ! threshold energy index - real(8) :: f ! interpolation factor - real(8) :: r ! pseudo-random number - real(8) :: elastic ! elastic cross section - real(8) :: capture ! (n,gamma) cross section - real(8) :: fission ! fission cross section - real(8) :: inelastic ! inelastic cross section - - micro_xs % use_ptable = .true. - - associate (urr => this % urr_data(i_temp)) - ! determine energy table - i_energy = 1 - do - if (E < urr % energy(i_energy + 1)) exit - i_energy = i_energy + 1 - end do - - ! determine interpolation factor on table - f = (E - urr % energy(i_energy)) / & - (urr % energy(i_energy + 1) - urr % energy(i_energy)) - - ! sample probability table using the cumulative distribution - - ! Random numbers for xs calculation are sampled from a separated stream. - ! This guarantees the randomness and, at the same time, makes sure we reuse - ! random number for the same nuclide at different temperatures, therefore - ! preserving correlation of temperature in probability tables. - call prn_set_stream(STREAM_URR_PTABLE) - r = future_prn(int(this % i_nuclide, 8)) - call prn_set_stream(STREAM_TRACKING) - - i_low = 1 - do - if (urr % prob(i_energy, URR_CUM_PROB, i_low) > r) exit - i_low = i_low + 1 - end do - i_up = 1 - do - if (urr % prob(i_energy + 1, URR_CUM_PROB, i_up) > r) exit - i_up = i_up + 1 - end do - - ! determine elastic, fission, and capture cross sections from probability - ! table - if (urr % interp == LINEAR_LINEAR) then - elastic = (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_low) + & - f * urr % prob(i_energy + 1, URR_ELASTIC, i_up) - fission = (ONE - f) * urr % prob(i_energy, URR_FISSION, i_low) + & - f * urr % prob(i_energy + 1, URR_FISSION, i_up) - capture = (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_low) + & - f * urr % prob(i_energy + 1, URR_N_GAMMA, i_up) - elseif (urr % interp == LOG_LOG) then - ! Get logarithmic interpolation factor - f = log(E / urr % energy(i_energy)) / & - log(urr % energy(i_energy + 1) / urr % energy(i_energy)) - - ! Calculate elastic cross section/factor - elastic = ZERO - if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then - elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, & - i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, & - i_up))) - end if - - ! Calculate fission cross section/factor - fission = ZERO - if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then - fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, & - i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, & - i_up))) - end if - - ! Calculate capture cross section/factor - capture = ZERO - if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then - capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, & - i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, & - i_up))) - end if - end if - - ! Determine treatment of inelastic scattering - inelastic = ZERO - if (urr % inelastic_flag > 0) then - ! Get index on energy grid and interpolation factor - i_energy = micro_xs % index_grid - f = micro_xs % interp_factor - - ! Determine inelastic scattering cross section - associate (rx => this % reactions(this % urr_inelastic)) - threshold = rx % xs_threshold(i_temp) - if (i_energy >= threshold) then - inelastic = (ONE - f) * rx % xs(i_temp, i_energy - threshold + 1) + & - f * rx % xs(i_temp, i_energy - threshold + 2) - end if - end associate - end if - - ! Multiply by smooth cross-section if needed - if (urr % multiply_smooth) then - call this % calculate_elastic_xs(micro_xs) - elastic = elastic * micro_xs % elastic - capture = capture * (micro_xs % absorption - micro_xs % fission) - fission = fission * micro_xs % fission - end if - - ! Check for negative values - if (elastic < ZERO) elastic = ZERO - if (fission < ZERO) fission = ZERO - if (capture < ZERO) capture = ZERO - - ! Set elastic, absorption, fission, and total cross sections. Note that the - ! total cross section is calculated as sum of partials rather than using the - ! table-provided value - micro_xs % elastic = elastic - micro_xs % absorption = capture + fission - micro_xs % fission = fission - micro_xs % total = elastic + inelastic + capture + fission - - ! Determine nu-fission cross section - if (this % fissionable) then - micro_xs % nu_fission = this % nu(E, EMISSION_TOTAL) * & - micro_xs % fission - end if - end associate - - end subroutine calculate_urr_xs - !=============================================================================== ! CHECK_DATA_VERSION checks for the right version of nuclear data within HDF5 ! files diff --git a/src/physics.cpp b/src/physics.cpp index 43fc2871b4..df2f519134 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -627,7 +627,7 @@ void scatter(Particle* p, int i_nuclide, int i_nuc_mat) // Calculate elastic cross section if it wasn't precalculated if (micro.elastic == CACHE_INVALID) { - nuc->calculate_elastic_xs(i_nuclide); + nuc->calculate_elastic_xs(); } double prob = micro.elastic - micro.thermal; diff --git a/src/urr.cpp b/src/urr.cpp index 7a33a20e8f..1ecc8ca8c6 100644 --- a/src/urr.cpp +++ b/src/urr.cpp @@ -1,25 +1,45 @@ #include "openmc/urr.h" +#include + namespace openmc { void UrrData::from_hdf5(hid_t group_id) { // Read interpolation and other flags - read_attribute(group_id, "interpolation", interp_); + 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; + } + read_attribute(group_id, "inelastic", inelastic_flag_); read_attribute(group_id, "absorption", absorption_flag_); - int i; - read_attribute(group_id, "multiply_smooth", i); - multiply_smooth_ = (i == 1); + int temp_multiply_smooth; + read_attribute(group_id, "multiply_smooth", temp_multiply_smooth); + multiply_smooth_ = (temp_multiply_smooth == 1); // read the enrgies at which tables exist hid_t dset = open_dataset(group_id, "energy"); hsize_t dims[1]; get_shape(dset, dims); close_dataset(dset); - n_energy_ = static_cast(dims[0]); - energy_ = xt::xtensor({n_energy_}, 0.); + // n_energy_ = static_cast(dims[0]); + energy_ = xt::xtensor({dims[0]}, 0.); read_dataset_as_shape(group_id, "energy", energy_); // Read URR tables @@ -27,14 +47,9 @@ UrrData::from_hdf5(hid_t group_id) hsize_t dims3[3]; get_shape(dset, dims3); close_dataset(dset); - n_prob_ = static_cast(dims3[0]); - xt::xarray temp({n_energy_, 6, n_prob_}); - read_dataset(group_id, "table", temp); - - prob_ = xt::xtensor({n_energy_, 6, n_prob_}, 0.); - prob_ = temp; - //TODO: swap 1st and last indices? - + xt::xarray temp_arr({dims3[0], dims3[1], dims3[2]}, 0.); + read_dataset(group_id, "table", temp_arr); + prob_ = temp_arr; } } \ No newline at end of file From d4c8c3d34e4f377ca25eb14bdbc062a309f497d5 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Thu, 6 Dec 2018 19:28:26 -0500 Subject: [PATCH 3/6] completed conversion of Urr to C++ --- CMakeLists.txt | 1 - include/openmc/nuclide.h | 3 +- include/openmc/urr.h | 1 + src/nuclide.cpp | 12 +++++-- src/nuclide_header.F90 | 73 ++++++---------------------------------- src/urr.cpp | 3 +- src/urr_header.F90 | 72 --------------------------------------- 7 files changed, 25 insertions(+), 140 deletions(-) delete mode 100644 src/urr_header.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index da55400c43..eb2224867b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -361,7 +361,6 @@ add_library(libopenmc SHARED src/timer_header.F90 src/tracking.F90 src/track_output.F90 - src/urr_header.F90 src/vector_header.F90 src/volume_calc.F90 src/volume_header.F90 diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 662a6fa7c2..67444a106f 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -175,10 +175,9 @@ 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 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); - } // namespace openmc #endif // OPENMC_NUCLIDE_H diff --git a/include/openmc/urr.h b/include/openmc/urr.h index acfe4d0590..36ea073324 100644 --- a/include/openmc/urr.h +++ b/include/openmc/urr.h @@ -20,6 +20,7 @@ public: 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 diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 64c983fd39..727968a4d8 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -576,9 +576,17 @@ void set_micro_xs() } extern "C" void -nuclide_calculate_urr_xs(const int i_nuclide, const int i_temp, const double E) +nuclide_calculate_urr_xs(const bool use_mp, const int i_nuclide, + const int i_temp, const double E) { - data::nuclides[i_nuclide - 1]->calculate_urr_xs(i_temp - 1, E); + Nuclide* nuc = data::nuclides[i_nuclide - 1].get(); + if (settings::urr_ptables_on && (nuc->urr_present_ && !use_mp)) { + if ((E > nuc->urr_data_[i_temp - 1].energy_(0)) && + (E < nuc->urr_data_[i_temp - 1].energy_( + nuc->urr_data_[i_temp - 1].n_energy_ - 1))) { + nuc->calculate_urr_xs(i_temp - 1, E); + } + } } } // namespace openmc diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index e05a703130..b556df665b 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -23,7 +23,6 @@ module nuclide_header use stl_vector, only: VectorInt, VectorReal use string use simulation_header, only: need_depletion_rx - use urr_header, only: UrrData implicit none @@ -76,11 +75,6 @@ module nuclide_header integer, allocatable :: index_fission(:) ! indices in reactions class(Function1D), allocatable :: total_nu - ! Unresolved resonance data - logical :: urr_present = .false. - integer :: urr_inelastic - type(UrrData), allocatable :: urr_data(:) - ! Multipole data logical :: mp_present = .false. type(MultipoleArray), pointer :: multipole => null() @@ -199,12 +193,13 @@ module nuclide_header type(C_PTR) :: path end function - subroutine nuclide_calculate_urr_xs_c(i_nuclide, i_temp, E) & + subroutine nuclide_calculate_urr_xs_c(use_mp, i_nuclide, i_temp, E) & bind(C, name='nuclide_calculate_urr_xs') - import C_INT, C_DOUBLE - integer(C_INT), value, intent(in) :: i_nuclide - integer(C_INT), value, intent(in) :: i_temp - real(C_DOUBLE), value, intent(in) :: E + 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 interface @@ -261,7 +256,7 @@ contains integer :: i integer :: i_closest integer :: n_temperature - integer(HID_T) :: urr_group, nu_group + integer(HID_T) :: nu_group integer(HID_T) :: energy_group, energy_dset integer(HID_T) :: kT_group integer(HID_T) :: rxs_group @@ -446,47 +441,6 @@ contains end do call close_group(rxs_group) - ! Read unresolved resonance probability tables if present - if (object_exists(group_id, 'urr')) then - this % urr_present = .true. - allocate(this % urr_data(n_temperature)) - - do i = 1, n_temperature - ! Get temperature as a string - temp_str = trim(to_str(temps_to_read % data(i))) // "K" - - ! Read probability tables for i-th temperature - urr_group = open_group(group_id, 'urr/' // trim(temp_str)) - call this % urr_data(i) % from_hdf5(urr_group) - call close_group(urr_group) - - ! Check for negative values - if (any(this % urr_data(i) % prob < ZERO) .and. master) then - call warning("Negative value(s) found on probability table & - &for nuclide " // this % name // " at " // trim(temp_str)) - end if - end do - - ! if the inelastic competition flag indicates that the inelastic cross - ! section should be determined from a normal reaction cross section, we - ! need to get the index of the reaction - if (n_temperature > 0) then - if (this % urr_data(1) % inelastic_flag > 0) then - do i = 1, size(this % reactions) - if (this % reactions(i) % MT == this % urr_data(1) % inelastic_flag) then - this % urr_inelastic = i - end if - end do - - ! Abort if no corresponding inelastic reaction was found - if (this % urr_inelastic == NONE) then - call fatal_error("Could not find inelastic reaction specified on & - &unresolved resonance probability table.") - end if - end if - end if - end if - ! Check for nu-total if (object_exists(group_id, 'total_nu')) then nu_group = open_group(group_id, 'total_nu') @@ -784,7 +738,7 @@ contains real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache - logical :: use_mp ! true if XS can be calculated with windowed multipole + logical(C_BOOL) :: use_mp ! true if XS can be calculated with windowed multipole integer :: i_temp ! index for temperature integer :: i_grid ! index on nuclide energy grid integer :: i_low ! lower logarithmic mapping index @@ -839,7 +793,7 @@ contains ! 1. physics.F90 - scatter - For inelastic scatter. ! 2. physics.F90 - sample_fission - For partial fissions. ! 3. tally.F90 - score_general - For tallying on MTxxx reactions. - ! 4. nuclide_header.F90 - calculate_urr_xs - For unresolved purposes. + ! 4. nuclide.h - calculate_urr_xs - For unresolved purposes. ! It is worth noting that none of these occur in the resolved ! resonance range, so the value here does not matter. index_temp is ! set to -1 to force a segfault in case a developer messes up and tries @@ -982,15 +936,10 @@ contains call calculate_sab_xs(this, i_sab, E, sqrtkT, sab_frac, micro_xs) end if + ! If the particle is in the unresolved resonance range and there are ! probability tables, we need to determine cross sections from the table - - if (urr_ptables_on .and. this % urr_present .and. .not. use_mp) then - if (E > this % urr_data(i_temp) % energy(1) .and. E < this % & - urr_data(i_temp) % energy(this % urr_data(i_temp) % n_energy)) then - call nuclide_calculate_urr_xs_c(this % i_nuclide, i_temp, E) - end if - end if + call nuclide_calculate_urr_xs_c(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 1ecc8ca8c6..1ed7bca394 100644 --- a/src/urr.cpp +++ b/src/urr.cpp @@ -27,6 +27,7 @@ UrrData::from_hdf5(hid_t group_id) interp_ = Interpolation::log_log; } + // read the metadata read_attribute(group_id, "inelastic", inelastic_flag_); read_attribute(group_id, "absorption", absorption_flag_); int temp_multiply_smooth; @@ -38,7 +39,7 @@ UrrData::from_hdf5(hid_t group_id) hsize_t dims[1]; get_shape(dset, dims); close_dataset(dset); - // n_energy_ = static_cast(dims[0]); + n_energy_ = static_cast(dims[0]); energy_ = xt::xtensor({dims[0]}, 0.); read_dataset_as_shape(group_id, "energy", energy_); diff --git a/src/urr_header.F90 b/src/urr_header.F90 deleted file mode 100644 index 0ec9e89deb..0000000000 --- a/src/urr_header.F90 +++ /dev/null @@ -1,72 +0,0 @@ -module urr_header - - use hdf5_interface, only: read_attribute, open_dataset, read_dataset, & - close_dataset, get_shape, HID_T, HSIZE_T - - implicit none - -!=============================================================================== -! URRDATA contains probability tables for the unresolved resonance range. -!=============================================================================== - - type UrrData - integer :: n_energy ! # of incident neutron energies - integer :: n_prob ! # of probabilities - integer :: interp ! inteprolation (2=lin-lin, 5=log-log) - integer :: inelastic_flag ! inelastic competition flag - integer :: absorption_flag ! other absorption flag - logical :: multiply_smooth ! multiply by smooth cross section? - real(8), allocatable :: energy(:) ! incident energies - real(8), allocatable :: prob(:,:,:) ! actual probabibility tables - contains - procedure :: from_hdf5 => urr_from_hdf5 - end type UrrData - -contains - - subroutine urr_from_hdf5(this, group_id) - class(UrrData), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - - integer :: i, j, k - integer(HID_T) :: energy - integer(HID_T) :: table - integer(HSIZE_T) :: dims(1) - integer(HSIZE_T) :: dims3(3) - real(8), allocatable :: temp(:,:,:) - - ! Read interpolation and other flags - call read_attribute(this % interp, group_id, 'interpolation') - call read_attribute(this % inelastic_flag, group_id, 'inelastic') - call read_attribute(this % absorption_flag, group_id, 'absorption') - call read_attribute(i, group_id, 'multiply_smooth') - this % multiply_smooth = (i == 1) - - ! Read energies at which tables exist - energy = open_dataset(group_id, 'energy') - call get_shape(energy, dims) - this % n_energy = int(dims(1), 4) - allocate(this % energy(this % n_energy)) - call read_dataset(this % energy, energy) - call close_dataset(energy) - - ! Read URR tables - table = open_dataset(group_id, 'table') - call get_shape(table, dims3) - this % n_prob = int(dims3(1), 4) - allocate(temp(this % n_prob, 6, this % n_energy)) - call read_dataset(temp, table) - call close_dataset(table) - - ! Swap first and last indices - allocate(this % prob(this % n_energy, 6, this % n_prob)) - do i = 1, this % n_energy - do j = 1, 6 - do k = 1, this % n_prob - this % prob(i, j, k) = temp(k, j, i) - end do - end do - end do - end subroutine urr_from_hdf5 - -end module urr_header From 999415be6a840200290cff9b5e561670e0c7b1a2 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Thu, 6 Dec 2018 19:40:10 -0500 Subject: [PATCH 4/6] removed a todo from physics_mg --- include/openmc/physics_mg.h | 16 ++++----------- src/physics_mg.cpp | 40 ++++++++++++++++++++----------------- src/tracking.F90 | 7 +++---- 3 files changed, 29 insertions(+), 34 deletions(-) diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index e5648c4a52..8febad19a5 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -10,24 +10,18 @@ namespace openmc { -//TODO: Remove material_xs parameters when they reside on -// the C-side this should happen after materials, physics, input, and tallies -// are brought over - //! \brief samples particle behavior after a collision event. //! \param p Particle to operate on -//! \param material_xs The cross section cache for the current material extern "C" void -collision_mg(Particle* p, const MaterialMacroXS* material_xs); +collision_mg(Particle* p); //! \brief samples a reaction type. //! //! Note that there is special logic when suvival biasing is turned on since //! fission and disappearance are treated implicitly. //! \param p Particle to operate on -//! \param material_xs The cross section cache for the current material void -sample_reaction(Particle* p, const MaterialMacroXS* material_xs); +sample_reaction(Particle* p); //! \brief Samples the scattering event //! \param p Particle to operate on @@ -40,16 +34,14 @@ scatter(Particle* p); //! \param bank_array The particle bank to populate //! \param size_bank Number of particles currently in the bank //! \param bank_array_size Allocated size of the bank -//! \param material_xs The cross section cache for the current material void create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, - int64_t bank_array_size, const MaterialMacroXS* material_xs); + int64_t bank_array_size); //! \brief Handles an absorption event //! \param p Particle to operate on -//! \param material_xs The cross section cache for the current material void -absorption(Particle* p, const MaterialMacroXS* material_xs); +absorption(Particle* p); } // namespace openmc #endif // OPENMC_PHYSICS_MG_H diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 202e5abe30..1d475aaa74 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -22,13 +22,13 @@ namespace openmc { void -collision_mg(Particle* p, const MaterialMacroXS* material_xs) +collision_mg(Particle* p) { // Add to the collision counter for the particle p->n_collision++; // Sample the reaction type - sample_reaction(p, material_xs); + sample_reaction(p); // Display information about collision if ((settings::verbosity >= 10) || (simulation::trace)) { @@ -39,7 +39,7 @@ collision_mg(Particle* p, const MaterialMacroXS* material_xs) } void -sample_reaction(Particle* p, const MaterialMacroXS* material_xs) +sample_reaction(Particle* p) { // Create fission bank sites. Note that while a fission reaction is sampled, // it never actually "happens", i.e. the weight of the particle does not @@ -48,19 +48,20 @@ sample_reaction(Particle* p, const MaterialMacroXS* material_xs) if (model::materials[p->material - 1]->fissionable) { if (settings::run_mode == RUN_MODE_EIGENVALUE) { - create_fission_sites(p, simulation::fission_bank.data(), &simulation::n_bank, - simulation::fission_bank.size(), material_xs); + create_fission_sites( + p, simulation::fission_bank.data(), &simulation::n_bank, + simulation::fission_bank.size()); } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && (settings::create_fission_neutrons)) { create_fission_sites(p, p->secondary_bank, &(p->n_secondary), - MAX_SECONDARY, material_xs); + MAX_SECONDARY); } } // If survival biasing is being used, the following subroutine adjusts the // weight of the particle. Otherwise, it checks to see if absorption occurs. - if (material_xs->absorption > 0.) { - absorption(p, material_xs); + if (simulation::material_xs.absorption > 0.) { + absorption(p); } else { p->absorb_wgt = 0.; } @@ -102,7 +103,7 @@ scatter(Particle* p) void create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, - int64_t bank_array_size, const MaterialMacroXS* material_xs) + int64_t bank_array_size) { // TODO: Heat generation from fission @@ -111,8 +112,8 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p->wgt / simulation::keff * weight * material_xs->nu_fission / - material_xs->total; + double nu_t = p->wgt / simulation::keff * weight * + simulation::material_xs.nu_fission / simulation::material_xs.total; // Sample the number of neutrons produced int nu = static_cast(nu_t); @@ -200,11 +201,12 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, } void -absorption(Particle* p, const MaterialMacroXS* material_xs) +absorption(Particle* p) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p->absorb_wgt = p->wgt * material_xs->absorption / material_xs->total; + p->absorb_wgt = p->wgt * + simulation::material_xs.absorption / simulation::material_xs.total; // Adjust weight of particle by the probability of absorption p->wgt -= p->absorb_wgt; @@ -212,13 +214,15 @@ absorption(Particle* p, const MaterialMacroXS* material_xs) // Score implicit absorpion estimate of keff #pragma omp atomic - global_tally_absorption += p->absorb_wgt * material_xs->nu_fission / - material_xs->absorption; + global_tally_absorption += p->absorb_wgt * + simulation::material_xs.nu_fission / + simulation::material_xs.absorption; } else { - if (material_xs->absorption > prn() * material_xs->total) { + if (simulation::material_xs.absorption > + prn() * simulation::material_xs.total) { #pragma omp atomic - global_tally_absorption += p->wgt * material_xs->nu_fission / - material_xs->absorption; + global_tally_absorption += p->wgt * simulation::material_xs.nu_fission / + simulation::material_xs.absorption; p->alive = false; p->event = EVENT_ABSORB; } diff --git a/src/tracking.F90 b/src/tracking.F90 index 48db0cf6df..d345d7cc0b 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -33,10 +33,9 @@ module tracking implicit none interface - subroutine collision_mg(p, material_xs) bind(C) - import Particle, C_DOUBLE, MaterialMacroXS + subroutine collision_mg(p) bind(C) + import Particle, C_DOUBLE type(Particle), intent(inout) :: p - type(MaterialMacroXS), intent(in) :: material_xs end subroutine collision_mg end interface @@ -234,7 +233,7 @@ contains if (run_CE) then call collision(p) else - call collision_mg(p, material_xs) + call collision_mg(p) end if ! Score collision estimator tallies -- this is done after a collision From af54b60c2eaeebe914a13f502cc20f04795e4d66 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Fri, 7 Dec 2018 20:38:24 -0500 Subject: [PATCH 5/6] Fixing @paulromano comments --- include/openmc/hdf5_interface.h | 1 - include/openmc/nuclide.h | 6 +-- include/openmc/urr.h | 17 ++++--- src/nuclide.cpp | 81 +++++++++++++++------------------ src/nuclide_header.F90 | 7 ++- src/urr.cpp | 20 +------- 6 files changed, 52 insertions(+), 80 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 313b841e20..0c940342c7 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 67444a106f..455dab5fc7 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 36ea073324..6764254425 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 727968a4d8..f7ef87b56c 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 b556df665b..296a4a6a27 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 1ed7bca394..4a6643bfc8 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_); From 3cfd9dff9db5d54252192bf49d48e2affe9adf53 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Mon, 10 Dec 2018 20:01:03 -0500 Subject: [PATCH 6/6] Added xtensor read_dataset --- include/openmc/hdf5_interface.h | 28 ++++++++++++++++++++++++++++ src/urr.cpp | 21 ++++++--------------- 2 files changed, 34 insertions(+), 15 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 0c940342c7..5b5d2a9991 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -271,6 +271,34 @@ void read_dataset(hid_t obj_id, const char* name, xt::xarray& arr, bool indep close_dataset(dset); } + +template +void read_dataset(hid_t obj_id, const char* name, xt::xtensor& arr, bool indep=false) +{ + // Open dataset and read array + hid_t dset = open_dataset(obj_id, name); + + // Get shape of dataset + std::vector hsize_t_shape = object_shape(dset); + close_dataset(dset); + + // cast from hsize_t to size_t + std::vector shape(hsize_t_shape.size()); + for (int i = 0; i < shape.size(); i++) { + shape[i] = static_cast(hsize_t_shape[i]); + } + + // Allocate new xarray to read data into + xt::xarray xarr(shape); + + // Read data from the dataset + read_dataset(obj_id, name, xarr); + + // Copy into xtensor + arr = xarr; + +} + template void read_dataset_as_shape(hid_t obj_id, const char* name, xt::xtensor& arr, bool indep=false) diff --git a/src/urr.cpp b/src/urr.cpp index 4a6643bfc8..28ad2505f0 100644 --- a/src/urr.cpp +++ b/src/urr.cpp @@ -18,23 +18,14 @@ UrrData::UrrData(hid_t group_id) read_attribute(group_id, "multiply_smooth", temp_multiply_smooth); multiply_smooth_ = (temp_multiply_smooth == 1); - // read the enrgies at which tables exist - hid_t dset = open_dataset(group_id, "energy"); - hsize_t dims[1]; - get_shape(dset, dims); - close_dataset(dset); - n_energy_ = static_cast(dims[0]); - energy_ = xt::xtensor({dims[0]}, 0.); - read_dataset_as_shape(group_id, "energy", energy_); + // read the energies at which tables exist + read_dataset(group_id, "energy", energy_); + + // Set n_energy_ + n_energy_ = energy_.shape()[0]; // Read URR tables - dset = open_dataset(group_id, "table"); - hsize_t dims3[3]; - get_shape(dset, dims3); - close_dataset(dset); - xt::xarray temp_arr({dims3[0], dims3[1], dims3[2]}, 0.); - read_dataset(group_id, "table", temp_arr); - prob_ = temp_arr; + read_dataset(group_id, "table", prob_); } } \ No newline at end of file