Commits

Matthew Turk  committed c326360

Refactored ARTIOOctreeContainer's add function

This changes add(), which was previously designed to be called from Python with
arrays, to add_oct, a cdef'd function that adds individual Octs. Doug has
mentioned he has an idea for refactoring this more extensively.

  • Participants
  • Parent commits 6764800

Comments (0)

Files changed (3)

File yt/frontends/artio/_artio_caller.pyx

         oct_count = 0
         level = 0
         for iz in range(self.num_grid/2) :
-            pos[0,2] = iz*2+1
+            dpos[2] = iz*2+1
             for iy in range(self.num_grid/2) :
-                pos[0,1] = iy*2+1
+                dpos[1] = iy*2+1
                 for ix in range(self.num_grid/2) :
-                    pos[0,0]=ix*2+1
-                    oct_handler.add(self.cpu+1, level, self.num_grid, pos, self.domain_id) 
+                    dpos[0]=ix*2+1
+                    oct_handler.add_oct(self.cpu+1, level, dpos)
                     oct_count += 1
 
         # Now do real ART octs
                     pos[0,0] = dpos[0]
                     pos[0,1] = dpos[1]
                     pos[0,2] = dpos[2]
-                    oct_count += oct_handler.add(self.cpu+1, level+1, self.num_grid, pos, self.domain_id )
+                    oct_count += oct_handler.add_oct(self.cpu+1, level+1, dpos)
 
                 status = artio_grid_read_level_end( self.handle )
                 check_artio_status(status) 

File yt/geometry/oct_container.pxd

 
 cdef class ARTIOOctreeContainer(OctreeContainer):
     cdef OctAllocationContainer **domains
+    cdef int add_oct(self, int curdom, int curlevel,
+                     np.float64_t pp[3])
 
 cdef class RAMSESOctreeContainer(OctreeContainer):
     cdef OctAllocationContainer **domains

File yt/geometry/oct_container.pyx

     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
-    def add(self, int curdom, int curlevel, int ng,
-            np.ndarray[np.float64_t, ndim=2] pos,
-            int local_domain, int skip_boundary = 1):
+    cdef int add_oct(self, int curdom, int curlevel, 
+                     np.float64_t pp[3]):
         cdef int level, no, p, i, j, k, ind[3]
-        cdef int local = (local_domain == curdom)
         cdef Oct *cur, *next = NULL
-        cdef np.float64_t pp[3], cp[3], dds[3]
-        no = pos.shape[0] #number of octs
-        if curdom > self.max_domain: curdom = local_domain
-        if curdom < 1 : curdom = 1
+        cdef np.float64_t cp[3], dds[3]
         cdef OctAllocationContainer *cont = self.domains[curdom - 1]
         cdef int initial = cont.n_assigned
-        cdef int in_boundary = 0
-        for p in range(no): #oct number p out of no octs to allocate
-            #for every oct we're trying to add, find the 
-            #floating point unitary position on this level
-            in_boundary = 0
+        for i in range(3):
+            dds[i] = (self.DRE[i] + self.DLE[i])/self.nn[i]
  1. Douglas Rudd

    Sorry for the flood, but need to be careful in this code. Shouldn't this be DRE - DLE? Is the bug hidden because DLE is usually 0?

    I'm replacing the calculation of the true oct position with the following:

            for i in range(3):
                if pp[i] < self.DLE[i] or pp[i] > self.DRE[i] :
                    raise RuntimeError
                dds = (self.DRE[i] - self.DLE[i])/(<np.int64_t)self.nn[i]<<curlevel)                                        
                pos[i] = <np.int64_t> floor((pp[i]-self.DLE[i])/dds)
    

    any comments?

    1. Matthew Turk author

      It absolutely should be. This, like the + self.DLE[i] fix, was something I caught in RAMSES, but I didn't have the ARTIO stuff merged in at that time and so I missed this.

+            ind[i] = <np.int64_t> floor((pp[i]-self.DLE[i])/dds[i])
    1. Matthew Turk author

      In this case, no -- in fact, ind here is always going to be between 0 and self.nn[i]. ind also will be the child index once we traversed the root level, so it will be either 0 or 1 for level > 0.

          1. Douglas Rudd

            I like it. Kinda.

            The true shift needs to take into account the difference in levels, so

            for i in range(3) :
              if pos[i] < (cur.pos[i] << (curlevel-level)) :
                ind[i] = 0
              else :
                ind[i] = 1
            
+            cp[i] = (ind[i] + 0.5) * dds[i]
+        cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
+        if cur == NULL:
+            if curlevel != 0:
+                raise RuntimeError
+            cur = &cont.my_octs[cont.n_assigned] 
    1. Matthew Turk author

      Yes, I think you may be right. My connection is increasingly high latency for a little while so I may not be able to make this change immediately.

+            cur.parent = NULL
+            cur.level = 0
             for i in range(3):
-                pp[i] = pos[p, i] 
-                dds[i] = (self.DRE[i] + self.DLE[i])/self.nn[i]
-                ind[i] = <np.int64_t> floor((pp[i]-self.DLE[i])/dds[i])
-                cp[i] = (ind[i] + 0.5) * dds[i]
-                if ind[i] < 0 or ind[i] >= self.nn[i]:
-                    in_boundary = 1
-            if skip_boundary == in_boundary == 1: continue
-            cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
-            if cur == NULL:
-                if curlevel != 0:
-                    raise RuntimeError
-                cur = &cont.my_octs[cont.n_assigned] 
-                cur.parent = NULL
-                cur.level = 0
+                cur.pos[i] = ind[i]
+            cont.n_assigned += 1
+            self.nocts += 1
+            self.root_mesh[ind[0]][ind[1]][ind[2]] = cur
+        # Now we find the location we want
+        for level in range(curlevel):
+            # At every level, find the cell this oct lives inside
+            for i in range(3):
+                #as we get deeper, oct size halves
+                dds[i] = dds[i] / 2.0
+                if cp[i] > pp[i]: 
+                    ind[i] = 0
+                    cp[i] -= dds[i]/2.0
+                else:
+                    ind[i] = 1
+                    cp[i] += dds[i]/2.0
+            # Check if it has not been allocated
+            next = cur.children[ind[0]][ind[1]][ind[2]]
+            if next == NULL:
+                next = &cont.my_octs[cont.n_assigned]
+                cur.children[ind[0]][ind[1]][ind[2]] = next
+                cont.n_assigned += 1
+                next.parent = cur
                 for i in range(3):
-                    cur.pos[i] = ind[i]
-                cont.n_assigned += 1
+                    next.pos[i] = ind[i] + (cur.pos[i] << 1)
+                next.level = level + 1  
                 self.nocts += 1
-                self.root_mesh[ind[0]][ind[1]][ind[2]] = cur
-            # Now we find the location we want
-            for level in range(curlevel):
-                # At every level, find the cell this oct lives inside
-                for i in range(3):
-                    #as we get deeper, oct size halves
-                    dds[i] = dds[i] / 2.0
-                    if cp[i] > pp[i]: 
-                        ind[i] = 0
-                        cp[i] -= dds[i]/2.0
-                    else:
-                        ind[i] = 1
-                        cp[i] += dds[i]/2.0
-                # Check if it has not been allocated
-                next = cur.children[ind[0]][ind[1]][ind[2]]
-                if next == NULL:
-                    next = &cont.my_octs[cont.n_assigned]
-                    cur.children[ind[0]][ind[1]][ind[2]] = next
-                    cont.n_assigned += 1
-                    next.parent = cur
-                    for i in range(3):
-                        next.pos[i] = ind[i] + (cur.pos[i] << 1)
-                    next.level = level + 1  
-                    self.nocts += 1
-                cur = next
-            cur.domain = curdom 
-            if local == 1:
-                cur.ind = p
-            cur.level = curlevel
+            cur = next
+        cur.domain = curdom 
+        cur.ind = 1
+        cur.level = curlevel
         return cont.n_assigned - initial
     # ii:mask/art ; ci=ramses loop backward (k<-fast, j ,i<-slow) 
     # ii=0 000 art 000 ci 000