Commits

Matthew Turk committed fed5981

This is getting closer to a working ART implementation

  • Participants
  • Parent commits 834acd6

Comments (0)

Files changed (2)

File art_header.py

     #le = le/2**(root_level-1-level)-1
 
     #try without the -1
-    le = le/2**(root_level-2-level)-1
+    le = le.astype("float64")
+    le /= 2**(root_level - level + 1)
 
     #now read the hvars and vars arrays
     #we are looking for iOctCh

File test_art_octre.py

 from art_header import read_header
 from art_header import _count_art_octs
 from art_header import _read_art_level_info
+from yt.geometry.oct_container import RAMSESOctreeContainer
+import yt.geometry.selection_routines as sel
 
 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)*256
+NX = na.ones(3)*128 / 2
 LE = na.array([0.0, 0.0, 0.0], dtype='float64')
 RE = na.array([1.0, 1.0, 1.0], dtype='float64')
 
-nocts = 256**3
+nocts = NX.prod()
 pos = f.tell()
+#max_level = 5
 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, nLevel = _read_art_level_info(f, level_oct_offsets,level,root_level=15)
+    le, 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
 rpos = na.vstack([p.ravel() for p in rpos]).T
 print rpos.shape
 art_oct.add(1, 0, -1, rpos, 1)
+assert(art_oct.nocts == rpos.shape[0])
 
 f.seek(pos)
 for level in range(1, max_level):
+    o1 = art_oct.nocts
     # read level info into POSx
     # 1 is the CPU, -1 is not used, 1 is the domain ID
-    le, fln, nLevel = _read_art_level_info(f, level_oct_offsets,level,root_level=15)
+    le, 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):
-    fle = (le.astype("float64") / (NX*2**level)) * (RE - LE) + LE
-    print "LEVEL", level, fle.min(), fle.max()
+    fle = le.astype("float64") / (NX*2**(level - 1))
     art_oct.add(1, level, -1, fle, 1)
+    print "LEVEL", level, fle.min(), fle.max(), art_oct.nocts
+    o2 = art_oct.nocts
+    assert(o2 - o1 == fle.shape[0])
+print "AO", art_oct.nocts, nocts
 
 # Now let's test a selector
 class Empty(object):
 
 sp_obj = Empty()
 sp_obj.center = [0.5, 0.5, 0.5]
-sp_obj.radius = 0.1
+sp_obj.radius = 1.0
 sp = sel.SphereSelector(sp_obj)
 
-mask = na.zeros(sp.nocts, dtype="bool")
-
+mask = sp.select_octs(art_oct)
 cells = art_oct.count_cells(sp, mask)
 nc = cells.sum()
-print nc, sp.nocts
+print nc, art_oct.nocts
 # 1 is domain_id
-icoords = art_oct.icoords(1, mask, nc)
-ires = art_oct.ires(1, mask, nc)
+#icoords = art_oct.icoords(1, mask, nc)
+#ires = art_oct.ires(1, mask, nc)
 
-print icoords
-print ires
+#print icoords
+#print ires