crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Determine Tally.sum, sum_sq, mean, std_dev on demand. No more read_results()!

Browse source
This commit is contained in:
Paul Romano 2015-09-16 16:19:12 +07:00
parent aac5eb6f52
commit 97a7582c44
8 changed files with 67 additions and 95 deletions
Download patch file Download diff file Expand all files Collapse all files

47
openmc/statepoint.py
View file

@ -146,7 +146,6 @@ class StatePoint(object):
# Set flags for what data has been read
self._meshes_read = False
self._tallies_read = False
self._results_read = False
self._source_read = False
self._with_summary = False
@ -317,6 +316,8 @@ class StatePoint(object):
# Add mesh to the global dictionary of all Meshes
self._meshes[mesh_id] = mesh
self._meshes_read = True
return self._meshes
@property
@ -401,6 +402,7 @@ class StatePoint(object):
# Create Tally object and assign basic properties
tally = openmc.Tally(tally_key)
tally._statepoint = self
tally.estimator = ESTIMATOR_TYPES[estimator_type]
tally.num_realizations = n_realizations
@ -497,6 +499,8 @@ class StatePoint(object):
# Add Tally to the global dictionary of all Tallies
self._tallies[tally_key] = tally
self._tallies_read = True
return self._tallies
@property
@ -513,47 +517,6 @@ class StatePoint(object):
def with_summary(self):
return self._with_summary
def read_results(self):
"""Read tally results and store them in the ``tallies`` attribute. No results
are read when the statepoint is instantiated.
"""
base = 'tallies/tally '
# Read Tally results
if self.tallies_present:
# Iterate over and extract the results for all Tallies
for tally_key, tally in self.tallies.items():
# Compute the total number of bins for this Tally
num_tot_bins = tally.num_bins
# Extract Tally data from the file
data = self._f['{0}{1}/results'.format(base, tally_key)].value
sum = data['sum']
sum_sq = data['sum_sq']
# Define a routine to convert 0 to 1
def nonzero(val):
return 1 if not val else val
# Reshape the results arrays
new_shape = (nonzero(tally.num_filter_bins),
nonzero(tally.num_nuclides),
nonzero(tally.num_score_bins))
sum = np.reshape(sum, new_shape)
sum_sq = np.reshape(sum_sq, new_shape)
# Set the data for this Tally
tally.sum = sum
tally.sum_sq = sum_sq
# Indicate that Tally results have been read
self._results_read = True
def compute_ci(self, confidence=0.95):
"""Computes confidence intervals for each Tally bin.

10
openmc/summary.py
View file

@ -532,18 +532,11 @@ class Summary(object):
# Iterate over all Tallies
for tally_key in tally_keys:
tally_id = int(tally_key.strip('tally '))
subbase = '{0}{1}'.format(base, tally_id)
# Read Tally name metadata
name_size = self._f['{0}/name_size'.format(subbase)][...]
if (name_size > 0):
tally_name = self._f['{0}/name'.format(subbase)][...][0]
tally_name = tally_name.lstrip('[\'')
tally_name = tally_name.rstrip('\']')
else:
tally_name = ''
tally_name = self._f['{0}/name'.format(subbase)].value.decode()
# Create Tally object and assign basic properties
tally = openmc.Tally(tally_id, tally_name)
@ -560,7 +553,6 @@ class Summary(object):
# Initialize all Filters
for j in range(1, num_filters+1):
subsubbase = '{0}/filter {1}'.format(subbase, j)
# Read filter type (e.g., "cell", "energy", etc.)

79
openmc/tallies.py
View file

@ -103,6 +103,9 @@ class Tally(object):
self._with_batch_statistics = False
self._derived = False
self._statepoint = None
self._results_read = False
def __deepcopy__(self, memo):
existing = memo.get(id(self))
@ -121,6 +124,8 @@ class Tally(object):
clone._with_summary = self.with_summary
clone._with_batch_statistics = self.with_batch_statistics
clone._derived = self.derived
clone._statepoint = self._statepoint
clone._results_read = self._results_read
clone._filters = []
for filter in self.filters:
@ -259,24 +264,66 @@ class Tally(object):
@property
def sum(self):
if not self._statepoint:
return None
if not self._results_read:
# Extract Tally data from the file
data = self._statepoint._f['tallies/tally {0}/results'.format(
self.id)].value
sum = data['sum']
sum_sq = data['sum_sq']
# Define a routine to convert 0 to 1
def nonzero(val):
return 1 if not val else val
# Reshape the results arrays
new_shape = (nonzero(self.num_filter_bins),
nonzero(self.num_nuclides),
nonzero(self.num_score_bins))
sum = np.reshape(sum, new_shape)
sum_sq = np.reshape(sum_sq, new_shape)
# Set the data for this Tally
self._sum = sum
self._sum_sq = sum_sq
# Indicate that Tally results have been read
self._results_read = True
return self._sum
@property
def sum_sq(self):
if not self._statepoint:
return None
if not self._results_read:
# Force reading of sum and sum_sq
self.sum
return self._sum_sq
@property
def mean(self):
# Compute the mean if needed
if self._mean is None:
self.compute_mean()
if not self._statepoint:
return None
self._mean = self.sum / self.num_realizations
return self._mean
@property
def std_dev(self):
# Compute the standard deviation if needed
if self._std_dev is None:
self.compute_std_dev()
if not self._statepoint:
return None
n = self.num_realizations
self._std_dev = np.sqrt((self.sum_sq/n - self.mean**2)/(n - 1))
self.with_batch_statistics = True
return self._std_dev
@property
@ -456,30 +503,6 @@ class Tally(object):
self._nuclides.remove(nuclide)
def compute_mean(self):
"""Compute the sample mean for each bin in the tally"""
# Calculate sample mean
self._mean = self.sum / self.num_realizations
def compute_std_dev(self, t_value=1.0):
"""Compute the sample standard deviation for each bin in the tally
Parameters
----------
t_value : float, optional
Student's t-value applied to the uncertainty. Defaults to 1.0,
meaning the reported value is the sample standard deviation.
"""
# Calculate sample standard deviation
self.compute_mean()
self._std_dev = np.sqrt((self.sum_sq / self.num_realizations -
self.mean**2) / (self.num_realizations - 1))
self._std_dev *= t_value
self.with_batch_statistics = True
def __repr__(self):
string = 'Tally\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self.id)

1
tests/test_entropy/test_entropy.py
View file

@ -14,7 +14,6 @@ class EntropyTestHarness(TestHarness):
# Read the statepoint file.
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
sp = StatePoint(statepoint)
sp.read_results()
# Write out k-combined.
outstr = 'k-combined:\n'

11
tests/test_fixed_source/test_fixed_source.py
View file

@ -14,17 +14,16 @@ class FixedSourceTestHarness(TestHarness):
# Read the statepoint file.
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
sp = StatePoint(statepoint)
sp.read_results()
# Write out tally data.
outstr = ''
if self._tallies:
tally_num = 1
for tally_ind in sp._tallies:
tally = sp._tallies[tally_ind]
results = np.zeros((tally._sum.size*2, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
for tally_ind in sp.tallies:
tally = sp.tallies[tally_ind]
results = np.zeros((tally.sum.size*2, ))
results[0::2] = tally.sum.ravel()
results[1::2] = tally.sum_sq.ravel()
results = ['{0:12.6E}'.format(x) for x in results]
outstr += 'tally ' + str(tally_num) + ':\n'

1
tests/test_sourcepoint_batch/test_sourcepoint_batch.py
View file

@ -21,7 +21,6 @@ class SourcepointTestHarness(TestHarness):
# Read the statepoint file.
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
sp = StatePoint(statepoint)
sp.read_results()
# Get the eigenvalue information.
outstr = TestHarness._get_results(self)

1
tests/test_sourcepoint_interval/test_sourcepoint_interval.py
View file

@ -21,7 +21,6 @@ class SourcepointTestHarness(TestHarness):
# Read the statepoint file.
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
sp = StatePoint(statepoint)
sp.read_results()
# Get the eigenvalue information.
outstr = TestHarness._get_results(self)

12
tests/testing_harness.py
View file

@ -93,7 +93,6 @@ class TestHarness(object):
# Read the statepoint file.
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
sp = StatePoint(statepoint)
sp.read_results()
# Write out k-combined.
outstr = 'k-combined:\n'
@ -103,11 +102,11 @@ class TestHarness(object):
# Write out tally data.
if self._tallies:
tally_num = 1
for tally_ind in sp._tallies:
tally = sp._tallies[tally_ind]
results = np.zeros((tally._sum.size*2, ))
results[0::2] = tally._sum.ravel()
results[1::2] = tally._sum_sq.ravel()
for tally_ind in sp.tallies:
tally = sp.tallies[tally_ind]
results = np.zeros((tally.sum.size*2, ))
results[0::2] = tally.sum.ravel()
results[1::2] = tally.sum_sq.ravel()
results = ['{0:12.6E}'.format(x) for x in results]
outstr += 'tally ' + str(tally_num) + ':\n'
@ -204,7 +203,6 @@ class CMFDTestHarness(TestHarness):
# Read the statepoint file.
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
sp = StatePoint(statepoint)
sp.read_results()
# Write out the eigenvalue and tallies.
outstr = TestHarness._get_results(self)