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

BUG fix, indexing of nuclides in reactionrates arrays was being ignored + cleaned debugging lines + added plots to examples

Browse source
This commit is contained in:
guillaume 2018-06-03 21:21:39 -04:00
parent da49be46aa
commit 9de6108ac0
8 changed files with 72 additions and 82 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -73,27 +73,28 @@ results = openmc.deplete.ResultsList("depletion_results.h5")
# Obtain K_eff as a function of time
time, keff = results.get_eigenvalue()
# Plot eigenvalue as a function of time
# Obtain U235 concentration as a function of time
time, n_U235 = results.get_atoms('1', 'U235')
# Obtain Xe135 absorption as a function of time
time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
###############################################################################
# Generate plots
###############################################################################
plt.figure()
plt.plot(time/(24*60*60), keff, label="K-effective")
plt.xlabel("Time (days)")
plt.ylabel("Keff")
plt.show()
# Obtain U235 concentration as a function of time
time, n_U235 = results.get_atoms('1', 'U235')
print(time/(24*60*60))
print(n_U235)
plt.figure()
plt.plot(time/(24*60*60), n_U235, label="U 235")
plt.xlabel("Time (days)")
plt.ylabel("n U5 (-)")
plt.show()
# Obtain Xe135 absorption as a function of time
time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
plt.figure()
plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption")
plt.xlabel("Time (days)")

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

@ -14,7 +14,7 @@ particles = 1000
# Depletion simulation parameters
time_step = 1*24*60*60 # s
final_time = 4*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'
@ -136,7 +136,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
@ -144,13 +144,29 @@ results = openmc.deplete.ResultsList("depletion_results.h5")
# Obtain K_eff as a function of time
time, keff = results.get_eigenvalue()
# Obtain U235 concentration as a function of time
time, n_U235 = results.get_atoms('1', 'U235')
# Obtain Xe135 absorption as a function of time
time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
###############################################################################
# Generate plots
###############################################################################
plt.figure()
plt.plot(time/(24*60*60), keff, label="K-effective")
plt.xlabel("Time (days)")
plt.ylabel("Keff")
plt.show()
plt.figure()
plt.plot(time/(24*60*60), n_U235, label="U 235")
plt.xlabel("Time (days)")
plt.ylabel("n U5 (-)")
plt.show()
plt.figure()
plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption")
plt.xlabel("Time (days)")

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

@ -69,15 +69,20 @@ def cecm(operator, timesteps, power, print_out=True):
op_results = [operator(x[0], p)]
else:
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
# Get rates
op_results = [operator.prev_res[-1]]
op_results[0].rates = op_results[0].rates[0]
# Set first stage value of keff
op_results[0].k = op_results[0].k[0]
# Scale reaction rates by ratio of powers
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]
print(x)
print(op_results[0].rates)
op_results[0].rates[0] *= ratio_power[0]
# Deplete for first half of timestep
x_middle = deplete(chain, x[0], op_results[0], dt/2, print_out)

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

@ -38,7 +38,6 @@ def predictor(operator, timesteps, power, print_out=True):
"""
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
print(power)
# Generate initial conditions
with operator as vec:
@ -60,7 +59,6 @@ def predictor(operator, timesteps, power, print_out=True):
# Get beginning-of-timestep concentrations and reaction rates
# Avoid doing first transport run if already done in previous
# calculation
print("i", i, "sp i", i_res + i)
if i > 0 or operator.prev_res is None:
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], p)]
@ -68,62 +66,21 @@ def predictor(operator, timesteps, power, print_out=True):
# Create results, write to disk
Results.save(operator, x, op_results, [t, t + dt], p, i_res + i)
else:
print("Data", operator.prev_res[-1].data)
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
print(x)
x = [copy.deepcopy(vec)]
print(x)
# Get rates, indexed by mat_to_ind in previous results
# Get rates
op_results = [operator.prev_res[-1]]
nuc_to_ind_current = {nuc: i for i, nuc in \
enumerate(operator.number.burnable_nuclides)}
nuc_to_ind_res = [*operator.prev_res[-1].nuc_to_ind]
nuc_to_ind_res = {nuc: i for i, nuc in enumerate(nuc_to_ind_res)}
print(nuc_to_ind_current)
print(nuc_to_ind_current.keys())
print(nuc_to_ind_res)
print(operator.prev_res[-1].nuc_to_ind.keys())
match = [nuc_to_ind_res[nuc] for nuc in nuc_to_ind_current.keys()]
print(match)
match = [nuc_to_ind_current[nuc] for nuc in nuc_to_ind_res.keys()]
print(match)
match = [operator.prev_res[-1].nuc_to_ind[nuc] for nuc in nuc_to_ind_current.keys()]
print(match)
match = [nuc_to_ind_current[nuc] for nuc in operator.prev_res[-1].nuc_to_ind.keys()]
print(match)
match = range(9)
print(match)
sv = op_results[0].rates[0][0] ######
print("Index of nuclides in rr", sv.index_nuc)
print("Index of rxn in rr", sv.index_rx)
print("Index of mat in rr", sv.index_mat)
#x = [[operator.prev_res[-1].data[0][0][[nuc_to_ind_res[nuc] for nuc in nuc_to_ind_current.keys()]]]]
print(x)
op_results[0].rates = op_results[0].rates[0]
# Scale reaction rates by ratio of powers
power_res = operator.prev_res[-1].power
ratio_power = p / power_res
op_results[0].rates[0] *= ratio_power[0]
op_results = [operator(x[0], p)]
print("old", sv)
print("new", op_results[0].rates)
print(operator.prev_res[-1].nuc_to_ind)
# Deplete for full timestep
print("x[0]", x[0])
x_end = deplete(chain, x[0], op_results[0], dt, print_out)
print("xend", x_end)
# Advance time, update vector
t += dt
vec = copy.deepcopy(x_end)
@ -131,8 +88,6 @@ def predictor(operator, timesteps, power, print_out=True):
# Perform one last simulation
x = [copy.deepcopy(vec)]
op_results = [operator(x[0], power[-1])]
print("Final power", power[-1])
# Create results, write to disk
print("i" , len(timesteps), "sp i" , i_res + len(timesteps))
Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))

15
openmc/deplete/operator.py
View file

@ -124,8 +124,10 @@ class Operator(TransportOperator):
# Determine which nuclides have incident neutron data
self.nuclides_with_data = self._get_nuclides_with_data()
self._burnable_nucs = [nuc for nuc in self.nuclides_with_data
if nuc in self.chain]
# Select nuclides with data that are also in the chain
self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides
if nuc.name in self.nuclides_with_data]
# Extract number densities from the geometry / previous depletion run
self._extract_number(self.local_mats, volume, nuclides, self.prev_res)
@ -295,14 +297,13 @@ class Operator(TransportOperator):
mat_id = str(mat.id)
# Get nuclide lists from geometry and depletion results
depl_nuc = prev_res[-1].nuc_to_ind.keys()
depl_nuc = prev_res[-1].nuc_to_ind
geom_nuc_densities = mat.get_nuclide_atom_densities()
geom_nuc = {x[0] for x in geom_nuc_densities.values()}
# Merge lists of nuclides
nuc_set = set(depl_nuc) | geom_nuc
# Merge lists of nuclides, with the same order for every calculation
geom_nuc_densities.update(depl_nuc)
for nuclide in nuc_set:
for nuclide in geom_nuc_densities.keys():
if nuclide in depl_nuc:
concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1]
volume = prev_res[-1].volume[mat_id]

21
openmc/deplete/reaction_rates.py
View file

@ -22,6 +22,9 @@ class ReactionRates(np.ndarray):
Depletable nuclides
reactions : list of str
Transmutation reactions being tracked
from_results : boolean
If the reaction rates are loaded from results, indexing dictionnaries
need to be kept the same.
Attributes
----------
@ -47,16 +50,24 @@ class ReactionRates(np.ndarray):
# the __array_finalize__ method (discussed here:
# https://docs.scipy.org/doc/numpy/user/basics.subclassing.html)
def __new__(cls, local_mats, nuclides, reactions):
def __new__(cls, local_mats, nuclides, reactions, from_results=False):
# Create appropriately-sized zeroed-out ndarray
shape = (len(local_mats), len(nuclides), len(reactions))
obj = super().__new__(cls, shape)
obj[:] = 0.0
# Add mapping attributes
obj.index_mat = {mat: i for i, mat in enumerate(local_mats)}
obj.index_nuc = {nuc: i for i, nuc in enumerate(nuclides)}
obj.index_rx = {rx: i for i, rx in enumerate(reactions)}
# Add mapping attributes, keep same indexing if from depletion_results
if from_results:
obj.index_mat = local_mats
obj.index_nuc = nuclides
obj.index_rx = reactions
# Else, assumes that reaction rates are ordered the same way as
# the lists of local_mats, nuclides and reactions (or keys if these
# are dictionnaries)
else:
obj.index_mat = {mat: i for i, mat in enumerate(local_mats)}
obj.index_nuc = {nuc: i for i, nuc in enumerate(nuclides)}
obj.index_rx = {rx: i for i, rx in enumerate(reactions)}
return obj

6
openmc/deplete/results.py
View file

@ -165,6 +165,7 @@ class Results(object):
else:
kwargs = {}
# Write new file if first time step, else add to existing file
kwargs['mode'] = "w" if step == 0 else "a"
with h5py.File(filename, **kwargs) as handle:
@ -261,8 +262,6 @@ class Results(object):
comm.barrier()
print(self.nuc_to_ind)
# Grab handles
number_dset = handle["/number"]
rxn_dset = handle["/reaction rates"]
@ -307,7 +306,6 @@ class Results(object):
inds = [self.mat_to_hdf5_ind[mat] for mat in self.mat_to_ind]
low = min(inds)
high = max(inds)
print("indexes", inds)
for i in range(n_stages):
number_dset[index, i, low:high+1, :] = self.data[i, :, :]
rxn_dset[index, i, low:high+1, :, :] = self.rates[i][:, :, :]
@ -369,7 +367,7 @@ class Results(object):
results.rates = []
# Reconstruct reactions
for i in range(results.n_stages):
rate = ReactionRates(results.mat_to_ind, rxn_nuc_to_ind, rxn_to_ind)
rate = ReactionRates(results.mat_to_ind, rxn_nuc_to_ind, rxn_to_ind, True)
rate[:] = handle["/reaction rates"][step, i, :, :, :]
results.rates.append(rate)

3
tests/unit_tests/test_deplete_integrator.py
View file

@ -25,6 +25,9 @@ def test_results_save(run_in_tmpdir):
# Mock geometry
op = MagicMock()
# Avoid DummyOperator thinking it's doing a restart calculation
op.prev_res = None
vol_dict = {}
full_burn_list = []