Commits

Matthew Turk  committed b4dc04e

Attempting to deal with boundary issues. This is still not functional.

  • Participants
  • Parent commits 7038c8a
  • Branches geometry_handling

Comments (0)

Files changed (3)

File yt/frontends/ramses/data_structures.py

         self._read_amr_header()
 
     _hydro_offset = None
+    _level_count = None
+
+    @property
+    def level_count(self):
+        if self._level_count is not None: return self._level_count
+        self.hydro_offset
+        return self._level_count
 
     @property
     def hydro_offset(self):
         # It goes: level, CPU, 8-variable
         hydro_offset = na.zeros(self.amr_header['nlevelmax'], dtype='int64')
         hydro_offset -= 1
+        level_count = na.zeros(self.amr_header['nlevelmax'], dtype='int64')
         for level in range(self.amr_header['nlevelmax']):
             for cpu in range(self.amr_header['nboundary'] +
                              self.amr_header['ncpu']):
                            ('file_ncache', 1, 'I') )
                 hvals = fpu.read_attrs(f, header)
                 if hvals['file_ncache'] == 0: continue
-                if cpu + 1 == self.domain_id: hydro_offset[level] = f.tell()
+                assert(hvals['file_ilevel'] == level+1)
+                if cpu + 1 == self.domain_id:
+                    hydro_offset[level] = f.tell()
+                    level_count[level] = hvals['file_ncache']
                 fpu.skip(f, 8 * self.nvar)
         self._hydro_offset = hydro_offset
+        self._level_count = level_count
         return self._hydro_offset
 
     def _read_amr_header(self):
         fpu.skip(f)
         if hvals['nboundary'] > 0:
             fpu.skip(f, 2)
-            ngridbound = fpu.read_vector(f, 'i')
+            self.ngridbound = fpu.read_vector(f, 'i')
+        else:
+            self.ngridbound = 0
         free_mem = fpu.read_attrs(f, (('free_mem', 5, 'i'), ) )
         ordering = fpu.read_vector(f, 'c')
         fpu.skip(f, 4)
             if c < self.amr_header['ncpu']:
                 ng = self.amr_header['numbl'][l, c]
             else:
-                ng = ngridbound[c - self.amr_header['ncpu'] +
+                ng = self.ngridbound[c - self.amr_header['ncpu'] +
                                 self.amr_header['nboundary']*l]
             return ng
         for level in range(self.amr_header['nlevelmax']):
                 pos[:,0] = fpu.read_vector(f, "d")
                 pos[:,1] = fpu.read_vector(f, "d")
                 pos[:,2] = fpu.read_vector(f, "d")
-                pos *= self.pf.domain_width
+                #pos *= self.pf.domain_width
+                print pos.min(), pos.max()
                 #pos += self.parameter_file.domain_left_edge
                 #print pos.min(), pos.max()
-                parents = fpu.read_vector(f, "I")
-                fpu.skip(f, 6)
-                children = na.empty((ng, 8), dtype='int64')
-                for i in range(8):
-                    children[:,i] = fpu.read_vector(f, "I")
-                cpu_map = na.empty((ng, 8), dtype="int64")
-                for i in range(8):
-                    cpu_map[:,i] = fpu.read_vector(f, "I")
-                rmap = na.empty((ng, 8), dtype="int64")
-                for i in range(8):
-                    rmap[:,i] = fpu.read_vector(f, "I")
+                fpu.skip(f, 31)
+                #parents = fpu.read_vector(f, "I")
+                #fpu.skip(f, 6)
+                #children = na.empty((ng, 8), dtype='int64')
+                #for i in range(8):
+                #    children[:,i] = fpu.read_vector(f, "I")
+                #cpu_map = na.empty((ng, 8), dtype="int64")
+                #for i in range(8):
+                #    cpu_map[:,i] = fpu.read_vector(f, "I")
+                #rmap = na.empty((ng, 8), dtype="int64")
+                #for i in range(8):
+                #    rmap[:,i] = fpu.read_vector(f, "I")
                 # We don't want duplicate grids.
                 if cpu + 1 >= self.domain_id: 
                     assert(pos.shape[0] == ng)
-                    oct_handler.add(cpu + 1, level, ng, pos, ind, cpu_map,
-                                    int(cpu + 1 == self.domain_id))
+                    oct_handler.add(cpu + 1, level, ng, pos, 
+                                    self.domain_id)
 
     def select(self, selector):
         if id(selector) == self._last_selector_id:
         for level, offset in enumerate(self.domain.hydro_offset):
             if offset == -1: continue
             content.seek(offset)
+            nc = self.domain.level_count[level]
+            level_offset = 0
             temp = {}
+            for field in all_fields:
+                temp[field] = na.empty((nc, 8), dtype="float64")
             for i in range(8):
                 for field in all_fields:
                     if field not in fields:
                         #print "Skipping %s in %s : %s" % (field, level,
                         #        self.domain.domain_id)
                         fpu.skip(content)
-                        continue
                     else:
                         #print "Reading %s in %s : %s" % (field, level,
                         #        self.domain.domain_id)
-                        tt = fpu.read_vector(content, 'd') # cell 1
-                        if i == 0:
-                            temp[field] = na.empty((tt.shape[0], 8), dtype="float64")
-                        temp[field][:,i] = tt
-            filled, pos = oct_handler.fill_level(self.domain.domain_id, level,
-                                            tr, temp, self.mask, filled, pos)
+                        temp[field][:,i] = fpu.read_vector(content, 'd') # cell 1
+            oct_handler.fill_level(self.domain.domain_id, level,
+                                   tr, temp, self.mask, level_offset)
             #print "FILL (%s : %s) %s" % (self.domain.domain_id, level, filled)
         #print "DONE (%s) %s of %s" % (self.domain.domain_id, filled,
         #self.cell_count)
             self.parameter_file.domain_right_edge)
         mylog.debug("Allocating %s octs", total_octs)
         self.oct_handler.allocate_domains(
-            [dom.local_oct_count for dom in self.domains])
+            [dom.local_oct_count + dom.ngridbound
+             for dom in self.domains])
         #this actually reads every oct and loads it into the octree
         for dom in self.domains:
             dom._read_amr(self.oct_handler)
 
     def _detect_fields(self):
         # TODO: Add additional fields
-        self.field_list = [ "Density", "x-velocity", "y-velocity",
-	                        "z-velocity", "Pressure", "Metallicity"]
+        self.field_list = ( "Density", "x-velocity", "y-velocity",
+	                        "z-velocity", "Pressure", "Metallicity" )
     
     def _setup_classes(self):
         dd = self._get_data_reader_dict()
                 self.hilbert_indices[int(dom)] = (float(mi), float(ma))
         self.parameters.update(rheader)
         self.current_time = self.parameters['time'] * self.parameters['unit_t']
-        self.domain_right_edge = na.ones(3, dtype='float64') \
-                                           * rheader['boxlen']
         self.domain_left_edge = na.zeros(3, dtype='float64')
-        self.domain_dimensions = na.ones(3, dtype='int32') * 2
+        self.domain_dimensions = na.ones(3, dtype='int32') * 3
+        self.domain_right_edge = na.ones(3, dtype='float64') * self.domain_dimensions
         # This is likely not true, but I am not sure how to otherwise
         # distinguish them.
         mylog.warning("No current mechanism of distinguishing cosmological simulations in RAMSES!")

File yt/geometry/oct_container.pxd

     cdef int nn[3]
     cdef np.float64_t DLE[3], DRE[3]
     cdef public int nocts
-    cdef int max_domain
+    cdef public int max_domain
 
 cdef class RAMSESOctreeContainer(OctreeContainer):
     cdef OctAllocationContainer **domains

File yt/geometry/oct_container.pyx

 cdef class RAMSESOctreeContainer(OctreeContainer):
 
     def allocate_domains(self, domain_counts):
+        print "I AM ALLOCATING", domain_counts
         cdef int count, i
         cdef OctAllocationContainer *cur = self.cont
         assert(cur == NULL)
         self.domains = <OctAllocationContainer **> malloc(
             sizeof(OctAllocationContainer *) * len(domain_counts))
         for i,count in enumerate(domain_counts):
+            print "ALLOCATE", i, count
             cur = allocate_octs(count, cur)
             if self.cont == NULL: self.cont = cur
             self.domains[i] = cur
     @cython.cdivision(True)
     def add(self, int curdom, int curlevel, int ng,
             np.ndarray[np.float64_t, ndim=2] pos,
-            np.ndarray[np.int64_t, ndim=1] findices,
-            np.ndarray[np.int64_t, ndim=2] cpumap,
-            int local):
+            int local_domain):
         cdef int level, no, p, i, j, k, ind[3]
+        cdef int local = (local_domain == curdom)
         cdef Oct* cur = self.root_mesh[0][0][0]
         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
         cdef OctAllocationContainer *cont = self.domains[curdom - 1]
         # How do we bootstrap ourselves?
         for p in range(no):
                 dds[i] = (self.DRE[i] + self.DLE[i])/self.nn[i]
                 ind[i] = <np.int64_t> (pp[i]/dds[i])
                 cp[i] = (ind[i] + 0.5) * dds[i]
+            print ind[0], ind[1], ind[2]
             cur = self.root_mesh[ind[0]][ind[1]][ind[2]]
             if cur == NULL:
                 if curlevel != 0: raise RuntimeError
                     self.nocts += 1
                 cur = next
             cur.domain = curdom
-            if local == 1: cur.ind = p
+            if local == 1:
+                #print cur.local_ind - p
+                cur.ind = p
             cur.level = curlevel
 
     @cython.boundscheck(False)
         return coords
 
     def fill_level(self, int domain, int level, dest_fields, source_fields,
-                   np.ndarray[np.uint8_t, ndim=1, cast=True] mask,
-                   np.int64_t filled, np.int64_t pos):
+                   np.ndarray[np.uint8_t, ndim=1, cast=True] mask, int offset):
         cdef np.ndarray[np.float64_t, ndim=2] source
         cdef np.ndarray[np.float64_t, ndim=1] dest
         cdef OctAllocationContainer *dom = self.domains[domain - 1]
         cdef int i, j, k
         cdef int local_pos, local_filled
         for key in dest_fields:
-            local_filled = filled
+            local_filled = 0
             dest = dest_fields[key]
             source = source_fields[key]
-            # n ranges from the start of the octs in this domain to the end
-            # when we first iterate, it will be at pos = 0.
-            for n in range(dom.n_assigned):
+            for n in range(dom.n):
+                if mask[n + dom.offset] == 0: continue
                 o = &dom.my_octs[n]
-                if o.level != level: continue
-                if mask[n + dom.offset] == 0: continue
                 for i in range(2):
                     for j in range(2):
                         for k in range(2):
                             if o.children[i][j][k] != NULL: continue
                             #print dom.n_assigned, local_filled, dest.shape[0], n, pos, n-pos, source.shape[0], o.level, level
-                            pos = o.ind
-                            dest[local_filled] = source[pos, i*4+j*2+k]
+                            if o.level == level:
+                                #print level, o.ind, offset, o.ind-offset, source.shape[0]
+                                dest[local_filled] = source[o.ind,((k*2)+j)*2+i]
                             local_filled += 1
-        return local_filled, local_pos