Commits

Matthew Turk committed ed562fc

Updating read_art_octre.py to work with the data I have on hand

  • Participants
  • Parent commits d29f79a

Comments (0)

Files changed (1)

test_art_octre.py

 from yt.mods import *
-from yt.geometry.oct_container import RAMSESOctreeContainer
-import yt.geometry.selection_routines as sel
 from art_header import read_header
 from art_header import _count_art_octs
 from art_header import _read_art_level_info
 
-parameter_filename = '/u/cmoody3/data/art_snapshots/SFG1/10MpcBox_csf512_a0.330.d'
+parameter_filename = 'ART/10MpcBox_csf512_a0.150.d'
 root_grid_mask_offset,root_grid_offset,child_grid_offset,min_level,max_level = read_header(parameter_filename)
 f = open(parameter_filename,'rb')
 nhydrovars, inoll, level_oct_offsets, level_child_offsets =_count_art_octs(f, child_grid_offset,min_level, max_level)
 
-NX = na.ones(3)*128
+NX = na.ones(3)*256
 LE = na.array([0.0, 0.0, 0.0], dtype='float64')
 RE = na.array([1.0, 1.0, 1.0], dtype='float64')
 
-nocts = 0
-for level in range(max_level):
+nocts = 256**3
+pos = f.tell()
+for level in range(1, max_level+1):
     # read level info into POS
     # 1 is the CPU, -1 is not used, 1 is the domain ID
-    le,fln= _read_art_level_info(f, level_oct_offsets,level,root_level=15)
+    le, fln, nLevel = _read_art_level_info(f, level_oct_offsets,level,root_level=15)
     nocts += le.shape[0]
+    print level, nLevel
+print "NUMBER OF OCTS", nocts
 
 art_oct = RAMSESOctreeContainer(NX, LE, RE)
-art_oct.allocate_domains(na.ones(3)*nocts)
+art_oct.allocate_domains([nocts])
 
-for level in range(max_level):
+root_dx = (RE - LE) / NX
+LL = LE + root_dx/2.0
+RL = RE - root_dx/2.0
+rpos = na.mgrid[ LL[0]:RL[0]:NX[0]*1j,
+                 LL[1]:RL[1]:NX[1]*1j,
+                 LL[2]:RL[2]:NX[2]*1j ]
+rpos = na.vstack([p.ravel() for p in rpos]).T
+print rpos.shape
+art_oct.add(1, 0, -1, rpos, 1)
+
+f.seek(pos)
+for level in range(1, max_level):
     # read level info into POSx
     # 1 is the CPU, -1 is not used, 1 is the domain ID
-    le,fln= _read_art_level_info(f, level_oct_offsets,level,root_level=15)
+    le, fln, nLevel = _read_art_level_info(f, level_oct_offsets,level,root_level=15)
     # int curdom, int curlevel, int ng,
     #       np.ndarray[np.float64_t, ndim=2] pos,
     #        int local_domain):
-    art_oct.add(1, level, -1, le, 1)
+    fle = (le.astype("float64") / (NX*2**level)) * (RE - LE) + LE
+    print "LEVEL", level, fle.min(), fle.max()
+    art_oct.add(1, level, -1, fle, 1)
 
 # Now let's test a selector
 class Empty(object):