added restart in predictor, working on skipping initial step if already ran

This commit is contained in:
guillaume 2018-05-29 17:45:31 -04:00
parent ed65c15b30
commit 1673b2983a
5 changed files with 72 additions and 23 deletions

View file

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

View file

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

View file

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

View file

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

View file

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