diff --git a/docs/source/methods/depletion.rst b/docs/source/methods/depletion.rst index 1b131bed5..87a9976dd 100644 --- a/docs/source/methods/depletion.rst +++ b/docs/source/methods/depletion.rst @@ -257,3 +257,85 @@ choose one of two methods for estimating the heating rate, including: The method for normalization can be chosen through the ``normalization_mode`` argument to the :class:`openmc.deplete.CoupledOperator` class. + +-------------- +Transfer Rates +-------------- + +OpenMC allows continuous removal or feed of nuclides by adding an +extra transfer rate term to the depletion matrix. An application of this feature +is the chemical processing of Molten Salt Reactors (MSRs), where one can +model the removal of fission products or feeding fresh fuel into the system. + +A transfer rate as defined here is the rate at which nuclides are +continuously removed/fed from/to a material. + +.. note:: + + A transfer rate can be positive or negative, indicating removal or feed + respectively. + +Mathematically, it can be thought of as an additional term :math:`\mathbf{T}` +in the depletion equation that is proportional to the nuclide density, which can be written as: + +.. math:: + + \begin{aligned}\frac{dN_i(t)}{dt} = &\underbrace{\sum\limits_j f_{j\rightarrow i} + \int_0^\infty dE \; \sigma_j (E,t) \phi(E,t) N_j(t) - \int_0^\infty dE \; \sigma_i(E,t) + \phi(E,t) N_i(t)}_\textbf{R} \\ + &+ \underbrace{\sum_j \left [ \lambda_{j\rightarrow i} N_j(t) - \lambda_{i\rightarrow j} N_i(t) \right ]}_\textbf{D} \\ + &- \underbrace{t_i N_i(t)}_\textbf{T} \end{aligned} + +where the reaction term :math:`\mathbf{R}`, the decay term :math:`\mathbf{D}` +and the new transfer term :math:`\mathbf{T}` have been grouped together so that +:math:`\mathbf{A} = \mathbf{R}+\mathbf{D}-\mathbf{T}`. +The transfer rate coefficient :math:`t_i` defines the continuous transfer of the +nuclide :math:`i`, which behaves similar to radioactive decay. +:math:`t_i` can also be defined as the reciprocal of a cycle time +:math:`T_{cyc}`, intended as the time needed to process the whole inventory. + +Note that this formulation assumes homogeneous distribution of nuclide +:math:`i` throughout the material. + +A more rigorous description of removal rate and its implementation can be found +in the paper by `Hombourger +`_. + +The resulting burnup matrix can be solved with the same integration algorithms +that are used in the absence of the transfer term. + +.. note:: + + If no ``destination_material`` is specified, nuclides that are removed + or fed will not be tracked afterwards. + +Coupling materials +------------------ + +To keep track of removed nuclides or to feed nuclides from one depletable material +to another, the respective depletion equations have to be coupled. This can be +achieved by defining one block matrix, with diagonal blocks corresponding to +depletion matrices :math:`\mathbf{A_{ii}}`, where the index :math:`i` indicates +the depletable material id, and off-diagonal blocks corresponding to inter-material +coupling matrices :math:`\mathbf{T_{ij}}`, positioned so that that the indices :math:`i` and +:math:`j` indicate the nuclides receiving and losing materials, respectively. +The nuclide vectors are assembled together in one single vector and the resulting +system is solved with the same integration algorithms seen before. + +As an example, consider the case of two depletable materials and one +transfer defined from material 1 to material 2. The final system will look like: + +.. math:: + + \begin{aligned}\frac{d}{dt}\begin{pmatrix}\vec{N_1}\\ \vec{N_2}\end{pmatrix} &= + \begin{pmatrix}\mathbf{A_{11}} & \mathbf{0}\\ \mathbf{T_{21}} & \mathbf{A_{22 }} + \end{pmatrix} \begin{pmatrix}\vec{N_1}\\ \vec{N_2}\end{pmatrix} \end{aligned} + +where: + +:math:`\mathbf{A_{11}} = \mathbf{R_{11}}+\mathbf{D_{11}}-\mathbf{T_{21}}`, and + +:math:`\mathbf{A_{22}} = \mathbf{R_{22}}+\mathbf{D_{22}}`. + +Note that mass conservation is guaranteed by transferring the number +of atoms directly. diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index dd679a40a..dc9370ab5 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -54,7 +54,7 @@ provides the following transport operator classes: IndependentOperator The :class:`CoupledOperator` and :class:`IndependentOperator` classes must also -have some knowledge of how nuclides transmute and decay. This is handled by the +have some knowledge of how nuclides transmute and decay. This is handled by the :class:`Chain` class. Minimal Example @@ -195,14 +195,24 @@ total system energy. helpers.FissionYieldCutoffHelper helpers.FluxCollapseHelper -The :class:`openmc.deplete.IndependentOperator` uses inner classes subclassed +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. + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + transfer_rates.TransferRates + Intermediate Classes -------------------- Specific implementations of abstract base classes may utilize some of -the same methods and data structures. These methods and data are stored +the same methods and data structures. These methods and data are stored in intermediate classes. Methods common to tally-based implementation of :class:`FissionYieldHelper` diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index a16e98ca7..439e30305 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -262,7 +262,7 @@ Loading and Generating Microscopic Cross Sections ------------------------------------------------- As mentioned earlier, any transport code could be used to calculate one-group -microscopic cross sections. The :mod:`openmc.deplete` module provides the +microscopic cross sections. The :mod:`openmc.deplete` module provides the :class:`~openmc.deplete.MicroXS` class, which contains methods to read in pre-calculated cross sections from a ``.csv`` file or from data arrays:: @@ -355,9 +355,9 @@ Multiple Materials A transport-independent depletion simulation using ``source-rate`` normalization will calculate reaction rates for each material independently. This can be -useful for running many different cases of a particular scenario. A +useful for running many different cases of a particular scenario. A transport-independent depletion simulation using ``fission-q`` normalization -will sum the fission energy values across all materials into :math:`Q_i` in +will sum the fission energy values across all materials into :math:`Q_i` in Equation :math:numref:`fission-q`, and Equation :math:numref:`fission-q` provides the flux we use to calculate the reaction rates in each material. This can be useful for running a scenario with multiple depletable materials @@ -370,3 +370,67 @@ The values of the one-group microscopic cross sections passed to :class:`openmc.deplete.IndependentOperator` are fixed for the entire depletion simulation. This implicit assumption may produce inaccurate results for certain scenarios. + +Transfer Rates +============== + +Transfer rates define removal or feed of nuclides to or from one or more +depletable materials. This can be useful to model continuous fuel reprocessing, +online fission products separation, etc. + +Transfer rates are defined by calling the +:meth:`~openmc.deplete.abc.Integrator.add_transfer_rate()` method directly from +one of the Integrator classes:: + + ... + integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power) + integrator.add_transfer_rate(...) + +Defining transfer rates +----------------------- + +The :meth:`~openmc.deplete.abc.Integrator.add_transfer_rate()` method requires a +:class:`~openmc.Material` instance (alternatively, a material id or +the name) as the depletable material from which nuclides are processed, +a list of elements that share the same transfer rate, and a transfer rate itself. + +.. caution:: + + Make sure you set the transfer rate value with the right sign. + A positive transfer rate assumes removal, while a negative one assumes feed. + +The ``transfer_rate_units`` argument specifies the units for the transfer rate. +The default is `1/s`, but '1/min', '1/h', '1/d' and '1/a' are also valid +options. + +For example, to define continuous removal of xenon from one material with a +removal rate value of 0.1 s\ :sup:`-1` (or a cycle time of 10 s), you'd use:: + + mat1 = openmc.Material(material_id=1, name='fuel') + + ... + + integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power) + # by openmc.Material object + integrator.add_transfer_rate(mat1, ['Xe'], 0.1) + # or by material id + integrator.add_transfer_rate(1, ['Xe'], 0.1) + # or by material name + integrator.add_transfer_rate('fuel', ['Xe'], 0.1) + +Note that in this case the xenon isotopes that are removed will not be tracked. + +Defining a destination material +------------------------------- + +To transfer elements from one depletable material to another, the +``destination_material`` parameter needs to be passed to the +:meth:`~openmc.deplete.abc.Integrator.add_transfer_rate()` method. For example, +to transfer xenon from one material to another, you'd use:: + + ... + mat2 = openmc.Material(name='storage') + + ... + + integrator.add_transfer_rate(mat1, ['Xe'], 0.1, destination_material=mat2) diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 329ee8b52..8a9509e90 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -16,6 +16,7 @@ from .atom_number import * from .stepresult import * from .results import * from .integrators import * +from .transfer_rates import * from . import abc from . import cram from . import helpers diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 5b405b1d5..e871eccfe 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -24,6 +24,7 @@ from .stepresult import StepResult from .chain import Chain from .results import Results from .pool import deplete +from .transfer_rates import TransferRates __all__ = [ @@ -547,7 +548,6 @@ class Integrator(ABC): :attr:`solver`. .. versionadded:: 0.12 - Attributes ---------- operator : openmc.deplete.abc.TransportOperator @@ -574,7 +574,10 @@ class Integrator(ABC): * ``n1`` is a :class:`numpy.ndarray` of compositions at the next time step. Expected to be of the same shape as ``n0`` - .. versionadded:: 0.12 + transfer_rates : openmc.deplete.TransferRates + Instance of TransferRates class to perform continuous transfer during depletion + + .. versionadded:: 0.13.4 """ @@ -654,6 +657,8 @@ class Integrator(ABC): self.timesteps = asarray(seconds) self.source_rates = asarray(source_rates) + self.transfer_rates = None + if isinstance(solver, str): # Delay importing of cram module, which requires this file if solver == "cram48": @@ -704,7 +709,8 @@ class Integrator(ABC): def _timed_deplete(self, concs, rates, dt, matrix_func=None): start = time.time() results = deplete( - self._solver, self.chain, concs, rates, dt, matrix_func) + self._solver, self.chain, concs, rates, dt, matrix_func, + self.transfer_rates) return time.time() - start, results @abstractmethod @@ -835,6 +841,31 @@ class Integrator(ABC): self.operator.finalize() + def add_transfer_rate(self, material, elements, transfer_rate, + transfer_rate_units='1/s', destination_material=None): + """Add transfer rates to depletable material. + + Parameters + ---------- + material : openmc.Material or str or int + Depletable material + elements : list of str + List of strings of elements that share transfer rate + 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. + + """ + if self.transfer_rates is None: + self.transfer_rates = TransferRates(self.operator, self.operator.model) + + self.transfer_rates.set_transfer_rate(material, elements, transfer_rate, + transfer_rate_units, destination_material) @add_params class SIIntegrator(Integrator): diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 372ed3510..5e3157f90 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -698,6 +698,64 @@ class Chain: dict.update(matrix_dok, matrix) return matrix_dok.tocsr() + def form_rr_term(self, transfer_rates, materials): + """Function to form the transfer rate term matrices. + + .. versionadded:: 0.13.4 + + Parameters + ---------- + transfer_rates : openmc.deplete.TransferRates + Instance of openmc.deplete.TransferRates + materials : string or two-tuple of strings + Two cases are possible: + + 1) Material ID as string: + Nuclide transfer only. In this case the transfer rate terms will be + subtracted from the respective depletion matrix + + 2) Two-tuple of material IDs as strings: + Nuclide transfer from one material into another. + The pair is assumed to be + ``(destination_material, source_material)``, where + ``destination_material`` and ``source_material`` are the nuclide + receiving and losing materials, respectively. + The transfer rate terms get placed in the final matrix with indexing + position corresponding to the ID of the materials set. + + Returns + ------- + scipy.sparse.csr_matrix + Sparse matrix representing transfer term. + + """ + matrix = defaultdict(float) + + for i, nuclide in enumerate(self.nuclides): + element = re.split(r'\d+', nuclide.name)[0] + # Build transfer terms matrices + if isinstance(materials, str): + material = materials + if element in transfer_rates.get_elements(material): + matrix[i, i] = transfer_rates.get_transfer_rate(material, element) + else: + matrix[i, i] = 0.0 + #Build transfer terms matrices + elif isinstance(materials, tuple): + destination_material, material = materials + if transfer_rates.get_destination_material(material, element) == destination_material: + matrix[i, i] = transfer_rates.get_transfer_rate(material, element) + else: + warn(f'Material {destination_material} is not defined ' + f'as a destination material for Material {material}. ' + 'Setting transfer rate to 0.0') + matrix[i, i] = 0.0 + #Nothing else is allowed + n = len(self) + matrix_dok = sp.dok_matrix((n, n)) + dict.update(matrix_dok, matrix) + return matrix_dok.tocsr() + def get_branch_ratios(self, reaction="(n,gamma)"): """Return a dictionary with reaction branching ratios diff --git a/openmc/deplete/pool.py b/openmc/deplete/pool.py index aba60b558..fc006a884 100644 --- a/openmc/deplete/pool.py +++ b/openmc/deplete/pool.py @@ -4,6 +4,8 @@ Provided to avoid some circular imports """ from itertools import repeat, starmap from multiprocessing import Pool +from scipy.sparse import bmat +import numpy as np # Configurable switch that enables / disables the use of @@ -15,7 +17,8 @@ USE_MULTIPROCESSING = True NUM_PROCESSES = None -def deplete(func, chain, x, rates, dt, matrix_func=None, *matrix_args): +def deplete(func, chain, x, rates, dt, matrix_func=None, transfer_rates=None, + *matrix_args): """Deplete materials using given reaction rates for a specified time Parameters @@ -37,7 +40,11 @@ def deplete(func, chain, x, rates, dt, matrix_func=None, *matrix_args): ``fission_yields = {parent: {product: yield_frac}}`` Expected to return the depletion matrix required by ``func`` - matrix_args : Any, optional + transfer_rates : openmc.deplete.TransferRates, Optional + Object to perform continuous reprocessing. + + .. versionadded:: 0.13.4 + matrix_args: Any, optional Additional arguments passed to matrix_func Returns @@ -62,6 +69,51 @@ def deplete(func, chain, x, rates, dt, matrix_func=None, *matrix_args): matrices = map(matrix_func, repeat(chain), rates, fission_yields, *matrix_args) + if transfer_rates is not None: + # Calculate transfer rate terms as diagonal matrices + transfers = map(chain.form_rr_term, repeat(transfer_rates), + transfer_rates.burnable_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: + # 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 + } + + # Combine all matrices together in a single matrix of matrices + # to be solved in one go + n_rows = n_cols = len(transfer_rates.burnable_mats) + rows = [] + for row in range(n_rows): + cols = [] + for col in range(n_cols): + mat_pair = (transfer_rates.burnable_mats[row], + transfer_rates.burnable_mats[col]) + if row == col: + # Fill the diagonals with the Bateman matrices + cols.append(matrices[row]) + elif mat_pair in transfer_rates.index_transfer: + # Fill the off-diagonals with the transfer pair matrices + cols.append(transfer_pair[mat_pair]) + else: + cols.append(None) + + rows.append(cols) + matrix = bmat(rows) + + # Concatenate vectors of nuclides in one + x_multi = np.concatenate([xx for xx in x]) + x_result = func(matrix, x_multi, dt) + + # Split back the nuclide vector result into the original form + x_result = np.split(x_result, np.cumsum([len(i) for i in x])[:-1]) + + return x_result + inputs = zip(matrices, x, repeat(dt)) if USE_MULTIPROCESSING: diff --git a/openmc/deplete/transfer_rates.py b/openmc/deplete/transfer_rates.py new file mode 100644 index 000000000..c6f0cd94d --- /dev/null +++ b/openmc/deplete/transfer_rates.py @@ -0,0 +1,188 @@ +from numbers import Real + +from openmc.checkvalue import check_type, check_value +from openmc import Material +from openmc.data import ELEMENT_SYMBOL + + +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. + + An instance of this class can be passed directly to an instance of one of + the :class:`openmc.deplete.Integrator` classes. + + .. versionadded:: 0.13.4 + + 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. + + Attributes + ---------- + burnable_mats : list of str + All burnable material IDs. + transfer_rates : dict of str to dict + Container of transfer rates, elements and destination material + index_transfer : Set of pair of str + Pair of strings needed to build final matrix (destination_material, mat) + """ + + def __init__(self, operator, model): + + self.materials = model.materials + self.burnable_mats = operator.burnable_mats + + #initialize transfer rates container dict + self.transfer_rates = {mat: {} for mat in self.burnable_mats} + self.index_transfer = set() + + def _get_material_id(self, val): + """Helper method for getting material id from Material obj or name. + + Parameters + ---------- + val : openmc.Material or str or int representing material name/id + + Returns + ------- + material_id : str + + """ + if isinstance(val, Material): + check_value('Depeletable Material', str(val.id), self.burnable_mats) + val = val.id + + elif isinstance(val, str): + if val.isnumeric(): + 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]) + val = [mat.id for mat in self.materials if mat.name == val][0] + + elif isinstance(val, int): + check_value('Material ID', str(val), self.burnable_mats) + + return str(val) + + def get_transfer_rate(self, material, element): + """Return transfer rate for given material and element. + + Parameters + ---------- + material : openmc.Material or str or int + Depletable material + element : str + Element to get transfer rate value + + Returns + ------- + transfer_rate : float + Transfer rate value + + """ + material_id = self._get_material_id(material) + check_value('element', element, ELEMENT_SYMBOL.values()) + return self.transfer_rates[material_id][element][0] + + def get_destination_material(self, material, element): + """Return destination (or transfer) material for given material and + element, if defined. + + Parameters + ---------- + material : openmc.Material or str or int + Depletable material + element : str + Element that gets transferred to another material. + + Returns + ------- + destination_material_id : str + Depletable material ID to where the element gets transferred + + """ + material_id = self._get_material_id(material) + check_value('element', element, ELEMENT_SYMBOL.values()) + if element in self.transfer_rates[material_id]: + return self.transfer_rates[material_id][element][1] + + def get_elements(self, material): + """Extract removing elements for a given material + + Parameters + ---------- + material : openmc.Material or str or int + Depletable material + + Returns + ------- + elements : list + List of elements where transfer rates exist + + """ + material_id = self._get_material_id(material) + if material_id in self.transfer_rates: + return self.transfer_rates[material_id].keys() + + def set_transfer_rate(self, material, elements, transfer_rate, transfer_rate_units='1/s', + destination_material=None): + """Set element transfer rates in a depletable material. + + Parameters + ---------- + material : openmc.Material or str or int + Depletable material + elements : list of str + List of strings of elements that share transfer rate + 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. + + """ + material_id = self._get_material_id(material) + check_type('transfer_rate', transfer_rate, Real) + + if destination_material is not None: + destination_material_id = self._get_material_id(destination_material) + if len(self.burnable_mats) > 1: + check_value('destination_material', str(destination_material_id), + self.burnable_mats) + else: + raise ValueError(f'Transfer to material {destination_material_id} '\ + 'is set, but there is only one depletable material') + else: + destination_material_id = None + + if transfer_rate_units in ('1/s', '1/sec'): + unit_conv = 1 + elif transfer_rate_units in ('1/min', '1/minute'): + unit_conv = 60 + elif transfer_rate_units in ('1/h', '1/hr', '1/hour'): + unit_conv = 60*60 + elif transfer_rate_units in ('1/d', '1/day'): + unit_conv = 24*60*60 + elif transfer_rate_units in ('1/a', '1/year'): + unit_conv = 365.25*24*60*60 + else: + raise ValueError("Invalid transfer rate unit '{}'".format(transfer_rate_units)) + + for element in elements: + check_value('element', element, ELEMENT_SYMBOL.values()) + self.transfer_rates[material_id][element] = transfer_rate / unit_conv, destination_material_id + if destination_material_id is not None: + self.index_transfer.add((destination_material_id, material_id)) diff --git a/tests/regression_tests/deplete_with_transfer_rates/__init__.py b/tests/regression_tests/deplete_with_transfer_rates/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_feed.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_feed.h5 new file mode 100644 index 000000000..cb8d49bfc Binary files /dev/null and b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_feed.h5 differ 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 new file mode 100644 index 000000000..b20bbc645 Binary files /dev/null and b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_removal.h5 differ 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 new file mode 100644 index 000000000..d428cdd7c Binary files /dev/null and b/tests/regression_tests/deplete_with_transfer_rates/ref_depletion_with_transfer.h5 differ 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 new file mode 100644 index 000000000..95b7236f0 Binary files /dev/null and b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_feed.h5 differ diff --git a/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_removal.h5 b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_removal.h5 new file mode 100644 index 000000000..81c97c755 Binary files /dev/null and b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_only_removal.h5 differ 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 new file mode 100644 index 000000000..407c2114e Binary files /dev/null and b/tests/regression_tests/deplete_with_transfer_rates/ref_no_depletion_with_transfer.h5 differ diff --git a/tests/regression_tests/deplete_with_transfer_rates/test.py b/tests/regression_tests/deplete_with_transfer_rates/test.py new file mode 100644 index 000000000..7e20a7607 --- /dev/null +++ b/tests/regression_tests/deplete_with_transfer_rates/test.py @@ -0,0 +1,80 @@ +""" TransferRates depletion test suite """ + +from pathlib import Path + +import numpy as np +import pytest +import openmc +import openmc.deplete +from openmc.deplete import CoupledOperator + +@pytest.fixture +def model(): + f = openmc.Material(name="f") + f.add_element("U", 1, percent_type="ao", enrichment=4.25) + f.add_element("O", 2) + f.set_density("g/cc", 10.4) + + w = openmc.Material(name="w") + w.add_element("O", 1) + w.add_element("H", 2) + w.set_density("g/cc", 1.0) + w.depletable = True + + radii = [0.42, 0.45] + f.volume = np.pi * radii[0] ** 2 + w.volume = np.pi * (radii[1]**2 - radii[0]**2) + + materials = openmc.Materials([f, w]) + + surf_f = openmc.Sphere(r=radii[0]) + surf_w = openmc.Sphere(r=radii[1], boundary_type='reflective') + cell_f = openmc.Cell(fill=f, region=-surf_f) + cell_w = openmc.Cell(fill=w, region=+surf_f & -surf_w) + geometry = openmc.Geometry([cell_f, cell_w]) + + settings = openmc.Settings() + settings.particles = 100 + settings.inactive = 0 + settings.batches = 10 + + 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'), + (1e-5, None, 174.0, 'depletion_with_removal'), + (-1e-5, None, 174.0, 'depletion_with_feed'), + (-1e-5, 'w', 0.0, 'no_depletion_with_transfer'), + (1e-5, 'w', 174.0, 'depletion_with_transfer'), + ]) +def test_transfer_rates(run_in_tmpdir, model, rate, dest_mat, power, ref_result): + """Tests transfer_rates depletion class with transfer rates""" + + chain_file = Path(__file__).parents[2] / 'chain_simple.xml' + + transfer_elements = ['Xe'] + + op = CoupledOperator(model, chain_file) + integrator = openmc.deplete.PredictorIntegrator( + op, [1], power, timestep_units = 'd') + integrator.add_transfer_rate('f', transfer_elements, rate, + destination_material=dest_mat) + 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') + + # Load the reference/test results + res_ref = openmc.deplete.Results(path_reference) + res_test = openmc.deplete.Results(path_test) + + for mat in res_test[0].rates[0].index_mat: + for nuc in res_test[0].rates[0].index_nuc: + for rx in res_test[0].rates[0].index_rx: + y_test = res_test.get_reaction_rate(mat, nuc, rx)[1] / \ + res_test.get_atoms(mat, nuc)[1] + y_ref = res_ref.get_reaction_rate(mat, nuc, rx)[1] / \ + res_ref.get_atoms(mat, nuc)[1] + assert y_test == pytest.approx(y_ref, abs=1e-5) diff --git a/tests/unit_tests/test_deplete_transfer_rates.py b/tests/unit_tests/test_deplete_transfer_rates.py new file mode 100644 index 000000000..1572755c2 --- /dev/null +++ b/tests/unit_tests/test_deplete_transfer_rates.py @@ -0,0 +1,118 @@ +""" Tests for TransferRates class """ + +from pathlib import Path +from math import exp + +import pytest +import numpy as np + +import openmc +from openmc.deplete import CoupledOperator +from openmc.deplete.transfer_rates import TransferRates +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(): + f = openmc.Material(name="f") + f.add_element("U", 1, percent_type="ao", enrichment=4.25) + f.add_element("O", 2) + f.set_density("g/cc", 10.4) + + w = openmc.Material(name="w") + w.add_element("O", 1) + w.add_element("H", 2) + w.set_density("g/cc", 1.0) + w.depletable = True + + radii = [0.42, 0.45] + f.volume = np.pi * radii[0] ** 2 + w.volume = np.pi * (radii[1]**2 - radii[0]**2) + materials = openmc.Materials([f, w]) + + surf_f = openmc.Sphere(r=radii[0]) + surf_w = openmc.Sphere(r=radii[1], boundary_type='vacuum') + cell_f = openmc.Cell(fill=f, region=-surf_f) + cell_w = openmc.Cell(fill=w, region=+surf_f & -surf_w) + geometry = openmc.Geometry([cell_f, cell_w]) + + settings = openmc.Settings() + settings.particles = 1000 + settings.inactive = 10 + settings.batches = 50 + + return openmc.Model(geometry, materials, settings) + + +def test_get_set(model): + """Tests the get/set methods""" + + # create transfer rates for U and Xe + transfer_rates = {'U': 0.01, 'Xe': 0.1} + op = CoupledOperator(model, CHAIN_PATH) + transfer = TransferRates(op, model) + + # Test by Openmc material, material name and material id + material, dest_material = [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, + dest_material.id]: + for element, transfer_rate in transfer_rates.items(): + transfer.set_transfer_rate(material_input, [element], transfer_rate, + destination_material=dest_material_input) + assert transfer.get_transfer_rate(material_input, + element) == transfer_rate + assert transfer.get_destination_material(material_input, + element) == str(dest_material.id) + assert transfer.get_elements(material_input) == transfer_rates.keys() + + +@pytest.mark.parametrize("transfer_rate_units, unit_conv", [ + ('1/s', 1), + ('1/sec', 1), + ('1/min', _SECONDS_PER_MINUTE), + ('1/minute', _SECONDS_PER_MINUTE), + ('1/h', _SECONDS_PER_HOUR), + ('1/hr', _SECONDS_PER_HOUR), + ('1/hour', _SECONDS_PER_HOUR), + ('1/d', _SECONDS_PER_DAY), + ('1/day', _SECONDS_PER_DAY), + ('1/a', _SECONDS_PER_JULIAN_YEAR), + ('1/year', _SECONDS_PER_JULIAN_YEAR), + ]) +def test_units(transfer_rate_units, unit_conv, model): + """ Units testing""" + # create transfer rate Xe + element = 'Xe' + transfer_rate = 1e-5 + op = CoupledOperator(model, CHAIN_PATH) + transfer = TransferRates(op, model) + + transfer.set_transfer_rate('f', [element], transfer_rate * unit_conv, + transfer_rate_units=transfer_rate_units) + assert transfer.get_transfer_rate('f', element) == transfer_rate + + +def test_transfer(run_in_tmpdir, model): + """Tests transfer depletion class without neither reaction rates nor decay + but only transfer rates""" + + # create transfer rate for U + element = ['U'] + transfer_rate = 1e-5 + op = CoupledOperator(model, CHAIN_PATH) + integrator = openmc.deplete.PredictorIntegrator( + op, [1,1], 0.0, timestep_units = 'd') + integrator.add_transfer_rate('f', element, transfer_rate) + integrator.integrate() + + # Get number of U238 atoms from results + results = openmc.deplete.Results('depletion_results.h5') + _, atoms = results.get_atoms(model.materials[0], "U238") + + # Ensure number of atoms equal transfer decay + assert atoms[1] == pytest.approx(atoms[0]*exp(-transfer_rate*3600*24)) + assert atoms[2] == pytest.approx(atoms[1]*exp(-transfer_rate*3600*24))