Store OperatorResults.k as uncertainties.ufloat

This allows the propagation of the uncertainty through
the SIE iterations, and alows the expressions for updating
k to be preserved. The OperatorResult.k was previous stored
as a tuple of two floats. The interface into the depletion_results
file is not changed, as the Results.save method, which writes the
OperatorResult data, breaks the ufloat into k and uncertainties.

Documentation has been updated on the OperatorResult.

openmc/deplete/abc.py

@@ -20,8 +20,8 @@ Result of applying transport operator
 
 Parameters
 ----------
-k : float
-    Resulting eigenvalue
+k : uncertainties.ufloat
+    Resulting eigenvalue and standard deviation
 rates : openmc.deplete.ReactionRates
     Resulting reaction rates
 

openmc/deplete/integrator/abc.py

@@ -2,7 +2,9 @@ from copy import deepcopy
 from abc import ABC, abstractmethod
 from collections.abc import Iterable
 
-from openmc.deplete import Results
+from uncertainties import ufloat
+
+from openmc.deplete import Results, OperatorResult
 
 
 class Integrator(ABC):
@@ -96,12 +98,12 @@ class Integrator(ABC):
     def _get_bos_data_from_restart(self, step_index, step_power, bos_conc):
         res = self.operator.prev_res[-1]
         bos_conc = res.data[0]
-        res.rates = res.rates[0]
-        res.k = res.k[0]
+        rates = res.rates[0]
+        k = ufloat(res.k[0,0], res.k[0, 1])
 
         # Scale rates by ratio of powers
-        res.rates *= step_power / res.power[0]
-        return bos_conc, res
+        rates *= step_power / res.power[0]
+        return bos_conc, OperatorResult(k, rates)
 
     def _get_start_data(self):
         if self.operator.prev_res is None:

openmc/deplete/integrator/si_celi.py

@@ -3,6 +3,8 @@
 import copy
 from collections.abc import Iterable
 
+from uncertainties import ufloat
+
 from .cram import timed_deplete
 from ..results import Results
 from ..abc import OperatorResult
@@ -82,7 +84,7 @@ def si_celi(operator, timesteps, power=None, power_density=None,
             op_results[0].rates = op_results[0].rates[0]
 
             # Set first stage value of keff
-            op_results[0].k = op_results[0].k[0]
+            op_results[0].k = ufloat(*op_results[0].k[0])
 
             # Scale reaction rates by ratio of powers
             power_res = operator.prev_res[-1].power

openmc/deplete/integrator/si_leqi.py

@@ -4,6 +4,8 @@ import copy
 from collections.abc import Iterable
 from itertools import repeat
 
+from uncertainties import ufloat
+
 from .si_celi import si_celi_inner
 from .leqi import _leqi_f1, _leqi_f2, _leqi_f3, _leqi_f4
 from .cram import timed_deplete
@@ -84,7 +86,7 @@ def si_leqi(operator, timesteps, power=None, power_density=None,
             op_results[0].rates = op_results[0].rates[0]
 
             # Set first stage value of keff
-            op_results[0].k = op_results[0].k[0]
+            op_results[0].k = ufloat(*op_results[0].k[0])
 
             # Scale reaction rates by ratio of powers
             power_res = operator.prev_res[-1].power

openmc/deplete/operator.py

@@ -16,6 +16,7 @@ import xml.etree.ElementTree as ET
 
 import h5py
 import numpy as np
+from uncertainties import ufloat
 
 import openmc
 import openmc.capi
@@ -526,7 +527,7 @@ class Operator(TransportOperator):
         rates[:, :, :] = 0.0
 
         # Get k and uncertainty
-        k_combined = openmc.capi.keff()
+        k_combined = ufloat(*openmc.capi.keff())
 
         # Extract tally bins
         materials = self.burnable_mats

openmc/deplete/results.py

@@ -21,8 +21,8 @@ class Results(object):
 
     Attributes
     ----------
-    k : list of float
-        Eigenvalue for each substep.
+    k : list of (float, float)
+        Eigenvalue and uncertainty for each substep.
     time : list of float
         Time at beginning, end of step, in seconds.
     power : float
@@ -446,7 +446,7 @@ class Results(object):
             for mat_i in range(n_mat):
                 results[i, mat_i, :] = x[offset + i][mat_i][:]
 
-        results.k = [r.k for r in op_results]
+        results.k = [(r.k.nominal_value, r.k.std_dev) for r in op_results]
         results.rates = [r.rates for r in op_results]
         results.time = t
         results.power = power

tests/dummy_operator.py

@@ -1,5 +1,7 @@
 import numpy as np
 import scipy.sparse as sp
+from uncertainties import ufloat
+
 from openmc.deplete.reaction_rates import ReactionRates
 from openmc.deplete.abc import TransportOperator, OperatorResult
 
@@ -48,7 +50,7 @@ class DummyOperator(TransportOperator):
         reaction_rates[0, 1, 0] = vec[0][1]
 
         # Create a fake rates object
-        return OperatorResult(0.0, reaction_rates)
+        return OperatorResult(ufloat(0.0, 0.0), reaction_rates)
 
     @property
     def chain(self):

tests/unit_tests/test_deplete_integrator.py

@@ -11,6 +11,8 @@ import os
 from unittest.mock import MagicMock
 
 import numpy as np
+from uncertainties import ufloat
+
 from openmc.deplete import (ReactionRates, Results, ResultsList, comm,
                             OperatorResult)
 
@@ -74,8 +76,8 @@ def test_results_save(run_in_tmpdir):
     t1 = [0.0, 1.0]
     t2 = [1.0, 2.0]
 
-    op_result1 = [OperatorResult(k, rates) for k, rates in zip(eigvl1, rate1)]
-    op_result2 = [OperatorResult(k, rates) for k, rates in zip(eigvl2, rate2)]
+    op_result1 = [OperatorResult(ufloat(*k), rates) for k, rates in zip(eigvl1, rate1)]
+    op_result2 = [OperatorResult(ufloat(*k), rates) for k, rates in zip(eigvl2, rate2)]
     Results.save(op, x1, op_result1, t1, 0, 0)
     Results.save(op, x2, op_result2, t2, 0, 1)