[Yt-svn] yt: 3 new changesets
hg at spacepope.org
hg at spacepope.org
Tue Aug 10 10:49:56 PDT 2010
hg Repository: yt
details: yt/rev/479eeff90b61
changeset: 1921:479eeff90b61
user: Matthew Turk <matthewturk at gmail.com>
date:
Tue Aug 10 08:49:11 2010 -0700
description:
Midway through patch-coalescing inside the hierarchy. Then, data IO
coalescing.
hg Repository: yt
details: yt/rev/db8a27305d49
changeset: 1922:db8a27305d49
user: Matthew Turk <matthewturk at gmail.com>
date:
Tue Aug 10 10:03:09 2010 -0700
description:
Intermediate commit; grids now instantiated, but the form of this instantiation
may change. Coalescing is currently extremely inefficient.
hg Repository: yt
details: yt/rev/2760302d04cb
changeset: 1923:2760302d04cb
user: Matthew Turk <matthewturk at gmail.com>
date:
Tue Aug 10 10:49:52 2010 -0700
description:
Changed some things about the calculation of protosubgrids, made things much
much more efficient.
diffstat:
yt/lagos/BaseGridType.py | 10 +--
yt/lagos/HierarchyType.py | 120 ++++++++++++++++++++++++++++++++++------
yt/ramses_reader.pyx | 62 +++++++++++++-------
3 files changed, 145 insertions(+), 47 deletions(-)
diffs (250 lines):
diff -r 6a3c23ef3ae4 -r 2760302d04cb yt/lagos/BaseGridType.py
--- a/yt/lagos/BaseGridType.py Tue Aug 10 07:55:20 2010 -0700
+++ b/yt/lagos/BaseGridType.py Tue Aug 10 10:49:52 2010 -0700
@@ -722,16 +722,14 @@
class RAMSESGrid(AMRGridPatch):
_id_offset = 0
#__slots__ = ["_level_id", "stop_index"]
- def __init__(self, id, hierarchy, level, domain, grid_offset, parent_id,
- cm):
+ def __init__(self, id, hierarchy, level, locations):
AMRGridPatch.__init__(self, id, filename = hierarchy.hierarchy_filename,
hierarchy = hierarchy)
self.Level = level
- self.domain = domain
- self.grid_offset = grid_offset
- self._parent_id = parent_id
+ self._parent_id = -1
self._children_ids = []
- self._child_mask = cm.reshape((2,2,2), order="F")
+ #self._child_mask = cm.reshape((2,2,2), order="F")
+ self.locations = locations
def _del_child_mask(self):
return
diff -r 6a3c23ef3ae4 -r 2760302d04cb yt/lagos/HierarchyType.py
--- a/yt/lagos/HierarchyType.py Tue Aug 10 07:55:20 2010 -0700
+++ b/yt/lagos/HierarchyType.py Tue Aug 10 10:49:52 2010 -0700
@@ -28,6 +28,18 @@
import string, re, gc, time
import cPickle
from itertools import chain, izip
+try:
+ import yt.ramses_reader as ramses_reader
+except ImportError:
+ mylog.warning("Ramses Reader not imported")
+
+def num_deep_inc(f):
+ def wrap(self, *args, **kwargs):
+ self.num_deep += 1
+ rv = f(self, *args, **kwargs)
+ self.num_deep -= 1
+ return rv
+ return wrap
class AMRHierarchy(ObjectFindingMixin, ParallelAnalysisInterface):
float_type = 'float64'
@@ -1542,26 +1554,98 @@
self.object_types.sort()
def _count_grids(self):
- self._level_info = self.tree_proxy.count_zones()
- self.num_grids = sum(self._level_info)
+ # We have to do all the patch-coalescing here.
+ level_info = self.tree_proxy.count_zones()
+ num_ogrids = sum(level_info)
+ ogrid_left_edge = na.zeros((num_ogrids,3), dtype='float64')
+ ogrid_right_edge = na.zeros((num_ogrids,3), dtype='float64')
+ ogrid_levels = na.zeros((num_ogrids,1), dtype='int32')
+ ogrid_file_locations = na.zeros((num_ogrids,3), dtype='int64')
+ ochild_masks = na.zeros((num_ogrids, 8), dtype='int32')
+ self.tree_proxy.fill_hierarchy_arrays(
+ ogrid_left_edge, ogrid_right_edge,
+ ogrid_levels, ogrid_file_locations, ochild_masks)
+ # We now have enough information to run the patch coalescing
+ self.proto_grids = []
+ for level in xrange(len(level_info)):
+ if level_info[level] == 0: continue
+ ggi = (ogrid_levels == level).ravel()
+ left_index = ((ogrid_left_edge[ggi,:]) * (2.0**(level+1))).astype('int64')
+ right_index = left_index + 2
+ dims = na.ones((ggi.sum(), 3), dtype='int64') * 2
+ fl = ogrid_file_locations[ggi,:]
+ # Now our initial protosubgrid
+ initial_left = na.zeros(3, dtype='int64')
+ idims = na.ones(3, dtype='int64') * (2**(level+1))
+ #if level == 6: raise RuntimeError
+ psg = ramses_reader.ProtoSubgrid(initial_left, idims,
+ left_index, right_index, dims, fl)
+ self.proto_grids.append(self._recursive_patch_splitting(
+ psg, idims, initial_left,
+ left_index, right_index, dims, fl))
+ sums = na.zeros(3, dtype='int64')
+ if len(self.proto_grids[level]) == 1: continue
+ for g in self.proto_grids[level]:
+ sums += [s.sum() for s in g.sigs]
+ assert(na.all(sums == dims.prod(axis=1).sum()))
+ self.num_grids = sum(len(l) for l in self.proto_grids)
+
+ #num_deep = 0
+
+ #@num_deep_inc
+ def _recursive_patch_splitting(self, psg, dims, ind,
+ left_index, right_index, gdims, fl):
+ min_eff = 0.2 # This isn't always respected.
+ if psg.efficiency > min_eff or psg.efficiency < 0.0:
+ return [psg]
+ tt, ax, fp = psg.find_split()
+ #print " " * self.num_deep + "Got ax", ax, "fp", fp
+ dims_l = dims.copy()
+ dims_l[ax] = fp
+ li_l = ind.copy()
+ if na.any(dims_l <= 0): return [psg]
+ L = ramses_reader.ProtoSubgrid(
+ li_l, dims_l, left_index, right_index, gdims, fl)
+ #print " " * self.num_deep + "L", tt, L.efficiency
+ if L.efficiency > 1.0: raise RuntimeError
+ if L.efficiency <= 0.0: L = []
+ elif L.efficiency < min_eff:
+ L = self._recursive_patch_splitting(L, dims_l, li_l,
+ left_index, right_index, gdims, fl)
+ else:
+ L = [L]
+ dims_r = dims.copy()
+ dims_r[ax] -= fp
+ li_r = ind.copy()
+ li_r[ax] += fp
+ if na.any(dims_r <= 0): return [psg]
+ R = ramses_reader.ProtoSubgrid(
+ li_r, dims_r, left_index, right_index, gdims, fl)
+ #print " " * self.num_deep + "R", tt, R.efficiency
+ if R.efficiency > 1.0: raise RuntimeError
+ if R.efficiency <= 0.0: R = []
+ elif R.efficiency < min_eff:
+ R = self._recursive_patch_splitting(R, dims_r, li_r,
+ left_index, right_index, gdims, fl)
+ else:
+ R = [R]
+ return L + R
def _parse_hierarchy(self):
- self.grid_dimensions.flat[:] = 2
- grid_file_locations = na.zeros((self.num_grids, 3), dtype='int64')
- child_masks = na.zeros((self.num_grids, 8), dtype='int32')
- self.tree_proxy.fill_hierarchy_arrays(
- self.grid_left_edge, self.grid_right_edge,
- self.grid_levels, grid_file_locations, child_masks)
- ggi = 0
- gs = []
- gs = [self.grid(i, self, self.grid_levels[i,0],
- grid_file_locations[i,0],
- grid_file_locations[i,1],
- grid_file_locations[i,2],
- child_masks[i,:])
- for i in xrange(self.num_grids)]
- self.grids = na.array(gs, dtype='object')
- self.grid_file_locations = grid_file_locations
+ # We have important work to do
+ grids = []
+ gi = 0
+ for level, grid_list in enumerate(self.proto_grids):
+ for g in grid_list:
+ fl = g.grid_file_locations
+ props = g.get_properties()
+ self.grid_left_edge[gi,:] = props[0,:] / (2.0**(level+1))
+ self.grid_right_edge[gi,:] = props[1,:] / (2.0**(level+1))
+ self.grid_dimensions[gi,:] = props[2,:]
+ self.grid_levels[gi,:] = level
+ grids.append(self.grid(gi, self, level, fl))
+ gi += 1
+ self.grids = na.array(grids, dtype='object')
def _populate_grid_objects(self):
for gi,g in enumerate(self.grids):
diff -r 6a3c23ef3ae4 -r 2760302d04cb yt/ramses_reader.pyx
--- a/yt/ramses_reader.pyx Tue Aug 10 07:55:20 2010 -0700
+++ b/yt/ramses_reader.pyx Tue Aug 10 10:49:52 2010 -0700
@@ -604,6 +604,7 @@
cdef public np.float64_t efficiency
cdef public object sigs
cdef public object grid_file_locations
+ cdef public object dd
#@cython.boundscheck(False)
#@cython.wraparound(False)
@@ -683,37 +684,52 @@
if used == 1:
self.grid_file_locations.append(grid_file_locations[gi,:])
- for i in range(3): efficiency /= self.dimensions[i]
+ self.dd = np.ones(3, dtype='int64')
+ for i in range(3):
+ efficiency /= self.dimensions[i]
+ self.dd[i] = self.dimensions[i]
#print "Efficiency is %0.3e" % (efficiency)
self.efficiency = efficiency
- def find_split(self, int ax):
+ def find_split(self):
# First look for zeros
- cdef int i, center
- center = self.dimensions[i] / 2
+ cdef int i, center, ax
+ cdef np.ndarray[ndim=1, dtype=np.int64_t] axes
cdef np.int64_t strength, zcstrength, zcp
- cdef np.ndarray[np.int64_t] sig = self.sigs[ax]
- for i in range(self.dimensions[ax]):
- if sig[i] == 0:
- #print "zero: %s (%s)" % (i, self.dimensions[ax])
- return i
- cdef np.int64_t *sig2d = <np.int64_t *> calloc(
- sizeof(np.int64_t), self.dimensions[ax])
- sig2d[0] = sig2d[self.dimensions[ax]-1] = 0
- for i in range(1, self.dimensions[ax] - 1):
- sig2d[i] = (sig[i-1] - 2*sig[i] + sig[i+1])
+ axes = np.argsort(self.dd)[::-1]
+ cdef np.ndarray[np.int64_t] sig
+ for axi in range(3):
+ ax = axes[axi]
+ center = self.dimensions[ax] / 2
+ sig = self.sigs[ax]
+ for i in range(self.dimensions[ax]):
+ if sig[i] == 0 and i > 0 and i < self.dimensions[ax] - 1:
+ #print "zero: %s (%s)" % (i, self.dimensions[ax])
+ return 'zs', ax, i
zcstrength = 0
zcp = 0
- for i in range(1, self.dimensions[ax] - 1):
- if sig2d[i] * sig2d[i+1] <= 0:
- strength = labs(sig2d[i] - sig2d[i+1])
- if (strength > zcstrength) or \
- (strength == zcstrength and (abs(center - i) <
- abs(center - zcp))):
- zcstrength = strength
- zcp = i
+ zca = -1
+ cdef int temp
+ cdef np.int64_t *sig2d
+ for axi in range(3):
+ ax = axes[axi]
+ sig = self.sigs[ax]
+ sig2d = <np.int64_t *> malloc(sizeof(np.int64_t) * self.dimensions[ax])
+ sig2d[0] = sig2d[self.dimensions[ax]-1] = 0
+ for i in range(1, self.dimensions[ax] - 1):
+ sig2d[i] = sig[i-1] - 2*sig[i] + sig[i+1]
+ for i in range(1, self.dimensions[ax] - 1):
+ if sig2d[i] * sig2d[i+1] <= 0:
+ strength = labs(sig2d[i] - sig2d[i+1])
+ if (strength > zcstrength) or \
+ (strength == zcstrength and (abs(center - i) <
+ abs(center - zcp))):
+ zcstrength = strength
+ zcp = i
+ zca = ax
+ free(sig2d)
#print "zcp: %s (%s)" % (zcp, self.dimensions[ax])
- return zcp
+ return 'zc', ax, zcp
def get_properties(self):
cdef np.ndarray[np.int64_t, ndim=2] tr = np.empty((3,3), dtype='int64')
More information about the yt-svn
mailing list