From 72bf884d9bdfe62af065642207fe7fc6b79df9c5 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 2 Aug 2019 17:00:38 -0500 Subject: [PATCH] Pass fission yields to depletion matrix_func Takes a single set of fission yields and passes them as an additional argument to matrix_func: >>> A = matrix_func(chain, rates, fission_yields) Applied to cf4, epc_rk4, celi, and leqi functions. Assumes that fission yields will not change during a depletion event. This change is probably overshadowed by how much the reaction rates may change, but still worth pointing out. --- openmc/deplete/_matrix_funcs.py | 83 +++++++++++++++++---------------- openmc/deplete/chain.py | 11 ++--- openmc/deplete/cram.py | 21 +++++++-- tests/dummy_operator.py | 69 ++++++++++++++------------- 4 files changed, 99 insertions(+), 85 deletions(-) diff --git a/openmc/deplete/_matrix_funcs.py b/openmc/deplete/_matrix_funcs.py index 1cacbd2c06..a606d739f7 100644 --- a/openmc/deplete/_matrix_funcs.py +++ b/openmc/deplete/_matrix_funcs.py @@ -1,77 +1,78 @@ """Functions to form the special matrix for depletion""" -def celi_f1(chain, rates): - return (5 / 12 * chain.form_matrix(rates[0]) - + 1 / 12 * chain.form_matrix(rates[1])) +def celi_f1(chain, rates, fission_yields=None): + return (5 / 12 * chain.form_matrix(rates[0], fission_yields) + + 1 / 12 * chain.form_matrix(rates[1], fission_yields)) -def celi_f2(chain, rates): - return (1 / 12 * chain.form_matrix(rates[0]) - + 5 / 12 * chain.form_matrix(rates[1])) +def celi_f2(chain, rates, fission_yields=None): + return (1 / 12 * chain.form_matrix(rates[0], fission_yields) + + 5 / 12 * chain.form_matrix(rates[1], fission_yields)) -def cf4_f1(chain, rates): - return 1 / 2 * chain.form_matrix(rates) +def cf4_f1(chain, rates, fission_yields=None): + return 1 / 2 * chain.form_matrix(rates, fission_yields) -def cf4_f2(chain, rates): - return -1 / 2 * chain.form_matrix(rates[0]) + chain.form_matrix(rates[1]) +def cf4_f2(chain, rates, fission_yields=None): + return (-1 / 2 * chain.form_matrix(rates[0], fission_yields) + + chain.form_matrix(rates[1], fission_yields)) -def cf4_f3(chain, rates): - return (1 / 4 * chain.form_matrix(rates[0]) - + 1 / 6 * chain.form_matrix(rates[1]) - + 1 / 6 * chain.form_matrix(rates[2]) - - 1 / 12 * chain.form_matrix(rates[3])) +def cf4_f3(chain, rates, fission_yields=None): + return (1 / 4 * chain.form_matrix(rates[0], fission_yields) + + 1 / 6 * chain.form_matrix(rates[1], fission_yields) + + 1 / 6 * chain.form_matrix(rates[2], fission_yields) + - 1 / 12 * chain.form_matrix(rates[3], fission_yields)) -def cf4_f4(chain, rates): - return (-1 / 12 * chain.form_matrix(rates[0]) - + 1 / 6 * chain.form_matrix(rates[1]) - + 1 / 6 * chain.form_matrix(rates[2]) - + 1 / 4 * chain.form_matrix(rates[3])) +def cf4_f4(chain, rates, fission_yields=None): + return (-1 / 12 * chain.form_matrix(rates[0], fission_yields) + + 1 / 6 * chain.form_matrix(rates[1], fission_yields) + + 1 / 6 * chain.form_matrix(rates[2], fission_yields) + + 1 / 4 * chain.form_matrix(rates[3], fission_yields)) -def rk4_f1(chain, rates): - return 1 / 2 * chain.form_matrix(rates) +def rk4_f1(chain, rates, fission_yields=None): + return 1 / 2 * chain.form_matrix(rates, fission_yields) -def rk4_f4(chain, rates): - return (1 / 6 * chain.form_matrix(rates[0]) - + 1 / 3 * chain.form_matrix(rates[1]) - + 1 / 3 * chain.form_matrix(rates[2]) - + 1 / 6 * chain.form_matrix(rates[3])) +def rk4_f4(chain, rates, fission_yields=None): + return (1 / 6 * chain.form_matrix(rates[0], fission_yields) + + 1 / 3 * chain.form_matrix(rates[1], fission_yields) + + 1 / 3 * chain.form_matrix(rates[2], fission_yields) + + 1 / 6 * chain.form_matrix(rates[3], fission_yields)) -def leqi_f1(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) +def leqi_f1(chain, inputs, fission_yields): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) dt_l, dt = inputs[2], inputs[3] return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 -def leqi_f2(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) +def leqi_f2(chain, inputs, fission_yields=None): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) dt_l, dt = inputs[2], inputs[3] return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 -def leqi_f3(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - f3 = chain.form_matrix(inputs[2]) +def leqi_f3(chain, inputs, fission_yields=None): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) + f3 = chain.form_matrix(inputs[2], fission_yields) dt_l, dt = inputs[3], inputs[4] return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1 + (dt ** 2 + 6 * dt * dt_l + 5 * dt_l ** 2) / (12 * dt_l * (dt + dt_l)) * f2 + dt_l / (12 * (dt + dt_l)) * f3) -def leqi_f4(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - f3 = chain.form_matrix(inputs[2]) +def leqi_f4(chain, inputs, fission_yields=None): + f1 = chain.form_matrix(inputs[0], fission_yields) + f2 = chain.form_matrix(inputs[1], fission_yields) + f3 = chain.form_matrix(inputs[2], fission_yields) dt_l, dt = inputs[3], inputs[4] return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1 + (dt ** 2 + 2 * dt * dt_l + dt_l ** 2) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 865f533ceb..62158255a7 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -129,7 +129,6 @@ class Chain(object): self.nuclides = [] self.reactions = [] self.nuclide_dict = OrderedDict() - self._default_fsn_yields = None def __contains__(self, nuclide): return nuclide in self.nuclide_dict @@ -375,7 +374,7 @@ class Chain(object): clean_indentation(root_elem) tree.write(str(filename), encoding='utf-8') - def _build_default_yields(self): + def get_thermal_fission_yields(self): """Return dictionary {str: {str: float}}""" # Take lowest energy for back compatability # Should be removed by end of this feature @@ -383,8 +382,8 @@ class Chain(object): for nuc in self.nuclides: if len(nuc.yield_data) == 0: continue - _energy, yield_data = sorted(nuc.yield_data.items())[0] - out[nuc.name] = {prod: frac for prod, frac in yield_data} + yield_obj = nuc.yield_data[min(nuc.yield_energies)] + out[nuc.name] = dict(yield_obj) return out def form_matrix(self, rates, fission_yields=None): @@ -407,9 +406,7 @@ class Chain(object): reactions = set() if fission_yields is None: - if self._default_fsn_yields is None: - self._default_fsn_yields = self._build_default_yields() - fission_yields = self._default_fsn_yields + fission_yields = self.get_thermal_fission_yields() for i, nuc in enumerate(self.nuclides): diff --git a/openmc/deplete/cram.py b/openmc/deplete/cram.py index 41a4e4c710..d0c868004c 100644 --- a/openmc/deplete/cram.py +++ b/openmc/deplete/cram.py @@ -30,7 +30,9 @@ def deplete(chain, x, rates, dt, matrix_func=None): dt : float Time in [s] to deplete for maxtrix_func : Callable, optional - Function of two variables: ``chain`` and ``rates``. + Function to form the depletion matrix after calling + ``matrix_func(chain, rates, fission_yields)``, where + ``fission_yields = {parent: {product: yield_frac}}`` Expected to return the depletion matrix required by :func:`CRAM48`. @@ -40,9 +42,15 @@ def deplete(chain, x, rates, dt, matrix_func=None): Updated atom number vectors for each material """ + if not hasattr(chain, "fission_yields"): + fission_yields = repeat(chain.get_thermal_fission_yields()) + else: + fission_yields = chain.fission_yields + # Use multiprocessing pool to distribute work with Pool() as pool: - iters = zip(repeat(chain), x, rates, repeat(dt), repeat(matrix_func)) + iters = zip(repeat(chain), x, rates, repeat(dt), + fission_yields, repeat(matrix_func)) x_result = list(pool.starmap(_cram_wrapper, iters)) return x_result @@ -67,7 +75,7 @@ def timed_deplete(*args, **kwargs): return time.time() - start, results -def _cram_wrapper(chain, n0, rates, dt, matrix_func=None): +def _cram_wrapper(chain, n0, rates, dt, fission_yields, matrix_func=None): """Wraps depletion matrix creation / CRAM solve for multiprocess execution Parameters @@ -82,6 +90,9 @@ def _cram_wrapper(chain, n0, rates, dt, matrix_func=None): Time to integrate to. maxtrix_func : function, optional Function to form the depletion matrix + fission_yields : dict + Single-energy fission yields of the form + ``{parent: {product: fission_yield}}`` Returns ------- @@ -90,9 +101,9 @@ def _cram_wrapper(chain, n0, rates, dt, matrix_func=None): """ if matrix_func is None: - A = chain.form_matrix(rates) + A = chain.form_matrix(rates, fission_yields) else: - A = matrix_func(chain, rates) + A = matrix_func(chain, rates, fission_yields) return CRAM48(A, n0, dt) diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py index c23abef631..acfdc5b583 100644 --- a/tests/dummy_operator.py +++ b/tests/dummy_operator.py @@ -6,6 +6,42 @@ from openmc.deplete.reaction_rates import ReactionRates from openmc.deplete.abc import TransportOperator, OperatorResult +class TestChain(object): + + @staticmethod + def get_thermal_fission_yields(): + return None + + def form_matrix(self, rates, _fission_yields=None): + """Forms the f(y) matrix in y' = f(y)y. + + Nominally a depletion matrix, this is abstracted on the off chance + that the function f has nothing to do with depletion at all. + + Parameters + ---------- + rates : numpy.ndarray + Slice of reaction rates for a single material + _fission_yields : optional + Not used + + Returns + ------- + scipy.sparse.csr_matrix + Sparse matrix representing f(y). + """ + + y_1 = rates[0, 0] + y_2 = rates[1, 0] + + a11 = np.sin(y_2) + a12 = np.cos(y_1) + a21 = -np.cos(y_2) + a22 = np.sin(y_1) + + return sp.csr_matrix(np.array([[a11, a12], [a21, a22]])) + + class DummyOperator(TransportOperator): """This is a dummy operator class with no statistical uncertainty. @@ -21,6 +57,7 @@ class DummyOperator(TransportOperator): """ def __init__(self, previous_results=None): self.prev_res = previous_results + self.chain = TestChain() def __call__(self, vec, power, print_out=False): """Evaluates F(y) @@ -52,38 +89,6 @@ class DummyOperator(TransportOperator): # Create a fake rates object return OperatorResult(ufloat(0.0, 0.0), reaction_rates) - @property - def chain(self): - return self - - def form_matrix(self, rates): - """Forms the f(y) matrix in y' = f(y)y. - - Nominally a depletion matrix, this is abstracted on the off chance - that the function f has nothing to do with depletion at all. - - Parameters - ---------- - rates : numpy.ndarray - Slice of reaction rates for a single material - - Returns - ------- - scipy.sparse.csr_matrix - Sparse matrix representing f(y). - """ - - y_1 = rates[0, 0] - y_2 = rates[1, 0] - - mat = np.zeros((2, 2)) - a11 = np.sin(y_2) - a12 = np.cos(y_1) - a21 = -np.cos(y_2) - a22 = np.sin(y_1) - - return sp.csr_matrix(np.array([[a11, a12], [a21, a22]])) - @property def volume(self): """