[Yt-svn] yt: 2 new changesets
hg at spacepope.org
hg at spacepope.org
Mon Aug 16 21:34:29 PDT 2010
hg Repository: yt
details: yt/rev/863c57567671
changeset: 1947:863c57567671
user: Christopher Erick Moody <cemoody at ucsc.edu>
date:
Mon Aug 16 21:32:50 2010 -0700
description:
Debugging Gadget components
hg Repository: yt
details: yt/rev/fc315090030c
changeset: 1948:fc315090030c
user: Christopher Erick Moody <cemoody at ucsc.edu>
date:
Mon Aug 16 21:34:18 2010 -0700
description:
Debugged Gadget and merged
diffstat:
yt/lagos/BaseGridType.py | 15 +++-
yt/lagos/GadgetFields.py | 15 +++--
yt/lagos/HierarchyType.py | 116 ++++++++++++++++++++++++++++++--------
yt/lagos/OutputTypes.py | 12 ++-
yt/lagos/__init__.py | 1 +
yt/ramses_reader.pyx | 90 +++++++++++++++++++----------
6 files changed, 178 insertions(+), 71 deletions(-)
diffs (truncated from 578 to 300 lines):
diff -r f76765bcffe4 -r fc315090030c yt/lagos/BaseGridType.py
--- a/yt/lagos/BaseGridType.py Wed Aug 11 10:09:12 2010 -0700
+++ b/yt/lagos/BaseGridType.py Mon Aug 16 21:34:18 2010 -0700
@@ -27,6 +27,8 @@
#import yt.enki, gc
from yt.funcs import *
+import pdb
+
class AMRGridPatch(object):
_spatial = True
_num_ghost_zones = 0
@@ -362,7 +364,7 @@
mask[startIndex[0]:endIndex[0],
startIndex[1]:endIndex[1],
startIndex[2]:endIndex[2]] = tofill
-
+
def __generate_child_mask(self):
"""
Generates self.child_mask, which is zero where child grids exist (and
@@ -598,8 +600,9 @@
self.N = 0
self.Address = ''
self.NumberOfParticles = self.N
- self.ActiveDimensions = [0,0,0]
+ self.ActiveDimensions = na.array([0,0,0])
self._id_offset = 0
+ self.start_index = na.array([0,0,0])
for key,val in kwargs.items():
if key in dir(self):
@@ -609,6 +612,9 @@
#def __repr__(self):
# return "GadgetGrid_%04i" % (self.Address)
+ def get_global_startindex(self):
+ return self.start_index
+
def _prepare_grid(self):
#all of this info is already included in the snapshots
pass
@@ -621,12 +627,13 @@
# So first we figure out what the index is. We don't assume
# that dx=dy=dz , at least here. We probably do elsewhere.
id = self.id
- LE, RE = self.hierarchy.grid_left_edge[id,:], \
- self.hierarchy.grid_right_edge[id,:]
+ LE, RE = self.LeftEdge,self.RightEdge
self.dds = na.array((RE-LE)/self.ActiveDimensions)
if self.pf["TopGridRank"] < 2: self.dds[1] = 1.0
if self.pf["TopGridRank"] < 3: self.dds[2] = 1.0
self.data['dx'], self.data['dy'], self.data['dz'] = self.dds
+
+
class ChomboGrid(AMRGridPatch):
diff -r f76765bcffe4 -r fc315090030c yt/lagos/GadgetFields.py
--- a/yt/lagos/GadgetFields.py Wed Aug 11 10:09:12 2010 -0700
+++ b/yt/lagos/GadgetFields.py Mon Aug 16 21:34:18 2010 -0700
@@ -38,11 +38,14 @@
"Mass":"MASS"
}
-for f,v in translation_dict.items():
- if v not in GadgetFieldInfo:
- add_field(f, function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("POS")],
- units=r"\rm{cm}")
+#for f,v in translation_dict.items():
+# add_field(f, function=lambda a,b: None, take_log=True,
+# validators = [ValidateDataField(v)],
+# units=r"\rm{cm}")
+# add_field(v, function=lambda a,b: None, take_log=True,
+# validators = [ValidateDataField(v)],
+# units=r"\rm{cm}")
+
add_field("Density", function=lambda a,b: None, take_log=True,
validators = [ValidateDataField("POS")],
@@ -50,7 +53,7 @@
add_field("VEL", function=lambda a,b: None, take_log=True,
validators = [ValidateDataField("VEL")],
- units=r"\rm{cm}/\rm{s}")
+ units=r"")
add_field("ID", function=lambda a,b: None, take_log=True,
validators = [ValidateDataField("ID")],
diff -r f76765bcffe4 -r fc315090030c yt/lagos/HierarchyType.py
--- a/yt/lagos/HierarchyType.py Wed Aug 11 10:09:12 2010 -0700
+++ b/yt/lagos/HierarchyType.py Mon Aug 16 21:34:18 2010 -0700
@@ -1201,6 +1201,7 @@
self.ngrids = ngrids
self.grids = []
+
class GadgetHierarchy(AMRHierarchy):
grid = GadgetGrid
@@ -1219,7 +1220,7 @@
#example string:
#"(S'VEL'\np1\nS'ID'\np2\nS'MASS'\np3\ntp4\n."
#fields are surrounded with '
- fields_string=self._handle['root'].attrs['names']
+ fields_string=self._handle['root'].attrs['fieldnames']
#splits=fields_string.split("'")
#pick out the odd fields
#fields= [splits[j] for j in range(1,len(splits),2)]
@@ -1263,8 +1264,8 @@
kwargs = {}
kwargs['Address'] = loc
- kwargs['Children'] = [ch for ch in node.values()]
kwargs['Parent'] = parent
+ kwargs['Axis'] = self.pf._get_param('divideaxis',location=loc)
kwargs['Level'] = self.pf._get_param('level',location=loc)
kwargs['LeftEdge'] = self.pf._get_param('leftedge',location=loc)
kwargs['RightEdge'] = self.pf._get_param('rightedge',location=loc)
@@ -1274,22 +1275,31 @@
dx = self.pf._get_param('dx',location=loc)
dy = self.pf._get_param('dy',location=loc)
dz = self.pf._get_param('dz',location=loc)
- kwargs['ActiveDimensions'] = (dx,dy,dz)
+ divdims = na.array([1,1,1]
+ if not kwargs['IsLeaf']:
+ divdims[kwargs['Axis']] = 2
+ kwargs['ActiveDimensions'] = divdims
+ #Active dimensions:
+ #This is the number of childnodes, along with dimensiolaity
+ #ie, binary tree can be (2,1,1) but octree is (2,2,2)
+
+ idx+=1
+ #pdb.set_trace()
+ children = []
+ if not kwargs['IsLeaf']:
+ for child in node.values():
+ children,idx=self._walk_nodes(node,child,children,idx=idx)
+
+ kwargs['Children'] = children
grid = self.grid(idx,self.pf.parameter_filename,self,**kwargs)
- idx+=1
+ grids += children
grids += [grid,]
- #pdb.set_trace()
- if kwargs['IsLeaf']:
- return grids,idx
- else:
- for child in node.values():
- grids,idx=self._walk_nodes(node,child,grids,idx=idx)
return grids,idx
-
+
def _populate_grid_objects(self):
for g in self.grids:
g._prepare_grid()
- self.max_level = self._handle['root'].attrs['maxlevel']
+ self.max_level = cPickle.loads(self._handle['root'].attrs['maxlevel'])
def _setup_unknown_fields(self):
pass
@@ -1665,7 +1675,13 @@
self.pf["TopGridDimensions"],
ogrid_left_edge, ogrid_right_edge,
ogrid_levels, ogrid_file_locations, ochild_masks)
- import pdb;pdb.set_trace()
+ # Now we can rescale
+ mi, ma = ogrid_left_edge.min(), ogrid_right_edge.max()
+ DL = self.pf["DomainLeftEdge"]
+ DR = self.pf["DomainRightEdge"]
+ ogrid_left_edge = (ogrid_left_edge - mi)/(ma - mi) * (DR - DL) + DL
+ ogrid_right_edge = (ogrid_right_edge - mi)/(ma - mi) * (DR - DL) + DL
+ #import pdb;pdb.set_trace()
# We now have enough information to run the patch coalescing
self.proto_grids = []
for level in xrange(len(level_info)):
@@ -1673,19 +1689,65 @@
ggi = (ogrid_levels == level).ravel()
mylog.info("Re-gridding level %s: %s octree grids", level, ggi.sum())
nd = self.pf["TopGridDimensions"] * 2**level
- left_index = na.rint((ogrid_left_edge[ggi,:]) * nd).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') * nd
#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))
+ # We want grids that cover no more than MAX_EDGE cells in every direction
+ MAX_EDGE = 128
+ psgs = []
+ left_index = na.rint((ogrid_left_edge[ggi,:]) * nd).astype('int64')
+ right_index = left_index + 2
+ lefts = [na.mgrid[0:nd[i]:MAX_EDGE] for i in range(3)]
+ #lefts = zip(*[l.ravel() for l in lefts])
+ pbar = get_pbar("Re-gridding ", lefts[0].size)
+ min_ind = na.min(left_index, axis=0)
+ max_ind = na.max(right_index, axis=0)
+ for i,dli in enumerate(lefts[0]):
+ pbar.update(i)
+ if min_ind[0] > dli + nd[0]: continue
+ if max_ind[0] < dli: continue
+ idim = min(nd[0] - dli, MAX_EDGE)
+ gdi = ((dli <= right_index[:,0])
+ & (dli + idim >= left_index[:,0]))
+ if not na.any(gdi): continue
+ for dlj in lefts[1]:
+ if min_ind[1] > dlj + nd[1]: continue
+ if max_ind[1] < dlj: continue
+ idim = min(nd[1] - dlj, MAX_EDGE)
+ gdj = ((dlj <= right_index[:,1])
+ & (dlj + idim >= left_index[:,1])
+ & (gdi))
+ if not na.any(gdj): continue
+ for dlk in lefts[2]:
+ if min_ind[2] > dlk + nd[2]: continue
+ if max_ind[2] < dlk: continue
+ idim = min(nd[2] - dlk, MAX_EDGE)
+ gdk = ((dlk <= right_index[:,2])
+ & (dlk + idim >= left_index[:,2])
+ & (gdj))
+ if not na.any(gdk): continue
+ left = na.array([dli, dlj, dlk])
+ domain_left = left.ravel()
+ initial_left = na.zeros(3, dtype='int64') + domain_left
+ idims = na.ones(3, dtype='int64') * na.minimum(nd - domain_left, MAX_EDGE)
+ # We want to find how many grids are inside.
+ dleft_index = left_index[gdk,:]
+ dright_index = right_index[gdk,:]
+ ddims = dims[gdk,:]
+ dfl = fl[gdk,:]
+ psg = ramses_reader.ProtoSubgrid(initial_left, idims,
+ dleft_index, dright_index, ddims, dfl)
+ #print "Gridding from %s to %s + %s" % (
+ # initial_left, initial_left, idims)
+ if psg.efficiency <= 0: continue
+ self.num_deep = 0
+ psgs.extend(self._recursive_patch_splitting(
+ psg, idims, initial_left,
+ dleft_index, dright_index, ddims, dfl))
+ #psgs.extend([psg])
+ pbar.finish()
+ self.proto_grids.append(psgs)
sums = na.zeros(3, dtype='int64')
mylog.info("Final grid count: %s", len(self.proto_grids[level]))
if len(self.proto_grids[level]) == 1: continue
@@ -1694,12 +1756,16 @@
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 = 0
- #@num_deep_inc
+ @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.
+ min_eff = 0.1 # This isn't always respected.
+ if self.num_deep > 40:
+ # If we've recursed more than 100 times, we give up.
+ psg.efficiency = min_eff
+ return [psg]
if psg.efficiency > min_eff or psg.efficiency < 0.0:
return [psg]
tt, ax, fp = psg.find_split()
diff -r f76765bcffe4 -r fc315090030c yt/lagos/OutputTypes.py
--- a/yt/lagos/OutputTypes.py Wed Aug 11 10:09:12 2010 -0700
+++ b/yt/lagos/OutputTypes.py Mon Aug 16 21:34:18 2010 -0700
@@ -637,6 +637,8 @@
# generalization.
self.parameters["TopGridRank"] = 3
self.parameters["RefineBy"] = 2
+ self.parameters["DomainLeftEdge"] = self.leftedge
+ self.parameters["DomainRightEdge"] = self.rightedge
def _parse_parameter_file(self):
@@ -654,7 +656,7 @@
continue
val = fh['root'].attrs[kw]
if type(val)==type(''):
- try: val = cPickle.load(val)
+ try: val = cPickle.loads(val)
except: pass
#also, includes unit info
setattr(self,kw,val)
@@ -662,9 +664,8 @@
def _get_param(self,kw,location='/root'):
fh = h5py.File(self.parameter_filename)
val = fh[location].attrs[kw]
- if type(val)==type(''):
- try: val = cPickle.load(val)
- except: pass
More information about the yt-svn
mailing list