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>

openmc/deplete/abc.py

@@ -1013,6 +1013,36 @@ class Integrator(ABC):
             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 @@ loaded from an .xml file and all the nuclides are linked together.
 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 @@ class Chain:
         # 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 @@ class ExternalRates:
         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 @@ class TransferRates(ExternalRates):
             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.

tests/regression_tests/deplete_with_transfer_rates/test.py

@@ -12,7 +12,6 @@ from openmc.deplete import CoupledOperator
 from tests.regression_tests import config, assert_reaction_rates_equal, \
     assert_atoms_equal
 
-
 @pytest.fixture
 def model():
     openmc.reset_auto_ids()
@@ -54,6 +53,9 @@ def model():
     (-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'),
+    (0.0, None, 174.0, 'depletion_with_redox'),
+    (1e-5, None, 174.0, 'depletion_with_removal_and_redox'),
+    (1e-5, 'w', 174.0, 'depletion_with_transfer_and_redox'),
     ])
 def test_transfer_rates(run_in_tmpdir, model, rate, dest_mat, power, ref_result):
     """Tests transfer_rates depletion class with transfer rates"""
@@ -61,13 +63,18 @@ def test_transfer_rates(run_in_tmpdir, model, rate, dest_mat, power, ref_result)
     chain_file = Path(__file__).parents[2] / 'chain_simple.xml'
 
     transfer_elements = ['Xe']
+    os = {'I': -1, 'Xe':0, 'Cs': 1, 'Gd': 3, 'U': 4}
 
     op = CoupledOperator(model, chain_file)
     op.round_number = True
     integrator = openmc.deplete.PredictorIntegrator(
         op, [1], power, timestep_units = 'd')
-    integrator.add_transfer_rate('f', transfer_elements, rate,
-                                 destination_material=dest_mat)
+    if rate != 0.0:
+        integrator.add_transfer_rate('f', transfer_elements, rate,
+                                    destination_material=dest_mat)
+    if 'redox' in ref_result.split('_'):
+        integrator.add_redox('f', {'Gd157':1}, os)
+
     integrator.integrate()
 
     # Get path to test and reference results
@@ -83,9 +90,8 @@ def test_transfer_rates(run_in_tmpdir, model, rate, dest_mat, power, ref_result)
     res_ref = openmc.deplete.Results(path_reference)
     res_test = openmc.deplete.Results(path_test)
 
-    assert_atoms_equal(res_ref, res_test)
-    assert_reaction_rates_equal(res_ref, res_test)
-
+    assert_atoms_equal(res_ref, res_test, tol=1e-3)
+    assert_reaction_rates_equal(res_ref, res_test, tol=1e-3)
 
 @pytest.mark.parametrize("rate, power, ref_result", [
     (1e-1, 0.0, 'no_depletion_with_ext_source'),
@@ -118,5 +124,5 @@ def test_external_source_rates(run_in_tmpdir, model, rate, power, ref_result):
     res_ref = openmc.deplete.Results(path_reference)
     res_test = openmc.deplete.Results(path_test)
 
-    assert_atoms_equal(res_ref, res_test)
-    assert_reaction_rates_equal(res_ref, res_test)
+    assert_atoms_equal(res_ref, res_test, tol=1e-3)
+    assert_reaction_rates_equal(res_ref, res_test, tol=1e-3)

tests/unit_tests/test_deplete_transfer_rates.py

@@ -198,3 +198,35 @@ def test_transfer(run_in_tmpdir, model):
     # 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))
+
+@pytest.mark.parametrize("case_name, buffer, ox", [
+    ('redox', {'Gd157':1}, {'Gd': 3, 'U': 4}),
+    ('buffer_invalid', {'Gd158':1}, {'Gd': 3, 'U': 4}),
+    ('elm_invalid', {'Gd157':1}, {'Gb': 3, 'U': 4}),
+    ])
+def test_redox(case_name, buffer, ox, model):
+    op = CoupledOperator(model, CHAIN_PATH)
+    number_of_timesteps = 2
+    transfer = TransferRates(op, model.materials, number_of_timesteps)
+
+    # Test by Openmc material, material name and material id
+    material, dest_material, dest_material2 = [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]:
+
+            if case_name == 'buffer_invalid':
+                with pytest.raises(ValueError, match='Gd158 is not a valid '
+                                   'nuclide.'):
+                    transfer.set_redox(material_input, buffer, ox)
+
+            elif case_name == 'elm_invalid':
+                with pytest.raises(ValueError, match='Gb is not a valid '
+                                   'element.'):
+                    transfer.set_redox(material_input, buffer, ox)
+            else:
+                transfer.set_redox(material_input, buffer, ox)
+                mat_id = transfer._get_material_id(material_input)
+                assert transfer.redox[mat_id][0] == buffer
+                assert transfer.redox[mat_id][1] == ox