Add reactivity control to coupled transport-depletion analyses

openmc/deplete/__init__.py

@@ -17,6 +17,7 @@
 from .results import *
 from .integrators import *
 from .transfer_rates import *
+from .reactivity_control import *
 from . import abc
 from . import cram
 from . import helpers

openmc/deplete/abc.py

@@ -20,7 +20,7 @@
 import numpy as np
 from uncertainties import ufloat
 
-from openmc.checkvalue import check_type, check_greater_than, PathLike
+from openmc.checkvalue import check_value, check_type, check_greater_than, PathLike
 from openmc.mpi import comm
 from openmc import Material
 from .stepresult import StepResult
@@ -29,7 +29,12 @@
 from .pool import deplete
 from .reaction_rates import ReactionRates
 from .transfer_rates import TransferRates
-
+from openmc import Material, Cell
+from .reactivity_control import (
+    GeometricalCellReactivityController,
+    TemperatureCellReactivityController,
+    RefuelMaterialReactivityController
+)
 
 __all__ = [
     "OperatorResult", "TransportOperator",
@@ -563,6 +568,9 @@ class Integrator(ABC):
 
     transfer_rates : openmc.deplete.TransferRates
         Instance of TransferRates class to perform continuous transfer during depletion
+    reactivity_control : openmc.deplete.ReactivityController
+        Instance of ReactivityController class to perform reactivity control during
+        transport-depletion simulation.
 
         .. versionadded:: 0.14.0
 
@@ -653,6 +661,7 @@ def __init__(
         self.source_rates = np.asarray(source_rates)
 
         self.transfer_rates = None
+        self._reactivity_control = None
 
         if isinstance(solver, str):
             # Delay importing of cram module, which requires this file
@@ -701,6 +710,10 @@ def solver(self, func):
 
         self._solver = func
 
+    @property
+    def reactivity_control(self):
+        return self._reactivity_control
+
     def _timed_deplete(self, n, rates, dt, matrix_func=None):
         start = time.time()
         results = deplete(
@@ -787,6 +800,13 @@ def _get_start_data(self):
         return (self.operator.prev_res[-1].time[-1],
                 len(self.operator.prev_res) - 1)
 
+    def _get_bos_from_reactivity_control(self, step_index, bos_conc):
+        """Get BOS from reactivity control."""
+        x = deepcopy(bos_conc)
+        # Get new vector after keff criticality control
+        x, root = self._reactivity_control.search_for_keff(x, step_index)
+        return x, root
+
     def integrate(
             self,
             final_step: bool = True,
@@ -821,10 +841,18 @@ def integrate(
 
                 # Solve transport equation (or obtain result from restart)
                 if i > 0 or self.operator.prev_res is None:
+                    # Update geometry/material according to reactivity control
+                    if self._reactivity_control is not None and source_rate != 0.0:
+                        n, root = self._get_bos_from_reactivity_control(i, n)
+                    else:
+                        root = None
                     n, res = self._get_bos_data_from_operator(i, source_rate, n)
                 else:
                     n, res = self._get_bos_data_from_restart(source_rate, n)
-
+                    if self._reactivity_control:
+                        root = self.operator.prev_res[-1].reac_cont
+                    else:
+                        root = None
                 # Solve Bateman equations over time interval
                 proc_time, n_list, res_list = self(n, res.rates, dt, source_rate, i)
 
@@ -834,9 +862,8 @@ def integrate(
 
                 # Remove actual EOS concentration for next step
                 n = n_list.pop()
-
                 StepResult.save(self.operator, n_list, res_list, [t, t + dt],
-                                source_rate, self._i_res + i, proc_time, path)
+                                source_rate, self._i_res + i, proc_time, root, path)
 
                 t += dt
 
@@ -846,9 +873,13 @@ def integrate(
             # solve)
             if output and final_step and comm.rank == 0:
                 print(f"[openmc.deplete] t={t} (final operator evaluation)")
+            if self._reactivity_control is not None and source_rate != 0.0:
+                n, root = self._get_bos_from_reactivity_control(i+1, n)
+            else:
+                root = None
             res_list = [self.operator(n, source_rate if final_step else 0.0)]
             StepResult.save(self.operator, [n], res_list, [t, t],
-                         source_rate, self._i_res + len(self), proc_time)
+                         source_rate, self._i_res + len(self), proc_time, root)
             self.operator.write_bos_data(len(self) + self._i_res)
 
         self.operator.finalize()
@@ -887,6 +918,33 @@ def add_transfer_rate(
         self.transfer_rates.set_transfer_rate(material, components, transfer_rate,
                                       transfer_rate_units, destination_material)
 
+    def add_reactivity_control(self, obj, attr, **kwargs):
+        """Add reactivity control to integrator scheme.
+
+        Parameters
+        ----------
+        obj : openmc.Cell or openmc.Material object or id or str name
+            Cell or Materials identifier to where add reactivity control
+        attr : str
+            Attribute to specify the type of reactivity control. Accepted values
+            are: 'translation', 'rotation', 'temperature' for an openmc.Cell
+            object; 'refuel' for an openmc.Material object.
+        **kwargs
+            keyword arguments that are passed to ReactivityController.
+
+        """
+        check_value('attribute', attr, ('translation', 'rotation',
+                                        'temperature', 'refuel'))
+        if attr in ('translation', 'rotation'):
+            reactivity_control = GeometricalCellReactivityController
+        elif attr == 'temperature':
+            reactivity_control = TemperatureCellReactivityController
+        elif attr == 'refuel':
+            reactivity_control = RefuelMaterialReactivityController
+
+        self._reactivity_control = reactivity_control.from_params(obj, attr,
+                                                        self.operator,**kwargs)
+
 @add_params
 class SIIntegrator(Integrator):
     r"""Abstract class for the Stochastic Implicit Euler integrators
@@ -1046,13 +1104,13 @@ def integrate(
                 n = n_list.pop()
 
                 StepResult.save(self.operator, n_list, res_list, [t, t + dt],
-                             p, self._i_res + i, proc_time, path)
+                             p, self._i_res + i, proc_time, path=path)
 
                 t += dt
 
             # No final simulation for SIE, use last iteration results
             StepResult.save(self.operator, [n], [res_list[-1]], [t, t],
-                         p, self._i_res + len(self), proc_time, path)
+                         p, self._i_res + len(self), proc_time, path=path)
             self.operator.write_bos_data(self._i_res + len(self))
 
         self.operator.finalize()