Merge pull request #1726 from gridley/print_materials

Export of depletion results to openmc.Materials

openmc/deplete/results_list.py

@@ -6,8 +6,10 @@ import h5py
 import numpy as np
 
 from .results import Results, VERSION_RESULTS
-from openmc.checkvalue import check_filetype_version, check_value, check_type
-
+import openmc.checkvalue as cv
+from openmc.data.library import DataLibrary
+from openmc.material import Material, Materials
+from openmc.exceptions import DataError, InvalidArgumentError
 
 __all__ = ["ResultsList"]
 
@@ -34,7 +36,7 @@ class ResultsList(list):
             New instance of depletion results
         """
         with h5py.File(str(filename), "r") as fh:
-            check_filetype_version(fh, 'depletion results', VERSION_RESULTS[0])
+            cv.check_filetype_version(fh, 'depletion results', VERSION_RESULTS[0])
             new = cls()
 
             # Get number of results stored
@@ -78,8 +80,8 @@ class ResultsList(list):
             Concentration of specified nuclide in units of ``nuc_units``
 
         """
-        check_value("time_units", time_units, {"s", "d", "min", "h"})
-        check_value("nuc_units", nuc_units,
+        cv.check_value("time_units", time_units, {"s", "d", "min", "h"})
+        cv.check_value("nuc_units", nuc_units,
                     {"atoms", "atom/b-cm", "atom/cm3"})
 
         times = np.empty_like(self, dtype=float)
@@ -215,7 +217,7 @@ class ResultsList(list):
             1-D vector of time points
 
         """
-        check_type("time_units", time_units, str)
+        cv.check_type("time_units", time_units, str)
 
         times = np.fromiter(
             (r.time[0] for r in self),
@@ -269,9 +271,9 @@ class ResultsList(list):
         int
 
         """
-        check_type("time", time, numbers.Real)
-        check_type("atol", atol, numbers.Real)
-        check_type("rtol", rtol, numbers.Real)
+        cv.check_type("time", time, numbers.Real)
+        cv.check_type("atol", atol, numbers.Real)
+        cv.check_type("rtol", rtol, numbers.Real)
 
         times = self.get_times(time_units)
 
@@ -296,3 +298,74 @@ class ResultsList(list):
             "relative tolerances {} and {}.".format(
                 time, time_units, atol, rtol)
         )
+
+    def export_to_materials(self, burnup_index, nuc_with_data=None) -> Materials:
+        """Return openmc.Materials object based on results at a given step
+
+        Parameters
+        ----------
+        burn_index : int
+            Index of burnup step to evaluate. See also: get_step_where for
+            obtaining burnup step indices from other data such as the time.
+        nuc_with_data : Iterable of str, optional
+            Nuclides to include in resulting materials.
+            This can be specified if not all nuclides appearing in
+            depletion results have associated neutron cross sections, and
+            as such cannot be used in subsequent transport calculations.
+            If not provided, nuclides from the cross_sections element of
+            materials.xml will be used. If that element is not present,
+            nuclides from OPENMC_CROSS_SECTIONS will be used.
+
+        Returns
+        -------
+        mat_file : Materials
+            A modified Materials instance containing depleted material data
+            and original isotopic compositions of non-depletable materials
+        """
+        result = self[burnup_index]
+
+        # Only materials found in the original materials.xml file will be
+        # updated. If for some reason you have modified OpenMC to produce
+        # new materials as depletion takes place, this method will not
+        # work as expected and leave out that material.
+        mat_file = Materials.from_xml("materials.xml")
+
+        # Only nuclides with valid transport data will be written to
+        # the new materials XML file. The precedence of nuclides to select
+        # is first ones provided as a kwarg here, then ones specified
+        # in the materials.xml file if provided, then finally from
+        # the environment variable OPENMC_CROSS_SECTIONS.
+        if nuc_with_data:
+            cv.check_iterable_type('nuclide names', nuc_with_data, str)
+            available_cross_sections = nuc_with_data
+        else:
+            # select cross_sections.xml file to use
+            if mat_file.cross_sections:
+                this_library = DataLibrary.from_xml(path=mat_file.cross_sections)
+            else:
+                this_library = DataLibrary.from_xml()
+
+            # Find neutron libraries we have access to
+            available_cross_sections = set()
+            for lib in this_library.libraries:
+                if lib['type'] == 'neutron':
+                    available_cross_sections.update(lib['materials'])
+            if not available_cross_sections:
+                raise DataError('No neutron libraries found in cross_sections.xml')
+
+        # Overwrite material definitions, if they can be found in the depletion
+        # results, and save them to the new depleted xml file.
+        for mat in mat_file:
+            mat_id = str(mat.id)
+            if mat_id in result.mat_to_ind:
+                mat.volume = result.volume[mat_id]
+                for nuc in result.nuc_to_ind:
+                    if nuc not in available_cross_sections:
+                        continue
+                    atoms = result[0, mat_id, nuc]
+                    if atoms > 0.0:
+                        atoms_per_barn_cm = 1e-24 * atoms / mat.volume
+                        mat.remove_nuclide(nuc) # Replace if it's there
+                        mat.add_nuclide(nuc, atoms_per_barn_cm)
+
+        return mat_file

tests/regression_tests/deplete/test.py

@@ -4,6 +4,7 @@ from math import floor
 import shutil
 from pathlib import Path
 
+from difflib import unified_diff
 import numpy as np
 import pytest
 import openmc
@@ -130,3 +131,47 @@ def test_full(run_in_tmpdir, problem, multiproc):
 
     # Check that no additional tallies are loaded from the files
     assert np.all(n_tallies == 0)
+
+
+def test_depletion_results_to_material(run_in_tmpdir, problem):
+    """Checks openmc.Materials objects can be created from depletion results"""
+    # Load the reference/test results
+    path_reference = Path(__file__).with_name('test_reference.h5')
+    res_ref = openmc.deplete.ResultsList.from_hdf5(path_reference)
+
+    # Firstly need to export materials.xml file for the initial simulation state
+    geometry, lower_left, upper_right = problem
+    materials = openmc.Materials()
+    for mat in geometry.root_universe.get_all_materials().values():
+        materials.append(mat)
+    materials.export_to_xml()
+
+    # Export last step of depletion to its own openmc.Materials object,
+    # using only nuclides available in the current nuclear data library
+    last_step_materials = res_ref.export_to_materials(-1)
+
+    # Export final depletion  step materials to XML
+    output_xml_file = 'last_step_materials.xml'
+    last_step_materials.export_to_xml(path=output_xml_file)
+    with open(output_xml_file, 'r') as result_file:
+        result_file_lines = result_file.readlines()
+
+    # If updating results, do so and return. We write out the last-step
+    # depleted materials as an XML, and save the list of lines to diff.
+    reference_file = Path(__file__).with_name('last_step_reference_materials.xml')
+    if config['update']:
+        with open(reference_file, 'w') as ref_file:
+            ref_file.writelines(result_file_lines)
+        return
+
+    # Check text of final depletion point material XML matches reference
+    with open(reference_file) as ref_file:
+        reference_lines = ref_file.readlines()
+    diff_vs_expected = unified_diff(reference_lines, result_file_lines)
+
+    # Check all lines match, printing errors along the way
+    success = True 
+    for line in diff_vs_expected:
+        success = False
+        print(line.rstrip())
+    assert success