diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 126d6b97e8..5c71b3f5b3 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -43,8 +43,9 @@ class StatePoint(object): Shannon entropy of fission source at each batch gen_per_batch : int Number of fission generations per batch - global_tallies : ndarray - Global tallies and their uncertainties + global_tallies : ndarray of compound datatype + Global tallies for k-effective estimates and leakage. The compound + datatype has fields 'name', 'sum', 'sum_sq', 'mean', and 'std_dev'. k_combined : list Combined estimator for k-effective and its uncertainty k_col_abs : float @@ -102,6 +103,7 @@ class StatePoint(object): self._meshes_read = False self._tallies_read = False self._with_summary = False + self._global_tallies = None def close(self): self._f.close() @@ -177,8 +179,24 @@ class StatePoint(object): @property def global_tallies(self): - data = self._f['global_tallies'].value - return np.column_stack((data['sum'], data['sum_sq'])) + if self._global_tallies is None: + data = self._f['global_tallies'].value + gt = np.zeros_like(data, dtype=[ + ('name', 'a14'), ('sum', 'f8'), ('sum_sq', 'f8'), + ('mean', 'f8'), ('std_dev', 'f8')]) + gt['name'] = ['k-collision', 'k-absorption', 'k-tracklength', + 'leakage'] + gt['sum'] = data['sum'] + gt['sum_sq'] = data['sum_sq'] + + # Calculate mean and sample standard deviation of mean + n = self.n_realizations + gt['mean'] = gt['sum']/n + gt['std_dev'] = np.sqrt((gt['sum_sq']/n - gt['mean']**2)/(n - 1)) + + self._global_tallies = gt + + return self._global_tallies @property def k_cmfd(self): @@ -460,60 +478,6 @@ class StatePoint(object): def with_summary(self): return self._with_summary - def compute_ci(self, confidence=0.95): - """Computes confidence intervals for each Tally bin. - - This method is equivalent to calling compute_stdev(...) when the - confidence is known as opposed to its corresponding t value. - - Parameters - ---------- - confidence : float, optional - Confidence level. Defaults to 0.95. - - """ - - # Determine significance level and percentile for two-sided CI - alpha = 1 - confidence - percentile = 1 - alpha/2 - - # Calculate t-value - t_value = scipy.stats.t.ppf(percentile, self._n_realizations - 1) - self.compute_stdev(t_value) - - def compute_stdev(self, t_value=1.0): - """Computes the sample mean and the standard deviation of the mean - for each Tally bin. - - 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. - - """ - - # Determine number of realizations - n = self._n_realizations - - # Calculate the standard deviation for each global tally - for i in range(len(self._global_tallies)): - - # Get sum and sum of squares - s, s2 = self._global_tallies[i] - - # Calculate sample mean and replace value - s /= n - self._global_tallies[i, 0] = s - - # Calculate standard deviation - if s != 0.0: - self._global_tallies[i, 1] = t_value * np.sqrt((s2 / n - s**2) / (n-1)) - - # Calculate sample mean and standard deviation for user-defined Tallies - for tally_id, tally in self.tallies.items(): - tally.compute_std_dev(t_value) - def get_tally(self, scores=[], filters=[], nuclides=[], name=None, id=None, estimator=None): """Finds and returns a Tally object with certain properties. diff --git a/tests/test_fixed_source/test_fixed_source.py b/tests/test_fixed_source/test_fixed_source.py index 0595ea1dba..c3bd34856d 100644 --- a/tests/test_fixed_source/test_fixed_source.py +++ b/tests/test_fixed_source/test_fixed_source.py @@ -30,9 +30,10 @@ class FixedSourceTestHarness(TestHarness): outstr += '\n'.join(results) + '\n' tally_num += 1 + gt = sp.global_tallies outstr += 'leakage:\n' - outstr += '{0:12.6E}'.format(sp.global_tallies[3][0]) + '\n' - outstr += '{0:12.6E}'.format(sp.global_tallies[3][1]) + '\n' + outstr += '{0:12.6E}'.format(gt[gt['name'] == b'leakage'][0]['sum']) + '\n' + outstr += '{0:12.6E}'.format(gt[gt['name'] == b'leakage'][0]['sum_sq']) + '\n' return outstr