Add properties to settings w/ documentation, c++ loading of filename, and python round-trip test (#3808)

Co-authored-by: Patrick C Shriwise <pshriwise@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Author: 3 people

docs/source/io_formats/settings.rst

@@ -555,6 +555,18 @@ generator during generation of colors in plots.
 
   *Default*: 1
 
+.. _properties_file:
+
+-----------------------------
+``<properties_file>`` Element
+-----------------------------
+
+  The ``properties_file`` element has no attributes and contains the path to a
+  properties HDF5 file to load cell temperatures/densities and material
+  densities.
+
+  *Default*: None
+
 ---------------------
 ``<ptables>`` Element
 ---------------------

include/openmc/settings.h

@@ -115,6 +115,8 @@ extern std::string path_sourcepoint;      //!< path to a source file
 extern std::string path_statepoint;       //!< path to a statepoint file
 extern std::string weight_windows_file;   //!< Location of weight window file to
                                           //!< load on simulation initialization
+extern std::string properties_file;       //!< Location of properties file to
+                                          //!< load on simulation initialization
 
 // This is required because the c_str() may not be the first thing in
 // std::string. Sometimes it is, but it seems libc++ may not be like that

openmc/lib/core.py

@@ -693,8 +693,8 @@ def __init__(self, model=None, cwd=None, **init_kwargs):
         self.model = model
 
         # Determine MPI intercommunicator
-        self.init_kwargs.setdefault('intracomm', comm)
-        self.comm = self.init_kwargs['intracomm']
+        self.comm = self.init_kwargs.get('intracomm') or comm
+        self.init_kwargs['intracomm'] = self.comm
 
     def __enter__(self):
         """Initialize the OpenMC library in a temporary directory."""

openmc/settings.py

@@ -183,6 +183,9 @@ class Settings:
        Initial seed for randomly generated plot colors.
     ptables : bool
         Determine whether probability tables are used.
+    properties_file : PathLike
+        Location of the properties file to load cell temperatures/densities
+        and material densities
     random_ray : dict
         Options for configuring the random ray solver. Acceptable keys are:
 
@@ -409,6 +412,7 @@ def __init__(self, **kwargs):
         self._atomic_relaxation = None
         self._plot_seed = None
         self._ptables = None
+        self._properties_file = None
         self._uniform_source_sampling = None
         self._seed = None
         self._stride = None
@@ -1067,6 +1071,18 @@ def temperature(self, temperature: dict):
 
         self._temperature = temperature
 
+    @property
+    def properties_file(self) -> PathLike | None:
+        return self._properties_file
+
+    @properties_file.setter
+    def properties_file(self, value: PathLike | None):
+        if value is None:
+            self._properties_file = None
+        else:
+            cv.check_type('properties file', value, PathLike)
+            self._properties_file = input_path(value)
+
     @property
     def trace(self) -> Iterable:
         return self._trace
@@ -1772,6 +1788,12 @@ def _create_temperature_subelements(self, root):
                 else:
                     element.text = str(value)
 
+    def _create_properties_file_element(self, root):
+        if self.properties_file is not None:
+            element = ET.Element("properties_file")
+            element.text = str(self.properties_file)
+            root.append(element)
+
     def _create_trace_subelement(self, root):
         if self._trace is not None:
             element = ET.SubElement(root, "trace")
@@ -2284,6 +2306,11 @@ def _temperature_from_xml_element(self, root):
         if text is not None:
             self.temperature['multipole'] = text in ('true', '1')
 
+    def _properties_file_from_xml_element(self, root):
+        text = get_text(root, 'properties_file')
+        if text is not None:
+            self.properties_file = text
+
     def _trace_from_xml_element(self, root):
         text = get_elem_list(root, "trace", int)
         if text is not None:
@@ -2522,6 +2549,7 @@ def to_xml_element(self, mesh_memo=None):
         self._create_ifp_n_generation_subelement(element)
         self._create_tabular_legendre_subelements(element)
         self._create_temperature_subelements(element)
+        self._create_properties_file_element(element)
         self._create_trace_subelement(element)
         self._create_track_subelement(element)
         self._create_ufs_mesh_subelement(element, mesh_memo)
@@ -2639,6 +2667,7 @@ def from_xml_element(cls, elem, meshes=None):
         settings._ifp_n_generation_from_xml_element(elem)
         settings._tabular_legendre_from_xml_element(elem)
         settings._temperature_from_xml_element(elem)
+        settings._properties_file_from_xml_element(elem)
         settings._trace_from_xml_element(elem)
         settings._track_from_xml_element(elem)
         settings._ufs_mesh_from_xml_element(elem, meshes)

src/finalize.cpp

@@ -140,6 +140,7 @@ int openmc_finalize()
   settings::temperature_multipole = false;
   settings::temperature_range = {0.0, 0.0};
   settings::temperature_tolerance = 10.0;
+  settings::properties_file.clear();
   settings::trigger_on = false;
   settings::trigger_predict = false;
   settings::trigger_batch_interval = 1;

src/initialize.cpp

@@ -118,6 +118,10 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
   if (!read_model_xml())
     read_separate_xml_files();
 
+  if (!settings::properties_file.empty()) {
+    openmc_properties_import(settings::properties_file.c_str());
+  }
+
   // Reset locale to previous state
   if (std::setlocale(LC_ALL, prev_locale.c_str()) == NULL) {
     fatal_error("Cannot reset locale.");

src/settings.cpp

@@ -97,6 +97,7 @@ std::string path_sourcepoint;
 std::string path_statepoint;
 const char* path_statepoint_c {path_statepoint.c_str()};
 std::string weight_windows_file;
+std::string properties_file;
 
 int32_t n_inactive {0};
 int32_t max_lost_particles {10};
@@ -751,6 +752,14 @@ void read_settings_xml(pugi::xml_node root)
     }
   }
 
+  // read properties from file
+  if (check_for_node(root, "properties_file")) {
+    properties_file = get_node_value(root, "properties_file");
+    if (!file_exists(properties_file)) {
+      fatal_error(fmt::format("File '{}' does not exist.", properties_file));
+    }
+  }
+
   // Particle trace
   if (check_for_node(root, "trace")) {
     auto temp = get_node_array<int64_t>(root, "trace");

tests/unit_tests/test_settings.py

@@ -1,8 +1,17 @@
+from pathlib import Path
+
+import h5py
+import pytest
+
 import openmc
+import openmc.lib
 import openmc.stats
 
 
 def test_export_to_xml(run_in_tmpdir):
+
+    tmp_properties_file = 'properties_test.h5'
+
     s = openmc.Settings(run_mode='fixed source', batches=1000, seed=17)
     s.generations_per_batch = 10
     s.inactive = 100
@@ -22,7 +31,7 @@ def test_export_to_xml(run_in_tmpdir):
     s.surf_source_read = {'path': 'surface_source_1.h5'}
     s.surf_source_write = {'surface_ids': [2], 'max_particles': 200}
     s.surface_grazing_ratio = 0.7
-    s.surface_grazing_cutoff = 0.1     
+    s.surface_grazing_cutoff = 0.1
     s.confidence_intervals = True
     s.ptables = True
     s.plot_seed = 100
@@ -45,6 +54,7 @@ def test_export_to_xml(run_in_tmpdir):
     s.tabular_legendre = {'enable': True, 'num_points': 50}
     s.temperature = {'default': 293.6, 'method': 'interpolation',
                      'multipole': True, 'range': (200., 1000.)}
+    s.properties_file = tmp_properties_file
     s.trace = (10, 1, 20)
     s.track = [(1, 1, 1), (2, 1, 1)]
     s.ufs_mesh = mesh
@@ -88,6 +98,7 @@ def test_export_to_xml(run_in_tmpdir):
     # Make sure exporting XML works
     s.export_to_xml()
 
+
     # Generate settings from XML
     s = openmc.Settings.from_xml()
     assert s.run_mode == 'fixed source'
@@ -111,7 +122,7 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.surf_source_read['path'].name == 'surface_source_1.h5'
     assert s.surf_source_write == {'surface_ids': [2], 'max_particles': 200}
     assert s.surface_grazing_ratio == 0.7
-    assert s.surface_grazing_cutoff == 0.1     
+    assert s.surface_grazing_cutoff == 0.1
     assert s.confidence_intervals
     assert s.ptables
     assert s.plot_seed == 100
@@ -134,6 +145,7 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.tabular_legendre == {'enable': True, 'num_points': 50}
     assert s.temperature == {'default': 293.6, 'method': 'interpolation',
                              'multipole': True, 'range': [200., 1000.]}
+    assert s.properties_file == Path(tmp_properties_file)
     assert s.trace == [10, 1, 20]
     assert s.track == [(1, 1, 1), (2, 1, 1)]
     assert isinstance(s.ufs_mesh, openmc.RegularMesh)
@@ -178,3 +190,59 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.max_secondaries == 1_000_000
     assert s.source_rejection_fraction == 0.01
     assert s.free_gas_threshold == 800.0
+
+
+def test_properties_file_load(tmp_path, mpi_intracomm):
+    model = openmc.examples.pwr_assembly()
+
+    # Session 1: export a structurally valid properties file via the C++ API,
+    # then collect the cell/material structure so we can patch it with h5py.
+    cell_instances = {}   # {cell_id: n_instances} — material cells only
+    mat_densities = {}    # {mat_id: original atom/b-cm density}
+
+    props_path = tmp_path / 'properties.h5'
+    with openmc.lib.TemporarySession(model, intracomm=mpi_intracomm):
+        openmc.lib.export_properties(str(props_path))
+        for cell_id, cell in openmc.lib.cells.items():
+            try:
+                cell.fill  # raises NotImplementedError for non-material cells
+                cell_instances[cell_id] = cell.num_instances
+            except NotImplementedError:
+                pass
+        for mat_id, mat in openmc.lib.materials.items():
+            mat_densities[mat_id] = mat.get_density('atom/b-cm')
+
+    assert any(n > 1 for n in cell_instances.values())
+
+    # Patch the exported properties file overwriting temperatures
+    # with per-instance values and scale material atom densities.
+    density_factor = 0.75
+    with h5py.File(props_path, 'r+') as f:
+        cells_grp = f['geometry/cells']
+        for cell_id, n in cell_instances.items():
+            cell_grp = cells_grp[f'cell {cell_id}']
+            del cell_grp['temperature']
+            cell_grp.create_dataset(
+                'temperature', data=[500.0 + 5.0 * i for i in range(n)]
+            )
+
+        for mat_id, orig_density in mat_densities.items():
+            f['materials'][f'material {mat_id}'].attrs['atom_density'] = \
+                orig_density * density_factor
+
+    # now apply the newly patched properties file using the settings
+    # and load the model again, checking that the new temperature and
+    # density values match those in the new file
+    model.settings.properties_file = props_path
+
+    with openmc.lib.TemporarySession(model, intracomm=mpi_intracomm):
+        for cell_id, n in cell_instances.items():
+            cell = openmc.lib.cells[cell_id]
+            for i in range(n):
+                assert cell.get_temperature(i) == pytest.approx(500.0 + 5.0 * i)
+
+        for mat_id, orig_density in mat_densities.items():
+            mat = openmc.lib.materials[mat_id]
+            assert mat.get_density('atom/b-cm') == pytest.approx(
+                orig_density * density_factor, rel=1e-5
+            )