From 71e36db418a61525525cda64617a9f13d8ee3f0d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 16 Nov 2018 16:22:22 -0600 Subject: [PATCH] Store reactions as a vector on Nuclide --- include/openmc/nuclide.h | 8 ++- include/openmc/reaction.h | 2 - src/nuclide.cpp | 138 ++++++++++++++++++++++++++++++++++++-- src/nuclide_header.F90 | 32 +++++---- src/reaction.cpp | 8 --- src/reaction_header.F90 | 31 +++------ 6 files changed, 164 insertions(+), 55 deletions(-) diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index c86ba6334..bc3c2a1cc 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -5,11 +5,13 @@ #define OPENMC_NUCLIDE_H #include +#include // for unique_ptr #include #include #include "openmc/constants.h" +#include "openmc/reaction.h" namespace openmc { @@ -20,7 +22,7 @@ namespace openmc { class Nuclide { public: // Constructors - Nuclide(hid_t group); + Nuclide(hid_t group, const double* temperature, int n); // Data members std::string name_; //! Name of nuclide, e.g. "U235" @@ -28,7 +30,7 @@ public: int A_; //! Mass number int metastable_; //! Metastable state double awr_; //! Atomic weight ratio - + std::vector> reactions_; //! Reactions }; //============================================================================== @@ -42,7 +44,7 @@ namespace data { extern std::array energy_min; extern std::array energy_max; -extern std::vector nuclides; +extern std::vector> nuclides; } // namespace data diff --git a/include/openmc/reaction.h b/include/openmc/reaction.h index d44babd90..7f20e7a95 100644 --- a/include/openmc/reaction.h +++ b/include/openmc/reaction.h @@ -51,8 +51,6 @@ std::string reaction_name(int mt); //============================================================================== extern "C" { - Reaction* reaction_from_hdf5(hid_t group, int* temperatures, int n); - void reaction_delete(Reaction* rx); int reaction_mt(Reaction* rx); double reaction_q_value(Reaction* rx); bool reaction_scatter_in_cm(Reaction* rx); diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 82022159f..629cd86f5 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1,6 +1,14 @@ #include "openmc/nuclide.h" +#include "openmc/container_util.h" +#include "openmc/error.h" #include "openmc/hdf5_interface.h" +#include "openmc/message_passing.h" +#include "openmc/settings.h" +#include "openmc/string_utils.h" + +#include // for sort +#include // for to_string, stoi namespace openmc { @@ -13,7 +21,7 @@ namespace data { std::array energy_min {0.0, 0.0}; std::array energy_max {INFTY, INFTY}; -std::vector nuclides; +std::vector> nuclides; } // namespace data @@ -21,7 +29,7 @@ std::vector nuclides; // Nuclide implementation //============================================================================== -Nuclide::Nuclide(hid_t group) +Nuclide::Nuclide(hid_t group, const double* temperature, int n) { // Get name of nuclide from group, removing leading '/' name_ = object_name(group).substr(1); @@ -30,6 +38,120 @@ Nuclide::Nuclide(hid_t group) read_attribute(group, "A", A_); read_attribute(group, "metastable", metastable_); read_attribute(group, "atomic_weight_ratio", awr_); + + // Determine temperatures available + hid_t kT_group = open_group(group, "kTs"); + auto dset_names = dataset_names(kT_group); + std::vector temps_available; + for (const auto& name : dset_names) { + double T; + read_dataset(kT_group, name.c_str(), T); + temps_available.push_back(T / K_BOLTZMANN); + } + std::sort(temps_available.begin(), temps_available.end()); + close_group(kT_group); + + // If only one temperature is available, revert to nearest temperature + if (temps_available.size() == 1 && settings::temperature_method == TEMPERATURE_INTERPOLATION) { + if (mpi::master) { + warning("Cross sections for " + name_ + " are only available at one " + "temperature. Reverting to nearest temperature method."); + } + settings::temperature_method = TEMPERATURE_NEAREST; + } + + // 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; + double T_min = settings::temperature_range[0]; + double T_max = settings::temperature_range[1]; + if (T_max > 0.0) { + for (auto T : temps_available) { + if (T_min <= T && T <= T_max) { + temps_to_read.push_back(std::round(T)); + } + } + } + + switch (settings::temperature_method) { + case TEMPERATURE_NEAREST: + // Find nearest temperatures + for (int i = 0; i < n; ++i) { + double T_desired = temperature[i]; + + // Determine closest temperature + double min_delta_T = INFTY; + double T_actual; + for (auto T : temps_available) { + double delta_T = std::abs(T - T_desired); + if (delta_T < min_delta_T) { + T_actual = T; + min_delta_T = delta_T; + } + } + + if (std::abs(T_actual - T_desired) < settings::temperature_tolerance) { + if (!contains(temps_to_read, std::round(T_actual))) { + temps_to_read.push_back(std::round(T_actual)); + + // Write warning for resonance scattering data if 0K is not available + if (std::abs(T_actual - T_desired) > 0 && T_desired == 0 && mpi::master) { + warning(name_ + " does not contain 0K data needed for resonance " + "scattering options selected. Using data at " + std::to_string(T_actual) + + " K instead."); + } + } + } else { + fatal_error("Nuclear data library does not contain cross sections for " + + name_ + " at or near " + std::to_string(T_desired) + " K."); + } + } + 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 (int i = 0; i < n; ++i) { + double T_desired = temperature[i]; + + bool found_pair = false; + for (int j = 0; j < temps_available.size() - 1; ++j) { + if (temps_available[j] <= T_desired && T_desired < temps_available[j + 1]) { + int T_j = std::round(temps_available[j]); + int T_j1 = std::round(temps_available[j+1]); + if (!contains(temps_to_read, T_j)) { + temps_to_read.push_back(T_j); + } + if (!contains(temps_to_read, T_j1)) { + temps_to_read.push_back(T_j1); + } + found_pair = true; + } + } + + if (!found_pair) { + fatal_error("Nuclear data library does not contain cross sections for " + + name_ +" at temperatures that bound " + std::to_string(T_desired) + " K."); + } + } + break; + } + + // Sort temperatures to read + std::sort(temps_to_read.begin(), temps_to_read.end()); + + // Read reactions + hid_t rxs_group = open_group(group, "reactions"); + for (auto name : group_names(rxs_group)) { + if (starts_with(name, "reaction_")) { + hid_t rx_group = open_group(rxs_group, name.c_str()); + reactions_.push_back(std::make_unique(rx_group, temps_to_read)); + close_group(rx_group); + } + } + close_group(rxs_group); + } //============================================================================== @@ -43,9 +165,17 @@ set_particle_energy_bounds(int particle, double E_min, double E_max) data::energy_max[particle - 1] = E_max; } -extern "C" void nuclide_from_hdf5_c(hid_t group) +extern "C" Nuclide* nuclide_from_hdf5_c(hid_t group, const double* temperature, int n) { - data::nuclides.emplace_back(group); + data::nuclides.push_back(std::make_unique(group, temperature, n)); + return data::nuclides.back().get(); } +extern "C" Reaction* nuclide_reaction(Nuclide* nuc, int i_rx) +{ + return nuc->reactions_[i_rx-1].get(); +} + +extern "C" void nuclides_clear() { data::nuclides.clear(); } + } // namespace openmc diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index aaf619d8c..3fb19b60f 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -105,6 +105,8 @@ module nuclide_header class(Function1D), allocatable :: fission_q_prompt ! fragments and prompt neutrons, gammas class(Function1D), allocatable :: fission_q_recov ! fragments, neutrons, gammas, betas + type(C_PTR) :: ptr + contains procedure :: assign_0K_elastic_scattering procedure :: clear => nuclide_clear @@ -312,18 +314,6 @@ contains class(Nuclide), intent(inout) :: this ! The Nuclide object to clear integer :: i - interface - subroutine reaction_delete(rx) bind(C) - import C_PTR - type(C_PTR), value :: rx - end subroutine reaction_delete - end interface - - do i = 1, size(this % reactions) - call reaction_delete(this % reactions(i) % ptr) - end do - deallocate(this % reactions) - if (associated(this % multipole)) deallocate(this % multipole) end subroutine nuclide_clear @@ -364,14 +354,17 @@ contains type(VectorInt) :: index_inelastic_scatter interface - subroutine nuclide_from_hdf5_c(group) bind(C) - import HID_T + function nuclide_from_hdf5_c(group, temperature, n) result(ptr) bind(C) + import HID_T, C_DOUBLE, C_INT, C_PTR integer(HID_T), value :: group - end subroutine + real(C_DOUBLE), intent(in) :: temperature(*) + integer(C_INT), value :: n + type(C_PTR) :: ptr + end function end interface ! Read data on C++ side - call nuclide_from_hdf5_c(group_id) + this % ptr = nuclide_from_hdf5_c(group_id, temperature % data(1), temperature % size()) ! Get name of nuclide from group this % name = get_name(group_id) @@ -523,7 +516,8 @@ contains rx_group = open_group(rxs_group, 'reaction_' // trim(& zero_padded(MTs % data(i), 3))) - call this % reactions(i) % from_hdf5(rx_group, temps_to_read) + ! Set pointer for each reaction + call this % reactions(i) % init(this % ptr, i) ! Check for 0K elastic scattering if (this % reactions(i) % MT == 2) then @@ -1583,6 +1577,9 @@ contains interface subroutine library_clear() bind(C) end subroutine + + subroutine nuclides_clear() bind(C) + end subroutine end interface ! Deallocate cross section data, listings, and cache @@ -1592,6 +1589,7 @@ contains call nuclides(i) % clear() end do deallocate(nuclides) + call nuclides_clear() end if n_nuclides = 0 diff --git a/src/reaction.cpp b/src/reaction.cpp index 12880bad5..918d61fea 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -275,14 +275,6 @@ std::string reaction_name(int mt) // Fortran compatibility functions //============================================================================== -Reaction* reaction_from_hdf5(hid_t group, int* temperatures, int n) -{ - std::vector temps {temperatures, temperatures + n}; - return new Reaction{group, temps}; -} - -void reaction_delete(Reaction* rx) { delete rx; } - int reaction_mt(Reaction* rx) { return rx->mt_; } double reaction_q_value(Reaction* rx) { return rx->q_value_; } diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index 0762846ef..9d57db5cc 100644 --- a/src/reaction_header.F90 +++ b/src/reaction_header.F90 @@ -22,7 +22,7 @@ module reaction_header logical(C_BOOL) :: scatter_in_cm ! scattering system in center-of-mass? logical(C_BOOL) :: redundant ! redundant reaction? contains - procedure :: from_hdf5 + procedure :: init procedure :: mt_ procedure :: q_value_ procedure :: scatter_in_cm_ @@ -40,11 +40,10 @@ module reaction_header end type Reaction interface - function reaction_from_hdf5(group, temperatures, n) result(ptr) bind(C) - import C_PTR, HID_T, C_INT - integer(HID_T), value :: group - integer(C_INT), intent(in) :: temperatures - integer(C_INT), value :: n + function nuclide_reaction(nuc_ptr, i_rx) result(ptr) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: nuc_ptr + integer(C_INT), value :: i_rx type(C_PTR) :: ptr end function @@ -148,27 +147,17 @@ module reaction_header contains - subroutine from_hdf5(this, group_id, temperatures) + subroutine init(this, nuc_ptr, i_rx) class(Reaction), intent(inout) :: this - integer(HID_T), intent(in) :: group_id - type(VectorInt), intent(in) :: temperatures + type(C_PTR), intent(in) :: nuc_ptr + integer(C_INT), intent(in) :: i_rx - integer(C_INT) :: dummy - integer(C_INT) :: n - - n = temperatures % size() - if (n > 0) then - this % ptr = reaction_from_hdf5(group_id, temperatures % data(1), n) - else - ! In this case, temperatures % data(1) doesn't exist, so we just pass a - ! dummy value - this % ptr = reaction_from_hdf5(group_id, dummy, n) - end if + this % ptr = nuclide_reaction(nuc_ptr, i_rx) this % MT = reaction_mt(this % ptr) this % Q_value = reaction_q_value(this % ptr) this % scatter_in_cm = reaction_scatter_in_cm(this % ptr) this % redundant = reaction_redundant(this % ptr) - end subroutine from_hdf5 + end subroutine function mt_(this) result(mt) class(Reaction), intent(in) :: this