From 7382b5d1c8b00342335b39159a00667ca8768b05 Mon Sep 17 00:00:00 2001 From: Lorenzo Chierici Date: Mon, 12 May 2025 22:48:21 +0200 Subject: [PATCH] External transfer rates source term (#3088) Co-authored-by: Paul Romano --- docs/source/pythonapi/deplete.rst | 5 +- openmc/deplete/abc.py | 110 ++++-- openmc/deplete/chain.py | 70 +++- openmc/deplete/integrators.py | 56 +-- openmc/deplete/pool.py | 59 ++- openmc/deplete/transfer_rates.py | 355 ++++++++++++++---- .../ref_depletion_with_ext_source.h5 | Bin 0 -> 37328 bytes .../ref_depletion_with_feed.h5 | Bin 37328 -> 37328 bytes .../ref_depletion_with_removal.h5 | Bin 37328 -> 37328 bytes .../ref_depletion_with_transfer.h5 | Bin 37328 -> 37328 bytes .../ref_no_depletion_only_feed.h5 | Bin 37328 -> 37328 bytes .../ref_no_depletion_only_removal.h5 | Bin 37328 -> 37328 bytes .../ref_no_depletion_with_ext_source.h5 | Bin 0 -> 37328 bytes .../ref_no_depletion_with_transfer.h5 | Bin 37328 -> 37328 bytes .../deplete_with_transfer_rates/test.py | 45 ++- .../test_deplete_external_source_rates.py | 170 +++++++++ .../unit_tests/test_deplete_transfer_rates.py | 65 ++-- 17 files changed, 751 insertions(+), 184 deletions(-) create mode 100644 tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_ext_source.h5 create mode 100644 tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_ext_source.h5 create mode 100644 tests/unit_tests/test_deplete_external_source_rates.py diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 3967f1a18..f112cf8cc 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -206,14 +206,15 @@ total system energy. The :class:`openmc.deplete.IndependentOperator` uses inner classes subclassed from those listed above to perform similar calculations. -The following classes are used to define transfer rates to model continuous -removal or feed of nuclides during depletion. +The following classes are used to define external source rates or transfer rates +to model continuous removal or feed of nuclides during depletion. .. autosummary:: :toctree: generated :nosignatures: :template: myclass.rst + transfer_rates.ExternalSourceRates transfer_rates.TransferRates Intermediate Classes diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index ab087b98e..c304b8dc6 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -29,7 +29,7 @@ from .results import Results, _SECONDS_PER_MINUTE, _SECONDS_PER_HOUR, \ _SECONDS_PER_DAY, _SECONDS_PER_JULIAN_YEAR from .pool import deplete from .reaction_rates import ReactionRates -from .transfer_rates import TransferRates +from .transfer_rates import TransferRates, ExternalSourceRates __all__ = [ @@ -607,9 +607,14 @@ class Integrator(ABC): next time step. Expected to be of the same shape as ``n0`` transfer_rates : openmc.deplete.TransferRates - Instance of TransferRates class to perform continuous transfer during depletion + Transfer rates for the depletion system used to model continuous + removal/feed between materials. .. versionadded:: 0.14.0 + external_source_rates : openmc.deplete.ExternalSourceRates + External source rates for the depletion system. + + .. versionadded:: 0.15.3 """ @@ -686,6 +691,7 @@ class Integrator(ABC): self.source_rates = np.asarray(source_rates) self.transfer_rates = None + self.external_source_rates = None if isinstance(solver, str): # Delay importing of cram module, which requires this file @@ -731,11 +737,11 @@ class Integrator(ABC): self._solver = func - def _timed_deplete(self, n, rates, dt, matrix_func=None): + def _timed_deplete(self, n, rates, dt, i=None, matrix_func=None): start = time.time() results = deplete( - self._solver, self.chain, n, rates, dt, matrix_func, - self.transfer_rates) + self._solver, self.chain, n, rates, dt, i, matrix_func, + self.transfer_rates, self.external_source_rates) return time.time() - start, results @abstractmethod @@ -885,13 +891,14 @@ class Integrator(ABC): self.operator.finalize() def add_transfer_rate( - self, - material: Union[str, int, Material], - components: Sequence[str], - transfer_rate: float, - transfer_rate_units: str = '1/s', - destination_material: Optional[Union[str, int, Material]] = None - ): + self, + material: str | int | Material, + components: Sequence[str], + transfer_rate: float, + transfer_rate_units: str = '1/s', + timesteps: Sequence[int] | None = None, + destination_material: str | int | Material | None = None + ): """Add transfer rates to depletable material. Parameters @@ -905,18 +912,79 @@ class Integrator(ABC): transfer_rate : float Rate at which elements are transferred. A positive or negative values set removal of feed rates, respectively. - destination_material : openmc.Material or str or int, Optional - Destination material to where nuclides get fed. transfer_rate_units : {'1/s', '1/min', '1/h', '1/d', '1/a'} Units for values specified in the transfer_rate argument. 's' means seconds, 'min' means minutes, 'h' means hours, 'a' means Julian years. + timesteps : list of int, optional + List of timestep indices where to set external source rates. + Defaults to None, which means the external source rate is set for + all timesteps. + destination_material : openmc.Material or str or int, Optional + Destination material to where nuclides get fed. """ if self.transfer_rates is None: - self.transfer_rates = TransferRates(self.operator, self.operator.model) + if hasattr(self.operator, 'model'): + materials = self.operator.model.materials + elif hasattr(self.operator, 'materials'): + materials = self.operator.materials + self.transfer_rates = TransferRates( + self.operator, materials, len(self.timesteps)) + + if self.external_source_rates is not None and destination_material: + raise ValueError('Currently is not possible to set a transfer rate ' + 'with destination matrial in combination with ' + 'external source rates.') + + self.transfer_rates.set_transfer_rate( + material, components, transfer_rate, transfer_rate_units, + timesteps, destination_material) + + def add_external_source_rate( + self, + material: str | int | Material, + composition: dict[str, float], + rate: float, + rate_units: str = 'g/s', + timesteps: Sequence[int] | None = None + ): + """Add external source rates to depletable material. + + Parameters + ---------- + material : openmc.Material or str or int + Depletable material + composition : dict of str to float + External source rate composition vector, where key can be an element + or a nuclide and value the corresponding weight percent. + rate : float + External source rate in units of mass per time. A positive or + negative value corresponds to a feed or removal rate, respectively. + units : {'g/s', 'g/min', 'g/h', 'g/d', 'g/a'} + Units for values specified in the `rate` argument. 's' for seconds, + 'min' for minutes, 'h' for hours, 'a' for Julian years. + timesteps : list of int, optional + List of timestep indices where to set external source rates. + Defaults to None, which means the external source rate is set for + all timesteps. + + """ + if self.external_source_rates is None: + if hasattr(self.operator, 'model'): + materials = self.operator.model.materials + elif hasattr(self.operator, 'materials'): + materials = self.operator.materials + self.external_source_rates = ExternalSourceRates( + self.operator, materials, len(self.timesteps)) + + if self.transfer_rates is not None and self.transfer_rates.index_transfer: + raise ValueError('Currently is not possible to set an external ' + 'source rate in combination with transfer rates ' + 'with destination matrial.') + + self.external_source_rates.set_external_source_rate( + material, composition, rate, rate_units, timesteps) - self.transfer_rates.set_transfer_rate(material, components, transfer_rate, - transfer_rate_units, destination_material) @add_params class SIIntegrator(Integrator): @@ -1047,10 +1115,10 @@ class SIIntegrator(Integrator): return inherited def integrate( - self, - output: bool = True, - path: PathLike = "depletion_results.h5" - ): + self, + output: bool = True, + path: PathLike = "depletion_results.h5" + ): """Perform the entire depletion process across all steps Parameters diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 4998a2c3d..9f234683a 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -703,7 +703,7 @@ class Chain: # Return CSC representation instead of DOK return matrix.tocsc() - def form_rr_term(self, tr_rates, mats): + def form_rr_term(self, tr_rates, current_timestep, mats): """Function to form the transfer rate term matrices. .. versionadded:: 0.14.0 @@ -712,6 +712,8 @@ class Chain: ---------- tr_rates : openmc.deplete.TransferRates Instance of openmc.deplete.TransferRates + current_timestep : int + Current timestep index mats : string or two-tuple of strings Two cases are possible: @@ -740,32 +742,74 @@ class Chain: for i, nuc in enumerate(self.nuclides): elm = re.split(r'\d+', nuc.name)[0] - # Build transfer terms matrices + # Build transfer terms (nuclide transfer only) if isinstance(mats, str): mat = mats - components = tr_rates.get_components(mat) + components = tr_rates.get_components(mat, current_timestep) + if not components: + break if elm in components: - matrix[i, i] = sum(tr_rates.get_transfer_rate(mat, elm)) + matrix[i, i] = sum( + tr_rates.get_external_rate(mat, elm, current_timestep)) elif nuc.name in components: - matrix[i, i] = sum(tr_rates.get_transfer_rate(mat, nuc.name)) + matrix[i, i] = sum( + tr_rates.get_external_rate(mat, nuc.name, current_timestep)) else: matrix[i, i] = 0.0 - #Build transfer terms matrices + + # Build transfer terms (transfer from one material into another) elif isinstance(mats, tuple): dest_mat, mat = mats - if dest_mat in tr_rates.get_destination_material(mat, elm): - dest_mat_idx = tr_rates.get_destination_material(mat, elm).index(dest_mat) - matrix[i, i] = tr_rates.get_transfer_rate(mat, elm)[dest_mat_idx] - elif dest_mat in tr_rates.get_destination_material(mat, nuc.name): - dest_mat_idx = tr_rates.get_destination_material(mat, nuc.name).index(dest_mat) - matrix[i, i] = tr_rates.get_transfer_rate(mat, nuc.name)[dest_mat_idx] + components = tr_rates.get_components(mat, current_timestep, dest_mat) + if elm in components: + matrix[i, i] = tr_rates.get_external_rate( + mat, elm, current_timestep, dest_mat)[0] + elif nuc.name in components: + matrix[i, i] = tr_rates.get_external_rate( + mat, nuc.name, current_timestep, dest_mat)[0] else: matrix[i, i] = 0.0 - #Nothing else is allowed # Return CSC instead of DOK return matrix.tocsc() + def form_ext_source_term(self, ext_source_rates, current_timestep, mat): + """Function to form the external source rate term vectors. + + .. versionadded:: 0.15.3 + + Parameters + ---------- + ext_source_rates : openmc.deplete.ExternalSourceRates + Instance of openmc.deplete.ExternalSourceRates + current_timestep : int + Current timestep index + mat : string + Material id + + Returns + ------- + scipy.sparse.csc_matrix + Sparse vector representing external source term. + + """ + if not ext_source_rates.get_components(mat, current_timestep): + return + # Use DOK as intermediate representation + n = len(self) + vector = sp.dok_matrix((n, 1)) + + for i, nuc in enumerate(self.nuclides): + # Build source term vector + if nuc.name in ext_source_rates.get_components(mat, current_timestep): + vector[i] = sum(ext_source_rates.get_external_rate( + mat, nuc.name, current_timestep)) + else: + vector[i] = 0.0 + + # Return CSC instead of DOK + return vector.tocsc() + def get_branch_ratios(self, reaction="(n,gamma)"): """Return a dictionary with reaction branching ratios diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index 50810c88a..000cb2c41 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -40,8 +40,8 @@ class PredictorIntegrator(Integrator): Time in [s] for the entire depletion interval source_rate : float Power in [W] or source rate in [neutron/sec] - _i : int or None - Iteration index. Not used + _i : int, optional + Current iteration count. Not used Returns ------- @@ -54,7 +54,7 @@ class PredictorIntegrator(Integrator): with predictor """ - proc_time, n_end = self._timed_deplete(n, rates, dt) + proc_time, n_end = self._timed_deplete(n, rates, dt, _i) return proc_time, [n_end], [] @@ -106,12 +106,12 @@ class CECMIntegrator(Integrator): Eigenvalue and reaction rates from transport simulations """ # deplete across first half of interval - time0, n_middle = self._timed_deplete(n, rates, dt / 2) + time0, n_middle = self._timed_deplete(n, rates, dt / 2, _i) res_middle = self.operator(n_middle, source_rate) # deplete across entire interval with BOS concentrations, # MOS reaction rates - time1, n_end = self._timed_deplete(n, res_middle.rates, dt) + time1, n_end = self._timed_deplete(n, res_middle.rates, dt, _i) return time0 + time1, [n_middle, n_end], [res_middle] @@ -172,26 +172,26 @@ class CF4Integrator(Integrator): """ # Step 1: deplete with matrix 1/2*A(y0) time1, n_eos1 = self._timed_deplete( - n_bos, bos_rates, dt, matrix_func=cf4_f1) + n_bos, bos_rates, dt, _i, matrix_func=cf4_f1) res1 = self.operator(n_eos1, source_rate) # Step 2: deplete with matrix 1/2*A(y1) time2, n_eos2 = self._timed_deplete( - n_bos, res1.rates, dt, matrix_func=cf4_f1) + n_bos, res1.rates, dt, _i, matrix_func=cf4_f1) res2 = self.operator(n_eos2, source_rate) # Step 3: deplete with matrix -1/2*A(y0)+A(y2) list_rates = list(zip(bos_rates, res2.rates)) time3, n_eos3 = self._timed_deplete( - n_eos1, list_rates, dt, matrix_func=cf4_f2) + n_eos1, list_rates, dt, _i, matrix_func=cf4_f2) res3 = self.operator(n_eos3, source_rate) # Step 4: deplete with two matrix exponentials list_rates = list(zip(bos_rates, res1.rates, res2.rates, res3.rates)) time4, n_inter = self._timed_deplete( - n_bos, list_rates, dt, matrix_func=cf4_f3) + n_bos, list_rates, dt, _i, matrix_func=cf4_f3) time5, n_eos5 = self._timed_deplete( - n_inter, list_rates, dt, matrix_func=cf4_f4) + n_inter, list_rates, dt, _i, matrix_func=cf4_f4) return (time1 + time2 + time3 + time4 + time5, [n_eos1, n_eos2, n_eos3, n_eos5], @@ -249,17 +249,17 @@ class CELIIntegrator(Integrator): simulation """ # deplete to end using BOS rates - proc_time, n_ce = self._timed_deplete(n_bos, rates, dt) + proc_time, n_ce = self._timed_deplete(n_bos, rates, dt, _i) res_ce = self.operator(n_ce, source_rate) # deplete using two matrix exponentials list_rates = list(zip(rates, res_ce.rates)) time_le1, n_inter = self._timed_deplete( - n_bos, list_rates, dt, matrix_func=celi_f1) + n_bos, list_rates, dt, _i, matrix_func=celi_f1) time_le2, n_end = self._timed_deplete( - n_inter, list_rates, dt, matrix_func=celi_f2) + n_inter, list_rates, dt, _i, matrix_func=celi_f2) return proc_time + time_le1 + time_le1, [n_ce, n_end], [res_ce] @@ -316,20 +316,20 @@ class EPCRK4Integrator(Integrator): """ # Step 1: deplete with matrix A(y0) / 2 - time1, n1 = self._timed_deplete(n, rates, dt, matrix_func=rk4_f1) + time1, n1 = self._timed_deplete(n, rates, dt, _i, matrix_func=rk4_f1) res1 = self.operator(n1, source_rate) # Step 2: deplete with matrix A(y1) / 2 - time2, n2 = self._timed_deplete(n, res1.rates, dt, matrix_func=rk4_f1) + time2, n2 = self._timed_deplete(n, res1.rates, dt, _i, matrix_func=rk4_f1) res2 = self.operator(n2, source_rate) # Step 3: deplete with matrix A(y2) - time3, n3 = self._timed_deplete(n, res2.rates, dt) + time3, n3 = self._timed_deplete(n, res2.rates, dt, _i) res3 = self.operator(n3, source_rate) # Step 4: deplete with matrix built from weighted rates list_rates = list(zip(rates, res1.rates, res2.rates, res3.rates)) - time4, n4 = self._timed_deplete(n, list_rates, dt, matrix_func=rk4_f4) + time4, n4 = self._timed_deplete(n, list_rates, dt, _i, matrix_func=rk4_f4) return (time1 + time2 + time3 + time4, [n1, n2, n3, n4], [res1, res2, res3]) @@ -414,9 +414,9 @@ class LEQIIntegrator(Integrator): self._prev_rates, bos_res.rates, repeat(prev_dt), repeat(dt))) time1, n_inter = self._timed_deplete( - n_bos, le_inputs, dt, matrix_func=leqi_f1) + n_bos, le_inputs, dt, i, matrix_func=leqi_f1) time2, n_eos0 = self._timed_deplete( - n_inter, le_inputs, dt, matrix_func=leqi_f2) + n_inter, le_inputs, dt, i, matrix_func=leqi_f2) res_inter = self.operator(n_eos0, source_rate) @@ -425,9 +425,9 @@ class LEQIIntegrator(Integrator): repeat(prev_dt), repeat(dt))) time3, n_inter = self._timed_deplete( - n_bos, qi_inputs, dt, matrix_func=leqi_f3) + n_bos, qi_inputs, dt, i, matrix_func=leqi_f3) time4, n_eos1 = self._timed_deplete( - n_inter, qi_inputs, dt, matrix_func=leqi_f4) + n_inter, qi_inputs, dt, i, matrix_func=leqi_f4) # store updated rates self._prev_rates = copy.deepcopy(bos_res.rates) @@ -478,7 +478,7 @@ class SICELIIntegrator(SIIntegrator): Eigenvalue and reaction rates from intermediate transport simulations """ - proc_time, n_eos = self._timed_deplete(n_bos, bos_rates, dt) + proc_time, n_eos = self._timed_deplete(n_bos, bos_rates, dt, _i) n_inter = copy.deepcopy(n_eos) # Begin iteration @@ -494,9 +494,9 @@ class SICELIIntegrator(SIIntegrator): list_rates = list(zip(bos_rates, res_bar.rates)) time1, n_inter = self._timed_deplete( - n_bos, list_rates, dt, matrix_func=celi_f1) + n_bos, list_rates, dt, _i, matrix_func=celi_f1) time2, n_inter = self._timed_deplete( - n_inter, list_rates, dt, matrix_func=celi_f2) + n_inter, list_rates, dt, _i, matrix_func=celi_f2) proc_time += time1 + time2 # end iteration @@ -560,9 +560,9 @@ class SILEQIIntegrator(SIIntegrator): inputs = list(zip(self._prev_rates, bos_rates, repeat(prev_dt), repeat(dt))) proc_time, n_inter = self._timed_deplete( - n_bos, inputs, dt, matrix_func=leqi_f1) + n_bos, inputs, dt, i, matrix_func=leqi_f1) time1, n_eos = self._timed_deplete( - n_inter, inputs, dt, matrix_func=leqi_f2) + n_inter, inputs, dt, i, matrix_func=leqi_f2) proc_time += time1 n_inter = copy.deepcopy(n_eos) @@ -580,9 +580,9 @@ class SILEQIIntegrator(SIIntegrator): inputs = list(zip(self._prev_rates, bos_rates, res_bar.rates, repeat(prev_dt), repeat(dt))) time1, n_inter = self._timed_deplete( - n_bos, inputs, dt, matrix_func=leqi_f3) + n_bos, inputs, dt, i, matrix_func=leqi_f3) time2, n_inter = self._timed_deplete( - n_inter, inputs, dt, matrix_func=leqi_f4) + n_inter, inputs, dt, i, matrix_func=leqi_f4) proc_time += time1 + time2 return proc_time, [n_eos, n_inter], [res_bar] diff --git a/openmc/deplete/pool.py b/openmc/deplete/pool.py index 27ecaa4dd..03b050af3 100644 --- a/openmc/deplete/pool.py +++ b/openmc/deplete/pool.py @@ -5,7 +5,7 @@ Provided to avoid some circular imports from itertools import repeat, starmap from multiprocessing import Pool -from scipy.sparse import bmat +from scipy.sparse import bmat, hstack, vstack, csc_matrix import numpy as np from openmc.mpi import comm @@ -40,8 +40,8 @@ def _distribute(items): return items[j:j + chunk_size] j += chunk_size -def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None, - *matrix_args): +def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None, + transfer_rates=None, external_source_rates=None, *matrix_args): """Deplete materials using given reaction rates for a specified time Parameters @@ -58,15 +58,21 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None, Reaction rates (from transport operator) dt : float Time in [s] to deplete for + current_timestep : int + Current timestep index maxtrix_func : callable, optional 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`` transfer_rates : openmc.deplete.TransferRates, Optional - Object to perform continuous reprocessing. + Transfer rates for continuous removal/feed. .. versionadded:: 0.14.0 + external_source_rates : openmc.deplete.ExternalSourceRates, Optional + External source rates for continuous removal/feed. + + .. versionadded:: 0.15.3 matrix_args: Any, optional Additional arguments passed to matrix_func @@ -93,15 +99,17 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None, matrices = map(matrix_func, repeat(chain), rates, fission_yields, *matrix_args) - if transfer_rates is not None: + if (transfer_rates is not None and + current_timestep in transfer_rates.external_timesteps): # Calculate transfer rate terms as diagonal matrices transfers = map(chain.form_rr_term, repeat(transfer_rates), - transfer_rates.local_mats) + repeat(current_timestep), transfer_rates.local_mats) + # Subtract transfer rate terms from Bateman matrices matrices = [matrix - transfer for (matrix, transfer) in zip(matrices, transfers)] - if len(transfer_rates.index_transfer) > 0: + if current_timestep in transfer_rates.index_transfer: # Gather all on comm.rank 0 matrices = comm.gather(matrices) n = comm.gather(n) @@ -112,10 +120,12 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None, n = [n_elm for n_mat in n for n_elm in n_mat] # Calculate transfer rate terms as diagonal matrices - transfer_pair = { - mat_pair: chain.form_rr_term(transfer_rates, mat_pair) - for mat_pair in transfer_rates.index_transfer - } + transfer_pair = {} + for mat_pair in transfer_rates.index_transfer[current_timestep]: + transfer_matrix = chain.form_rr_term(transfer_rates, + current_timestep, + mat_pair) + transfer_pair[mat_pair] = transfer_matrix # Combine all matrices together in a single matrix of matrices # to be solved in one go @@ -129,7 +139,7 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None, if row == col: # Fill the diagonals with the Bateman matrices cols.append(matrices[row]) - elif mat_pair in transfer_rates.index_transfer: + elif mat_pair in transfer_rates.index_transfer[current_timestep]: # Fill the off-diagonals with the transfer pair matrices cols.append(transfer_pair[mat_pair]) else: @@ -155,6 +165,25 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None, return n_result + if (external_source_rates is not None and + current_timestep in external_source_rates.external_timesteps): + # Calculate external source term vectors + sources = map(chain.form_ext_source_term, repeat(external_source_rates), + repeat(current_timestep), external_source_rates.local_mats) + + # stack vector column at the end of the matrix + matrices = [ + hstack([matrix, source]) + for matrix, source in zip(matrices, sources) + ] + + # Add a last row of zeroes to the matrices and append 1 to the last row + # of the nuclide vectors + for i, matrix in enumerate(matrices): + if not np.equal(*matrix.shape): + matrices[i] = vstack([matrix, csc_matrix([0]*matrix.shape[1])]) + n[i] = np.append(n[i], 1.0) + inputs = zip(matrices, n, repeat(dt)) if USE_MULTIPROCESSING: @@ -163,4 +192,10 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None, else: n_result = list(starmap(func, inputs)) + # Remove extra value at the end of the nuclide vectors + if (external_source_rates is not None and + current_timestep in external_source_rates.external_timesteps): + external_source_rates.reformat_nuclide_vectors(n) + external_source_rates.reformat_nuclide_vectors(n_result) + return n_result diff --git a/openmc/deplete/transfer_rates.py b/openmc/deplete/transfer_rates.py index 01c9b2e35..4c28c7d15 100644 --- a/openmc/deplete/transfer_rates.py +++ b/openmc/deplete/transfer_rates.py @@ -1,31 +1,31 @@ +from collections import defaultdict from numbers import Real import re +from typing import Sequence + +import numpy as np from openmc.checkvalue import check_type, check_value from openmc import Material -from openmc.data import ELEMENT_SYMBOL +from openmc.data import ELEMENT_SYMBOL, isotopes, AVOGADRO, atomic_mass +from .results import _SECONDS_PER_MINUTE, _SECONDS_PER_HOUR, \ + _SECONDS_PER_DAY, _SECONDS_PER_JULIAN_YEAR -class TransferRates: - """Class for defining continuous removals and feeds. - Molten Salt Reactors (MSRs) benefit from continuous reprocessing, - which removes fission products and feeds fresh fuel into the system. MSRs - inspired the development of this class. +class ExternalRates: + """External rates class for defining addition terms of depletion equation. - An instance of this class can be passed directly to an instance of one of - the :class:`openmc.deplete.Integrator` classes. - - .. versionadded:: 0.14.0 + .. versionadded:: 0.15.3 Parameters ---------- operator : openmc.TransportOperator Depletion operator - model : openmc.Model - OpenMC model containing materials and geometry. If using - :class:`openmc.deplete.CoupledOperator`, the model must also contain - a :class:`opnemc.Settings` object. + materials : openmc.Materials + OpenMC materials. + number_of_timesteps : int + Total number of depletion timesteps Attributes ---------- @@ -33,22 +33,25 @@ class TransferRates: All burnable material IDs. local_mats : list of str All burnable material IDs being managed by a single process - transfer_rates : dict of str to dict - Container of transfer rates, components (elements and/or nuclides) and - destination material - index_transfer : Set of pair of str - Pair of strings needed to build final matrix (destination_material, mat) + number_of_timesteps : int + Total number of depletion timesteps + external_rates : dict of str to dict + Container of timesteps, external rates, components (elements and/or + nuclides) and optionally destination material + external_timesteps : list of int + Container of all timesteps indeces with an external rate defined. """ - def __init__(self, operator, model): + def __init__(self, operator, materials, number_of_timesteps): - self.materials = model.materials + self.materials = materials self.burnable_mats = operator.burnable_mats self.local_mats = operator.local_mats + self.number_of_timesteps = number_of_timesteps - #initialize transfer rates container dict - self.transfer_rates = {mat: {} for mat in self.burnable_mats} - self.index_transfer = set() + # initialize transfer rates container dict + self.external_rates = {mat: defaultdict(list) for mat in self.burnable_mats} + self.external_timesteps = [] def _get_material_id(self, val): """Helper method for getting material id from Material obj or name. @@ -71,7 +74,7 @@ class TransferRates: check_value('Material ID', str(val), self.burnable_mats) else: check_value('Material name', val, - [mat.name for mat in self.materials if mat.depletable]) + [mat.name for mat in self.materials if mat.depletable]) val = [mat.id for mat in self.materials if mat.name == val][0] elif isinstance(val, int): @@ -79,7 +82,13 @@ class TransferRates: return str(val) - def get_transfer_rate(self, material, component): + def get_external_rate( + self, + material: str | int | Material, + component: str, + timestep: int, + destination_material: str | int | Material | None = None + ): """Return transfer rate for given material and element. Parameters @@ -88,62 +97,114 @@ class TransferRates: Depletable material component : str Element or nuclide to get transfer rate value + timestep : int + Current timestep index + destination_material : openmc.Material or str or int, Optional + Destination material to where nuclides get fed Returns ------- - transfer_rate : list of floats - Transfer rate values + external_rate : list of floats + External rate values """ material_id = self._get_material_id(material) check_type('component', component, str) - return [i[0] for i in self.transfer_rates[material_id][component]] - - def get_destination_material(self, material, component): - """Return destination material for given material and - component, if defined. - - Parameters - ---------- - material : openmc.Material or str or int - Depletable material - component : str - Element or nuclide that gets transferred to another material. - - Returns - ------- - destination_material_id : list of str - Depletable material ID to where the element or nuclide gets - transferred - - """ - material_id = self._get_material_id(material) - check_type('component', component, str) - if component in self.transfer_rates[material_id]: - return [i[1] for i in self.transfer_rates[material_id][component]] + if destination_material is not None: + dest_mat_id = self._get_material_id(destination_material) + return [i[1] for i in self.external_rates[material_id][component] + if timestep in i[0] and dest_mat_id == i[2]] else: - return [] + return [i[1] for i in self.external_rates[material_id][component] + if timestep in i[0]] - def get_components(self, material): - """Extract removing elements and/or nuclides for a given material + def get_components(self, material, timestep, destination_material=None): + """Extract removing elements and/or nuclides for a given material at a + given timestep Parameters ---------- material : openmc.Material or str or int Depletable material + timestep : int + Current timestep index + destination_material : openmc.Material or str or int, Optional + Destination material to where nuclides get fed Returns ------- - elements : list - List of elements and nuclides where transfer rates exist + components : list + List of elements or nuclides with external rates set at a given + timestep """ material_id = self._get_material_id(material) - if material_id in self.transfer_rates: - return self.transfer_rates[material_id].keys() + if destination_material is not None: + dest_mat_id = self._get_material_id(destination_material) + else: + dest_mat_id = None + + all_components = [] + if material_id in self.external_rates: + mat_components = self.external_rates[material_id] + + for component in mat_components: + if dest_mat_id: + # check for both timestep and destination material ids + if np.isin(timestep, [val[0] for val in mat_components[component]]) and \ + np.isin(dest_mat_id, [val[2] for val in mat_components[component]]): + all_components.append(component) + else: + # check only for timesteps + if np.isin(timestep, [val[0] for val in mat_components[component]]): + all_components.append(component) + return all_components + + +class TransferRates(ExternalRates): + """Class for defining continuous removals and feeds. + + Molten Salt Reactors (MSRs) benefit from continuous reprocessing, + which removes fission products and feeds fresh fuel into the system. MSRs + inspired the development of this class. + + An instance of this class can be passed directly to an instance of one of + the :class:`openmc.deplete.Integrator` classes. + + .. versionadded:: 0.14.0 + + Parameters + ---------- + operator : openmc.TransportOperator + Depletion operator + materials : openmc.Materials + OpenMC materials. + number_of_timesteps : int + Total number of depletion timesteps + + Attributes + ---------- + burnable_mats : list of str + All burnable material IDs. + local_mats : list of str + All burnable material IDs being managed by a single process + external_rates : dict of str to dict + Container of timesteps, transfer rates, components (elements and/or + nuclides) and destination material + external_timesteps : list of int + Container of all timesteps indeces with an external rate defined. + index_transfer : Set of pair of str + Pair of strings needed to build final matrix (destination_material, mat) + """ + + def __init__(self, operator, materials, number_of_timesteps): + super().__init__(operator, materials, number_of_timesteps) + self.index_transfer = defaultdict(list) + self.chain_nuclides = [nuc.name for nuc in operator.chain.nuclides] def set_transfer_rate(self, material, components, transfer_rate, - transfer_rate_units='1/s', destination_material=None): + transfer_rate_units='1/s', timesteps=None, + destination_material=None): """Set element and/or nuclide transfer rates in a depletable material. Parameters @@ -157,11 +218,14 @@ class TransferRates: transfer_rate : float Rate at which elements and/or nuclides are transferred. A positive or negative value corresponds to a removal or feed rate, respectively. - destination_material : openmc.Material or str or int, Optional - Destination material to where nuclides get fed. transfer_rate_units : {'1/s', '1/min', '1/h', '1/d', '1/a'} Units for values specified in the transfer_rate argument. 's' for seconds, 'min' for minutes, 'h' for hours, 'a' for Julian years. + timesteps : list of int, Optional + List of timestep indeces where to set transfer rates. + Default to None means the transfer rate is set for all timesteps. + destination_material : openmc.Material or str or int, Optional + Destination material to where nuclides get fed. """ material_id = self._get_material_id(material) @@ -183,19 +247,25 @@ class TransferRates: if transfer_rate_units in ('1/s', '1/sec'): unit_conv = 1 elif transfer_rate_units in ('1/min', '1/minute'): - unit_conv = 60 + unit_conv = _SECONDS_PER_MINUTE elif transfer_rate_units in ('1/h', '1/hr', '1/hour'): - unit_conv = 60*60 + unit_conv = _SECONDS_PER_HOUR elif transfer_rate_units in ('1/d', '1/day'): - unit_conv = 24*60*60 + unit_conv = _SECONDS_PER_DAY elif transfer_rate_units in ('1/a', '1/year'): - unit_conv = 365.25*24*60*60 + unit_conv = _SECONDS_PER_JULIAN_YEAR else: - raise ValueError('Invalid transfer rate unit ' - f'"{transfer_rate_units}"') + raise ValueError(f'Invalid transfer rate unit "{transfer_rate_units}"') + + if timesteps is not None: + for timestep in timesteps: + check_value('timestep', timestep, range(self.number_of_timesteps)) + timesteps = np.array(timesteps) + else: + timesteps = np.arange(self.number_of_timesteps) for component in components: - current_components = self.transfer_rates[material_id].keys() + current_components = self.external_rates[material_id].keys() split_component = re.split(r'\d+', component) element = split_component[0] if element not in ELEMENT_SYMBOL.values(): @@ -219,11 +289,144 @@ class TransferRates: f'where element {element} already has ' 'a transfer rate.') - if component in self.transfer_rates[material_id]: - self.transfer_rates[material_id][component].append( - (transfer_rate / unit_conv, destination_material_id)) - else: - self.transfer_rates[material_id][component] = [ - (transfer_rate / unit_conv, destination_material_id)] + self.external_rates[material_id][component].append( + (timesteps, transfer_rate/unit_conv, destination_material_id)) + if destination_material_id is not None: - self.index_transfer.add((destination_material_id, material_id)) + for timestep in timesteps: + self.index_transfer[timestep].append( + (destination_material_id, material_id)) + + self.external_timesteps = np.unique(np.concatenate( + [self.external_timesteps, timesteps])) + + +class ExternalSourceRates(ExternalRates): + """Class for defining external source rates. + + An instance of this class can be passed directly to an instance of one of + the :class:`openmc.deplete.Integrator` classes. + + .. versionadded:: 0.15.3 + + Parameters + ---------- + operator : openmc.TransportOperator + Depletion operator + materials : openmc.Materials + OpenMC materials. + number_of_timesteps : int + Total number of depletion timesteps + + Attributes + ---------- + burnable_mats : list of str + All burnable material IDs. + local_mats : list of str + All burnable material IDs being managed by a single process + external_timesteps : list of int + Container of all timesteps indeces with an external rate defined. + external_rates : dict of str to dict + Container of timesteps external source rates, and components + (elements and/or nuclides) + """ + + def reformat_nuclide_vectors(self, vectors): + """Remove last element of nuclide vector that was added for handling + external source rates by the depletion solver. + + Parameters + ---------- + vectors : list of array + List of nuclides vector to reformat + + """ + for mat_index, i in enumerate(self.local_mats): + if self.external_rates[i]: + vectors[mat_index] = vectors[mat_index][:-1] + + def set_external_source_rate( + self, + material: str | int | Material, + composition: dict[str, float], + rate: float, + rate_units: str = 'g/s', + timesteps: Sequence[int] | None = None + ): + """Set element and/or nuclide composition vector external source rates + to a depletable material. + + Parameters + ---------- + material : openmc.Material or str or int + Depletable material + composition : dict of str to float + External source rate composition vector, where key can be an element + or a nuclide and value the corresponding weight percent. + rate : float + External source rate in units of mass per time. A positive or + negative value corresponds to a feed or removal rate, respectively. + units : {'g/s', 'g/min', 'g/h', 'g/d', 'g/a'} + Units for values specified in the `rate` argument. 's' for seconds, + 'min' for minutes, 'h' for hours, 'a' for Julian years. + timesteps : list of int, optional + List of timestep indices where to set external source rates. Default + to None, which means the external source rate is set for all + timesteps. + + """ + + material_id = self._get_material_id(material) + check_type('rate', rate, Real) + check_type('composition', composition, dict, str) + + if rate_units in ('g/s', 'g/sec'): + unit_conv = 1 + elif rate_units in ('g/min', 'g/minute'): + unit_conv = _SECONDS_PER_MINUTE + elif rate_units in ('g/h', 'g/hr', 'g/hour'): + unit_conv = _SECONDS_PER_HOUR + elif rate_units in ('g/d', 'g/day'): + unit_conv = _SECONDS_PER_DAY + elif rate_units in ('g/a', 'g/year'): + unit_conv = _SECONDS_PER_JULIAN_YEAR + else: + raise ValueError(f'Invalid external source rate unit "{rate_units}"') + + if timesteps is not None: + for timestep in timesteps: + check_value('timestep', timestep, range(self.number_of_timesteps)) + timesteps = np.asarray(timesteps) + else: + timesteps = np.arange(self.number_of_timesteps) + + components = composition.keys() + percents = composition.values() + norm_percents = [float(i) / sum(percents) for i in percents] + + atoms_per_nuc = {} + for component, percent in zip(components, norm_percents): + split_component = re.split(r'\d+', component) + element = split_component[0] + if element not in ELEMENT_SYMBOL.values(): + raise ValueError(f'{component} is not a valid nuclide or element.') + + if len(split_component) == 1: + if not isotopes(component): + raise ValueError(f'Cannot add element {component} ' + 'as it is not naturally abundant. ' + 'Specify a nuclide vector instead.') + for nuc, frac in isotopes(component): + atoms_per_nuc[nuc] = (rate / atomic_mass(nuc) * AVOGADRO * + frac * percent / unit_conv) + + else: + atoms_per_nuc[component] = (rate / atomic_mass(component) * + AVOGADRO * percent / unit_conv) + + for nuc, val in atoms_per_nuc.items(): + self.external_rates[material_id][nuc].append((timesteps, val, None)) + + self.external_timesteps = np.unique(np.concatenate( + [self.external_timesteps, timesteps] + )) diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_ext_source.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_ext_source.h5 new file mode 100644 index 0000000000000000000000000000000000000000..aa09e1bcf4e352f94798f2cc266718b5f77fcec1 GIT binary patch literal 37328 zcmeHPeQX@X6`%9j;Rr4{Lr9&psWuQ1$AQKsrg1;czVt4sX(_QF5KGR-Iq@Yv;V=AE z7b?{Sr!j2-8;a@_TA89yr&5h6X`n#GRRv)oikwuI(n_TlRi#F)RF^6!RTI+Po%h@I zjlFldjW2P-yFYIB&AfT@=C|+nc4lw)9;^y)`q;u{3x)CV`9!Ypn}_&1OP}!pSpymE z5jyZt*+=DYjy&+_7LB$-s`rz9%_}QLsea4t8#fA(BVZpTHH-cQ@~E6Ds&#GS6EVN} z>rzLcYU74#V~&Uy#91Vm-I*6qykb1Wb-t)F2lfb`$k(I0E85l`YiSmKpH3sQMLOp3 zNM;Zn^az(9@C9YQAboMU-rpzo#2TX=UvG_)y0}_bEYJn$*GF4*1${)VE!y7M*wIdk z9$mxG0pH;NDdJlp@ksK`AIV^LVDCsj>$YA9P_}%@}I|@9|)*<06YrP^#NXjC*ULM-Ig*B@FB>BZF-Fou>-4EA)6|a0zVouB>7U=62pNKWrM~$32WBsB`KVQhq z=Mk5Q#d^Q4md4H|eU(z6RnibaQFetWJFDc`-QU=C%h8ulR66CQ63VZfmlVKRoD~75 zymVaQIr7C^caUL#@F2<>iK&{fTf zK38692!KzcG%@gHj28G#an?{KJM?`CLU3d2TLw!r6#xN=J zIu(d0YA%*nqKEc}7-yK|G4jJtrb|9951>~= za`#eMUr2g`%aZltBqu&8>ow97pGwv%BDu;^fa3ftaa~pMvhRysWv)}}I$D~v=FX;F z(KZchR%1p^>sDHZ>V1Re#69||)mZPgMeBCU+h=th(P@^_*>TEOQ*_=qUzsrTR~iJI z^3@>45xOHk%vHVu&(wTl^TQp_^2lM^jys;gjuX#nmdk;dXNV9kQz78QGeL2LzS-iL zdcC&!;f`my z949Qzj(`)-$|;V}D_cA}Mc0Kw8qem3JD$GvaQp?(4&+PAdcC|wR*O4tT4(-Zua}0JA)60$)kLL&Q5yZ=U&K(#!Uxa!7Lo0S74iV% zs{70m(i^9BG4#}RvV1yCh;lLvJ@q;tBzr@Ym!PNSsVdU*t(6DRDb7ckAf6#|~IXR8@yq(&;`kSmqo74LxdHS<_MthDo z$28#+nrYdalQ+1@JTjY262QHw=J?2s zkvC^SVvBUzaBMn0hG|^p^utV7RMM#W-K0Ne*%qs7#Cyukox2-j^-+k{^O&fRk9Z2J zUTf10=#A0X-e_}IU1KM8?T8s~9|_Uk(%H5`~;?`(PI`IGba zcvGK0r}LHh`O!>gyg9igawbXfJMFwgY^Sq$G!qYrphwcpif6e$9co^p9DIW6tssDQb9Nf=8C$Z0-|1=)-r9WvG zPy_D+1Zm9b%G84wQFD>9zE|I?X&}9lP-=S-YSzA(T5rC$8tZ4zEdK}7pP}Rahh+!! z0}s=gnR+By1N~4UzhNCSEg+LSvcw250*nA7zz8q`i~u9R2rvSS03&cg5OCT@YS_qS z=*1B2C(%DsGd0Y-okU<4QeMt~7u1Q-EEfDvE>7=iOZ z!0CR1{~T?`>33!5W8d%8eP|%L1^OuY{y;K5EHMI%03*N%FanGKBftnS0`mz0xBG|j zbKEzKI_wvCpD>@IfPG;E7y(9r5nu!u0Y)G_0#573bkBt)Xx#`s+&5N7j619wp^yG* zKb9E*Mt~7u1Q-EEfDvE><`)8)t{W%QtQW@|)``46oL}diy!1AAW z9~xP8G<3MEzGD5EK={nw`qzghdc$)y=D?$G74A5+r~>xZUU{bP-koIm{IT0#wdBa@zx7Uams;{)j};gEpetm_-=19X zn+NjCE&0(U+a9UB|3*tb;(vC{SNHv^^7z5Ty$@}E!t%q!tMFvdEaOqLto_qi<+uEy z9g(+x5qf;<$P+Jw4~DihOiaAD{>9J-#SdM3w&+0Ug`$@mR!?jX9eDUhzubScG4$?p z&s<%4S9SQCzg~RDsei`9+rOV*d(X+zji%i>_^R18ysWqP8Vj%0W#73o@!M)k?tkgs z;sbmB`vnX&Q1#^6-fLr4`(KK+$gkQTv*eSNqaVC6_JAc%tStEPeQ(_#cGO@vOiF}bYS7K+lB+x@3+)GUfdmcrKfD~bvvu9ei-j!m+v3h_w=jp{IT+P Ucl%fGx&FKI{NlW#oOXWy2Me7taXSOqlGuvgUC!_3S$G-Ub z>PcrL;#Tt8F+hPS-`|d^MQ`AIlLrRv;c*jeVe(3^wOIarkbukYSUJnf{M1r7f93Lp z55gWQ!ujSi3S71r8`|-%_xWe{^^fg=bB6cV>8RWNTqwCp`-5Y>-QuzXv9b^U*fLur zwLe^U)i(71y?BKl1-Jn_UruH@P{wQzbHF;c6>c%v3*Zj&=RdOU=w%hS`j!62H(v_c z59fRRdEWA&!Vu0seo0GkdYLnnufxF5`hx9XGmiyaL8w?vzkQP~oWD(Ay|!=bM|;mD zlTY<4>ArS3Vt@IQOj&ZjfpgqqU4doi^qr^FGp>wm^mo?h>M_q;tK+PwX{*|1a>{A? zWXTEXhIW5`ZSN}QgBuvGerM^9V~8N0el6erozFzLk2)rAoFHGnY}V?<#(Q+&x*n#u zN`E_LVz)JC)|5kx|7>sGX$`Os*Rb1@__2D{KS#SI+>dqI&iuEn;8^OO zt(QTB0wr0;_XFx>R|lfGn& z5uCsMY-cgY88J>MM8&(C#?Xa9#tDeQIn<|Z#b*n3A$UNupP1r&||h;P(U delta 693 zcmcbxnCZe|rVXd;7#TL7wv%TBQm^gj1F6Z~j{kw2ZYO>aXSOqlGuvgUC!@?{$G-Ub zZIjMO#I5AFV}JsCzP}w+i{8Nbjt>ml!{a8{!sIPqYq9+MAOV+uv2vD|`KhIF{_f=q zAA~(rg!4UT6u4|LHndY)@AJ>@>mS>D=M3+!(^0qku~2fA_6Ns$yNhK9Vr3uxu@$sP zYJa%ys;$ZYd+`cA3UC8Xy`0Q)pp4lb=73vnE8JqT7r-4<#(!kp(aS1u^%MP%Z@v_? zAI{JI^StFlg&~~3)nlHyR>xUO(^j?3eJ$VqozFzLkET!FI6=OC`mEK9jrZulbv;gT zmHu|h#BN^BtSN^W|Jk0o(;8qOu3@()@niL@e~xzT+>dqI&iuD6=2-ghbN*G^DWCW6 zt(QTBMM8Z_Ra#Xa9#tDI9hB<|Z#b*vCdpUNupP1r&||wBXb{ diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_removal.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_removal.h5 index 4b24bed52ad30a270df65704cecf8aef59078557..271af24103cda8f8c4886da551422a824dcde360 100644 GIT binary patch delta 720 zcmcbxnCZe|rVXd;7@0Powv%TBQm^gj1F6Z~j{kw2ZYO>aXSOqlGuvgUC!_3S$G-Ub z`blRb;#Tt8F+hO|-`|d^MQ`AI!v_ZK;c*jeVe-PSwOIarkbukYUOCIl{M1r7fAR8# z55gWQ!ue)13S71r8`}L}=kw3*>mS=a=M3+!(^0qky-;$M_6Ns$yBTE%Vr3uxvHfSB z)c$bYRok%t_u>_L6yOFdcsZHnKpC?=%mF*xR=CAvFMvD9ng7VTqnB0S>QDF|-+U=( zKb-IX=XuMA3PU*m?j1EDPz77M!#20J_n|UnY3QEOd`t6%^;ry)v>$QDbKiW@R zHu+S)l5Tw95&LVOWXh8J4V)i-%?pV7q2jzS|7oOSos~0-Rkiu-Zz|4;nzpKKCa0XH zPnMjJUa#TK<#D4A5yZl)U3P48Lfz{@F?$ zH(nv%Fo^RFe%OS1TEZ2iO+31{_^BzJzm#RxO@U1x?VT4* LUNupP1r)>paXSOqlGuvgUC!@?{$G-Ub zU6am8#I5AFV}JrDzP}w+i{8NbE)NXa!{a8{!sM-AYq9+MAOV+uy>gb9`KhIF{=Vf4 zAA~(rg!8>-6u4|LHnh`P@AJ>@>mS?u=M3+!(^0qku~2fA_6Ns$yDMb}Vr3uxu@$jM zYJa%ys;$}od+`cA3UC9?yqwH(pp4lb=74){E8JqT7r-4<#eZbo(aS1u_0#;1Z@v_? zAI{JE^StFlg&~~3>ynn>^fG5CUx$IA>jm4vW*!T;f^%Xq{q{||aQ=FM_1eCzAMF#C zPCnJIq<^@FkP;p+E|1{FE&dT|ZWwrV2Zz|5pnzpKKCa0WM zOqQIGUT^5m<#D4A5yXP4U3P48Lfz{@Jn| zH(?*JrcAgLsoqrpqj?rSKpQy)n=1-F!rp9Q#pUwCd0i zxc=&VFURQ>MsWVhITbwGu>nv%Fo;VGe%OS1TEZ2yO+31{_^BzJ-^nuTrog6;_FYRR LubQaD0t#XPxH8vc diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_transfer.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_transfer.h5 index 51173f778b945f5c666bd9a8ed9ace75f47eb0ec..25365f6f3a486e7cc6ffa3ce844d4ed3fa17a95c 100644 GIT binary patch delta 747 zcmcbxnCZe|rVXd;7@0Powv%TBQm^gj1F6Z~j{kw2ZYO>aXSOqlGuvgUC!_3S$G-Ub znn`CQ;#Tt8F+hPC-`|d^MQ`AIod*W(;c*jeVe)dXwOIarkbujtT{+9k{M1r7fA#W( z55gWQ!ubv}3S71r8``n1_xWe{^^fg=bB6cV>8RWNS}3_n`-5Y>-Q2PRv9b^U*#0+9 zYJa%ys%^^ud+`cA3UCAVznsi+pp4lb=72SBE8JqT7r-54$A4tq(aS1u^=ti)Z@v_? zAI?ws^StFlg&~}O?UI(@^fG5CUx$HV!V9*8%{&%x1(U>L`t6%^;ruNE>$QDbKiW$y zn|!KYIYqL_{>mqrvgCdP=lfst0-}DXI4{V58tGVP<^0dG+I;pm6=y|FTh%s`Q%=*? zFN}Km^|q7b6BRc^s6l{FfSF{d#K&8>4(#Mfvfj7Ezoh^ zJvnfKzTxvU{%g7!d~kQHYY!LjzWo}`UoVu~79}td9uV%6H%^eRSKe)LTyU`-TvY;- zPRfT!JG-#7kbTw*|JmB}?}#v$(6rkf@axu!eU5h5#G9{8-1^UU!k!BcORitDz1ih? zCPfMne=FouU#G=p=3n$G%*2jC9K&N)0GkkJIr z-*Q~&^`-5uaQ?N71A)5=Ea7~~OTuLjZomy>@MwKGX-dW?`+quQC{jF delta 747 zcmcbxnCZe|rVXd;7#TL7wv%TBQm^gj1F6Z~j{kw2ZYO>aXSOqlGuvgUC!@?{$G-Ub zZIjMO#I5AFV}JsCzP}w+i{8NbP7e&)!{a8{!sLx#Yq9+MAOV+uwsMx2`KhIF{@&#a zAA~(rg!A2J6u4|LHndY;@AJ>@>mS>D=M3+!(^0qku~2fA_6Ns$yGvyUVr3uxu@$sP zYJa%ys;%+=d+`cA3UC8Xy`0Q)pp4lb=75`SE8JqT7r-4<#(!kp(aS1u_5J?GH(v_c z59eq8dEWA&!Vu2ia!E^YdYLnnufxF5@`CMPGmiya!DX?Se)}d}IDfssdTrmm zC!gw9PMKS1f8~=*S#rOD^WCp`0Z~6xoR{Z6jdZNDa%QrsHlO`X#aTad^^4A~US^?dzuUT8Iwe*<{c4K~%u6}<9%?si`r-6R;OhNN3v`^f zP7a))Z}=*W|C(+FAKV>V+QS9BZ@-4~4+-VAMF~uV2ZYt+jT7YS6?a=47hJ3dS5?fU zlky?b&Mqb`WS{lIf41iQJ0i>_H0`zq{JOPbpQGJ#@#bq2xBjzTwCBRZlIz!OuXK5y zNs&UtpTN-%vElm=!GBgNIa7VsLb$D*YDrrIjNzdb&5)j;rnCOy0k{KVat==jWHf>E zR~;96eQCQZoPRLmK;W(dOE~}SMd7jsH{b>`c(uNqG$rGc{af|Pt0pQl1A`v`ms{(o diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_feed.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_feed.h5 index d53adc7d8413c904289119fcb922be662b396f96..d080b447093018f62b3eced413db81fa606d3338 100644 GIT binary patch delta 77 zcmcbxnCZe|rVXd;7@0Powv%TBQm^gj1F6Z~j{kw2ZYO>aXSOqlGuvgUC!_3U$G+8! Z%msdjCnvhfPrM+sSz*E&MiwC79spTR99sYY delta 77 zcmcbxnCZe|rVXd;7#TL7wv%TBQm^gj1F6Z~j{kw2ZYO>aXSOqlGuvgUC!@?}$G+8! Z%zKV6nw;n=KkaXSOqlGuvgUC!_3U$G+8! Z%!S%VCnvhfPrM+sSz*E&MiwC79spL192)=t delta 77 zcmcbxnCZe|rVXd;7#TL7wv%TBQm^gj1F6Z~j{kw2ZYO>aXSOqlGuvgUC!@?}$G+8! Z%u{q0O-^)`pLjuNv%-Wmj4VLDJpfd>91s8i diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_ext_source.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_ext_source.h5 new file mode 100644 index 0000000000000000000000000000000000000000..f3b1b171ebcab9d99fb7b73f5b8273ac72187971 GIT binary patch literal 37328 zcmeHPO>7%Q6rQz{)@=i%rIb)nv=vA-|4y5>0o2x|PDw8`f`WWN;56Abmhwk_N)Q~X z3Q*tx=`AM|2_dBiP!2sHaq$7998)Bu9>Xoia;TI`Rc80SFB?xhcD#vG(>@|;Ck?^OOLkx&h)IXJ}Eef}(5rxw)FpW_3E%SkS0V|3tXjecJX zw>PtXH%<+6+2UNA29urKSB(F8G$OhaG{=j?-QK8sRVSvs}8i zE8orjiuvHzO{z1P*ii}9q`mu!H!_+Z8dS{*O{1WP409Y}2Gc=%IDNnulue4|)eb#A zp)Te7y|K@SJr-9_YQf#Q0{LEVSPSSG-4Sneynk$z1syHn%YbjNf0_B#!aO3rIaa<+ zQ$#SL0j!boVMhRb%ko0uPwnhKSbr-Y9Oe@KnT9}><@@iR&p4`4P4KvMro4h6C}@%9 z)8*xelS(Zx4`8Rs3#_y71nGm^64w)ZsJ)%1OP?%1?}O^KU3V9EkwsH%7{vi)jROzL zewnvE5HRBac$DVr1FQs3z(=%4Eh!Gf_*$ye?d_xq^Ps~uu2aKJ#d6bJ|LyUpmi)iR z;regl)ud9v!W^|t=PwDBAMEw~n7bx_F_zz?DDpXKyV|SAT^Z^hAJAFKZub;VkydTH zRNH+c&h~tH@!YN7|2bugODT@85|<3XrCJgJTU@$r@Eq}Cv*Hr?WX6?h`2@S7^69V1 zeIMPpcT@AJbiQxqE0y-w>8A*0m(iZ!gs!tB`qKEYApkzj@xs8Dc~08NmQ&9a%O$zq zxgr3;`Ydr>HSr4Xi|q;W)b6pN0e5hG;G#F;B4_niWXoGQ4b@}Q!NMJ#wfgz)h}Ydi zx6iuAyj7Yi>9NJDJAB?GUIj4AR|N!Y@oJj=2*VK{HY;9%XJ)(!`4Nq0jcjnJCmPS7 z$Huel0h&m71`m<490E3;DfT0btrgGA>vhPFXgq6RgF`*hcm_Q-o)ui0NO%Shk+K{D zHlFpl#iMG)Gc(@}`4Nq0usYNejc3qf<5}_rk`SH&DJjb#VB^`s^TnfT#j~k{G|wE* z7xE(-&+v0{s3#iFpvT6usY7J9@T}NRq*fXM8_(L=k1(oMJX_}LLJQ9q@*^70_~+9N z-v@yn8_$+sB4mVT@DM4>Azy40eWkqVWuRY&Az5qebIbWy-A3W!qWZiuS=>YSZ`^*%}EpT28IWteDr_+>bXT^{+uk&fvH^XrWa%P#vft&65#ygBJ&cDPJ zASPug1Z;6-+Tba848JO+tC}|O!Hf&l@&R`F@9G7|+~DB7*am7q{s(!~cd*fT&)e1? z*Ke}=eYW>Y8uhb$et#?;3sfYO8)%Nl8mG?&M?rTWqO1P29KtsXQJMtbYw^eiT~1?Q z-kSG+O*~MqgbU=V$A>mX@z}P-3uLr!*lK)y)ws;*!)h}sQfht==*KKa^4 zr$688p}O)sCK~7wPqotSj%WqC)5~A>2CsDYkMq#6y#Mx*QlmrTBR$?HzEm)E_P_py`r7kH>v!Vi&!6plReyfmz?gWfAw?@G^50eGB^Mf9!=tr$5P}ZTZN;DA zkH7)->et`G>yHnAO6!7e^4ceAxZ299*NX@I4S1lw0X$~E2=VaM)hQly{>=Y901naD zA6xu7T>mN_OjbT=TPT6|0n$8Y*QxT8C%HRmj=SBr%gfb!ub)3V zCH@a)s!GF~zo!8h2Od^ZbLWR*35-LD_=Y@Yl|vQ}l$szQ2nYg#fFK|U2m*qDARq_` z0)oJnAYj`^y4c87X~hihCo#TSHCM3_<#{n4rMa&}UG|q%o+uF)1Ox#=KoAfF1OY)n z5D)|e0YN|z5Ck@WfbD)l^E!1F_;+Q76itcqH&~i7uS~C6P9~J(M xYPB>SfAiT-KW9Gd%3SL^aQ?f03z=)javzW1ZS7bqDJ@%u0M08)p4dOX{{cmRix2<+ literal 0 HcmV?d00001 diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_transfer.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_transfer.h5 index 469149f6742736085bebd22499c4eb2a45f519ed..b0ba997407f2471adcb32897dd08378bfe312814 100644 GIT binary patch delta 77 zcmcbxnCZe|rVXd;7@0Powv%TBQm^gj1F6Z~j{kw2ZYO>aXSOqlGuvgUC!_3U$G+8! Z%nw!-O-^)`pLjuNv%-Wmj4VLDJpg559Z>)P delta 77 zcmcbxnCZe|rVXd;7#TL7wv%TBQm^gj1F6Z~j{kw2ZYO>aXSOqlGuvgUC!@?}$G+8! Z%nfS`CnvhfPrM+sSz*E&MiwC79spV!9Ebn_ diff --git a/tests/regression_tests/deplete_with_transfer_rates/test.py b/tests/regression_tests/deplete_with_transfer_rates/test.py index 10d60866a..4ad009d22 100644 --- a/tests/regression_tests/deplete_with_transfer_rates/test.py +++ b/tests/regression_tests/deplete_with_transfer_rates/test.py @@ -1,8 +1,7 @@ -""" TransferRates depletion test suite """ +""" ExternalRates depletion test suite """ from pathlib import Path import shutil -import sys import numpy as np import pytest @@ -11,7 +10,7 @@ import openmc.deplete from openmc.deplete import CoupledOperator from tests.regression_tests import config, assert_reaction_rates_equal, \ - assert_atoms_equal, assert_same_mats + assert_atoms_equal @pytest.fixture @@ -44,9 +43,11 @@ def model(): settings.particles = 100 settings.inactive = 0 settings.batches = 10 + settings.seed = 1 return openmc.Model(geometry, materials, settings) + @pytest.mark.parametrize("rate, dest_mat, power, ref_result", [ (1e-5, None, 0.0, 'no_depletion_only_removal'), (-1e-5, None, 0.0, 'no_depletion_only_feed'), @@ -83,6 +84,40 @@ def test_transfer_rates(run_in_tmpdir, model, rate, dest_mat, power, ref_result) res_ref = openmc.deplete.Results(path_reference) res_test = openmc.deplete.Results(path_test) - assert_same_mats(res_ref, res_test) - assert_atoms_equal(res_ref, res_test, 1e-6) + assert_atoms_equal(res_ref, res_test) + assert_reaction_rates_equal(res_ref, res_test) + + +@pytest.mark.parametrize("rate, power, ref_result", [ + (1e-1, 0.0, 'no_depletion_with_ext_source'), + (1e-1, 174., 'depletion_with_ext_source'), +]) +def test_external_source_rates(run_in_tmpdir, model, rate, power, ref_result): + """Tests external_rates depletion class with external source rates""" + + chain_file = Path(__file__).parents[2] / 'chain_simple.xml' + + external_source_vector = {'U': 1} + + op = CoupledOperator(model, chain_file) + op.round_number = True + integrator = openmc.deplete.PredictorIntegrator( + op, [1], power, timestep_units='d') + integrator.add_external_source_rate('f', external_source_vector, rate) + integrator.integrate() + + # Get path to test and reference results + path_test = op.output_dir / 'depletion_results.h5' + path_reference = Path(__file__).with_name(f'ref_{ref_result}.h5') + + # If updating results, do so and return + if config['update']: + shutil.copyfile(str(path_test), str(path_reference)) + return + + # Load the reference/test results + res_ref = openmc.deplete.Results(path_reference) + res_test = openmc.deplete.Results(path_test) + + assert_atoms_equal(res_ref, res_test) assert_reaction_rates_equal(res_ref, res_test) diff --git a/tests/unit_tests/test_deplete_external_source_rates.py b/tests/unit_tests/test_deplete_external_source_rates.py new file mode 100644 index 000000000..a8cf9dde5 --- /dev/null +++ b/tests/unit_tests/test_deplete_external_source_rates.py @@ -0,0 +1,170 @@ +""" Tests for ExternalSourceRates class """ + +from pathlib import Path + +import pytest +import numpy as np +import re + +import openmc +from openmc.data import AVOGADRO, atomic_mass +from openmc.deplete import CoupledOperator +from openmc.deplete.transfer_rates import ExternalSourceRates +from openmc.deplete.abc import (_SECONDS_PER_MINUTE, _SECONDS_PER_HOUR, + _SECONDS_PER_DAY, _SECONDS_PER_JULIAN_YEAR) + +CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" + +@pytest.fixture +def model(): + openmc.reset_auto_ids() + f = openmc.Material(name="f") + f.add_element("U", 1, enrichment=4.25) + f.add_element("O", 2) + f.set_density("g/cm3", 10.4) + + w = openmc.Material(name="w") + w.add_element("O", 1) + w.add_element("H", 2) + w.set_density("g/cm3", 1.0) + w.depletable = True + + # material just to test multiple destination material + h = openmc.Material(name="h") + h.add_element("He", 1) + h.set_density("g/cm3", 1.78e-4) + + radii = [0.42, 0.45] + f.volume = np.pi * radii[0] ** 2 + w.volume = np.pi * (radii[1]**2 - radii[0]**2) + h.volume = 1 + materials = openmc.Materials([f, w, h]) + + surf_f = openmc.Sphere(r=radii[0]) + surf_w = openmc.Sphere(r=radii[1], boundary_type='vacuum') + surf_h = openmc.Sphere(x0=10, r=1, boundary_type='vacuum') + cell_f = openmc.Cell(fill=f, region=-surf_f) + cell_w = openmc.Cell(fill=w, region=+surf_f & -surf_w) + cell_h = openmc.Cell(fill=h, region=-surf_h) + geometry = openmc.Geometry([cell_f, cell_w, cell_h]) + + settings = openmc.Settings() + settings.particles = 1000 + settings.inactive = 10 + settings.batches = 50 + + return openmc.Model(geometry, materials, settings) + + +@pytest.mark.parametrize( +"case_name, external_source_vectors, external_source_rate, timesteps", [ + ('elements', [{'U': 0.9, 'Xe': 0.1}], 1, None), + ('nuclides', [{'I135': 0.1, 'Gd156': 0.3, 'Gd157': 0.6}], 1, None), + ('nuclides_elements', [{'I135': 0.01, 'Gd156': 0.1, 'Gd157': 0.01, 'U': 0.8, + 'Xe': 0.08}], 1, None), + ('elements_nuclides', [{'U': 0.78, 'Xe': 0.1, 'I135': 0.01, 'Gd156': 0.1, + 'Gd157': 0.01}], 1, None), + ('multiple_vectors', [{'U': 1.}, {'Xe': 1}], 1, None), + ('timesteps', [{'U': 0.9, 'Xe': 0.1}], 1, [1]), + ('rates_invalid_1', [{'Gb': 1.}], 1, None), + ('rates_invalid_2', [{'Pu': 1.}], 1, None) + ]) +def test_get_set(model, case_name, external_source_vectors, external_source_rate, + timesteps): + """Tests the get/set methods""" + + op = CoupledOperator(model, CHAIN_PATH) + number_of_timesteps = 2 + transfer = ExternalSourceRates(op, model.materials, number_of_timesteps) + + if timesteps is None: + timesteps = np.arange(number_of_timesteps) + + # Test by Openmc material, material name and material id + material= [m for m in model.materials if m.depletable][0] + + for material_input in [material, material.name, material.id]: + for external_source_vector in external_source_vectors: + if case_name == 'rates_invalid_1': + with pytest.raises(ValueError, match='Gb is not a valid ' + 'nuclide or element.'): + transfer.set_external_source_rate(material_input, + external_source_vector, + external_source_rate) + elif case_name == 'rates_invalid_2': + with pytest.raises(ValueError, match='Cannot add element Pu'): + transfer.set_external_source_rate(material_input, + external_source_vector, + external_source_rate) + else: + transfer.set_external_source_rate(material_input, + external_source_vector, + external_source_rate, + timesteps=timesteps) + for component, percent in external_source_vector.items(): + split_component = re.split(r'\d+', component) + if len(split_component) == 1: + for nuc, frac in openmc.data.isotopes(component): + val = external_source_rate * percent * frac * \ + AVOGADRO / atomic_mass(nuc) + assert transfer.get_external_rate( + material_input, nuc, timesteps)[0] == pytest.approx(val) + else: + val = external_source_rate * percent * AVOGADRO / atomic_mass(component) + assert transfer.get_external_rate( + material_input, component, timesteps)[0] == pytest.approx(val) + + assert np.all(transfer.external_timesteps == timesteps) + + +@pytest.mark.parametrize("units, unit_conv", [ + ('g/s', 1), + ('g/sec', 1), + ('g/min', _SECONDS_PER_MINUTE), + ('g/minute', _SECONDS_PER_MINUTE), + ('g/h', _SECONDS_PER_HOUR), + ('g/hr', _SECONDS_PER_HOUR), + ('g/hour', _SECONDS_PER_HOUR), + ('g/d', _SECONDS_PER_DAY), + ('g/day', _SECONDS_PER_DAY), + ('g/a', _SECONDS_PER_JULIAN_YEAR), + ('g/year', _SECONDS_PER_JULIAN_YEAR), + ]) +def test_units(units, unit_conv, model): + """ Units testing""" + # create external rate Xe + components = ['Xe135', 'U235'] + external_source_rate = 1.0 + number_of_timesteps = 2 + op = CoupledOperator(model, CHAIN_PATH) + transfer = ExternalSourceRates(op, model.materials, number_of_timesteps) + timesteps = np.arange(number_of_timesteps) + + for component in components: + rate = external_source_rate * unit_conv * atomic_mass(component) / AVOGADRO + transfer.set_external_source_rate('f', {component: 1}, rate, rate_units=units) + assert transfer.get_external_rate( + 'f', component, timesteps)[0] == pytest.approx(external_source_rate) + + +def test_external_source(run_in_tmpdir, model): + """Tests external source depletion class without neither reaction rates nor + decay but only external source rates""" + # create transfer rate for U + vector = {'U235': 1} + external_source = 10 # grams + op = CoupledOperator(model, CHAIN_PATH) + integrator = openmc.deplete.PredictorIntegrator( + op, [1, 1], 0.0, timestep_units = 'd') + integrator.add_external_source_rate('f', vector, external_source/(24*3600)) + integrator.integrate() + + # Get number of U238 atoms from results + results = openmc.deplete.Results('depletion_results.h5') + _, atoms = results.get_atoms(model.materials[0], "U235") + + # Ensure number of atoms equal external source + assert atoms[1] - atoms[0] == pytest.approx( + external_source * AVOGADRO / atomic_mass('U235')) + assert atoms[2] - atoms[1] == pytest.approx( + external_source * AVOGADRO / atomic_mass('U235')) diff --git a/tests/unit_tests/test_deplete_transfer_rates.py b/tests/unit_tests/test_deplete_transfer_rates.py index 140777cd6..a3228e9fb 100644 --- a/tests/unit_tests/test_deplete_transfer_rates.py +++ b/tests/unit_tests/test_deplete_transfer_rates.py @@ -16,6 +16,7 @@ CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml" @pytest.fixture def model(): + openmc.reset_auto_ids() f = openmc.Material(name="f") f.add_element("U", 1, percent_type="ao", enrichment=4.25) f.add_element("O", 2) @@ -54,32 +55,35 @@ def model(): return openmc.Model(geometry, materials, settings) -@pytest.mark.parametrize("case_name, transfer_rates", [ - ('elements', {'U': 0.01, 'Xe': 0.1}), - ('nuclides', {'I135': 0.01, 'Gd156': 0.1, 'Gd157': 0.01}), +@pytest.mark.parametrize("case_name, transfer_rates, timesteps", [ + ('elements', {'U': 0.01, 'Xe': 0.1}, None), + ('nuclides', {'I135': 0.01, 'Gd156': 0.1, 'Gd157': 0.01}, None), ('nuclides_elements', {'I135': 0.01, 'Gd156': 0.1, 'Gd157': 0.01, 'U': 0.01, - 'Xe': 0.1}), + 'Xe': 0.1}, None), ('elements_nuclides', {'U': 0.01, 'Xe': 0.1, 'I135': 0.01, 'Gd156': 0.1, - 'Gd157': 0.01}), + 'Gd157': 0.01}, None), ('multiple_transfer', {'U': 0.01, 'Xe': 0.1, 'I135': 0.01, 'Gd156': 0.1, - 'Gd157': 0.01}), - ('rates_invalid_1', {'Gd': 0.01, 'Gd157': 0.01, 'Gd156': 0.01}), - ('rates_invalid_2', {'Gd156': 0.01, 'Gd157': 0.01, 'Gd': 0.01}), - ('rates_invalid_3', {'Gb156': 0.01}), - ('rates_invalid_4', {'Gb': 0.01}) + 'Gd157': 0.01}, None), + ('timesteps', {'U': 0.01, 'Xe': 0.1}, [1]), + ('rates_invalid_1', {'Gd': 0.01, 'Gd157': 0.01, 'Gd156': 0.01}, None), + ('rates_invalid_2', {'Gd156': 0.01, 'Gd157': 0.01, 'Gd': 0.01}, None), + ('rates_invalid_3', {'Gb156': 0.01}, None), + ('rates_invalid_4', {'Gb': 0.01}, None) ]) -def test_get_set(model, case_name, transfer_rates): +def test_get_set(model, case_name, transfer_rates, timesteps): """Tests the get/set methods""" - - openmc.reset_auto_ids() op = CoupledOperator(model, CHAIN_PATH) - transfer = TransferRates(op, model) + number_of_timesteps = 2 + transfer = TransferRates(op, model.materials, number_of_timesteps) + + if timesteps is None: + timesteps = np.arange(number_of_timesteps) # Test by Openmc material, material name and material id material, dest_material, dest_material2 = [m for m in model.materials if m.depletable] for material_input in [material, material.name, material.id]: - for dest_material_input in [dest_material, dest_material.name, + for dest_material_input in [None, dest_material, dest_material.name, dest_material.id]: if case_name == 'rates_invalid_1': with pytest.raises(ValueError, match='Cannot add transfer ' @@ -117,14 +121,21 @@ def test_get_set(model, case_name, transfer_rates): for component, transfer_rate in transfer_rates.items(): transfer.set_transfer_rate(material_input, [component], transfer_rate, + timesteps=timesteps, destination_material=\ dest_material_input) - assert transfer.get_transfer_rate( - material_input, component)[0] == transfer_rate - assert transfer.get_destination_material( - material_input, component)[0] == str(dest_material.id) - assert transfer.get_components(material_input) == \ - transfer_rates.keys() + assert transfer.get_external_rate( + material_input, component, timesteps, + dest_material_input)[0] == transfer_rate + assert np.all(transfer.external_timesteps == timesteps) + + if timesteps is not None: + for timestep in timesteps: + assert transfer.get_components(material_input, timestep, + dest_material_input) == list(transfer_rates.keys()) + else: + assert transfer.get_components(material_input, timesteps, + dest_material_input) == list(transfer_rates.keys()) if case_name == 'multiple_transfer': for dest2_material_input in [dest_material2, dest_material2.name, @@ -135,10 +146,8 @@ def test_get_set(model, case_name, transfer_rates): destination_material=\ dest2_material_input) for id, dest_mat in zip([0,1],[dest_material,dest_material2]): - assert transfer.get_transfer_rate( - material_input, component)[id] == transfer_rate - assert transfer.get_destination_material( - material_input, component)[id] == str(dest_mat.id) + assert transfer.get_external_rate( + material_input, component, timesteps)[0] == transfer_rate @pytest.mark.parametrize("transfer_rate_units, unit_conv", [ ('1/s', 1), @@ -158,13 +167,15 @@ def test_units(transfer_rate_units, unit_conv, model): # create transfer rate Xe components = ['Xe', 'U235'] transfer_rate = 1e-5 + number_of_timesteps = 2 op = CoupledOperator(model, CHAIN_PATH) - transfer = TransferRates(op, model) + transfer = TransferRates(op, model.materials, number_of_timesteps) for component in components: transfer.set_transfer_rate('f', [component], transfer_rate * unit_conv, transfer_rate_units=transfer_rate_units) - assert transfer.get_transfer_rate('f', component)[0] == transfer_rate + for timestep in range(transfer.number_of_timesteps): + assert transfer.get_external_rate('f', component, timestep)[0] == transfer_rate def test_transfer(run_in_tmpdir, model):