Merge pull request #1946 from paulromano/depletion-bugfixes

Depletion-related bug fixes

examples/pincell_depletion/restart_depletion.py

@@ -39,9 +39,11 @@ settings.entropy_mesh = entropy_mesh
 #                   Initialize and run depletion calculation
 ###############################################################################
 
+model = openmc.Model(geometry=geometry, settings=settings)
+
 # Create depletion "operator"
-chain_file = './chain_simple.xml'
-op = openmc.deplete.Operator(geometry, settings, chain_file, previous_results)
+chain_file = 'chain_simple.xml'
+op = openmc.deplete.Operator(model, chain_file, previous_results)
 
 # Perform simulation using the predictor algorithm
 time_steps = [1.0, 1.0, 1.0, 1.0, 1.0]  # days
@@ -70,21 +72,20 @@ time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
 ###############################################################################
 
 days = 24*60*60
-plt.figure()
-plt.plot(time/days, keff, label="K-effective")
-plt.xlabel("Time (days)")
-plt.ylabel("Keff")
+fig, ax = plt.subplots()
+ax.errorbar(time/days, keff[:, 0], keff[:, 1], label="K-effective")
+ax.set_xlabel("Time [d]")
+ax.set_ylabel("Keff")
 plt.show()
 
-plt.figure()
-plt.plot(time/days, n_U235, label="U 235")
-plt.xlabel("Time (days)")
-plt.ylabel("n U5 (-)")
+fig, ax = plt.subplots()
+ax.plot(time/days, n_U235, label="U235")
+ax.set_xlabel("Time [d]")
+ax.set_ylabel("U235 atoms")
 plt.show()
 
-plt.figure()
-plt.plot(time/days, Xe_capture, label="Xe135 capture")
-plt.xlabel("Time (days)")
-plt.ylabel("RR (-)")
+fig, ax = plt.subplots()
+ax.plot(time/days, Xe_capture, label="Xe135 capture")
+ax.set_xlabel("Time [d]")
+ax.set_ylabel("Xe135 capture rate")
 plt.show()
-plt.close('all')

examples/pincell_depletion/run_depletion.py

@@ -87,7 +87,7 @@ settings.entropy_mesh = entropy_mesh
 model = openmc.Model(geometry=geometry, settings=settings)
 
 # Create depletion "operator"
-chain_file = './chain_simple.xml'
+chain_file = 'chain_simple.xml'
 op = openmc.deplete.Operator(model, chain_file)
 
 # Perform simulation using the predictor algorithm
@@ -117,21 +117,20 @@ time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
 ###############################################################################
 
 days = 24*60*60
-plt.figure()
-plt.plot(time/days, keff, label="K-effective")
-plt.xlabel("Time (days)")
-plt.ylabel("Keff")
+fig, ax = plt.subplots()
+ax.errorbar(time/days, keff[:, 0], keff[:, 1], label="K-effective")
+ax.set_xlabel("Time [d]")
+ax.set_ylabel("Keff")
 plt.show()
 
-plt.figure()
-plt.plot(time/days, n_U235, label="U235")
-plt.xlabel("Time (days)")
-plt.ylabel("n U5 (-)")
+fig, ax = plt.subplots()
+ax.plot(time/days, n_U235, label="U235")
+ax.set_xlabel("Time [d]")
+ax.set_ylabel("U235 atoms")
 plt.show()
 
-plt.figure()
-plt.plot(time/days, Xe_capture, label="Xe135 capture")
-plt.xlabel("Time (days)")
-plt.ylabel("RR (-)")
+fig, ax = plt.subplots()
+ax.plot(time/days, Xe_capture, label="Xe135 capture")
+ax.set_xlabel("Time [d]")
+ax.set_ylabel("Xe135 capture rate")
 plt.show()
-plt.close('all')

openmc/deplete/operator.py

@@ -309,7 +309,7 @@ class Operator(TransportOperator):
         if normalization_mode == "fission-q":
             self._normalization_helper = ChainFissionHelper()
         elif normalization_mode == "energy-deposition":
-            score = "heating" if settings.photon_transport else "heating-local"
+            score = "heating" if self.settings.photon_transport else "heating-local"
             self._normalization_helper = EnergyScoreHelper(score)
         else:
             self._normalization_helper = SourceRateHelper()
@@ -642,7 +642,7 @@ class Operator(TransportOperator):
         if mfile.exists():
             tree = ET.parse(str(mfile))
             xs = tree.find('cross_sections')
-            if xs is not None:
+            if xs is not None and self.materials.cross_sections is None:
                 self.materials.cross_sections = xs.text
 
         self.materials.export_to_xml()

openmc/deplete/results_list.py

@@ -361,6 +361,7 @@ class ResultsList(list):
             mat_id = str(mat.id)
             if mat_id in result.mat_to_ind:
                 mat.volume = result.volume[mat_id]
+                mat.set_density('sum')
                 for nuc in result.nuc_to_ind:
                     if nuc not in available_cross_sections:
                         continue

openmc/model/model.py

@@ -1,4 +1,5 @@
 from collections.abc import Iterable
+from functools import lru_cache
 import os
 from pathlib import Path
 from numbers import Integral
@@ -88,30 +89,6 @@ class Model:
         if plots is not None:
             self.plots = plots
 
-        # Store dictionaries to the materials and cells by ID and names
-        if materials is None:
-            mats = self.geometry.get_all_materials().values()
-        else:
-            mats = self.materials
-        cells = self.geometry.get_all_cells()
-        # Get the ID maps
-        self._materials_by_id = {mat.id: mat for mat in mats}
-        self._cells_by_id = {cell.id: cell for cell in cells.values()}
-
-        # Get the names maps, but since names are not unique, store a list for
-        # each name key. In this way when the user requests a change by a name,
-        # the change will be applied to all of the same name.
-        self._cells_by_name = {}
-        for cell in cells.values():
-            if cell.name not in self._cells_by_name:
-                self._cells_by_name[cell.name] = set()
-            self._cells_by_name[cell.name].add(cell)
-        self._materials_by_name = {}
-        for mat in mats:
-            if mat.name not in self._materials_by_name:
-                self._materials_by_name[mat.name] = set()
-            self._materials_by_name[mat.name].add(mat)
-
     @property
     def geometry(self):
         return self._geometry
@@ -140,6 +117,50 @@ class Model:
         except ImportError:
             return False
 
+    @property
+    @lru_cache(maxsize=None)
+    def _materials_by_id(self):
+        """Dictionary mapping material ID --> material"""
+        if self.materials is None:
+            mats = self.geometry.get_all_materials().values()
+        else:
+            mats = self.materials
+        return {mat.id: mat for mat in mats}
+
+    @property
+    @lru_cache(maxsize=None)
+    def _cells_by_id(self):
+        """Dictionary mapping cell ID --> cell"""
+        cells = self.geometry.get_all_cells()
+        return {cell.id: cell for cell in cells.values()}
+
+    @property
+    @lru_cache(maxsize=None)
+    def _cells_by_name(self):
+        # Get the names maps, but since names are not unique, store a set for
+        # each name key. In this way when the user requests a change by a name,
+        # the change will be applied to all of the same name.
+        result = {}
+        for cell in self.geometry.get_all_cells().values():
+            if cell.name not in result:
+                result[cell.name] = set()
+            result[cell.name].add(cell)
+        return result
+
+    @property
+    @lru_cache(maxsize=None)
+    def _materials_by_name(self):
+        if self.materials is None:
+            mats = self.geometry.get_all_materials().values()
+        else:
+            mats = self.materials
+        result = {}
+        for mat in mats:
+            if mat.name not in result:
+                result[mat.name] = set()
+            result[mat.name].add(mat)
+        return result
+
     @geometry.setter
     def geometry(self, geometry):
         check_type('geometry', geometry, openmc.Geometry)