From 5f3022989cd4b75c94e1315e5932815eeb11481d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 14 Aug 2018 22:21:41 -0500 Subject: [PATCH] Initial version of ThermalScattering class --- CMakeLists.txt | 1 + src/distribution.h | 4 + src/secondary_correlated.h | 24 ++- src/settings.cpp | 9 +- src/settings.h | 9 +- src/thermal.cpp | 353 +++++++++++++++++++++++++++++++++++++ src/thermal.h | 78 ++++++++ 7 files changed, 468 insertions(+), 10 deletions(-) create mode 100644 src/thermal.cpp create mode 100644 src/thermal.h diff --git a/CMakeLists.txt b/CMakeLists.txt index dc9a1231ba..efbffcb58c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -418,6 +418,7 @@ add_library(libopenmc SHARED src/state_point.cpp src/string_functions.cpp src/surface.cpp + src/thermal.cpp src/xml_interface.cpp src/xsdata.cpp) set_target_properties(libopenmc PROPERTIES diff --git a/src/distribution.h b/src/distribution.h index 84819e3f7c..815248a2e6 100644 --- a/src/distribution.h +++ b/src/distribution.h @@ -107,6 +107,10 @@ public: //! Sample a value from the distribution //! \return Sampled value double sample() const; + + // x property + std::vector& x() { return x_; } + const std::vector& x() const { return x_; } private: std::vector x_; //!< tabulated independent variable std::vector p_; //!< tabulated probability density diff --git a/src/secondary_correlated.h b/src/secondary_correlated.h index 3b6aae9c19..1e3eab0ab5 100644 --- a/src/secondary_correlated.h +++ b/src/secondary_correlated.h @@ -21,14 +21,6 @@ namespace openmc { class CorrelatedAngleEnergy : public AngleEnergy { public: - explicit CorrelatedAngleEnergy(hid_t group); - - //! Sample distribution for an angle and energy - //! \param[in] E_in Incoming energy in [eV] - //! \param[out] E_out Outgoing energy in [eV] - //! \param[out] mu Outgoing cosine with respect to current direction - void sample(double E_in, double& E_out, double& mu) const; -private: //! Outgoing energy/angle at a single incoming energy struct CorrTable { int n_discrete; //!< Number of discrete lines @@ -39,6 +31,22 @@ private: std::vector angle; //!< Angle distribution }; + explicit CorrelatedAngleEnergy(hid_t group); + + //! Sample distribution for an angle and energy + //! \param[in] E_in Incoming energy in [eV] + //! \param[out] E_out Outgoing energy in [eV] + //! \param[out] mu Outgoing cosine with respect to current direction + void sample(double E_in, double& E_out, double& mu) const; + + // energy property + std::vector& energy() { return energy_; } + const std::vector& energy() const { return energy_; } + + // distribution property + std::vector& distribution() { return distribution_; } + const std::vector& distribution() const { return distribution_; } +private: int n_region_; //!< Number of interpolation regions std::vector breakpoints_; //!< Breakpoints between regions std::vector interpolation_; //!< Interpolation laws diff --git a/src/settings.cpp b/src/settings.cpp index 8ecc04a790..8ca07a5276 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -1,5 +1,6 @@ #include "settings.h" +#include "constants.h" #include "error.h" #include "openmc.h" #include "string_utils.h" @@ -20,6 +21,12 @@ std::string path_multipole; std::string path_output; std::string path_source; +int temperature_method {TEMPERATURE_NEAREST}; +bool temperature_multipole {false}; +double temperature_tolerance {10.0}; +double temperature_default {293.6}; +std::array temperature_range {0.0, 0.0}; + //============================================================================== // Functions //============================================================================== @@ -67,4 +74,4 @@ void read_settings(pugi::xml_node* root) } } -} // namespace openmc \ No newline at end of file +} // namespace openmc diff --git a/src/settings.h b/src/settings.h index c203aafde7..ea4dfebdc1 100644 --- a/src/settings.h +++ b/src/settings.h @@ -4,6 +4,7 @@ //! \file settings.h //! \brief Settings for OpenMC +#include #include #include "pugixml.hpp" @@ -31,6 +32,12 @@ extern std::string path_multipole; extern std::string path_output; extern std::string path_source; +extern int temperature_method; +extern bool temperature_multipole; +extern double temperature_tolerance; +extern double temperature_default; +extern std::array temperature_range; + //============================================================================== //! Read settings from XML file //! \param[in] root XML node for @@ -40,4 +47,4 @@ extern "C" void read_settings(pugi::xml_node* root); } // namespace openmc -#endif // OPENMC_SETTINGS_H \ No newline at end of file +#endif // OPENMC_SETTINGS_H diff --git a/src/thermal.cpp b/src/thermal.cpp new file mode 100644 index 0000000000..b4818d0802 --- /dev/null +++ b/src/thermal.cpp @@ -0,0 +1,353 @@ +#include "thermal.h" + +#include // for sort, move +#include // for round +#include // for stringstream + +#include "xtensor/xarray.hpp" +#include "xtensor/xbuilder.hpp" +#include "xtensor/xmath.hpp" +#include "xtensor/xsort.hpp" +#include "xtensor/xtensor.hpp" +#include "xtensor/xview.hpp" + +#include "constants.h" +#include "error.h" +#include "random_lcg.h" +#include "search.h" +#include "secondary_correlated.h" +#include "settings.h" + +namespace openmc { + +//============================================================================== +// ThermalScattering implementation +//============================================================================== + +ThermalScattering::ThermalScattering(hid_t group, const std::vector& temperature, + int method, double tolerance, const double* minmax) +{ + // Get name of table from group + name_ = object_name(group); + + // Get rid of leading '/' + name_ = name_.substr(1); + + read_attribute(group, "atomic_weight_ratio", awr_); + read_attribute(group, "nuclides", nuclides_); + std::string sec_mode; + read_attribute(group, "secondary_mode", sec_mode); + if (sec_mode == "equal") { + secondary_mode_ = SAB_SECONDARY_EQUAL; + } else if (sec_mode == "skewed") { + secondary_mode_ = SAB_SECONDARY_SKEWED; + } else if (sec_mode == "continuous") { + secondary_mode_ = SAB_SECONDARY_CONT; + } + + // Read temperatures + hid_t kT_group = open_group(group, "kTs"); + + // Determine temperatures available + auto dset_names = dataset_names(kT_group); + auto n = dset_names.size(); + auto temps_available = xt::empty({n}); + for (int i = 0; i < dset_names.size(); ++i) { + // Read temperature value + double T; + read_dataset(kT_group, dset_names[i].data(), T); + temps_available[i] = T / K_BOLTZMANN; + } + std::sort(temps_available.begin(), temps_available.end()); + + // Determine actual temperatures to read -- start by checking whether a + // temperature range was given, in which case all temperatures in the range + // are loaded irrespective of what temperatures actually appear in the model + std::vector temps_to_read; + if (minmax[1] > 0.0) { + for (const auto& T : temps_available) { + if (minmax[0] <= T && T <= minmax[1]) { + temps_to_read.push_back(std::round(T)); + } + } + } + + switch (method) { + case TEMPERATURE_NEAREST: + // Determine actual temperatures to read + for (const auto& T : temperature) { + + auto i_closest = xt::argmin(xt::abs(temps_available - T))[0]; + auto temp_actual = temps_available[i_closest]; + if (std::fabs(temp_actual - T) < tolerance) { + if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temp_actual)) + == temps_to_read.end()) { + temps_to_read.push_back(std::round(temp_actual)); + } + } else { + std::stringstream msg; + msg << "Nuclear data library does not contain cross sections for " + << name_ << " at or near " << std::round(T) << " K."; + fatal_error(msg); + } + } + break; + + case TEMPERATURE_INTERPOLATION: + // If temperature interpolation or multipole is selected, get a list of + // bounding temperatures for each actual temperature present in the model + for (const auto& T : temperature) { + bool found = false; + for (int j = 0; j < temps_available.size() - 1; ++j) { + if (temps_available[j] <= T && T < temps_available[j + 1]) { + int T_j = std::round(temps_available[j]); + int T_j1 = std::round(temps_available[j + 1]); + if (std::find(temps_to_read.begin(), temps_to_read.end(), T_j) == temps_to_read.end()) { + temps_to_read.push_back(T_j); + } + if (std::find(temps_to_read.begin(), temps_to_read.end(), T_j1) == temps_to_read.end()) { + temps_to_read.push_back(T_j1); + } + found = true; + } + } + if (!found) { + std::stringstream msg; + msg << "Nuclear data library does not contain cross sections for " + << name_ << " at temperatures that bound " << std::round(T) << " K."; + fatal_error(msg); + } + } + } + + // Sort temperatures to read + std::sort(temps_to_read.begin(), temps_to_read.end()); + + auto n_temperature = temps_to_read.size(); + kTs_.reserve(n_temperature); + data_.reserve(n_temperature); + + for (auto T : temps_to_read) { + // Get temperature as a string + std::string temp_str = std::to_string(T) + "K"; + + // Read exact temperature value + double kT; + read_dataset(kT_group, temp_str.data(), kT); + kTs_.push_back(kT); + + // Open group for temperature i + hid_t T_group = open_group(group, temp_str.data()); + data_.emplace_back(T_group, secondary_mode_); + close_group(group); + } + + close_group(kT_group); +} + +void +ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, + double* elastic, double* inelastic) +{ + // Determine temperature for S(a,b) table + double kT = sqrtkT*sqrtkT; + int i; + if (temperature_method == TEMPERATURE_NEAREST) { + // If using nearest temperature, do linear search on temperature + for (i = 0; i < kTs_.size(); ++i) { + if (abs(kTs_[i] - kT) < K_BOLTZMANN*temperature_tolerance) { + break; + } + } + } else { + // Find temperatures that bound the actual temperature + for (i = 0; i < kTs_.size() - 1; ++i) { + if (kTs_[i] <= kT && kT < kTs_[i+1]) { + break; + } + } + + // Randomly sample between temperature i and i+1 + double f = (kT - kTs_[i]) / (kTs_[i+1] - kTs_[i]); + if (f > prn()) ++i; + } + + // Set temperature index + *i_temp = i; + + // Get pointer to S(a,b) table + auto& sab = data_[i]; + + // Get index and interpolation factor for inelastic grid + int i_grid; + double f; + if (E < sab.inelastic_e_in_.front()) { + i_grid = 0; + f = 0.0; + } else { + auto& E_in = sab.inelastic_e_in_; + i_grid = lower_bound_index(E_in.begin(), E_in.end(), E); + f = (E - E_in[i_grid]) / (E_in[i_grid+1] - E_in[i_grid]); + } + + // Calculate S(a,b) inelastic scattering cross section + auto& xs = sab.inelastic_sigma_; + *inelastic = (1.0 - f) * xs[i_grid] + f * xs[i_grid + 1]; + + // Check for elastic data + if (E < sab.threshold_elastic_) { + // Determine whether elastic scattering is given in the coherent or + // incoherent approximation. For coherent, the cross section is + // represented as P/E whereas for incoherent, it is simply P + + auto& E_in = sab.elastic_e_in_; + + if (sab.elastic_mode_ == SAB_ELASTIC_EXACT) { + if (E < E_in.front()) { + // If energy is below that of the lowest Bragg peak, the elastic + // cross section will be zero + *elastic = 0.0; + } else { + i_grid = lower_bound_index(E_in.begin(), E_in.end(), E); + *elastic = sab.elastic_P_[i_grid] / E; + } + } else { + // Determine index on elastic energy grid + if (E < E_in.front()) { + i_grid = 0; + } else { + i_grid = lower_bound_index(E_in.begin(), E_in.end(), E); + } + + // Get interpolation factor for elastic grid + f = (E - E_in[i_grid])/(E_in[i_grid+1] - E_in[i_grid]); + + // Calculate S(a,b) elastic scattering cross section + auto& xs = sab.elastic_P_; + *elastic = (1.0 - f) * xs[i_grid] + f * xs[i_grid + 1]; + } + } else { + // No elastic data + *elastic = 0.0; + } +} + +//============================================================================== +// ThermalData implementation +//============================================================================== + +ThermalData::ThermalData(hid_t group, int secondary_mode) +{ + // Coherent elastic data + if (object_exists(group, "elastic")) { + // Read cross section data + hid_t elastic_group = open_group(group, "elastic"); + + // Read elastic cross section + xt::xarray temp; + hid_t dset = open_dataset(elastic_group, "xs"); + read_dataset(dset, temp); + + // Get view on energies and cross section/probability values + auto E_in = xt::view(temp, 0); + auto P = xt::view(temp, 1); + + // Set cross section data and type + std::copy(E_in.begin(), E_in.end(), std::back_inserter(elastic_e_in_)); + std::copy(P.begin(), P.end(), std::back_inserter(elastic_P_)); + n_elastic_e_in_ = elastic_e_in_.size(); + + // Determine elastic type + std::string type; + read_attribute(dset, "type", type); + if (type == "tab1") { + elastic_mode_ = SAB_ELASTIC_DISCRETE; + } else if (type == "bragg") { + elastic_mode_ = SAB_ELASTIC_EXACT; + } + close_dataset(dset); + + // Set elastic threshold + threshold_elastic_ = elastic_e_in_.back(); + + // Read angle distribution + if (elastic_mode_ != SAB_ELASTIC_EXACT) { + xt::xarray mu_out; + read_dataset(elastic_group, "mu_out", mu_out); + elastic_mu_ = mu_out; + } + + close_group(elastic_group); + } + + // Inelastic data + if (object_exists(group, "inelastic")) { + // Read type of inelastic data + hid_t inelastic_group = open_group(group, "inelastic"); + + // Read cross section data + xt::xarray temp; + read_dataset(inelastic_group, "xs", temp); + + // Get view of inelastic cross section and energy grid + auto E_in = xt::view(temp, 0); + auto xs = xt::view(temp, 1); + + // Set cross section data + std::copy(E_in.begin(), E_in.end(), std::back_inserter(inelastic_e_in_)); + std::copy(xs.begin(), xs.end(), std::back_inserter(inelastic_sigma_)); + n_inelastic_e_in_ = inelastic_e_in_.size(); + + // Set inelastic threshold + threshold_inelastic_ = inelastic_e_in_.back(); + + if (secondary_mode != SAB_SECONDARY_CONT) { + // Read energy distribution + xt::xarray E_out; + read_dataset(inelastic_group, "energy_out", E_out); + inelastic_e_out_ = E_out; + + // Read angle distribution + xt::xarray mu_out; + read_dataset(inelastic_group, "mu_out", mu_out); + inelastic_mu_ = mu_out; + } else { + // Read correlated angle-energy distribution + CorrelatedAngleEnergy dist {inelastic_group}; + + // Convert to S(a,b) native format + for (const auto& edist : dist.distribution()) { + // Create temporary distribution + DistEnergySab d; + + // Copy outgoing energy distribution + d.n_e_out = edist.e_out.size(); + d.e_out = edist.e_out; + d.e_out_pdf = edist.p; + d.e_out_cdf = edist.c; + + for (int j = 0; j < d.n_e_out; ++j) { + auto adist = dynamic_cast(edist.angle[j].get()); + if (adist) { + // On first pass, allocate space for angles + if (j == 0) { + auto n_mu = adist->x().size(); + n_inelastic_mu_ = n_mu; + d.mu = xt::empty({d.n_e_out, n_mu}); + } + + // Copy outgoing angles + auto mu_j = xt::view(d.mu, j); + std::copy(adist->x().begin(), adist->x().end(), mu_j.begin()); + } + } + + inelastic_data_.push_back(std::move(d)); + } + } + + close_group(inelastic_group); + } +} + +} // namespace openmc diff --git a/src/thermal.h b/src/thermal.h new file mode 100644 index 0000000000..59f5c97ee7 --- /dev/null +++ b/src/thermal.h @@ -0,0 +1,78 @@ +#ifndef OPENMC_THERMAL_SCATTERING_H +#define OPENMC_THERMAL_SCATTERING_H + +#include +#include +#include + +#include "xtensor/xtensor.hpp" + +#include "hdf5_interface.h" + +namespace openmc { + +class ThermalData { +public: + ThermalData(hid_t group, int secondary_mode); +private: + struct DistEnergySab { + std::size_t n_e_out; + xt::xtensor e_out; + xt::xtensor e_out_pdf; + xt::xtensor e_out_cdf; + xt::xtensor mu; + }; + + // Threshold for thermal scattering treatment (usually ~4 eV) + double threshold_inelastic_; + double threshold_elastic_ {0.0}; + + // Inelastic scattering data + std::size_t n_inelastic_e_in_; // # of incoming E for inelastic + std::size_t n_inelastic_e_out_; // # of outgoing E for inelastic + std::size_t n_inelastic_mu_; // # of outgoing angles for inelastic + std::vector inelastic_e_in_; + std::vector inelastic_sigma_; + + // The following are used only if secondary_mode is 0 or 1 + xt::xtensor inelastic_e_out_; + xt::xtensor inelastic_mu_; + + // The following is used only if secondary_mode is 3 + // The different implementation is necessary because the continuous + // representation has a variable number of outgoing energy points for each + // incoming energy + std::vector inelastic_data_; // One for each Ein + + // Elastic scattering data + int elastic_mode_; // elastic mode (discrete/exact) + std::size_t n_elastic_e_in_; // # of incoming E for elastic + std::size_t n_elastic_mu_; // # of outgoing angles for elastic + std::vector elastic_e_in_; + std::vector elastic_P_; + xt::xtensor elastic_mu_; + + friend class ThermalScattering; +}; + +class ThermalScattering { +public: + ThermalScattering(hid_t group, const std::vector& temperature, int method, + double tolerance, const double* minmax); + + void calculate_xs(double E, double sqrtkT, int* i_temp, double* elastic, + double* inelastic); + + std::string name_; // name of table, e.g. "c_H_in_H2O" + double awr_; // weight of nucleus in neutron masses + std::vector kTs_; // temperatures in eV (k*T) + std::vector nuclides_; // List of valid nuclides + int secondary_mode_; // secondary mode (equal/skewed/continuous) + + // cross sections and distributions at each temperature + std::vector data_; +}; + +} // namespace openmc + +#endif // OPENMC_THERMAL_SCATTERING_H