restarting killed depletion simulations improperly saving h5 data

[lewisgross1296]

Bug Description

SEE TL;DR at end for concise summary

In #3272, I added the capability for depletion runs to continue an existing run, given that the initial timesteps/source rates found in depletion_results.h5 match those provided to the integrator requesting a continue run.

I found that the continue feature runs "successfully" once but that a second attempt at continuing the same simulation causes the error handling to think you've improperly defined a continue run, i.e. hitting this block.

openmc/openmc/deplete/abc.py

Lines 666 to 671 in 820648d

	if not np.array_equal(prev_timesteps, timesteps[:num_prev]):
	raise ValueError(
	"You are attempting to continue a run in which the previous timesteps "
	"do not have the same initial timesteps as those provided to the "
	"Integrator. Please make sure you are using the correct timesteps."
	)

This happened despite using identical data in the first, second and third run... perplexing. I was getting this issue every time I attempted a second continue.

I think I've narrowed down the issue to the writing of the depletion_results.h5 file. I modified the Integrator class source code locally to print out what the class thought all the data was and discovered that one of the timesteps in the depletion_results.h5 did not match the actual dt supplied at that step.

For example, here are the timesteps passed to the integrator in my depletion run script

This simulation will deplete at 225000.0 W with the following time steps (units of days)
[1.0, 1.0, 1.0, 1.0, 1.0, 5.0, 5.0, 5.0, 15.0, 15.0, 15.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0]

When attempting the second continue run, I printed some timestep checks inside of openmc/deplete/abc.py Here's the printing of prev_timesteps declared here

openmc/openmc/deplete/abc.py

Lines 658 to 661 in 820648d

	# Get timesteps and source rates from previous results
	prev_times = operator.prev_res.get_times(timestep_units)
	prev_source_rates = operator.prev_res.get_source_rates()
	prev_timesteps = np.diff(prev_times)

and timesteps[:num_prev] from the snippet above with the if condition

prev_timesteps = 
 [  1.   1.   1.   1.   1.   5.   5.   5.  15.  15.  15.  60.  60. 120.  60.  60.  60.  60.  60.  60.  60.  60.  60.  60.  60.  60.  60.  60.]
timesteps[:num_prev] = 
 [1.0, 1.0, 1.0, 1.0, 1.0, 5.0, 5.0, 5.0, 15.0, 15.0, 15.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0, 60.0]

Notice that there is a 120 in the timesteps from depletion_results.h5 where there is a 60 in the timesteps provided to the integrator.

Looking at the HPC, I noticed coincidentally that the 120 corresponds exactly to the last step of the first run. This makes me suspect that the first continue run mistakenly increments the time at that step in the first write to h5 of the continue.

Steps to Reproduce

The fastest way to reproduce would be to launch a depletion simulation with two or more steps. Then, launch a continue run (with the same first steps) and two additional steps. If you launch another properly defined continue run and any number of new steps, the error message for attempting a continue run with invalid steps should occur. I can add an MWE for this soon

Environment

Initially found on an HPC environment using OpenMC v0.15.2. I recreated locally with this commit

[lewisgross1296]

I've made an MWE (see zip) for debugging purposes that corresponds with this branch on my fork. The commit on that branch is only python changes, so it should be easy to grab the changed openmc/deplete/abc.py and use it locally (no recompile required).

Right now, the MWE produces this output at the failure point

timestep checks
 num_prev = 4
 prev_timesteps = 
 [0. 2. 1. 1.]
 initial_timesteps_provided_to_integrator = 
 [1, 1, 1, 1]
Traceback (most recent call last):
  File "/home/lgross/restart_failure/mwe/deplete.py", line 54, in <module>
    integrator = openmc.deplete.CECMIntegrator(
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/lgross/openmc/openmc/deplete/abc.py", line 673, in __init__
    raise ValueError(
ValueError: You are attempting to continue a run in which the previous timesteps do not have the same initial timesteps as those provided to the Integrator. Please make sure you are using the correct timesteps.

mwe.zip
Steps to reproduce:

1) python build_model.py
2) python deplete.py

I can attach the captured output from the depletion runs too if that would be useful

[lewisgross1296]

So as a sanity check, I tried to duplicate the simulation without continue_timesteps = True. It looks like any restarted depletion simulation corrupts the h5 results (EDIT: more accurately would be to say it only writes the current steps to h5 and doesn't copy the data from before the restart and instead writes 0s). I used the same XML as above but restarted by just supplying time_steps = [1] * 2 three total times. Here's what prints from each of the h5 files produced from depletion

TESTING EACH H5 FLIE
initial results get_times()
[0. 1. 2.]


first restart get_times()
[0. 0. 2. 3. 4.]


second restart get_times()
[0. 0. 0. 0. 4. 5. 6.]

What I would expect is

initial results get_times()
[0. 1. 2.]

first restart get_times()
[0. 1. 2. 3. 4.]

second restart get_times()
[0. 1. 2. 3. 4. 5. 6.]

I think something deeper is wrong with the h5 writing and restarted simulations and isn't actually caused by #3272 since the symptom appears in the absence of continue runs.

Here's the file in case anyone else wants to try

import openmc
import openmc.deplete

# load model
model = openmc.Model.from_xml()
power = 225000 # W

# first depletion solve
time_steps = [1] * 2
operator = openmc.deplete.CoupledOperator( model , chain_file="casl_chain.xml")
integrator = openmc.deplete.CECMIntegrator(
    operator,
    time_steps,
    power=power,
    timestep_units="d",
)
print("BEFORE FIRST SOLVE")
integrator.integrate(final_step=True, output=True)
print("\n")
print("\n")
print("\n")
print("FIRST DEPLETION SOLVE COMPLETE")

# add two more 1 day timestpes and rerun
time_steps = [1] * 2
initial_results = (
    openmc.deplete.Results("depletion_results.h5")
)
operator = openmc.deplete.CoupledOperator(
    model, prev_results=initial_results, chain_file="casl_chain.xml"
)
integrator = openmc.deplete.CECMIntegrator(
    operator,
    time_steps,
    power=power,
    timestep_units="d")
print("\n")
print("\n")
print("\n")
print("BEFORE SECOND SOLVE")
integrator.integrate(final_step=True, output=True, path="first_restart.h5")
print("SECOND DEPLETION SOLVE COMPLETE")

# add two more 1 day timestpes and rerun
time_steps = [1] * 2
first_restart = (
    openmc.deplete.Results("first_restart.h5")
)
operator = openmc.deplete.CoupledOperator(
    model, prev_results=first_restart, chain_file="casl_chain.xml"
)
integrator = openmc.deplete.CECMIntegrator(
    operator,
    time_steps,
    power=power,
    timestep_units="d",
)
print("\n")
print("\n")
print("\n")
print("BEFORE THIRD SOLVE")
integrator.integrate(final_step=True, output=True, path="second_restart.h5")
print("THIRD DEPLETION SOLVE COMPLETE")



final_results = (
    openmc.deplete.Results("second_restart.h5")
)


print("TESTING EACH H5 FLIE")
print("initial results get_times()")
print(initial_results.get_times())
print("\n")
print("first restart get_times()")
print(first_restart.get_times())
print("\n")
print("second restart get_times()")
print(final_results.get_times())

[lewisgross1296]

Still trying to figure out the reason for these observed behavior(s).

While using different names for the depletion files results in 0s for the points simulated by earlier runs, I've noticed when a depletion simulation re-uses the same depletion_results.h5 file, all the time steps appear as expected.

I recreated the above test with a single depletion_results.h5 file and the continue_timesteps = True feature worked. If other's want to test, this file runs with no errors

import openmc
import openmc.deplete

# load model
model = openmc.Model.from_xml()
power = 225000 # W

# first depletion solve
time_steps = [1] * 2
operator = openmc.deplete.CoupledOperator( model , chain_file="casl_chain.xml")
integrator = openmc.deplete.CECMIntegrator(
    operator,
    time_steps,
    power=power,
    timestep_units="d",
)
integrator.integrate(final_step=True, output=True)
print("\n")
print("\n")
print("\n")
print("FIRST DEPLETION SOLVE COMPLETE")

# add two more 1 day timestpes and rerun
time_steps = [1] * 4
prev_results = (
    openmc.deplete.Results("depletion_results.h5")
)
print("initial results get_times()")
print(prev_results.get_times())
print(prev_results.get_keff())
operator = openmc.deplete.CoupledOperator(
    model, prev_results=prev_results, chain_file="casl_chain.xml"
)
integrator = openmc.deplete.CECMIntegrator(
    operator,
    time_steps,
    power=power,
    timestep_units="d",
    continue_timesteps=True
)
print("\n")
print("\n")
print("\n")
print("BEFORE SECOND SOLVE")
integrator.integrate(final_step=True, output=True)
print("SECOND DEPLETION SOLVE COMPLETE")

# add two more 1 day timestpes and rerun
time_steps = [1] * 6
prev_results = (
    openmc.deplete.Results("depletion_results.h5")
)
print("first continue results get_times(), get_keff")
print(prev_results.get_times())
print(prev_results.get_keff())
operator = openmc.deplete.CoupledOperator(
    model, prev_results=prev_results, chain_file="casl_chain.xml"
)
integrator = openmc.deplete.CECMIntegrator(
    operator,
    time_steps,
    power=power,
    timestep_units="d",
    continue_timesteps=True
)
print("\n")
print("\n")
print("\n")
print("BEFORE THIRD SOLVE")
integrator.integrate(final_step=True, output=True)
print("THIRD DEPLETION SOLVE COMPLETE")
final_results = openmc.deplete.Results("depletion_results.h5")
print("second continue results get_times(), get_keff")
print(prev_results.get_times())
print(prev_results.get_keff())

and prints

THIRD DEPLETION SOLVE COMPLETE
second continue results get_times(), get_keff
[0. 1. 2. 3. 4.]
(array([     0.,  86400., 172800., 259200., 345600.]), array([[1.26682373, 0.00566476],
       [1.2076182 , 0.00432229],
       [1.15074937, 0.00514333],
       [1.08728389, 0.00420871],
       [1.03888102, 0.00545256]]))

At this point, I'm not really sure how to keep investigating. There's no "error" being thrown at the point in the simulation that is going wrong (i.e. the point at which data is written or read improperly), so if anyone has ideas for how to keep debugging this, I would appreciate them.

[gonuke]

@lewisgross1296 and I may have found a solution - as little as 1 character change (plus some comments) - stay tuned.

[lewisgross1296]

TL;DR

There are two behaviors above, but one is not responsible for any of the observed issues. Basically, writing restarts to different filenames creates restarts that do not contain the data from their prev_results (they are stored as zeros). I will make a second issue for this, but it's a low priority task, since you can get around it by only using depletion_results.h5 as your filename. (Perhaps a good task for an intern to learn more about OpenMC).

The main issue has to do with which time a restarted simulation grabs from hdf5. It turns out that when a simulation is interrupted, the restart/continue actually grabs the end of timestep time, when it needs to grab the beginning of timestep time (since the step didn't really finish and needs to be rerun). This isn't a problem for uninterrupted simulations (continue feature works as of rn for those) because the final timestep is the final simulation time and continue/restarts from there grab the correct time.

I will open a PR that 1) pushes some tests that expose the bad behavior and will fail with the current code, and then 2) push the fix to show that it actually fixes a problem.