crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

addressed paulromano's review, added support for restart in cecm integrator

Browse source
This commit is contained in:
guillaume 2018-06-02 14:02:14 -04:00
parent f25cdb93f2
commit a85db74b2b
7 changed files with 44 additions and 23 deletions
Download patch file Download diff file Expand all files Collapse all files

14
examples/python/pincell_depletion/restart_depletion.py
View file

@ -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()

2
examples/python/pincell_depletion/run_depletion.py
View file

@ -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

36
openmc/deplete/integrator/cecm.py
View file

@ -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))

7
openmc/deplete/integrator/predictor.py
View file

@ -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))

4
openmc/deplete/operator.py
View file

@ -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:

3
openmc/statepoint.py
View file

@ -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

1
tests/dummy_operator.py
View file

@ -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)