DISPATCH
dispatch_utils.py
1 
2 # Pythn 2/3 compatibility
3 from __future__ import print_function
4 import sys
5 if sys.version_info[0]!=2:
6  from importlib import reload
7 
8 import numpy as np
9 import os
10 import re
11 from time import time
12 
13 import dispatch_data as dd
14 reload(dd)
15 
16 def dispatch_utils():
17  '''
18  Main methods: patches, cache
19  patches() Return list of objects for all 'data/run' patches
20  cache() Adds cached arrays to the object
21  Typical use:
22  pp=patches(iout,'run')
23  pp=[p for p in pp if select(p)]
24  pp=cache(pp)
25  for p in pp:
26  ... do something with the data in p.d, p.ux, ...
27  Other methods:
28  patch_files() Return list of patch file root names
29  p=patch() Return one patch object
30  '''
31  print( "Type dispatch_utils?<RETURN> for help text, or dispach_utils(\
32  for pop-up help in Canopy")
33 
34 def snapshot (iout=0, run='', data='../data', overlap=0.01, verbose=0):
35  snap=dd.Snapshot (iout, run, data, overlap=overlap, verbose=verbose)
36  snap.iout=iout
37  snap.run=run
38  snap.data=data
39  snap.dict={}
40  for p in snap.patches:
41  snap.dict[p.id]=p
42  return snap
43 
44 class limits():
45  ''' Return extent limits in x,y,z,values '''
46  def __init__(self,pp,ix=None,iy=None,iz=None,iv=0,verbose=0):
47  self.xlim=np.zeros(2)
48  self.ylim=np.zeros(2)
49  self.zlim=np.zeros(2)
50  self.vlim=np.zeros(2)
51  p=pp[0]
52  xr=p.active_corners[:,0]
53  self.xlim[0]=xr[0]
54  self.xlim[1]=xr[1]
55  yr=p.active_corners[:,1]
56  self.ylim[0]=yr[0]
57  self.ylim[1]=yr[1]
58  zr=p.active_corners[:,2]
59  self.zlim[0]=zr[0]
60  self.zlim[1]=zr[1]
61  if ix:
62  self.vlim[0]=np.min(p.data[ix,:,:,iv])
63  self.vlim[1]=np.max(p.data[ix,:,:,iv])
64  if iy:
65  self.vlim[0]=np.min(p.data[:,iy,:,iv])
66  self.vlim[1]=np.max(p.data[:,iy,:,iv])
67  if iz:
68  self.vlim[0]=np.min(p.data[:,:,iz,iv])
69  self.vlim[1]=np.max(p.data[:,:,iz,iv])
70  for p in pp:
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])
80  if ix:
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]))
83  if iy:
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]))
86  if iz:
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]))
89  if verbose:
90  print('xlim:',self.xlim)
91  print('ylim:',self.ylim)
92  print('zlim:',self.zlim)
93  print('vlim:',self.vlim)
94 
95 def fsplit(name):
96  res=name.split('_')
97  if np.size(res)==2:
98  return res[0],res[1]
99  else:
100  res=name.split('/')
101  return res[1],res[0]
102 
103 def patch_files(iout=-1, run='', data='../data', verbose=0):
104  '''
105  Return a list of patch file root names, and count patches and snapshots,
106  unless verbose is negative
107  '''
108  if verbose>0: print('patch_files: run={} data={} iout={}'.format(run,data,iout))
109  rundir=data+'/'+run+'/'
110  if iout>=0:
111  dir='{:05d}'.format(iout)
112  if os.path.isdir(rundir+dir):
113  dirs=[dir]
114  else:
115  dirs=[]
116  else:
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)]
119  nsnaps=np.size(dirs)
120  files=[]
121  for dir in dirs:
122  if verbose>1: print(rundir+dir)
123  ff=np.sort(os.listdir(rundir+dir))
124  ext='.txt'
125  ff=[rundir+dir+'/'+f.replace('.txt','') for f in ff if f.endswith(ext)]
126  if np.size(ff)>0:
127  npatch=np.size(ff)
128  for f in ff:
129  files.append(f)
130  if iout>=0 and np.size(files)>0:
131  if verbose>0: print(npatch,'patches found')
132  return files
133 
134  if verbose>0:
135  if iout<0:
136  print(nsnaps,'snapshots with',npatch,'patches each found')
137  else:
138  print(npatch,'patches found')
139  return files
140 
141 def patch (id=2, iout=0, run='', data='../data', overlap=0.01, verbose=0, **kwargs):
142  '''
143  Return patch structure for a single patch, identified by id and iout
144  '''
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)
148  else:
149  files = patch_files(iout,run,data)
150  p=None
151  for file in files:
152  with open(file+'.txt','r') as fd:
153  success = dd.Search (fd, id, verbose)
154  if (success):
155  p=dd.Patch(file,fd=fd,overlap=overlap,verbose=verbose)
156  if not p:
157  print ('not found')
158  return p
159 
160 def is_inside(point,p,verbose=0):
161  ''' Return True/False depending on if the point is inside the patch
162  '''
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]
167  if verbose>0:
168  print('{:03d} {} {:8.1f} {:8.1f}'.format(p.id,p.centre,left1,left2))
169  return left.all() and right.all()
170 
171 def count_inside(point,p,verbose=0):
172  ''' Return the number of dimension in which point is inside the patch
173  '''
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]
178  if verbose>0:
179  print('{:03d} {} {:8.1f} {:8.1f}'.format(p.id,p.centre,left1,left2))
180  count=0
181  for i in range(3):
182  if left[i] and right[i]:
183  count+=1
184  return count
185 
186 def patch_at(point,pp,verbose=0):
187  ''' Find the patch at a given point
188  '''
189  level=-1
190  p1=None
191  for p in pp:
192  if is_inside(point,p,verbose):
193  if p.level>level:
194  p1=p
195  level=p.level
196  return p1
197 
198 def patches_along(point,pp,dir=0,verbose=0):
199  ''' Get the patches along a given direction through a point
200  '''
201  pt=np.copy(point)
202  if verbose>0:
203  rfmt = lambda x: "%8.2f" % x
204  dict = {'float':rfmt}
205  np.set_printoptions(formatter=dict)
206  p=patch_at(pt,pp)
207  if p:
208  p1=p
209  while p:
210  if verbose>1:
211  print('A {:3d} {}{}'.format(p.id,p.centre,p.size))
212  pt[dir]=p.centre[dir]-p.size[dir]*0.6
213  p=patch_at(pt,pp)
214  if p:
215  p1=p
216  p=p1
217  out=[]
218  while p:
219  out.append(p)
220  if verbose>1:
221  print('B {:3d}{}{}'.format(p.id,p.centre,p.size))
222  pt[dir]=p.centre[dir]+p.size[dir]*0.6
223  p=patch_at(pt,pp)
224  if verbose>0:
225  for p in out:
226  print('id: {:4d} position: {}'.format(p.id,p.centre))
227  else:
228  out=None
229  np.set_printoptions(formatter=None)
230  return out
231 
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)
235 
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)
239 
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)
243 
244 def indices_and_weights(point,p,iv=0):
245  '''
246  Return indices and interpolation weights for a point in a patch p
247  '''
248  x0=p.x[0]
249  y0=p.y[0]
250  z0=p.z[0]
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)
262  weights=p-indices
263  return indices, weights
264 
265 def values_along(point,pp,dir=0,iv=0,var=None,verbose=0,all=0):
266  '''
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
269  to axis values
270  '''
271  patches=patches_along(point,pp,dir=dir,verbose=verbose)
272  ss=[]
273  ff=[]
274  for p in patches:
275  ii,w=indices_and_weights(point,p,iv)
276  if all:
277  m=np.zeros(3,dtype=np.int32)
278  n=p.gn
279  else:
280  m=p.nghost
281  n=p.n
282  if p.ioformat==2 or p.ioformat==4:
283  n=n-1
284  if var:
285  p.vars()
286  data=p.var[var]
287  else:
288  data=p.data[:,:,:,iv]
289  if dir==0:
290  j=ii[1]
291  k=ii[2]
292  if verbose>2:
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]
298  if 'p1' in p.varidx:
299  if iv==p.varidx['p1']:
300  ss.append(p.xs[i])
301  else:
302  ss.append(p.x[i])
303  else:
304  ss.append(p.x[i])
305  ff.append(f)
306  elif dir==1:
307  i=ii[0]
308  k=ii[2]
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]
313  if 'p2' in p.varidx:
314  if iv==p.varidx['p2']:
315  ss.append(p.ys[j])
316  else:
317  ss.append(p.y[j])
318  else:
319  ss.append(p.y[j])
320  ff.append(f)
321  else:
322  i=ii[0]
323  j=ii[1]
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]
328  if 'p3' in p.varidx:
329  if iv==p.varidx['p3']:
330  ss.append(p.zs[k])
331  else:
332  ss.append(p.z[k])
333  else:
334  ss.append(p.z[k])
335  ff.append(f)
336  ss=np.array(ss)
337  ff=np.array(ff)
338  return ss,ff
339 
340 def values_in(point,p,dir=0,iv=0,verbose=0,all=0):
341  '''
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
344  to axis values
345  '''
346  ss=[]
347  ff=[]
348  ii,w=indices_and_weights(point,p,iv)
349  if verbose:
350  print('{:3d} {} {} {:5.1f} {:4.1f} {:4.1f}'.\
351  format(p.id,iv,ii,w[0],w[1],w[2]))
352  if all:
353  m=np.zeros(3,dtype=np.int32)
354  n=p.gn
355  else:
356  m=p.nghost
357  n=p.n
358  if dir==0:
359  j=ii[1]
360  k=ii[2]
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]
365  if 'p1' in p.varidx:
366  if iv==p.varidx['p1']:
367  ss.append(p.xs[i])
368  else:
369  ss.append(p.x[i])
370  else:
371  ss.append(p.x[i])
372  ff.append(f)
373  elif dir==1:
374  i=ii[0]
375  k=ii[2]
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]
380  if 'p2' in p.varidx:
381  if iv==p.varidx['p2']:
382  ss.append(p.ys[j])
383  else:
384  ss.append(p.y[j])
385  else:
386  ss.append(p.y[j])
387  ff.append(f)
388  else:
389  i=ii[0]
390  j=ii[1]
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]
395  if 'p3' in p.varidx:
396  if iv==p.varidx['p3']:
397  ss.append(p.zs[k])
398  else:
399  ss.append(p.z[k])
400  else:
401  ss.append(p.z[k])
402  ff.append(f)
403  ss=np.array(ss)
404  ff=np.array(ff)
405  return ss,ff
406 
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)
410  return s,f
411 
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)
415  return s,f
416 
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)
420  return s,f
421 
422 def patches (iout=-1, run='', data='../data', verbose=0, x=None, y=None, z=None, same=True, overlap=0.01):
423  '''
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().
428 
429  pp = patches('run')
430  time_select = [p for p in pp if p.iout==iout]
431  patch_select = [p for p in pp if p.id==id]
432  '''
433  if verbose: print('patches: run={} data={} iout={}'.format(run,data,iout))
434  #files = patch_files (run=run, data=data, iout=iout, verbose=verbose)
435  snap = snapshot(iout,run,data,overlap=overlap,verbose=verbose)
436  pp=[]
437  nghost=-1
438  #for f in files:
439  for s in snap.patches:
440  #s=dd.Patch(f,verbose=verbose)
441  if same:
442  if nghost==-1:
443  nghost=s.nghost[0]
444  if s.nghost[0]==nghost:
445  pp.append(s)
446  else:
447  pp.append(s)
448  if x:
449  pp=[p for p in pp if x>=p.active_corners[0,0] and x<p.active_corners[1,0]]
450  if y:
451  pp=[p for p in pp if y>=p.active_corners[0,1] and y<p.active_corners[1,1]]
452  if z:
453  pp=[p for p in pp if z>=p.active_corners[0,2] and z<p.active_corners[1,2]]
454  return pp
455 
456 def corners(pp,active=1):
457  '''Corners of patch collection'''
458  if active:
459  c=np.array(pp[0].active_corners)
460  else:
461  c=np.array(pp[0].corners)
462  for p in pp:
463  if active:
464  co=np.array(p.active_corners)
465  else:
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])
469  return c
470 
471 def rotate(d,e):
472  '''Rotate data and extent to landscape format'''
473  d=d.transpose()
474  e=[e[2],e[3],e[0],e[1]]
475  return d,e
476 
477 def cache(pp):
478  '''
479  Cache actual data into the list of structures -- typically for a selected subset
480  '''
481  for p in pp:
482  p.cache()
483 
485  def __init__(self,pp,r=10,dr=2.0,verbose=0):
486  ''' return dicts with shell values '''
487  if verbose:
488  start=time()
489  self.var={}
490  self.r=r
491  self.dr=dr
492  p=pp[0]
493  p.vars()
494  for key in p.var:
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)
498  npatch=0
499  if verbose:
500  print ('vars:',np.sort(self.var.keys()))
501  for p in pp:
502  r2l=(r-dr)**2
503  r2u=(r+dr)**2
504  c2=p.corner_radii()
505  if c2.min()<r and c2.max()>r:
506  if verbose:
507  if verbose>1:
508  print(p.id)
509  npatch+=1
510  p.cache()
511  p.stripit()
512  p.vars()
513  rr=np.meshgrid(p.x,p.y,p.z)
514  p.var['x']=rr[0]
515  p.var['y']=rr[1]
516  p.var['z']=rr[2]
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)
521  w=np.where(cmp<0.0)
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)
525  if verbose:
526  print (npatch,'patches used',time()-start,'sec')
527 
528  def radial_components(self):
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']
532 
533  def angles(self):
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'])
538 
540  def __init__(self,pp,radius=10,nds=2,verbose=0):
541  ''' return a dict with shell data (slower) '''
542  if verbose:
543  start=time()
544  self.var={}
545  p=pp[0]
546  p.vars()
547  for key in p.var:
548  self.var[key]=[]
549  self.var['r2']=[]
550  npatch=0
551  for p in pp:
552  r2l=(radius-nds*p.ds[0])**2
553  r2u=(radius+nds*p.ds[0])**2
554  c2=p.corner_radii()
555  if c2.min()<radius and c2.max()>radius:
556  if verbose:
557  if verbose>1:
558  print(p.id)
559  npatch+=1
560  p.cache()
561  p.stripit()
562  p.vars()
563  r2=sum(np.array(np.meshgrid(p.x,p.y,p.z))**2,0)
564  p.var['r2']=r2
565  cmp=(r2-r2l)*(r2-r2u)
566  w=np.where(cmp<0.0)
567  for key in self.var.keys():
568  vv=p.var[key][w[0],w[1],w[2]]
569  for v in vv:
570  self.var[key].append(v)
571  if verbose:
572  print (npatch,'patches used',time()-start,'sec')
573 
574 def mass_fluxes(pp,radius=100,nds=2):
575  ''' return the mass fluxes at given radius '''
576  s=shell_values(pp,radius=radius,nds=nds)
577  p=pp[0]
578  ix=p.varidx['p1']
579  iy=p.varidx['p2']
580  iz=p.varidx['p3']
581  f=[]
582  dv=[]
583  for pt in s:
584  f.append(pt.dv*(pt.x*pt.data[ix]++pt.y*pt.data[iy]+pt.z*pt.data[iz]))
585  dv.append(pt.dv)
586  return np.array(f),np.array(dv)
587 
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)
592 
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")
596  bad=True
597  while bad:
598  if id:
599  while bad:
600  print ('checking',it,)
601  p=patch(id,it,run,data)
602  bad=any(np.isnan(p.data[:,:,:,0]))
603  if bad:
604  it=it-1
605  print ("bad")
606  pp=patches(it,run,data)
607  bad=False
608  for p in pp:
609  if any(np.isnan(p.data[:,:,:,0])):
610  bad=True
611  id=p.id
612  break
613  if bad:
614  print("bad")
615  else:
616  print("good")
617  print("first good:",it)
618  it=it+1
619  print ("checking",it)
620  pp=patches(it,run,data)
621  n=0
622  for p in pp:
623  if np.isnan(p.data[:,:,:,0]).any():
624  n+=1
625  print ("snapshot",it,"has",n,"bad patches")
626 
627 def fmax_rt(p,kappa):
628  '''
629  Return a data block, containing the f_max(RT) values. input is a patch
630  and kappa in cgs units.
631 
632  p.fmax_rt = du.fmax_rt(p,kappa)
633  '''
634  from scaling import scaling, cgs
635  sc = scaling(cgs)
636  k=kappa*sc.m/sc.l**2
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)
639  return fmax
640 
641 def umax_rt(p, courant, cdtd):
642  '''
643  Return the umax value. Input is a patch, the courant number and cdtd
644 
645  p.umax_rt = du.umax_rt(p)
646  '''
647  from scaling import scaling, cgs
648  sc = scaling(cgs)
649  u_max = sc.stefan*np.pi*16./3.*p.fmax_rt.max()*courant/cdtd
650  return u_max
651 
652 def minloc(a):
653  ''' return the 3-D index of the maximum '''
654  i = a.argmin()
655  return np.unravel_index(i,np.shape(a))
656 
657 def maxloc(a):
658  ''' return the 3-D index of the maximum '''
659  i = a.argmax()
660  return np.unravel_index(i,np.shape(a))
661 
662 def first_and_last(run,data,verbose=0):
663  id=2
664  found=False
665  for it in range(100,9999):
666  file=data+'/'+run+'/{:05d}/00002.txt'.format(it)
667  try:
668  fd=open(file,'rb')
669  exists=True
670  except:
671  exists=False
672  if not found and exists:
673  p=patch(id,it,run,data)
674  if verbose:
675  print ('found first snapshot',it,'at t=',p.time)
676  found=True
677  it0=it
678  t0=p.time
679  elif found and not exists:
680  p=patch(id,it-1,run,data)
681  if verbose:
682  print ('found last snapshot',it-1,'at t=',p.time)
683  it1=it
684  t1=p.time
685  break
686  class it:
687  first=it0
688  last=it1
689  time=(t1-t0)/max(it1-1-it0,1)
690  return it()
def __init__(self, pp, r=10, dr=2.0, verbose=0)
def __init__(self, pp, radius=10, nds=2, verbose=0)