Depletion

A transport calculation answers a question about one instant. Depletion asks what happens over months: U-235 is consumed, plutonium grows in from U-238, fission products accumulate and some of them absorb neutrons voraciously, and the reactivity of the fuel drifts down until the cycle ends.

The governing equations are the Bateman equations — one ordinary differential equation per nuclide, saying that its concentration changes by production from decay and transmutation of others, minus its own decay and absorption. What makes them hard is not the arithmetic but the coefficients: every transmutation term needs a reaction rate, which is a cross section folded with the neutron flux, and the flux depends on the composition that is being solved for. The system is coupled and nonlinear.

OpenMC splits the work in two, and the vocabulary follows the split. An operator runs transport and reports one-group reaction rates for the current composition. An integrator advances the composition over a timestep using those rates, deciding how many transport solves to spend and how to combine them. You build both, hand the operator to the integrator, and call integrate().

A working transport model is not yet a depletion model. Three additions are required, and two of them produce errors that are much easier to understand once you know what they are asking for.

The first is marking materials as depletable. Only materials with depletable = True are tracked, which is what keeps a moderator or a steel structure from acquiring a fission product inventory it has no business having.

The second is a volume, and this one surprises people. Transport never needs to know how large a material is — it works in atom densities, and the geometry supplies the shape. The Bateman equations are written in absolute numbers of atoms, so depletion needs to convert a density into an inventory, and that conversion is the volume. Without it OpenMC cannot proceed and will tell you so. The value is in cubic centimeters, and for anything harder than a cylinder, openmc.VolumeCalculation will compute it stochastically.

import openmc

fuel = openmc.Material(name='UO2')
fuel.set_density('g/cm3', 10.4)
fuel.add_nuclide('U235', 0.045)
fuel.add_nuclide('U238', 0.955)
fuel.add_element('O', 2.0)

fuel.depletable = True
fuel.volume = 3.14159 * 0.4096**2 * 366.0   # cm3 — required, not optional

The third is a depletion chain, which is a separate data file from the cross sections. It holds decay constants, branching ratios, and fission product yields — everything about how nuclides turn into other nuclides, none of which lives in a transport library. Chains are distributed from openmc.org, and they come in sizes: a full chain carries the better part of two thousand nuclides, while the CASL chain is reduced to a few hundred chosen for light-water reactor work and runs considerably faster.

# Once, globally — the operator picks this up by default
openmc.config['chain_file'] = '/path/to/chain_casl_pwr.xml'

# Or explicitly, per operator
operator = openmc.deplete.CoupledOperator(model, chain_file='chain_casl_pwr.xml')

A material in OpenMC is a composition, and a single material object can fill hundreds of cells. That is exactly what you want for transport, where all 264 fuel pins of an assembly are identical and sharing one object saves memory. In depletion it produces a result that is quietly wrong.

Depletion works on materials, not on cells. If one material fills every pin, there is one inventory for all of them, depleting at the flux-weighted average rate. The pins at the assembly edge, which see a softer spectrum and burn differently from those at the center, are averaged away — and so is the whole point of computing a burnup distribution.

Differentiating them creates a separate material per instance. The volume of the original is divided among the new materials, which is why the volume you set has to be the volume of one pin's worth of fuel or the total, consistently with what you asked for.

# Split each multiply-instantiated depletable material into one material
# per cell instance. The original volume is divided equally by default.
operator = openmc.deplete.CoupledOperator(
    model,
    diff_burnable_mats=True,
    diff_volume_method='divide equally',   # or 'match cell'
)

The cost is real: 264 depletion zones means 264 Bateman solves per step and 264 sets of reaction rate tallies, and a full-core model differentiated pin by pin is a serious calculation. Differentiate what you need resolved and leave the rest lumped.

Reaction rates from a transport solve are per source particle. Turning them into rates per second requires knowing how hard the reactor is actually running, which is what the power argument is for. Give it as power in watts, or as power_density in watts per gram of heavy metal — the second is often more natural, since 30 to 40 W/gHM is a recognizable pressurized water reactor and a total wattage for a single pin is not.

# Constant power, six timesteps of 30 days
integrator = openmc.deplete.CECMIntegrator(
    operator,
    timesteps=[30.0] * 6,
    power=None,
    power_density=35.0,      # W per gram of heavy metal
    timestep_units='d',
)
integrator.integrate()

Power may also vary. Passing a list the same length as the timesteps gives a power history, and a zero-power step is how you model a shutdown — decay continues, transmutation stops, and that is how you follow the Xe-135 peak after a trip or cool spent fuel in a pond.

The normalization_mode setting decides how the energy release is computed. The default fission-q uses fission Q values from the chain, attributing the full recoverable energy to the fission event. energy-deposition instead tallies where energy is actually deposited, which is more faithful and requires the transport model to carry the photons that make it more faithful. For a fixed-source problem — an accelerator target, a fusion blanket — there is no fission power to normalize to, and source-rate takes neutrons per second instead.

Timesteps are where depletion calculations are won or lost, and the right length is set by the fastest thing you care about rather than by tidy divisions of the cycle.

Xe-135 is the usual constraint at startup. It is the strongest neutron absorber in the core, it is produced mostly through the decay of I-135 rather than directly, and the pair reach equilibrium over roughly forty hours. A first timestep of thirty days steps straight over that transient, and the reactivity history you get will be missing the several-thousand-pcm dip that a real core goes through in its first two days. Begin with steps measured in hours, lengthen through days, and only then move to the thirty-day steps that carry the cycle.

Burnable absorbers impose the same discipline for a different reason. Gadolinium and boron are consumed quickly and shield themselves while they last, so their absorption changes sharply over the first few thousand megawatt-days. A step long enough to average across that change will misplace the reactivity peak.

# A schedule shaped by the physics rather than by the calendar
timesteps = (
    [0.25, 0.5, 1.0, 2.0, 4.0]     # days: catch the Xe-135 transient
    + [10.0] * 3                   # early gadolinium burnout
    + [30.0] * 15                  # the rest of the cycle
)

integrator = openmc.deplete.CECMIntegrator(
    operator, timesteps, power_density=35.0, timestep_units='d'
)

Whether a schedule is fine enough is an empirical question with a straightforward answer: halve the steps and see whether the answer moves. If it does, the original schedule was too coarse, and no amount of statistical precision within each step will fix that.

Every integrator solves the same equations; they differ in how many transport solves they spend per timestep and what they do with them. Since the transport solve dominates the cost entirely, that count is the price list.

PredictorIntegrator spends one. It takes the reaction rates at the start of the step and holds them constant across it, which is first-order accurate and systematically wrong in a knowable direction — rates are evaluated on fuel that has not yet burned. CECMIntegrator spends two, evaluating rates at the start and again at the midpoint, and is the sensible default: the second solve buys most of the accuracy available. CELIIntegrator and LEQIIntegrator also spend two, with different interpolation schemes, and LEQI reuses information from the previous step, which helps when step lengths vary. CF4Integrator and EPCRK4Integrator spend four for fourth-order accuracy, worth it when a reference solution is the goal. The stochastic implicit integrators, SICELIIntegrator and SILEQIIntegrator, iterate to self-consistency and cost the most.

The choice interacts with step length. A better integrator with long steps and a simple integrator with short steps can reach similar accuracy for similar total cost, and the second is easier to debug. Start with CECMIntegrator and a schedule you have tested by halving.

The run writes depletion_results.h5, containing k and the full inventory of every depletable material at every step. Note that concentrations come back in whichever units you ask for — atoms by default, which is a raw count and not what you want on a plot against a density axis.

import matplotlib.pyplot as plt
import openmc.deplete

results = openmc.deplete.Results('depletion_results.h5')

# k over the cycle. Column 0 is the value, column 1 its standard deviation.
time, keff = results.get_keff(time_units='d')

# Inventories. nuc_units defaults to 'atoms' — ask for a density explicitly.
_, u235 = results.get_atoms('1', 'U235', nuc_units='atom/b-cm', time_units='d')
_, pu239 = results.get_atoms('1', 'Pu239', nuc_units='atom/b-cm', time_units='d')
_, xe135 = results.get_atoms('1', 'Xe135', nuc_units='atom/b-cm', time_units='d')

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))

ax1.errorbar(time, keff[:, 0], yerr=keff[:, 1], marker='o')
ax1.axhline(1.0, color='k', ls=':')
ax1.set(xlabel='time (days)', ylabel='k-effective')

ax2.semilogy(time, u235, label='U-235')
ax2.semilogy(time, pu239, label='Pu-239')
ax2.semilogy(time, xe135, label='Xe-135')
ax2.set(xlabel='time (days)', ylabel='atom/b-cm')
ax2.legend()

fig.tight_layout()
fig.savefig('depletion.png', dpi=150)

Three features tell you the calculation is behaving. U-235 falls almost linearly, because consumption is proportional to flux and the flux is held roughly constant by the power normalization. Pu-239 rises steeply and then flattens as its own destruction catches up with its production from U-238 — a fresh pressurized water reactor pin ends its life getting a third or more of its fissions from plutonium that did not exist at loading. Xe-135 jumps within the first two days and then sits flat, which is the equilibrium the fine early steps were there to resolve.

Burnup is conventionally reported in megawatt-days per kilogram of initial heavy metal, and the operator has already computed the denominator: operator.heavy_metal is the initial heavy metal inventory in grams. Multiplying power by elapsed time and dividing gives the axis that lets you compare against published depletion benchmarks, where a discharged pressurized water reactor assembly is somewhere around 45 to 55 MWd/kgU.

A long depletion is worth restarting rather than repeating. Passing the previous results to a new operator through prev_results continues from the last state, which is what you want when a cycle finishes and you would like to extend it, or when a job hits a wall clock limit. For plotting the flux and reaction rate distributions at a particular burnup, the statepoint files written along the way behave exactly like those from a standalone run.