diff --git a/python/cardinal_operator.py b/python/cardinal_operator.py
new file mode 100644
index 000000000..debb427b0
--- /dev/null
+++ b/python/cardinal_operator.py
@@ -0,0 +1,301 @@
+import os
+import copy
+import logging
+from numbers import Integral
+import shlex
+
+from uncertainties import ufloat
+
+import openmc
+import openmc.lib
+
+# openmc
+import openmc.deplete
+from openmc.deplete import CoupledOperator
+from openmc.deplete.abc import OperatorResult
+
+# moose
+from MooseControl import MooseControl
+
+logger = logging.getLogger('CardinalOperator')
+
+class CardinalOperator(CoupledOperator):
+
+ def __init__(self, *args, **kwargs):
+ self._control = None
+ self.cardinal_cmd = None
+ super().__init__(*args, **kwargs)
+
+ def __del__(self):
+ # Ensure the Cardinal process is terminated
+ # when this object is deleted
+ self.stop_cardinal()
+
+ @staticmethod
+ def material_userobject_args(mat):
+ uo_path = f'UserObjects/openmc_mat{mat.id}'
+ return [f'{uo_path}/type=OpenMCNuclideDensities',
+ f'{uo_path}/material_id={mat.id}',
+ f'{uo_path}/names="{" ".join(mat.nuclides)}"',
+ f'{uo_path}/densities="{" ".join([str(d) for d in mat.densities.flat])}"']
+
+ @staticmethod
+ def filter_userobject_args(id, domain_type, domain_ids):
+ filter_path = f'UserObjects/openmc_filter{id}'
+ return [f'{filter_path}/type=OpenMCDomainFilterEditor',
+ f'{filter_path}/create_filter=true',
+ f'{filter_path}/filter_id="{id}"',
+ f'{filter_path}/filter_type="{domain_type}"',
+ f'{filter_path}/bins="{" ".join([str(d) for d in domain_ids])}"']
+
+ @staticmethod
+ def tally_userobject_args(id):
+ tally_path = f'UserObjects/openmc_tally{id}'
+ return [f'{tally_path}/type=OpenMCTallyEditor',
+ f'{tally_path}/create_tally=true',
+ f'{tally_path}/tally_id={id}',
+ f'{tally_path}/scores=""',
+ f'{tally_path}/nuclides=""',
+ f'{tally_path}/filter_ids=""',
+ f'{tally_path}/multiply_density=false']
+
+ def start_cardinal(self):
+ """Builds and starts the MooseControl object that runs cardinal
+ """
+ if self.cardinal_cmd is None:
+ raise Exception('cardinal_cmd was not set')
+ cardinal_cmd = shlex.split(self.cardinal_cmd)
+
+ # add server information
+ control_name = 'web_server'
+ control_path = f'Controls/{control_name}'
+
+ cardinal_cmd += [f'{control_path}/type=WebServerControl',
+ f'{control_path}/execute_on="TIMESTEP_BEGIN TIMESTEP_END"']
+ # add material user objects
+ for mat in openmc.lib.materials.values():
+ cardinal_cmd += self.material_userobject_args(mat)
+ # add tally user objects
+ for tally in openmc.lib.tallies.values():
+ cardinal_cmd += self.tally_userobject_args(tally.id)
+
+ for filter in openmc.lib.filters.values():
+ if isinstance(filter, (openmc.lib.CellFilter, openmc.lib.MaterialFilter, openmc.lib.UniverseFilter)):
+ bins = [bin.id for bin in filter.bins]
+ elif isinstance(filter, openmc.lib.MeshFilter):
+ bins = [filter.mesh.id]
+ else:
+ continue
+ cardinal_cmd += self.filter_userobject_args(filter.id, filter.filter_type, bins)
+
+ with open('server.i', 'w') as fh:
+ fh.write('\n'.join(cardinal_cmd))
+
+ self._control = MooseControl(moose_command=cardinal_cmd,
+ moose_control_name=control_name)
+ self._control.initialize()
+
+ def stop_cardinal(self):
+ """Stops the Cardinal background process"""
+ if self._control is not None:
+ self._control.kill()
+
+ def run_cardinal(self):
+ """Execute a Cardinal run and wait for it to finish
+ """
+ logger.info('Running cardinal')
+
+ # Wait for cardinal to be a the beginning of a solve
+ self._control.wait('TIMESTEP_BEGIN')
+
+ # Start the solve
+ self._control.setContinue()
+
+ # Wait for cardinal to finish the solve
+ self._control.wait('TIMESTEP_END')
+
+ return f'statepoint.{self.model.settings.batches}.h5'
+
+ def update_cardinal_materials(self):
+ """Pass material compositions from the local in-memory OpenMC model to
+ Cardinal's OpenMC model
+ """
+ logger.info('Waiting to update carindal materials')
+
+ # Wait for cardinal to be at any flag. On the first solve, we'll be
+ # sitting at TIMESTEP_BEGIN. But after any other solves, we'll be
+ # at TIMESTEP_END from the previous timestep
+ current_flag = self._control.wait()
+
+ # If we're at TIMESTEP_END, we need to continue on
+ if current_flag == 'TIMESTEP_END':
+ self._control.setContinue()
+
+ # Wait for cardinal to be ready at the beginning of a timestep
+ self._control.wait('TIMESTEP_BEGIN')
+
+ logger.info('Updating cardinal materials')
+ for m in openmc.lib.materials.values():
+ uo_path = f'UserObjects/openmc_mat{m.id}'
+ logger.info(f'Updating material {m.id} via {uo_path}')
+
+ names_path = f'{uo_path}/names'
+ self._control.setControllableVectorString(names_path, m.nuclides)
+
+ densities_path = f'{uo_path}/densities'
+ densities = list(m.densities.flat)
+ self._control.setControllableVectorReal(densities_path, densities)
+
+ break # hack for only sending the first
+
+ def _update_materials_and_nuclides(self, vec):
+ """Update the number density, material compositions, and nuclide
+ lists in helper objects
+
+ Parameters
+ ----------
+ vec : list of numpy.ndarray
+ Total atoms.
+
+ """
+ # Update the number densities regardless of the source rate
+ self.number.set_density(vec)
+ self._update_materials()
+
+ # Prevent OpenMC from complaining about re-creating tallies
+ openmc.reset_auto_ids()
+
+ # Update tally nuclides data in preparation for transport solve
+ nuclides = self._get_reaction_nuclides()
+ print(f"Nuclides: {nuclides}")
+ self._rate_helper.nuclides = nuclides
+ self._normalization_helper.nuclides = nuclides
+ self._yield_helper.update_tally_nuclides(nuclides)
+
+ def update_cardinal_tallies(self):
+ """Update the set of nuclides in the depletion tallies in the Cardinal
+ OpenMC model
+ """
+ logger.info('Waiting to update cardinal tally nuclides')
+
+ if self._control is None:
+ self.start_cardinal()
+
+ # Wait for cardinal to be at any flag
+ current_flag = self._control.wait()
+
+ # If we're at timestep end, keep going
+ if current_flag == 'TIMESTEP_END':
+ self._control.setContinue()
+
+ # Wait for cardinal to be ready at the beginning of a timestep
+ self._control.wait('TIMESTEP_BEGIN')
+
+ logger.info('Updating cardinal tally nuclides')
+ for t in openmc.lib.tally.tallies.values():
+ uo_path = f'UserObjects/openmc_tally{t.id}'
+ logger.info(f'Updating tally {t.id} via {uo_path}')
+ nuclides = [n for n in t.nuclides]
+ names_path = f'{uo_path}/nuclides'
+ self._control.setControllableVectorString(names_path, nuclides)
+
+ scores = [s for s in t.scores]
+ self._control.setControllableVectorString(f'{uo_path}/scores', scores)
+
+ filter_ids = [f.id for f in t.filters]
+ self._control.setControllableVectorString(f'{uo_path}/filter_ids', filter_ids)
+
+ break
+
+ def load_cardinal_results(self, sp_file):
+ """Load results from the statepoint file generated by the Cardinal
+ execution
+ """
+ with openmc.lib.quiet_dll(False):
+ openmc.lib.simulation_init()
+ openmc.lib.statepoint_load(sp_file)
+ openmc.lib.simulation_finalize()
+ logger.info('Statepoint loaded')
+
+ def initial_condition(self):
+ """Performs setup, stars the Cardinal problem (with a server enabled) and returns initial condition.
+
+ Parameters
+ ----------
+ materials : list of openmc.lib.Material
+ list of materials
+
+ Returns
+ -------
+ list of numpy.ndarray
+ Total density for initial conditions.
+
+ """
+
+ n = super().initial_condition()
+ self._update_materials_and_nuclides(n)
+
+ # Ensure that tally data is written to the statepoint file
+ for t in openmc.lib.tally.tallies.values():
+ t.writable = True
+
+ if self._control is None:
+ self.start_cardinal()
+
+ return n
+
+ def __call__(self, vec, source_rate):
+ """Runs a simulation.
+
+ Simulation will abort under the following circumstances:
+
+ 1) No energy is computed using OpenMC tallies.
+
+ Parameters
+ ----------
+ vec : list of numpy.ndarray
+ Total atoms to be used in function.
+ source_rate : float
+ Power in [W] or source rate in [neutron/sec]
+
+ Returns
+ -------
+ openmc.deplete.OperatorResult
+ Eigenvalue and reaction rates resulting from transport operator
+
+ """
+ # Reset results in OpenMC
+ openmc.lib.reset()
+
+ # start Cardinal if it hasn't been started already
+ if self._control is None:
+ self.start_cardinal()
+
+ self._update_materials_and_nuclides(vec)
+
+ # Update materials in Cardinal
+ self.update_cardinal_materials()
+
+ # If the source rate is zero, return zero reaction rates without running
+ # a transport solve
+ if source_rate == 0.0:
+ rates = self.reaction_rates.copy()
+ rates.fill(0.0)
+ return OperatorResult(ufloat(0.0, 0.0), rates)
+
+ # EXECUTE CARDINAL
+ self.update_cardinal_tallies()
+ sp_file = self.run_cardinal()
+
+ # LOAD RESULTS FROM STATEPOINT FILE GENERATED BY CARDINAL
+ self.load_cardinal_results(sp_file)
+
+ # Extract results
+ rates = self._calculate_reaction_rates(source_rate)
+
+ # Get k and uncertainty
+ keff = ufloat(*openmc.lib.keff())
+
+ op_result = OperatorResult(keff, rates)
+
+ return copy.deepcopy(op_result)
diff --git a/test/tests/depletion/chain_simple.xml b/test/tests/depletion/chain_simple.xml
new file mode 100644
index 000000000..c2e50a370
--- /dev/null
+++ b/test/tests/depletion/chain_simple.xml
@@ -0,0 +1,49 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 2.53000e-02
+
+ Gd157 Gd156 I135 Xe135 Xe136 Cs135
+ 1.093250e-04 2.087260e-04 2.780820e-02 6.759540e-03 2.392300e-02 4.356330e-05
+
+
+
+
+
+
+ 2.53000e-02
+
+ Gd157 Gd156 I135 Xe135 Xe136 Cs135
+ 6.142710e-5 1.483250e-04 0.0292737 0.002566345 0.0219242 4.9097e-6
+
+
+
+
+
+
+ 2.53000e-02
+
+ Gd157 Gd156 I135 Xe135 Xe136 Cs135
+ 4.141120e-04 7.605360e-04 0.0135457 0.00026864 0.0024432 3.7100E-07
+
+
+
+
diff --git a/test/tests/depletion/depletion.py b/test/tests/depletion/depletion.py
new file mode 100755
index 000000000..40d7ed018
--- /dev/null
+++ b/test/tests/depletion/depletion.py
@@ -0,0 +1,54 @@
+#!/usr/bin/env python3
+import os, shlex, sys
+from MooseControl import MooseControl
+
+# All of these imports are awful and we'll figure
+# out a better way to do it before merging
+try:
+ import openmc
+except:
+ openmc_contrib = os.path.abspath(os.path.join(os.path.dirname(__file__), '..', '..', '..', 'contrib', 'openmc'))
+ if os.path.exists(openmc_contrib):
+ sys.path.append(openmc_contrib)
+ import openmc
+
+
+import openmc.lib
+import openmc.deplete
+from openmc.deplete import CoupledOperator
+
+try:
+ from cardinal_operator import CardinalOperator
+except:
+ cardinal_contrib = os.path.abspath(os.path.join(os.path.dirname(__file__), '..', '..', '..', 'python'))
+ if os.path.exists(cardinal_contrib):
+ sys.path.append(cardinal_contrib)
+ from cardinal_operator import CardinalOperator
+
+if __name__ == '__main__':
+ # Get the command we should run
+ RUNAPP_COMMAND = os.environ.get('RUNAPP_COMMAND')
+ if RUNAPP_COMMAND is None:
+ sys.exit('Missing expected command variable RUNAPP_COMMAND')
+
+ # We need to automate the removal of this state because
+ # it'll definitely break things if we don't
+ files = ['tallies.out', 'tallies.xml', 'statepoint.40.h5', 'summary.h5']
+ for file in files:
+ if os.path.exists(file):
+ os.remove(file)
+
+ model = openmc.Model.from_xml()
+
+ operator = CardinalOperator(model, 'chain_simple.xml')
+ operator.cardinal_cmd = RUNAPP_COMMAND
+
+ timesteps = [0.1] * 2
+ timestep_units = 'd'
+ power = 174 # Wth (pincell)
+
+ integrator = openmc.deplete.PredictorIntegrator(operator,
+ timesteps,
+ power,
+ timestep_units=timestep_units)
+ integrator.integrate()
diff --git a/test/tests/depletion/geometry.xml b/test/tests/depletion/geometry.xml
new file mode 100644
index 000000000..bf6e32ea6
--- /dev/null
+++ b/test/tests/depletion/geometry.xml
@@ -0,0 +1,14 @@
+
+
+ |
+ |
+ |
+ |
+
+
+
+
+
+
+
+
diff --git a/test/tests/depletion/gold/openmc_out.e b/test/tests/depletion/gold/openmc_out.e
new file mode 100644
index 000000000..7dd732e78
Binary files /dev/null and b/test/tests/depletion/gold/openmc_out.e differ
diff --git a/test/tests/depletion/materials.xml b/test/tests/depletion/materials.xml
new file mode 100644
index 000000000..90f605eeb
--- /dev/null
+++ b/test/tests/depletion/materials.xml
@@ -0,0 +1,53 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/test/tests/depletion/openmc.i b/test/tests/depletion/openmc.i
new file mode 100644
index 000000000..c024bc1b4
--- /dev/null
+++ b/test/tests/depletion/openmc.i
@@ -0,0 +1,69 @@
+[Mesh]
+ [file]
+ type = FileMeshGenerator
+ file = gold/openmc_out.e
+ []
+[]
+
+[AuxVariables]
+ [cell_temperature]
+ family = MONOMIAL
+ order = CONSTANT
+ []
+[]
+
+[AuxKernels]
+ [cell_temperature]
+ type = CellTemperatureAux
+ variable = cell_temperature
+ []
+[]
+
+[Problem]
+ type = OpenMCCellAverageProblem
+ verbose = true
+ power = ${fparse 3000e6 / 273 / (17 * 17)}
+ cell_level = 0
+ temperature_blocks = '1 2 3'
+
+
+ [Tallies]
+ [heat_source]
+ type = CellTally
+ name = heat_source
+ block = '2 3'
+ []
+ []
+
+ volume_calculation = vol
+[]
+
+[UserObjects]
+ [vol]
+ type = OpenMCVolumeCalculation
+ n_samples = 100000
+ []
+[]
+
+[Executioner]
+ type = Transient
+ num_steps = 3
+[]
+
+[Outputs]
+ exodus = true
+[]
+
+[Postprocessors]
+ [heat_source]
+ type = ElementIntegralVariablePostprocessor
+ variable = heat_source
+ []
+ [max_tally_rel_err]
+ type = TallyRelativeError
+ []
+ [max_heat_source]
+ type = ElementExtremeValue
+ variable = heat_source
+ []
+[]
diff --git a/test/tests/depletion/settings.xml b/test/tests/depletion/settings.xml
new file mode 100644
index 000000000..5c6d2e045
--- /dev/null
+++ b/test/tests/depletion/settings.xml
@@ -0,0 +1,13 @@
+
+
+ eigenvalue
+ 1000
+ 5
+ 2
+
+
+ -0.62992 -0.62992 -1.0 0.62992 0.62992 1.0
+
+
+ 100
+
diff --git a/test/tests/depletion/tests b/test/tests/depletion/tests
new file mode 100644
index 000000000..9bd258cdc
--- /dev/null
+++ b/test/tests/depletion/tests
@@ -0,0 +1,9 @@
+[Tests]
+ [test]
+ type = Exodiff
+ input = openmc.i
+ exodiff = 'openmc_out.e'
+ command_proxy = depletion.py
+ required_objects = 'OpenMCCellAverageProblem'
+ []
+[]