From a85db74b2b1fc624606e5d67d91cc200b738e1e0 Mon Sep 17 00:00:00 2001 From: guillaume Date: Sat, 2 Jun 2018 14:02:14 -0400 Subject: [PATCH] addressed paulromano's review, added support for restart in cecm integrator --- .../pincell_depletion/restart_depletion.py | 14 ++++---- .../python/pincell_depletion/run_depletion.py | 2 +- openmc/deplete/integrator/cecm.py | 36 ++++++++++++++++--- openmc/deplete/integrator/predictor.py | 7 ++-- openmc/deplete/operator.py | 4 +-- openmc/statepoint.py | 3 -- tests/dummy_operator.py | 1 - 7 files changed, 44 insertions(+), 23 deletions(-) diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py index a99565235..cfca2adf0 100644 --- a/examples/python/pincell_depletion/restart_depletion.py +++ b/examples/python/pincell_depletion/restart_depletion.py @@ -26,11 +26,8 @@ power = 174 # W/cm, for 2D simulations only (use W for 3D) # Load geometry from statepoint statepoint = 'statepoint.100.h5' -sp = openmc.StatePoint(statepoint) -geometry = sp.summary.geometry - -# Close statepoint and summary files to be able to write over them -sp.close() +with openmc.StatePoint(statepoint) as sp: + geometry = sp.summary.geometry # Load previous depletion results previous_results = openmc.deplete.ResultsList("depletion_results.h5") @@ -60,8 +57,8 @@ settings_file.entropy_mesh = entropy_mesh # Initialize and run depletion calculation ############################################################################### -op = openmc.deplete.Operator(geometry, settings_file, chain_file, \ - previous_results) +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) @@ -78,7 +75,8 @@ 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.plot(time/(24*60*60), keff, label="K-effective") plt.xlabel("Time (days)") plt.ylabel("Keff") plt.show() +plt.close() diff --git a/examples/python/pincell_depletion/run_depletion.py b/examples/python/pincell_depletion/run_depletion.py index 03e25fc0e..6600bec06 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.predictor(op, time_steps, power) +openmc.deplete.integrator.cecm(op, time_steps, power) ############################################################################### # Read depletion calculation results diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index b2c766e99..0ab430226 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -47,11 +47,36 @@ def cecm(operator, timesteps, power, print_out=True): # Generate initial conditions with operator as vec: chain = operator.chain - t = 0.0 + + # Initialize time + if operator.prev_res is None: + t = 0.0 + else: + t = operator.prev_res[-1].time[-1] + + # Initialize starting index for saving results + if operator.prev_res is None: + i_res = 0 + else: + i_res = len(operator.prev_res) + for i, (dt, p) in enumerate(zip(timesteps, power)): - # Get beginning-of-timestep reaction rates + # Get beginning-of-timestep concentrations x = [copy.deepcopy(vec)] - op_results = [operator(x[0], p)] + + # Get beginning-of-timestep reaction rates + # Avoid doing first transport run if already done in previous + # calculation + if i > 0 or operator.prev_res is 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 = ratio_power[0] * op_results[0].rates[0] + op_results[0].k = op_results[0].k[0] # Deplete for first half of timestep x_middle = deplete(chain, x[0], op_results[0], dt/2, print_out) @@ -61,10 +86,11 @@ def cecm(operator, timesteps, power, print_out=True): op_results.append(operator(x_middle, p)) # Deplete for full timestep using beginning-of-step materials + # and middle-of-timestep reaction rates x_end = deplete(chain, x[0], op_results[1], dt, print_out) # Create results, write to disk - Results.save(operator, x, op_results, [t, t + dt], p, i) + Results.save(operator, x, op_results, [t, t + dt], p, i_res + i) # Advance time, update vector t += dt @@ -75,4 +101,4 @@ def cecm(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], p, len(timesteps)) + Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps)) diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index ea9880294..4711b7f85 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -66,13 +66,14 @@ def predictor(operator, timesteps, power, print_out=True): op_results = [operator(x[0], p)] # Create results, write to disk - Results.save(operator, x, op_results, [t, t + dt], p, i + i_res) + Results.save(operator, x, op_results, [t, t + dt], p, i_res + i) else: power_res = operator.prev_res[-1].power + print(power_res) ratio_power = p / power_res op_results = [operator.prev_res[-1]] - op_results[0].rates = ratio_power * op_results[0].rates[0] + op_results[0].rates = ratio_power[0] * op_results[0].rates[0] # Deplete for full timestep x_end = deplete(chain, x[0], op_results[0], dt, print_out) @@ -86,4 +87,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], p, len(timesteps) + i_res) + Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps)) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 71031be02..9ffe74063 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -296,10 +296,10 @@ class Operator(TransportOperator): # Get nuclide lists from geometry and depletion results depl_nuc = prev_res[-1].nuc_to_ind.keys() geom_nuc_densities = mat.get_nuclide_atom_densities() - geom_nuc = [x[0] for x in list(geom_nuc_densities.values())] + geom_nuc = {x[0] for x in geom_nuc_densities.values()} # Merge lists of nuclides - nuc_set = set(depl_nuc) | set(geom_nuc) + nuc_set = set(depl_nuc) | geom_nuc for nuclide in nuc_set: if nuclide in depl_nuc: diff --git a/openmc/statepoint.py b/openmc/statepoint.py index fcc076b94..a200e900e 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -153,9 +153,6 @@ class StatePoint(object): if self._summary is not None: self._summary._f.close() - def close(self): - self.__exit__() - @property def cmfd_on(self): return self._f.attrs['cmfd_on'] > 0 diff --git a/tests/dummy_operator.py b/tests/dummy_operator.py index 867d30dd6..fd230635d 100644 --- a/tests/dummy_operator.py +++ b/tests/dummy_operator.py @@ -19,7 +19,6 @@ class DummyOperator(TransportOperator): """ def __init__(self): self.prev_res = None - pass def __call__(self, vec, power, print_out=False): """Evaluates F(y)