[Yt-svn] yt-commit r588 - branches/yt-generalization/yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Fri Jun 20 16:06:01 PDT 2008
Author: mturk
Date: Fri Jun 20 16:05:57 2008
New Revision: 588
URL: http://yt.spacepope.org/changeset/588
Log:
Split the hierarchy up a bit. I am still astounded by how ugly this code --
much of it untouched since 2006 -- is.
Passes all of the unit tests.
Modified:
branches/yt-generalization/yt/lagos/HierarchyType.py
Modified: branches/yt-generalization/yt/lagos/HierarchyType.py
==============================================================================
--- branches/yt-generalization/yt/lagos/HierarchyType.py (original)
+++ branches/yt-generalization/yt/lagos/HierarchyType.py Fri Jun 20 16:05:57 2008
@@ -36,72 +36,47 @@
}
class AMRHierarchy:
- pass
-
-class EnzoHierarchy(AMRHierarchy):
- eiTopGrid = None
- _strip_path = False
- @time_execution
- def __init__(self, pf, data_style=None):
- """
- This is the grid structure as Enzo sees it, with some added bonuses.
- It's primarily used as a class factory, to generate data objects and
- access grids.
-
- It should never be created directly -- you should always access it via
- calls to an affiliated :class:`~yt.lagos.EnzoStaticOutput`.
-
- On instantiation, it processes the hierarchy and generates the grids.
- """
- # Expect filename to be the name of the parameter file, not the
- # hierarchy
- self.hierarchy_filename = os.path.abspath(pf.parameter_filename) \
- + ".hierarchy"
- self.boundary_filename = os.path.abspath(pf.parameter_filename) \
- + ".boundary"
- self.directory = os.path.dirname(self.hierarchy_filename)
+ def __init__(self, pf):
self.parameter_file = weakref.proxy(pf)
- self.__data_file = None
- # Now we search backwards from the end of the file to find out how many
- # grids we have, which allows us to preallocate memory
- self.__hierarchy_lines = open(self.hierarchy_filename).readlines()
- if len(self.__hierarchy_lines) == 0:
- raise IOError(-1,"File empty", self.hierarchy_filename)
- self.__hierarchy_string = open(self.hierarchy_filename).read()
- for line in reversed(self.__hierarchy_lines):
- if line.startswith("Grid ="):
- self.num_grids = int(line.split("=")[-1])
- break
- self.data_style = data_style
- self.__guess_data_style()
+ self._data_file = None
+ self._setup_classes()
+ self._initialize_grids()
- # For some reason, r8 seems to want Float64
- if self.parameters.has_key("CompilerPrecision") \
- and self.parameters["CompilerPrecision"] == "r4":
- EnzoFloatType = 'float32'
- else:
- EnzoFloatType = 'float64'
+ # For use with derived quantities depending on centers
+ # Although really, I think perhaps we should take a closer look
+ # at how "center" is used.
+ self.center = None
+ self.bulkVelocity = None
+ self._initialize_level_stats()
+
+ mylog.debug("Initializing data file")
+ self._initialize_data_file()
+ mylog.debug("Populating hierarchy")
+ self._populate_hierarchy()
+ mylog.debug("Done populating hierarchy")
+
+ def _initialize_grids(self):
mylog.debug("Allocating memory for %s grids", self.num_grids)
- self.__setup_classes()
self.gridDimensions = na.zeros((self.num_grids,3), 'int32')
self.gridStartIndices = na.zeros((self.num_grids,3), 'int32')
self.gridEndIndices = na.zeros((self.num_grids,3), 'int32')
- self.gridLeftEdge = na.zeros((self.num_grids,3), EnzoFloatType)
- self.gridRightEdge = na.zeros((self.num_grids,3), EnzoFloatType)
+ self.gridLeftEdge = na.zeros((self.num_grids,3), self.float_type)
+ self.gridRightEdge = na.zeros((self.num_grids,3), self.float_type)
self.gridLevels = na.zeros((self.num_grids,1), 'int32')
- self.gridDxs = na.zeros((self.num_grids,1), EnzoFloatType)
- self.gridDys = na.zeros((self.num_grids,1), EnzoFloatType)
- self.gridDzs = na.zeros((self.num_grids,1), EnzoFloatType)
+ self.gridDxs = na.zeros((self.num_grids,1), self.float_type)
+ self.gridDys = na.zeros((self.num_grids,1), self.float_type)
+ self.gridDzs = na.zeros((self.num_grids,1), self.float_type)
self.gridTimes = na.zeros((self.num_grids,1), 'float64')
self.gridNumberOfParticles = na.zeros((self.num_grids,1))
mylog.debug("Done allocating")
-
+ mylog.debug("Creating grid objects")
self.grids = na.array([self.grid(i+1) for i in xrange(self.num_grids)])
- mylog.debug("Done creating grid objects")
self.gridReverseTree = [-1] * self.num_grids
self.gridTree = [ [] for i in range(self.num_grids)]
+ mylog.debug("Done creating grid objects")
+ def _initialize_level_stats(self):
# Now some statistics:
# 0 = number of grids
# 1 = number of cells
@@ -113,79 +88,8 @@
self.level_stats['numgrids'] = [0 for i in range(MAXLEVEL)]
self.level_stats['numcells'] = [0 for i in range(MAXLEVEL)]
- # For use with derived quantities depending on centers
- # Although really, I think perhaps we should take a closer look
- # at how "center" is used.
- self.center = None
- self.bulkVelocity = None
-
- mylog.debug("Initializing data file")
- self.__initialize_data_file()
- mylog.debug("Populating hierarchy")
- self.__populate_hierarchy()
- mylog.debug("Done populating hierarchy")
-
- def __guess_data_style(self):
- if self.data_style: return
- for i in xrange(len(self.__hierarchy_lines)-1,0,-1):
- line = self.__hierarchy_lines[i]
- if line.startswith("BaryonFileName") or \
- line.startswith("FileName "):
- testGrid = line.split("=")[-1].strip().rstrip()
- break
- if testGrid[0] != os.path.sep:
- testGrid = os.path.join(self.directory, testGrid)
- if not os.path.exists(testGrid):
- testGrid = os.path.join(self.directory,
- os.path.basename(testGrid))
- mylog.debug("Your data uses the annoying hardcoded path.")
- self._strip_path = True
- try:
- a = SD.SD(testGrid)
- self.data_style = 4
- mylog.debug("Detected HDF4")
- except:
- a = tables.openFile(testGrid, 'r')
- for b in a.iterNodes("/"):
- c = "%s" % (b)
- break
- if c.startswith("/Grid"):
- mylog.debug("Detected packed HDF5")
- self.data_style = 6
- else:
- mylog.debug("Detected unpacked HDF5")
- self.data_style = 5
- a.close()
- def __setup_filemap(self, grid):
- if not self.data_style == 6:
- return
- self.cpu_map[grid.filename].append(grid)
-
- def __setup_classes(self):
- dd = { 'readDataFast' : _data_style_funcs[self.data_style][0],
- 'readAllData' : _data_style_funcs[self.data_style][1],
- 'getFields' : _data_style_funcs[self.data_style][2],
- 'readDataSlice' : _data_style_funcs[self.data_style][3],
- '_read_data' : _data_style_funcs[self.data_style][0],
- '_read_all_data' : _data_style_funcs[self.data_style][1],
- '_read_field_names' : _data_style_funcs[self.data_style][2],
- '_read_data_slice' : _data_style_funcs[self.data_style][3],
- '_read_exception' : _data_style_funcs[self.data_style][4](),
- 'pf' : self.parameter_file, # Already weak
- 'hierarchy': weakref.proxy(self) }
- self.grid = classobj("EnzoGrid",(EnzoGridBase,), dd)
- self.proj = classobj("AMRProj",(AMRProjBase,), dd)
- self.slice = classobj("AMRSlice",(AMRSliceBase,), dd)
- self.region = classobj("AMRRegion",(AMRRegionBase,), dd)
- self.covering_grid = classobj("AMRCoveringGrid",(AMRCoveringGrid,), dd)
- self.smoothed_covering_grid = classobj("AMRSmoothedCoveringGrid",(AMRSmoothedCoveringGrid,), dd)
- self.sphere = classobj("AMRSphere",(AMRSphereBase,), dd)
- self.cutting = classobj("AMRCuttingPlane",(AMRCuttingPlaneBase,), dd)
- self.ray = classobj("AMROrthoRay",(AMROrthoRayBase,), dd)
- self.disk = classobj("AMRCylinder",(AMRCylinderBase,), dd)
-
- def __initialize_data_file(self):
+ def _initialize_data_file(self):
if not ytcfg.getboolean('lagos','serialize'): return
fn = os.path.join(self.directory,"%s.yt" % self["CurrentTimeIdentifier"])
if ytcfg.getboolean('lagos','onlydeserialize'):
@@ -193,16 +97,16 @@
else:
mode = 'a'
try:
- self.__data_file = tables.openFile(fn, mode)
+ self._data_file = tables.openFile(fn, mode)
my_name = self.get_data("/","MyName")
if my_name is None:
self.save_data(str(self.parameter_file), "/", "MyName")
else:
if str(my_name.read())!=str(self.parameter_file):
- self.__data_file.close()
- self.__data_file = None
+ self._data_file.close()
+ self._data_file = None
except:
- self.__data_file = None
+ self._data_file = None
pass
def save_data(self, array, node, name, set_attr=None, force=False):
@@ -211,36 +115,337 @@
described by *node* and *name*. If data file does not exist, it throws
no error and simply does not save.
"""
- if self.__data_file is None: return
+ if self._data_file is None: return
if force:
try:
- node_loc = self.__data_file.getNode(node)
+ node_loc = self._data_file.getNode(node)
if name in node_loc:
- self.__data_file.removeNode(node, name, recursive=True)
+ self._data_file.removeNode(node, name, recursive=True)
except tables.exceptions.NoSuchNodeError:
pass
- arr = self.__data_file.createArray(node, name, array, createparents=True)
+ arr = self._data_file.createArray(node, name, array, createparents=True)
if set_attr is not None:
for i, j in set_attr.items(): arr.setAttr(i,j)
- self.__data_file.flush()
+ self._data_file.flush()
def get_data(self, node, name):
"""
Return the dataset with a given *name* located at *node* in the
datafile.
"""
- if self.__data_file == None:
+ if self._data_file == None:
+ return None
+ try:
+ return self._data_file.getNode(node, name)
+ except tables.exceptions.NoSuchNodeError:
return None
+
+ def _close_data_file(self):
+ if self._data_file:
+ self._data_file.close()
+ del self._data_file
+ self._data_file = None
+
+ def _setup_classes(self, dd):
+ self.proj = classobj("AMRProj",(AMRProjBase,), dd)
+ self.slice = classobj("AMRSlice",(AMRSliceBase,), dd)
+ self.region = classobj("AMRRegion",(AMRRegionBase,), dd)
+ self.covering_grid = classobj("AMRCoveringGrid",(AMRCoveringGrid,), dd)
+ self.smoothed_covering_grid = classobj("AMRSmoothedCoveringGrid",(AMRSmoothedCoveringGrid,), dd)
+ self.sphere = classobj("AMRSphere",(AMRSphereBase,), dd)
+ self.cutting = classobj("AMRCuttingPlane",(AMRCuttingPlaneBase,), dd)
+ self.ray = classobj("AMROrthoRay",(AMROrthoRayBase,), dd)
+ self.disk = classobj("AMRCylinder",(AMRCylinderBase,), dd)
+
+ def __deserialize_hierarchy(self, harray):
+ mylog.debug("Cached entry found.")
+ self.gridDimensions[:] = harray[:,0:3]
+ self.gridStartIndices[:] = harray[:,3:6]
+ self.gridEndIndices[:] = harray[:,6:9]
+ self.gridLeftEdge[:] = harray[:,9:12]
+ self.gridRightEdge[:] = harray[:,12:15]
+ self.gridLevels[:] = harray[:,15:16]
+ self.gridTimes[:] = harray[:,16:17]
+ self.gridNumberOfParticles[:] = harray[:,17:18]
+
+ def _get_data_reader_dict(self):
+ dd = { 'readDataFast' : _data_style_funcs[self.data_style][0],
+ 'readAllData' : _data_style_funcs[self.data_style][1],
+ 'getFields' : _data_style_funcs[self.data_style][2],
+ 'readDataSlice' : _data_style_funcs[self.data_style][3],
+ '_read_data' : _data_style_funcs[self.data_style][0],
+ '_read_all_data' : _data_style_funcs[self.data_style][1],
+ '_read_field_names' : _data_style_funcs[self.data_style][2],
+ '_read_data_slice' : _data_style_funcs[self.data_style][3],
+ '_read_exception' : _data_style_funcs[self.data_style][4](),
+ 'pf' : self.parameter_file, # Already weak
+ 'hierarchy': weakref.proxy(self) }
+ return dd
+
+ def select_grids(self, level):
+ """
+ Returns an array of grids at *level*.
+ """
+ return self.grids[self._select_level(level)]
+
+ def get_smallest_dx(self):
+ """
+ Returns (in code units) the smallest cell size in the simulation.
+ """
+ return self.gridDxs.min()
+
+ def print_stats(self):
+ """
+ Prints out (stdout) relevant information about the simulation
+ """
+ for i in xrange(MAXLEVEL):
+ if (self.level_stats['numgrids'][i]) == 0:
+ break
+ print "% 3i\t% 6i\t% 11i" % \
+ (i, self.level_stats['numgrids'][i], self.level_stats['numcells'][i])
+ dx = self.gridDxs[self.levelIndices[i][0]]
+ print "-" * 28
+ print " \t% 6i\t% 11i" % (self.level_stats['numgrids'].sum(), self.level_stats['numcells'].sum())
+ print "\n"
+ try:
+ print "z = %0.8f" % (self["CosmologyCurrentRedshift"])
+ except:
+ pass
+ t_s = self["InitialTime"] * self["Time"]
+ print "t = %0.8e = %0.8e s = %0.8e years" % \
+ (self["InitialTime"], \
+ t_s, t_s / (365*24*3600.0) )
+ print "\nSmallest Cell:"
+ u=[]
+ for item in self.parameter_file.units.items():
+ u.append((item[1],item[0]))
+ u.sort()
+ for unit in u:
+ print "\tWidth: %0.3e %s" % (dx*unit[0], unit[1])
+
+ def find_point(self, coord):
+ """
+ Returns the (objects, indices) of grids containing an (x,y,z) point
+ """
+ mask=na.ones(self.num_grids)
+ for i in xrange(len(coord)):
+ na.choose(na.greater(self.gridLeftEdge[:,i],coord[i]), (mask,0), mask)
+ na.choose(na.greater(self.gridRightEdge[:,i],coord[i]), (0,mask), mask)
+ ind = na.where(mask == 1)
+ return self.grids[ind], ind
+
+ def find_ray_grids(self, coord, axis):
+ """
+ Returns the (objects, indices) of grids that an (x,y) ray intersects
+ along *axis*
+ """
+ # Let's figure out which grids are on the slice
+ mask=na.ones(self.num_grids)
+ # So if gRE > coord, we get a mask, if not, we get a zero
+ # if gLE > coord, we get a zero, if not, mask
+ # Thus, if the coordinate is between the two edges, we win!
+ na.choose(na.greater(self.gridRightEdge[:,x_dict[axis]],coord[0]),(0,mask),mask)
+ na.choose(na.greater(self.gridLeftEdge[:,x_dict[axis]],coord[0]),(mask,0),mask)
+ na.choose(na.greater(self.gridRightEdge[:,y_dict[axis]],coord[1]),(0,mask),mask)
+ na.choose(na.greater(self.gridLeftEdge[:,y_dict[axis]],coord[1]),(mask,0),mask)
+ ind = na.where(mask == 1)
+ return self.grids[ind], ind
+
+ def find_slice_grids(self, coord, axis):
+ """
+ Returns the (objects, indices) of grids that a slice intersects along
+ *axis*
+ """
+ # Let's figure out which grids are on the slice
+ mask=na.ones(self.num_grids)
+ # So if gRE > coord, we get a mask, if not, we get a zero
+ # if gLE > coord, we get a zero, if not, mask
+ # Thus, if the coordinate is between the edges, we win!
+ #ind = na.where( na.logical_and(self.gridRightEdge[:,axis] > coord, \
+ #self.gridLeftEdge[:,axis] < coord))
+ na.choose(na.greater(self.gridRightEdge[:,axis],coord),(0,mask),mask)
+ na.choose(na.greater(self.gridLeftEdge[:,axis],coord),(mask,0),mask)
+ ind = na.where(mask == 1)
+ return self.grids[ind], ind
+
+ def find_sphere_grids(self, center, radius):
+ """
+ Returns objects, indices of grids within a sphere
+ """
+ centers = (self.gridRightEdge + self.gridLeftEdge)/2.0
+ long_axis = na.maximum.reduce(self.gridRightEdge - self.gridLeftEdge, 1)
+ t = centers - center
+ dist = na.sqrt(t[:,0]**2+t[:,1]**2+t[:,2]**2)
+ gridI = na.where(na.logical_and((self.gridDxs<=radius)[:,0],(dist < (radius + long_axis))) == 1)
+ return self.grids[gridI], gridI
+
+ def get_box_grids(self, leftEdge, rightEdge):
+ """
+ Gets back all the grids between a left edge and right edge
+ """
+ gridI = na.where((na.all(self.gridRightEdge > leftEdge, axis=1)
+ & na.all(self.gridLeftEdge < rightEdge, axis=1)) == True)
+ return self.grids[gridI], gridI
+
+ @time_execution
+ def find_max(self, field, finestLevels = True):
+ """
+ Returns (value, center) of location of maximum for a given field.
+ """
+ if finestLevels:
+ gI = na.where(self.gridLevels >= self.maxLevel - NUMTOCHECK)
+ else:
+ gI = na.where(self.gridLevels >= 0) # Slow but pedantic
+ maxVal = -1e100
+ for grid in self.grids[gI[0]]:
+ mylog.debug("Checking %s (level %s)", grid.id, grid.Level)
+ val, coord = grid.find_max(field)
+ if val > maxVal:
+ maxCoord = coord
+ maxVal = val
+ maxGrid = grid
+ mc = na.array(maxCoord)
+ pos=maxGrid.get_position(mc)
+ pos[0] += 0.5*maxGrid.dx
+ pos[1] += 0.5*maxGrid.dx
+ pos[2] += 0.5*maxGrid.dx
+ mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s %s", \
+ maxVal, pos[0], pos[1], pos[2], maxGrid, maxGrid.Level, mc)
+ self.center = pos
+ # This probably won't work for anyone else
+ self.bulkVelocity = (maxGrid["x-velocity"][maxCoord], \
+ maxGrid["y-velocity"][maxCoord], \
+ maxGrid["z-velocity"][maxCoord])
+ self.parameters["Max%sValue" % (field)] = maxVal
+ self.parameters["Max%sPos" % (field)] = "%s" % (pos)
+ return maxVal, pos
+
+ findMax = find_max
+
+ @time_execution
+ def find_min(self, field):
+ """
+ Returns (value, center) of location of minimum for a given field
+ """
+ gI = na.where(self.gridLevels >= 0) # Slow but pedantic
+ minVal = 1e100
+ for grid in self.grids[gI[0]]:
+ mylog.debug("Checking %s (level %s)", grid.id, grid.Level)
+ val, coord = grid.find_min(field)
+ if val < minVal:
+ minCoord = coord
+ minVal = val
+ minGrid = grid
+ mc = na.array(minCoord)
+ pos=minGrid.get_position(mc)
+ pos[0] += 0.5*minGrid.dx
+ pos[1] += 0.5*minGrid.dx
+ pos[2] += 0.5*minGrid.dx
+ mylog.info("Min Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s", \
+ minVal, pos[0], pos[1], pos[2], minGrid, minGrid.Level)
+ self.center = pos
+ # This probably won't work for anyone else
+ self.binkVelocity = (minGrid["x-velocity"][minCoord], \
+ minGrid["y-velocity"][minCoord], \
+ minGrid["z-velocity"][minCoord])
+ self.parameters["Min%sValue" % (field)] = minVal
+ self.parameters["Min%sPos" % (field)] = "%s" % (pos)
+ return minVal, pos
+
+ def _get_parameters(self):
+ return self.parameter_file.parameters
+ parameters=property(_get_parameters)
+
+ def __getitem__(self, item):
+ return self.parameter_file[item]
+
+
+class EnzoHierarchy(AMRHierarchy):
+ eiTopGrid = None
+ _strip_path = False
+ @time_execution
+ def __init__(self, pf, data_style=None):
+ """
+ This is the grid structure as Enzo sees it, with some added bonuses.
+ It's primarily used as a class factory, to generate data objects and
+ access grids.
+
+ It should never be created directly -- you should always access it via
+ calls to an affiliated :class:`~yt.lagos.EnzoStaticOutput`.
+
+ On instantiation, it processes the hierarchy and generates the grids.
+ """
+ # Expect filename to be the name of the parameter file, not the
+ # hierarchy
+ self.data_style = data_style
+ self.hierarchy_filename = os.path.abspath(pf.parameter_filename) \
+ + ".hierarchy"
+ self.__hierarchy_lines = open(self.hierarchy_filename).readlines()
+ if len(self.__hierarchy_lines) == 0:
+ raise IOError(-1,"File empty", self.hierarchy_filename)
+ self.boundary_filename = os.path.abspath(pf.parameter_filename) \
+ + ".boundary"
+ self.directory = os.path.dirname(self.hierarchy_filename)
+ # Now we search backwards from the end of the file to find out how many
+ # grids we have, which allows us to preallocate memory
+ self.__hierarchy_string = open(self.hierarchy_filename).read()
+ for line in reversed(self.__hierarchy_lines):
+ if line.startswith("Grid ="):
+ self.num_grids = int(line.split("=")[-1])
+ break
+ self.__guess_data_style()
+ # For some reason, r8 seems to want Float64
+ if pf.has_key("CompilerPrecision") \
+ and pf["CompilerPrecision"] == "r4":
+ self.float_type = 'float32'
+ else:
+ self.float_type = 'float64'
+
+ AMRHierarchy.__init__(self, pf)
+
+ del self.__hierarchy_string, self.__hierarchy_lines
+
+ def _setup_classes(self):
+ dd = self._get_data_reader_dict()
+ self.grid = classobj("EnzoGrid",(EnzoGridBase,), dd)
+ AMRHierarchy._setup_classes(self, dd)
+
+ def __guess_data_style(self):
+ if self.data_style: return
+ for i in xrange(len(self.__hierarchy_lines)-1,0,-1):
+ line = self.__hierarchy_lines[i]
+ if line.startswith("BaryonFileName") or \
+ line.startswith("FileName "):
+ testGrid = line.split("=")[-1].strip().rstrip()
+ break
+ if testGrid[0] != os.path.sep:
+ testGrid = os.path.join(self.directory, testGrid)
+ if not os.path.exists(testGrid):
+ testGrid = os.path.join(self.directory,
+ os.path.basename(testGrid))
+ mylog.debug("Your data uses the annoying hardcoded path.")
+ self._strip_path = True
try:
- return self.__data_file.getNode(node, name)
- except tables.exceptions.NoSuchNodeError:
- return None
+ a = SD.SD(testGrid)
+ self.data_style = 4
+ mylog.debug("Detected HDF4")
+ except:
+ a = tables.openFile(testGrid, 'r')
+ for b in a.iterNodes("/"):
+ c = "%s" % (b)
+ break
+ if c.startswith("/Grid"):
+ mylog.debug("Detected packed HDF5")
+ self.data_style = 6
+ else:
+ mylog.debug("Detected unpacked HDF5")
+ self.data_style = 5
+ a.close()
- def _close_data_file(self):
- if self.__data_file:
- self.__data_file.close()
- del self.__data_file
- self.__data_file = None
+ def __setup_filemap(self, grid):
+ if not self.data_style == 6:
+ return
+ self.cpu_map[grid.filename].append(grid)
def __del__(self):
self._close_data_file()
@@ -254,34 +459,9 @@
del self.gridReverseTree
del self.gridLeftEdge, self.gridRightEdge
del self.gridLevels, self.gridStartIndices, self.gridEndIndices
- del self.gridTimes, self.__hierarchy_string, self.__hierarchy_lines
+ del self.gridTimes
del self.gridTree
- def __deserialize_hierarchy(self, harray):
- mylog.debug("Cached entry found.")
- self.gridDimensions[:] = harray[:,0:3]
- self.gridStartIndices[:] = harray[:,3:6]
- self.gridEndIndices[:] = harray[:,6:9]
- self.gridLeftEdge[:] = harray[:,9:12]
- self.gridRightEdge[:] = harray[:,12:15]
- self.gridLevels[:] = harray[:,15:16]
- self.gridTimes[:] = harray[:,16:17]
- self.gridNumberOfParticles[:] = harray[:,17:18]
- del harray
- mylog.debug("Copied to local array.")
- # Now get the baryon filenames
- mylog.debug("Getting baryon filenames")
- re_BaryonFileName = constructRegularExpressions("BaryonFileName",('s'))
- fn_results = re.findall(re_BaryonFileName, self.__hierarchy_string)
- if len(fn_results):
- self.__set_all_filenames(fn_results)
- return
- re_BaryonFileName = constructRegularExpressions("FileName",('s'))
- fn_results = re.findall(re_BaryonFileName, self.__hierarchy_string)
- self.__set_all_filenames(fn_results)
- # This is pretty bad, but we do it to save a significant amount of time
- # in larger runs.
-
def __set_all_filenames(self, fns):
if self._strip_path:
for fnI, fn in enumerate(fns):
@@ -358,6 +538,22 @@
self.save_data(allArrays, "/","Hierarchy")
del allArrays
+ def __obtain_filenames(self):
+ mylog.debug("Copied to local array.")
+ # This needs to go elsewhere:
+ # Now get the baryon filenames
+ mylog.debug("Getting baryon filenames")
+ re_BaryonFileName = constructRegularExpressions("BaryonFileName",('s'))
+ fn_results = re.findall(re_BaryonFileName, self.__hierarchy_string)
+ if len(fn_results):
+ self.__set_all_filenames(fn_results)
+ return
+ re_BaryonFileName = constructRegularExpressions("FileName",('s'))
+ fn_results = re.findall(re_BaryonFileName, self.__hierarchy_string)
+ self.__set_all_filenames(fn_results)
+ # This is pretty bad, but we do it to save a significant amount of time
+ # in larger runs.
+
def __setup_grid_tree(self):
mylog.debug("No cached tree found, creating")
self.grids[0].Level = 0 # Bootstrap
@@ -389,7 +585,7 @@
self.save_data(self.gridLevels, "/", "Levels")
@time_execution
- def __populate_hierarchy(self):
+ def _populate_hierarchy(self):
"""
Instantiates all of the grid objects, with their appropriate
parameters. This is the work-horse.
@@ -404,6 +600,7 @@
self.__deserialize_hierarchy(harray)
else:
self.__parse_hierarchy_file()
+ self.__obtain_filenames()
treeArray = self.get_data("/", "Tree")
if treeArray == None:
self.__setup_grid_tree()
@@ -439,7 +636,7 @@
self.level_stats[level]['numgrids'] = na.where(self.gridLevels==level)[0].size
li = na.where(self.gridLevels[:,0] == level)
self.level_stats[level]['numcells'] = ad[li].prod(axis=1).sum()
- self.levelIndices[level] = self.__select_level(level)
+ self.levelIndices[level] = self._select_level(level)
self.levelNum[level] = len(self.levelIndices[level])
mylog.debug("Hierarchy fully populated.")
@@ -516,181 +713,11 @@
- def __select_level(self, level):
+ def _select_level(self, level):
# We return a numarray of the indices of all the grids on a given level
indices = na.where(self.gridLevels[:,0] == level)[0]
return indices
- def select_grids(self, level):
- """
- Returns an array of grids at *level*.
- """
- return self.grids[self.__select_level(level)]
-
- def get_smallest_dx(self):
- """
- Returns (in code units) the smallest cell size in the simulation.
- """
- return self.gridDxs.min()
-
- def print_stats(self):
- """
- Prints out (stdout) relevant information about the simulation
- """
- for i in xrange(MAXLEVEL):
- if (self.level_stats['numgrids'][i]) == 0:
- break
- print "% 3i\t% 6i\t% 11i" % \
- (i, self.level_stats['numgrids'][i], self.level_stats['numcells'][i])
- dx = self.gridDxs[self.levelIndices[i][0]]
- print "-" * 28
- print " \t% 6i\t% 11i" % (self.level_stats['numgrids'].sum(), self.level_stats['numcells'].sum())
- print "\n"
- try:
- print "z = %0.8f" % (self["CosmologyCurrentRedshift"])
- except:
- pass
- t_s = self["InitialTime"] * self["Time"]
- print "t = %0.8e = %0.8e s = %0.8e years" % \
- (self["InitialTime"], \
- t_s, t_s / (365*24*3600.0) )
- print "\nSmallest Cell:"
- u=[]
- for item in self.parameter_file.units.items():
- u.append((item[1],item[0]))
- u.sort()
- for unit in u:
- print "\tWidth: %0.3e %s" % (dx*unit[0], unit[1])
-
- def find_point(self, coord):
- """
- Returns the (objects, indices) of grids containing an (x,y,z) point
- """
- mask=na.ones(self.num_grids)
- for i in xrange(len(coord)):
- na.choose(na.greater(self.gridLeftEdge[:,i],coord[i]), (mask,0), mask)
- na.choose(na.greater(self.gridRightEdge[:,i],coord[i]), (0,mask), mask)
- ind = na.where(mask == 1)
- return self.grids[ind], ind
-
- def find_ray_grids(self, coord, axis):
- """
- Returns the (objects, indices) of grids that an (x,y) ray intersects
- along *axis*
- """
- # Let's figure out which grids are on the slice
- mask=na.ones(self.num_grids)
- # So if gRE > coord, we get a mask, if not, we get a zero
- # if gLE > coord, we get a zero, if not, mask
- # Thus, if the coordinate is between the two edges, we win!
- na.choose(na.greater(self.gridRightEdge[:,x_dict[axis]],coord[0]),(0,mask),mask)
- na.choose(na.greater(self.gridLeftEdge[:,x_dict[axis]],coord[0]),(mask,0),mask)
- na.choose(na.greater(self.gridRightEdge[:,y_dict[axis]],coord[1]),(0,mask),mask)
- na.choose(na.greater(self.gridLeftEdge[:,y_dict[axis]],coord[1]),(mask,0),mask)
- ind = na.where(mask == 1)
- return self.grids[ind], ind
-
- def find_slice_grids(self, coord, axis):
- """
- Returns the (objects, indices) of grids that a slice intersects along
- *axis*
- """
- # Let's figure out which grids are on the slice
- mask=na.ones(self.num_grids)
- # So if gRE > coord, we get a mask, if not, we get a zero
- # if gLE > coord, we get a zero, if not, mask
- # Thus, if the coordinate is between the edges, we win!
- #ind = na.where( na.logical_and(self.gridRightEdge[:,axis] > coord, \
- #self.gridLeftEdge[:,axis] < coord))
- na.choose(na.greater(self.gridRightEdge[:,axis],coord),(0,mask),mask)
- na.choose(na.greater(self.gridLeftEdge[:,axis],coord),(mask,0),mask)
- ind = na.where(mask == 1)
- return self.grids[ind], ind
-
- def find_sphere_grids(self, center, radius):
- """
- Returns objects, indices of grids within a sphere
- """
- centers = (self.gridRightEdge + self.gridLeftEdge)/2.0
- long_axis = na.maximum.reduce(self.gridRightEdge - self.gridLeftEdge, 1)
- t = centers - center
- dist = na.sqrt(t[:,0]**2+t[:,1]**2+t[:,2]**2)
- gridI = na.where(na.logical_and((self.gridDxs<=radius)[:,0],(dist < (radius + long_axis))) == 1)
- return self.grids[gridI], gridI
-
- def get_box_grids(self, leftEdge, rightEdge):
- """
- Gets back all the grids between a left edge and right edge
- """
- gridI = na.where((na.all(self.gridRightEdge > leftEdge, axis=1)
- & na.all(self.gridLeftEdge < rightEdge, axis=1)) == True)
- return self.grids[gridI], gridI
-
- @time_execution
- def find_max(self, field, finestLevels = True):
- """
- Returns (value, center) of location of maximum for a given field.
- """
- if finestLevels:
- gI = na.where(self.gridLevels >= self.maxLevel - NUMTOCHECK)
- else:
- gI = na.where(self.gridLevels >= 0) # Slow but pedantic
- maxVal = -1e100
- for grid in self.grids[gI[0]]:
- mylog.debug("Checking %s (level %s)", grid.id, grid.Level)
- val, coord = grid.find_max(field)
- if val > maxVal:
- maxCoord = coord
- maxVal = val
- maxGrid = grid
- mc = na.array(maxCoord)
- pos=maxGrid.get_position(mc)
- pos[0] += 0.5*maxGrid.dx
- pos[1] += 0.5*maxGrid.dx
- pos[2] += 0.5*maxGrid.dx
- mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s %s", \
- maxVal, pos[0], pos[1], pos[2], maxGrid, maxGrid.Level, mc)
- self.center = pos
- # This probably won't work for anyone else
- self.bulkVelocity = (maxGrid["x-velocity"][maxCoord], \
- maxGrid["y-velocity"][maxCoord], \
- maxGrid["z-velocity"][maxCoord])
- self.parameters["Max%sValue" % (field)] = maxVal
- self.parameters["Max%sPos" % (field)] = "%s" % (pos)
- return maxVal, pos
-
- findMax = find_max
-
- @time_execution
- def find_min(self, field):
- """
- Returns (value, center) of location of minimum for a given field
- """
- gI = na.where(self.gridLevels >= 0) # Slow but pedantic
- minVal = 1e100
- for grid in self.grids[gI[0]]:
- mylog.debug("Checking %s (level %s)", grid.id, grid.Level)
- val, coord = grid.find_min(field)
- if val < minVal:
- minCoord = coord
- minVal = val
- minGrid = grid
- mc = na.array(minCoord)
- pos=minGrid.get_position(mc)
- pos[0] += 0.5*minGrid.dx
- pos[1] += 0.5*minGrid.dx
- pos[2] += 0.5*minGrid.dx
- mylog.info("Min Value is %0.5e at %0.16f %0.16f %0.16f in grid %s at level %s", \
- minVal, pos[0], pos[1], pos[2], minGrid, minGrid.Level)
- self.center = pos
- # This probably won't work for anyone else
- self.binkVelocity = (minGrid["x-velocity"][minCoord], \
- minGrid["y-velocity"][minCoord], \
- minGrid["z-velocity"][minCoord])
- self.parameters["Min%sValue" % (field)] = minVal
- self.parameters["Min%sPos" % (field)] = "%s" % (pos)
- return minVal, pos
-
@time_execution
def export_particles_pb(self, filename, filter = 1, indexboundary = 0, fields = None, scale=1.0):
"""
@@ -740,7 +767,7 @@
f=open(filename,"w")
for l in xrange(self.maxLevel):
f.write("add object g%s = l%s\n" % (l,l))
- ind = self.__select_level(l)
+ ind = self._select_level(l)
for i in ind:
f.write("add box -n %s -l %s %s,%s %s,%s %s,%s\n" % \
(i+1, self.gridLevels.ravel()[i],
@@ -749,14 +776,6 @@
self.gridLeftEdge[i,2], self.gridRightEdge[i,2]))
- """@todo: Fix these!"""
- def _get_parameters(self):
- return self.parameter_file.parameters
- parameters=property(_get_parameters)
-
- def __getitem__(self, item):
- return self.parameter_file[item]
-
scanf_regex = {}
scanf_regex['e'] = r"[-+]?\d+\.?\d*?|\.\d+[eE][-+]?\d+?"
scanf_regex['g'] = scanf_regex['e']
More information about the yt-svn
mailing list