Material addition

openmc/deplete/abc.py

@@ -943,6 +943,9 @@ def add_density_function(self, mats, density_func, oxidation_states):
     def add_redox(self, mat, buffer, oxidation_states):
         self.transfer_rates.set_redox(mat, buffer, oxidation_states)
 
+    def add_material(self, mat, value, mat_vector, timestep, quantity='grams'):
+        self.batchwise.add_material(mat, value, mat_vector, timestep,
+                                    quantity)
 @add_params
 class SIIntegrator(Integrator):
     r"""Abstract class for the Stochastic Implicit Euler integrators

openmc/deplete/batchwise.py

@@ -126,6 +126,8 @@ def __init__(self, operator, model, bracket, bracket_limit,
 
         self.print_iterations = print_iterations
         self.search_for_keff_output = search_for_keff_output
+        # List of materials to add as single action in time
+        self.materials_to_add = {}
 
     @property
     def bracketed_method(self):
@@ -384,7 +386,7 @@ def _update_densities(self):
                 else:
                     print(f'Mat ID: {mat_id}, with comp: {mol_comp}, returned nan density')
 
-    def _update_materials(self, x):
+    def _update_materials(self, x, step_index):
         """
         Update number density and material compositions in OpenMC on all processes.
         If density_treatment is set to 'constant-density'
@@ -399,6 +401,11 @@ def _update_materials(self, x):
         """
         self.operator.number.set_density(x)
 
+        if self.materials_to_add:
+            for (mat, t), values in self.materials_to_add.items():
+                if t == step_index:
+                    x = self._add_material(x, mat, values)
+
         if self.density_functions:
             self._update_densities()
             self._update_volumes()
@@ -425,6 +432,8 @@ def _update_materials(self, x):
                 # Update densities on C API side
                 openmc.lib.materials[int(mat_id)].set_densities(nuclides, densities)
 
+        return x
+
     def _update_x_and_set_volumes(self, x, volumes):
         """
         Update x vector with new volumes, before assign them to AtomNumber.
@@ -506,7 +515,6 @@ def _get_molar_composition(self, mat_mol_id):
 
         return mol_comp
 
-
     def _get_redox(self, mat_rx_id):
         """
         Parameters
@@ -552,6 +560,60 @@ def _balance_redox(self, x):
                 x[mat_idx][nuc_idx] = number_i[mat_id, nuc]
         return x
 
+    def add_material(self, mat, value, mat_vector, timestep, units='grams'):
+        """
+        Add material entries to self.materials_to_add dictionary for later use,
+        before converting grams into atoms.
+        Parameters
+        ----------
+        mat : openmc.Material or str or int
+            Material identifier
+        value : float
+            Material to add in units of 'units'
+        mat_vector : dict of float
+            Nuclide composition of material to add where keys are the nuclides
+            and values the fractions.
+        units : str
+            Units of material value.
+            Default : grams
+        """
+        # convert grams in atoms
+        self.materials_to_add[(mat,timestep)] = dict()
+        for nuc,frac in mat_vector.items():
+            self.materials_to_add[(mat,timestep)][nuc] = \
+                    frac * value / atomic_mass(nuc) * AVOGADRO
+
+    def _add_material(self, x, mat, values):
+        """
+        Private method to add materials as number of atoms to number and x.
+        Parameters
+        ----------
+        x : list of numpy.ndarray
+            Total atoms concentrations
+        mat : openmc.Material or str or int
+            Material identifier
+        values : dict of floats
+            Numebr of atoms per nuclide to add, where keys are nuclides and
+            values number of atoms.
+        Return
+        ------
+        x : list of numpy.ndarray
+            Total atoms concentrations
+        """
+        number_i = self.operator.number
+        #get material id
+        mats = self._get_materials([mat])
+
+        for mat_idx, mat_id in enumerate(self.local_mats):
+            if mat_id in mats:
+                for nuc, val in values.items():
+                    number_i[mat_id, nuc] += val
+                    #Now updates x vector
+                    if nuc in number_i.burnable_nuclides:
+                        nuc_idx = number_i.index_nuc[nuc]
+                        x[mat_idx][nuc_idx] += val
+        return x
+
 class BatchwiseCell(Batchwise):
     """Abstract class holding batch wise cell-based functions.
 
@@ -732,7 +794,7 @@ def search_for_keff(self, x, step_index):
         # Update all material densities from concentration vectors
         #before performing the search_for_keff. This is needed before running
         # the  transport equations in the search_for_keff algorithm
-        super()._update_materials(x)
+        x = super()._update_materials(x, step_index)
 
         # Calculate new cell attribute
         root = super()._search_for_keff(val)
@@ -741,8 +803,9 @@ def search_for_keff(self, x, step_index):
         self._set_cell_attrib(root)
 
         # if at least one of the cell materials is depletable, calculate new
-        # volume and update x and number accordingly
-        # new volume
+        # volume and update x and number accordingly.
+        # Alternatively, if material has been edded x needs to be updated
+        # TODO: improve this
         if self.cell_materials:
             volumes = self._calculate_volumes()
             x = super()._update_x_and_set_volumes(x, volumes)
@@ -1121,7 +1184,7 @@ def search_for_keff(self, x, step_index):
         """
         # Update AtomNumber with new conc vectors x. Materials are also updated
         # even though they are also re-calculated when running the search_for_kef
-        super()._update_materials(x)
+        x = super()._update_materials(x, step_index)
 
         # Solve search_for_keff and find new value
         root = super()._search_for_keff(0)
@@ -1872,6 +1935,9 @@ def set_density_function(self, mats, density_func, oxidation_states):
         for bw in self.bw_list:
             bw.set_density_function(mats, density_func, oxidation_states)
 
+    def add_material(self, mat, value, mat_vector, timestep, units):
+        self.bw_geom.add_material(mat, value, mat_vector, timestep, units)
+
     def update_from_restart(self, x, root):
         """
         This is for a restart. TODO update abc class
@@ -1977,6 +2043,9 @@ def set_density_function(self, mats, density_func, oxidation_states):
         for bw in self.bw_list:
             bw.set_density_function(mats, density_func, oxidation_states)
 
+    def add_material(self, mat, value, mat_vector, timestep, units):
+        self.bw_geom.add_material(mat, value, mat_vector, timestep, units)
+
     def update_from_restart(self, x, root):
         """
         This is for a restart. TODO update abc class