From 8be229ee624b6254a98027c65b50fe395d41b323 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 20 Dec 2018 14:50:53 -0600 Subject: [PATCH] Start converting full Nuclide object to C++ (calculate_xs, init_grid) --- include/openmc/endf.h | 2 + include/openmc/nuclide.h | 54 +++++-- src/endf.cpp | 19 +++ src/nuclide.cpp | 338 +++++++++++++++++++++++++++++++++++++-- 4 files changed, 386 insertions(+), 27 deletions(-) diff --git a/include/openmc/endf.h b/include/openmc/endf.h index a7030d0355..170a647d50 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -4,6 +4,7 @@ #ifndef OPENMC_ENDF_H #define OPENMC_ENDF_H +#include #include #include "hdf5.h" @@ -98,6 +99,7 @@ private: std::vector factors_; //!< Partial sums of structure factors [eV-b] }; +std::unique_ptr read_function(hid_t group, const char* name); } // namespace openmc diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 455dab5fc7..1ee9ceba09 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -97,6 +97,17 @@ public: // Constructors Nuclide(hid_t group, const double* temperature, int n, int i_nuclide); + //! Initialize logarithmic grid for energy searches + //! \param E_min Minimum energy in [eV] + //! \param E_max Maximum energy in [eV] + //! \param M Number of equally log-spaced bins + void init_grid(double E_min, double E_max, int M); + + void calculate_xs(int i_sab, double E, int i_log_union, + double sqrtkT, double sab_frac); + + void calculate_sab_xs(int i_sab, double E, double sqrtkT, double sab_frac); + // Methods double nu(double E, EmissionMode mode, int group=0) const; void calculate_elastic_xs() const; @@ -107,23 +118,29 @@ public: //! \brief Determines cross sections in the unresolved resonance range //! from probability tables. - void calculate_urr_xs(int i_temp, double E); + void calculate_urr_xs(int i_temp, 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 - int i_nuclide_; //! Index in the nuclides array + 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 + 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 + // Temperature dependent cross section data + std::vector kTs_; //!< temperatures in eV (k*T) + std::vector grid_; //!< Energy grid at each temperature + std::vector> xs_; //!< Cross sections at each temperature + + // Fission data + 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 + std::unique_ptr fission_q_prompt_; //!< Prompt fission energy release + std::unique_ptr fission_q_recov_; //!< Recoverable fission energy release // Resonance scattering information bool resonant_ {false}; @@ -136,11 +153,18 @@ public: int urr_inelastic_ {C_NONE}; std::vector urr_data_; - std::vector> reactions_; //! Reactions + std::vector> reactions_; //!< Reactions + std::array reaction_index_; //!< Index of each reaction std::vector index_inelastic_scatter_; private: void create_derived(); + + static constexpr int XS_TOTAL {0}; + static constexpr int XS_ABSORPTION {1}; + static constexpr int XS_FISSION {2}; + static constexpr int XS_NU_FISSION {3}; + static constexpr int XS_PHOTON_PROD {4}; }; //============================================================================== diff --git a/src/endf.cpp b/src/endf.cpp index db1332591e..86f08b8a86 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -76,6 +76,25 @@ bool is_inelastic_scatter(int mt) } } +std::unique_ptr +read_function(hid_t group, const char* name) +{ + hid_t dset = open_dataset(group, name); + std::string func_type; + read_attribute(dset, "type", func_type); + std::unique_ptr func; + if (func_type == "Tabulated1D") { + func = std::make_unique(dset); + } else if (func_type == "Polynomial") { + func = std::make_unique(dset); + } else { + throw std::runtime_error{"Unknown function type " + func_type + + " for dataset " + object_name(dset)}; + } + close_dataset(dset); + return func; +} + //============================================================================== // Polynomial implementation //============================================================================== diff --git a/src/nuclide.cpp b/src/nuclide.cpp index f7ef87b56c..188847677f 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -8,8 +8,12 @@ #include "openmc/random_lcg.h" #include "openmc/search.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/string_utils.h" +#include "xtensor/xbuilder.hpp" +#include "xtensor/xview.hpp" + #include // for sort #include // for to_string, stoi @@ -234,36 +238,74 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide) } } - // Check for nu-total + // Check for total nu data if (object_exists(group, "total_nu")) { // Read total nu data hid_t nu_group = open_group(group, "total_nu"); - hid_t nu_dset = open_dataset(nu_group, "yield"); - std::string func_type; - read_attribute(nu_dset, "type", func_type); - if (func_type == "Tabulated1D") { - total_nu_ = std::make_unique(nu_dset); - } else if (func_type == "Polynomial") { - total_nu_ = std::make_unique(nu_dset); - } - close_dataset(nu_dset); + total_nu_ = read_function(nu_group, "yield"); close_group(nu_group); } + // Read fission energy release data if present + if (object_exists(group, "fission_energy_release")) { + hid_t fer_group = open_group(group, "fission_energy_release"); + fission_q_prompt_ = read_function(fer_group, "q_prompt"); + fission_q_recov_ = read_function(fer_group, "q_recoverable"); + close_group(fer_group); + } + this->create_derived(); } void Nuclide::create_derived() { + for (const auto& grid : grid_) { + // Allocate and initialize cross section + std::array shape {grid.energy.size(), 5}; + xs_.emplace_back(shape, 0.0); + } + + reaction_index_.fill(-1); for (int i = 0; i < reactions_.size(); ++i) { const auto& rx {reactions_[i]}; + // Set entry in direct address table for reaction + reaction_index_[rx->mt_] = i; + for (int t = 0; t < kTs_.size(); ++t) { + // TODO: off-by-one + int j = rx->xs_[t].threshold - 1; + int n = rx->xs_[t].value.size(); + auto xs = xt::adapt(rx->xs_[t].value); + + for (const auto& p : rx->products_) { + if (p.particle_ == ParticleType::photon) { + auto pprod = xt::view(xs_[t], xt::range(j, j+n), XS_PHOTON_PROD); + for (int k = 0; k < n; ++k) { + double E = grid_[t].energy[k+j]; + pprod[k] += xs[k] * (*p.yield_)(E); + } + } + } + // Skip redundant reactions if (rx->redundant_) continue; + // Add contribution to total cross section + auto total = xt::view(xs_[t], xt::range(j,j+n), XS_TOTAL); + total += xs; + + // Add contribution to absorption cross section + auto absorption = xt::view(xs_[t], xt::range(j,j+n), XS_ABSORPTION); + if (is_disappearance(rx->mt_)) { + absorption += xs; + } + if (is_fission(rx->mt_)) { fissionable_ = true; + auto fission = xt::view(xs_[t], xt::range(j,j+n), XS_FISSION); + fission += xs; + absorption += xs; // Keep track of fission reactions if (t == 0) { @@ -283,6 +325,18 @@ void Nuclide::create_derived() } } + // Calculate nu-fission cross section + for (int t = 0; t < kTs_.size(); ++t) { + if (fissionable_) { + int n = grid_[t].energy.size(); + for (int i = 0; i < n; ++i) { + double E = grid_[t].energy[i]; + xs_[t](i, XS_NU_FISSION) = nu(E, EmissionMode::total) + * xs_[t](i, XS_FISSION); + } + } + } + if (settings::res_scat_on) { // Determine if this nuclide should be treated as a resonant scatterer if (!settings::res_scat_nuclides.empty()) { @@ -327,6 +381,33 @@ void Nuclide::create_derived() } } +void Nuclide::init_grid(double E_min, double E_max, int M) +{ + // Determine equal-logarithmic energy spacing + double spacing = std::log(E_max/E_min)/M; + + // Create equally log-spaced energy grid + auto umesh = xt::linspace(0.0, M*spacing, M+1); + + for (auto& grid : grid_) { + // Resize array for storing grid indices + grid.grid_index.resize(M + 1); + + // Determine corresponding indices in nuclide grid to energies on + // equal-logarithmic grid + int j = 0; + for (int k = 0; k <= M; ++k) { + while (std::log(grid.energy[j]/E_min) <= umesh(k)) { + // Ensure that for isotopes where maxval(grid.energy) << E_max that + // there are no out-of-bounds issues. + if (j == grid.energy.size()) break; + ++j; + } + grid.grid_index[k] = j; + } + } +} + double Nuclide::nu(double E, EmissionMode mode, int group) const { if (!fissionable_) return 0.0; @@ -404,7 +485,240 @@ 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(int i_temp, double E) +void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, + double sqrtkT, double sab_frac) +{ + auto& micro_xs = simulation::micro_xs[i_nuclide_]; + + // Initialize cached cross sections to zero + micro_xs.elastic = CACHE_INVALID; + micro_xs.thermal = 0.0; + micro_xs.thermal_elastic = 0.0; + + // Check to see if there is multipole data present at this energy + bool use_mp = false; + if (nuclide_wmp_present(i_nuclide_)) { + use_mp = (E >= nuclide_wmp_emin(i_nuclide_) && + E <= nuclide_wmp_emax(i_nuclide_)); + } + + int i_temp = -1; + + // Evaluate multipole or interpolate + if (use_mp) { + // Call multipole kernel + double sig_s, sig_a, sig_f; + //multipole_eval(this % multipole, E, sqrtkT, sig_s, sig_a, sig_f) + + micro_xs.total = sig_s + sig_a; + micro_xs.elastic = sig_s; + micro_xs.absorption = sig_a; + micro_xs.fission = sig_f; + micro_xs.nu_fission = fissionable_ ? + sig_f * this->nu(E, EmissionMode::total) : 0.0; + + if (simulation::need_depletion_rx) { + // Only non-zero reaction is (n,gamma) + micro_xs.reaction[0] = sig_a - sig_f; + + // Set all other reaction cross sections to zero + for (int i = 1; i < DEPLETION_RX.size(); ++i) { + micro_xs.reaction[i] = 0.0; + } + } + + // Ensure these values are set + // Note, the only time either is used is in one of 4 places: + // 1. physics.cpp - scatter - For inelastic scatter. + // 2. physics.cpp - sample_fission - For partial fissions. + // 3. tally.F90 - score_general - For tallying on MTxxx reactions. + // 4. nuclide.cpp - 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 + // to use it with multipole. + micro_xs.index_temp = -1; + micro_xs.index_grid = 0; + micro_xs.interp_factor = 0.0; + + } else { + // Find the appropriate temperature index. + double kT = sqrtkT*sqrtkT; + double f; + switch (settings::temperature_method) { + case TEMPERATURE_NEAREST: + { + double max_diff = INFTY; + for (int t = 0; t < kTs_.size(); ++t) { + if (std::abs(kTs_[t] - kT) < max_diff) { + i_temp = t; + } + } + } + break; + + case TEMPERATURE_INTERPOLATION: + // Find temperatures that bound the actual temperature + for (i_temp = 0; i_temp < kTs_.size() - 1; ++i_temp) { + if (kTs_[i_temp] <= kT && kT < kTs_[i_temp + 1]) break; + } + + // Randomly sample between temperature i and i+1 + f = (kT - kTs_[i_temp]) / (kTs_[i_temp + 1] - kTs_[i_temp]); + if (f > prn()) ++i_temp; + break; + } + + // Determine the energy grid index using a logarithmic mapping to + // reduce the energy range over which a binary search needs to be + // performed + + const auto& grid {grid_[i_temp]}; + const auto& xs {xs_[i_temp]}; + + int i_grid; + if (E < grid.energy.front()) { + i_grid = 0; + } else if (E > grid.energy.back()) { + i_grid = grid.energy.size() - 2; + } else { + // Determine bounding indices based on which equal log-spaced + // interval the energy is in + int i_low = grid.grid_index[i_log_union]; + int i_high = grid.grid_index[i_log_union + 1] + 1; + + // Perform binary search over reduced range + i_grid = lower_bound_index(&grid.energy[i_low], &grid.energy[i_high], E) - 1; + } + + // check for rare case where two energy points are the same + if (grid.energy[i_grid] == grid.energy[i_grid + 1]) ++i_grid; + + // calculate interpolation factor + f = (E - grid.energy[i_grid]) / + (grid.energy[i_grid + 1]- grid.energy[i_grid]); + + micro_xs.index_temp = i_temp; + micro_xs.index_grid = i_grid; + micro_xs.interp_factor = f; + + // Calculate microscopic nuclide total cross section + micro_xs.total = (1.0 - f)*xs(i_grid,XS_TOTAL) + + f*xs(i_grid + 1,XS_TOTAL); + + // Calculate microscopic nuclide absorption cross section + micro_xs.absorption = (1.0 - f)*xs(i_grid,XS_ABSORPTION) + + f*xs(i_grid + 1,XS_ABSORPTION); + + if (fissionable_) { + // Calculate microscopic nuclide total cross section + micro_xs.fission = (1.0 - f)*xs(i_grid,XS_FISSION) + + f*xs(i_grid + 1,XS_FISSION); + + // Calculate microscopic nuclide nu-fission cross section + micro_xs.nu_fission = (1.0 - f)*xs(i_grid,XS_NU_FISSION) + + f*xs(i_grid + 1,XS_NU_FISSION); + } else { + micro_xs.fission = 0.0; + micro_xs.nu_fission = 0.0; + } + + // Calculate microscopic nuclide photon production cross section + micro_xs.photon_prod = (1.0 - f)*xs(i_grid,XS_PHOTON_PROD) + + f*xs(i_grid + 1,XS_PHOTON_PROD); + + // Depletion-related reactions + if (simulation::need_depletion_rx) { + // Initialize all reaction cross sections to zero + for (double& xs_i : micro_xs.reaction) { + xs_i = 0.0; + } + + for (int j = 0; j < DEPLETION_RX.size(); ++j) { + // If reaction is present and energy is greater than threshold, set the + // reaction xs appropriately + int i_rx = reaction_index_[DEPLETION_RX[j]]; + + const auto& rx = reactions_[i_rx]; + const auto& rx_xs = rx->xs_[i_temp].value; + + if (i_rx >= 0) { + // Physics says that (n,gamma) is not a threshold reaction, so we don't + // need to specifically check its threshold index + if (j == 0) { + micro_xs.reaction[0] = (1.0 - f)*rx_xs[i_grid] + + f*rx_xs[i_grid + 1]; + continue; + } + + // TODO: off-by-one + int threshold = rx->xs_[i_temp].threshold - 1; + if (i_grid >= threshold) { + micro_xs.reaction[j] = (1.0 - f)*rx_xs[i_grid - threshold] + + f*rx_xs[i_grid - threshold + 1]; + } else if (j >= 3) { + // One can show that the the threshold for (n,(x+1)n) is always + // higher than the threshold for (n,xn). Thus, if we are below + // the threshold for, e.g., (n,2n), there is no reason to check + // the threshold for (n,3n) and (n,4n). + break; + } + } + } + } + } + + // Initialize sab treatment to false + micro_xs.index_sab = -1; + micro_xs.sab_frac = 0.0; + + // Initialize URR probability table treatment to false + micro_xs.use_ptable = false; + + // If there is S(a,b) data for this nuclide, we need to set the sab_scatter + // and sab_elastic cross sections and correct the total and elastic cross + // sections. + + if (i_sab >= 0) this->calculate_sab_xs(i_sab, E, sqrtkT, sab_frac); + + // If the particle is in the unresolved resonance range and there are + // probability tables, we need to determine cross sections from the table + this->calculate_urr_xs(i_temp, E); + + micro_xs.last_E = E; + micro_xs.last_sqrtkT = sqrtkT; +} + +void Nuclide::calculate_sab_xs(int i_sab, double E, double sqrtkT, double sab_frac) +{ + auto& micro {simulation::micro_xs[i_nuclide_]}; + + // Set flag that S(a,b) treatment should be used for scattering + micro.index_sab = i_sab; + + // Calculate the S(a,b) cross section + int i_temp; + double elastic; + double inelastic; + //data::sab_tables[i_sab]->calculate_xs(E, sqrtkT, &i_temp, &elastic, &inelastic); + + // Store the S(a,b) cross sections. + micro.thermal = sab_frac * (elastic + inelastic); + micro.thermal_elastic = sab_frac * elastic; + + // Calculate free atom elastic cross section + this->calculate_elastic_xs(); + + // Correct total and elastic cross sections + micro.total = micro.total + micro.thermal - sab_frac*micro.elastic; + micro.elastic = micro.thermal + (1.0 - sab_frac)*micro.elastic; + + // Save temperature index and thermal fraction + micro.index_temp_sab = i_temp; + micro.sab_frac = sab_frac; +} + +void Nuclide::calculate_urr_xs(int i_temp, double E) const { auto& micro = simulation::micro_xs[i_nuclide_]; micro.use_ptable = true; @@ -414,7 +728,7 @@ void Nuclide::calculate_urr_xs(int i_temp, double E) // Determine the energy table int i_energy = 0; - while(E >= urr.energy_(i_energy + 1)) {++i_energy;}; + while (E >= urr.energy_(i_energy + 1)) {++i_energy;}; // Sample the probability table using the cumulative distribution