How to include burnup in coupled calculations?

[LeonDong1998]

Hello everyone!
I have encountered a difficult problem. I want to implement the coupled calculation of MOOSE and openmc burnup using the same time scale. But it seems that burnup calculation can only be done through integrator calculation, and burnup settings cannot be generated into XML files.
I'd like to express my gratitude in advance!

[aprilnovak]

Hi @LeonDong1998, we are very close to having this natively supported in Cardinal (you can find details here: #921). In this PR, you can find an in-memory approach for iteratively running (i) Cardinal multiphysics statepoint to predict temperatures and densities and (ii) OpenMC depletion calculations to update composition. Our team is still currently reviewing the capability and hope to have this merged within the next year but at the moment there are a few other developments higher on the task list.

We would greatly welcome you using this PR for your needs - the functionality should be present, but we do not yet have any documentation.

There is an alternative (equivalent) approach where you can iteratively run Cardinal and OpenMC depletion calculations, but using a file-based approach. With this alternative option, you can (i) run Cardinal as you normally do, outputting a properties.h5 file which contains temperature and density; (ii) loading this file into an OpenMC depletion script that runs just 1 depletion step and outputs a statepoint file with the new composition. You can then repeat these two steps until the end of the depletion cycle, loading data from files. I've attached an example script which shows how to use this.

import sys, os

# This decreases the number of warnings that get dumped to console during the build process.
import warnings
warnings.filterwarnings("ignore")

# You'll need lib and deplete at a minimum if you want to do "in memory" multiphysics
import openmc
import openmc.lib
import openmc.deplete

# I use these for post-processing. Feel free to remove them
import pandas as pd
import numpy as np

# Function to build the boronated model.
def build_boronated_model(boron_ppm):
  # This script pulls in a model from an xml file, so we need to reset the ID mixins whenever the new model is built to
  # prevent them from getting scrambled (multiphysics data is mapped by id + instance)
  openmc.Material.reset_ids()
  openmc.Cell.reset_ids()
  openmc.Surface.reset_ids()
  openmc.Universe.reset_ids()
  openmc.PlotBase.reset_ids()

  # Import the model from the model.xml file used by Cardinal to ensure that IDs map to the multiphysics data
  m = openmc.Model.from_model_xml('../quarter_core_simplified_mp/model.xml')

  # Swap the materials out for the depleted materials (if they exist).
  if os.path.isfile('./depletion_results.h5'):
    print('Loading depleted materials from the previous step.')
    result_mats = openmc.deplete.Results("./depletion_results.h5").export_to_materials(burnup_index=1)
    m.materials = result_mats

  # Modify materials. This is necessary for BEAVRS - modify as required for your own use case!
  for i in range(len(m.materials)):
    # Set bornable absorbers to depletable.
    if m.materials[i].id == 15:
      m.materials[i].depletable = True

    # Change the boron ppm in the coolant.
    if m.materials[i].id == 16:
      mat_name = m.materials[i].name
      m.materials[i] = openmc.model.borated_water(boron_ppm, temperature=550, temp_unit='F', pressure=2250, press_unit='psi')
      m.materials[i].id = 16
      m.materials[i].name = mat_name

  # Some more practical settings for prototyping.
  m.settings.inactive = 200
  m.settings.batches = 1200
  m.settings.particles = 10000

  # Initialize OpenMC lib and use it to load temperatures and densities. Running model.import_properties is
  # too slow for very large geometries with many distribcells.
  # TODO: can we get away with doing this once somehow?
  print('Loading cell densities / temperatures.')
  with openmc.lib.quiet_dll(False):
    m.init_lib()
    openmc.lib.import_properties('../quarter_core_simplified_mp/properties.h5')

  return m

# Function to search for k_eff and return the critical boron ppm.
def criticality_search():
  crit_ppm, guesses, keffs = openmc.search_for_keff(build_boronated_model, bracket=[0., 1000.],
                                                    tol=1e-2, print_iterations=True)
  return crit_ppm

# Function to run a single depletion step.
def single_depletion_step(dt, power, previous_model):
  print('Depleting model from the previous timestep.')
  # We really should be running with distribmats, but auto distribmats don't appear to work for the BEAVRS model (to big?)
  # previous_model.differentiate_depletable_mats(diff_volume_method='divide equally')
  op = openmc.deplete.CoupledOperator(previous_model, chain_file='./chain_casl_pwr.xml')
  openmc.deplete.CECMIntegrator(op, [dt], [power], timestep_units='d').integrate(final_step=False)

  print('Searching for critical boron.')
  return criticality_search()

def main():
  # Settings for depletion. Change these as required
  # openmc.deplete.pool.USE_MULTIPROCESSING = False
  openmc.deplete.pool.NUM_PROCESSES = 100

  # Load the power history and timesteps.
  history = pd.read_csv('./simplified_power_history.csv')
  # Convert percent full power to quarter-core power.
  power   = 0.25 * 3411e6 * history['power_history'].to_numpy() / 100.0
  delta_t = history['dt'].to_numpy()

  # Compute the initial condition with a criticality search.
  boron_letdown = [criticality_search()]
  depleted_model = build_boronated_model(boron_letdown[0])

  # Run depletion with a criticality search one step at a time.
  for i in range(len(power)):
    # Save the boron letdown curve to disk as we're running.
    np.savetxt('./boron_letdown_curve.csv', np.array(boron_letdown))
    # Depletion followed by a criticality search.
    boron_letdown.append(single_depletion_step(delta_t[i], power[i], depleted_model))
    # Regenerate the model using the given boron concentration
    # to use for the next depletion step.
    depleted_model = build_boronated_model(boron_letdown[-1])

  print('Boron letdown curve:')
  print(np.array(boron_letdown))
  np.savetxt('./boron_letdown_curve.csv', np.array(boron_letdown))
  openmc.lib.finalize()

if __name__ == "__main__":
    main()