diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py index 6f30c68ea..928c45a28 100644 --- a/examples/python/pincell_depletion/restart_depletion.py +++ b/examples/python/pincell_depletion/restart_depletion.py @@ -1,6 +1,7 @@ import openmc import openmc.deplete import numpy as np +import matplotlib.pyplot as plt ############################################################################### # Simulation Input File Parameters @@ -13,7 +14,7 @@ particles = 1000 # Depletion simulation parameters time_step = 1*24*60*60 # s -final_time = 15*24*60*60 # s +final_time = 5*24*60*60 # s time_steps = np.full(final_time // time_step, time_step) chain_file = './chain_simple.xml' @@ -28,14 +29,15 @@ statepoint = 'statepoint.100.h5' sp = openmc.StatePoint(statepoint) geometry = sp.summary.geometry -# Load previous delpletion results +# Close statepoint and summary files to be able to write over them +sp.summary._f.close() +sp._f.close() + +# Load previous depletion results previous_results = openmc.deplete.ResultsList("depletion_results.h5") -# Reload volumes into geometry -previous_results.transfer_volumes(geometry) - ############################################################################### -# Set transport calculation settings +# Transport calculation settings ############################################################################### # Instantiate a Settings object, set all runtime parameters, and export to XML @@ -59,7 +61,8 @@ settings_file.entropy_mesh = entropy_mesh # Initialize and run depletion calculation ############################################################################### -op = openmc.deplete.Operator(geometry, settings_file, chain_file) +op = openmc.deplete.Operator(geometry, settings_file, chain_file, \ + previous_results) # Perform simulation using the predictor algorithm openmc.deplete.integrator.predictor(op, time_steps, power) @@ -73,6 +76,13 @@ results = openmc.deplete.ResultsList("depletion_results.h5") # Obtain K_eff as a function of time time, keff = results.get_eigenvalue() + +print(time/24/60/60) +print(keff) -# Obtain U235 concentration as a function of time -time, n_U235 = results.get_atoms('1', 'U235') +# Plot eigenvalue as a function of time +plt.figure() +plt.plot(time/24/60/60, keff, label="K-effective") +plt.xlabel("Time (day)") +plt.ylabel("Keff") +plt.show() diff --git a/examples/python/pincell_depletion/run_depletion.py b/examples/python/pincell_depletion/run_depletion.py index a17bf5fd9..03e25fc0e 100644 --- a/examples/python/pincell_depletion/run_depletion.py +++ b/examples/python/pincell_depletion/run_depletion.py @@ -13,7 +13,7 @@ particles = 1000 # Depletion simulation parameters time_step = 1*24*60*60 # s -final_time = 15*24*60*60 # s +final_time = 5*24*60*60 # s time_steps = np.full(final_time // time_step, time_step) chain_file = './chain_simple.xml' @@ -94,7 +94,7 @@ root.add_cells([fuel, gap, clad, water]) geometry = openmc.Geometry(root) ############################################################################### -# Volumes of depletable materials +# Set volumes of depletable materials ############################################################################### # Compute cell areas @@ -105,7 +105,7 @@ area[fuel] = np.pi * fuel_or.coefficients['R'] ** 2 uo2.volume = area[fuel] ############################################################################### -# Transport calculation settings +# Transport calculation settings ############################################################################### # Instantiate a Settings object, set all runtime parameters, and export to XML @@ -135,7 +135,7 @@ op = openmc.deplete.Operator(geometry, settings_file, chain_file) openmc.deplete.integrator.predictor(op, time_steps, power) ############################################################################### -# Read depletion calculation results +# Read depletion calculation results ############################################################################### # Open results file diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index 9be992c16..d3b38aa3a 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -42,17 +42,44 @@ def predictor(operator, timesteps, power, print_out=True): # Generate initial conditions with operator as vec: chain = operator.chain - t = 0.0 - for i, (dt, p) in enumerate(zip(timesteps, power)): - # Get beginning-of-timestep reaction rates - x = [copy.deepcopy(vec)] - op_results = [operator(x[0], p)] - # Create results, write to disk - Results.save(operator, x, op_results, [t, t + dt], i) + # Initialize time + if operator.prev_res == None: + t = 0.0 + else: + t = operator.prev_res.get_eigenvalue()[0][-1] + print("Time", t/24/60/60) + + # Initialize starting index for saving results + if operator.prev_res == None: + i_res = 0 + else: + i_res = len(operator.prev_res.get_eigenvalue()[0]) + print(i_res) + + for i, (dt, p) in enumerate(zip(timesteps, power)): + + # Avoid doing first run if already done in previous calculation + if i > 0 or operator.prev_res == None or p != p_end: + # Get beginning-of-timestep reaction rates + x = [copy.deepcopy(vec)] + op_results = [operator(x[0], p)] + + # Create results, write to disk + Results.save(operator, x, op_results, [t, t + dt], i + i_res) + + else: + x = [copy.deepcopy(vec)] + op_results = operator.prev_res[0] + print(op_results) + op_results = [operator(x[0], p)] + print(op_results) # Deplete for full timestep + #print(x[0], op_results[0]) x_end = deplete(chain, x[0], op_results[0], dt, print_out) + print("Time", t/24/60/60) + print("Step", i + i_res) # Advance time, update vector t += dt @@ -63,4 +90,4 @@ def predictor(operator, timesteps, power, print_out=True): op_results = [operator(x[0], power[-1])] # Create results, write to disk - Results.save(operator, x, op_results, [t, t], len(timesteps)) + Results.save(operator, x, op_results, [t, t], len(timesteps) + i_res) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index b40ff6344..ebd025f94 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -66,6 +66,8 @@ class Operator(TransportOperator): chain_file : str, optional Path to the depletion chain XML file. Defaults to the :envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists. + prev_results : ResultsList, optional + Results from the previous depletion calculation. Defaults to None Attributes ---------- @@ -94,14 +96,25 @@ class Operator(TransportOperator): All burnable material IDs local_mats : list of str All burnable material IDs being managed by a single process + prev_res : ResultsList + Results from the previous depletion run """ - def __init__(self, geometry, settings, chain_file=None): + def __init__(self, geometry, settings, chain_file=None, prev_results=None): super().__init__(chain_file) self.round_number = False self.settings = settings self.geometry = geometry + if prev_results != None: + # Reload volumes into geometry + prev_results.transfer_volumes(geometry) + + # Store previous results in operator + self.prev_res = prev_results + else: + self.prev_res = None + # Clear out OpenMC, create task lists, distribute openmc.reset_auto_ids() self.burnable_mats, volume, nuclides = self._get_burnable_mats() diff --git a/src/summary.F90 b/src/summary.F90 index de3d20a17..45c92062d 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -440,7 +440,6 @@ contains if (m % depletable) then call write_attribute(material_group, "depletable", 1) - call write_attribute(material_group, "volume", m % volume) else call write_attribute(material_group, "depletable", 0) end if