3 from __future__
import print_function
5 if sys.version_info[0]!=2:
6 from importlib
import reload
13 import dispatch_data
as dd
18 Main methods: patches, cache 19 patches() Return list of objects for all 'data/run' patches 20 cache() Adds cached arrays to the object 22 pp=patches(iout,'run') 23 pp=[p for p in pp if select(p)] 26 ... do something with the data in p.d, p.ux, ... 28 patch_files() Return list of patch file root names 29 p=patch() Return one patch object 31 print(
"Type dispatch_utils?<RETURN> for help text, or dispach_utils(\ 32 for pop-up help in Canopy")
34 def snapshot (iout=0, run='', data='../data', overlap=0.01, verbose=0):
35 snap=dd.Snapshot (iout, run, data, overlap=overlap, verbose=verbose)
40 for p
in snap.patches:
45 ''' Return extent limits in x,y,z,values ''' 46 def __init__(self,pp,ix=None,iy=None,iz=None,iv=0,verbose=0):
52 xr=p.active_corners[:,0]
55 yr=p.active_corners[:,1]
58 zr=p.active_corners[:,2]
62 self.
vlim[0]=np.min(p.data[ix,:,:,iv])
63 self.
vlim[1]=np.max(p.data[ix,:,:,iv])
65 self.
vlim[0]=np.min(p.data[:,iy,:,iv])
66 self.
vlim[1]=np.max(p.data[:,iy,:,iv])
68 self.
vlim[0]=np.min(p.data[:,:,iz,iv])
69 self.
vlim[1]=np.max(p.data[:,:,iz,iv])
71 xr=p.active_corners[:,0]
72 self.
xlim[0]=np.minimum(self.
xlim[0],xr[0])
73 self.
xlim[1]=np.maximum(self.
xlim[1],xr[1])
74 yr=p.active_corners[:,1]
75 self.
ylim[0]=np.minimum(self.
ylim[0],yr[0])
76 self.
ylim[1]=np.maximum(self.
ylim[1],yr[1])
77 zr=p.active_corners[:,2]
78 self.
zlim[0]=np.minimum(self.
zlim[0],zr[0])
79 self.
zlim[1]=np.maximum(self.
zlim[1],zr[1])
81 self.
vlim[0]=np.minimum(self.
vlim[0],np.min(p.data[ix,:,:,iv]))
82 self.
vlim[1]=np.maximum(self.
vlim[1],np.max(p.data[ix,:,:,iv]))
84 self.
vlim[0]=np.minimum(self.
vlim[0],np.min(p.data[:,iy,:,iv]))
85 self.
vlim[1]=np.maximum(self.
vlim[1],np.max(p.data[:,iy,:,iv]))
87 self.
vlim[0]=np.minimum(self.
vlim[0],np.min(p.data[:,:,iz,iv]))
88 self.
vlim[1]=np.maximum(self.
vlim[1],np.max(p.data[:,:,iz,iv]))
90 print(
'xlim:',self.
xlim)
91 print(
'ylim:',self.
ylim)
92 print(
'zlim:',self.
zlim)
93 print(
'vlim:',self.
vlim)
103 def patch_files(iout=-1, run='', data='../data', verbose=0):
105 Return a list of patch file root names, and count patches and snapshots, 106 unless verbose is negative 108 if verbose>0: print(
'patch_files: run={} data={} iout={}'.format(run,data,iout))
109 rundir=data+
'/'+run+
'/' 111 dir=
'{:05d}'.format(iout)
112 if os.path.isdir(rundir+dir):
117 dirs=[d
for d
in np.sort(os.listdir(rundir))
if os.path.isdir(rundir+d)]
118 dirs=[d
for d
in dirs
if re.match(
"^[0-9][0-9][0-9][0-9][0-9]$",d)]
122 if verbose>1: print(rundir+dir)
123 ff=np.sort(os.listdir(rundir+dir))
125 ff=[rundir+dir+
'/'+f.replace(
'.txt',
'')
for f
in ff
if f.endswith(ext)]
130 if iout>=0
and np.size(files)>0:
131 if verbose>0: print(npatch,
'patches found')
136 print(nsnaps,
'snapshots with',npatch,
'patches each found')
138 print(npatch,
'patches found')
141 def patch (id=2, iout=0, run='', data='../data', overlap=0.01, verbose=0, **kwargs):
143 Return patch structure for a single patch, identified by id and iout 145 name=data+
'/'+run+
'/{:05d}/{:05d}'.format(iout,id)
146 if os.path.isfile(name+
'.txt'):
147 return dd.Patch(name,overlap=overlap,verbose=verbose)
149 files = patch_files(iout,run,data)
152 with open(file+
'.txt',
'r') as fd: 153 success = dd.Search (fd, id, verbose) 155 p=dd.Patch(file,fd=fd,overlap=overlap,verbose=verbose)
160 def is_inside(point,p,verbose=0):
161 ''' Return True/False depending on if the point is inside the patch 163 left=point>=p.active_corners[0]
164 right=point<=p.active_corners[1]
165 left1=p.centre[0]-p.size[0]*0.5
166 left2=p.active_corners[0][0]
168 print(
'{:03d} {} {:8.1f} {:8.1f}'.format(p.id,p.centre,left1,left2))
169 return left.all()
and right.all()
171 def count_inside(point,p,verbose=0):
172 ''' Return the number of dimension in which point is inside the patch 174 left=point>=p.active_corners[0]
175 right=point<=p.active_corners[1]
176 left1=p.centre[0]-p.size[0]*0.5
177 left2=p.active_corners[0][0]
179 print(
'{:03d} {} {:8.1f} {:8.1f}'.format(p.id,p.centre,left1,left2))
182 if left[i]
and right[i]:
186 def patch_at(point,pp,verbose=0):
187 ''' Find the patch at a given point 192 if is_inside(point,p,verbose):
198 def patches_along(point,pp,dir=0,verbose=0):
199 ''' Get the patches along a given direction through a point 203 rfmt =
lambda x:
"%8.2f" % x
204 dict = {
'float':rfmt}
205 np.set_printoptions(formatter=dict)
211 print(
'A {:3d} {}{}'.format(p.id,p.centre,p.size))
212 pt[dir]=p.centre[dir]-p.size[dir]*0.6
221 print(
'B {:3d}{}{}'.format(p.id,p.centre,p.size))
222 pt[dir]=p.centre[dir]+p.size[dir]*0.6
226 print(
'id: {:4d} position: {}'.format(p.id,p.centre))
229 np.set_printoptions(formatter=
None)
232 def patches_along_x(point,pp,verbose=0):
233 ''' Get the patches along the x direction through a point ''' 234 patches_along(point,pp,dir=0,verbose=verbose)
236 def patches_along_y(point,pp,verbose=0):
237 ''' Get the patches along the y direction through a point ''' 238 patches_along(point,pp,dir=1,verbose=verbose)
240 def along_z(point,pp,verbose=0):
241 ''' Get the patches along the z direction through a point ''' 242 patches_along(point,pp,dir=2,verbose=verbose)
244 def indices_and_weights(point,p,iv=0):
246 Return indices and interpolation weights for a point in a patch p 251 if 'p1' in p.varidx.keys():
252 if iv==p.varidx[
'p1']: x0=p.xs[0]
253 if iv==p.varidx[
'p2']: y0=p.ys[0]
254 if iv==p.varidx[
'p3']: z0=p.zs[0]
255 if 'b1' in p.varidx.keys():
256 if iv==p.varidx[
'b1']: x0=p.xs[0]
257 if iv==p.varidx[
'b2']: y0=p.ys[0]
258 if iv==p.varidx[
'b3']: z0=p.zs[0]
259 corner=np.array([x0,y0,z0])
260 p=(point-corner)/p.ds
261 indices=np.array(p,dtype=np.int32)
263 return indices, weights
265 def values_along(point,pp,dir=0,iv=0,var=None,verbose=0,all=0):
267 Return s,f(s) with s the coordinates and f the values in the iv 268 slot of data, taken along the direction v -- so far restricted 271 patches=patches_along(point,pp,dir=dir,verbose=verbose)
275 ii,w=indices_and_weights(point,p,iv)
277 m=np.zeros(3,dtype=np.int32)
282 if p.ioformat==2
or p.ioformat==4:
288 data=p.data[:,:,:,iv]
293 print(
'id: {:3d} iy: {:5.2f} iz: {:5.2f}'.format(p.id,j+w[1],k+w[2]))
294 for i
in range(m[0],m[0]+n[0]):
295 f1=data[i,j,k ]*(1-w[1])+data[i,j+1,k ]*w[1]
296 f2=data[i,j,k+1]*(1-w[1])+data[i,j+1,k+1]*w[1]
297 f=f1*(1-w[2])+f2*w[2]
299 if iv==p.varidx[
'p1']:
309 for j
in range(m[1],m[1]+n[1]):
310 f1=data[i,j,k ]*(1-w[0])+data[i+1,j,k ]*w[0]
311 f2=data[i,j,k+1]*(1-w[0])+data[i+1,j,k+1]*w[0]
312 f=f1*(1-w[2])+f2*w[2]
314 if iv==p.varidx[
'p2']:
324 for k
in range(m[2],m[2]+n[2]):
325 f1=data[i,j ,k]*(1-w[0])+data[i+1,j ,k]*w[0]
326 f2=data[i,j+1,k]*(1-w[0])+data[i+1,j+1,k]*w[0]
327 f=f1*(1-w[1])+f2*w[1]
329 if iv==p.varidx[
'p3']:
340 def values_in(point,p,dir=0,iv=0,verbose=0,all=0):
342 Return s,f(s) with s the coordinates and f the values in the iv 343 slot of data, taken along the direction v -- so far restricted 348 ii,w=indices_and_weights(point,p,iv)
350 print(
'{:3d} {} {} {:5.1f} {:4.1f} {:4.1f}'.\
351 format(p.id,iv,ii,w[0],w[1],w[2]))
353 m=np.zeros(3,dtype=np.int32)
361 for i
in range(m[0],m[0]+n[0]):
362 f1=p.data[i,j,k ,iv]*(1-w[1])+p.data[i,j+1,k, iv]*w[1]
363 f2=p.data[i,j,k+1,iv]*(1-w[1])+p.data[i,j+1,k+1,iv]*w[1]
364 f=f1*(1-w[2])+f2*w[2]
366 if iv==p.varidx[
'p1']:
376 for j
in range(m[1],m[1]+n[1]):
377 f1=p.data[i,j,k ,iv]*(1-w[0])+p.data[i+1,j,k ,iv]*w[0]
378 f2=p.data[i,j,k+1,iv]*(1-w[0])+p.data[i+1,j,k+1,iv]*w[0]
379 f=f1*(1-w[2])+f2*w[2]
381 if iv==p.varidx[
'p2']:
391 for k
in range(m[2],m[2]+n[2]):
392 f1=p.data[i,j ,k,iv]*(1-w[0])+p.data[i+1,j ,k,iv]*w[0]
393 f2=p.data[i,j+1,k,iv]*(1-w[0])+p.data[i+1,j+1,k,iv]*w[0]
394 f=f1*(1-w[1])+f2*w[1]
396 if iv==p.varidx[
'p3']:
407 def values_along_x(point,p,iv=0,verbose=0):
408 ''' Get the iv values along the x direction through a point ''' 409 s,f=values_along(point,p,dir=0,iv=iv,verbose=verbose)
412 def values_along_y(point,p,iv=0,verbose=0):
413 ''' Get the iv values along the y direction through a point ''' 414 s,f=values_along(point,p,dir=1,iv=iv,verbose=verbose)
417 def values_along_z(point,p,iv=0,verbose=0):
418 ''' Get the iv values along the z direction through a point ''' 419 s,f=values_along(point,p,dir=2,iv=iv,verbose=verbose)
422 def patches (iout=-1, run='', data='../data', verbose=0, x=None, y=None, z=None, same=True, overlap=0.01):
424 Return a list of structures s={id,iout,time,data,...}, where id,iout 425 are the integer patch and snapshot numbers, time is in code units, 426 data is a memory mapping of the patch data, and ... are a number of 427 other patch properties returned by ThisPatch(). 430 time_select = [p for p in pp if p.iout==iout] 431 patch_select = [p for p in pp if p.id==id] 433 if verbose: print(
'patches: run={} data={} iout={}'.format(run,data,iout))
435 snap = snapshot(iout,run,data,overlap=overlap,verbose=verbose)
439 for s
in snap.patches:
444 if s.nghost[0]==nghost:
449 pp=[p
for p
in pp
if x>=p.active_corners[0,0]
and x<p.active_corners[1,0]]
451 pp=[p
for p
in pp
if y>=p.active_corners[0,1]
and y<p.active_corners[1,1]]
453 pp=[p
for p
in pp
if z>=p.active_corners[0,2]
and z<p.active_corners[1,2]]
456 def corners(pp,active=1):
457 '''Corners of patch collection''' 459 c=np.array(pp[0].active_corners)
461 c=np.array(pp[0].corners)
464 co=np.array(p.active_corners)
466 co=np.array(p.corners)
467 for i
in range(3): c[0,i]=min(c[0,i],co[0,i])
468 for i
in range(3): c[1,i]=max(c[1,i],co[1,i])
472 '''Rotate data and extent to landscape format''' 474 e=[e[2],e[3],e[0],e[1]]
479 Cache actual data into the list of structures -- typically for a selected subset 486 ''' return dicts with shell values ''' 495 self.
var[key]=np.array([],dtype=np.float32)
496 for key
in (
'r','dv','x','y','z'):
497 self.
var[key]=np.array([],dtype=np.float32)
500 print (
'vars:',np.sort(self.
var.keys()))
505 if c2.min()<r
and c2.max()>r:
513 rr=np.meshgrid(p.x,p.y,p.z)
517 r2=sum(np.array(rr)**2,0)
518 p.var[
'r']=np.sqrt(r2) 519 p.var['dv']=np.ones(np.shape(r2))*np.product(p.ds,dtype=np.float32)
520 cmp=(r2-r2l)*(r2-r2u)
522 for key
in self.
var.keys():
523 var=p.var[key][w[0],w[1],w[2]]
524 self.
var[key]=np.append(self.
var[key],var)
526 print (npatch,
'patches used',time()-start,
'sec')
529 ''' compute radial components of mass flux and velocity ''' 530 self.
var[
'pr']=(self.
var[
'p1']*self.
var[
'x']+self.
var[
'p2']*self.
var[
'y']+self.
var[
'p3']*self.
var[
'z'])/self.
var[
'r'] 531 self.var['ur']=(self.var['u1']*self.var['x']+self.var['u2']*self.var['y']+self.var['u3']*self.var['z'])/self.var['
r'] 534 ''' compute latitude and longitude ''' 535 self.
var[
'mu']=self.
var[
'z']/self.
var[
'r'] 536 self.var['lat']=np.arcsin(self.
var[
'mu'])
537 self.
var[
'lon']=np.arctan2(self.
var[
'x'],self.
var[
'y'])
541 ''' return a dict with shell data (slower) ''' 552 r2l=(radius-nds*p.ds[0])**2
553 r2u=(radius+nds*p.ds[0])**2
555 if c2.min()<radius
and c2.max()>radius:
563 r2=sum(np.array(np.meshgrid(p.x,p.y,p.z))**2,0)
565 cmp=(r2-r2l)*(r2-r2u)
567 for key
in self.
var.keys():
568 vv=p.var[key][w[0],w[1],w[2]]
570 self.
var[key].append(v)
572 print (npatch,
'patches used',time()-start,
'sec')
574 def mass_fluxes(pp,radius=100,nds=2):
575 ''' return the mass fluxes at given radius ''' 584 f.append(pt.dv*(pt.x*pt.data[ix]++pt.y*pt.data[iy]+pt.z*pt.data[iz]))
586 return np.array(f),np.array(dv)
588 def mass_flux(pp,radius=100,nds=2):
589 ''' return the net mass flux at given radius ''' 590 f,dv=mass_fluxes(pp,radius=radius,nds=nds)
591 return np.sum(f)/np.sum(dv)
593 def find_nan (id=0,it=0,run='',data='../data'):
594 ''' check for NaN ''' 595 print (
"checking snapshot",it,
"of "+data+run+
" for NaN")
600 print (
'checking',it,)
601 p=patch(id,it,run,data)
602 bad=any(np.isnan(p.data[:,:,:,0]))
606 pp=patches(it,run,data)
609 if any(np.isnan(p.data[:,:,:,0])):
617 print(
"first good:",it)
619 print (
"checking",it)
620 pp=patches(it,run,data)
623 if np.isnan(p.data[:,:,:,0]).any():
625 print (
"snapshot",it,
"has",n,
"bad patches")
627 def fmax_rt(p,kappa):
629 Return a data block, containing the f_max(RT) values. input is a patch 630 and kappa in cgs units. 632 p.fmax_rt = du.fmax_rt(p,kappa) 634 from scaling
import scaling, cgs
637 fmax = (p.pg/p.d)**4/p.pg*(k*p.d*p.dx.min())
638 fmax = fmax/(1.+(k*p.d*p.dx.min())**2)
641 def umax_rt(p, courant, cdtd):
643 Return the umax value. Input is a patch, the courant number and cdtd 645 p.umax_rt = du.umax_rt(p) 647 from scaling
import scaling, cgs
649 u_max = sc.stefan*np.pi*16./3.*p.fmax_rt.max()*courant/cdtd
653 ''' return the 3-D index of the maximum ''' 655 return np.unravel_index(i,np.shape(a))
658 ''' return the 3-D index of the maximum ''' 660 return np.unravel_index(i,np.shape(a))
662 def first_and_last(run,data,verbose=0):
665 for it
in range(100,9999):
666 file=data+
'/'+run+
'/{:05d}/00002.txt'.format(it)
672 if not found
and exists:
673 p=patch(id,it,run,data)
675 print (
'found first snapshot',it,
'at t=',p.time)
679 elif found
and not exists:
680 p=patch(id,it-1,run,data)
682 print (
'found last snapshot',it-1,
'at t=',p.time)
689 time=(t1-t0)/max(it1-1-it0,1)
def __init__(self, pp, r=10, dr=2.0, verbose=0)
def __init__(self, pp, radius=10, nds=2, verbose=0)
def radial_components(self)