Store reactions as a vector on Nuclide

include/openmc/nuclide.h

@@ -5,11 +5,13 @@
 #define OPENMC_NUCLIDE_H
 
 #include <array>
+#include <memory> // for unique_ptr
 #include <vector>
 
 #include <hdf5.h>
 
 #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<std::unique_ptr<Reaction>> reactions_; //! Reactions
 };
 
 //==============================================================================
@@ -42,7 +44,7 @@ namespace data {
 extern std::array<double, 2> energy_min;
 extern std::array<double, 2> energy_max;
 
-extern std::vector<Nuclide> nuclides;
+extern std::vector<std::unique_ptr<Nuclide>> nuclides;
 
 } // namespace data
 

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);

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 <algorithm> // for sort
+#include <string> // for to_string, stoi
 
 namespace openmc {
 
@@ -13,7 +21,7 @@ namespace data {
 std::array<double, 2> energy_min {0.0, 0.0};
 std::array<double, 2> energy_max {INFTY, INFTY};
 
-std::vector<Nuclide> nuclides;
+std::vector<std::unique_ptr<Nuclide>> nuclides;
 
 } // namespace data
 
@@ -21,7 +29,7 @@ std::vector<Nuclide> 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<double> 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<int> 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<Reaction>(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<Nuclide>(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

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
 

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<int> 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_; }

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