325 lines
11 KiB
Python
325 lines
11 KiB
Python
#!/usr/bin/env python
|
|
#-*- coding: utf-8 -*-
|
|
"""
|
|
Generate nice graphs from mesh-tally 1 mctal files.
|
|
|
|
At the current stage only simple text processing is done that doesn't work with
|
|
multiple tallies or energy dependent meshes, yet.
|
|
"""
|
|
|
|
import gzip
|
|
from numpy import *
|
|
from pylab import *
|
|
from matplotlib.colors import LogNorm, BoundaryNorm, ListedColormap
|
|
from matplotlib.ticker import LogFormatterMathtext, LogLocator
|
|
from mpl_toolkits.axes_grid1 import make_axes_locatable
|
|
|
|
COLORMAP='tab20c'
|
|
SCALING=1.0
|
|
|
|
# generate colormap to emphasize radiation levels properly
|
|
Rcref=['#3182BD', '#6BAED6', '#9ECAE1', '#C6DBEF',
|
|
'#31A354', '#74C476', '#A1D99B', '#C7E9C0',
|
|
'#756BB1', '#9E9AC8', '#BCBDDC', '#DADAEB',
|
|
'#E6550D', '#FD8D3C', '#FDAE6B', '#FDD0A2',
|
|
'#000000', '#111111', '#222222', '#333333',
|
|
'#444444', '#555555', '#666666', '#777777',
|
|
'#888888', '#999999', '#aaaaaa', '#bbbbbb',
|
|
'#cccccc', '#dddddd', '#eeeeee', '#ffffff',
|
|
]
|
|
Rbounds=[
|
|
1e-7, 2.5e-7, 5e-7, 7.5e-7,
|
|
1e-6, 2.5e-6, 5e-6, 7.5e-6,
|
|
1e-5, 2.5e-5, 5e-5, 7.5e-5,
|
|
1e-4, 2.5e-4, 5e-4, 7.5e-4,
|
|
1e-3, 2.5e-3, 5e-3, 7.5e-3,
|
|
1e-2, 2.5e-2, 5e-2, 7.5e-2,
|
|
1e-1, 2.5e-1, 5e-1, 7.5e-1,
|
|
1, 2.5, 5, 7.5, 10
|
|
]
|
|
Rbounds=[log10(bi) for bi in Rbounds]
|
|
Rcmap=ListedColormap(Rcref)
|
|
Rnorm=BoundaryNorm(Rbounds, Rcmap.N)
|
|
|
|
XY_PLANES=[0.]
|
|
XZ_PLANES=[0.]
|
|
YZ_PLANES=[1010, 1100, 1300, 1550,
|
|
2050, 2300, 2900, 3500]
|
|
|
|
try:
|
|
# image data generated from MCNP plot output for geometry lines
|
|
geometry=load('/home/glavic_a/estia/media/data and simulation/shielding/e2-estia_baseline_v07-0-gb8bc945_lines.npz')
|
|
except:
|
|
geometry=None
|
|
|
|
def read_tally(fname, use_tally=None):
|
|
if fname.endswith('.gz'):
|
|
txt=gzip.open(fname, 'r').read()
|
|
else:
|
|
txt=open(fname, 'r').read()
|
|
|
|
tally_ranges=[]
|
|
tidx=0
|
|
while 'tally' in txt[tidx+1:]:
|
|
tidx=txt.index('tally ', tidx+1)
|
|
tnr=int(txt[tidx:].split(None,3)[1])
|
|
if len(tally_ranges)>0:
|
|
tally_ranges[-1][2]=tidx-1
|
|
tally_ranges.append([tnr, tidx, -1])
|
|
if use_tally is None:
|
|
txt=txt[tally_ranges[0][1]:tally_ranges[0][2]]
|
|
else:
|
|
tidx=[ti[0] for ti in tally_ranges].index(use_tally)
|
|
txt=txt[tally_ranges[tidx][1]:tally_ranges[tidx][2]]
|
|
header, data=txt.split('vals')
|
|
|
|
|
|
# get dimensions of data
|
|
fstart=header.index('\nf ')
|
|
hi=header[fstart+3:].splitlines()[0]
|
|
N,nE,nX,nY,nZ=map(int, hi.strip().split())
|
|
|
|
grid=header[fstart+3+len(hi):header.index('\nd ')]
|
|
grid=array(grid.replace('\n', '').strip().split(), dtype=float)
|
|
|
|
x=grid[:nX+1]
|
|
y=grid[nX+1:nX+nY+2]
|
|
z=grid[nX+nY+2:nX+nY+nZ+3]
|
|
|
|
data=array(data.replace('\n','').strip().split(), dtype=float)
|
|
I,dI=data.reshape(nZ, nY, nX, 2).transpose(3,2,1,0)
|
|
|
|
return x,y,z,I,dI
|
|
|
|
def plot_xy(outfile, res, zidx=5):
|
|
x,y,z,I,dI=res
|
|
xs=x[-1]-x[0]
|
|
ys=y[-1]-y[0]
|
|
fig=figure(figsize=(xs*0.00394+0.5, ys*0.00394/0.95+1.0))
|
|
rel_w, rel_h= (xs*0.00394/(xs*0.00394+0.5),
|
|
ys*0.00394/(ys*0.00394/0.95+1.0))
|
|
position=[(1.-rel_w)/2., (1.-rel_h)-0.05, rel_w, rel_h]
|
|
ax=axes(position, aspect='equal')
|
|
p=ax.pcolormesh(x/100., y/100.,
|
|
log10(I[:,:,zidx].T+1e-8), norm=Rnorm, cmap=Rcmap)
|
|
xlabel('x [m]')
|
|
ylabel('y [m]')
|
|
title('XY plane at z=%.1fm'%((z[zidx]+z[zidx+1])/200.))
|
|
if geometry is not None:
|
|
ax.pcolormesh(geometry['x'], geometry['y'],
|
|
log10(geometry['lines']*0+2e-3),
|
|
norm=Rnorm, cmap=Rcmap
|
|
)
|
|
|
|
divider = make_axes_locatable(ax)
|
|
cax = divider.append_axes("bottom", size="5%", pad=0.6)
|
|
cbar=fig.colorbar(p, orientation='horizontal', spacing='proportional',
|
|
label='Dose rate [Sv/h]', cax=cax,
|
|
#fraction=0.1, pad=0.1,
|
|
#aspect=50.,
|
|
format="10$^{%i}$")
|
|
#fig.savefig(outfile)
|
|
#pc=p.get_facecolors()
|
|
#pc[:,3]=where(I[:,:,zidx].T.flatten()>0, 0.9-0.25*dI[:,:,zidx].T.flatten(), 0.)
|
|
#ax.add_patch(Circle((-0.151,0.021), radius=5.50, color='black', fill=False, lw=2))
|
|
#ax.add_patch(Circle((-0.151,0.021), radius=11.50, color='black', fill=False, lw=2))
|
|
#ax.add_patch(Circle((-0.151,0.021), radius=15.00, color='black', fill=False, lw=2))
|
|
## shielding walls outside bunker
|
|
#ax.add_line(Line2D([21.12-0.151, 21.12-0.151], [-1.5,1.5], color='black', lw=2))
|
|
#ax.add_line(Line2D([14.845,19.426], [-0.541,-0.713], color='black', lw=2))
|
|
#ax.add_line(Line2D([19.426,19.405], [-0.713,-1.248], color='black', lw=2))
|
|
#ax.add_line(Line2D([19.405,21.205], [-1.248,-1.460], color='black', lw=2))
|
|
#ax.add_line(Line2D([14.723, 22.740], [2.232, 3.084], color='black', lw=2))
|
|
#ax.add_line(Line2D([14.801, 20.894], [-1.210, -1.724], color='black', lw=2))
|
|
fig.savefig(outfile, transparent=False, dpi=300)
|
|
|
|
def plot_xz(outfile, res, yidx=5):
|
|
x,y,z,I,dI=res
|
|
xs=x[-1]-x[0]
|
|
ys=z[-1]-z[0]
|
|
fig=figure(figsize=(xs*0.00394+0.5, ys*0.00394/0.95+1.0))
|
|
rel_w, rel_h= (xs*0.00394/(xs*0.00394+0.5),
|
|
ys*0.00394/(ys*0.00394/0.95+1.0))
|
|
position=[(1.-rel_w)/2., (1.-rel_h)-0.08, rel_w, rel_h]
|
|
ax=axes(position, aspect='equal')
|
|
p=ax.pcolormesh(x/100., z/100.,
|
|
log10(I[:,yidx].T+1e-8), norm=Rnorm, cmap=Rcmap)
|
|
xlabel('x [m]')
|
|
ylabel('z [m]')
|
|
ypos=(y[yidx]+y[yidx+1])/200.
|
|
title('XZ plane at y=%.1fm'%ypos)
|
|
|
|
divider = make_axes_locatable(ax)
|
|
cax = divider.append_axes("bottom", size="5%", pad=0.6)
|
|
cbar=fig.colorbar(p, orientation='horizontal', spacing='proportional',
|
|
label='Dose rate [Sv/h]', cax=cax,
|
|
#fraction=0.1, pad=0.1,
|
|
#aspect=50.,
|
|
format="10$^{%i}$")
|
|
#fig.savefig(outfile)
|
|
#pc=p.get_facecolors()
|
|
#pc[:,3]=where(I[:,yidx].T.flatten()>0, 0.9-0.25*dI[:,yidx].T.flatten(), 0.)
|
|
#ax.plot([sqrt(5.5**2-ypos**2)-0.151, sqrt(5.5**2-ypos**2)-0.151],
|
|
#[z.min()/100., z.max()/100.], color='black', lw=2)
|
|
#ax.plot([sqrt(11.5**2-ypos**2)-0.151, sqrt(11.5**2-ypos**2)-0.151],
|
|
#[z.min()/100., z.max()/100.], color='black', lw=2)
|
|
#ax.plot([sqrt(15.0**2-ypos**2)-0.151, sqrt(15.0**2-ypos**2)-0.151],
|
|
#[z.min()/100., z.max()/100.], color='black', lw=2)
|
|
fig.savefig(outfile, transparent=False, dpi=300)
|
|
|
|
def plot_yz(outfile, res, xidx=5):
|
|
x,y,z,I,dI=res
|
|
xs=y[-1]-y[0]
|
|
ys=z[-1]-z[0]
|
|
fig=figure(figsize=(xs*0.00394+0.5, ys*0.00394/0.95+1.0))
|
|
rel_w, rel_h= (xs*0.00394/(xs*0.00394+0.5),
|
|
ys*0.00394/(ys*0.00394/0.95+1.0))
|
|
position=[(1.-rel_w)/2., (1.-rel_h)-0.05, rel_w, rel_h]
|
|
ax=axes(position, aspect='equal')
|
|
p=ax.pcolormesh(y/100., z/100.,
|
|
log10(I[xidx].T+1e-8), norm=Rnorm, cmap=Rcmap)
|
|
xlim(y.max()/100., y.min()/100.)
|
|
xlabel('y [m]')
|
|
ylabel('z [m]')
|
|
title('YZ plane at x=%.1fm'%((x[xidx]+x[xidx+1])/200.))
|
|
|
|
divider = make_axes_locatable(ax)
|
|
cax = divider.append_axes("bottom", size="5%", pad=0.6)
|
|
cbar=fig.colorbar(p, orientation='horizontal', spacing='proportional',
|
|
label='Dose rate [Sv/h]', cax=cax,
|
|
#fraction=0.1, pad=0.1,
|
|
#aspect=50.,
|
|
format="10$^{%i}$")
|
|
#fig.savefig(outfile)
|
|
#pc=p.get_facecolors()
|
|
#pc[:,3]=where(I[xidx].T.flatten()>0, 0.9-0.25*dI[xidx].T.flatten(), 0.)
|
|
fig.savefig(outfile, transparent=False, dpi=300)
|
|
|
|
def get_name(fname):
|
|
# extract tag-name and version info from filename
|
|
stripped=os.path.basename(fname)[9:].split('.')[0]
|
|
tag, version=stripped.split('-', 1)
|
|
return tag, version
|
|
|
|
|
|
if __name__=='__main__':
|
|
import sys, os
|
|
if '-o' in sys.argv:
|
|
idx=sys.argv.index('-o')
|
|
sys.argv.pop(idx)
|
|
prefix=sys.argv.pop(idx)
|
|
else:
|
|
prefix='estia_dosemap'
|
|
if '-x' in sys.argv:
|
|
idx=sys.argv.index('-x')
|
|
sys.argv.pop(idx)
|
|
output_path=os.path.join('.', 'images', sys.argv.pop(idx))
|
|
else:
|
|
output_path=None
|
|
maxadd=[]
|
|
while '-m' in sys.argv:
|
|
idx=sys.argv.index('-m')
|
|
sys.argv.pop(idx)
|
|
maxadd.append(sys.argv.pop(idx))
|
|
print sys.argv[1]
|
|
if ':' in sys.argv[1]:
|
|
fname, tally=sys.argv[1].split(':',2)
|
|
tally=int(tally)
|
|
else:
|
|
fname=sys.argv[1]
|
|
tally=None
|
|
res=read_tally(fname, use_tally=tally)
|
|
x,y,z,I,dI=res
|
|
base_points=array([pi.flatten() for pi in
|
|
meshgrid((y[:-1]+y[1:])/2.,
|
|
(x[:-1]+x[1:])/2.,
|
|
(z[:-1]+z[1:])/2.)])[[1,0,2]].T
|
|
I*=SCALING
|
|
dI*=SCALING
|
|
fpath=os.path.join('.', 'images', get_name(fname)[0]+'-'+get_name(fname)[1])
|
|
if tally is not None:
|
|
fpath+='[%s]'%tally
|
|
if len(sys.argv)>2:
|
|
for fi in sys.argv[2:]:
|
|
print fi
|
|
if ':' in fi:
|
|
fi, tally=fi.split(':',2)
|
|
tally=int(tally)
|
|
else:
|
|
tally=None
|
|
resi=read_tally(fi, use_tally=tally)
|
|
xi,yi,zi,Ii,dIi=resi
|
|
Ii*=SCALING
|
|
dI*=SCALING
|
|
if Ii.shape!=I.shape:
|
|
print "different shape, interpolate to existing grid..."
|
|
from scipy.interpolate import interpn
|
|
interp=interpn(((xi[:-1]+xi[1:])/2.,
|
|
(yi[:-1]+yi[1:])/2.,
|
|
(zi[:-1]+zi[1:])/2.), Ii, base_points,
|
|
fill_value=0., bounds_error=False)
|
|
dinterp=interpn(((xi[:-1]+xi[1:])/2.,
|
|
(yi[:-1]+yi[1:])/2.,
|
|
(zi[:-1]+zi[1:])/2.), dIi, base_points,
|
|
fill_value=0., bounds_error=False)
|
|
#print(xi.mean()-x.mean(), yi.mean()-y.mean(),
|
|
#zi.mean()-z.mean(), Ii.max(), interp.max())
|
|
Ii=interp.reshape(*I.shape)
|
|
dIi=dinterp.reshape(*I.shape)
|
|
I+=Ii
|
|
dI=sqrt(dI**2+dIi**2)
|
|
res=(x,y,z,I,dI)
|
|
fpath+='+'+get_name(fi)[1]
|
|
if tally is not None:
|
|
fpath+='[%s]'%tally
|
|
for fi in maxadd:
|
|
print "create maximum of other data and", fi
|
|
if ':' in fi:
|
|
fi, tally=fi.split(':',2)
|
|
tally=int(tally)
|
|
else:
|
|
tally=None
|
|
resi=read_tally(fi, use_tally=tally)
|
|
xi,yi,zi,Ii,dIi=resi
|
|
Ii*=SCALING
|
|
dI*=SCALING
|
|
if Ii.shape!=I.shape:
|
|
print "different shape, interpolate to existing grid..."
|
|
from scipy.interpolate import interpn
|
|
interp=interpn(((xi[:-1]+xi[1:])/2.,
|
|
(yi[:-1]+yi[1:])/2.,
|
|
(zi[:-1]+zi[1:])/2.), Ii, base_points,
|
|
fill_value=0., bounds_error=False)
|
|
dinterp=interpn(((xi[:-1]+xi[1:])/2.,
|
|
(yi[:-1]+yi[1:])/2.,
|
|
(zi[:-1]+zi[1:])/2.), dIi, base_points,
|
|
fill_value=0., bounds_error=False)
|
|
Ii=interp.reshape(*I.shape)
|
|
dIi=dinterp.reshape(*I.shape)
|
|
dI=where(I>Ii, dI, dIi)
|
|
I=maximum(I, Ii)
|
|
res=(x,y,z,I,dI)
|
|
fpath+='+'+get_name(fi)[1]
|
|
if tally is not None:
|
|
fpath+='[%s]'%tally
|
|
if output_path is not None:
|
|
fpath=output_path
|
|
os.mkdir(fpath)
|
|
fpre=os.path.join(fpath, prefix)
|
|
for zi in XY_PLANES:
|
|
zidx=where(z>zi)[0][0]
|
|
print fpre+'_xy_%+05i.png'%zi
|
|
figure(figsize=(12,8))
|
|
plot_xy(fpre+'_xy_%+05i.png'%zi, res, zidx=zidx)
|
|
for yi in XZ_PLANES:
|
|
yidx=where(y>yi)[0][0]
|
|
print fpre+'_xz_%+05i.png'%yi
|
|
figure(figsize=(12,8))
|
|
plot_xz(fpre+'_xz_%+05i.png'%yi, res, yidx=yidx)
|
|
for xi in YZ_PLANES:
|
|
xidx=where(x>xi)[0][0]
|
|
print fpre+'_yz_%+05i.png'%xi
|
|
figure(figsize=(8,10))
|
|
plot_yz(fpre+'_yz_%+05i.png'%xi, res, xidx=xidx)
|
|
|