Merge branch 'openmc-dev:develop' into ttb

Author: GuySten

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

openmc/deplete/abc.py

@@ -29,7 +29,7 @@
     _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 @@ def __init__(
         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 @@ def solver(self, func):
 
         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 @@ def integrate(
         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 @@ def add_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.
+        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 @@ def _get_bos_data_from_operator(self, step_index, step_power, n_bos):
         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

openmc/deplete/chain.py

@@ -703,7 +703,7 @@ def form_matrix(self, rates, fission_yields=None):
         # 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 @@ def form_rr_term(self, tr_rates, mats):
         ----------
         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 @@ def form_rr_term(self, tr_rates, mats):
 
         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