diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py index 9ed20485c..f9324a8cb 100644 --- a/examples/python/pincell_depletion/restart_depletion.py +++ b/examples/python/pincell_depletion/restart_depletion.py @@ -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)") diff --git a/examples/python/pincell_depletion/run_depletion.py b/examples/python/pincell_depletion/run_depletion.py index 0c58b27d2..a2a5e0e7f 100644 --- a/examples/python/pincell_depletion/run_depletion.py +++ b/examples/python/pincell_depletion/run_depletion.py @@ -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)") diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index be274bf8f..0b17b49b5 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -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) diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index 980eb0cde..a45017813 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -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)) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 24e8843ab..a0b2417c8 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -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] diff --git a/openmc/deplete/reaction_rates.py b/openmc/deplete/reaction_rates.py index fddb88b19..cea2f1997 100644 --- a/openmc/deplete/reaction_rates.py +++ b/openmc/deplete/reaction_rates.py @@ -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 diff --git a/openmc/deplete/results.py b/openmc/deplete/results.py index 259ba6f64..6170a4464 100644 --- a/openmc/deplete/results.py +++ b/openmc/deplete/results.py @@ -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) diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index c19f1c5f6..8cade45b1 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -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 = []