getting initial nuclide densities from depletion results, other choice would be to update python densities, but depletion results are needed anyway

examples/python/pincell_depletion/restart_depletion.py

@@ -40,7 +40,7 @@ previous_results = openmc.deplete.ResultsList("depletion_results.h5")
 #                      Transport calculation settings
 ###############################################################################
 
-# Instantiate a Settings object, set all runtime parameters, and export to XML
+# Instantiate a Settings object, set all runtime parameters
 settings_file = openmc.Settings()
 settings_file.batches = batches
 settings_file.inactive = inactive

openmc/deplete/integrator/predictor.py

@@ -55,7 +55,7 @@ def predictor(operator, timesteps, power, print_out=True):
         else:
             i_res = len(operator.prev_res.get_eigenvalue()[0])
 
-        #TODO : Get last time step power from previous results, and run a 
+        #TODO : Get last time step power from previous results, and run a
         # new calculation if different from power at the first time step
         # If no TH coupling, just re-scale rates by ratio of power
 
@@ -68,7 +68,6 @@ def predictor(operator, timesteps, power, print_out=True):
             # Avoid doing first run if already done in previous calculation
             if i > 0 or operator.prev_res == None:
                 op_results = [operator(x[0], p)]
-
             else:
                 op_results = [operator.prev_res[-1]]
                 op_results[0].rates = op_results[0].rates[0]

openmc/deplete/operator.py

@@ -77,7 +77,7 @@ class Operator(TransportOperator):
         OpenMC settings object
     dilute_initial : float
         Initial atom density to add for nuclides that are zero in initial
-        condition to ensure they exist in the decay chain.  Only done for
+        condition to ensure they exist in the decay chain. Only done for
         nuclides with reaction rates. Defaults to 1.0e3.
     output_dir : pathlib.Path
         Path to output directory to save results.
@@ -112,6 +112,9 @@ class Operator(TransportOperator):
 
             # Store previous results in operator
             self.prev_res = prev_results
+
+            # Get number densities from previous results
+            #self.number = prev_results
         else:
             self.prev_res = None
 
@@ -125,8 +128,9 @@ class Operator(TransportOperator):
         self._burnable_nucs = [nuc for nuc in self.nuclides_with_data
                                if nuc in self.chain]
 
-        # Extract number densities from the geometry
-        self._extract_number(self.local_mats, volume, nuclides)
+        # Extract number densities from the geometry / previous depletion run
+        self._extract_number(self.local_mats, volume, nuclides, \
+                                                self.prev_res)
 
         # Create reaction rates array
         self.reaction_rates = ReactionRates(
@@ -225,7 +229,7 @@ class Operator(TransportOperator):
 
         return burnable_mats, volume, nuclides
 
-    def _extract_number(self, local_mats, volume, nuclides):
+    def _extract_number(self, local_mats, volume, nuclides, prev_res=None):
         """Construct AtomNumber using geometry
 
         Parameters
@@ -244,10 +248,18 @@ class Operator(TransportOperator):
             for nuc in self._burnable_nucs:
                 self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial)
 
-        # Now extract the number densities and store
-        for mat in self.geometry.get_all_materials().values():
-            if str(mat.id) in local_mats:
-                self._set_number_from_mat(mat)
+        # Now extract and store the number densities
+        # From the geometry if no previous depletion results
+        if prev_res == None:
+            for mat in self.geometry.get_all_materials().values():
+                if str(mat.id) in local_mats:
+                    self._set_number_from_mat(mat)
+
+        # Else from previous depletion results
+        else:
+            for mat in self.geometry.get_all_materials().values():
+                if str(mat.id) in local_mats:
+                    self._set_number_from_results(mat, prev_res)
 
     def _set_number_from_mat(self, mat):
         """Extracts material and number densities from openmc.Material
@@ -262,6 +274,23 @@ class Operator(TransportOperator):
 
         for nuclide, density in mat.get_nuclide_atom_densities().values():
             number = density * 1.0e24
+            print(nuclide)
+            self.number.set_atom_density(mat_id, nuclide, number)
+
+    def _set_number_from_results(self, mat, prev_res):
+        """Extracts material and number densities from previous results
+
+        Parameters
+        ----------
+        mat : openmc.Material
+            The material to read from
+
+        """
+        mat_id = str(mat.id)
+
+        for nuclide, density in mat.get_nuclide_atom_densities().values():
+            number = density * 1.0e24
+            print(nuclide, number)
             self.number.set_atom_density(mat_id, nuclide, number)
 
     def initial_condition(self):
@@ -325,9 +354,13 @@ class Operator(TransportOperator):
                                       " is negative (density = ", val, " at/barn-cm)")
                             number_i[mat, nuc] = 0.0
 
+                # Update densities on C API side
                 mat_internal = openmc.capi.materials[int(mat)]
                 mat_internal.set_densities(nuclides, densities)
 
+                #TODO Update densities on the Python side, otherwise the
+                # summary.h5 file contains densities at the first time step
+
     def _generate_materials_xml(self):
         """Creates materials.xml from self.number.