crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #447 from paulromano/new-hdf5-interface

Browse source 
Complete revamp of binary output featuring new HDF5 interface
This commit is contained in:
Sterling Harper 2015-09-14 16:46:30 -04:00
commit f8bfa401cb
46 changed files with 4404 additions and 7379 deletions
Download patch file Download diff file Expand all files Collapse all files

375
scripts/openmc-tally-convergence
View file

@ -1,375 +0,0 @@
#!/usr/bin/env python
# This program takes OpenMC statepoint binary files and creates a variety of
# outputs from them which should provide the user with an idea of the
# convergence behavior of all the tallies and filters defined by the user in
# tallies.xml. The program can directly plot the value and errors of each
# tally, filter, score combination; it can save these plots to a file; and
# it can also save the data used in these plots to a CSV file for importing in
# to other plotting packages such as Excel, gnuplot, MathGL, or Veusz.
# To use the program, run this program from the working directory of the openMC
# problem to analyze.
# The USER OPTIONS block below provides four options for the user to set:
# fileType, printxs, showImg, and savetoCSV. See the options block for more
# information.
from __future__ import print_function
from math import sqrt, pow
from glob import glob
import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
from openmc.statepoint import StatePoint
##################################### USER OPTIONS
# Set filetype (the file extension desired, without the period.)
# Options are backend dependent, but most backends support png, pdf, ps, eps
# and svg. Write "none" if no saved files are desired.
fileType = "none"
# Set if cross-sections or reaction rates are desired printxs = True means X/S
printxs = False
# Set if the figures should be displayed to screen or not (True means show)
showImg = False
# Save to CSV for use in more advanced plotting programs like GNUPlot, MathGL
savetoCSV = True
##################################### END USER OPTIONS
## Find if tallies.xml exists.
#if glob('./tallies.xml') != None:
# # It exists
# tallyData = talliesXML('tallies.xml')
#else:
# # It does not exist.
# tallyData = None
# Find all statepoints in this directory.
files = glob('./statepoint.*.binary')
fileNums = []
begin = 13
# Arrange the file list in increasing batch order
for i in range(len(files)):
end = files[i].find(".binary")
fileNums.append(int(files[i][begin:end]))
fileNums.sort()
# Re-make filenames
files = []
for i in range(len(fileNums)):
files.append("./statepoint." + str(fileNums[i]) + ".binary")
# Initialize arrays as needed
mean = [None for x in range(len(files))]
uncert = [None for x in range(len(files))]
scoreType = [None for x in range(len(files))]
active_batches = [None for x in range(len(files))]
for i_batch in range(len(files)):
# Get filename
batch_filename = files[i_batch]
# Create StatePoint object
sp = StatePoint(batch_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!")
# Increase the dimensionality of our main variables
mean[i_batch] = [None for x in range(len(sp.tallies))]
uncert[i_batch] = [None for x in range(len(sp.tallies))]
scoreType[i_batch] = [None for x in range(len(sp.tallies))]
# Loop over all tallies
for i_tally, 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)
# Store the batch count
active_batches[i_batch] = n
# Resize the 2nd dimension
mean[i_batch][i_tally] = [None for x in range(t.total_filter_bins)]
uncert[i_batch][i_tally] = [None for x in range(t.total_filter_bins)]
scoreType[i_batch][i_tally] = [None for x in range(t.total_filter_bins)]
for i_filter in range(t.total_filter_bins):
# Resize the 3rd dimension
mean[i_batch][i_tally][i_filter] = [None for x in range(t.n_nuclides)]
uncert[i_batch][i_tally][i_filter] = [None for x in range(t.n_nuclides)]
scoreType[i_batch][i_tally][i_filter] = [None for x in range(t.n_nuclides)]
print(t.total_filter_bins,t.n_nuclides)
for i_nuclide in range(t.n_nuclides):
mean[i_batch][i_tally][i_filter][i_nuclide] = \
[None for x in range(t.n_scores)]
uncert[i_batch][i_tally][i_filter][i_nuclide] = \
[None for x in range(t.n_scores)]
scoreType[i_batch][i_tally][i_filter][i_nuclide] = \
[None for x in range(t.n_scores)]
for i_score in range(t.n_scores):
scoreType[i_batch][i_tally][i_filter][i_nuclide][i_score] = \
t.scores[i_score]
s, s2 = sp._get_double(2)
s /= n
mean[i_batch][i_tally][i_filter][i_nuclide][i_score] = s
if s != 0.0:
relative_error = t_value*sqrt((s2/n - s*s)/(n-1))/s
else:
relative_error = 0.0
uncert[i_batch][i_tally][i_filter][i_nuclide][i_score] = relative_error
# Reorder the data lists in to a list order more conducive for plotting:
# The indexing should be: [tally][filter][score][batch]
meanPlot = [None for x in range(len(mean[0]))] # Set to the number of tallies
uncertPlot = [None for x in range(len(mean[0]))] # Set to the number of tallies
absUncertPlot = [None for x in range(len(mean[0]))] # Set to number of tallies
filterLabel = [None for x in range(len(mean[0]))] #Set to the number of tallies
fluxLoc = [None for x in range(len(mean[0]))] # Set to the number of tallies
printxs = [False for x in range(len(mean[0]))] # Set to the number of tallies
# Get and set the correct sizes for the rest of the dimensions
for i_tally in range(len(meanPlot)):
# Set 2nd (score) dimension
meanPlot[i_tally] = [None for x in range(len(mean[0][i_tally]))]
uncertPlot[i_tally] = [None for x in range(len(mean[0][i_tally]))]
absUncertPlot[i_tally] = [None for x in range(len(mean[0][i_tally]))]
filterLabel[i_tally] = [None for x in range(len(mean[0][i_tally]))]
# Initialize flux location so it will be -1 if not found
fluxLoc[i_tally] = -1
for i_filter in range(len(meanPlot[i_tally])):
# Set 3rd (filter) dimension
meanPlot[i_tally][i_filter] = \
[None for x in range(len(mean[0][i_tally][i_filter]))]
uncertPlot[i_tally][i_filter] = \
[None for x in range(len(mean[0][i_tally][i_filter]))]
absUncertPlot[i_tally][i_filter] = \
[None for x in range(len(mean[0][i_tally][i_filter]))]
filterLabel[i_tally][i_filter] = \
[None for x in range(len(mean[0][i_tally][i_filter]))]
for i_nuclide in range(len(meanPlot[i_tally][i_filter])):
# Set 4th (nuclide)) dimension
meanPlot[i_tally][i_filter][i_nuclide] = \
[None for x in range(len(mean[0][i_tally][i_filter][i_nuclide]))]
uncertPlot[i_tally][i_filter][i_nuclide] = \
[None for x in range(len(mean[0][i_tally][i_filter][i_nuclide]))]
absUncertPlot[i_tally][i_filter][i_nuclide] = \
[None for x in range(len(mean[0][i_tally][i_filter][i_nuclide]))]
for i_score in range(len(meanPlot[i_tally][i_filter][i_nuclide])):
# Set 5th (batch) dimension
meanPlot[i_tally][i_filter][i_nuclide][i_score] = \
[None for x in range(len(mean))]
uncertPlot[i_tally][i_filter][i_nuclide][i_score] = \
[None for x in range(len(mean))]
absUncertPlot[i_tally][i_filter][i_nuclide][i_score] = \
[None for x in range(len(mean))]
# Get filterLabel (this should be moved to its own function)
#??? How to do?
# Set flux location if found
# all batches and all tallies will have the same score ordering, hence
# the 0's in the 1st, 3rd, and 4th dimensions.
if scoreType[0][i_tally][0][0][i_score] == 'flux':
fluxLoc[i_tally] = i_score
# Set printxs array according to the printxs input
if printxs:
for i_tally in range(len(fluxLoc)):
if fluxLoc[i_tally] != -1:
printxs[i_tally] = True
# Now rearrange the data as suitable, and perform xs conversion if necessary
for i_batch in range(len(mean)):
for i_tally in range(len(mean[i_batch])):
for i_filter in range(len(mean[i_batch][i_tally])):
for i_nuclide in range(len(mean[i_batch][i_tally][i_filter])):
for i_score in range(len(mean[i_batch][i_tally][i_filter][i_nuclide])):
if (printxs[i_tally] and \
((scoreType[0][i_tally][i_filter][i_nuclide][i_score] != 'flux') and \
(scoreType[0][i_tally][i_filter][i_nuclide][i_score] != 'current'))):
# Perform rate to xs conversion
# mean is mean/fluxmean
meanPlot[i_tally][i_filter][i_nuclide][i_score][i_batch] = \
mean[i_batch][i_tally][i_filter][i_nuclide][i_score] / \
mean[i_batch][i_tally][i_filter][i_nuclide][fluxLoc[i_tally]]
# Update the relative uncertainty via error propagation
uncertPlot[i_tally][i_filter][i_nuclide][i_score][i_batch] = \
sqrt(pow(uncert[i_batch][i_tally][i_filter][i_nuclide][i_score],2) \
+ pow(uncert[i_batch][i_tally][i_filter][i_nuclide][fluxLoc[i_tally]],2))
else:
# Do not perform rate to xs conversion
meanPlot[i_tally][i_filter][i_nuclide][i_score][i_batch] = \
mean[i_batch][i_tally][i_filter][i_nuclide][i_score]
uncertPlot[i_tally][i_filter][i_nuclide][i_score][i_batch] = \
uncert[i_batch][i_tally][i_filter][i_nuclide][i_score]
# Both have the same absolute uncertainty calculation
absUncertPlot[i_tally][i_filter][i_nuclide][i_score][i_batch] = \
uncert[i_batch][i_tally][i_filter][i_nuclide][i_score] * \
mean[i_batch][i_tally][i_filter][i_nuclide][i_score]
# Set plotting constants
xLabel = "Batches"
xLabel = xLabel.title() # not necessary for now, but is left in to handle if
# the previous line changes
# Begin plotting
for i_tally in range(len(meanPlot)):
# Set tally string (placeholder until I put tally labels in statePoint)
tallyStr = "Tally " + str(i_tally + 1)
for i_filter in range(len(meanPlot[i_tally])):
# Set filter string
filterStr = "Filter " + str(i_filter + 1)
for i_nuclide in range(len(meanPlot[i_tally][i_filter])):
nuclideStr = "Nuclide " + str(i_nuclide + 1)
for i_score in range(len(meanPlot[i_tally][i_filter][i_nuclide])):
# Set score string
scoreStr = scoreType[i_batch][i_tally][i_filter][i_nuclide][i_score]
scoreStr = scoreStr.title()
if (printxs[i_tally] and ((scoreStr != 'Flux') and \
(scoreStr != 'Current'))):
scoreStr = scoreStr + "-XS"
# set Title
title = "Convergence of " + scoreStr + " in " + tallyStr + " for "\
+ filterStr + " and " + nuclideStr
# set yLabel
yLabel = scoreStr
yLabel = yLabel.title()
# Set saving filename
fileName = "tally_" + str(i_tally + 1) + "_" + scoreStr + \
"_filter_" + str(i_filter+1) + "_nuclide_" + str(i_nuclide+1) \
+ "." + fileType
REfileName = "tally_" + str(i_tally + 1) + "_" + scoreStr + \
"RE_filter_" + str(i_filter+1) + "_nuclide_" + str(i_nuclide+1) \
+ "." + fileType
# Plot mean with absolute error bars
plt.errorbar(active_batches, \
meanPlot[i_tally][i_filter][i_nuclide][i_score][:], \
absUncertPlot[i_tally][i_filter][i_nuclide][i_score][:],fmt='o-',aa=True)
plt.xlabel(xLabel)
plt.ylabel(yLabel)
plt.title(title)
if (fileType != 'none'):
plt.savefig(fileName)
if showImg:
plt.show()
plt.clf()
# Plot relative uncertainty
plt.plot(active_batches, \
uncertPlot[i_tally][i_filter][i_nuclide][i_score][:],'o-',aa=True)
plt.xlabel(xLabel)
plt.ylabel("Relative Error of " + yLabel)
plt.title("Relative Error of " + title)
if (fileType != 'none'):
plt.savefig(REfileName)
if showImg:
plt.show()
plt.clf()
if savetoCSV:
# This block loops through each tally, and for each tally:
# Creates a new file
# Writes the scores and filters for that tally in csv format.
# The columns will be: batches,then for each filter: all the scores
# The rows, of course, are the data points per batch.
for i_tally in range(len(meanPlot)):
# Set tally string (placeholder until I put tally labels in statePoint)
tallyStr = "Tally " + str(i_tally + 1)
CSV_filename = "./tally" + str(i_tally+1)+".csv"
# Open the file
f = open(CSV_filename, 'w')
# Write the header line
lineText = "Batches"
for i_filter in range(len(meanPlot[i_tally])):
# Set filter string
filterStr = "Filter " + str(i_filter + 1)
for i_nuclide in range(len(meanPlot[i_tally][i_filter])):
nuclideStr = "Nuclide " + str(i_nuclide + 1)
for i_score in range(len(meanPlot[i_tally][i_filter][i_nuclide])):
# Set the title
scoreStr = scoreType[i_batch][i_tally][i_filter][i_nuclide][i_score]
scoreStr = scoreStr.title()
if (printxs[i_tally] and ((scoreStr != 'Flux') and \
(scoreStr != 'Current'))):
scoreStr = scoreStr + "-XS"
# set header
headerText = scoreStr + " for " + filterStr + " for " + nuclideStr
lineText = lineText + "," + headerText + \
",Abs Unc of " + headerText + \
",Rel Unc of " + headerText
f.write(lineText + "\n")
# Write the data lines, each row is a different batch
for i_batch in range(len(meanPlot[i_tally][0][0][0])):
lineText = repr(active_batches[i_batch])
for i_filter in range(len(meanPlot[i_tally])):
for i_nuclide in range(len(meanPlot[i_tally][i_filter])):
for i_score in range(len(meanPlot[i_tally][i_filter][i_nuclide])):
fieldText = \
repr(meanPlot[i_tally][i_filter][i_nuclide][i_score][i_batch]) + \
"," + \
repr(absUncertPlot[i_tally][i_filter][i_nuclide][i_score][i_batch]) +\
"," + \
repr(uncertPlot[i_tally][i_filter][i_nuclide][i_score][i_batch])
lineText = lineText + "," + fieldText
f.write(lineText + "\n")

47
scripts/openmc-track-to-vtk
View file

@ -18,6 +18,7 @@ Usage information can be obtained by running 'track.py --help':
import os
import argparse
import h5py
import struct
import vtk
@ -39,56 +40,26 @@ def main():
# Parse commandline arguments.
args = _parse_args()
# Check input file extensions.
for fname in args.input:
if not (fname.endswith('.h5') or fname.endswith('.binary')):
raise ValueError("Input file names must either end with '.h5' or"
"'.binary'.")
# Make sure that the output filename ends with '.pvtp'.
if not args.out:
args.out = 'tracks.pvtp'
elif not args.out.endswith('.pvtp'):
args.out += '.pvtp'
# Import HDF library if HDF files are present
for fname in args.input:
if fname.endswith('.h5'):
import h5py
break
# Initialize data arrays and offset.
points = vtk.vtkPoints()
cells = vtk.vtkCellArray()
point_offset = 0
for fname in args.input:
# Write coordinate values to points array.
if fname.endswith('.binary'):
track = open(fname, 'rb')
# Determine number of particles and tracks/particle
n_particles = struct.unpack('i', track.read(4))[0]
n_coords = struct.unpack('i'*n_particles, track.read(4*n_particles))
coords = []
for i in range(n_particles):
# Read coordinates for each particle
coords.append([struct.unpack('ddd', track.read(24))
for j in range(n_coords[i])])
# Add coordinates to points data
for triplet in coords[i]:
points.InsertNextPoint(triplet)
else:
track = h5py.File(fname)
n_particles = track['n_particles'].value[0]
n_coords = track['n_coords']
coords = []
for i in range(n_particles):
coords.append(track['coordinates_' + str(i + 1)].value)
for j in range(n_coords[i]):
points.InsertNextPoint(coords[i][j,:])
track = h5py.File(fname)
n_particles = track['n_particles'].value
n_coords = track['n_coords']
coords = []
for i in range(n_particles):
coords.append(track['coordinates_' + str(i + 1)].value)
for j in range(n_coords[i]):
points.InsertNextPoint(coords[i][j,:])
for i in range(n_particles):
# Create VTK line and assign points to line.