ESS-Estia/Estia-MCNP
1
0
Fork 0
Code Actions Packages Releases Activity

New improved plotting and add overview drawing svg and complete gimp file

Browse Source
This commit is contained in:
glavic_a
2019-05-03 18:31:33 +02:00
parent 83511136e1
commit 6e282839c3
3 changed files with 159 additions and 53 deletions
Download Patch File Download Diff File Expand all files Collapse all files

200
plot_tally.py
View File

@@ -12,6 +12,7 @@ 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
@@ -40,6 +41,11 @@ 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]
def read_tally(fname, use_tally=None):
if fname.endswith('.gz'):
txt=gzip.open(fname, 'r').read()
@@ -81,80 +87,102 @@ def read_tally(fname, use_tally=None):
def plot_xy(outfile, res, zidx=5):
x,y,z,I,dI=res
fig=gcf()
ax=gca()
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)
ax.set_aspect('equal')
fig.colorbar(p, orientation='horizontal', shrink=0.9, spacing='proportional',
label='Dose rate [Sv/h]',
format="10$^{%i}$")
xlabel('x [m]')
ylabel('y [m]')
title('XY plane at z=%.1fm'%((z[zidx]+z[zidx+1])/200.))
fig.subplots_adjust()
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))
if outfile is not None:
fig.savefig(outfile, transparent=False, dpi=300)
#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
fig=gcf()
ax=gca()
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)
ax.set_aspect('equal')
fig.colorbar(p, orientation='horizontal', shrink=0.9, spacing='proportional',
label='Dose rate [Sv/h]',
format="10$^{%i}$")
xlabel('x [m]')
ylabel('z [m]')
ypos=(y[yidx]+y[yidx+1])/200.
title('XZ plane at y=%.1fm'%ypos)
fig.subplots_adjust()
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)
if outfile is not None:
fig.savefig(outfile, transparent=False, dpi=300)
#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
fig=gcf()
ax=gca()
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)
ax.set_aspect('equal')
fig.colorbar(p, orientation='horizontal', shrink=0.9, spacing='proportional',
label='Dose rate [Sv/h]',
format="10$^{%i}$")
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.))
fig.subplots_adjust()
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.)
if outfile is not None:
fig.savefig(outfile, transparent=False, dpi=300)
fig.savefig(outfile, transparent=False, dpi=300)
def get_name(fname):
# extract tag-name and version info from filename
@@ -171,6 +199,17 @@ if __name__=='__main__':
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)
@@ -180,7 +219,12 @@ if __name__=='__main__':
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
@@ -195,25 +239,75 @@ if __name__=='__main__':
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)
xi, yi, zi=res[3].shape
for zidx in range(zi):
print fpre+'_xy_%02i.png'%zidx
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_%02i.png'%zidx, res, zidx=zidx)
for yidx in range(yi):
print fpre+'_xz_%02i.png'%yidx
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_%02i.png'%yidx, res, yidx=yidx)
for xidx in range(xi):
print fpre+'_yz_%02i.png'%xidx
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_%02i.png'%xidx, res, xidx=xidx)
plot_yz(fpre+'_yz_%+05i.png'%xi, res, xidx=xidx)