mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
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.
This commit is contained in:
Andrew Johnson
parent
f3e7c7c623
commit
72bf884d9b
4 changed files with 99 additions and 85 deletions
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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):
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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):
|
||||
"""
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue