depletion-thermochemistry: Redox control transfer rates (#2783)

Co-authored-by: Gavin Ridley <gavin.keith.ridley@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Author: 3 people

openmc/deplete/abc.py

@@ -1013,6 +1013,36 @@ def add_external_source_rate(
             material, composition, rate, rate_units, timesteps)
 
 
+    def add_redox(self, material, buffer, oxidation_states, timesteps=None):
+        """Add redox control to depletable material.
+
+        Parameters
+        ----------
+        material : openmc.Material or str or int
+            Depletable material
+        buffer : dict
+            Dictionary of buffer nuclides used to maintain redox balance. Keys
+            are nuclide names (strings) and values are their respective
+            fractions (float) that collectively sum to 1.
+        oxidation_states : dict
+            User-defined oxidation states for elements. Keys are element symbols
+            (e.g., 'H', 'He'), and values are their corresponding oxidation
+            states as integers (e.g., +1, 0).
+        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.transfer_rates is None:
+            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))
+
+        self.transfer_rates.set_redox(material, buffer, oxidation_states, timesteps)
+
 @add_params
 class SIIntegrator(Integrator):
     r"""Abstract class for the Stochastic Implicit Euler integrators

openmc/deplete/chain.py

@@ -7,6 +7,7 @@
 from io import StringIO
 from itertools import chain
 import math
+import numpy as np
 import re
 from collections import defaultdict, namedtuple
 from collections.abc import Mapping, Iterable
@@ -714,6 +715,59 @@ def setval(i, j, val):
         # Return CSC representation instead of DOK
         return sp.csc_matrix((vals, (rows, cols)), shape=(n, n))
 
+    def add_redox_term(self, matrix, buffer, oxidation_states):
+        """Adds a redox term to the depletion matrix from data contained in
+        the matrix itself and a few user-inputs.
+
+        The redox term to add to the buffer nuclide :math:`N_j` can be written
+        as: :math:`\frac{dN_j(t)}{dt} =
+                \cdots - \frac{1}{OS_j}\sum_i N_i a_{ij} \cdot OS_i `
+
+        where :math:`OS` is the oxidation states vector and `a_{ij}` the
+        corresponding term in the Bateman matrix.
+
+        Parameters
+        ----------
+        matrix : scipy.sparse.csc_matrix
+            Sparse matrix representing depletion
+        buffer : dict
+            Dictionary of buffer nuclides used to maintain anoins net balance.
+            Keys are nuclide names (strings) and values are their respective
+            fractions (float) that collectively sum to 1.
+        oxidation_states : dict
+            User-defined oxidation states for elements. Keys are element symbols
+            (e.g., 'H', 'He'), and values are their corresponding oxidation
+            states as integers (e.g., +1, 0).
+        Returns
+        -------
+        matrix : scipy.sparse.csc_matrix
+            Sparse matrix with redox term added
+        """
+        # Elements list with the same size as self.nuclides
+        elements = [re.split(r'\d+', nuc.name)[0] for nuc in self.nuclides]
+
+        # Match oxidation states with all elements and add 0 if not data
+        os = np.array([oxidation_states[elm] if elm in oxidation_states else 0
+                       for elm in elements])
+
+        # Buffer idx with nuclide index as value
+        buffer_idx = {nuc: self.nuclide_dict[nuc] for nuc in buffer}
+        array = matrix.toarray()
+        redox_change = np.array([])
+
+        # calculate the redox array
+        for i in range(len(self)):
+            # Net redox impact of reaction: multiply the i-th column of the
+            # depletion matrix by the oxidation states
+            redox_change = np.append(redox_change, sum(array[:, i]*os))
+
+        # Subtract redox vector to the buffer nuclides in the matrix scaling by
+        # their respective oxidation states
+        for nuc, idx in buffer_idx.items():
+            array[idx] -= redox_change * buffer[nuc] / os[idx]
+
+        return sp.csc_matrix(array)
+
     def form_rr_term(self, tr_rates, current_timestep, mats):
         """Function to form the transfer rate term matrices.
 

openmc/deplete/pool.py

@@ -109,6 +109,13 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
         matrices = [matrix - transfer for (matrix, transfer) in zip(matrices,
                                                                     transfers)]
 
+        if transfer_rates.redox:
+            for mat_idx, mat_id in enumerate(transfer_rates.local_mats):
+                if mat_id in transfer_rates.redox:
+                    matrices[mat_idx] = chain.add_redox_term(matrices[mat_idx],
+                                                transfer_rates.redox[mat_id][0],
+                                                transfer_rates.redox[mat_id][1])
+
         if current_timestep in transfer_rates.index_transfer:
             # Gather all on comm.rank 0
             matrices = comm.gather(matrices)
@@ -125,6 +132,12 @@ def deplete(func, chain, n, rates, dt, current_timestep=None, matrix_func=None,
                     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])
                     transfer_pair[mat_pair] = transfer_matrix
 
                 # Combine all matrices together in a single matrix of matrices

openmc/deplete/transfer_rates.py

@@ -49,9 +49,10 @@ def __init__(self, operator, materials, number_of_timesteps):
         self.local_mats = operator.local_mats
         self.number_of_timesteps = number_of_timesteps
 
-        # initialize transfer rates container dict
+        #initialize transfer rates container dict
         self.external_rates = {mat: defaultdict(list) for mat in self.burnable_mats}
         self.external_timesteps = []
+        self.redox = {}
 
     def _get_material_id(self, val):
         """Helper method for getting material id from Material obj or name.
@@ -300,6 +301,46 @@ def set_transfer_rate(self, material, components, transfer_rate,
             self.external_timesteps = np.unique(np.concatenate(
                     [self.external_timesteps, timesteps]))
 
+    def set_redox(self, material, buffer, oxidation_states, timesteps=None):
+        """Add redox control to depletable material.
+
+        Parameters
+        ----------
+        material : openmc.Material or str or int
+            Depletable material
+        buffer : dict
+            Dictionary of buffer nuclides used to maintain redox balance.
+            Keys are nuclide names (strings) and values are their respective
+            fractions (float) that collectively sum to 1.
+        oxidation_states : dict
+            User-defined oxidation states for elements.
+            Keys are element symbols (e.g., 'H', 'He'), and values are their
+            corresponding oxidation states as integers (e.g., +1, 0).
+        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.
+
+        """
+        material_id = self._get_material_id(material)
+        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)
+        #Check nuclides in buffer exist
+        for nuc in buffer:
+            if nuc not in self.chain_nuclides:
+                raise ValueError(f'{nuc} is not a valid nuclide.')
+        # Checks element in oxidation states exist
+        for elm in oxidation_states:
+            if elm not in ELEMENT_SYMBOL.values():
+                raise ValueError(f'{elm} is not a valid element.')
+
+        self.redox[material_id] =  (buffer, oxidation_states)
+        self.external_timesteps = np.unique(np.concatenate(
+                    [self.external_timesteps, timesteps]))
 
 class ExternalSourceRates(ExternalRates):
     """Class for defining external source rates.