continue_timesteps producing h5 files with incomplete information when used to restart a simulation that did not specify final_step = True in Integrator

[lewisgross1296]

Bug Description

I'm running a model that needs to do one depletion step at a time, as I'm loading new settings.properties_files and changing a rotation in my model to follow a specific control drum history on each step. Since I'm running one depletion step at a time, I'm aware of two options

1. Overwrite the XML and depletion_results.h5 file at each step (and perhaps save the old one for postprocessing)
2. Use the continue_timesteps feature with one new additional timestep for each depletion run to produce on depletion_results.h5 file that should have all the data I want in one file.

I decided to go option 2 to minimize the number of output files and File I/O required during my run. Interestingly, my script "runs" as is and the output stream seems to indicate that transport/depletion is happening. However, when I examine the h5 output files, they do not contain all the information I'd hoped to see. For example,

>>> import openmc.deplete
>>> r = openmc.deplete.Results("depletion_results.h5")
>>> t,k = r.get_keff()
>>> k
array([[1.00072744e+00, 3.34656867e-04],
       [0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00],
... more zeros ...
       [0.00000000e+00, 0.00000000e+00]])
>>> s0 = openmc.StatePoint("openmc_simulation_n0.h5")
>>> s0.keff
1.0007274367827947+/-0.00033465686722658917
>>> s1 = openmc.StatePoint("openmc_simulation_n1.h5")
>>> s1.keff
0.9999761313651004+/-inf

The first StatePoint seems to be okay but the subsequent ones are missing sigma data. There must be some interaction between h5 writing and using continue_timesteps that is causing this behavior. I don't believe the actual transport or nuclide updates are producing wrong info. Feels more like an hdf5 write bug.

Steps to Reproduce

My model with the above results is a bit complicated, so I slightly modified the pincell_depletion example to use the same looping structure that my code ran with. Interestingly, the same issues were present with the pincell_depletion case, which indicates a bug, as opposed to a model error or dependency error, to me.

from math import pi

import openmc
import openmc.deplete
import os

###############################################################################
#                              Define materials
###############################################################################

# Instantiate some Materials and register the appropriate Nuclides
uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment')
uo2.set_density('g/cm3', 10.29769)
uo2.add_element('U', 1., enrichment=2.4)
uo2.add_element('O', 2.)

helium = openmc.Material(name='Helium for gap')
helium.set_density('g/cm3', 0.001598)
helium.add_element('He', 2.4044e-4)

zircaloy = openmc.Material(name='Zircaloy 4')
zircaloy.set_density('g/cm3', 6.55)
zircaloy.add_element('Sn', 0.014, 'wo')
zircaloy.add_element('Fe', 0.00165, 'wo')
zircaloy.add_element('Cr', 0.001, 'wo')
zircaloy.add_element('Zr', 0.98335, 'wo')

borated_water = openmc.Material(name='Borated water')
borated_water.set_density('g/cm3', 0.740582)
borated_water.add_element('B', 4.0e-5)
borated_water.add_element('H', 5.0e-2)
borated_water.add_element('O', 2.4e-2)
borated_water.add_s_alpha_beta('c_H_in_H2O')

###############################################################################
#                             Create geometry
###############################################################################

# Define surfaces
pitch = 1.25984
fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR')
clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
box = openmc.model.RectangularPrism(pitch, pitch, boundary_type='reflective')

# Define cells
fuel = openmc.Cell(fill=uo2, region=-fuel_or)
gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
water = openmc.Cell(fill=borated_water, region=+clad_or & -box)

# Define overall geometry
geometry = openmc.Geometry([fuel, gap, clad, water])

###############################################################################
#                     Set volumes of depletable materials
###############################################################################

# Set material volume for depletion. For 2D simulations, this should be an area.
uo2.volume = pi * fuel_or.r**2

###############################################################################
#                     Transport calculation settings
###############################################################################

# Instantiate a Settings object, set all runtime parameters, and export to XML
settings = openmc.Settings()
settings.batches = 100
settings.inactive = 10
settings.particles = 1000

# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
settings.source = openmc.IndependentSource(
    space=uniform_dist, constraints={'fissionable': True})

###############################################################################
#                   Initialize and run depletion calculation
###############################################################################

model = openmc.Model(geometry=geometry, settings=settings)

# Perform simulation using the predictor algorithm
time_steps = [1.0, 1.0, 1.0, 1.0, 1.0]  # days
power = 174  # W/cm, for 2D simulations only (use W for 3D)

t_in_prog = []
for dt in time_steps:
    t_in_prog.append(dt)
    if len(t_in_prog) == 1:
        ct = False
        pr = None
    else:
        ct = True
        if os.path.isfile("depletion_results.h5"):
            pr = openmc.deplete.Results("depletion_results.h5")
    powers = [power] * len(t_in_prog)
    op = openmc.deplete.CoupledOperator(model, prev_results=pr, chain_file='chain_simple.xml')
    openmc.deplete.CECMIntegrator(
        op, t_in_prog, powers, timestep_units="d", continue_timesteps=ct).integrate(final_step=False)

###############################################################################
#                    Read depletion calculation results
###############################################################################

# Open results file
results = openmc.deplete.Results("depletion_results.h5")

# Obtain K_eff as a function of time
time, keff = results.get_keff(time_units='d')

print("testing depletion_results.h5")
print(time)
print(keff)

steps = [0,1,2,3,4,5]
for s in steps:
    sp = openmc.StatePoint(f"openmc_simulation_n{s}.h5")
    print(sp.keff)

Environment

Using this commit

     Commit Hash | 111eb770668d2f8e39d7fd4524a08df8bf9819f5

Tested both locally and on a cluster, so pretty sure it's not a dependency or environment issue.

[lewisgross1296]

A few observations, when the above MWE is rerun with only changing final_step=False -> final_step=True, all the data appears to be present.

We discussed today that it's actually implicitly required that if you want to use continue_timesteps = True, then you would need to have run a final transport solve on the last step. This is because the continue_timesteps option starts the next depletion step assuming that data for the previous step is complete.

@pshriwise suggested that we can enforce this by adding a new variable to the Results class that indicates whether final_step was True or False. If we attempt to do a continue_timesteps run (in this logic block) and final_step was originally False in the last simulation, then we can have OpenMC run the extra transport it would need for continuity of depletion.

For now, everything should work as long as final_step=True is used with continue runs. I'd like to get a PR going to prevent users from using these options incorrectly, but it might be a little before I get a chance.

Hopefully, I'll get a chance before v0.15.4 release