diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d3785aa6a..5fb11848d 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -261,15 +261,12 @@ set_target_properties(${program} PROPERTIES #=============================================================================== install(TARGETS ${program} RUNTIME DESTINATION bin) -install(PROGRAMS utils/statepoint_cmp.py +install(PROGRAMS utils/plot_mesh_tally.py DESTINATION bin - RENAME statepoint_cmp) + RENAME plot_mesh_tally) install(PROGRAMS utils/statepoint_histogram.py DESTINATION bin RENAME statepoint_histogram) -install(PROGRAMS utils/statepoint_meshplot.py - DESTINATION bin - RENAME statepoint_meshplot) install(PROGRAMS utils/update_inputs.py DESTINATION bin RENAME update_inputs) diff --git a/src/state_point.F90 b/src/state_point.F90 index fb9402e21..5cbd93677 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -307,7 +307,7 @@ contains MOMENT_LOOP: do j = 1, tally % n_user_score_bins select case(tally % score_bins(k)) case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) - moment_name = 'P' // to_str(tally % moment_order(k)) + moment_name = 'P' // trim(to_str(tally % moment_order(k))) call sp % write_data(moment_name, "order" // trim(to_str(k)), & group="tallies/tally " // trim(to_str(tally % id)) // & "/moments") diff --git a/src/utils/particle_restart.py b/src/utils/openmc/particle_restart.py similarity index 100% rename from src/utils/particle_restart.py rename to src/utils/openmc/particle_restart.py diff --git a/src/utils/openmc/statepoint.py b/src/utils/openmc/statepoint.py index 765346c12..8896fd6d5 100644 --- a/src/utils/openmc/statepoint.py +++ b/src/utils/openmc/statepoint.py @@ -1,5 +1,6 @@ import copy import struct +import sys import numpy as np import scipy.stats @@ -7,6 +8,9 @@ import scipy.stats import openmc from openmc.constants import * +if sys.version > '3': + long = int + class SourceSite(object): @@ -111,6 +115,17 @@ class StatePoint(object): def tallies(self): return self._tallies + @property + def tallies_present(self): + return self._tallies_present + + @property + def global_tallies(self): + return self._global_tallies + + @property + def n_realizations(self): + return self._n_realizations def _read_metadata(self): @@ -372,7 +387,8 @@ class StatePoint(object): filter.num_bins = n_bins if FILTER_TYPES[filter_type] == 'mesh': - filter.mesh = self._meshes[bins] + key = self._mesh_keys[self._mesh_ids.index(bins)] + filter.mesh = self._meshes[key] # Add Filter to the Tally tally.add_filter(filter) diff --git a/src/utils/openmc/tallies.py b/src/utils/openmc/tallies.py index 810611b26..b5817a60c 100644 --- a/src/utils/openmc/tallies.py +++ b/src/utils/openmc/tallies.py @@ -1124,10 +1124,10 @@ class Tally(object): # Convert (x,y,z) to a single bin -- this is similar to # subroutine mesh_indices_to_bin in openmc/src/mesh.F90. - value = ((filter_bins[i][0] - 1) * ny * nz + - (filter_bins[i][1] - 1) * nz + - (filter_bins[i][2] - 1)) - filter_index += value * test_filter._stride + val = ((filter_bins[i][0] - 1) * ny * nz + + (filter_bins[i][1] - 1) * nz + + (filter_bins[i][2] - 1)) + filter_index += val * test_filter._stride # Filter bins for distribcell are the "IDs" of each unique placement # of the Cell in the Geometry (integers starting at 0) diff --git a/src/utils/plot_mesh_tally.py b/src/utils/plot_mesh_tally.py index b78305ebb..5c46d927c 100755 --- a/src/utils/plot_mesh_tally.py +++ b/src/utils/plot_mesh_tally.py @@ -11,7 +11,7 @@ from matplotlib.figure import Figure import matplotlib.pyplot as plt import numpy as np -from statepoint import * +from openmc.statepoint import StatePoint if sys.version_info[0] < 3: import Tkinter as tk @@ -130,14 +130,13 @@ class MeshPlotter(tk.Frame): selectedTally = self.datafile.tallies[tally_id] # Get mesh for selected tally - self.mesh = self.datafile.meshes[ - selectedTally.filters['mesh'].bins[0] - 1] + self.mesh = selectedTally.filters_by_name['mesh'].mesh # Get mesh dimensions self.nx, self.ny, self.nz = self.mesh.dimension # Repopulate comboboxes baesd on current basis selection - text = self.basisBox['values'][self.basisBox.current()] + text = self.basisBox.get() if text == 'xy': self.axialBox['values'] = [str(i+1) for i in range(self.nz)] elif text == 'yz': @@ -164,7 +163,8 @@ class MeshPlotter(tk.Frame): # create a label/combobox for each filter in selected tally count = 0 - for filterType in selectedTally.filters: + for f in selectedTally.filters: + filterType = f.type if filterType == 'mesh': continue count += 1 @@ -176,7 +176,6 @@ class MeshPlotter(tk.Frame): self.filterBoxes[filterType] = combobox # Set combobox items - f = selectedTally.filters[filterType] if filterType in ['energyin', 'energyout']: combobox['values'] = ['{0} to {1}'.format(*f.bins[i:i+2]) for i in range(f.length)] @@ -198,7 +197,7 @@ class MeshPlotter(tk.Frame): def redraw(self, event=None): basis = self.basisBox.current() + 1 axial_level = self.axialBox.current() + 1 - is_mean = self.meanBox.current() + mbvalue = self.meanBox.get() # Get selected tally tally_id = self.meshTallies[self.tallyBox.current()] @@ -206,71 +205,45 @@ class MeshPlotter(tk.Frame): # Create spec_list spec_list = [] - for f in selectedTally.filters.values(): + for f in selectedTally.filters: if f.type == 'mesh': + mesh_filter = f continue index = self.filterBoxes[f.type].current() - spec_list.append((f.type, index)) + spec_list.append((f, index)) - # Take is_mean and convert it to an index of the score - score_loc = is_mean - if score_loc > 1: - score_loc = 1 - - text = self.basisBox['values'][self.basisBox.current()] + text = self.basisBox.get() if text == 'xy': - matrix = np.zeros((self.nx, self.ny)) - for i in range(self.nx): - for j in range(self.ny): - matrix[i, j] = self.datafile.get_value(tally_id, - spec_list + [('mesh', (i + 1, j + 1, axial_level))], - self.scoreBox.current())[score_loc] - # Calculate relative uncertainty from absolute, if requested - if is_mean == 2: - # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(tally_id, - spec_list + [('mesh', (i + 1, j + 1, axial_level))], - self.scoreBox.current())[0] - if mean_val > 0.0: - matrix[i, j] = matrix[i, j] / mean_val - else: - matrix[i, j] = 0.0 - + dims = (self.nx, self.ny) elif text == 'yz': - matrix = np.zeros((self.ny, self.nz)) - for i in range(self.ny): - for j in range(self.nz): - matrix[i, j] = self.datafile.get_value(tally_id, - spec_list + [('mesh', (axial_level, i + 1, j + 1))], - self.scoreBox.current())[score_loc] - # Calculate relative uncertainty from absolute, if requested - if is_mean == 2: - # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(tally_id, - spec_list + [('mesh', (axial_level, i + 1, j + 1))], - self.scoreBox.current())[0] - if mean_val > 0.0: - matrix[i, j] = matrix[i, j] / mean_val - else: - matrix[i, j] = 0.0 - + dims = (self.ny, self.nz) else: - matrix = np.zeros((self.nx, self.nz)) - for i in range(self.nx): - for j in range(self.nz): - matrix[i, j] = self.datafile.get_value(tally_id, - spec_list + [('mesh', (i + 1, axial_level, j + 1))], - self.scoreBox.current())[score_loc] - # Calculate relative uncertainty from absolute, if requested - if is_mean == 2: - # Take care to handle zero means when normalizing - mean_val = self.datafile.get_value(tally_id, - spec_list + [('mesh', (i + 1, axial_level, j + 1))], - self.scoreBox.current())[0] - if mean_val > 0.0: - matrix[i, j] = matrix[i, j] / mean_val - else: - matrix[i, j] = 0.0 + dims = (self.nx, self.nz) + matrix = np.zeros(dims) + for i in range(dims[0]): + for j in range(dims[1]): + if text == 'xy': + meshtuple = (i + 1, j + 1, axial_level) + elif text == 'yz': + meshtuple = (axial_level, i + 1, j + 1) + else: + meshtuple = (i + 1, axial_level, j + 1) + filters, filter_bins = zip(*spec_list + [ + (mesh_filter, meshtuple)]) + mean = selectedTally.get_value( + self.scoreBox.get(), filters, filter_bins) + stdev = selectedTally.get_value( + self.scoreBox.get(), filters, filter_bins, + value='std_dev') + if mbvalue == 'Mean': + matrix[i, j] = mean + elif mbvalue == 'Absolute uncertainty': + matrix[i, j] = stdev + else: + if mean > 0.: + matrix[i, j] = stdev/mean + else: + matrix[i, j] = 0. # Clear the figure self.fig.clear() @@ -292,13 +265,14 @@ class MeshPlotter(tk.Frame): # Create StatePoint object and read in data self.datafile = StatePoint(filename) self.datafile.read_results() - self.datafile.generate_stdev() + self.datafile.compute_stdev() # Find which tallies are mesh tallies self.meshTallies = [] - for itally, tally in enumerate(self.datafile.tallies): - if 'mesh' in tally.filters: + for itally, tally in self.datafile.tallies.iteritems(): + if any([f.type == 'mesh' for f in tally.filters]): self.meshTallies.append(itally) + tally.filters_by_name = {f.type: f for f in tally.filters} if not self.meshTallies: messagebox.showerror("Invalid StatePoint File", diff --git a/src/utils/statepoint_cmp.py b/src/utils/statepoint_cmp.py deleted file mode 100755 index 7ad284370..000000000 --- a/src/utils/statepoint_cmp.py +++ /dev/null @@ -1,86 +0,0 @@ -#!/usr/bin/env python - -import sys - -from numpy.testing import assert_allclose, assert_equal - -from openmc.statepoint import StatePoint - -if len(sys.argv) > 2: - path1 = sys.argv[1] - path2 = sys.argv[2] -else: - raise - -# Create StatePoint objects -sp1 = StatePoint(path1) -sp2 = StatePoint(path2) - -# Read tally results -sp1.read_results() -sp2.read_results() - -# Compare header information -assert sp1.revision == sp2.revision -assert sp1.version == sp2.version -assert sp1.seed == sp2.seed -assert sp1.run_mode == sp2.run_mode -assert sp1.n_particles == sp2.n_particles -assert sp1.n_batches == sp2.n_batches - -if sp1.run_mode == 2: - # Compare criticality information - assert sp1.n_inactive == sp2.n_inactive - assert sp1.gen_per_batch == sp2.gen_per_batch - assert sp1.current_batch == sp2.current_batch - - # Compare keff results - assert_allclose(sp1.k_batch, sp2.k_batch) - - # Compare entropy results - assert_allclose(sp1.entropy, sp2.entropy) - -# Compare global tallies -assert_allclose(sp1.global_tallies, sp2.global_tallies) - -# Compare meshes -assert len(sp1.meshes) == len(sp2.meshes) -for m1, m2 in zip(sp1.meshes, sp2.meshes): - assert m1.type == m2.type - assert_equal(m1.dimension, m2.dimension) - assert_equal(m1.lower_left, m2.lower_left) - assert_equal(m1.upper_right, m2.upper_right) - assert_equal(m1.width, m2.width) - -# Compare tallies -assert len(sp1.tallies) == len(sp2.tallies) -for t1, t2 in zip(sp1.tallies, sp2.tallies): - # Compare size of tallies - assert t1.total_score_bins == t2.total_score_bins - assert t1.total_filter_bins == t2.total_filter_bins - - # Compare filters - assert len(t1.filters) == len(t2.filters) - for f1, f2 in zip(t1.filters.values(), t2.filters.values()): - assert f1.type == f2.type - assert f1.length == f2.length - assert_equal(f1.bins, f2.bins) - - # Compare nuclide and score bins - assert_equal(t1.nuclides, t2.nuclides) - assert_equal(t1.scores, t2.scores) - - # Compare tally results - assert_allclose(t1.results, t2.results) - -# If criticality, compare source sites -if sp1.run_mode == 2: - sp1.read_source() - sp2.read_source() - - assert len(sp1.source) == len(sp2.source) - for s1, s2 in zip(sp1.source, sp2.source): - assert s1.weight == s2.weight - assert_allclose(s1.xyz, s2.xyz) - assert_allclose(s1.uvw, s2.uvw) - assert s1.E == s2.E diff --git a/src/utils/statepoint_histogram.py b/src/utils/statepoint_histogram.py index ed7b2682b..26e8bdae6 100755 --- a/src/utils/statepoint_histogram.py +++ b/src/utils/statepoint_histogram.py @@ -13,60 +13,31 @@ from openmc.statepoint import StatePoint # Get filename filename = argv[1] -# Determine axial level -axial_level = int(argv[2]) if len(argv) > 2 else 1 -score = int(argv[3]) if len(argv) > 3 else 1 - # Create StatePoint object sp = StatePoint(filename) - -# Read number of realizations for global tallies -sp.n_realizations = sp._get_int()[0] - -# Read global tallies -n_global_tallies = sp._get_int()[0] -sp.global_tallies = np.array(sp._get_double(2*n_global_tallies)) -sp.global_tallies.shape = (n_global_tallies, 2) - -# Flag indicating if tallies are present -tallies_present = sp._get_int()[0] +sp.read_results() +sp.compute_ci() # Check if tallies are present -if not tallies_present: +if not sp.tallies_present: raise Exception("No tally data in state point!") # Loop over all tallies -for i, t in enumerate(sp.tallies): - # Calculate t-value for 95% two-sided CI - n = t.n_realizations - t_value = scipy.stats.t.ppf(0.975, n - 1) +for i, t in sp.tallies.items(): + # Determine relative error and fraction of bins with less than 1% half-width + # of CI + n_bins = t.mean.size + relative_error = t.std_dev[t.mean > 0.] / t.mean[t.mean > 0.] + fraction = float(sum(relative_error < 0.01))/n_bins - # Create lists for tallies - mean = [] - uncertainties = [] - nonzero = [] - - n_bins = t.total_score_bins * t.total_filter_bins - - # Loop over filter/score bins - for j in range(n_bins): - # Read sum and sum-squared - s, s2 = sp._get_double(2) - s /= n - if s != 0.0: - relative_error = t_value*sqrt((s2/n - s*s)/(n-1))/s - uncertainties.append(relative_error) - - # Display number of tallies with less than 1% CIs - fraction = (sum([1.0 if re < 0.01 else 0.0 for re in uncertainties]) - / n_bins) - print("Tally " + str(i+1)) + # Display results + print("Tally " + str(i)) print(" Fraction under 1% = {0}".format(fraction)) - print(" Min relative error = {0}".format(min(uncertainties))) - print(" Max relative error = {0}".format(max(uncertainties))) + print(" Min relative error = {0}".format(min(relative_error))) + print(" Max relative error = {0}".format(max(relative_error))) print(" Non-scoring bins = {0}".format( - 1.0 - float(len(uncertainties))/n_bins)) + 1.0 - float(relative_error.size)/n_bins)) # Plot histogram - plt.hist(uncertainties, 100) + plt.hist(relative_error, 100) plt.show() diff --git a/src/utils/statepoint_meshplot.py b/src/utils/statepoint_meshplot.py deleted file mode 100755 index ba09ef5e6..000000000 --- a/src/utils/statepoint_meshplot.py +++ /dev/null @@ -1,95 +0,0 @@ -#!/usr/bin/env python - -from __future__ import print_function, division - -from sys import argv -from math import sqrt - -import numpy as np -import scipy.stats -import matplotlib.pyplot as plt - -from openmc.statepoint import StatePoint - -# Get filename -filename = argv[1] - -# Determine axial level -axial_level = int(argv[2]) if len(argv) > 2 else 1 -score = int(argv[3]) if len(argv) > 3 else 1 - -# Create StatePoint object -sp = StatePoint(filename) - -# Read number of realizations for global tallies -sp.n_realizations = sp._get_int()[0] - -# Read global tallies -n_global_tallies = sp._get_int()[0] -sp.global_tallies = np.array(sp._get_double(2*n_global_tallies)) -sp.global_tallies.shape = (n_global_tallies, 2) - -# Flag indicating if tallies are present -tallies_present = sp._get_int()[0] - -# Check if tallies are present -if not tallies_present: - raise Exception("No tally data in state point!") - -# Loop over all tallies -print("Reading data...") -for t in sp.tallies: - # Calculate t-value for 95% two-sided CI - n = t.n_realizations - t_value = scipy.stats.t.ppf(0.975, n - 1) - - n_bins = t.total_score_bins * t.total_filter_bins - - # Check for mesh - if 'mesh' not in t.filters: - continue - - # Get Mesh object and determine size - m = sp.meshes[t.filters['mesh'].bins[0] - 1] - if len(m.dimension) == 2: - nx, ny = m.dimension - nz = 1 - else: - nx, ny, nz = m.dimension - - # Calculate number of score bins - ns = t.total_score_bins * t.total_filter_bins // (nx*ny*nz) - assert n_bins == nx*ny*nz*ns - - # Create lists for tallies - mean = np.zeros((nx, ny)) - error = np.zeros((nx, ny)) - criteria = np.zeros((nx, ny)) - - # Determine starting position for data - start = sp._f.tell() + (axial_level-1)*ns*16 + (score - 1)*16 - - for x in range(nx): - for y in range(ny): - # Seek to position of data - sp._f.seek(start + x*ny*nz*ns*16 + y*nz*ns*16) - - # Read sum and sum-squared - s, s2 = sp._get_double(2) - s /= n - mean[x, y] = s - if s != 0.0: - error[x, y] = t_value*sqrt((s2/n - s*s)/(n-1))/s - criteria[x, y] = 1.0 if error[x, y] < 0.05 else 0.0 - - # Make figure - print("Making colorplot...") - plt.imshow(mean, interpolation="nearest") - plt.colorbar() - plt.xlim((0, nx)) - plt.ylim((0, ny)) - plt.xticks(np.linspace(0, nx, 5), - np.linspace(m.lower_left[0], m.upper_right[0], 5)) - plt.yticks(np.linspace(0, ny, 5), - np.linspace(m.lower_left[1], m.upper_right[1], 5)) - plt.show() diff --git a/tests/test_particle_restart_eigval/results.py b/tests/test_particle_restart_eigval/results.py index ae446b855..2ff1ef2b8 100644 --- a/tests/test_particle_restart_eigval/results.py +++ b/tests/test_particle_restart_eigval/results.py @@ -5,7 +5,7 @@ import sys sys.path.insert(0, '../../src/utils') # import particle restart -import particle_restart as pr +import openmc.particle_restart as pr # read in particle restart file if len(sys.argv) > 1: diff --git a/tests/test_particle_restart_fixed/results.py b/tests/test_particle_restart_fixed/results.py index a5bac1da9..6c9a40d5a 100644 --- a/tests/test_particle_restart_fixed/results.py +++ b/tests/test_particle_restart_fixed/results.py @@ -4,8 +4,8 @@ import sys sys.path.insert(0, '../../src/utils') -# import particle restart -import particle_restart as pr +# import particle restart +import openmc.particle_restart as pr # read in particle restart file if len(sys.argv) > 1: @@ -15,8 +15,8 @@ else: # set up output string outstr = '' - -# write out properties + +# write out properties outstr += 'current batch:\n' outstr += "{0:12.6E}\n".format(p.current_batch) outstr += 'current gen:\n'