diff --git a/CMakeLists.txt b/CMakeLists.txt index 55f4d7f4ee..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 @@ -417,6 +416,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/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/hdf5_interface.h b/include/openmc/hdf5_interface.h index 313b841e20..5b5d2a9991 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -272,6 +272,33 @@ void read_dataset(hid_t obj_id, const char* name, xt::xarray& arr, bool indep } +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/include/openmc/nuclide.h b/include/openmc/nuclide.h index 9a99f71e82..455dab5fc7 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 { @@ -23,58 +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_; - - 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 @@ -132,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(int i_temp, 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 //============================================================================== @@ -164,7 +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(bool use_mp, int i_nuclide, + int i_temp, double E); } // namespace openmc 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/include/openmc/urr.h b/include/openmc/urr.h new file mode 100644 index 0000000000..6764254425 --- /dev/null +++ b/include/openmc/urr.h @@ -0,0 +1,33 @@ +//! \brief UrrData information for the unresolved resonance treatment + +#ifndef OPENMC_URR_H +#define OPENMC_URR_H + +#include "xtensor/xtensor.hpp" + +#include "openmc/constants.h" +#include "openmc/hdf5_interface.h" + +namespace openmc { + +//============================================================================== +//! UrrData contains probability tables for the unresolved resonance range. +//============================================================================== + +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 + + //! \brief Load the URR data from the provided HDF5 group + explicit UrrData(hid_t group_id); +}; + +} // namespace openmc + +#endif // OPENMC_URR_H diff --git a/src/nuclide.cpp b/src/nuclide.cpp index ea1e213d57..f7ef87b56c 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"); @@ -191,6 +193,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_.reserve(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_.emplace_back(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 @@ -324,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; @@ -361,6 +404,134 @@ 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) +{ + auto& micro = simulation::micro_xs[i_nuclide_]; + micro.use_ptable = true; + + // Create a shorthand for the URR data + const auto& urr = urr_data_[i_temp]; + + // Determine the energy table + int i_energy = 0; + 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 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); + //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 (urr.prob_(i_energy, URR_CUM_PROB, i_low) <= r) {++i_low;}; + + int i_up = 0; + 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 + 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 //============================================================================== @@ -374,7 +545,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(); } @@ -395,4 +567,17 @@ void set_micro_xs() } } +extern "C" void +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)) { + 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 becb865240..296a4a6a27 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() @@ -198,6 +192,14 @@ module nuclide_header character(kind=C_CHAR), intent(in) :: name(*) type(C_PTR) :: path end function + + 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 end interface contains @@ -253,7 +255,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 @@ -438,47 +440,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') @@ -776,7 +737,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 @@ -831,7 +792,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 @@ -974,15 +935,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 calculate_urr_xs(this, i_temp, E, micro_xs) - end if - end if + call nuclide_calculate_urr_xs(use_mp, this % i_nuclide, i_temp, E) micro_xs % last_E = E micro_xs % last_sqrtkT = sqrtkT @@ -1246,152 +1202,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 4d8c01cf35..27d0e7b81f 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/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 diff --git a/src/urr.cpp b/src/urr.cpp new file mode 100644 index 0000000000..28ad2505f0 --- /dev/null +++ b/src/urr.cpp @@ -0,0 +1,31 @@ +#include "openmc/urr.h" + +#include + +namespace openmc { + +UrrData::UrrData(hid_t group_id) +{ + // Read interpolation and other flags + int interp_temp; + read_attribute(group_id, "interpolation", interp_temp); + interp_ = static_cast(interp_temp); + + // read the metadata + read_attribute(group_id, "inelastic", inelastic_flag_); + read_attribute(group_id, "absorption", absorption_flag_); + int temp_multiply_smooth; + read_attribute(group_id, "multiply_smooth", temp_multiply_smooth); + multiply_smooth_ = (temp_multiply_smooth == 1); + + // read the energies at which tables exist + read_dataset(group_id, "energy", energy_); + + // Set n_energy_ + n_energy_ = energy_.shape()[0]; + + // Read URR tables + read_dataset(group_id, "table", prob_); +} + +} \ No newline at end of file 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