Hi! I am having trouble with depletion calculations. I created a typical cubic reactor and I’m trying to perform depletion calculations. Below I have attached some of the important settings and geometry I’ve done; we are running 120 batches. The fuel I’m running the simulations for contains a great percentage Uranium-235 but is contaminated by other myriad of isotopes all over the periodic table. My main issue is that when I first run OpenMC I get a supercritical value (~1.10), but when I perform depletion calculations I get two files with subcritical values (0.62, 0.53: i.e., the numbers are not consistent). And even the graph below shows that our initial keff is past 1.6, not even being consistent with the initial OpenMC run. I have had a graduate student mentor look at it and we cannot pinpoint the issue. Some questions we have asked ourselves are: why is the initial keff changing when we change the volume? Why is the steady-state keff different from the initial burnup core keff, which should be identical? Manipulating our inputs, why is the keff only consistent when the fuel comprises only three nuclides (O16, U235, U238)?

We might think there is an issue with how Jupyter performs the depletion calculations, but we’re not sure. Any help is greatly appreciated and your luck will triple!! We’ve been dealing with this issue for a while with no clear solution. :'^(

The depletion looks like (I would’ve posted a screenshot but this post only allowed one image):

```
#DEPLETION
import openmc.deplete
model = openmc.Model(geometry=geometry, materials=materials, settings=settings)
operator = openmc.deplete.CoupledOperator(model, "/Users/kdoll/openmc/lib/chain_endfb71_pwr.xml")
#We will take depletion step sizes of 30days
# and instruct OpenMC to re-run a transport simulation every 30 days until we have modeled the problem over a six month cycle.
#The depletion interface expects the time to be given in seconds, so we will have to convert.
total_power = 3400*10**3 #MWthermal
time_steps = [300] * 1
integrator = openmc.deplete.PredictorIntegrator(operator, time_steps, total_power, timestep_units='d')
integrator.integrate()
```