From bdec9759236b02bc1a35641922494ae7377955c7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 9 Jan 2022 13:11:27 -0600 Subject: [PATCH 1/5] Fix material density in ResultsList.export_to_materials --- examples/pincell_depletion/restart_depletion.py | 8 +++++--- examples/pincell_depletion/run_depletion.py | 4 ++-- openmc/deplete/results_list.py | 1 + 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/examples/pincell_depletion/restart_depletion.py b/examples/pincell_depletion/restart_depletion.py index 0bbab57a58..6a38b4bd98 100644 --- a/examples/pincell_depletion/restart_depletion.py +++ b/examples/pincell_depletion/restart_depletion.py @@ -39,9 +39,11 @@ settings.entropy_mesh = entropy_mesh # Initialize and run depletion calculation ############################################################################### +model = openmc.Model(geometry=geometry, settings=settings) + # Create depletion "operator" -chain_file = './chain_simple.xml' -op = openmc.deplete.Operator(geometry, settings, chain_file, previous_results) +chain_file = 'chain_simple.xml' +op = openmc.deplete.Operator(model, chain_file, previous_results) # Perform simulation using the predictor algorithm time_steps = [1.0, 1.0, 1.0, 1.0, 1.0] # days @@ -71,7 +73,7 @@ time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') days = 24*60*60 plt.figure() -plt.plot(time/days, keff, label="K-effective") +plt.plot(time/days, keff[:, 0], label="K-effective") plt.xlabel("Time (days)") plt.ylabel("Keff") plt.show() diff --git a/examples/pincell_depletion/run_depletion.py b/examples/pincell_depletion/run_depletion.py index 72e71afe5a..1343e8b593 100644 --- a/examples/pincell_depletion/run_depletion.py +++ b/examples/pincell_depletion/run_depletion.py @@ -87,7 +87,7 @@ settings.entropy_mesh = entropy_mesh model = openmc.Model(geometry=geometry, settings=settings) # Create depletion "operator" -chain_file = './chain_simple.xml' +chain_file = 'chain_simple.xml' op = openmc.deplete.Operator(model, chain_file) # Perform simulation using the predictor algorithm @@ -118,7 +118,7 @@ time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') days = 24*60*60 plt.figure() -plt.plot(time/days, keff, label="K-effective") +plt.plot(time/days, keff[:, 0], label="K-effective") plt.xlabel("Time (days)") plt.ylabel("Keff") plt.show() diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index 2dc486feb3..a9a0e52119 100644 --- a/openmc/deplete/results_list.py +++ b/openmc/deplete/results_list.py @@ -361,6 +361,7 @@ class ResultsList(list): mat_id = str(mat.id) if mat_id in result.mat_to_ind: mat.volume = result.volume[mat_id] + mat.set_density('sum') for nuc in result.nuc_to_ind: if nuc not in available_cross_sections: continue From 66835f465f51c26c1efdde7912e4fad13fdea042 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 9 Jan 2022 13:14:25 -0600 Subject: [PATCH 2/5] Fix undefined variable in Operator.__init__ --- openmc/deplete/operator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index d727b26317..fce0670e51 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -309,7 +309,7 @@ class Operator(TransportOperator): if normalization_mode == "fission-q": self._normalization_helper = ChainFissionHelper() elif normalization_mode == "energy-deposition": - score = "heating" if settings.photon_transport else "heating-local" + score = "heating" if self.settings.photon_transport else "heating-local" self._normalization_helper = EnergyScoreHelper(score) else: self._normalization_helper = SourceRateHelper() From 141a32487a64d9b49097ad31ee4b567106fa7769 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 9 Jan 2022 13:16:02 -0600 Subject: [PATCH 3/5] Don't overwrite Operator.materials.cross_sections if already defined --- openmc/deplete/operator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index fce0670e51..7e3b34d0c1 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -642,7 +642,7 @@ class Operator(TransportOperator): if mfile.exists(): tree = ET.parse(str(mfile)) xs = tree.find('cross_sections') - if xs is not None: + if xs is not None and self.materials.cross_sections is None: self.materials.cross_sections = xs.text self.materials.export_to_xml() From 2ccbf0669bf63f8adc096427e6ccec242bb30b4b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 9 Jan 2022 15:47:08 -0600 Subject: [PATCH 4/5] Update pincell depletion example plots --- .../pincell_depletion/restart_depletion.py | 25 +++++++++---------- examples/pincell_depletion/run_depletion.py | 25 +++++++++---------- 2 files changed, 24 insertions(+), 26 deletions(-) diff --git a/examples/pincell_depletion/restart_depletion.py b/examples/pincell_depletion/restart_depletion.py index 6a38b4bd98..fd0a4a8dde 100644 --- a/examples/pincell_depletion/restart_depletion.py +++ b/examples/pincell_depletion/restart_depletion.py @@ -72,21 +72,20 @@ time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') ############################################################################### days = 24*60*60 -plt.figure() -plt.plot(time/days, keff[:, 0], label="K-effective") -plt.xlabel("Time (days)") -plt.ylabel("Keff") +fig, ax = plt.subplots() +ax.errorbar(time/days, keff[:, 0], keff[:, 1], label="K-effective") +ax.set_xlabel("Time [d]") +ax.set_ylabel("Keff") plt.show() -plt.figure() -plt.plot(time/days, n_U235, label="U 235") -plt.xlabel("Time (days)") -plt.ylabel("n U5 (-)") +fig, ax = plt.subplots() +ax.plot(time/days, n_U235, label="U235") +ax.set_xlabel("Time [d]") +ax.set_ylabel("U235 atoms") plt.show() -plt.figure() -plt.plot(time/days, Xe_capture, label="Xe135 capture") -plt.xlabel("Time (days)") -plt.ylabel("RR (-)") +fig, ax = plt.subplots() +ax.plot(time/days, Xe_capture, label="Xe135 capture") +ax.set_xlabel("Time [d]") +ax.set_ylabel("Xe135 capture rate") plt.show() -plt.close('all') diff --git a/examples/pincell_depletion/run_depletion.py b/examples/pincell_depletion/run_depletion.py index 1343e8b593..45dfd2a875 100644 --- a/examples/pincell_depletion/run_depletion.py +++ b/examples/pincell_depletion/run_depletion.py @@ -117,21 +117,20 @@ time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') ############################################################################### days = 24*60*60 -plt.figure() -plt.plot(time/days, keff[:, 0], label="K-effective") -plt.xlabel("Time (days)") -plt.ylabel("Keff") +fig, ax = plt.subplots() +ax.errorbar(time/days, keff[:, 0], keff[:, 1], label="K-effective") +ax.set_xlabel("Time [d]") +ax.set_ylabel("Keff") plt.show() -plt.figure() -plt.plot(time/days, n_U235, label="U235") -plt.xlabel("Time (days)") -plt.ylabel("n U5 (-)") +fig, ax = plt.subplots() +ax.plot(time/days, n_U235, label="U235") +ax.set_xlabel("Time [d]") +ax.set_ylabel("U235 atoms") plt.show() -plt.figure() -plt.plot(time/days, Xe_capture, label="Xe135 capture") -plt.xlabel("Time (days)") -plt.ylabel("RR (-)") +fig, ax = plt.subplots() +ax.plot(time/days, Xe_capture, label="Xe135 capture") +ax.set_xlabel("Time [d]") +ax.set_ylabel("Xe135 capture rate") plt.show() -plt.close('all') From 5395d32fde9971ffe359211b2557db51b61833fe Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 10 Jan 2022 07:39:29 -0600 Subject: [PATCH 5/5] Fix use of _{materials,cells}_by_id and _by_name attributes in Model --- openmc/model/model.py | 69 ++++++++++++++++++++++++++++--------------- 1 file changed, 45 insertions(+), 24 deletions(-) diff --git a/openmc/model/model.py b/openmc/model/model.py index d253bfb5b0..998992c903 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1,4 +1,5 @@ from collections.abc import Iterable +from functools import lru_cache import os from pathlib import Path from numbers import Integral @@ -88,30 +89,6 @@ class Model: if plots is not None: self.plots = plots - # Store dictionaries to the materials and cells by ID and names - if materials is None: - mats = self.geometry.get_all_materials().values() - else: - mats = self.materials - cells = self.geometry.get_all_cells() - # Get the ID maps - self._materials_by_id = {mat.id: mat for mat in mats} - self._cells_by_id = {cell.id: cell for cell in cells.values()} - - # Get the names maps, but since names are not unique, store a list for - # each name key. In this way when the user requests a change by a name, - # the change will be applied to all of the same name. - self._cells_by_name = {} - for cell in cells.values(): - if cell.name not in self._cells_by_name: - self._cells_by_name[cell.name] = set() - self._cells_by_name[cell.name].add(cell) - self._materials_by_name = {} - for mat in mats: - if mat.name not in self._materials_by_name: - self._materials_by_name[mat.name] = set() - self._materials_by_name[mat.name].add(mat) - @property def geometry(self): return self._geometry @@ -140,6 +117,50 @@ class Model: except ImportError: return False + @property + @lru_cache(maxsize=None) + def _materials_by_id(self): + """Dictionary mapping material ID --> material""" + if self.materials is None: + mats = self.geometry.get_all_materials().values() + else: + mats = self.materials + return {mat.id: mat for mat in mats} + + @property + @lru_cache(maxsize=None) + def _cells_by_id(self): + """Dictionary mapping cell ID --> cell""" + cells = self.geometry.get_all_cells() + return {cell.id: cell for cell in cells.values()} + + @property + @lru_cache(maxsize=None) + def _cells_by_name(self): + # Get the names maps, but since names are not unique, store a set for + # each name key. In this way when the user requests a change by a name, + # the change will be applied to all of the same name. + result = {} + for cell in self.geometry.get_all_cells().values(): + if cell.name not in result: + result[cell.name] = set() + result[cell.name].add(cell) + return result + + @property + @lru_cache(maxsize=None) + def _materials_by_name(self): + if self.materials is None: + mats = self.geometry.get_all_materials().values() + else: + mats = self.materials + result = {} + for mat in mats: + if mat.name not in result: + result[mat.name] = set() + result[mat.name].add(mat) + return result + @geometry.setter def geometry(self, geometry): check_type('geometry', geometry, openmc.Geometry)