Merge remote-tracking branch 'upstream/develop' into squash-distribcells

This commit is contained in:
Will Boyd 2015-04-30 23:31:04 -04:00
commit bf89189273
11 changed files with 86 additions and 309 deletions

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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