Add mesh plotting and run with WW+ExpT and tallies at E01+E03

e2-estia.i

@@ -732,6 +732,24 @@ tmesh
 endmd
 fm1      8.7374e15 
 c
+fc44     neutron flux outside bunker wall E01 and E03
+f44:n    708 709
+e44      1e-10 1.58e-10 2.51e-10 3.98e-10 6.31e-10
+         1e-9 1.58e-9 2.51e-9 3.98e-9 6.31e-9 
+         1e-8 1.58e-8 2.51e-8 3.98e-8 6.31e-8 
+         1e-7 1.58e-7 2.51e-7 3.98e-7 6.31e-7 
+         1e-6 1.58e-6 2.51e-6 3.98e-6 6.31e-6 
+         1e-5 1.58e-5 2.51e-5 3.98e-5 6.31e-5 
+         1e-4 1.58e-4 2.51e-4 3.98e-4 6.31e-4 
+         1e-3 1.58e-3 2.51e-3 3.98e-3 6.31e-3 
+         1e-2 1.58e-2 2.51e-2 3.98e-2 6.31e-2 
+         1e-1 1.58e-1 2.51e-1 3.98e-1 6.31e-1 
+         1e-0   1.58    2.51    3.98    6.31 
+         1e+1 1.58e+1 2.51e+1 3.98e+1 6.31e+1 
+         1e+2 1.58e+2 2.51e+2 3.98e+2 6.31e+2
+         1e+3 1.58e+3 2.51e+3
+fq44 e f
+fm44      8.7374e15 
 c 
 c ========================================== 
 c
@@ -762,7 +780,7 @@ c -------------------------------------------------------
 c --------------- PHYSICS CARDS --------------------------
 c -------------------------------------------------------
 c dbcn 375642321 j j j j j j j j j j j 15291711
-nps 1e9
+nps 3e9
 c histp
 mode n p n h / d t s a
 c calculate volume for tally cells (spheres)

plot_tally.py

@@ -0,0 +1,127 @@
+#!/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
+
+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, Imin=1e-6, Imax=1e4, cticks=None):
+  x,y,z,I,dI=res
+  fig=gcf()
+  ax=gca()
+  p=ax.pcolormesh(x/100., y/100., I[:,:,zidx].T, norm=LogNorm(Imin, Imax), cmap='gist_ncar')
+  ax.set_aspect('equal')
+  fig.colorbar(p, orientation='horizontal', shrink=0.9, ticks=cticks, label='Dose rate [Sv/h]')
+  xlabel('x [m]')
+  ylabel('y [m]')
+  title('XY plane at z=%.1fm'%((z[zidx]+z[zidx+1])/200.))
+  #clabel('I [Sv/h]')
+  fig.subplots_adjust()
+  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,0), radius=5.50, color='black', fill=False, lw=2))
+  ax.add_patch(Circle((0,0), radius=11.50, color='black', fill=False, lw=2))
+  ax.add_patch(Circle((0,0), radius=15.00, color='black', fill=False, lw=2))
+  fig.savefig(outfile, transparent=True, dpi=300)
+
+def plot_xz(outfile, res, yidx=5, Imin=1e-6, Imax=1e4, cticks=None):
+  x,y,z,I,dI=res
+  fig=gcf()
+  ax=gca()
+  p=ax.pcolormesh(x/100., z/100., I[:,yidx].T, norm=LogNorm(Imin, Imax), cmap='gist_ncar')
+  ax.set_aspect('equal')
+  fig.colorbar(p, orientation='horizontal', shrink=0.9, ticks=cticks, label='Dose rate [Sv/h]')
+  xlabel('x [m]')
+  ylabel('z [m]')
+  title('XZ plane at y=%.1fm'%((y[yidx]+y[yidx+1])/200.))
+  #clabel('I [Sv/h]')
+  fig.subplots_adjust()
+  fig.savefig(outfile)
+  pc=p.get_facecolors()
+  pc[:,3]=where(I[:,yidx].T.flatten()>0, 0.9-0.25*dI[:,yidx].T.flatten(), 0.)
+  fig.savefig(outfile, transparent=True, dpi=300)
+
+def plot_yz(outfile, res, xidx=5, Imin=1e-6, Imax=1e4, cticks=None):
+  x,y,z,I,dI=res
+  fig=gcf()
+  ax=gca()
+  p=ax.pcolormesh(y/100., z/100., I[xidx].T, norm=LogNorm(Imin, Imax), cmap='gist_ncar')
+  ax.set_aspect('equal')
+  fig.colorbar(p, orientation='horizontal', shrink=0.9, ticks=cticks, label='Dose rate [Sv/h]')
+  xlabel('y [m]')
+  ylabel('z [m]')
+  title('YZ plane at x=%.1fm'%((x[xidx]+x[xidx+1])/200.))
+  #clabel('I [Sv/h]')
+  fig.subplots_adjust()
+  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=True, dpi=300)
+
+if __name__=='__main__':
+  import sys, os
+  res=read_tally(sys.argv[1], use_tally=1)
+  x,y,z,I,dI=res
+  I*=3600.
+  fpre=os.path.join('.', 'images', os.path.basename(sys.argv[1])[:-4])
+  xi, yi, zi=res[3].shape
+  for zidx in range(zi):
+    print fpre+'_xy_%02i.png'%zidx
+    figure(figsize=(12,8))
+    plot_xy(fpre+'_xy_%02i.png'%zidx, res, zidx=zidx, cticks=[1e-6, 1e-3, 1, 1e3])
+  for yidx in range(yi):
+    print fpre+'_xz_%02i.png'%yidx
+    figure(figsize=(12,8))
+    plot_xz(fpre+'_xz_%02i.png'%yidx, res, yidx=yidx, cticks=[1e-6, 1e-3, 1, 1e3])
+  for xidx in range(xi):
+    print fpre+'_yz_%02i.png'%xidx
+    figure(figsize=(8,10))
+    plot_yz(fpre+'_yz_%02i.png'%xidx, res, xidx=xidx, cticks=[1e-6, 1e-3, 1, 1e3])
+