Simplify bindings for sab_tables

include/openmc/thermal.h

@@ -103,8 +103,7 @@ private:
 
 class ThermalScattering {
 public:
-  ThermalScattering(hid_t group, const std::vector<double>& temperature, int method,
-                    double tolerance, const double* minmax);
+  ThermalScattering(hid_t group, const std::vector<double>& temperature);
 
   //! Determine inelastic/elastic cross section at given energy
   //!

src/input_xml.F90

@@ -2023,8 +2023,10 @@ contains
     integer :: i_nuclide
     integer :: i_element
     integer :: i_sab
+    integer(C_INT) :: n
     integer(HID_T) :: file_id
     integer(HID_T) :: group_id
+    real(C_DOUBLE) :: dummy
     logical :: mp_found     ! if windowed multipole libraries were found
     character(MAX_WORD_LEN) :: name
     character(MAX_FILE_LEN) :: filename
@@ -2034,7 +2036,6 @@ contains
 
     allocate(nuclides(n_nuclides))
     allocate(elements(n_elements))
-    allocate(sab_tables(n_sab_tables))
     if (photon_transport .and. electron_treatment == ELECTRON_TTB) then
       allocate(ttb(n_materials))
     end if
@@ -2183,8 +2184,14 @@ contains
 
           ! Read S(a,b) data from HDF5
           group_id = open_group(file_id, name)
-          call sab_tables(i_sab) % from_hdf5(group_id, sab_temps(i_sab), &
-               temperature_method, temperature_tolerance, temperature_range)
+          n = sab_temps(i_sab) % size()
+          if (n > 0) then
+            call sab_from_hdf5(group_id, sab_temps(i_sab) % data(1), n)
+          else
+            ! In this case, data(1) doesn't exist, so we just pass a dummy value
+            call sab_from_hdf5(group_id, dummy, n)
+          end if
+
           call close_group(group_id)
           call file_close(file_id)
 

src/material_header.F90

@@ -12,7 +12,7 @@ module material_header
   use sab_header
   use simulation_header, only: log_spacing
   use stl_vector, only: VectorReal, VectorInt
-  use string, only: to_str, to_f_string
+  use string, only: to_str, to_f_string, to_c_string
 
   implicit none
 
@@ -251,6 +251,7 @@ contains
     integer :: j            ! index over nuclides in material
     integer :: k            ! index over S(a,b) tables in material
     integer :: m            ! position for sorting
+    integer :: i_sab
     integer :: temp_nuclide ! temporary value for sorting
     integer :: temp_table   ! temporary value for sorting
     real(8) :: temp_frac    ! temporary value for sorting
@@ -266,16 +267,15 @@ contains
       ! to search through the list of nuclides for one which has a matching
       ! name
       found = .false.
-      associate (sab => sab_tables(this % i_sab_tables(k)))
-        FIND_NUCLIDE: do j = 1, size(this % nuclide)
-          if (sab % has_nuclide(nuclides(this % nuclide(j)) % name)) then
-            call i_sab_tables % push_back(this % i_sab_tables(k))
-            call i_sab_nuclides % push_back(j)
-            call sab_fracs % push_back(this % sab_fracs(k))
-            found = .true.
-          end if
-        end do FIND_NUCLIDE
-      end associate
+      i_sab = this % i_sab_tables(k)
+      FIND_NUCLIDE: do j = 1, size(this % nuclide)
+        if (sab_has_nuclide(i_sab, to_c_string(nuclides(this % nuclide(j)) % name))) then
+          call i_sab_tables % push_back(i_sab)
+          call i_sab_nuclides % push_back(j)
+          call sab_fracs % push_back(this % sab_fracs(k))
+          found = .true.
+        end if
+      end do FIND_NUCLIDE
 
       ! Check to make sure S(a,b) table matched a nuclide
       if (.not. found) then
@@ -422,7 +422,7 @@ contains
 
           ! If particle energy is greater than the highest energy for the
           ! S(a,b) table, then don't use the S(a,b) table
-          if (p % E > sab_tables(i_sab) % threshold()) then
+          if (p % E > sab_threshold(i_sab)) then
             i_sab = 0
           end if
 

src/nuclide_header.F90

@@ -18,11 +18,9 @@ module nuclide_header
   use message_passing
   use random_lcg,             only: prn, future_prn, prn_set_stream
   use reaction_header,        only: Reaction
-  use sab_header,             only: SAlphaBeta, sab_tables
   use settings
   use stl_vector,             only: VectorInt, VectorReal
   use string
-  use simulation_header,      only: need_depletion_rx
 
   implicit none
 
@@ -1102,5 +1100,4 @@ contains
     E = nuclides(i_nuclide + 1) % multipole % E_max
   end function
 
-
 end module nuclide_header

src/output.F90

@@ -13,7 +13,6 @@ module output
   use mgxs_interface
   use nuclide_header
   use particle_header, only: LocalCoord, Particle
-  use sab_header,      only: SAlphaBeta
   use settings
   use simulation_header
   use surface_header,  only: surfaces

src/physics.F90

@@ -13,7 +13,6 @@ module physics
                                     thick_target_bremsstrahlung
   use physics_common
   use random_lcg,             only: prn
-  use sab_header,             only: sab_tables
   use settings
   use simulation_header
   use tally_header

src/sab_header.F90

@@ -4,100 +4,43 @@ module sab_header
 
   use dict_header, only: DictCharInt
   use hdf5_interface
-  use stl_vector, only: VectorReal
-  use string, only: to_c_string
 
   implicit none
   private
 
-  public :: free_memory_sab
-
-!===============================================================================
-! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off
-! of light isotopes such as water, graphite, Be, etc
-!===============================================================================
-
-  type, public :: SAlphaBeta
-    type(C_PTR) :: ptr
-  contains
-    procedure :: from_hdf5
-    procedure :: has_nuclide
-    procedure :: threshold
-  end type SAlphaBeta
+  public :: free_memory_sab, sab_from_hdf5, sab_has_nuclide, sab_threshold
 
   ! S(a,b) tables
-  type(SAlphaBeta), public, allocatable, target :: sab_tables(:)
-  integer(C_INT), public, bind(C) :: n_sab_tables
+  integer, public :: n_sab_tables
   type(DictCharInt), public :: sab_dict
 
   interface
-    function sab_from_hdf5(group_id, temperature, n, method, &
-         tolerance, minmax) result(ptr) bind(C)
-      import HID_T, C_DOUBLE, C_INT, C_PTR
+    subroutine sab_from_hdf5(group_id, temperature, n) bind(C)
+      import HID_T, C_DOUBLE, C_INT
       integer(HID_T), value :: group_id
       real(C_DOUBLE), intent(in) :: temperature
       integer(C_INT), value :: n
-      integer(C_INT), value :: method
-      real(C_DOUBLE), value :: tolerance
-      real(C_DOUBLE), intent(in) :: minmax(2)
-      type(C_PTR) :: ptr
-    end function
+    end subroutine
 
     subroutine sab_clear() bind(C)
     end subroutine
 
-    function sab_has_nuclide(ptr, name) result(val) bind(C)
-      import C_PTR, C_CHAR, C_BOOL
-      type(C_PTR), value :: ptr
+    function sab_has_nuclide(i_sab, name) result(val) bind(C)
+      import C_INT, C_CHAR, C_BOOL
+      integer(C_INT), value :: i_sab
       character(kind=C_CHAR), intent(in) :: name(*)
       logical(C_BOOL) :: val
     end function
 
-    function sab_threshold(ptr) result(threshold) bind(C)
-      import C_PTR, C_double
-      type(C_PTR), value :: ptr
+    function sab_threshold(i_sab) result(threshold) bind(C)
+      import C_INT, C_DOUBLE
+      integer(C_INT), value :: i_sab
       real(C_DOUBLE) :: threshold
     end function
   end interface
 
 contains
 
-  subroutine from_hdf5(this, group_id, temperature, method, &
-       tolerance, minmax)
-    class(SAlphaBeta), intent(inout) :: this
-    integer(HID_T),    intent(in)    :: group_id
-    type(VectorReal),  intent(in)    :: temperature ! list of temperatures
-    integer(C_INT),    intent(in)    :: method
-    real(C_DOUBLE),    intent(in)    :: tolerance
-    real(C_DOUBLE),    intent(in)    :: minmax(2)
-
-    real(C_DOUBLE) :: dummy
-    integer(C_INT) :: n
-
-    n = temperature % size()
-    if (n > 0) then
-      this % ptr = sab_from_hdf5(group_id, temperature % data(1), n, method, tolerance, minmax)
-    else
-      ! In this case, temperatures % data(1) doesn't exist, so we just pass a
-      ! dummy value
-      this % ptr = sab_from_hdf5(group_id, dummy, n, method, tolerance, minmax)
-    end if
-  end subroutine from_hdf5
-
-  function has_nuclide(this, name) result(val)
-    class(SAlphaBeta), intent(in) :: this
-    character(len=*), intent(in) :: name
-    logical(C_BOOL) :: val
-
-    val = sab_has_nuclide(this % ptr, to_c_string(name))
-  end function
-
-  function threshold(this)
-    class(SAlphaBeta), intent(in) :: this
-    real(C_DOUBLE) :: threshold
-    threshold = sab_threshold(this % ptr)
-  end function
-
 !===============================================================================
 ! FREE_MEMORY_SAB deallocates global arrays defined in this module
 !===============================================================================
@@ -105,7 +48,6 @@ contains
   subroutine free_memory_sab()
     n_sab_tables = 0
     call sab_clear()
-    if (allocated(sab_tables)) deallocate(sab_tables)
     call sab_dict % clear()
   end subroutine free_memory_sab
 

src/thermal.cpp

@@ -32,8 +32,7 @@ std::vector<std::unique_ptr<ThermalScattering>> thermal_scatt;
 // ThermalScattering implementation
 //==============================================================================
 
-ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& temperature,
-                                     int method, double tolerance, const double* minmax)
+ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& temperature)
 {
   // Get name of table from group
   name_ = object_name(group);
@@ -73,22 +72,23 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& tem
   // 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;
-  if (minmax[1] > 0.0) {
+  if (settings::temperature_range[1] > 0.0) {
     for (const auto& T : temps_available) {
-      if (minmax[0] <= T && T <= minmax[1]) {
+      if (settings::temperature_range[0] <= T &&
+          T <= settings::temperature_range[1]) {
         temps_to_read.push_back(std::round(T));
       }
     }
   }
 
-  switch (method) {
+  switch (settings::temperature_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::abs(temp_actual - T) < tolerance) {
+      if (std::abs(temp_actual - T) < settings::temperature_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));
@@ -591,29 +591,27 @@ ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
 //==============================================================================
 
 extern "C" ThermalScattering*
-sab_from_hdf5(hid_t group, const double* temperature, int n,
-    int method, double tolerance, const double* minmax)
+sab_from_hdf5(hid_t group, const double* temperature, int n)
 {
   // Convert temperatures to a vector
   std::vector<double> T {temperature, temperature + n};
 
   // Create new object and return it
-  data::thermal_scatt.push_back(std::make_unique<ThermalScattering>(
-    group, T, method, tolerance, minmax));
+  data::thermal_scatt.push_back(std::make_unique<ThermalScattering>(group, T));
 
   return data::thermal_scatt.back().get();
 }
 
 extern "C" void sab_clear() { data::thermal_scatt.clear(); }
 
-extern "C" bool sab_has_nuclide(ThermalScattering* data, const char* name)
+extern "C" bool sab_has_nuclide(int i_sab, const char* name)
 {
-  return data->has_nuclide(name);
+  return data::thermal_scatt[i_sab - 1]->has_nuclide(name);
 }
 
-extern "C" double sab_threshold(ThermalScattering* data)
+extern "C" double sab_threshold(int i_sab)
 {
-  return data->threshold();
+  return data::thermal_scatt[i_sab - 1]->threshold();
 }
 
 } // namespace openmc