From 98957fb47350d8b012236f00cc89a030c38e4590 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 26 Jul 2019 10:31:24 -0500 Subject: [PATCH] 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 | 4 ++-- openmc/deplete/integrator/abc.py | 12 +++++++----- openmc/deplete/integrator/si_celi.py | 4 +++- openmc/deplete/integrator/si_leqi.py | 4 +++- openmc/deplete/operator.py | 3 ++- openmc/deplete/results.py | 6 +++--- tests/dummy_operator.py | 4 +++- tests/unit_tests/test_deplete_integrator.py | 6 ++++-- 8 files changed, 27 insertions(+), 16 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 1be3e27ff9..5e2f1dc64f 100644 --- a/openmc/deplete/abc.py +++ b/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 diff --git a/openmc/deplete/integrator/abc.py b/openmc/deplete/integrator/abc.py index 89015ea2df..07151d13b4 100644 --- a/openmc/deplete/integrator/abc.py +++ b/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: diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py index 26087c1f2c..872863ef22 100644 --- a/openmc/deplete/integrator/si_celi.py +++ b/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 diff --git a/openmc/deplete/integrator/si_leqi.py b/openmc/deplete/integrator/si_leqi.py index 7d01782bae..86b184038f 100644 --- a/openmc/deplete/integrator/si_leqi.py +++ b/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 diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index a38c4e5016..e9c41093aa 100644 --- a/openmc/deplete/operator.py +++ b/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 diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index 3771f9683b..a2462b80a1 100644 --- a/openmc/deplete/results.py +++ b/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 diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py index c66070ad40..23dc22288a 100644 --- a/tests/dummy_operator.py +++ b/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): diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index 40090264ae..c50ac66a94 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/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)