TransferRates (old msr continuous capabilities)

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
+<https://doi.org/10.1016/j.anucene.2020.107504>`_.
+
+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.

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`

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)

openmc/deplete/__init__.py

@@ -16,6 +16,7 @@
 from .stepresult import *
 from .results import *
 from .integrators import *
+from .transfer_rates import *
 from . import abc
 from . import cram
 from . import helpers

openmc/deplete/abc.py

@@ -24,6 +24,7 @@
 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 @@ def __init__(self, operator, timesteps, power=None, power_density=None,
         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 @@ def solver(self, func):
     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 @@ def integrate(self, final_step=True, output=True):
 
         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):

openmc/deplete/chain.py

@@ -698,6 +698,64 @@ def form_matrix(self, rates, fission_yields=None):
         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