Fixed bugs in tally sparsification which broke MGXS library

This commit is contained in:
wbinventor@gmail.com 2015-12-13 11:14:01 -05:00
parent e85b659977
commit beb3189d8d
3 changed files with 197 additions and 202 deletions

View file

@ -73,6 +73,9 @@ class Library(object):
name : str, optional
Name of the multi-group cross section library. Used as a label to
identify tallies in OpenMC 'tallies.xml' file.
sparse : bool
Whether or not the Library's tallies use SciPy's LIL sparse matrix
format for compressed data storage
"""
@ -92,6 +95,7 @@ class Library(object):
self._all_mgxs = OrderedDict()
self._sp_filename = None
self._keff = None
self._sparse = False
self.name = name
self.openmc_geometry = openmc_geometry
@ -119,6 +123,7 @@ class Library(object):
clone._all_mgxs = self.all_mgxs
clone._sp_filename = self._sp_filename
clone._keff = self._keff
clone._sparse = self.sparse
clone._all_mgxs = OrderedDict()
for domain in self.domains:
@ -208,6 +213,10 @@ class Library(object):
def keff(self):
return self._keff
@property
def sparse(self):
return self._sparse
@openmc_geometry.setter
def openmc_geometry(self, openmc_geometry):
cv.check_type('openmc_geometry', openmc_geometry, openmc.Geometry)
@ -286,6 +295,28 @@ class Library(object):
cv.check_type('tally trigger', tally_trigger, openmc.Trigger)
self._tally_trigger = tally_trigger
@sparse.setter
def sparse(self, sparse):
"""Convert tally data from NumPy arrays to SciPy list of lists (LIL)
sparse matrices, and vice versa.
This property may be used to reduce the amount of data in memory during
tally data processing. The tally data will be stored as SciPy LIL
matrices internally within the Tally object. All tally data access
properties and methods will return data as a dense NumPy array.
"""
cv.check_type('sparse', sparse, bool)
# Sparsify or densify each MGXS in the Library
for domain in self.domains:
for mgxs_type in self.mgxs_types:
mgxs = self.get_mgxs(domain, mgxs_type)
mgxs.sparse = self.sparse
self._sparse = sparse
def build_library(self):
"""Initialize MGXS objects in each domain and for each reaction type
in the library.
@ -381,6 +412,7 @@ class Library(object):
for mgxs_type in self.mgxs_types:
mgxs = self.get_mgxs(domain, mgxs_type)
mgxs.load_from_statepoint(statepoint)
mgxs.sparse = self.sparse
def get_mgxs(self, domain, mgxs_type):
"""Return the MGXS object for some domain and reaction rate type.

View file

@ -258,6 +258,29 @@ class MGXS(object):
cv.check_type('tally trigger', tally_trigger, openmc.Trigger)
self._tally_trigger = tally_trigger
@sparse.setter
def sparse(self, sparse):
"""Convert tally data from NumPy arrays to SciPy list of lists (LIL)
sparse matrices, and vice versa.
This property may be used to reduce the amount of data in memory during
tally data processing. The tally data will be stored as SciPy LIL
matrices internally within the Tally object. All tally data access
properties and methods will return data as a dense NumPy array.
"""
cv.check_type('sparse', sparse, bool)
# Sparsify or densify the derived MGXS tally and its base tallies
if self.xs_tally:
self.xs_tally.sparse = sparse
for tally_name in self.tallies:
self.tallies[tally_name].sparse = sparse
self._sparse = sparse
@staticmethod
def get_mgxs(mgxs_type, domain=None, domain_type=None,
energy_groups=None, by_nuclide=False, name=''):
@ -512,10 +535,7 @@ class MGXS(object):
# Remove NaNs which may have resulted from divide-by-zero operations
self._xs_tally._mean = np.nan_to_num(self.xs_tally.mean)
self._xs_tally._std_dev = np.nan_to_num(self.xs_tally.std_dev)
# Sparsify the MGXS derived tally
if self.sparse:
self.xs_tally.sparsify()
self.xs_tally.sparse = self.sparse
def load_from_statepoint(self, statepoint):
"""Extracts tallies in an OpenMC StatePoint with the data needed to
@ -578,11 +598,9 @@ class MGXS(object):
estimator=tally.estimator)
sp_tally = sp_tally.get_slice(tally.scores, filters,
filter_bins, tally.nuclides)
sp_tally.sparse = self.sparse
self.tallies[tally_type] = sp_tally
if self.sparse:
sp_tally.sparsify()
# Compute the cross section from the tallies
self.compute_xs()
@ -732,6 +750,7 @@ class MGXS(object):
# Clone this MGXS to initialize the condensed version
condensed_xs = copy.deepcopy(self)
condensed_xs.sparse = False
condensed_xs.energy_groups = coarse_groups
# Build energy indices to sum across
@ -781,11 +800,7 @@ class MGXS(object):
# Compute the energy condensed multi-group cross section
condensed_xs.compute_xs()
# Sparsify the condensed MGXS
if self.sparse:
condensed_xs.sparsify()
condensed_xs.sparse = self.sparse
return condensed_xs
def get_subdomain_avg_xs(self, subdomains='all'):
@ -828,8 +843,9 @@ class MGXS(object):
# Clone this MGXS to initialize the subdomain-averaged version
avg_xs = copy.deepcopy(self)
avg_xs.sparse = False
# If domain is distribcell, make the new domain 'cell'
# If domain is distribcell, make the new domain 'cell'
if self.domain_type == 'distribcell':
avg_xs.domain_type = 'cell'
@ -868,11 +884,7 @@ class MGXS(object):
# Compute the subdomain-averaged multi-group cross section
avg_xs.compute_xs()
# Sparsify the subdomain-averaged MGXS
if self.sparse:
avg_xs.sparsify()
avg_xs.sparse = self.sparse
return avg_xs
def print_xs(self, subdomains='all', nuclides='all', xs_type='macro'):
@ -2103,59 +2115,6 @@ class Chi(MGXS):
super(Chi, self).compute_xs()
def sparsify(self):
"""Convert tally data from NumPy arrays to SciPy list of lists (LIL)
sparse matrices.
This method may be used to reduce the amount of data in memory during
tally data processing. The tally data will be stored as SciPy LIL
matrices internally within the Tally object. All tally data access
properties and methods will return data as a dense NumPy array.
See also
--------
MGXS.densify(), Tally.sparsify(), Tally.densify()
"""
# Sparsify the derived MGXS tally
if self.xs_tally:
self.xs_tally.sparsify()
# Sparsify the MGXS' base tallies
for tally_name in self.tallies:
self.tallies[tally_name].sparsify()
self._sparse = True
def densify(self):
"""Convert tally data from SciPy list of lists (LIL) sparse matrices
to NumPy arrrays.
This method may be used to restore a sparse tally to its original
state. This method will have no effect on the state of the tally
unless the Tally.sparse() method has been called.
See also
--------
MGXS.sparsify(), Tally.sparsify(), Tally.densify()
"""
# If the MGXS is already dense, simply return
if not self.sparse:
return
# Densify the derived MGXS tally
if self.xs_tally:
self.xs_tally.densify()
# Densify the MGXS' base tallies
for tally_name in self.tallies:
self.tallies[tally_name].densify()
self._sparse = False
def get_xs(self, groups='all', subdomains='all', nuclides='all',
xs_type='macro', order_groups='increasing', value='mean'):
"""Returns an array of the fission spectrum.

View file

@ -326,7 +326,8 @@ class Tally(object):
# Convert NumPy arrays to SciPy sparse LIL matrices
if self.sparse:
self._sum = sps.lil_matrix(self._sum.flatten(), self._sum.shape)
self._sum = \
sps.lil_matrix(self._sum.flatten(), self._sum.shape)
self._sum_sq = \
sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape)
@ -536,6 +537,48 @@ class Tally(object):
cv.check_type('sum_sq', sum_sq, Iterable)
self._sum_sq = sum_sq
@sparse.setter
def sparse(self, sparse):
"""Convert tally data from NumPy arrays to SciPy list of lists (LIL)
sparse matrices, and vice versa.
This property may be used to reduce the amount of data in memory during
tally data processing. The tally data will be stored as SciPy LIL
matrices internally within the Tally object. All tally data access
properties and methods will return data as a dense NumPy array.
"""
cv.check_type('sparse', sparse, bool)
# Convert NumPy arrays to SciPy sparse LIL matrices
if sparse and not self.sparse:
if self._sum is not None:
self._sum = \
sps.lil_matrix(self._sum.flatten(), self._sum.shape)
if self._sum_sq is not None:
self._sum_sq = \
sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape)
if self._mean is not None:
self._mean = \
sps.lil_matrix(self._mean.flatten(), self._mean.shape)
if self._std_dev is not None:
self._std_dev = \
sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape)
self._sparse = True
# Convert SciPy sparse LIL matrices to NumPy arrays
elif not sparse and self.sparse:
if self._sum is not None:
self._sum = np.reshape(self._sum.toarray(), self.shape)
if self._sum_sq is not None:
self._sum_sq = np.reshape(self._sum_sq.toarray(), self.shape)
if self._mean is not None:
self._mean = np.reshape(self._mean.toarray(), self.shape)
if self._std_dev is not None:
self._std_dev = np.reshape(self._std_dev.toarray(), self.shape)
self._sparse = False
def remove_score(self, score):
"""Remove a score from the tally
@ -1149,67 +1192,6 @@ class Tally(object):
return data
def sparsify(self):
"""Convert tally data from NumPy arrays to SciPy list of lists (LIL)
sparse matrices.
This method may be used to reduce the amount of data in memory during
tally data processing. The tally data will be stored as SciPy LIL
matrices internally within the Tally object. All tally data access
properties and methods will return data as a dense NumPy array.
See also
--------
Tally.densify()
"""
# Convert NumPy arrays to SciPy sparse LIL matrices
if self._sum is not None:
self._sum = \
sps.lil_matrix(self._sum.flatten(), self._sum.shape)
if self._sum_sq is not None:
self._sum_sq = \
sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape)
if self._mean is not None:
self._mean = \
sps.lil_matrix(self._mean.flatten(), self._mean.shape)
if self._std_dev is not None:
self._std_dev = \
sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape)
self._sparse = True
def densify(self):
"""Convert tally data from SciPy list of lists (LIL) sparse matrices
to NumPy arrrays.
This method may be used to restore a sparse tally to its original
state. This method will have no effect on the state of the tally
unless the Tally.sparse() method has been called.
See also
--------
Tally.sparsify()
"""
# If the tally is already dense, simply return
if not self.sparse:
return
# Convert SciPy sparse LIL matrices to NumPy arrays
if self._sum is not None:
self._sum = np.reshape(self._sum.toarray(), self.shape)
if self._sum_sq is not None:
self._sum_sq = np.reshape(self._sum_sq.toarray(), self.shape)
if self._mean is not None:
self._mean = np.reshape(self._mean.toarray(), self.shape)
if self._std_dev is not None:
self._std_dev = np.reshape(self._std_dev.toarray(), self.shape)
self._sparse = False
def get_pandas_dataframe(self, filters=True, nuclides=True,
scores=True, summary=None):
"""Build a Pandas DataFrame for the Tally data.
@ -1626,6 +1608,9 @@ class Tally(object):
self_copy = copy.deepcopy(self)
other_copy = copy.deepcopy(other)
self_copy.sparse = False
other_copy.sparse = False
# Align the tally data based on desired hybrid product
data = self_copy._align_tally_data(other_copy, filter_product,
nuclide_product, score_product)
@ -1691,7 +1676,8 @@ class Tally(object):
else:
all_nuclides = [self_copy.nuclides, other_copy.nuclides]
for self_nuclide, other_nuclide in itertools.product(*all_nuclides):
new_nuclide = CrossNuclide(self_nuclide, other_nuclide, binary_op)
new_nuclide = \
CrossNuclide(self_nuclide, other_nuclide, binary_op)
new_tally.add_nuclide(new_nuclide)
# Add scores to the new tally
@ -1700,7 +1686,8 @@ class Tally(object):
for self_score in self_copy.scores:
new_tally.add_score(self_score)
else:
new_tally.num_score_bins = self_copy.num_score_bins * other_copy.num_score_bins
new_tally.num_score_bins = \
self_copy.num_score_bins * other_copy.num_score_bins
all_scores = [self_copy.scores, other_copy.scores]
for self_score, other_score in itertools.product(*all_scores):
new_score = CrossScore(self_score, other_score, binary_op)
@ -1717,7 +1704,6 @@ class Tally(object):
def _align_tally_data(self, other, filter_product, nuclide_product,
score_product):
"""Aligns data from two tallies for tally arithmetic.
This is a helper method to construct a dict of dicts of the "aligned"
data arrays from each tally for tally arithmetic. The method analyzes
the filters, scores and nuclides in both tallies and determines how to
@ -1726,7 +1712,6 @@ class Tally(object):
'tile' and 'repeat' operations to the new data arrays such that all
possible combinations of the data in each tally's bins will be made
when the arithmetic operation is applied to the arrays.
Parameters
----------
other : Tally
@ -1740,24 +1725,18 @@ class Tally(object):
score_product : str
The type of product (tensor or entrywise) to be performed between
score data.
Returns
-------
dict
A dictionary of dictionaries to "aligned" 'mean' and 'std. dev'
NumPy arrays for each tally's data.
"""
# Use dense NumPy arrays for data alignment operations
self_sparse = self.sparse
other_sparse = other.sparse
self.densify()
other.densify()
# Get the set of filters that each tally is missing
other_missing_filters = set(self.filters).difference(set(other.filters))
self_missing_filters = set(other.filters).difference(set(self.filters))
other_missing_filters = \
set(self.filters).difference(set(other.filters))
self_missing_filters = \
set(other.filters).difference(set(self.filters))
# Add filters present in self but not in other to other
for filter in other_missing_filters:
@ -1784,30 +1763,40 @@ class Tally(object):
# Repeat and tile the data by nuclide in preparation for performing
# the tensor product across nuclides.
if nuclide_product == 'tensor':
self._mean = np.repeat(self.mean, other.num_nuclides, axis=1)
self._std_dev = np.repeat(self.std_dev, other.num_nuclides, axis=1)
other._mean = np.tile(other.mean, (1, self.num_nuclides, 1))
other._std_dev = np.tile(other.std_dev, (1, self.num_nuclides, 1))
self._mean = \
np.repeat(self.mean, other.num_nuclides, axis=1)
self._std_dev = \
np.repeat(self.std_dev, other.num_nuclides, axis=1)
other._mean = \
np.tile(other.mean, (1, self.num_nuclides, 1))
other._std_dev = \
np.tile(other.std_dev, (1, self.num_nuclides, 1))
# Add nuclides to each tally such that each tally contains the complete
# set of nuclides necessary to perform an entrywise product. New nuclides
# added to a tally will have all their scores set to zero.
# set of nuclides necessary to perform an entrywise product. New
# nuclides added to a tally will have all their scores set to zero.
else:
# Get the set of nuclides that each tally is missing
other_missing_nuclides = set(self.nuclides).difference(set(other.nuclides))
self_missing_nuclides = set(other.nuclides).difference(set(self.nuclides))
other_missing_nuclides = \
set(self.nuclides).difference(set(other.nuclides))
self_missing_nuclides = \
set(other.nuclides).difference(set(self.nuclides))
# Add nuclides present in self but not in other to other
for nuclide in other_missing_nuclides:
other._mean = np.insert(other.mean, other.num_nuclides, 0, axis=1)
other._std_dev = np.insert(other.std_dev, other.num_nuclides, 0, axis=1)
other._mean = \
np.insert(other.mean, other.num_nuclides, 0, axis=1)
other._std_dev = \
np.insert(other.std_dev, other.num_nuclides, 0, axis=1)
other.add_nuclide(nuclide)
# Add nuclides present in other but not in self to self
for nuclide in self_missing_nuclides:
self._mean = np.insert(self.mean, self.num_nuclides, 0, axis=1)
self._std_dev = np.insert(self.std_dev, self.num_nuclides, 0, axis=1)
self._mean = \
np.insert(self.mean, self.num_nuclides, 0, axis=1)
self._std_dev = \
np.insert(self.std_dev, self.num_nuclides, 0, axis=1)
self.add_nuclide(nuclide)
# Align other nuclides with self nuclides
@ -1821,30 +1810,40 @@ class Tally(object):
# Repeat and tile the data by score in preparation for performing
# the tensor product across scores.
if score_product == 'tensor':
self._mean = np.repeat(self.mean, other.num_score_bins, axis=2)
self._std_dev = np.repeat(self.std_dev, other.num_score_bins, axis=2)
other._mean = np.tile(other.mean, (1, 1, self.num_score_bins))
other._std_dev = np.tile(other.std_dev, (1, 1, self.num_score_bins))
self._mean = \
np.repeat(self.mean, other.num_score_bins, axis=2)
self._std_dev = \
np.repeat(self.std_dev, other.num_score_bins, axis=2)
other._mean = \
np.tile(other.mean, (1, 1, self.num_score_bins))
other._std_dev = \
np.tile(other.std_dev, (1, 1, self.num_score_bins))
# Add scores to each tally such that each tally contains the complete set
# of scores necessary to perform an entrywise product. New scores added
# to a tally will be set to zero.
# Add scores to each tally such that each tally contains the complete
# set of scores necessary to perform an entrywise product. New scores
# added to a tally will be set to zero.
else:
# Get the set of scores that each tally is missing
other_missing_scores = set(self.scores).difference(set(other.scores))
self_missing_scores = set(other.scores).difference(set(self.scores))
other_missing_scores = \
set(self.scores).difference(set(other.scores))
self_missing_scores = \
set(other.scores).difference(set(self.scores))
# Add scores present in self but not in other to other
for score in other_missing_scores:
other._mean = np.insert(other.mean, other.num_score_bins, 0, axis=2)
other._std_dev = np.insert(other.std_dev, other.num_score_bins, 0, axis=2)
other._mean = \
np.insert(other.mean, other.num_score_bins, 0, axis=2)
other._std_dev = \
np.insert(other.std_dev, other.num_score_bins, 0, axis=2)
other.add_score(score)
# Add scores present in other but not in self to self
for score in self_missing_scores:
self._mean = np.insert(self.mean, self.num_score_bins, 0, axis=2)
self._std_dev = np.insert(self.std_dev, self.num_score_bins, 0, axis=2)
self._mean = \
np.insert(self.mean, self.num_score_bins, 0, axis=2)
self._std_dev = \
np.insert(self.std_dev, self.num_score_bins, 0, axis=2)
self.add_score(score)
# Align other scores with self scores
@ -1867,12 +1866,6 @@ class Tally(object):
filter.stride = stride
stride *= filter.num_bins
# Restore tally operands to sparse storage
if self.sparse:
self.sparsify()
if other.sparse:
other.sparsify()
# Deep copy the mean and std dev data
self_mean = copy.deepcopy(self.mean)
self_std_dev = copy.deepcopy(self.std_dev)
@ -1886,6 +1879,7 @@ class Tally(object):
data['other']['mean'] = other_mean
data['self']['std. dev.'] = self_std_dev
data['other']['std. dev.'] = other_std_dev
return data
def swap_filters(self, filter1, filter2, inplace=False):
@ -2169,6 +2163,10 @@ class Tally(object):
if isinstance(other, Tally):
new_tally = self.hybrid_product(other, binary_op='+')
# If both tally operands were sparse, sparsify the new tally
if self.sparse and other.sparse:
new_tally.sparse = True
elif isinstance(other, Real):
new_tally = Tally(name='derived')
new_tally._derived = True
@ -2188,9 +2186,8 @@ class Tally(object):
for score in self.scores:
new_tally.add_score(score)
# If original tally operands were sparse, sparsify the sliced tally
if self.sparse and other.sparse:
new_tally.sparsify()
# If this tally operand is sparse, sparsify the new tally
new_tally.sparse = self.sparse
else:
msg = 'Unable to add "{0}" to Tally ID="{1}"'.format(other, self.id)
@ -2242,6 +2239,10 @@ class Tally(object):
if isinstance(other, Tally):
new_tally = self.hybrid_product(other, binary_op='-')
# If both tally operands were sparse, sparsify the new tally
if self.sparse and other.sparse:
new_tally.sparse = True
elif isinstance(other, Real):
new_tally = Tally(name='derived')
new_tally._derived = True
@ -2260,9 +2261,8 @@ class Tally(object):
for score in self.scores:
new_tally.add_score(score)
# If original tally operands were sparse, sparsify the sliced tally
if self.sparse and other.sparse:
new_tally.sparsify()
# If this tally operand is sparse, sparsify the new tally
new_tally.sparse = self.sparse
else:
msg = 'Unable to subtract "{0}" from Tally ' \
@ -2315,6 +2315,10 @@ class Tally(object):
if isinstance(other, Tally):
new_tally = self.hybrid_product(other, binary_op='*')
# If original tally operands were sparse, sparsify the new tally
if self.sparse and other.sparse:
new_tally.sparse = True
elif isinstance(other, Real):
new_tally = Tally(name='derived')
new_tally._derived = True
@ -2333,9 +2337,8 @@ class Tally(object):
for score in self.scores:
new_tally.add_score(score)
# If original tally operands were sparse, sparsify the sliced tally
if self.sparse and other.sparse:
new_tally.sparsify()
# If this tally operand is sparse, sparsify the new tally
new_tally.sparse = self.sparse
else:
msg = 'Unable to multiply Tally ID="{0}" ' \
@ -2388,6 +2391,10 @@ class Tally(object):
if isinstance(other, Tally):
new_tally = self.hybrid_product(other, binary_op='/')
# If original tally operands were sparse, sparsify the new tally
if self.sparse and other.sparse:
new_tally.sparse = True
elif isinstance(other, Real):
new_tally = Tally(name='derived')
new_tally._derived = True
@ -2406,9 +2413,8 @@ class Tally(object):
for score in self.scores:
new_tally.add_score(score)
# If original tally operands were sparse, sparsify the sliced tally
if self.sparse and other.sparse:
new_tally.sparsify()
# If this tally operand is sparse, sparsify the new tally
new_tally.sparse = self.sparse
else:
msg = 'Unable to divide Tally ID="{0}" ' \
@ -2464,6 +2470,10 @@ class Tally(object):
if isinstance(power, Tally):
new_tally = self.hybrid_product(power, binary_op='^')
# If original tally operands were sparse, sparsify the new tally
if self.sparse and other.sparse:
new_tally.sparse = True
elif isinstance(power, Real):
new_tally = Tally(name='derived')
new_tally._derived = True
@ -2484,8 +2494,7 @@ class Tally(object):
new_tally.add_score(score)
# If original tally was sparse, sparsify the exponentiated tally
if self.sparse:
new_tally.sparsify()
new_tally.sparse = self.sparse
else:
msg = 'Unable to raise Tally ID="{0}" to ' \
@ -2648,6 +2657,7 @@ class Tally(object):
raise ValueError(msg)
new_tally = copy.deepcopy(self)
new_tally.sparse = False
if self.sum is not None:
new_sum = self.get_values(scores, filters, filter_bins,
@ -2663,7 +2673,7 @@ class Tally(object):
new_tally._mean = new_mean
if self.std_dev is not None:
new_std_dev = self.get_values(scores, filters, filter_bins,
nuclides, 'std_dev')
nuclides, 'std_dev')
new_tally._std_dev = new_std_dev
# SCORES
@ -2725,9 +2735,7 @@ class Tally(object):
stride *= filter.num_bins
# If original tally was sparse, sparsify the sliced tally
if self.sparse:
new_tally.sparsify()
new_tally.sparse = self.sparse
return new_tally
def summation(self, scores=[], filter_type=None,
@ -2831,9 +2839,7 @@ class Tally(object):
tally_sum.add_filter(filters[-1])
# If original tally was sparse, sparsify the tally summation
if self.sparse:
tally_sum.sparsify()
tally_sum.sparse = self.sparse
return tally_sum
def diagonalize_filter(self, new_filter):
@ -2909,9 +2915,7 @@ class Tally(object):
stride *= filter.num_bins
# If original tally was sparse, sparsify the diagonalized tally
if self.sparse:
new_tally.sparsify
new_tally.sparse = self.sparse
return new_tally