Add Material::init_thermal()

include/openmc/material.h

@@ -34,6 +34,13 @@ extern std::unordered_map<int32_t, int32_t> material_map;
 class Material
 {
 public:
+  // Types
+  struct ThermalTable {
+    int index_table; //!< Index of table in data::thermal_scatt
+    int index_nuclide; //!< Index in nuclide_
+    double fraction; //!< How often to use table
+  };
+
   // Constructors
   Material() {};
   explicit Material(pugi::xml_node material_node);
@@ -44,6 +51,9 @@ public:
   //! Initialize bremsstrahlung data
   void init_bremsstrahlung();
 
+  //! Initialize thermal scattering mapping
+  void init_thermal();
+
   // Data
   int32_t id_; //!< Unique ID
   std::string name_; //!< Name of material
@@ -57,10 +67,8 @@ public:
   bool has_isotropic_nuclides_ {false};
   std::vector<bool> p0_; //!< Indicate which nuclides are to be treated with iso-in-lab scattering
 
-  // S(a,b) data
-  std::vector<int> i_sab_nuclides_; //!< Indices of nuclides that have S(a,b) table
-  std::vector<int> i_sab_tables_; //!< Corresponding indices in data::thermal_scatt
-  std::vector<double> sab_fracs_; //!< How often to use table
+  // Thermal scattering tables
+  std::vector<ThermalTable> thermal_tables_;
 
   //! \brief Default temperature for cells containing this material.
   //!

include/openmc/thermal.h

@@ -4,6 +4,7 @@
 #include <cstddef>
 #include <memory>
 #include <string>
+#include <unordered_map>
 #include <vector>
 
 #include "xtensor/xtensor.hpp"
@@ -36,6 +37,7 @@ class ThermalScattering;
 
 namespace data {
 extern std::vector<std::unique_ptr<ThermalScattering>> thermal_scatt;
+extern std::unordered_map<std::string, int> thermal_scatt_map;
 }
 
 //==============================================================================

src/material.cpp

@@ -1,6 +1,6 @@
 #include "openmc/material.h"
 
-#include <algorithm> // for min, max
+#include <algorithm> // for min, max, sort
 #include <cmath>
 #include <iterator>
 #include <string>
@@ -11,6 +11,7 @@
 #include "xtensor/xview.hpp"
 
 #include "openmc/cross_sections.h"
+#include "openmc/container_util.h"
 #include "openmc/error.h"
 #include "openmc/math_functions.h"
 #include "openmc/nuclide.h"
@@ -145,7 +146,7 @@ Material::Material(pugi::xml_node node)
     }
   } else {
     // Create list of nuclides based on those specified
-    for (auto node_nuc  : node.children("nuclide")) {
+    for (auto node_nuc : node.children("nuclide")) {
       // Check for empty name on nuclide
       if (!check_for_node(node_nuc, "name")) {
         fatal_error("No name specified on nuclide in material "
@@ -253,9 +254,7 @@ Material::Material(pugi::xml_node node)
       // Apply isotropic-in-lab treatment to specified nuclides
       for (const auto& nuc : iso_lab) {
         for (int j = 0; j < names.size(); ++j) {
-          if (names[j] == nuc) {
-            p0_[j] = true;
-          }
+          if (names[j] == nuc) p0_[j] = true;
         }
       }
     }
@@ -279,6 +278,100 @@ Material::Material(pugi::xml_node node)
   if (check_for_node(node, "volume")) {
     volume_ = std::stod(get_node_value(node, "volume"));
   }
+
+  // =======================================================================
+  // READ AND PARSE <sab> TAG FOR THERMAL SCATTERING DATA
+  if (settings::run_CE) {
+    // Loop over <sab> elements
+
+    std::vector<std::string> sab_names;
+    for (auto node_sab : node.children("sab")) {
+      // Determine name of thermal scattering table
+      if (!check_for_node(node_sab, "name")) {
+        fatal_error("Need to specify <name> for thermal scattering table.");
+      }
+      std::string name = get_node_value(node_sab, "name");
+      sab_names.push_back(name);
+
+      // Read the fraction of nuclei affected by this thermal scattering table
+      double fraction = 1.0;
+      if (check_for_node(node_sab, "fraction")) {
+        fraction = std::stod(get_node_value(node_sab, "fraction"));
+      }
+
+      // Check that the thermal scattering table is listed in the
+      // cross_sections.xml file
+      LibraryKey key {Library::Type::thermal, name};
+      if (data::library_map.find(key) == data::library_map.end()) {
+        fatal_error("Could not find thermal scattering data " + name +
+          " in cross_sections.xml file.");
+      }
+
+      // Determine index of thermal scattering data in global
+      // data::thermal_scatt array
+      int index_table;
+      if (data::thermal_scatt_map.find(name) == data::thermal_scatt_map.end()) {
+        index_table = data::thermal_scatt_map.size();
+        data::thermal_scatt_map[name] = index_table;
+      } else {
+        index_table = data::thermal_scatt_map[name];
+      }
+
+      // Add entry to thermal tables vector. For now, we put the nuclide index
+      // as zero since we don't know which nuclides the table is being applied
+      // to yet (this is assigned in init_thermal)
+      thermal_tables_.push_back({index_table, 0, fraction});
+    }
+  }
+}
+
+void Material::init_thermal()
+{
+  std::vector<ThermalTable> tables;
+
+  for (const auto& table : thermal_tables_) {
+    // In order to know which nuclide the S(a,b) table applies to, we need
+    // to search through the list of nuclides for one which has a matching
+    // name
+    bool found = false;
+    for (int j = 0; j < nuclide_.size(); ++j) {
+      const auto& name {data::nuclides[nuclide_[j]]->name_};
+      if (contains(data::thermal_scatt[table.index_table]->nuclides_, name)) {
+        tables.push_back({table.index_table, j, table.fraction});
+        found = true;
+      }
+    }
+
+    // Check to make sure thermal scattering table matched a nuclide
+    if (!found) {
+      fatal_error("Thermal scattering table " + data::thermal_scatt[
+        table.index_table]->name_  + " did not match any nuclide on material "
+        + std::to_string(id_));
+    }
+  }
+
+  // Make sure each nuclide only appears in one table.
+  for (int j = 0; j < tables.size(); ++j) {
+    for (int k = j+1; k < tables.size(); ++k) {
+      if (tables[j].index_nuclide == tables[k].index_nuclide) {
+        int index = nuclide_[tables[j].index_nuclide];
+        auto name = data::nuclides[index]->name_;
+        fatal_error(name + " in material " + std::to_string(id_) + " was found "
+          "in multiple thermal scattering tables. Each nuclide can appear in "
+          "only one table per material.");
+      }
+    }
+  }
+
+  // If there are multiple S(a,b) tables, we need to make sure that the
+  // entries in i_sab_nuclides are sorted or else they won't be applied
+  // correctly in the cross_section module.
+  std::sort(tables.begin(), tables.end(), [](ThermalTable a, ThermalTable b) {
+    return a.index_nuclide < b.index_nuclide;
+  });
+
+  // Update the list of thermal tables
+  thermal_tables_ = tables;
 }
 
 void Material::init_bremsstrahlung()
@@ -495,7 +588,7 @@ void Material::calculate_neutron_xs(const Particle& p) const
   int i_grid = std::log(p.E/data::energy_min[neutron])/simulation::log_spacing;
 
   // Determine if this material has S(a,b) tables
-  bool check_sab = (i_sab_tables_.size() > 0);
+  bool check_sab = (thermal_tables_.size() > 0);
 
   // Initialize position in i_sab_nuclides
   int j = 0;
@@ -511,10 +604,11 @@ void Material::calculate_neutron_xs(const Particle& p) const
     // Check if this nuclide matches one of the S(a,b) tables specified.
     // This relies on i_sab_nuclides being in sorted order
     if (check_sab) {
-      if (i == i_sab_nuclides_[j]) {
+      const auto& sab {thermal_tables_[j]};
+      if (i == sab.index_nuclide) {
         // Get index in sab_tables
-        i_sab = i_sab_tables_[j];
-        sab_frac = sab_fracs_[j];
+        i_sab = sab.index_table;
+        sab_frac = sab.fraction;
 
         // If particle energy is greater than the highest energy for the
         // S(a,b) table, then don't use the S(a,b) table
@@ -524,7 +618,7 @@ void Material::calculate_neutron_xs(const Particle& p) const
         ++j;
 
         // Don't check for S(a,b) tables if there are no more left
-        if (j == i_sab_tables_.size()) check_sab = false;
+        if (j == thermal_tables_.size() - 1) check_sab = false;
       }
     }
 
@@ -551,6 +645,7 @@ void Material::calculate_neutron_xs(const Particle& p) const
     double atom_density = atom_density_(i);
 
     // Add contributions to cross sections
+    const auto& macro {simulation::material_xs};
     simulation::material_xs.total += atom_density * micro.total;
     simulation::material_xs.absorption += atom_density * micro.absorption;
     simulation::material_xs.fission += atom_density * micro.fission;
@@ -716,6 +811,11 @@ extern "C" {
     mat->init_bremsstrahlung();
   }
 
+  void material_calculate_xs_c(Material* mat, const Particle* p)
+  {
+    mat->calculate_xs(*p);
+  }
+
   void extend_materials_c(int32_t n)
   {
     model::materials.reserve(model::materials.size() + n);

src/material_header.F90

@@ -360,18 +360,15 @@ contains
     class(Material), intent(in) :: this
     type(Particle),  intent(in) :: p
 
-    ! Set all material macroscopic cross sections to zero
-    material_xs % total           = ZERO
-    material_xs % absorption      = ZERO
-    material_xs % fission         = ZERO
-    material_xs % nu_fission      = ZERO
-
-    if (p % type == NEUTRON) then
-      call this % calculate_neutron_xs(p)
-    elseif (p % type == PHOTON) then
-      call this % calculate_photon_xs(p)
-    end if
+    interface
+      subroutine material_calculate_xs_c(ptr, p) bind(C)
+        import C_PTR, Particle
+        type(C_PTR), value :: ptr
+        type(Particle), intent(in) :: p
+      end subroutine
+    end interface
 
+    call material_calculate_xs_c(this % ptr, p)
   end subroutine material_calculate_xs
 
 !===============================================================================

src/physics.cpp

@@ -447,7 +447,7 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat)
   double prob = 0.0;
   while (prob < cutoff) {
     // Check to make sure that a nuclide was sampled
-    if (*i_nuc_mat > n) {
+    if (*i_nuc_mat >= n) {
       p->write_restart();
       fatal_error("Did not sample any nuclide during collision.");
     }

src/thermal.cpp

@@ -26,6 +26,7 @@ namespace openmc {
 
 namespace data {
 std::vector<std::unique_ptr<ThermalScattering>> thermal_scatt;
+std::unordered_map<std::string, int> thermal_scatt_map;
 }
 
 //==============================================================================