mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Merge pull request #154 from bhermanmit/master
Added Method for Tally statepoint binary extraction
This commit is contained in:
Paul Romano
commit
b988630357
2 changed files with 262 additions and 198 deletions
|
|
@ -1,206 +1,160 @@
|
|||
#!/usr/bin/python
|
||||
# Filename: eigenfunction_rms.py
|
||||
|
||||
import h5py
|
||||
# import packages
|
||||
import statepoint
|
||||
import numpy as np
|
||||
import cPickle
|
||||
import matplotlib.pyplot as plt
|
||||
import os
|
||||
import sys
|
||||
#
|
||||
class EigenFunction:
|
||||
#
|
||||
'''Represents the reference eigenfunction'''
|
||||
reference = 0.
|
||||
#
|
||||
def __init__(self,data):
|
||||
#
|
||||
'''Initializes the eigenfunction'''
|
||||
self.function = 0.
|
||||
self.data = data
|
||||
self.meanfunction = 0.
|
||||
#
|
||||
def read_hdf5(self,h5_file,cycle):
|
||||
#
|
||||
'''Read data from HDF5 file'''
|
||||
self.cycle = cycle
|
||||
f = h5py.File(h5_file,'r')
|
||||
group = '/cycle'+str(cycle)+'/'+self.data
|
||||
dataset = f[group]
|
||||
self.function = np.empty(dataset.shape,dataset.dtype)
|
||||
dataset.read_direct(self.function)
|
||||
self.function = self.function
|
||||
self.iamref = 'F'
|
||||
#
|
||||
def set_reference(self):
|
||||
#
|
||||
'''Sets instance to be reference calc'''
|
||||
self.iamref = 'T'
|
||||
EigenFunction.reference = self.function
|
||||
#
|
||||
def compute_rms(self):
|
||||
#
|
||||
'''Computes RMS value'''
|
||||
Np = self.function.size
|
||||
Np = 41772
|
||||
tmp = (self.meanfunction - EigenFunction.reference)**2
|
||||
tmp2 = tmp.sum()
|
||||
self.rms = np.sqrt((1.0/float(Np))*tmp2)
|
||||
#
|
||||
def read_runs(runpath,hdfile,cycle_start,cycle_end,run_start,run_end,data):
|
||||
runlist = []
|
||||
tmp = EigenFunction(data)
|
||||
tmp.read_hdf5(runpath+str(run_start)+'/'+hdfile,cycle_start)
|
||||
indices = tmp.function.shape # extent of all dimensions
|
||||
i = cycle_start
|
||||
while i <= cycle_end:
|
||||
meantmp = EigenFunction(data) # init mean object
|
||||
j = run_start
|
||||
runs = np.zeros((run_end - run_start + 1,indices[0],indices[1],indices[2],indices[3])) # init runs array
|
||||
while j <= run_end:
|
||||
tmp.read_hdf5(runpath+str(j)+'/'+hdfile,i) # read hdf5 file
|
||||
runs[j-run_start] = tmp.function # put function into runs
|
||||
if i == cycle_start:
|
||||
print runs[j-run_start,0,150,150,0]
|
||||
j += 1
|
||||
meantmp.meanfunction = np.average(runs, axis=0) # compute the mean
|
||||
meantmp.function = runs
|
||||
runlist.append(meantmp)
|
||||
print 'Read in from path: '+runpath+' Cycle: '+str(i)
|
||||
i += 1
|
||||
return runlist
|
||||
|
||||
def create_reference(runpath,hdfile,cycle,run_start,run_end,data):
|
||||
def main(tally_id, score_id, batch_start, batch_end, name):
|
||||
|
||||
# calculate reference solution
|
||||
print 'Calculating Reference solution...'
|
||||
tmp = EigenFunction(data)
|
||||
tmp.read_hdf5(runpath+str(run_start)+'/'+hdfile,cycle) # load first eigenfunction
|
||||
print 'Read in: '+runpath+str(1)+hdfile+' '+str(tmp.function[0,150,150,0])
|
||||
indices = tmp.function.shape # extent of all dimensions
|
||||
ref = np.zeros((run_end,indices[0],indices[1],indices[2],indices[3])) # initialize ref array
|
||||
ref[0] = tmp.function # set the first run in ref
|
||||
i = run_start + 1
|
||||
while i <= run_end: # begin loop around all runs
|
||||
tmp.read_hdf5(runpath+str(i)+'/'+hdfile,cycle)
|
||||
print 'Read in: '+runpath+str(i)+hdfile+' '+str(tmp.function[0,150,150,0])
|
||||
ref[i - run_start] = tmp.function
|
||||
i += 1
|
||||
meanref = np.average(ref, axis=0) # compute average of all runs
|
||||
EigenFunction.reference = meanref # set to global space in EigenFunction instances
|
||||
return meanref
|
||||
# read in statepoint header data
|
||||
sp = statepoint.StatePoint('statepoint.ref.binary')
|
||||
|
||||
def compute_rms(runlist):
|
||||
# read in results
|
||||
sp.read_results()
|
||||
|
||||
# calculate rms array
|
||||
print 'Calculating rms...'
|
||||
rms = np.zeros(len(runlist))
|
||||
i = 0
|
||||
while i < len(runlist):
|
||||
runlist[i].compute_rms()
|
||||
rms[i] = runlist[i].rms
|
||||
i += 1
|
||||
return rms
|
||||
# extract reference mean
|
||||
mean_ref = extract_mean(sp, tally_id, score_id)
|
||||
|
||||
def plot_rms(rms):
|
||||
print 'Generating plot...'
|
||||
ax = plt.subplot(111)
|
||||
size = rms.shape[0]
|
||||
x = np.linspace(1,size,size)*1e6
|
||||
y = (rms[0]/x[0]**(-0.5))*x**(-0.5)
|
||||
plt.loglog(x,rms*100,'b+')
|
||||
plt.loglog(x,y*100,'g--')
|
||||
ax.xaxis.grid(True,'minor')
|
||||
ax.yaxis.grid(True,'minor')
|
||||
ax.xaxis.grid(True,'major',linewidth=2)
|
||||
ax.yaxis.grid(True,'major',linewidth=2)
|
||||
plt.xlabel('# of Total Neutron Histories (active cycles)')
|
||||
plt.ylabel('RMS Error [%]')
|
||||
plt.legend(('1 million (10 runs)','Ideal Error 1 mil'))
|
||||
return
|
||||
# write gnuplot file
|
||||
write_src_gnuplot('testsrc_pin','Pin mesh',mean_ref,np.size(mean_ref,0))
|
||||
|
||||
# preallocate arrays
|
||||
hists = np.zeros(batch_end - batch_start + 1)
|
||||
norms = np.zeros(batch_end - batch_start + 1)
|
||||
|
||||
i = batch_start
|
||||
while i <= batch_end:
|
||||
|
||||
# process statepoint
|
||||
sp = statepoint.StatePoint('statepoint.'+str(i)+'.binary')
|
||||
sp.read_results()
|
||||
|
||||
# extract mean
|
||||
mean = extract_mean(sp, tally_id, score_id)
|
||||
|
||||
# calculate L2 norm
|
||||
norm = np.linalg.norm(mean - mean_ref)
|
||||
|
||||
# get history information
|
||||
n_inactive = sp.n_inactive
|
||||
current_batch = sp.current_batch
|
||||
n_particles = sp.n_particles
|
||||
gen_per_batch = sp.gen_per_batch
|
||||
n_histories = (current_batch - n_inactive)*n_particles*gen_per_batch
|
||||
|
||||
# batch in vectors
|
||||
hists[i - batch_start] = n_histories
|
||||
norms[i - batch_start] = norm
|
||||
|
||||
# print
|
||||
print 'Batch: '+str(i)+' Histories: '+str(n_histories)+' Norm: '+str(norm)
|
||||
|
||||
i += 1
|
||||
|
||||
# write out gnuplot file
|
||||
write_norm_gnuplot(name,hists,norms,np.size(hists))
|
||||
|
||||
def extract_mean(sp, tally_id,score_id):
|
||||
|
||||
# extract results
|
||||
results = sp.extract_results(tally_id,score_id)
|
||||
|
||||
# extract means and copy
|
||||
mean = results['mean'].copy()
|
||||
|
||||
# reshape and integrate over energy
|
||||
mean = mean.reshape(results['bin_max'],order='F')
|
||||
mean = np.sum(mean,0)
|
||||
mean = np.sum(mean,0)
|
||||
mean = mean/mean.sum()*(mean > 1.e-8).sum()
|
||||
|
||||
return mean
|
||||
|
||||
def write_norm_gnuplot(path,xdat,ydat,size):
|
||||
|
||||
# Header String for GNUPLOT
|
||||
headerstr = """#!/usr/bin/env gnuplot
|
||||
|
||||
set terminal pdf enhanced
|
||||
set output '{output}'
|
||||
set ylabel 'L-2 norm'
|
||||
set xlabel 'Histories'
|
||||
set log x
|
||||
set log y
|
||||
""".format(output=path+'.pdf')
|
||||
|
||||
# Write out the plot string
|
||||
pltstr = "plot '-' using 1:2 with lines"
|
||||
|
||||
# Write out the data string
|
||||
i = 0
|
||||
datastr = ''
|
||||
while i < size:
|
||||
datastr = datastr + '{0} {1}\n'.format(xdat[i],ydat[i])
|
||||
i += 1
|
||||
|
||||
# Concatenate all
|
||||
outstr = headerstr + '\n' + pltstr + '\n' + datastr
|
||||
|
||||
# Write File
|
||||
with open(path+".plot",'w') as f:
|
||||
f.write(outstr)
|
||||
|
||||
# Run GNUPLOT
|
||||
os.system('gnuplot ' + path+".plot")
|
||||
|
||||
def write_src_gnuplot(path,name,src,size):
|
||||
|
||||
# Header String for GNUPLOT
|
||||
headerstr = """#!/usr/bin/env gnuplot
|
||||
|
||||
set terminal pdf enhanced
|
||||
set output '{output}'
|
||||
set palette defined (0 '#000090', 1 '#000fff', 2 '#0090ff', 3 '#0fffee', 4 '#90ff70', 5 '#ffee00', 6 '#ff7000', 7 '#ee0000', 8 '#7f0000')
|
||||
set view map
|
||||
set size ratio -1
|
||||
set lmargin at screen 0.10
|
||||
set rmargin at screen 0.90
|
||||
set bmargin at screen 0.15
|
||||
set tmargin at screen 0.90
|
||||
unset xtics
|
||||
unset ytics
|
||||
set title '{title}'""".format(output=path+'.pdf',title=name)
|
||||
|
||||
# Write out the plot string
|
||||
pltstr = "splot '-' matrix with image "
|
||||
|
||||
# Write out the data string
|
||||
i = 0
|
||||
datastr = ''
|
||||
while i < size:
|
||||
j = 0
|
||||
while j < size:
|
||||
datastr = datastr + '{0} '.format(src[i,j][0])
|
||||
j += 1
|
||||
datastr = datastr + '\n'
|
||||
i += 1
|
||||
|
||||
# replace all nan with zero
|
||||
datastr = datastr.replace('nan','0.0')
|
||||
|
||||
# Concatenate all
|
||||
outstr = headerstr + '\n' + pltstr + '\n' + datastr
|
||||
|
||||
# Write File
|
||||
with open(path+".plot",'w') as f:
|
||||
f.write(outstr)
|
||||
|
||||
# Run GNUPLOT
|
||||
os.system('gnuplot ' + path+".plot")
|
||||
|
||||
def plot_source(source):
|
||||
plt.figure()
|
||||
X = np.linspace(1,272,272)
|
||||
Y = np.linspace(1,272,272)
|
||||
Y,X = np.meshgrid(Y,X)
|
||||
plt.contourf(X,Y,source[0,:,:,0],100)
|
||||
plt.colorbar()
|
||||
plt.xlabel('Mesh Cell in x-direction')
|
||||
plt.ylabel('Mesh Cell in y-direction')
|
||||
plt.title('OPR Converged Fission Source Distribution')
|
||||
plt.show()
|
||||
#
|
||||
if __name__ == "__main__":
|
||||
|
||||
if sys.argv[1] == 'restart':
|
||||
|
||||
# load in data
|
||||
print 'Loading input...'
|
||||
filein = open('rms.out','r')
|
||||
output = cPickle.load(filein)
|
||||
filein.close()
|
||||
rms = output['1milrms']
|
||||
meanref = output['ref']
|
||||
EigenFunction.reference = meanref
|
||||
|
||||
# plot rms
|
||||
plot_rms(rms)
|
||||
|
||||
# plot mean source distribution
|
||||
plot_source(meanref)
|
||||
|
||||
elif sys.argv[1] == 'interactive':
|
||||
|
||||
# load in data
|
||||
print 'Loading input...'
|
||||
filein = open('rms.out','r')
|
||||
output = cPickle.load(filein)
|
||||
filein.close()
|
||||
rms = output['1milrms']
|
||||
meanref = output['ref']
|
||||
EigenFunction.reference = meanref
|
||||
|
||||
# pop an interactive python shell
|
||||
from IPython import embed
|
||||
embed()
|
||||
|
||||
else:
|
||||
|
||||
# calculate reference solution
|
||||
runpath = '/media/Backup/opr_runs/64mil/run'
|
||||
hdfile = 'output.h5'
|
||||
cycle = 210
|
||||
run_start = 1
|
||||
run_end = 4
|
||||
data = 'openmc_src'
|
||||
meanref = create_reference(runpath,hdfile,cycle,run_start,run_end,data)
|
||||
|
||||
# calculate rms for 1 million case
|
||||
runpath = '/media/Backup/opr_runs/1mil/run'
|
||||
hdfile = 'output.h5'
|
||||
cycle_start = 201
|
||||
cycle_end = 840
|
||||
run_start = 1
|
||||
run_end = 10
|
||||
data = 'openmc_src'
|
||||
onemil = read_runs(runpath,hdfile,cycle_start,cycle_end,run_start,run_end,data)
|
||||
|
||||
# calculate rms array
|
||||
rms = compute_rms(onemil)
|
||||
|
||||
# write out numpy array to binary file
|
||||
print 'Writing output...'
|
||||
output = {}
|
||||
output.update({'1milrms':rms})
|
||||
output.update({'ref':meanref})
|
||||
fileout = open('rms.out','wb')
|
||||
cPickle.dump(output,fileout)
|
||||
fileout.close()
|
||||
|
||||
# plot rms
|
||||
plot_rms(rms)
|
||||
|
||||
# plot mean source distribution
|
||||
plot_source(meanref)
|
||||
|
||||
|
||||
tally_id = int(sys.argv[1])
|
||||
score_id = sys.argv[2]
|
||||
batch_start = int(sys.argv[3])
|
||||
batch_end = int(sys.argv[4])
|
||||
name = sys.argv[5]
|
||||
main(tally_id, score_id, batch_start, batch_end, name)
|
||||
|
|
|
|||
|
|
@ -183,10 +183,11 @@ class StatePoint(object):
|
|||
self.k_batch = self._get_double(self.current_batch, path='k_batch')
|
||||
self.entropy = self._get_double(
|
||||
self.current_batch*self.gen_per_batch, path='entropy')
|
||||
self.k_col_abs = self._get_double(path='k_col_abs')[0]
|
||||
self.k_col_tra = self._get_double(path='k_col_tra')[0]
|
||||
self.k_abs_tra = self._get_double(path='k_abs_tra')[0]
|
||||
self.k_combined = self._get_double(2, path='k_combined')
|
||||
if self.revision >= 8:
|
||||
self.k_col_abs = self._get_double(path='k_col_abs')[0]
|
||||
self.k_col_tra = self._get_double(path='k_col_tra')[0]
|
||||
self.k_abs_tra = self._get_double(path='k_abs_tra')[0]
|
||||
self.k_combined = self._get_double(2, path='k_combined')
|
||||
|
||||
# Read number of meshes
|
||||
n_meshes = self._get_int(path='tallies/n_meshes')[0]
|
||||
|
|
@ -439,6 +440,115 @@ class StatePoint(object):
|
|||
# has been called already.
|
||||
return t.results[filter_index, score_index]
|
||||
|
||||
def extract_results(self, tally_id, score_str):
|
||||
"""Returns a tally results dictionary given a tally_id and score string.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
tally_id : int
|
||||
Index for the tally in StatePoint.tallies list
|
||||
|
||||
score_str : string
|
||||
Corresponds to the string entered for a score in tallies.xml.
|
||||
For a flux score extraction it would be 'score'
|
||||
|
||||
"""
|
||||
|
||||
# get tally
|
||||
try:
|
||||
tally = self.tallies[tally_id-1]
|
||||
except:
|
||||
print 'Tally does not exist'
|
||||
return
|
||||
|
||||
# get the score index if it is present
|
||||
try:
|
||||
idx = tally.scores.index(score_str)
|
||||
except ValueError:
|
||||
print 'Score does not exist'
|
||||
print tally.scores
|
||||
return
|
||||
|
||||
# create numpy array for mean and 95% CI
|
||||
n_bins = len(tally.results)
|
||||
n_filters = len(tally.filters)
|
||||
n_scores = len(tally.scores)
|
||||
meanv = np.zeros(n_bins)
|
||||
unctv = np.zeros(n_bins)
|
||||
filters = np.zeros((n_bins,n_filters))
|
||||
filtmax = np.zeros(n_filters+1)
|
||||
meshmax = np.zeros(4)
|
||||
filtmax[0] = 1
|
||||
meshmax[0] = 1
|
||||
|
||||
# get number of realizations
|
||||
n = tally.n_realizations
|
||||
|
||||
# get t-value
|
||||
t_value = scipy.stats.t.ppf(0.975, n - 1)
|
||||
|
||||
# calculate mean
|
||||
meanv = tally.results[:,idx,0]
|
||||
meanv = meanv / n
|
||||
|
||||
# calculate 95% two-sided CI
|
||||
unctv = tally.results[:,idx,1]
|
||||
unctv = t_value*np.sqrt((unctv/n - meanv*meanv)/(n-1))/meanv
|
||||
|
||||
# create output dictionary
|
||||
data = {'mean':meanv,'CI95':unctv}
|
||||
|
||||
# get bounds of filter bins
|
||||
for akey in tally.filters.keys():
|
||||
idx = tally.filters.keys().index(akey)
|
||||
filtmax[n_filters - idx] = tally.filters[akey].length
|
||||
|
||||
# compute bin info
|
||||
for i in range(n_filters):
|
||||
|
||||
# compute indices for filter combination
|
||||
filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) %
|
||||
np.prod(filtmax[0:i+2]))/(np.prod(filtmax[0:i+1]))) + 1
|
||||
|
||||
# append in dictionary bin with filter
|
||||
data.update({tally.filters.keys()[n_filters - i - 1]:
|
||||
filters[:,n_filters - i - 1]})
|
||||
|
||||
# check for mesh
|
||||
if tally.filters.keys()[n_filters - i - 1] == 'mesh':
|
||||
dims = list(self.meshes[tally.filters['mesh'].bins[0] - 1].dimension)
|
||||
dims.reverse()
|
||||
dims = np.asarray(dims)
|
||||
if score_str == 'current':
|
||||
dims += 1
|
||||
meshmax[1:4] = dims
|
||||
mesh_bins = np.zeros((n_bins,3))
|
||||
mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) %
|
||||
np.prod(meshmax[0:2]))/(np.prod(meshmax[0:1]))) + 1
|
||||
mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) %
|
||||
np.prod(meshmax[0:3]))/(np.prod(meshmax[0:2]))) + 1
|
||||
mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) %
|
||||
np.prod(meshmax[0:4]))/(np.prod(meshmax[0:3]))) + 1
|
||||
data.update({'mesh':zip(mesh_bins[:,0],mesh_bins[:,1],
|
||||
mesh_bins[:,2])})
|
||||
i += 1
|
||||
|
||||
# add in maximum bin filters and order
|
||||
b = tally.filters.keys()
|
||||
b.reverse()
|
||||
filtmax = list(filtmax[1:])
|
||||
try:
|
||||
idx = b.index('mesh')
|
||||
filtmax[idx] = np.max(mesh_bins[:,2])
|
||||
filtmax.insert(idx,np.max(mesh_bins[:,1]))
|
||||
filtmax.insert(idx,np.max(mesh_bins[:,0]))
|
||||
|
||||
except ValueError:
|
||||
pass
|
||||
data.update({'bin_order':b,'bin_max':filtmax})
|
||||
|
||||
return data
|
||||
|
||||
def _get_data(self, n, typeCode, size):
|
||||
return list(struct.unpack('={0}{1}'.format(n,typeCode),
|
||||
self._f.read(n*size)))
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue