#!/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)