[yt-svn] commit/yt-3.0: MatthewTurk: RAMSES data now reads in correctly, constructing the correct number of Octs and
Bitbucket
commits-noreply at bitbucket.org
Fri Aug 10 15:08:11 PDT 2012
1 new commit in yt-3.0:
https://bitbucket.org/yt_analysis/yt-3.0/changeset/f75f82fca0a8/
changeset: f75f82fca0a8
branch: yt-3.0
user: MatthewTurk
date: 2012-08-11 00:08:02
summary: RAMSES data now reads in correctly, constructing the correct number of Octs and
so on. Not all data region selections work; 3D seems to, but 2D is not working
correctly as of yet. Single domain can do IO, multi cannot yet.
affected #: 5 files
diff -r e5e4c9adfee291745cfcb3aa7dab109874e46de4 -r f75f82fca0a8087ca05a2cd2269fa4d351f76767 yt/frontends/ramses/data_structures.py
--- a/yt/frontends/ramses/data_structures.py
+++ b/yt/frontends/ramses/data_structures.py
@@ -85,9 +85,11 @@
f = open(self.hydro_fn, "rb")
fpu.skip(f, 6)
# It goes: level, CPU, 8-variable
- hydro_offset = na.zeros(self.amr_header['nlevelmax'], dtype='int64')
+ min_level = self.pf.min_level
+ n_levels = self.amr_header['nlevelmax'] - min_level
+ hydro_offset = na.zeros(n_levels, dtype='int64')
hydro_offset -= 1
- level_count = na.zeros(self.amr_header['nlevelmax'], dtype='int64')
+ level_count = na.zeros(n_levels, dtype='int64')
for level in range(self.amr_header['nlevelmax']):
for cpu in range(self.amr_header['nboundary'] +
self.amr_header['ncpu']):
@@ -96,9 +98,9 @@
hvals = fpu.read_attrs(f, header)
if hvals['file_ncache'] == 0: continue
assert(hvals['file_ilevel'] == level+1)
- if cpu + 1 == self.domain_id:
- hydro_offset[level] = f.tell()
- level_count[level] = hvals['file_ncache']
+ if cpu + 1 == self.domain_id and level >= min_level:
+ hydro_offset[level - min_level] = f.tell()
+ level_count[level - min_level] = hvals['file_ncache']
fpu.skip(f, 8 * self.nvar)
self._hydro_offset = hydro_offset
self._level_count = level_count
@@ -125,7 +127,7 @@
# Now we iterate over each level and each CPU.
self.amr_header = hvals
self.amr_offset = f.tell()
- self.local_oct_count = hvals['numbl'][:, self.domain_id - 1].sum()
+ self.local_oct_count = hvals['numbl'][self.pf.min_level:, self.domain_id - 1].sum()
def _read_amr(self, oct_handler):
"""Open the oct file, read in octs level-by-level.
@@ -139,8 +141,8 @@
f = cStringIO.StringIO()
f.write(fb.read())
f.seek(0)
- mylog.debug("Reading domain AMR % 4i (%0.3e)",
- self.domain_id, self.local_oct_count)
+ mylog.debug("Reading domain AMR % 4i (%0.3e, %0.3e)",
+ self.domain_id, self.local_oct_count, self.ngridbound)
def _ng(c, l):
if c < self.amr_header['ncpu']:
ng = self.amr_header['numbl'][l, c]
@@ -148,6 +150,8 @@
ng = self.ngridbound[c - self.amr_header['ncpu'] +
self.amr_header['nboundary']*l]
return ng
+ min_level = self.pf.min_level
+ total = 0
for level in range(self.amr_header['nlevelmax']):
# Easier if do this 1-indexed
for cpu in range(self.amr_header['nboundary'] + self.amr_header['ncpu']):
@@ -155,16 +159,13 @@
ng = _ng(cpu, level)
if ng == 0: continue
ind = fpu.read_vector(f, "I").astype("int64")
- #print level, cpu, ind.min(), ind.max(), ind.size
fpu.skip(f, 2)
pos = na.empty((ng, 3), dtype='float64')
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
- print pos.min(), pos.max()
#pos += self.parameter_file.domain_left_edge
- #print pos.min(), pos.max()
fpu.skip(f, 31)
#parents = fpu.read_vector(f, "I")
#fpu.skip(f, 6)
@@ -178,9 +179,12 @@
#for i in range(8):
# rmap[:,i] = fpu.read_vector(f, "I")
# We don't want duplicate grids.
- if cpu + 1 >= self.domain_id:
+ # Note that we're adding *grids*, not individual cells.
+ if level >= min_level and cpu + 1 >= self.domain_id:
assert(pos.shape[0] == ng)
- oct_handler.add(cpu + 1, level, ng, pos,
+ if cpu + 1 == self.domain_id:
+ total += ng
+ oct_handler.add(cpu + 1, level - min_level, ng, pos,
self.domain_id)
def select(self, selector):
@@ -227,6 +231,7 @@
fields = [f for ft, f in fields]
tr = {}
filled = pos = 0
+ min_level = self.domain.pf.min_level
for field in fields:
tr[field] = na.zeros(self.cell_count, 'float64')
for level, offset in enumerate(self.domain.hydro_offset):
@@ -247,7 +252,8 @@
#print "Reading %s in %s : %s" % (field, level,
# self.domain.domain_id)
temp[field][:,i] = fpu.read_vector(content, 'd') # cell 1
- oct_handler.fill_level(self.domain.domain_id, level,
+ oct_handler.fill_level(self.domain.domain_id,
+ level - min_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,
@@ -275,7 +281,7 @@
#this merely allocates space for the oct tree
#and nothing else
self.oct_handler = RAMSESOctreeContainer(
- self.domains[0].amr_header['nx'],
+ self.parameter_file.domain_dimensions/2,
self.parameter_file.domain_left_edge,
self.parameter_file.domain_right_edge)
mylog.debug("Allocating %s octs", total_octs)
@@ -285,6 +291,8 @@
#this actually reads every oct and loads it into the octree
for dom in self.domains:
dom._read_amr(self.oct_handler)
+ #for dom in self.domains:
+ # self.oct_handler.check(dom.domain_id)
def _detect_fields(self):
# TODO: Add additional fields
@@ -399,7 +407,13 @@
for i in range(11): read_rhs(float)
f.readline()
read_rhs(str)
- # Now we read the hilber indices
+ # This next line deserves some comment. We specify a min_level that
+ # corresponds to the minimum level in the RAMSES simulation. RAMSES is
+ # one-indexed, but it also does refer to the *oct* dimensions -- so
+ # this means that a levelmin of 1 would have *1* oct in it. So a
+ # levelmin of 2 would have 8 octs at the root mesh level.
+ self.min_level = rheader['levelmin'] - 1
+ # Now we read the hilbert indices
self.hilbert_indices = {}
if rheader['ordering type'] == "hilbert":
f.readline() # header
@@ -409,8 +423,9 @@
self.parameters.update(rheader)
self.current_time = self.parameters['time'] * self.parameters['unit_t']
self.domain_left_edge = na.zeros(3, dtype='float64')
- self.domain_dimensions = na.ones(3, dtype='int32') * 3
- self.domain_right_edge = na.ones(3, dtype='float64') * self.domain_dimensions
+ self.domain_dimensions = na.ones(3, dtype='int32') * \
+ 2**(self.min_level+1)
+ self.domain_right_edge = na.ones(3, dtype='float64')
# 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!")
@@ -419,7 +434,7 @@
self.omega_lambda = rheader["omega_l"]
self.omega_matter = rheader["omega_m"]
self.hubble_constant = rheader["H0"]
- self.max_level = rheader['levelmax']
+ self.max_level = rheader['levelmax'] - rheader['levelmin']
@classmethod
def _is_valid(self, *args, **kwargs):
diff -r e5e4c9adfee291745cfcb3aa7dab109874e46de4 -r f75f82fca0a8087ca05a2cd2269fa4d351f76767 yt/geometry/oct_container.pxd
--- a/yt/geometry/oct_container.pxd
+++ b/yt/geometry/oct_container.pxd
@@ -31,7 +31,7 @@
np.int64_t local_ind
np.int64_t domain # (opt) addl int index
np.int64_t pos[3] # position in ints
- np.uint8_t level
+ np.int8_t level
Oct *children[2][2][2]
Oct *parent
diff -r e5e4c9adfee291745cfcb3aa7dab109874e46de4 -r f75f82fca0a8087ca05a2cd2269fa4d351f76767 yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -46,6 +46,8 @@
else:
n_cont.offset = prev.offset + prev.n
n_cont.my_octs = <Oct *> malloc(sizeof(Oct) * n_octs)
+ if n_cont.my_octs == NULL:
+ raise MemoryError
n_cont.n = n_octs
n_cont.n_assigned = 0
for n in range(n_octs):
@@ -53,6 +55,7 @@
oct.parent = NULL
oct.ind = oct.domain = -1
oct.local_ind = -1
+ oct.level = -1
for i in range(2):
for j in range(2):
for k in range(2):
@@ -150,7 +153,7 @@
o = &cur.my_octs[oi - cur.offset]
for i in range(3):
# This gives the *grid* width for this level
- dx[i] = base_dx[i] / (2 << o.level)
+ dx[i] = base_dx[i] / (1 << o.level)
pos[i] = self.DLE[i] + o.pos[i]*dx[i] + dx[i]/4.0
dx[i] = dx[i] / 2.0 # This is now the *offset*
for i in range(2):
@@ -169,15 +172,13 @@
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.max_domain = len(domain_counts) # 1-indexed
self.domains = <OctAllocationContainer **> malloc(
sizeof(OctAllocationContainer *) * len(domain_counts))
- for i,count in enumerate(domain_counts):
- print "ALLOCATE", i, count
+ for i, count in enumerate(domain_counts):
cur = allocate_octs(count, cur)
if self.cont == NULL: self.cont = cur
self.domains[i] = cur
@@ -205,6 +206,25 @@
count[cur.my_octs[i - cur.offset].domain - 1] += 1
return count
+ def check(self, int curdom):
+ cdef int dind, pi
+ cdef Oct oct
+ cdef OctAllocationContainer *cont = self.domains[curdom - 1]
+ cdef int nbad = 0
+ for pi in range(cont.n_assigned):
+ oct = cont.my_octs[pi]
+ for i in range(2):
+ for j in range(2):
+ for k in range(2):
+ if oct.children[i][j][k] != NULL and \
+ oct.children[i][j][k].level != oct.level + 1:
+ if curdom == 61:
+ print pi, oct.children[i][j][k].level,
+ print oct.level
+ nbad += 1
+ print "DOMAIN % 3i HAS % 9i BAD OCTS (%s / %s / %s)" % (curdom, nbad,
+ cont.n - cont.n_assigned, cont.n_assigned, cont.n)
+
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@@ -213,11 +233,12 @@
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 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 > self.max_domain: curdom = local_domain
cdef OctAllocationContainer *cont = self.domains[curdom - 1]
+ cdef int initial = cont.n_assigned
# How do we bootstrap ourselves?
for p in range(no):
#for every oct we're trying to add find the
@@ -227,10 +248,11 @@
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
+ if curlevel != 0:
+ raise RuntimeError
+ assert(cont.n_assigned < cont.n)
cur = &cont.my_octs[cont.n_assigned]
cur.local_ind = cont.offset + cont.n_assigned
cur.parent = NULL
@@ -240,6 +262,7 @@
cont.n_assigned += 1
self.nocts += 1
self.root_mesh[ind[0]][ind[1]][ind[2]] = cur
+ assert(cur.level == curlevel)
# Now we find the location we want
# Note that RAMSES I think 1-findiceses levels, but we don't.
for level in range(curlevel):
@@ -255,25 +278,26 @@
ind[i] = 1
cp[i] += dds[i]/2.0
# Check if it has not been allocated
- #print ind[0], ind[1], ind[2]
next = cur.children[ind[0]][ind[1]][ind[2]]
if next == NULL:
- cur.children[ind[0]][ind[1]][ind[2]] = \
- &cont.my_octs[cont.n_assigned]
- next = cur.children[ind[0]][ind[1]][ind[2]]
+ if cont.n_assigned >= cont.n:
+ raise AssertionError
+ next = &cont.my_octs[cont.n_assigned]
+ cont.n_assigned += 1
+ cur.children[ind[0]][ind[1]][ind[2]] = next
next.local_ind = cont.offset + cont.n_assigned
next.parent = cur
for i in range(3):
next.pos[i] = ind[i] + (cur.pos[i] << 1)
next.level = level + 1
- cont.n_assigned += 1
self.nocts += 1
+ assert(next.level == curlevel)
cur = next
cur.domain = curdom
if local == 1:
- #print cur.local_ind - p
cur.ind = p
cur.level = curlevel
+ return cont.n_assigned - initial
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -346,7 +370,7 @@
for i in range(3):
# This is the base_dx, but not the base distance from the center
# position. Note that the positions will also all be offset by
- # dx/2.0.
+ # dx/2.0. This is also for *oct grids*, not cells.
base_dx[i] = (self.DRE[i] - self.DLE[i])/self.nn[i]
ci = 0
for oi in range(cur.n):
@@ -354,7 +378,7 @@
o = &cur.my_octs[oi]
for i in range(3):
# This gives the *grid* width for this level
- dx[i] = base_dx[i] / (2 << o.level)
+ dx[i] = base_dx[i] / (1 << o.level)
pos[i] = self.DLE[i] + o.pos[i]*dx[i] + dx[i]/4.0
dx[i] = dx[i] / 2.0 # This is now the *offset*
for i in range(2):
@@ -387,8 +411,6 @@
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
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
diff -r e5e4c9adfee291745cfcb3aa7dab109874e46de4 -r f75f82fca0a8087ca05a2cd2269fa4d351f76767 yt/geometry/oct_geometry_handler.py
--- a/yt/geometry/oct_geometry_handler.py
+++ b/yt/geometry/oct_geometry_handler.py
@@ -53,7 +53,8 @@
"""
Returns (in code units) the smallest cell size in the simulation.
"""
- raise NotImplementedError
+ return (self.parameter_file.domain_width /
+ (2**self.max_level)).min()
def convert(self, unit):
return self.parameter_file.conversion_factors[unit]
diff -r e5e4c9adfee291745cfcb3aa7dab109874e46de4 -r f75f82fca0a8087ca05a2cd2269fa4d351f76767 yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -172,6 +172,7 @@
for j in range(octree.nn[1]):
pos[2] = octree.DLE[2] + dds[2]/2.0
for k in range(octree.nn[2]):
+ if octree.root_mesh[i][j][k] == NULL: continue
self.recursively_select_octs(
octree.root_mesh[i][j][k],
pos, dds, mask)
@@ -390,6 +391,7 @@
cdef int select_grid(self, np.float64_t left_edge[3],
np.float64_t right_edge[3]) nogil:
cdef np.float64_t box_center, relcenter, closest, dist, edge
+ return 1
cdef int id
if (left_edge[0] <= self.center[0] <= right_edge[0] and
left_edge[1] <= self.center[1] <= right_edge[1] and
Repository URL: https://bitbucket.org/yt_analysis/yt-3.0/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list