Changeset 39

Timestamp:

10/19/10 02:43:42 (19 months ago)

Author:

tabel

Message:

Speed up flash reading. Debugged enclosed gas mass and masses in flash support. Optimized Aake Nordlunds powerspectrum routine (factor of 4 faster now).

Location:

trunk/TOOLS

Files:

• AMRslice.pro (modified) (4 diffs)
• find_and_set_parent_grids.pro (modified) (3 diffs)
• flash2/read_all_flash2_grid_hdf5.pro (modified) (4 diffs)
• get_data.pro (modified) (5 diffs)
• get_units.pro (modified) (1 diff)
• grid_in_cube.pro (modified) (2 diffs)
• powerspectrum.pro (modified) (3 diffs)
• print_for_xmgr.pro (modified) (1 diff)

trunk/TOOLS/AMRslice.pro

@@ -616 +616 @@
            step = 1
            exists = ''
-           while ((step lt 100) and (exists[0] eq '')) do begin
+           while ((step lt 1000) and (exists[0] eq '')) do begin
 ;              findS = sp[0:ns-2] + STRING(max([LONG(timenum)-step,0]), FORMAT='(i4.4)')
               if regular then $
-                 findS = sp[0:ns-2] + STRING(max([LONG(timenum)-step,0]), FORMAT='(i4.4)') $
+                 findS = sp[0:((ns-2) > 0)] + STRING(max([LONG(timenum)-step,0]), FORMAT='(i4.4)') $
               else $
                  findS = sp[0:ns-1] + STRING(max([LONG(timenum)-step,0]), FORMAT='(i3.3)') 
@@ -674 +674 @@
            step = 1
            exists = ''
-           while ((step lt 100) and (exists[0] eq '')) do begin
+           while ((step lt 1000) and (exists[0] eq '')) do begin
               if regular then $
-                 findS = sp[0:ns-2] + STRING(LONG(timenum)+step, FORMAT='(i4.4)') $
+                 findS = sp[0:((ns-2) > 0)] + STRING(LONG(timenum)+step, FORMAT='(i4.4)') $
               else $
                  findS = sp[0:ns-1] + STRING(LONG(timenum)+step, FORMAT='(i3.3)') 
@@ -1206 +1206 @@
        construct_cube, datacube, fieldname
        print, '<data^2> =', mean(datacube^2)
-;       if vel.myturn eq 0 then ps = PowerSpectrum_3D(datacube,
-;       obin=obin) else $
+;       if vel.myturn eq 0 then ps = PowerSpectrum_3D(datacube, obin=obin)  $
        if vel.myturn eq 0 then $
-          power3d, reform(datacube[*,*,*]), $
-                   wavenumbers=obin, spectrum=ps,average=1,debug=verbose $
+          power3d, reform(datacube[*,*,*]), wavenumbers=obin, spectrum=ps,average=1,debug=verbose $
        else begin
 ;          ps1 = PowerSpectrum_3D(reform(datacube[*,*,*,0]), obin=obin)
@@ -2083 +2081 @@
   wItem = WIDGET_BUTTON(wMenu, VALUE='Color Map', UVALUE='CM_BUTTON', ACCELERATOR='Ctrl+C')
   wItem= WIDGET_BUTTON(wMenu, VALUE='Change Min and Max', UVALUE='SET_BUTTON', ACCELERATOR='Ctrl+Shift+M')
-  wItem = WIDGET_BUTTON(wMenu, VALUE='Data info', UVALUE='INFO_BUTTON', ACCELERATOR='Ctrl+Shift+I')
+  wItem = WIDGET_BUTTON(wMenu, VALUE='Dataset info', UVALUE='INFO_BUTTON', ACCELERATOR='Ctrl+Shift+I')
   wItem = WIDGET_BUTTON(wMenu, VALUE='Hierarchy file', UVALUE='VIEW_HIERARCHY')
   wItem = WIDGET_BUTTON(wMenu, $

trunk/TOOLS/find_and_set_parent_grids.pro

@@ -6 +6 @@
   num_of_levels = N_elements(count) 
 ;breaki
+  sind = sort(grid_info.num)
   for i=1, num_of_levels-1 DO BEGIN
      print, 'level: ', i
@@ -17 +18 @@
 ; set parent start and end indeces 
         IF (grid_info[cg].parent GT 0) THEN BEGIN
-           pi = where(grid_info[b].num eq grid_info[cg].parent, tcount) ; all parents
+;           pi = where(grid_info[b].num eq grid_info[cg].parent, tcount) ; search in all possible parents
+           pi = sind[grid_info[cg].parent-1]
+           tcount = 1
+;           print, grid_info[b[pi[0]]].num, grid_info[cg].parent, pi[0], b[pi[0]], 
  ;          print, tcount
            if tcount gt 0 then begin
-              pi  = b[pi]
+;              pi  = b[pi]
               p = grid_info[pi]
               c = grid_info[cg]
@@ -37 +41 @@
 
 function old_set_parent_start_end_indices, grid_info
-  print, 'start set_parent_start_end_indices'
+  print, 'start old_set_parent_start_end_indices'
   b =  0
   count = histogram(grid_info.level,reverse_indices=r)

trunk/TOOLS/flash2/read_all_flash2_grid_hdf5.pro

@@ -93 +93 @@
            dd = double(cg.right_edge-cg.left_edge)/ra
            d = ra
+           dims = N_elements(d)
            this_stride= (stride_defined) ? reform(stride[cnt,*]) : (intarr(n_elements(d)) + 1)
            this_start = (start_defined) ? reform(start[cnt,*]) : intarr(n_elements(d))
@@ -138 +139 @@
                ratio = this_count/this_stride > 1
 ;                   if total(this_start lt d) ne 3 then continue
-               dims = N_elements(d)
+               dims = N_elements(d)-1
 ;               print, dims
                if n_elements(d) ne 1 then begin
-                  H5S_SELECT_HYPERSLAB, dataspace_id, [this_start[0:dims-2],cg.num-1], [ratio[0:dims-2],1], STRIDE=[this_stride[0:dims-2],1], /RESET  
+                  H5S_SELECT_HYPERSLAB, dataspace_id, [this_start[0:dims-1],cg.num-1], [ratio[0:dims-1],1], STRIDE=[this_stride[0:dims-1],1], /RESET  
                   memory_space_id = H5S_CREATE_SIMPLE([ratio,1])  
                   data = H5D_READ(dsetid, FILE_SPACE=dataspace_id, MEMORY_SPACE=memory_space_id)  
@@ -154 +155 @@
 ;breaki
             ENDELSE ; not special
-            
+            d = d[0:dims-1]
             all[cnt].file_name   = gident
             all[cnt].np_total    = N_elements(data)
             s = size(data)
             if (s[0] eq 2) then s[3] = 1
+            
             all[cnt].li[*]       = this_start[*]
             all[cnt].ri[*]       = (s[1:3]+this_start-1)[*]
@@ -166 +168 @@
             if fields[i] eq 'particle massdensity' then data *= product(dd)
 ;            print,  product(dd[where(d gt 0.)]),  'dd:', dd[where(d gt 0.)], 'where ',where(d gt 0.)
+;            breaki
             if fields[i] eq 'mass' then data *= product(dd[where(d gt 0.)])
             all[cnt].Data[i]     = ptr_new(data, /no_copy)
+     if check_cancel() then goto, enough
           ENDFOR                 ; loop over patches
 
      ENDFOR                     ; loop over fields
-
+enough:
      if notspecial then H5D_CLOSE, dsetid           
      H5F_CLOSE, fid
-     if check_cancel() then goto, enough
+
   ENDFOR                        ; loop over unique files
 ; return information in the order grids and files were specified
   all = all[revind]
-enough:
+
   RETURN
 END

trunk/TOOLS/get_data.pro

@@ -9 +9 @@
 pro pack_as_list, use_g, a
 ; pack all data into 1D arrays 
-   zero = 0.D
+   zero = -1e30
    zero_solution_under_subgrid, use_g, a, ZERO=zero
 
@@ -153 +153 @@
 end
 
+function fn, aa, name
+@common_blocks.inc
+ind = where(strtrim(aa[0].data_fields,2) eq strtrim(name,2))
+if ind[0] eq -1 then begin
+   pind = where(name eq *parti.names)
+   if pind[0] ne -1 then return, parti.d[pind]
+   print, 'no field with name >'+name+'< found in ', aa
+   return, 0
+endif
+
+Num = N_elements(aa)
+res = (*aa[0].data[ind[0]])
+for i=1UL,Num-1 do res = [res,  (*aa[i].data[ind[0]])]
+return,  res
+end
+
 function enclosed_gas_mass, cen,x,y,z, mass
 x = x-cen[0]
@@ -162 +178 @@
 ;breaki
 return,Menc
-end
-
-function fn, aa, name
-@common_blocks.inc
-ind = where(strtrim(aa[0].data_fields,2) eq strtrim(name,2))
-if ind[0] eq -1 then begin
-   pind = where(name eq *parti.names)
-   if pind[0] ne -1 then return, parti.d[pind]
-   print, 'no field with name >'+name+'< found in ', aa
-   return, 0
-endif
-
-Num = N_elements(aa)
-res = (*aa[0].data[ind[0]])
-for i=1L,Num-1 do res = [res,  (*aa[i].data[ind[0]])]
-return,  res
 end
 
@@ -190 +190 @@
 *a[0].data[i] = Menc
 help, menc
+mi = MIN(Menc, max=ma)
+print, 'min, max:', mi, ma
 ;breaki
 return
@@ -610 +612 @@
 ;breaki
 
-print, 'RRAAAAAAAAAALLL'
-help, ral
+;print, 'RRAAAAAAAAAALLL'
+;help, ral
    for j=0L,N_elements(a)-1 DO BEGIN
       heap_free, a[j].data[Nvar:*] 

trunk/TOOLS/get_units.pro

@@ -163 +163 @@
         unitfactor = densityunit*lengthunit^3/1.989d33
         unitlabel = textoidl('M')+'!D!9n!X!N'
+        print, 'massunit:', unitfactor
      end
      (tlabel eq 'Dynamical time (gas)'): Begin

trunk/TOOLS/grid_in_cube.pro

@@ -18 +18 @@
   delta_distance = (grid_info.Right_edge - grid_info.Left_edge)/(index_range)
 
-  suff_res = delta_distance[0,*]
+;  suff_res = delta_distance[0,*]
   suff_res = where(delta_distance[0,*] ge b)
 ;  print, suff_res
-  print, fix(total(suff_res gt 0)), ' grids have sufficient resolution for cube' 
-  if total(suff_res gt 0) lt 1 THEN BEGIN
+  print, N_elements(suff_res), ' grids have sufficient resolution for cube' 
+  if suff_res[0] lt 0 THEN BEGIN
       grid_in_cube = grid_info[0]
       RETURN, grid_in_cube
@@ -40 +40 @@
           (Left_edge[2,*]   lt max_right[2])          )
 
-  print, 'there are ', total(help),' grids at least partially in cube'
+  print, 'there are ', total(help),' grids of sufficient resolution at least partially in cube'
   if total(help) gt 1 then grid_in = new_i(where(help gt 0))
   return, grid_in

trunk/TOOLS/powerspectrum.pro

@@ -368 +368 @@
 if n_elements(average) ne 1 then average=spectrum
 
+hist = histogram(rad, reverse_indices=r, min=0.+o-0.5, max=imax+o-0.499, nbins=imax, locations=obin)
+;breaki
 if keyword_set(debug) then print,'looping over shells'
 t0=systime(1)
 for i=0L,imax-1 do begin                                ; cut the corners?
-  ii=where((rad ge i+o-.5) and (rad lt i+o+.5),nw)      ; indices into shell
-  if keyword_set(debug) and i le 10 then print,i,nw
+;  ii=where((rad ge i+o-.5) and (rad lt i+o+.5),nw)      ; indices into shell
+   ii = [r[r[i]:r[i+1]-1]]      ;Tom Abel 2010 use this to speed up routine use line above for original
+   nw = N_elements(ii)
+  if keyword_set(debug) and i le 12 then print,i,nw; , N_elements(ij)
   if s(0) eq 3 then begin
     tmp=ta(ii) & tmp=abs(temporary(tmp))^2              ; conserve memory
@@ -380 +384 @@
   end
   ka(i)=aver(rad(ii))                                   ; average wave number
+;  print, ka(i), aver(rad[ij])
   if keyword_set(kindex) then ka(i)=i+plo
   if n_elements(average) ne 1 then begin                ; average array or not?
@@ -401 +406 @@
 
 if n_elements(c) gt 0 then begin
-  power3d,c,aver=aver,compensate=compensate,kindex=kindex,sp=sp1,/noplot,offset=o,potsdam=potsdam
+  power3d,c,aver=aver,compensate=compensate,kindex=kindex,sp=sp1,/noplot,offset=o,potsdam=potsdam,debug=debug
   spectrum=spectrum+sp1
 end
 if n_elements(b) gt 0 then begin
-  power3d,b,aver=aver,compensate=compensate,kindex=kindex,sp=sp1,/noplot,offset=o,potsdam=potsdam
+  power3d,b,aver=aver,compensate=compensate,kindex=kindex,sp=sp1,/noplot,offset=o,potsdam=potsdam,debug=debug
   spectrum=spectrum+sp1
 end

trunk/TOOLS/print_for_xmgr.pro

@@ -1 @@
-pro print_for_xmgr, x,y, set=set, file=filename
+ pro print_for_xmgr, x,y, set=set, file=filename
 ;
 ; print out x and y values for to be read by xmgr
 ;
-
+Nbins = N_elements(x)
+if N_elements(set) eq 0 then set=0
 if N_elements(filename) gt 0 then begin
    openw, lun, filename, /get_lun, /APPEND