diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py index cfca2adf0..9ed20485c 100644 --- a/examples/python/pincell_depletion/restart_depletion.py +++ b/examples/python/pincell_depletion/restart_depletion.py @@ -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') diff --git a/examples/python/pincell_depletion/run_depletion.py b/examples/python/pincell_depletion/run_depletion.py index 03e25fc0e..0c58b27d2 100644 --- a/examples/python/pincell_depletion/run_depletion.py +++ b/examples/python/pincell_depletion/run_depletion.py @@ -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') diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index a80f7ee26..be274bf8f 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -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) diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index 6918b222e..980eb0cde 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -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)) diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index f31f5672f..259ba6f64 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -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][:, :, :]