depletion-thermochemistry: Redox control transfer rates

openmc/deplete/abc.py

@@ -847,6 +847,24 @@ def add_transfer_rate(self, material, components, transfer_rate,
         self.transfer_rates.set_transfer_rate(material, components, transfer_rate,
                                       transfer_rate_units, destination_material)
 
+    def add_redox(self, material, buffer, oxidation_states):
+        """Add redox control to depletable material.
+
+        Parameters
+        ----------
+        material : openmc.Material or str or int
+            Depletable material
+        buffer : dict
+            Dictionary of buffer nuclides to be added to keep redox constant,
+            where keys are nuclide names and values fractions to 1.
+        oxidation_states : dict
+            User-defined oxidation states for elements.
+        """
+        if self.transfer_rates is None:
+            self.transfer_rates = TransferRates(self.operator, self.operator.model)
+
+        self.transfer_rates.set_redox(material, buffer, oxidation_states)
+
 @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 os
 import re
 from collections import defaultdict, namedtuple
@@ -686,6 +687,65 @@ def form_matrix(self, rates, fission_yields=None):
         dict.update(matrix_dok, matrix)
         return matrix_dok.tocsr()
 
+    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_b` can be written
+        as: :math:`\frac{dN_b(t)}{dt} =
+                \cdots + \frac{1}{Ox_b}\sum_i N_i\left( L_{ii}Ox_i -
+                \sum_j G_{i\rightarrow j } Ox_j\right)`
+
+        where :math:`Ox_b` and :math:`Ox_j` are the oxidation states for the
+        corresponding buffer elmenent and j-th nuclide.
+        The first term in the right hand side represent the losses in the
+        diagonal terms of the Bateman matrix, for each nuclide :math:`i`, and
+        the second one the gains in the off-diagonal terms, multiplied by their
+        respective oxidation states.
+
+        Parameters
+        ----------
+        matrix : scipy.sparse.csr_matrix
+            Sparse matrix representing depletion
+        buffer : dict
+            Dictionary of buffer nuclides to be added to keep redox constant,
+            where keys are nuclide names and values fractions to 1.
+        oxidation_states : dict
+            User-defined oxidation states for elements.
+        Returns
+        -------
+        matrix : scipy.sparse.csr_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
+        ox = 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 = np.array([])
+
+        # calculate the redox array
+        for i in range(len(self)):
+            # Gains: multiply each column of the depletion matrix by the oxidation states
+            # vector, excluding the i-th terms
+            gains = np.concatenate((array[:i,i], array[i+1:,i])) * np.delete(ox, i)
+            # Loss: multiply the i-th term by its corresponding oxidation state value
+            loss = array[i,i] * ox[i]
+            # Calculte the redox term
+            redox = np.append(redox, loss + sum(gains))
+
+        # 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 * buffer[nuc] / ox[idx]
+
+        return sp.dok_matrix(array)
+
     def form_rr_term(self, tr_rates, mats):
         """Function to form the transfer rate term matrices.
 

openmc/deplete/pool.py

@@ -97,10 +97,18 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None,
         # Calculate transfer rate terms as diagonal matrices
         transfers = map(chain.form_rr_term, repeat(transfer_rates),
                         transfer_rates.local_mats)
+
         # Subtract transfer rate terms from Bateman matrices
         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 len(transfer_rates.index_transfer) > 0:
             # Gather all on comm.rank 0
             matrices = comm.gather(matrices)
@@ -112,10 +120,16 @@ def deplete(func, chain, n, rates, dt, matrix_func=None, transfer_rates=None,
                 n = [n_elm for n_mat in n for  n_elm in n_mat]
 
                 # Calculate transfer rate terms as diagonal matrices
-                transfer_pair = {
-                    mat_pair: chain.form_rr_term(transfer_rates, mat_pair)
-                    for mat_pair in transfer_rates.index_transfer
-                }
+                transfer_pair = dict()
+                for mat_pair in transfer_rates.index_transfer:
+                    transfer_matrix = chain.form_rr_term(transfer_rates, 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
                 # to be solved in one go

openmc/deplete/transfer_rates.py

@@ -45,10 +45,12 @@ def __init__(self, operator, model):
         self.materials = model.materials
         self.burnable_mats = operator.burnable_mats
         self.local_mats = operator.local_mats
+        self.chain_nuclides = [nuc.name for nuc in operator.chain.nuclides]
 
         #initialize transfer rates container dict
         self.transfer_rates = {mat: {} for mat in self.burnable_mats}
         self.index_transfer = set()
+        self.redox = dict()
 
     def _get_material_id(self, val):
         """Helper method for getting material id from Material obj or name.
@@ -227,3 +229,28 @@ def set_transfer_rate(self, material, components, transfer_rate,
                     (transfer_rate / unit_conv, destination_material_id)]
             if destination_material_id is not None:
                 self.index_transfer.add((destination_material_id, material_id))
+
+    def set_redox(self, material, buffer, oxidation_states):
+        """Add redox control to depletable material.
+
+        Parameters
+        ----------
+        material : openmc.Material or str or int
+            Depletable material
+        buffer : dict
+            Dictionary of buffer nuclides to be added to keep redox constant,
+            where keys are nuclide names and values fractions to 1.
+        oxidation_states : dict
+            User-defined oxidation states for elements.
+        """
+        material_id = self._get_material_id(material)
+        #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)

tests/regression_tests/deplete_with_transfer_rates/test.py

@@ -40,9 +40,9 @@ def model():
     geometry = openmc.Geometry([cell_f, cell_w])
 
     settings = openmc.Settings()
-    settings.particles = 500
+    settings.particles = 100
     settings.inactive = 0
-    settings.batches = 2
+    settings.batches = 10
 
     return openmc.Model(geometry, materials, settings)
 
@@ -54,19 +54,27 @@ 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"""
 
     chain_file = Path(__file__).parents[2] / 'chain_simple.xml'
 
     transfer_elements = ['Xe']
+    ox = {'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}, ox)
     integrator.integrate()
 
     # Get path to test and reference results
@@ -82,5 +90,5 @@ 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, 1e-4)
+    assert_atoms_equal(res_ref, res_test, 1e-6)
     assert_reaction_rates_equal(res_ref, res_test)

tests/unit_tests/test_deplete_transfer_rates.py

@@ -186,3 +186,34 @@ 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)
+    transfer = TransferRates(op, model)
+
+    # 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