From 563e3436e462fbf96049dff205e4367bacb5cfdd Mon Sep 17 00:00:00 2001 From: guillaume Date: Tue, 29 May 2018 19:09:38 -0400 Subject: [PATCH] getting initial nuclide densities from depletion results, other choice would be to update python densities, but depletion results are needed anyway --- .../pincell_depletion/restart_depletion.py | 2 +- openmc/deplete/integrator/predictor.py | 3 +- openmc/deplete/operator.py | 49 ++++++++++++++++--- 3 files changed, 43 insertions(+), 11 deletions(-) diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py index 928c45a28..7bb365fe4 100644 --- a/examples/python/pincell_depletion/restart_depletion.py +++ b/examples/python/pincell_depletion/restart_depletion.py @@ -40,7 +40,7 @@ previous_results = openmc.deplete.ResultsList("depletion_results.h5") # Transport calculation settings ############################################################################### -# Instantiate a Settings object, set all runtime parameters, and export to XML +# Instantiate a Settings object, set all runtime parameters settings_file = openmc.Settings() settings_file.batches = batches settings_file.inactive = inactive diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index 7bc6b9546..f55393b54 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -55,7 +55,7 @@ def predictor(operator, timesteps, power, print_out=True): else: i_res = len(operator.prev_res.get_eigenvalue()[0]) - #TODO : Get last time step power from previous results, and run a + #TODO : Get last time step power from previous results, and run a # new calculation if different from power at the first time step # If no TH coupling, just re-scale rates by ratio of power @@ -68,7 +68,6 @@ def predictor(operator, timesteps, power, print_out=True): # Avoid doing first run if already done in previous calculation if i > 0 or operator.prev_res == None: op_results = [operator(x[0], p)] - else: op_results = [operator.prev_res[-1]] op_results[0].rates = op_results[0].rates[0] diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 44c57f160..49474a925 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -77,7 +77,7 @@ class Operator(TransportOperator): OpenMC settings object dilute_initial : float Initial atom density to add for nuclides that are zero in initial - condition to ensure they exist in the decay chain. Only done for + condition to ensure they exist in the decay chain. Only done for nuclides with reaction rates. Defaults to 1.0e3. output_dir : pathlib.Path Path to output directory to save results. @@ -112,6 +112,9 @@ class Operator(TransportOperator): # Store previous results in operator self.prev_res = prev_results + + # Get number densities from previous results + #self.number = prev_results else: self.prev_res = None @@ -125,8 +128,9 @@ class Operator(TransportOperator): self._burnable_nucs = [nuc for nuc in self.nuclides_with_data if nuc in self.chain] - # Extract number densities from the geometry - self._extract_number(self.local_mats, volume, nuclides) + # Extract number densities from the geometry / previous depletion run + self._extract_number(self.local_mats, volume, nuclides, \ + self.prev_res) # Create reaction rates array self.reaction_rates = ReactionRates( @@ -225,7 +229,7 @@ class Operator(TransportOperator): return burnable_mats, volume, nuclides - def _extract_number(self, local_mats, volume, nuclides): + def _extract_number(self, local_mats, volume, nuclides, prev_res=None): """Construct AtomNumber using geometry Parameters @@ -244,10 +248,18 @@ class Operator(TransportOperator): for nuc in self._burnable_nucs: self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial) - # Now extract the number densities and store - for mat in self.geometry.get_all_materials().values(): - if str(mat.id) in local_mats: - self._set_number_from_mat(mat) + # Now extract and store the number densities + # From the geometry if no previous depletion results + if prev_res == None: + for mat in self.geometry.get_all_materials().values(): + if str(mat.id) in local_mats: + self._set_number_from_mat(mat) + + # Else from previous depletion results + else: + for mat in self.geometry.get_all_materials().values(): + if str(mat.id) in local_mats: + self._set_number_from_results(mat, prev_res) def _set_number_from_mat(self, mat): """Extracts material and number densities from openmc.Material @@ -262,6 +274,23 @@ class Operator(TransportOperator): for nuclide, density in mat.get_nuclide_atom_densities().values(): number = density * 1.0e24 + print(nuclide) + self.number.set_atom_density(mat_id, nuclide, number) + + def _set_number_from_results(self, mat, prev_res): + """Extracts material and number densities from previous results + + Parameters + ---------- + mat : openmc.Material + The material to read from + + """ + mat_id = str(mat.id) + + for nuclide, density in mat.get_nuclide_atom_densities().values(): + number = density * 1.0e24 + print(nuclide, number) self.number.set_atom_density(mat_id, nuclide, number) def initial_condition(self): @@ -325,9 +354,13 @@ class Operator(TransportOperator): " is negative (density = ", val, " at/barn-cm)") number_i[mat, nuc] = 0.0 + # Update densities on C API side mat_internal = openmc.capi.materials[int(mat)] mat_internal.set_densities(nuclides, densities) + #TODO Update densities on the Python side, otherwise the + # summary.h5 file contains densities at the first time step + def _generate_materials_xml(self): """Creates materials.xml from self.number.