diff --git a/examples/python/pincell_depletion/run_depletion.py b/examples/python/pincell_depletion/run_depletion.py index 6600bec06..03e25fc0e 100644 --- a/examples/python/pincell_depletion/run_depletion.py +++ b/examples/python/pincell_depletion/run_depletion.py @@ -132,7 +132,7 @@ settings_file.entropy_mesh = entropy_mesh op = openmc.deplete.Operator(geometry, settings_file, chain_file) # Perform simulation using the predictor algorithm -openmc.deplete.integrator.cecm(op, time_steps, power) +openmc.deplete.integrator.predictor(op, time_steps, power) ############################################################################### # Read depletion calculation results diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index 0ab430226..a80f7ee26 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -61,16 +61,15 @@ def cecm(operator, timesteps, power, print_out=True): i_res = len(operator.prev_res) for i, (dt, p) in enumerate(zip(timesteps, power)): - # Get beginning-of-timestep concentrations - x = [copy.deepcopy(vec)] - - # Get beginning-of-timestep reaction rates + # Get beginning-of-timestep concentrations and reaction rates # Avoid doing first transport run if already done in previous # calculation if i > 0 or operator.prev_res is None: + x = [copy.deepcopy(vec)] op_results = [operator(x[0], p)] else: + x = [operator.prev_res[-1].data[0]] power_res = operator.prev_res[-1].power ratio_power = p / power_res diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index 4711b7f85..6918b222e 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -38,6 +38,7 @@ def predictor(operator, timesteps, power, print_out=True): """ if not isinstance(power, Iterable): power = [power]*len(timesteps) + print(power) # Generate initial conditions with operator as vec: @@ -56,20 +57,18 @@ def predictor(operator, timesteps, power, print_out=True): i_res = len(operator.prev_res) - 1 for i, (dt, p) in enumerate(zip(timesteps, power)): - # Get beginning-of-timestep concentrations - x = [copy.deepcopy(vec)] - - # Get beginning-of-timestep reaction rates + # Get beginning-of-timestep concentrations and reaction rates # Avoid doing first transport run if already done in previous # calculation if i > 0 or operator.prev_res is None: + x = [copy.deepcopy(vec)] op_results = [operator(x[0], p)] # Create results, write to disk Results.save(operator, x, op_results, [t, t + dt], p, i_res + i) else: + x = operator.prev_res[-1].data power_res = operator.prev_res[-1].power - print(power_res) ratio_power = p / power_res op_results = [operator.prev_res[-1]] diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 9ffe74063..24e8843ab 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -134,6 +134,7 @@ class Operator(TransportOperator): self.reaction_rates = ReactionRates( self.local_mats, self._burnable_nucs, self.chain.reactions) + def __call__(self, vec, power, print_out=True): """Runs a simulation. diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index 1d0dc3a0d..f31f5672f 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -5,6 +5,7 @@ Contains results generation and saving capabilities. from collections import OrderedDict import copy +from warnings import warn import numpy as np import h5py @@ -395,8 +396,18 @@ class Results(object): # Get indexing terms vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info() - # Create results + # For a restart calculation, limit number of stages saved to meet the + # format of the hdf5 file stages = len(x) + offset = 0 + if op.prev_res is not None and op.prev_res[0].n_stages < stages: + offset = stages - op.prev_res[0].n_stages + stages = min(stages, op.prev_res[0].n_stages) + warn("Number of restart integrator stages saved limited by initial" + " depletion integrator choice to {}" + .format(op.prev_res[0].n_stages)) + + # Create results results = Results() results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages) @@ -404,7 +415,7 @@ class Results(object): for i in range(stages): for mat_i in range(n_mat): - results[i, mat_i, :] = x[i][mat_i][:] + results[i, mat_i, :] = x[offset + i][mat_i][:] results.k = [r.k for r in op_results] results.rates = [r.rates for r in op_results] diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py index fd230635d..c66070ad4 100644 --- a/tests/dummy_operator.py +++ b/tests/dummy_operator.py @@ -17,8 +17,8 @@ class DummyOperator(TransportOperator): y_2(1.5) ~ 3.1726475740397628 """ - def __init__(self): - self.prev_res = None + def __init__(self, previous_results=None): + self.prev_res = previous_results def __call__(self, vec, power, print_out=False): """Evaluates F(y) diff --git a/tests/unit_tests/test_deplete_predictor.py b/tests/unit_tests/test_deplete_predictor.py index 50803e508..b39cc7dca 100644 --- a/tests/unit_tests/test_deplete_predictor.py +++ b/tests/unit_tests/test_deplete_predictor.py @@ -10,7 +10,7 @@ from tests import dummy_operator def test_predictor(run_in_tmpdir): - """Integral regression test of integrator algorithm using predictor/corrector""" + """Integral regression test of integrator algorithm using predictor""" op = dummy_operator.DummyOperator() op.output_dir = "test_integrator_regression"