crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 21:25:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

debugging rxn rates

Browse source
This commit is contained in:
guillaume 2018-06-02 22:47:24 -04:00
parent b2c720b296
commit da49be46aa
5 changed files with 92 additions and 7 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -72,11 +72,31 @@ 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
plt.figure()
plt.plot(time/(24*60*60), keff, label="K-effective")
plt.xlabel("Time (days)")
plt.ylabel("Keff")
plt.show()
plt.close()
# 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)")
plt.ylabel("RR (-)")
plt.show()
plt.close('all')

13
examples/python/pincell_depletion/run_depletion.py
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 = 5*24*60*60 # s
final_time = 4*24*60*60 # s
time_steps = np.full(final_time // time_step, time_step)
chain_file = './chain_simple.xml'
@ -146,3 +147,13 @@ 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)')
plt.figure()
plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption")
plt.xlabel("Time (days)")
plt.ylabel("RR (-)")
plt.show()
plt.close('all')

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

@ -76,6 +76,8 @@ def cecm(operator, timesteps, power, print_out=True):
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)
# Deplete for first half of timestep
x_middle = deplete(chain, x[0], op_results[0], dt/2, print_out)

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

@ -60,6 +60,7 @@ 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)]
@ -67,16 +68,62 @@ 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:
x = operator.prev_res[-1].data
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
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)
# 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.prev_res[-1]]
op_results[0].rates = ratio_power[0] * op_results[0].rates[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)
@ -84,6 +131,8 @@ 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))

3
openmc/deplete/results.py
View file

@ -261,6 +261,8 @@ class Results(object):
comm.barrier()
print(self.nuc_to_ind)
# Grab handles
number_dset = handle["/number"]
rxn_dset = handle["/reaction rates"]
@ -305,6 +307,7 @@ 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][:, :, :]