Merge pull request #1279 from drewejohnson/enh-dep-eig-unc

Store uncertainty on k through depletion

docs/source/io_formats/depletion_results.rst

@@ -12,9 +12,10 @@ The current version of the depletion results file format is 1.0.
              - **version** (*int[2]*) -- Major and minor version of the
                statepoint file format.
 
-:Datasets: - **eigenvalues** (*double[][]*) -- k-eigenvalues at each
+:Datasets: - **eigenvalues** (*double[][][2]*) -- k-eigenvalues at each
              time/stage. This array has shape (number of timesteps, number of
-             stages).
+             stages, value). The last axis contains the eigenvalue and the 
+             associated uncertainty
            - **number** (*double[][][][]*) -- Total number of atoms. This array
              has shape (number of timesteps, number of stages, number of
              materials, number of nuclides).

openmc/deplete/operator.py

@@ -525,7 +525,8 @@ class Operator(TransportOperator):
         rates = self.reaction_rates
         rates[:, :, :] = 0.0
 
-        k_combined = openmc.capi.keff()[0]
+        # Get k and uncertainty
+        k_combined = openmc.capi.keff()
 
         # Extract tally bins
         materials = self.burnable_mats

openmc/deplete/results.py

@@ -237,8 +237,8 @@ class Results(object):
                               chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn),
                               dtype='float64')
 
-        handle.create_dataset("eigenvalues", (1, n_stages),
-                              maxshape=(None, n_stages), dtype='float64')
+        handle.create_dataset("eigenvalues", (1, n_stages, 2),
+                              maxshape=(None, n_stages, 2), dtype='float64')
 
         handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
 

openmc/deplete/results_list.py

@@ -91,11 +91,13 @@ class ResultsList(list):
         time : numpy.ndarray
             Array of times in [s]
         eigenvalue : numpy.ndarray
-            k-eigenvalue at each time
+            k-eigenvalue at each time. Column 0
+            contains the eigenvalue, while column
+            1 contains the associated uncertainty
 
         """
         time = np.empty_like(self, dtype=float)
-        eigenvalue = np.empty_like(self, dtype=float)
+        eigenvalue = np.empty((len(self), 2), dtype=float)
 
         # Get time/eigenvalue at each point
         for i, result in enumerate(self):

tests/unit_tests/test_deplete_integrator.py

@@ -64,8 +64,9 @@ def test_results_save(run_in_tmpdir):
         r1[:] = np.random.rand(2, 2, 2)
 
     # Create global terms
-    eigvl1 = np.random.rand(stages)
-    eigvl2 = np.random.rand(stages)
+    # Col 0: eig, Col 1: uncertainty
+    eigvl1 = np.random.rand(stages, 2)
+    eigvl2 = np.random.rand(stages, 2)
 
     eigvl1 = comm.bcast(eigvl1, root=0)
     eigvl2 = comm.bcast(eigvl2, root=0)

tests/unit_tests/test_deplete_resultslist.py

@@ -44,6 +44,8 @@ def test_get_eigenvalue(res):
 
     t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
     k_ref = [1.16984322, 1.19097427, 1.03012572, 1.20045627]
+    u_ref = [0.0375587, 0.0347639, 0.07216021, 0.02839642]
 
     np.testing.assert_allclose(t, t_ref)
-    np.testing.assert_allclose(k, k_ref)
+    np.testing.assert_allclose(k[:, 0], k_ref)
+    np.testing.assert_allclose(k[:, 1], u_ref)