Implement OpenMCDepcode.read_depleted_materials()

- Add helper functions
- Add tests for read_depleted_materials() and helper functions
- add `saltproc_runtime_ref` directory in `openmc_data`
- remove unneeded import statements
- make importable time variables in app.py
- add assertion for `depcode_input[ouput_path]` in test_app.py

saltproc/app.py

@@ -27,6 +27,11 @@
 DAY_UNITS = ('d', 'day')
 YEAR_UNITS = ('a', 'year', 'yr')
 
+_SECONDS_PER_DAY = 60 * 60 * 24
+_MINUTES_PER_DAY = 60 * 24
+_HOURS_PER_DAY = 24
+_DAYS_PER_YEAR = 365.25
+
 def run():
     """ Inititializes main run"""
     threads, saltproc_input = parse_arguments()
@@ -303,8 +308,8 @@ def _process_main_input_reactor_params(reactor_input,
                                                      depletion_timesteps,
                                                      codename)
 
-    reactor_input['depletion_timesteps'] = list(depletion_timesteps)
-    reactor_input['power_levels'] = list(power_levels)
+    reactor_input['depletion_timesteps'] = depletion_timesteps.tolist()
+    reactor_input['power_levels'] = power_levels.tolist()
 
     return reactor_input
 
@@ -344,13 +349,13 @@ def _scale_depletion_timesteps(timestep_units, depletion_timesteps, codename):
     # serpent base timestep units are days or mwd/kg
     if not(timestep_units in DAY_UNITS) and timestep_units.lower() != 'mwd/kg' and codename == 'serpent':
         if timestep_units in SECOND_UNITS:
-            depletion_timesteps /= 60 * 60 * 24
+            depletion_timesteps /= _SECONDS_PER_DAY
         elif timestep_units in MINUTE_UNITS:
-            depletion_timesteps /= 60 * 24
+            depletion_timesteps /= _MINUTES_PER_DAY
         elif timestep_units in HOUR_UNITS:
-            depletion_timesteps /= 24
+            depletion_timesteps /= _HOURS_PER_DAY
         elif timestep_units in YEAR_UNITS:
-            depletion_timesteps *= 365.25
+            depletion_timesteps *= _DAYS_PER_YEAR
         else:
             raise IOError(f'Unrecognized time unit: {timestep_units}')
 

saltproc/openmc_depcode.py

@@ -4,13 +4,17 @@
 import re
 import json
 from pathlib import Path
+import numpy as np
 
 from pyne import nucname as pyname
 from pyne import serpent
 import openmc
 
 from saltproc import Materialflow
 from saltproc.abc import Depcode
+from openmc.deplete.abc import _SECONDS_PER_DAY
+from openmc.deplete import Results
+from openmc.data import atomic_mass
 
 class OpenMCDepcode(Depcode):
     """Interface for running depletion steps in OpenMC, as well as obtaining
@@ -90,6 +94,27 @@ def __init__(self,
              'settings': str((output_path / 'settings.xml').resolve())}
         self.runtime_matfile = str((output_path / 'materials.xml').resolve())
 
+        self._check_for_material_names(self.runtime_matfile)
+
+    def _check_for_material_names(self, filename):
+        """Checks that all materials in the material file
+        have a name.
+
+        Parameters
+        ----------
+        filename : str
+           Path to a materials.xml file
+
+        """
+
+        materials = openmc.Materials.from_xml(filename)
+
+        material_id_to_name = {}
+        for material in materials:
+            if material.name == '':
+                raise ValueError(f"Material {material.id} has no name.")
+
+
     def read_step_metadata(self):
         """Reads OpenMC's depletion step metadata and stores it in the
         :class:`OpenMCDepcode` object's :attr:`step_metadata` attribute.
@@ -124,6 +149,103 @@ def read_depleted_materials(self, read_at_end=False):
                 :class:`Materialflow` object holding composition and properties.
 
         """
+        # Determine moment in depletion step to read data from
+        if read_at_end:
+            moment = 1
+        else:
+            moment = 0
+
+        results_file = Path(self.output_path / 'depletion_results.h5')
+        depleted_materials = {}
+        results = Results(results_file)
+        #self.days = results['DAYS'][moment]
+        depleted_openmc_materials = results.export_to_materials(moment)
+        if read_at_end:
+            starting_openmc_materials = results.export_to_materials(0)
+        else:
+            # placeholder for starting materials
+            starting_openmc_materials = np.zeros(len(depleted_openmc_materials))
+
+        openmc_materials = zip(starting_openmc_materials, depleted_openmc_materials)
+
+        depleted_materials = {}
+        for starting_material, depleted_material in openmc_materials:
+            nucvec = self.create_mass_percents_dictionary(depleted_material)
+            name = depleted_material.name
+            depleted_materials[name] = Materialflow(nucvec)
+            depleted_materials[name].mass = depleted_material.get_mass()
+            depleted_materials[name].vol = depleted_material.volume
+            depleted_materials[name].density = material.get_mass() / material.volume
+
+            if read_at_end:
+                starting_heavy_metal_mass = starting_material.fissionable_mass
+                depleted_heavy_metal_mass = depleted_material.fissionable_mass
+                power = resuts[1].source_rate
+                days = np.diff(results[1].time)[0] / _SECONDS_PER_DAY
+                burnup = power * days / (starting_heavy_metal_mass - depleted_heavy_metal_mass)
+            else:
+                burnup = 0
+            depleted_materials[name].burnup = burnup
+        return depleted_materials
+
+    def _create_mass_percents_dictionary(self, mat):
+        """Creates a dicitonary with nuclide codes
+        in zzaaam formate as keys, and material composition
+        in mass percent as values
+
+        Parameters
+        ----------
+        mat : openmc.Material
+            A material
+
+        Returns
+        -------
+            mass_dict : dict of int to float
+        """
+        at_percents = []
+        nucs = []
+        at_mass = []
+        for nuc, pt, tp in mat.nuclides:
+            nucs.append(nuc)
+            at_percents.append(pt)
+            at_mass.append(atomic_mass(nuc))
+
+        at_percents = np.array(at_percents)
+        at_mass = np.array(at_mass)
+
+        mass_percents = at_percents*at_mass / np.dot(at_percents, at_mass)
+        zai = list(map(openmc.data.zam, nucs))
+        zam = list(map(self._z_a_m_to_zam, zai))
+
+        return dict(zip(zam, mass_percents))
+
+    def _z_a_m_to_zam(self,z_a_m_tuple):
+        """Helper function for :func:`_create_mass_percents_dictionary`.
+        Converts an OpenMC (Z,A,M) tuple into a zzaaam nuclide code
+
+        Parameters
+        ----------
+        z_a_m_tuple : 3-tuple of int
+            (z, a, m)
+
+        Returns
+        -------
+            zam : int
+            Nuclide code in zzaaam format
+
+        """
+        z, a, m = z_a_m_tuple
+        zam = 1000 * z
+        zam += a
+        if m > 0 and (z != 95 and a != 242):
+            zam += 300 + 100 * m
+        elif z == 95 and a == 242:
+            if m == 0:
+              zam = 95642
+            else:
+              zam = 95242
+        return zam
+
 
     def run_depletion_step(self, mpi_args=None, threads=None):
         """Runs a depletion step in OpenMC as a subprocess

tests/integration_tests/file_interface_openmc/test.py

@@ -1,30 +1,20 @@
+"""Test OpenMC file interface"""
 import json
+from os import remove
 from pathlib import Path
 
 import numpy as np
 import pytest
 import openmc
 
-from saltproc import Reactor
-
 
 @pytest.fixture
 def geometry_switch(scope='module'):
     path = Path(__file__).parents[2]
     return (path / 'openmc_data' / 'geometry_switch.xml')
 
 
-@pytest.fixture
-def msr(scope='module'):
-    reactor = Reactor(volume=1.0,
-                      power_levels=[1.250E+09, 1.250E+09, 5.550E+09],
-                      depletion_timesteps=[111.111, 2101.9, 3987.5],
-                      timestep_type='cumulative',
-                      timestep_units='d')
-    return reactor
-
-
-def test_write_runtime_input(openmc_depcode, msr):
+def test_write_runtime_input(openmc_depcode, openmc_reactor):
     # OpenMC
     input_materials = openmc.Materials.from_xml(
         openmc_depcode.template_input_file_path['materials'])
@@ -37,7 +27,7 @@ def test_write_runtime_input(openmc_depcode, msr):
     input_surfaces = input_geometry.get_all_surfaces()
     input_universes = input_geometry.get_all_universes()
 
-    openmc_depcode.write_runtime_input(msr,
+    openmc_depcode.write_runtime_input(openmc_reactor,
                                        0,
                                        False)
     # Load in the runtime_ objects
@@ -71,19 +61,19 @@ def test_write_runtime_input(openmc_depcode, msr):
     del input_materials, input_geometry
 
 
-def test_write_depletion_settings(openmc_depcode, msr):
+def test_write_depletion_settings(openmc_depcode, openmc_reactor):
     """
     Unit test for `Depcodeopenmc_depcode.write_depletion_settings`
     """
-    openmc_depcode.write_depletion_settings(msr, 0)
+    openmc_depcode.write_depletion_settings(openmc_reactor, 0)
     with open(openmc_depcode.output_path / 'depletion_settings.json') as f:
         j = json.load(f)
         assert Path(j['directory']).resolve() == Path(
             openmc_depcode.output_path)
-        assert j['timesteps'][0] == msr.depletion_timesteps[0]
+        assert j['timesteps'][0] == openmc_reactor.depletion_timesteps[0]
         assert j['operator_kwargs']['chain_file'] == \
             openmc_depcode.chain_file_path
-        assert j['integrator_kwargs']['power'] == msr.power_levels[0]
+        assert j['integrator_kwargs']['power'] == openmc_reactor.power_levels[0]
         assert j['integrator_kwargs']['timestep_units'] == 'd'
 
 
@@ -297,3 +287,30 @@ def _check_none_or_iterable_of_ndarray_equal(object1, object2):
             _check_none_or_iterable_of_ndarray_equal(subobject1, subobject2)
     else:
         assert object1 == object2
+
+def test_write_runtime_files(openmc_depcode, openmc_reactor):
+    ref_mats = openmc_depcode.read_depleted_materials(True)
+
+    # update_depletable_materials
+    openmc_depcode.update_depletable_materials(mats, 12.0)
+    file = openmc_depcode.runtime_matfile
+    test_mats = openmc.Materials.from_xml(openmc_depcode.runtime_matfile)
+
+    # compare material objects
+    ...
+
+    remove(openmc_depcode.runtime_matfile)
+
+    # write_runtime_input
+    openmc_depcode.write_runtime_input(openmc_reactor,
+                                        0,
+                                        False)
+
+    settings_file = openmc_depcode.runtime_inputfile['settings']
+    geometry_file = openmc_depcode.runtime_inputfile['geometry']
+
+    # compare settings and geometry files
+    ...
+
+    # switch_to_next_geometry
+    ...

tests/integration_tests/file_interface_serpent/test.py

@@ -5,8 +5,6 @@
 import numpy as np
 import pytest
 
-from saltproc import Reactor
-
 
 @pytest.fixture
 def geometry_switch(scope='module'):