Use openmc_load_nuclide consistently

include/openmc/capi.h

@@ -55,7 +55,7 @@ extern "C" {
   bool openmc_is_statepoint_batch();
   int openmc_legendre_filter_get_order(int32_t index, int* order);
   int openmc_legendre_filter_set_order(int32_t index, int order);
-  int openmc_load_nuclide(const char* name);
+  int openmc_load_nuclide(const char* name, const double* temps, int n);
   int openmc_material_add_nuclide(int32_t index, const char name[], double density);
   int openmc_material_get_densities(int32_t index, const int** nuclides, const double** densities, int* n);
   int openmc_material_get_id(int32_t index, int32_t* id);

openmc/lib/nuclide.py

@@ -1,5 +1,5 @@
 from collections.abc import Mapping
-from ctypes import c_int, c_char_p, POINTER, c_size_t
+from ctypes import c_int, c_double, c_char_p, POINTER, c_size_t
 from weakref import WeakValueDictionary
 
 from ..exceptions import DataError, AllocationError
@@ -14,7 +14,7 @@ __all__ = ['Nuclide', 'nuclides', 'load_nuclide']
 _dll.openmc_get_nuclide_index.argtypes = [c_char_p, POINTER(c_int)]
 _dll.openmc_get_nuclide_index.restype = c_int
 _dll.openmc_get_nuclide_index.errcheck = _error_handler
-_dll.openmc_load_nuclide.argtypes = [c_char_p]
+_dll.openmc_load_nuclide.argtypes = [c_char_p, POINTER(c_double), c_int]
 _dll.openmc_load_nuclide.restype = c_int
 _dll.openmc_load_nuclide.errcheck = _error_handler
 _dll.openmc_nuclide_name.argtypes = [c_int, POINTER(c_char_p)]
@@ -32,7 +32,7 @@ def load_nuclide(name):
         Name of the nuclide, e.g. 'U235'
 
     """
-    _dll.openmc_load_nuclide(name.encode())
+    _dll.openmc_load_nuclide(name.encode(), None, 0)
 
 
 class Nuclide(_FortranObject):

src/cross_sections.cpp

@@ -1,5 +1,6 @@
 #include "openmc/cross_sections.h"
 
+#include "openmc/capi.h"
 #include "openmc/constants.h"
 #include "openmc/container_util.h"
 #ifdef DAGMC
@@ -207,56 +208,11 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
       // If we've already read this nuclide, skip it
       if (already_read.find(name) != already_read.end()) continue;
 
-      LibraryKey key {Library::Type::neutron, name};
-      int idx = data::library_map[key];
-      std::string& filename = data::libraries[idx].path_;
+      const auto& temps = nuc_temps[i_nuc];
+      int err = openmc_load_nuclide(name.c_str(), temps.data(), temps.size());
+      if (err < 0) throw std::runtime_error{openmc_err_msg};
 
-      write_message("Reading " + name + " from " + filename, 6);
-
-      // Open file and make sure version is sufficient
-      hid_t file_id = file_open(filename, 'r');
-      check_data_version(file_id);
-
-      // Read nuclide data from HDF5
-      hid_t group = open_group(file_id, name.c_str());
-      int i_nuclide = data::nuclides.size();
-      data::nuclides.push_back(std::make_unique<Nuclide>(
-        group, nuc_temps[i_nuc]));
-
-      close_group(group);
-      file_close(file_id);
-
-      // Add name and alias to dictionary
       already_read.insert(name);
-
-      // Check if elemental data has been read, if needed
-      std::string element = to_element(name);
-      if (settings::photon_transport) {
-        if (already_read.find(element) == already_read.end()) {
-          // Read photon interaction data from HDF5 photon library
-          LibraryKey key {Library::Type::photon, element};
-          int idx = data::library_map[key];
-          std::string& filename = data::libraries[idx].path_;
-          write_message("Reading " + element + " from " + filename, 6);
-
-          // Open file and make sure version is sufficient
-          hid_t file_id = file_open(filename, 'r');
-          check_data_version(file_id);
-
-          // Read element data from HDF5
-          hid_t group = open_group(file_id, element.c_str());
-          data::elements.emplace_back(group);
-
-          close_group(group);
-          file_close(file_id);
-
-          // Add element to set
-          already_read.insert(element);
-        }
-      }
-
-      // Read multipole file into the appropriate entry on the nuclides array
-      if (settings::temperature_multipole) read_multipole_data(i_nuclide);
     }
   }
 

src/material.cpp

@@ -943,7 +943,7 @@ void Material::set_densities(const std::vector<std::string>& name,
   for (gsl::index i = 0; i < n; ++i) {
     const auto& nuc {name[i]};
     if (data::nuclide_map.find(nuc) == data::nuclide_map.end()) {
-      int err = openmc_load_nuclide(nuc.c_str());
+      int err = openmc_load_nuclide(nuc.c_str(), nullptr, 0);
       if (err < 0) throw std::runtime_error{openmc_err_msg};
     }
 
@@ -1053,7 +1053,7 @@ void Material::add_nuclide(const std::string& name, double density)
   }
 
   // If nuclide wasn't found, extend nuclide/density arrays
-  int err = openmc_load_nuclide(name.c_str());
+  int err = openmc_load_nuclide(name.c_str(), nullptr, 0);
   if (err < 0) throw std::runtime_error{openmc_err_msg};
 
   // Append new nuclide/density

src/nuclide.cpp

@@ -927,10 +927,12 @@ nuclides_size()
 // C API
 //==============================================================================
 
-extern "C" int openmc_load_nuclide(const char* name)
+extern "C" int openmc_load_nuclide(const char* name, const double* temps, int n)
 {
-  if (data::nuclide_map.find(name) == data::nuclide_map.end()) {
-    const auto& it = data::library_map.find({Library::Type::neutron, name});
+  if (data::nuclide_map.find(name) == data::nuclide_map.end() ||
+      data::nuclide_map.at(name) >= data::elements.size()) {
+    LibraryKey key {Library::Type::neutron, name};
+    const auto& it = data::library_map.find(key);
     if (it == data::library_map.end()) {
       set_errmsg("Nuclide '" + std::string{name} + "' is not present in library.");
       return OPENMC_E_DATA;
@@ -947,7 +949,7 @@ extern "C" int openmc_load_nuclide(const char* name)
 
     // Read nuclide data from HDF5
     hid_t group = open_group(file_id, name);
-    std::vector<double> temperature;
+    std::vector<double> temperature{temps, temps + n};
     data::nuclides.push_back(std::make_unique<Nuclide>(group, temperature));
 
     close_group(group);
@@ -960,10 +962,17 @@ extern "C" int openmc_load_nuclide(const char* name)
     // Read elemental data, if necessary
     if (settings::photon_transport) {
       auto element = to_element(name);
-      if (data::element_map.find(element) == data::element_map.end()) {
+      if (data::element_map.find(element) == data::element_map.end() ||
+          data::element_map.at(element) >= data::elements.size()) {
         // Read photon interaction data from HDF5 photon library
         LibraryKey key {Library::Type::photon, element};
-        int idx = data::library_map[key];
+        const auto& it = data::library_map.find(key);
+        if (it == data::library_map.end()) {
+          set_errmsg("Element '" + std::string{element} + "' is not present in library.");
+          return OPENMC_E_DATA;
+        }
+
+        int idx = it->second;
         const auto& filename = data::libraries[idx].path_;
         write_message("Reading " + element + " from " + filename, 6);