From 8726744ef2931c46dbdbe121f5dec85decf4becc Mon Sep 17 00:00:00 2001 From: guillaume Date: Tue, 29 May 2018 23:57:31 -0400 Subject: [PATCH] added support for changing power at beginning of restart without recomputing rates (no TH feedback) --- .../pincell_depletion/restart_depletion.py | 4 ++-- openmc/deplete/integrator/predictor.py | 10 ++++++---- openmc/deplete/results.py | 19 ++++++++++++++++++- 3 files changed, 26 insertions(+), 7 deletions(-) diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py index e37174ca8..65f73a3da 100644 --- a/examples/python/pincell_depletion/restart_depletion.py +++ b/examples/python/pincell_depletion/restart_depletion.py @@ -18,7 +18,7 @@ final_time = 5*24*60*60 # s time_steps = np.full(final_time // time_step, time_step) chain_file = './chain_simple.xml' -power = 174 # W/cm, for 2D simulations only (use W for 3D) +power = 180 # W/cm, for 2D simulations only (use W for 3D) ############################################################################### # Load previous simulation results @@ -80,6 +80,6 @@ time, keff = results.get_eigenvalue() # Plot eigenvalue as a function of time plt.figure() plt.plot(time/24/60/60, keff, label="K-effective") -plt.xlabel("Time (day)") +plt.xlabel("Time (days)") plt.ylabel("Keff") plt.show() diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index a0ef56dcf..ed6db4efb 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -60,7 +60,6 @@ def predictor(operator, timesteps, power, print_out=True): # If no TH coupling, just re-scale rates by ratio of power for i, (dt, p) in enumerate(zip(timesteps, power)): - # Get beginning-of-timestep concentrations x = [copy.deepcopy(vec)] @@ -69,11 +68,14 @@ def predictor(operator, timesteps, power, print_out=True): if i > 0 or operator.prev_res == None: op_results = [operator(x[0], p)] else: + power_res = operator.prev_res[-1].power + ratio_power = p / power_res + op_results = [operator.prev_res[-1]] - op_results[0].rates = op_results[0].rates[0] + op_results[0].rates = ratio_power * op_results[0].rates[0] # Create results, write to disk - Results.save(operator, x, op_results, [t, t + dt], i + i_res) + Results.save(operator, x, op_results, [t, t + dt], p, i + i_res) # Deplete for full timestep x_end = deplete(chain, x[0], op_results[0], dt, print_out) @@ -87,4 +89,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) + i_res) + Results.save(operator, x, op_results, [t, t], p, len(timesteps) + i_res) diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index ab740e61e..afdfffb74 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -24,6 +24,8 @@ class Results(object): Eigenvalue for each substep. time : list of float Time at beginning, end of step, in seconds. + power : float + Power during time step, in Watts n_mat : int Number of mats. n_nuc : int @@ -49,6 +51,7 @@ class Results(object): def __init__(self): self.k = None self.time = None + self.power = None self.rates = None self.volume = None @@ -237,6 +240,9 @@ class Results(object): handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64') + handle.create_dataset("power", (1, n_stages), maxshape=(None, n_stages), + dtype='float64') + def _to_hdf5(self, handle, index): """Converts results object into an hdf5 object. @@ -259,6 +265,7 @@ class Results(object): rxn_dset = handle["/reaction rates"] eigenvalues_dset = handle["/eigenvalues"] time_dset = handle["/time"] + power_dset = handle["/power"] # Get number of results stored number_shape = list(number_dset.shape) @@ -283,6 +290,10 @@ class Results(object): time_shape[0] = new_shape time_dset.resize(time_shape) + power_shape = list(power_dset.shape) + power_shape[0] = new_shape + power_dset.resize(power_shape) + # If nothing to write, just return if len(self.mat_to_ind) == 0: return @@ -300,6 +311,7 @@ class Results(object): eigenvalues_dset[index, i] = self.k[i] if comm.rank == 0: time_dset[index, :] = self.time + power_dset[index, :] = self.power @classmethod def from_hdf5(cls, handle, step): @@ -319,10 +331,12 @@ class Results(object): number_dset = handle["/number"] eigenvalues_dset = handle["/eigenvalues"] time_dset = handle["/time"] + power_dset = handle["/power"] results.data = number_dset[step, :, :, :] results.k = eigenvalues_dset[step, :] results.time = time_dset[step, :] + results.power = power_dset[step, :] # Reconstruct dictionaries results.volume = OrderedDict() @@ -359,7 +373,7 @@ class Results(object): return results @staticmethod - def save(op, x, op_results, t, step_ind): + def save(op, x, op_results, t, power, step_ind): """Creates and writes depletion results to disk Parameters @@ -372,6 +386,8 @@ class Results(object): Results of applying transport operator t : list of float Time indices. + power : float + Power during time step step_ind : int Step index. @@ -393,5 +409,6 @@ class Results(object): results.k = [r.k for r in op_results] results.rates = [r.rates for r in op_results] results.time = t + results.power = power results.export_to_hdf5("depletion_results.h5", step_ind)