Fix coupled external source rate and transfer rate with destination material (#3959)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

docs/source/methods/depletion.rst

@@ -339,3 +339,56 @@ where:
 
 Note that mass conservation is guaranteed by transferring the number
 of atoms directly.
+
+---------------------
+External Source Rates
+---------------------
+
+OpenMC allows the addition of external source rates to the depletion matrix.
+This is useful for modeling the feed or removal of fixed amounts of nuclides
+to or from a depletable material. Rates are specified as a mass flow (default
+units of [g/s]) and converted to [atom/s] source terms for the nuclides in the
+composition vector. A positive rate corresponds to feed and a negative rate to
+removal.
+
+Mathematically, this is represented as an external source term :math:`S_i` to
+the depletion equation:
+
+.. math::
+
+  \frac{dN_i}{dt} = \sum_j A_{ij} N_j + S_i
+
+The resulting linear system is non-homogeneous but can be recast in homogeneous
+form by augmenting the nuclide vector with a constant component equal to unity:
+
+.. math::
+
+  \frac{d}{dt}\begin{bmatrix}\mathbf{n}\\ 1\end{bmatrix} =
+  \begin{bmatrix}
+        \mathbf{A} & \mathbf{s}\\
+        \mathbf{0} & 0
+  \end{bmatrix}
+  \begin{bmatrix}
+        \mathbf{n}\\
+        1
+  \end{bmatrix}
+
+where :math:`\mathbf{s}` is the vector of external source rates in [atom/s]. The
+resulting system can be solved with the same integration algorithms that are
+used in the absence of the external source term.
+
+External source rates with transfer rates coupling materials
+------------------------------------------------------------
+
+In the presence of external source rates and coupled transfer rates between
+materials, the off-diagonal transfer blocks must match the dimensions of the
+corresponding diagonal blocks. Using the transfer example above, if an external
+source rate is applied to material 1 (the losing material), the coupling matrix
+:math:`\mathbf{T_{21}}` is extended with an additional column of zeroes so that
+its column count matches the augmented :math:`\mathbf{A_{11}}`. If the external
+source is applied to material 2 (the receiving material),
+:math:`\mathbf{T_{21}}` instead receives an additional row of zeroes so that its
+row count matches the augmented :math:`\mathbf{A_{22}}`.
+
+
+

docs/source/usersguide/depletion.rst

@@ -450,6 +450,67 @@ to transfer xenon from one material to another, you'd use::
 
     integrator.add_transfer_rate(mat1, ['Xe'], 0.1, destination_material=mat2)
 
+External Source Rates
+=====================
+
+External source rates define a fixed mass feed or removal of nuclides to or
+from a depletable material. Unlike transfer rates, which are proportional to
+the instantaneous nuclide inventory, external source rates add a constant
+source term to the depletion equations. This can be useful to model batch
+refueling, makeup fuel addition, or fixed-rate off-gas removal.
+
+External source rates are defined by calling the
+:meth:`~openmc.deplete.abc.Integrator.add_external_source_rate()` method
+directly from one of the Integrator classes::
+
+    ...
+    integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)
+    integrator.add_external_source_rate(...)
+
+Defining external source rates
+------------------------------
+
+The :meth:`~openmc.deplete.abc.Integrator.add_external_source_rate()` method
+requires a :class:`~openmc.Material` instance (or a material ID or the name) as
+the depletable material, a composition dictionary giving the relative weight
+fractions of elements and/or nuclides in the feed or removal stream, and a mass
+flow rate.
+
+.. caution::
+
+   Make sure you set the rate value with the right sign. A positive rate
+   corresponds to feed, while a negative rate corresponds to removal. This is
+   the opposite convention used for transfer rates.
+
+The ``rate_units`` argument specifies the units for the mass flow rate. The
+default is ``g/s``, but ``g/min``, ``g/h``, ``g/d``, and ``g/a`` are also valid
+options.
+
+For example, to feed U235 into a material at 10 g/day, you'd use::
+
+    mat = openmc.Material()
+    ...
+
+    integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)
+    composition = {'U235': 1.0}
+    integrator.add_external_source_rate(mat1, composition, 10, rate_units='g/d')
+
+Composition keys may be nuclides (e.g., ``'U235'``) or naturally abundant
+elements (e.g., ``'U'``). When an element is specified, the mass flow is
+distributed across its naturally occurring isotopes according to their natural
+abundances.
+
+The optional ``timesteps`` argument restricts the external source rate to
+specific depletion step indices. If omitted, the rate is applied at every step.
+
+Combining with transfer rates
+-----------------------------
+
+External source rates can be used together with transfer rates, including
+transfers between materials via ``destination_material``. See
+:ref:`methods_depletion` for the augmented-matrix formulation used when both
+features are active.
+
 Comparing to Other Codes
 ========================
 

openmc/deplete/abc.py

@@ -1028,11 +1028,6 @@ class Integrator(ABC):
             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)
@@ -1074,11 +1069,6 @@ class Integrator(ABC):
             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)
 

openmc/deplete/pool.py

@@ -6,10 +6,10 @@ from itertools import repeat, starmap
 from multiprocessing import Pool
 
 import numpy as np
-from scipy.sparse import hstack
+from scipy.sparse import hstack, vstack
 
 from openmc.mpi import comm
-from .._sparse_compat import block_array
+from .._sparse_compat import block_array, csc_array
 
 # Configurable switch that enables / disables the use of
 # multiprocessing routines during depletion
@@ -41,6 +41,29 @@ def _distribute(items):
             return items[j:j + chunk_size]
         j += chunk_size
 
+
+def _add_external_source(
+    matrices, n, chain, external_source_rates, current_timestep
+):
+    """Augment depletion matrices and nuclide vectors with external sources."""
+    sources = map(chain.form_ext_source_term, repeat(external_source_rates),
+                  repeat(current_timestep), external_source_rates.local_mats)
+    matrices = [
+        hstack([matrix, source])
+        for matrix, source in zip(matrices, sources)
+    ]
+    n_solve = [arr.copy() for arr in n]
+
+    # Homogenize the augmented matrices and nuclide vectors
+    for i, matrix in enumerate(matrices):
+        if matrix.shape[0] + 1 == matrix.shape[1]:
+            matrices[i] = vstack(
+                [matrix, csc_array((1, matrix.shape[1]))])
+            n_solve[i] = np.append(n_solve[i], 1.0)
+
+    return matrices, n_solve
+
+
 def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
             transfer_rates=None, external_source_rates=None, substeps=1,
             *matrix_args):
@@ -62,7 +85,7 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
         Time in [s] to deplete for
     current_timestep : int
         Current timestep index
-    maxtrix_func : callable, optional
+    matrix_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
@@ -87,7 +110,6 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
         list contains the number of [atom] of each nuclide.
 
     """
-
     fission_yields = chain.fission_yields
     if len(fission_yields) == 1:
         fission_yields = repeat(fission_yields[0])
@@ -103,8 +125,14 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
         matrices = map(matrix_func, repeat(chain), rates, fission_yields,
                        *matrix_args)
 
-    if (transfer_rates is not None and
-        current_timestep in transfer_rates.external_timesteps):
+    # Determine if transfer rates or external source rates are active
+    transfer_active = transfer_rates is not None and \
+        current_timestep in transfer_rates.external_timesteps
+    external_active = external_source_rates is not None and \
+        current_timestep in external_source_rates.external_timesteps
+
+    n_solve = n
+    if transfer_active:
         # Calculate transfer rate terms as diagonal matrices
         transfers = map(chain.form_rr_term, repeat(transfer_rates),
                         repeat(current_timestep), transfer_rates.local_mats)
@@ -120,10 +148,16 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
                                                 transfer_rates.redox[mat_id][0],
                                                 transfer_rates.redox[mat_id][1])
 
+        # Add external sources if present
+        if external_active:
+            matrices, n_solve = _add_external_source(
+                matrices, n, chain, external_source_rates, current_timestep)
+
+        # Set transfer rate terms with destination material if present
         if current_timestep in transfer_rates.index_transfer:
             # Gather all on comm.rank 0
             matrices = comm.gather(matrices)
-            n = comm.gather(n)
+            n = comm.gather(n_solve)
 
             if comm.rank == 0:
                 # Expand lists
@@ -132,20 +166,27 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
 
                 # Calculate transfer rate terms as diagonal matrices
                 transfer_pair = {}
-                for mat_pair in transfer_rates.index_transfer[current_timestep]:
+                for mat_pair in dict.fromkeys(transfer_rates.index_transfer[current_timestep]):
                     transfer_matrix = chain.form_rr_term(transfer_rates,
                                                          current_timestep,
                                                          mat_pair)
-
                     # check if destination material has a redox control
                     if mat_pair[0] in transfer_rates.redox:
                         transfer_matrix = chain.add_redox_term(transfer_matrix,
                                           transfer_rates.redox[mat_pair[0]][0],
                                           transfer_rates.redox[mat_pair[0]][1])
+                    # Add external source rates if present
+                    if external_active:
+                        if len(external_source_rates.get_components(mat_pair[0], current_timestep)) > 0:
+                            transfer_matrix = vstack([transfer_matrix,
+                                csc_array((1, transfer_matrix.shape[1]))])
+                        if len(external_source_rates.get_components(mat_pair[1], current_timestep)) > 0:
+                            transfer_matrix = hstack([transfer_matrix,
+                                csc_array((transfer_matrix.shape[0], 1))])
                     transfer_pair[mat_pair] = transfer_matrix
 
-                # Combine all matrices together in a single matrix of matrices
-                # to be solved in one go
+                # Combine all matrices together in a single block matrix of matrices
+                # to be solved on one rank
                 n_rows = n_cols = len(transfer_rates.burnable_mats)
                 rows = []
                 for row in range(n_rows):
@@ -171,37 +212,25 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
 
                 # Split back the nuclide vector result into the original form
                 n_result = np.split(n_result, np.cumsum([len(i) for i in n])[:-1])
-
             else:
                 n_result = None
 
-            # Braodcast result to other ranks
+            # Broadcast result to other MPI ranks and then distribute
             n_result = comm.bcast(n_result)
-            # Distribute results across MPI
             n_result = _distribute(n_result)
 
+            # Remove extra values based on the materials local to each rank
+            if external_active:
+                external_source_rates.reformat_nuclide_vectors(n_result)
+
             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)
+    # If only external source rates are present
+    elif external_active:
+        matrices, n_solve = _add_external_source(
+            matrices, n, chain, external_source_rates, current_timestep)
 
-        # 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):
-                matrix.resize(matrix.shape[1], matrix.shape[1])
-                n[i] = np.append(n[i], 1.0)
-
-    inputs = zip(matrices, n, repeat(dt), repeat(substeps))
+    inputs = zip(matrices, n_solve, repeat(dt), repeat(substeps))
 
     if USE_MULTIPROCESSING:
         with Pool(NUM_PROCESSES) as pool:
@@ -209,10 +238,8 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=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)
+    # Remove extra value at the end of the nuclide vectors if external source rates are present
+    if external_active:
         external_source_rates.reformat_nuclide_vectors(n_result)
 
     return n_result

tests/regression_tests/deplete_with_transfer_rates/test.py

@@ -126,3 +126,38 @@ def test_external_source_rates(run_in_tmpdir, model, rate, power, ref_result):
 
     assert_atoms_equal(res_ref, res_test, tol=1e-3)
     assert_reaction_rates_equal(res_ref, res_test, tol=1e-3)
+
+@pytest.mark.parametrize("external_source_rate, transfer_rate, power, ref_result", [
+    (1e-1, 1e-1, 174., 'depletion_with_ext_source_and_transfer'),
+])
+def test_external_source_rates_with_transfer_rates(run_in_tmpdir, model, external_source_rate, 
+    transfer_rate, power, ref_result):
+    """Tests external_rates depletion class with external source rates and transfer 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, external_source_rate)
+    integrator.add_transfer_rate('f', ['U235'], transfer_rate, destination_material='w')
+    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, tol=1e-3)
+    assert_reaction_rates_equal(res_ref, res_test, tol=1e-3)