mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
modified post processing script to work for any piece of output data
This commit is contained in:
Bryan Herman
parent
73fd4449df
commit
4ea142cc27
1 changed files with 141 additions and 79 deletions
|
|
@ -12,20 +12,22 @@ class EigenFunction:
|
|||
'''Represents the reference eigenfunction'''
|
||||
reference = 0.
|
||||
#
|
||||
def __init__(self):
|
||||
def __init__(self,data):
|
||||
#
|
||||
'''Initializes the eigenfunction'''
|
||||
self.function = 0.
|
||||
self.data = data
|
||||
#
|
||||
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)+'/openmc_src'
|
||||
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 * ((16)**2-5*2**2)*177
|
||||
self.iamref = 'F'
|
||||
#
|
||||
def set_reference(self):
|
||||
|
|
@ -38,102 +40,162 @@ class EigenFunction:
|
|||
#
|
||||
'''Computes RMS value'''
|
||||
Np = self.function.size
|
||||
tmp = (self.function - EigenFunction.reference)*(self.function - EigenFunction.reference)
|
||||
Np = ((16)**2-5*2**2)*177
|
||||
tmp = (self.function - EigenFunction.reference)**2
|
||||
tmp2 = tmp.sum()
|
||||
self.rms = np.sqrt((1.0/float(Np))*tmp2)
|
||||
#
|
||||
if sys.argv[1] != 'restart':
|
||||
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
|
||||
j += 1
|
||||
meantmp.function = np.average(runs, axis=0) # compute the mean
|
||||
runlist.append(meantmp)
|
||||
print 'Read in from path: '+runpath+' Cycle: '+str(i)
|
||||
i += 1
|
||||
return runlist
|
||||
|
||||
# calculate reference solution
|
||||
def create_reference(runpath,hdfile,cycle,run_start,run_end,data):
|
||||
|
||||
# calculate reference solution
|
||||
print 'Calculating Reference solution...'
|
||||
runpath = '/media/Backup/opr_runs/1mil/run' # the directory prefix
|
||||
hdfile = '/output.h5' # hdf5 file name
|
||||
cycle = 840 # cycle number to extract
|
||||
run_end = 25 # number of runs to look at
|
||||
tmp = EigenFunction()
|
||||
tmp.read_hdf5(runpath+str(1)+hdfile,cycle) # load first eigenfunction
|
||||
tmp = EigenFunction(data)
|
||||
tmp.read_hdf5(runpath+str(run_start)+'/'+hdfile,cycle) # load first eigenfunction
|
||||
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 = 2
|
||||
i = run_start + 1
|
||||
while i <= run_end: # begin loop around all runs
|
||||
tmp.read_hdf5(runpath+str(i)+hdfile,cycle)
|
||||
tmp.read_hdf5(runpath+str(i)+'/'+hdfile,cycle)
|
||||
print 'Read in: '+runpath+str(i)+hdfile
|
||||
ref[i-1] = tmp.function
|
||||
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
|
||||
|
||||
# calculate rms for 1 million case
|
||||
print 'Calculating 1million History case...'
|
||||
onemil = []
|
||||
cycle_start = 201
|
||||
cycle_end = 840
|
||||
run_end = 10
|
||||
i = cycle_start
|
||||
while i <= cycle_end:
|
||||
meantmp = EigenFunction() # init mean object
|
||||
j = 1
|
||||
runs = np.zeros((run_end,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-1] = tmp.function # put function into runs
|
||||
j += 1
|
||||
meantmp.function = np.average(runs, axis=0) # compute the mean
|
||||
onemil.append(meantmp)
|
||||
print 'Read in from path: '+runpath+' Cycle: '+str(i)
|
||||
i += 1
|
||||
def compute_rms(runlist):
|
||||
|
||||
# calculate rms array
|
||||
print 'Calculating rms...'
|
||||
rms = np.zeros((cycle_end - cycle_start + 1))
|
||||
rms = np.zeros(len(runlist))
|
||||
i = 0
|
||||
while i < len(onemil):
|
||||
onemil[i].compute_rms()
|
||||
rms[i] = onemil[i].rms
|
||||
while i < len(runlist):
|
||||
runlist[i].compute_rms()
|
||||
rms[i] = runlist[i].rms
|
||||
i += 1
|
||||
return rms
|
||||
|
||||
# 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()
|
||||
else:
|
||||
# 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
|
||||
def plot_rms(rms):
|
||||
print 'Generating plot...'
|
||||
ax = plt.subplot(111)
|
||||
size = rms.shape[0]
|
||||
x = np.linspace(1,size,size)*1e6
|
||||
y = (rms[size-1]/x[size-1]**(-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
|
||||
|
||||
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/1mil/run'
|
||||
hdfile = 'output.h5'
|
||||
cycle = 840
|
||||
run_start = 1
|
||||
run_end = 25
|
||||
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 = 1
|
||||
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)
|
||||
|
||||
# plot rms
|
||||
print 'Generating plot...'
|
||||
ax = plt.subplot(111)
|
||||
x = np.linspace(1,640,640)*1e6
|
||||
y = rms[0]/1e-3*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'))
|
||||
|
||||
# plot mean source distribution
|
||||
plt.figure()
|
||||
X = np.linspace(1,272,272)
|
||||
Y = np.linspace(1,272,272)
|
||||
Y,X = np.meshgrid(Y,X)
|
||||
plt.contourf(X,Y,meanref[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()
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue