From f207d4220afa0582942b2b8f5281e3151ebac429 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 24 Jan 2025 16:49:58 -0600 Subject: [PATCH] Adding methods to automatically apply results to existing Tally objects. (#2671) Co-authored-by: Paul Romano --- openmc/arithmetic.py | 52 +++---------- openmc/cell.py | 3 +- openmc/mgxs/mgxs.py | 2 +- openmc/model/model.py | 24 +++++- openmc/statepoint.py | 40 ++++++---- openmc/tallies.py | 104 ++++++++++++++++++++------ openmc/universe.py | 3 +- tests/regression_tests/source/test.py | 2 +- tests/unit_tests/test_tallies.py | 77 ++++++++++++++++++- 9 files changed, 223 insertions(+), 84 deletions(-) diff --git a/openmc/arithmetic.py b/openmc/arithmetic.py index 821014e36..92c42284c 100644 --- a/openmc/arithmetic.py +++ b/openmc/arithmetic.py @@ -93,9 +93,9 @@ class CrossNuclide: Parameters ---------- - left_nuclide : openmc.Nuclide or CrossNuclide + left_nuclide : str or CrossNuclide The left nuclide in the outer product - right_nuclide : openmc.Nuclide or CrossNuclide + right_nuclide : str or CrossNuclide The right nuclide in the outer product binary_op : str The tally arithmetic binary operator (e.g., '+', '-', etc.) used to @@ -103,9 +103,9 @@ class CrossNuclide: Attributes ---------- - left_nuclide : openmc.Nuclide or CrossNuclide + left_nuclide : str or CrossNuclide The left nuclide in the outer product - right_nuclide : openmc.Nuclide or CrossNuclide + right_nuclide : str or CrossNuclide The right nuclide in the outer product binary_op : str The tally arithmetic binary operator (e.g., '+', '-', etc.) used to @@ -134,7 +134,7 @@ class CrossNuclide: @left_nuclide.setter def left_nuclide(self, left_nuclide): cv.check_type('left_nuclide', left_nuclide, - (openmc.Nuclide, CrossNuclide, AggregateNuclide)) + (str, CrossNuclide, AggregateNuclide)) self._left_nuclide = left_nuclide @property @@ -144,7 +144,7 @@ class CrossNuclide: @right_nuclide.setter def right_nuclide(self, right_nuclide): cv.check_type('right_nuclide', right_nuclide, - (openmc.Nuclide, CrossNuclide, AggregateNuclide)) + (str, CrossNuclide, AggregateNuclide)) self._right_nuclide = right_nuclide @property @@ -159,26 +159,7 @@ class CrossNuclide: @property def name(self): - - string = '' - - # If the Summary was linked, the left nuclide is a Nuclide object - if isinstance(self.left_nuclide, openmc.Nuclide): - string += '(' + self.left_nuclide.name - # If the Summary was not linked, the left nuclide is the ZAID - else: - string += '(' + str(self.left_nuclide) - - string += ' ' + self.binary_op + ' ' - - # If the Summary was linked, the right nuclide is a Nuclide object - if isinstance(self.right_nuclide, openmc.Nuclide): - string += self.right_nuclide.name + ')' - # If the Summary was not linked, the right nuclide is the ZAID - else: - string += str(self.right_nuclide) + ')' - - return string + return f'({self.left_nuclide} {self.binary_op} {self.right_nuclide})' class CrossFilter: @@ -440,7 +421,7 @@ class AggregateNuclide: Parameters ---------- - nuclides : Iterable of str or openmc.Nuclide or CrossNuclide + nuclides : Iterable of str or CrossNuclide The nuclides included in the aggregation aggregate_op : str The tally aggregation operator (e.g., 'sum', 'avg', etc.) used @@ -448,7 +429,7 @@ class AggregateNuclide: Attributes ---------- - nuclides : Iterable of str or openmc.Nuclide or CrossNuclide + nuclides : Iterable of str or CrossNuclide The nuclides included in the aggregation aggregate_op : str The tally aggregation operator (e.g., 'sum', 'avg', etc.) used @@ -473,13 +454,7 @@ class AggregateNuclide: return str(other) == str(self) def __repr__(self): - - # Append each nuclide in the aggregate to the string - string = f'{self.aggregate_op}(' - names = [nuclide.name if isinstance(nuclide, openmc.Nuclide) - else str(nuclide) for nuclide in self.nuclides] - string += ', '.join(map(str, names)) + ')' - return string + return f'{self.aggregate_op}{self.name}' @property def nuclides(self): @@ -502,12 +477,9 @@ class AggregateNuclide: @property def name(self): - # Append each nuclide in the aggregate to the string - names = [nuclide.name if isinstance(nuclide, openmc.Nuclide) - else str(nuclide) for nuclide in self.nuclides] - string = '(' + ', '.join(map(str, names)) + ')' - return string + names = [str(nuclide) for nuclide in self.nuclides] + return '(' + ', '.join(map(str, names)) + ')' class AggregateFilter: diff --git a/openmc/cell.py b/openmc/cell.py index fe3939bbe..2e0f2ff13 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -403,9 +403,8 @@ class Cell(IDManagerMixin): if self._atoms is not None: volume = self.volume for name, atoms in self._atoms.items(): - nuclide = openmc.Nuclide(name) density = 1.0e-24 * atoms.n/volume # density in atoms/b-cm - nuclides[name] = (nuclide, density) + nuclides[name] = (name, density) else: raise RuntimeError( 'Volume information is needed to calculate microscopic ' diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 588b6f64e..1bc7f9387 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -956,7 +956,7 @@ class MGXS: self.xs_tally._nuclides = [] nuclides = self.get_nuclides() for nuclide in nuclides: - self.xs_tally.nuclides.append(openmc.Nuclide(nuclide)) + self.xs_tally.nuclides.append(nuclide) # Remove NaNs which may have resulted from divide-by-zero operations self.xs_tally._mean = np.nan_to_num(self.xs_tally.mean) diff --git a/openmc/model/model.py b/openmc/model/model.py index 1c261061a..80859c55d 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -14,7 +14,7 @@ import openmc import openmc._xml as xml from openmc.dummy_comm import DummyCommunicator from openmc.executor import _process_CLI_arguments -from openmc.checkvalue import check_type, check_value +from openmc.checkvalue import check_type, check_value, PathLike from openmc.exceptions import InvalidIDError from openmc.utility_funcs import change_directory @@ -604,7 +604,8 @@ class Model: def run(self, particles=None, threads=None, geometry_debug=False, restart_file=None, tracks=False, output=True, cwd='.', openmc_exec='openmc', mpi_args=None, event_based=None, - export_model_xml=True, **export_kwargs): + export_model_xml=True, apply_tally_results=False, + **export_kwargs): """Run OpenMC If the C API has been initialized, then the C API is used, otherwise, @@ -654,6 +655,11 @@ class Model: to True. .. versionadded:: 0.13.3 + apply_tally_results : bool + Whether to apply results of the final statepoint file to the + model's tally objects. + + .. versionadded:: 0.15.1 **export_kwargs Keyword arguments passed to either :meth:`Model.export_to_model_xml` or :meth:`Model.export_to_xml`. @@ -725,6 +731,10 @@ class Model: if mtime >= tstart: # >= allows for poor clock resolution tstart = mtime last_statepoint = sp + + if apply_tally_results: + self.apply_tally_results(last_statepoint) + return last_statepoint def calculate_volumes(self, threads=None, output=True, cwd='.', @@ -932,6 +942,16 @@ class Model: n_samples=n_samples, prn_seed=prn_seed ) + def apply_tally_results(self, statepoint: PathLike | openmc.StatePoint): + """Apply results from a statepoint to tally objects on the Model + + Parameters + ---------- + statepoint : PathLike or openmc.StatePoint + Statepoint file used to update tally results + """ + self.tallies.add_results(statepoint) + def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc', export_model_xml=True, **export_kwargs): """Creates plot images as specified by the Model.plots attribute diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 13aac0caf..5720f8cd0 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -448,14 +448,11 @@ class StatePoint: nuclide_names = group['nuclides'][()] # Add all nuclides to the Tally - for name in nuclide_names: - nuclide = openmc.Nuclide(name.decode().strip()) - tally.nuclides.append(nuclide) + tally.nuclides = [name.decode().strip() for name in nuclide_names] # Add the scores to the Tally scores = group['score_bins'][()] - for score in scores: - tally.scores.append(score.decode()) + tally.scores = [score.decode() for score in scores] # Add Tally to the global dictionary of all Tallies tally.sparse = self.sparse @@ -526,15 +523,16 @@ class StatePoint: def get_tally(self, scores=[], filters=[], nuclides=[], name=None, id=None, estimator=None, exact_filters=False, - exact_nuclides=False, exact_scores=False): + exact_nuclides=False, exact_scores=False, + multiply_density=None, derivative=None): """Finds and returns a Tally object with certain properties. This routine searches the list of Tallies and returns the first Tally found which satisfies all of the input parameters. NOTE: If any of the "exact" parameters are False (default), the input - parameters do not need to match the complete Tally specification and - may only represent a subset of the Tally's properties. If an "exact" + parameters do not need to match the complete Tally specification and may + only represent a subset of the Tally's properties. If an "exact" parameter is True then number of scores, filters, or nuclides in the parameters must precisely match those of any matching Tally. @@ -561,9 +559,15 @@ class StatePoint: to those in the matching Tally. If False (default), the nuclides in the parameters may be a subset of those in the matching Tally. exact_scores : bool - If True, the number of scores in the parameters must be identical - to those in the matching Tally. If False (default), the scores - in the parameters may be a subset of those in the matching Tally. + If True, the number of scores in the parameters must be identical to + those in the matching Tally. If False (default), the scores in the + parameters may be a subset of those in the matching Tally. Default + is None (no check). + multiply_density : bool, optional + If not None, the Tally must have the multiply_density attribute set + to the same value as this parameter. + derivative : openmc.TallyDerivative, optional + TallyDerivative object to match. Returns ------- @@ -591,17 +595,25 @@ class StatePoint: if id and id != test_tally.id: continue - # Determine if Tally has queried estimator - if estimator and estimator != test_tally.estimator: + # Determine if Tally has queried estimator, only move on to next tally + # if the estimator is both specified and the tally estimtor does not + # match + if estimator is not None and estimator != test_tally.estimator: continue # The number of filters, nuclides and scores must exactly match if exact_scores and len(scores) != test_tally.num_scores: continue - if exact_nuclides and len(nuclides) != test_tally.num_nuclides: + if exact_nuclides and nuclides and len(nuclides) != test_tally.num_nuclides: + continue + if exact_nuclides and not nuclides and test_tally.nuclides != ['total']: continue if exact_filters and len(filters) != test_tally.num_filters: continue + if derivative is not None and derivative != test_tally.derivative: + continue + if multiply_density is not None and multiply_density != test_tally.multiply_density: + continue # Determine if Tally has the queried score(s) if scores: diff --git a/openmc/tallies.py b/openmc/tallies.py index f39ffa078..b0d78d357 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1,3 +1,4 @@ +from __future__ import annotations from collections.abc import Iterable, MutableSequence import copy from functools import partial, reduce @@ -33,7 +34,7 @@ _NUCLIDE_CLASSES = (str, openmc.CrossNuclide, openmc.AggregateNuclide) _FILTER_CLASSES = (openmc.Filter, openmc.CrossFilter, openmc.AggregateFilter) # Valid types of estimators -ESTIMATOR_TYPES = ['tracklength', 'collision', 'analog'] +ESTIMATOR_TYPES = {'tracklength', 'collision', 'analog'} class Tally(IDManagerMixin): @@ -65,7 +66,9 @@ class Tally(IDManagerMixin): scores : list of str List of defined scores, e.g. 'flux', 'fission', etc. estimator : {'analog', 'tracklength', 'collision'} - Type of estimator for the tally + Type of estimator for the tally. If unset (None), OpenMC will automatically + select an appropriate estimator based on the tally filters and scores + with a preference for 'tracklength'. triggers : list of openmc.Trigger List of tally triggers num_scores : int @@ -131,6 +134,36 @@ class Tally(IDManagerMixin): self._sp_filename = None self._results_read = False + def __eq__(self, other): + if other.id != self.id: + return False + if other.name != self.name: + return False + # estimators are automatically set based on the tally filters and scores + # during OpenMC initialization if this value is None, so it is not + # considered a requirement for equivalence if it is unset on either + # tally as it implies that the user is allowing OpenMC to select an + # appropriate estimator. If the value is explicitly set on both tallies, + # then the values must match for the tallies to be considered equivalent. + if self.estimator is not None and other.estimator is not None and other.estimator != self.estimator: + return False + if other.filters != self.filters: + return False + # for tallies are loaded from statpoint files + # an empty nuclide list is equivalent to a list with 'total' + other_nuclides = other.nuclides.copy() + self_nuclides = self.nuclides.copy() + if 'total' in other_nuclides: + other_nuclides.remove('total') + if 'total' in self_nuclides: + self_nuclides.remove('total') + if other_nuclides != self_nuclides: + return False + for attr in {'scores', 'triggers', 'derivative', 'multiply_density'}: + if getattr(other, attr) != getattr(self, attr): + return False + return True + def __repr__(self): parts = ['Tally'] parts.append('{: <15}=\t{}'.format('ID', self.id)) @@ -266,7 +299,8 @@ class Tally(IDManagerMixin): @estimator.setter def estimator(self, estimator): - cv.check_value('estimator', estimator, ESTIMATOR_TYPES) + # allow the estimator to be set to None (let OpenMC choose the estimator at runtime) + cv.check_value('estimator', estimator, ESTIMATOR_TYPES | {None}) self._estimator = estimator @property @@ -304,8 +338,16 @@ class Tally(IDManagerMixin): # Open the HDF5 statepoint file with h5py.File(self._sp_filename, 'r') as f: + # Set number of realizations + group = f[f'tallies/tally {self.id}'] + self.num_realizations = int(group['n_realizations'][()]) + + # Update nuclides + nuclide_names = group['nuclides'][()] + self.nuclides = [name.decode().strip() for name in nuclide_names] + # Extract Tally data from the file - data = f[f'tallies/tally {self.id}/results'] + data = group['results'] sum_ = data[:, :, 0] sum_sq = data[:, :, 1] @@ -511,7 +553,7 @@ class Tally(IDManagerMixin): Parameters ---------- - nuclide : openmc.Nuclide + nuclide : str Nuclide to remove """ @@ -883,6 +925,21 @@ class Tally(IDManagerMixin): return element + def add_results(self, statepoint: cv.PathLike | openmc.StatePoint): + """Add results from the provided statepoint file to this tally instance + + .. versionadded:: 0.15.1 + + Parameters + ---------- + statepoint : openmc.PathLike or openmc.StatePoint + Statepoint used to update tally results + """ + if isinstance(statepoint, openmc.StatePoint): + self._sp_filename = statepoint._f.filename + else: + self._sp_filename = str(statepoint) + @classmethod def from_xml_element(cls, elem, **kwargs): """Generate tally object from an XML element @@ -1020,17 +1077,10 @@ class Tally(IDManagerMixin): in the Tally. """ - # Look for the user-requested nuclide in all of the Tally's Nuclides + # Look for the user-requested nuclide in all of the Tally's nuclides for i, test_nuclide in enumerate(self.nuclides): - # If the Summary was linked, then values are Nuclide objects - if isinstance(test_nuclide, openmc.Nuclide): - if test_nuclide.name == nuclide: - return i - - # If the Summary has not been linked, then values are ZAIDs - else: - if test_nuclide == nuclide: - return i + if test_nuclide == nuclide: + return i msg = (f'Unable to get the nuclide index for Tally since "{nuclide}" ' 'is not one of the nuclides') @@ -1357,9 +1407,7 @@ class Tally(IDManagerMixin): column_name = 'nuclide' for nuclide in self.nuclides: - if isinstance(nuclide, openmc.Nuclide): - nuclides.append(nuclide.name) - elif isinstance(nuclide, openmc.AggregateNuclide): + if isinstance(nuclide, openmc.AggregateNuclide): nuclides.append(nuclide.name) column_name = f'{nuclide.aggregate_op}(nuclide)' else: @@ -2350,7 +2398,8 @@ class Tally(IDManagerMixin): new_tally.name = self.name new_tally._mean = self.mean / other new_tally._std_dev = self.std_dev * np.abs(1. / other) - new_tally.estimator = self.estimator + if self.estimator is not None: + new_tally.estimator = self.estimator new_tally.with_summary = self.with_summary new_tally.num_realizations = self.num_realizations @@ -2643,8 +2692,8 @@ class Tally(IDManagerMixin): # Determine the nuclide indices from any of the requested nuclides for nuclide in self.nuclides: - if nuclide.name not in nuclides: - nuclide_index = self.get_nuclide_index(nuclide.name) + if nuclide not in nuclides: + nuclide_index = self.get_nuclide_index(nuclide) nuclide_indices.append(nuclide_index) # Loop over indices in reverse to remove excluded Nuclides @@ -3166,6 +3215,19 @@ class Tallies(cv.CheckedList): # Continue iterating from the first loop break + def add_results(self, statepoint: cv.PathLike | openmc.StatePoint): + """Add results from the provided statepoint file + + .. versionadded:: 0.15.1 + + Parameters + ---------- + statepoint : openmc.PathLike or openmc.StatePoint + Statepoint used to update tally results + """ + for tally in self: + tally.add_results(statepoint) + def _create_tally_subelements(self, root_element): for tally in self: root_element.append(tally.to_xml_element()) diff --git a/openmc/universe.py b/openmc/universe.py index 85ce6fd96..71ffa726b 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -598,9 +598,8 @@ class UniverseBase(ABC, IDManagerMixin): if self._atoms: volume = self.volume for name, atoms in self._atoms.items(): - nuclide = openmc.Nuclide(name) density = 1.0e-24 * atoms.n/volume # density in atoms/b-cm - nuclides[name] = (nuclide, density) + nuclides[name] = (name, density) else: raise RuntimeError( 'Volume information is needed to calculate microscopic cross ' diff --git a/tests/regression_tests/source/test.py b/tests/regression_tests/source/test.py index c1f4dccff..96efa77eb 100644 --- a/tests/regression_tests/source/test.py +++ b/tests/regression_tests/source/test.py @@ -11,7 +11,7 @@ class SourceTestHarness(PyAPITestHarness): super().__init__(*args, **kwargs) mat1 = openmc.Material(material_id=1, temperature=294) mat1.set_density('g/cm3', 4.5) - mat1.add_nuclide(openmc.Nuclide('U235'), 1.0) + mat1.add_nuclide('U235', 1.0) self._model.materials = openmc.Materials([mat1]) sphere = openmc.Sphere(surface_id=1, r=10.0, boundary_type='vacuum') diff --git a/tests/unit_tests/test_tallies.py b/tests/unit_tests/test_tallies.py index 544443312..f044df5d0 100644 --- a/tests/unit_tests/test_tallies.py +++ b/tests/unit_tests/test_tallies.py @@ -41,4 +41,79 @@ def test_xml_roundtrip(run_in_tmpdir): assert len(new_tally.triggers) == 1 assert new_tally.triggers[0].trigger_type == tally.triggers[0].trigger_type assert new_tally.triggers[0].threshold == tally.triggers[0].threshold - assert new_tally.triggers[0].scores == tally.triggers[0].scores \ No newline at end of file + assert new_tally.triggers[0].scores == tally.triggers[0].scores + assert new_tally.multiply_density == tally.multiply_density + + +def test_tally_equivalence(): + tally_a = openmc.Tally() + tally_b = openmc.Tally(tally_id=tally_a.id) + + tally_a.name = 'new name' + assert tally_a != tally_b + tally_b.name = tally_a.name + assert tally_a == tally_b + + assert tally_a == tally_b + ef_a = openmc.EnergyFilter([0.0, 0.1, 1.0, 10.0e6]) + ef_b = openmc.EnergyFilter([0.0, 0.1, 1.0, 10.0e6]) + + tally_a.filters = [ef_a] + assert tally_a != tally_b + tally_b.filters = [ef_b] + assert tally_a == tally_b + + tally_a.scores = ['flux', 'absorption', 'fission', 'scatter'] + assert tally_a != tally_b + tally_b.scores = ['flux', 'absorption', 'fission', 'scatter'] + assert tally_a == tally_b + + tally_a.nuclides = [] + tally_b.nuclides = [] + assert tally_a == tally_b + + tally_a.nuclides = ['total'] + assert tally_a == tally_b + + # a tally with an estimator set to None is equal to + # a tally with an estimator specified + tally_a.estimator = 'collision' + assert tally_a == tally_b + tally_b.estimator = 'collision' + assert tally_a == tally_b + + tally_a.multiply_density = False + assert tally_a != tally_b + tally_b.multiply_density = False + assert tally_a == tally_b + + trigger_a = openmc.Trigger('rel_err', 0.025) + trigger_b = openmc.Trigger('rel_err', 0.025) + + tally_a.triggers = [trigger_a] + assert tally_a != tally_b + tally_b.triggers = [trigger_b] + assert tally_a == tally_b + + +def test_tally_application(sphere_model, run_in_tmpdir): + # Create a tally with most possible gizmos + tally = openmc.Tally(name='test tally') + ef = openmc.EnergyFilter([0.0, 0.1, 1.0, 10.0e6]) + mesh = openmc.RegularMesh.from_domain(sphere_model.geometry, (2, 2, 2)) + mf = openmc.MeshFilter(mesh) + tally.filters = [ef, mf] + tally.scores = ['flux', 'absorption', 'fission', 'scatter'] + sphere_model.tallies = [tally] + + # run the simulation and apply results + sp_file = sphere_model.run(apply_tally_results=True) + with openmc.StatePoint(sp_file) as sp: + assert tally in sp.tallies.values() + sp_tally = sp.tallies[tally.id] + + # at this point the tally information regarding results should be the same + assert (sp_tally.mean == tally.mean).all() + assert (sp_tally.std_dev == tally.std_dev).all() + assert sp_tally.nuclides == tally.nuclides +