Load cell temperatures from HDF5 file indicated by new XML element

include/openmc/cell.h

@@ -48,6 +48,7 @@ class UniversePartitioner;
 namespace model {
 extern std::unordered_map<int32_t, int32_t> cell_map;
 extern vector<unique_ptr<Cell>> cells;
+extern std::string temperatures_from_h5;
 
 } // namespace model
 

openmc/geometry.py

@@ -34,19 +34,23 @@ class Geometry:
     surface_precision : int
         Number of decimal places to round to for comparing the coefficients of
         surfaces for considering them topologically equivalent.
+    temperatures_from_h5: str
+        File name to import cell temperatures
 
     """
 
     def __init__(
         self,
         root: openmc.UniverseBase | Iterable[openmc.Cell] | None = None,
         merge_surfaces: bool = False,
-        surface_precision: int = 10
+        surface_precision: int = 10,
+        temperature_from_h5: str | None = None
     ):
         self._root_universe = None
         self._offsets = {}
         self.merge_surfaces = merge_surfaces
         self.surface_precision = surface_precision
+        self.temperature_from_h5 = temperature_from_h5
         if root is not None:
             if isinstance(root, openmc.UniverseBase):
                 self.root_universe = root
@@ -89,6 +93,15 @@ def surface_precision(self, surface_precision):
         check_greater_than('surface_precision', surface_precision, 0)
         self._surface_precision = surface_precision
 
+    @property
+    def temperature_from_h5(self) -> str:
+        return self._temperature_from_h5
+
+    @temperature_from_h5.setter
+    def temperature_from_h5(self, temperature_from_h5):
+        check_type('temperature from h5', temperature_from_h5, str)
+        self._temperature_from_h5 = temperature_from_h5
+
     def add_volume_information(self, volume_calc):
         """Add volume information from a stochastic volume calculation.
 
@@ -132,6 +145,11 @@ def to_xml_element(self, remove_surfs=False) -> ET.Element:
 
         # Create XML representation
         element = ET.Element("geometry")
+        if self.temperature_from_h5 is not None:
+            element.set("temperature_from_h5", self.temperature_from_h5)
+        else:
+            element.set("temperature_from_h5", "")
+
         self.root_universe.create_xml_subelement(element)
 
         # Sort the elements in the file
@@ -199,6 +217,8 @@ def get_universe(univ_id):
                 universes[univ_id] = univ
             return universes[univ_id]
 
+        temperature_from_h5 = get_text(e, "temperature_from_h5")
+
         # Get surfaces
         surfaces = {}
         periodic = {}
@@ -256,7 +276,7 @@ def get_universe(univ_id):
         # child of any other object
         for u in universes.values():
             if not child_of[u]:
-                return cls(u)
+                return cls(u, temperature_from_h5=temperature_from_h5)
         else:
             raise ValueError('Error determining root universe.')
 

openmc/lib/core.py

@@ -304,7 +304,7 @@ def import_properties(filename):
     openmc.lib.export_properties
 
     """
-    _dll.openmc_properties_import(filename.encode())
+    _dll.openmc_properties_import(filename.encode(), True, True)
 
 
 def init(args=None, intracomm=None, output=True):

src/cell.cpp

@@ -33,6 +33,7 @@ namespace openmc {
 namespace model {
 std::unordered_map<int32_t, int32_t> cell_map;
 vector<unique_ptr<Cell>> cells;
+std::string temperatures_from_h5 = "";
 
 } // namespace model
 
@@ -407,30 +408,33 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
     }
   }
 
-  // Read the temperature element which may be distributed like materials.
-  if (check_for_node(cell_node, "temperature")) {
-    sqrtkT_ = get_node_array<double>(cell_node, "temperature");
-    sqrtkT_.shrink_to_fit();
+  if (model::temperature_from_h5 == "") {
 
-    // Make sure this is a material-filled cell.
-    if (material_.size() == 0) {
-      fatal_error(fmt::format(
-        "Cell {} was specified with a temperature but no material. Temperature"
-        "specification is only valid for cells filled with a material.",
-        id_));
-    }
+    // Read the temperature element which may be distributed like materials.
+    if (check_for_node(cell_node, "temperature")) {
+      sqrtkT_ = get_node_array<double>(cell_node, "temperature");
+      sqrtkT_.shrink_to_fit();
 
-    // Make sure all temperatures are non-negative.
-    for (auto T : sqrtkT_) {
-      if (T < 0) {
+      // Make sure this is a material-filled cell.
+      if (material_.size() == 0) {
         fatal_error(fmt::format(
-          "Cell {} was specified with a negative temperature", id_));
+          "Cell {} was specified with a temperature but no material. Temperature"
+          "specification is only valid for cells filled with a material.",
+          id_));
+      }
+
+      // Make sure all temperatures are non-negative.
+      for (auto T : sqrtkT_) {
+        if (T < 0) {
+          fatal_error(fmt::format(
+            "Cell {} was specified with a negative temperature", id_));
+        }
       }
-    }
 
-    // Convert to sqrt(k*T).
-    for (auto& T : sqrtkT_) {
-      T = std::sqrt(K_BOLTZMANN * T);
+      // Convert to sqrt(k*T).
+      for (auto& T : sqrtkT_) {
+        T = std::sqrt(K_BOLTZMANN * T);
+      }
     }
   }
 

src/geometry_aux.cpp

@@ -59,6 +59,9 @@ void read_geometry_xml(pugi::xml_node root)
   std::unordered_map<int, double> albedo_map;
   std::unordered_map<int, int> periodic_sense_map;
 
+  // confirm XML structure and whether we expect to find this here
+  model::temperature_from_h5 = get_node_value(root, "temperature_from_h5");
+
   read_surfaces(root, periodic_pairs, albedo_map, periodic_sense_map);
   read_cells(root);
   prepare_boundary_conditions(periodic_pairs, albedo_map, periodic_sense_map);
@@ -84,6 +87,10 @@ void read_geometry_xml(pugi::xml_node root)
 
   // if the root universe is DAGMC geometry, make sure the model is well-formed
   check_dagmc_root_univ();
+
+  if (model::temperature_from_h5 != "") {
+    openmc_properties_import(model::temperature_from_h5, true, false);
+  }
 }
 
 //==============================================================================

src/summary.cpp

@@ -190,7 +190,7 @@ extern "C" int openmc_properties_export(const char* filename)
   return 0;
 }
 
-extern "C" int openmc_properties_import(const char* filename)
+extern "C" int openmc_properties_import(const char* filename, const bool import_cells, const bool import_materials)
 {
   // Display output message
   auto msg = fmt::format("Importing properties from {}...", filename);
@@ -212,48 +212,51 @@ extern "C" int openmc_properties_import(const char* filename)
     return OPENMC_E_INVALID_ARGUMENT;
   }
 
-  // Make sure number of cells matches
-  auto geom_group = open_group(file, "geometry");
-  int32_t n;
-  read_attribute(geom_group, "n_cells", n);
-  if (n != openmc::model::cells.size()) {
+  if (import_cells) {
+    // Make sure number of cells matches
+    auto geom_group = open_group(file, "geometry");
+    int32_t n;
+    read_attribute(geom_group, "n_cells", n);
+    if (n != openmc::model::cells.size()) {
+      close_group(geom_group);
+      file_close(file);
+      set_errmsg(fmt::format(
+        "Number of cells in {} doesn't match current model.", filename));
+      return OPENMC_E_GEOMETRY;
+    }
+
+    // Read cell properties
+    auto cells_group = open_group(geom_group, "cells");
+    try {
+      for (const auto& c : model::cells) {
+        c->import_properties_hdf5(cells_group);
+      }
+    } catch (const std::exception& e) {
+      set_errmsg(e.what());
+      return OPENMC_E_UNASSIGNED;
+    }
+    close_group(cells_group);
     close_group(geom_group);
-    file_close(file);
-    set_errmsg(fmt::format(
-      "Number of cells in {} doesn't match current model.", filename));
-    return OPENMC_E_GEOMETRY;
   }
 
-  // Read cell properties
-  auto cells_group = open_group(geom_group, "cells");
-  try {
-    for (const auto& c : model::cells) {
-      c->import_properties_hdf5(cells_group);
+  if (import_materials) {
+    // Make sure number of cells matches
+    auto materials_group = open_group(file, "materials");
+    read_attribute(materials_group, "n_materials", n);
+    if (n != openmc::model::materials.size()) {
+      close_group(materials_group);
+      file_close(file);
+      set_errmsg(fmt::format(
+        "Number of materials in {} doesn't match current model.", filename));
+      return OPENMC_E_GEOMETRY;
     }
-  } catch (const std::exception& e) {
-    set_errmsg(e.what());
-    return OPENMC_E_UNASSIGNED;
-  }
-  close_group(cells_group);
-  close_group(geom_group);
 
-  // Make sure number of cells matches
-  auto materials_group = open_group(file, "materials");
-  read_attribute(materials_group, "n_materials", n);
-  if (n != openmc::model::materials.size()) {
+    // Read material properties
+    for (const auto& mat : model::materials) {
+      mat->import_properties_hdf5(materials_group);
+    }
     close_group(materials_group);
-    file_close(file);
-    set_errmsg(fmt::format(
-      "Number of materials in {} doesn't match current model.", filename));
-    return OPENMC_E_GEOMETRY;
   }
-
-  // Read material properties
-  for (const auto& mat : model::materials) {
-    mat->import_properties_hdf5(materials_group);
-  }
-  close_group(materials_group);
-
   // Terminate access to the file.
   file_close(file);
   return 0;

tests/unit_tests/test_lib.py

@@ -169,7 +169,7 @@ def test_properties_temperature(lib_init):
     assert cell.get_temperature() == pytest.approx(300.0)
 
     # Import properties and check that temperature is restored
-    openmc.lib.import_properties('properties.h5')
+    openmc.lib.import_properties('properties.h5', True, True)
     assert cell.get_temperature() == pytest.approx(200.0)
 
 
@@ -197,7 +197,7 @@ def test_properties_cell_density(lib_init):
     assert cell.get_density() == pytest.approx(3.0)
 
     # Import properties and check that density is restored
-    openmc.lib.import_properties('properties.h5')
+    openmc.lib.import_properties('properties.h5', True, True)
     assert cell.get_density() == pytest.approx(orig_density)
 
 
@@ -213,7 +213,7 @@ def test_properties_fail_cell(lib_init):
     # The number of cells was changed in the previous test, so the properties
     # file is no longer valid
     with pytest.raises(exc.GeometryError, match="Number of cells"):
-        openmc.lib.import_properties("properties.h5")
+        openmc.lib.import_properties("properties.h5", True, True)
 
 
 def test_material_mapping(lib_init):
@@ -267,7 +267,7 @@ def test_properties_density(lib_init):
     assert m.get_density() == pytest.approx(orig_density*2)
 
     # Import properties and check that density was restored
-    openmc.lib.import_properties('properties.h5')
+    openmc.lib.import_properties('properties.h5', True, True)
     assert m.get_density() == pytest.approx(orig_density)
 
     with pytest.raises(ValueError):
@@ -293,7 +293,7 @@ def test_properties_fail_material(lib_init):
     # The number of materials was changed in the previous test, so the properties
     # file is no longer valid
     with pytest.raises(exc.GeometryError, match="Number of materials"):
-        openmc.lib.import_properties("properties.h5")
+        openmc.lib.import_properties("properties.h5", True, True)
 
 
 def test_nuclide_mapping(lib_init):

tests/unit_tests/test_model.py

@@ -260,7 +260,7 @@ def test_import_properties(run_in_tmpdir, mpi_intracomm):
     openmc.lib.export_properties(output=False)
 
     # Import properties to existing model
-    model.import_properties("properties.h5")
+    model.import_properties("properties.h5", True, True)
 
     # Check to see that values are assigned to the C and python representations
     # First python