Adding methods to automatically apply results to existing Tally objects. (#2671)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
Patrick Shriwise 2025-01-24 16:49:58 -06:00 committed by GitHub
parent 7089780026
commit f207d4220a
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
9 changed files with 223 additions and 84 deletions

View file

@ -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:

View file

@ -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 '

View file

@ -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)

View file

@ -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

View file

@ -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:

View file

@ -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())

View file

@ -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 '

View file

@ -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')

View file

@ -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
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