[Yt-svn] yt-commit r950 - in trunk: tests yt yt/extensions yt/lagos yt/raven yt/reason
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Wed Nov 19 09:56:40 PST 2008
Author: mturk
Date: Wed Nov 19 09:56:37 2008
New Revision: 950
URL: http://yt.spacepope.org/changeset/950
Log:
Merging 2D and 1D back to trunk.
This includes the new yt-generalization FieldInfoContainer. DerivedFields will
disappear in the next commit.
Added:
trunk/yt/lagos/EnzoFields.py
- copied unchanged from r949, /branches/yt-non-3d/yt/lagos/EnzoFields.py
trunk/yt/lagos/FieldInfoContainer.py
- copied unchanged from r949, /branches/yt-non-3d/yt/lagos/FieldInfoContainer.py
trunk/yt/lagos/UniversalFields.py
- copied, changed from r949, /branches/yt-non-3d/yt/lagos/UniversalFields.py
Modified:
trunk/tests/test_lagos.py
trunk/yt/extensions/SpectralIntegrator.py
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/BaseGridType.py
trunk/yt/lagos/DataReadingFuncs.py
trunk/yt/lagos/DerivedFields.py
trunk/yt/lagos/HierarchyType.py
trunk/yt/lagos/OutputTypes.py
trunk/yt/lagos/ParallelTools.py
trunk/yt/lagos/Profiles.py
trunk/yt/lagos/__init__.py
trunk/yt/mods.py
trunk/yt/raven/Callbacks.py
trunk/yt/raven/PlotCollection.py
trunk/yt/raven/PlotTypes.py
trunk/yt/raven/_MPL.c
trunk/yt/raven/__init__.py
trunk/yt/reason/Functions.py
Modified: trunk/tests/test_lagos.py
==============================================================================
--- trunk/tests/test_lagos.py (original)
+++ trunk/tests/test_lagos.py Wed Nov 19 09:56:37 2008
@@ -2,7 +2,7 @@
Test that we can get outputs, and interact with them in some primitive ways.
"""
-# @TODO: Add unit test for deleting field from fieldInfo
+# @TODO: Add unit test for deleting field from FieldInfo
import unittest, glob, os.path, os, sys, StringIO
@@ -190,7 +190,8 @@
self.assertEqual(len(cid), 3)
-for field in yt.lagos.fieldInfo.values():
+for field_name in yt.lagos.FieldInfo:
+ field = yt.lagos.FieldInfo[field_name]
setattr(DataTypeTestingBase, "test%s" % field.name, _returnFieldFunction(field))
field = "Temperature"
Modified: trunk/yt/extensions/SpectralIntegrator.py
==============================================================================
--- trunk/yt/extensions/SpectralIntegrator.py (original)
+++ trunk/yt/extensions/SpectralIntegrator.py Wed Nov 19 09:56:37 2008
@@ -57,7 +57,7 @@
def add_frequency_bin_field(self, ev_min, ev_max):
"""
- Add a new field to the global fieldInfo dictionary, which is an
+ Add a new field to the FieldInfoContainer, which is an
integrated bin from *ev_min* to *ev_max*.
Returns the name of the new field.
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Wed Nov 19 09:56:37 2008
@@ -318,16 +318,16 @@
temp = grid[field]
def _generate_field(self, field):
- if fieldInfo.has_key(field):
+ if self.pf.field_info.has_key(field):
# First we check the validator
try:
- fieldInfo[field].check_available(self)
+ self.pf.field_info[field].check_available(self)
except NeedsGridType, ngt_exception:
# We leave this to be implementation-specific
self._generate_field_in_grids(field, ngt_exception.ghost_zones)
return False
else:
- self[field] = fieldInfo[field](self)
+ self[field] = self.pf.field_info[field](self)
return True
else: # Can't find the field, try as it might
raise exceptions.KeyError(field)
@@ -508,16 +508,16 @@
def _generate_field(self, field):
- if fieldInfo.has_key(field):
+ if self.pf.field_info.has_key(field):
# First we check the validator
try:
- fieldInfo[field].check_available(self)
+ self.pf.field_info[field].check_available(self)
except NeedsGridType, ngt_exception:
# We leave this to be implementation-specific
self._generate_field_in_grids(field, ngt_exception.ghost_zones)
return False
else:
- self[field] = fieldInfo[field](self)
+ self[field] = self.pf.field_info[field](self)
return True
else: # Can't find the field, try as it might
raise exceptions.KeyError(field)
@@ -638,16 +638,13 @@
self['py'] = t[:,1]
self['pz'] = t[:,2]
self['pdx'] = t[:,3]
- self['pdy'] = t[:,3]
- self['pdz'] = t[:,3]
+ self['pdy'] = t[:,4]
+ self['pdz'] = t[:,3] # Does not matter!
# Now we set the *actual* coordinates
self[axis_names[x_dict[self.axis]]] = t[:,0]
self[axis_names[y_dict[self.axis]]] = t[:,1]
self[axis_names[self.axis]] = t[:,2]
- self['dx'] = t[:,3]
- self['dy'] = t[:,3]
- self['dz'] = t[:,3]
self.ActiveDimensions = (t.shape[0], 1, 1)
@@ -657,7 +654,8 @@
def _generate_grid_coords(self, grid):
xaxis = x_dict[self.axis]
yaxis = y_dict[self.axis]
- wantedIndex = int(((self.coord-grid.LeftEdge[self.axis])/grid.dx))
+ ds, dx, dy = grid['dds'][self.axis], grid['dds'][xaxis], grid['dds'][yaxis]
+ wantedIndex = int(((self.coord-grid.LeftEdge[self.axis])/ds))
sl = [slice(None), slice(None), slice(None)]
sl[self.axis] = slice(wantedIndex, wantedIndex + 1)
#sl.reverse()
@@ -668,31 +666,33 @@
cmI = na.indices((nx,ny))
xind = cmI[0,:].ravel()
xpoints = na.ones(cm[0].shape, 'float64')
- xpoints *= xind[cm]*grid.dx+(grid.LeftEdge[xaxis] + 0.5*grid.dx)
+ xpoints *= xind[cm]*dx+(grid.LeftEdge[xaxis] + 0.5*dx)
yind = cmI[1,:].ravel()
ypoints = na.ones(cm[0].shape, 'float64')
- ypoints *= yind[cm]*grid.dx+(grid.LeftEdge[yaxis] + 0.5*grid.dx)
+ ypoints *= yind[cm]*dy+(grid.LeftEdge[yaxis] + 0.5*dy)
zpoints = na.ones(xpoints.shape, 'float64') * self.coord
- dx = na.ones(xpoints.shape, 'float64') * grid.dx/2.0
- t = na.array([xpoints, ypoints, zpoints, dx]).swapaxes(0,1)
+ dx = na.ones(xpoints.shape, 'float64') * dx/2.0
+ dy = na.ones(xpoints.shape, 'float64') * dy/2.0
+ t = na.array([xpoints, ypoints, zpoints, dx, dy]).swapaxes(0,1)
return t
@restore_grid_state
def _get_data_from_grid(self, grid, field):
# So what's our index of slicing? This is what we need to figure out
# first, so we can deal with our data in the fastest way.
- wantedIndex = int(((self.coord-grid.LeftEdge[self.axis])/grid.dx))
+ dx = grid['dds'][self.axis]
+ wantedIndex = int(((self.coord-grid.LeftEdge[self.axis])/dx))
sl = [slice(None), slice(None), slice(None)]
sl[self.axis] = slice(wantedIndex, wantedIndex + 1)
sl = tuple(sl)
- if fieldInfo.has_key(field) and fieldInfo[field].particle_type:
+ if self.pf.field_info.has_key(field) and self.pf.field_info[field].particle_type:
return grid[field]
- elif field in fieldInfo and fieldInfo[field].not_in_all:
+ elif field in self.pf.field_info and self.pf.field_info[field].not_in_all:
dv = grid[field][sl]
elif not grid.has_key(field):
conv_factor = 1.0
- if fieldInfo.has_key(field):
- conv_factor = fieldInfo[field]._convert_function(self)
+ if self.pf.field_info.has_key(field):
+ conv_factor = self.pf.field_info[field]._convert_function(self)
dv = self._read_data_slice(grid, field, self.axis, wantedIndex) * conv_factor
else:
dv = grid[field]
@@ -782,9 +782,7 @@
D += (x * self._norm_vec[0]).reshape(ss[0],1,1)
D += (y * self._norm_vec[1]).reshape(1,ss[1],1)
D += (z * self._norm_vec[2]).reshape(1,1,ss[2])
- diag_dist = na.sqrt(grid.dx**2.0
- + grid.dy**2.0
- + grid.dz**2.0)
+ diag_dist = na.sqrt(na.sum(grid['dds']**2.0))
cm = (na.abs(D) <= 0.5*diag_dist) # Boolean
return cm
@@ -808,7 +806,7 @@
return na.array(coords).swapaxes(0,1)
def _get_data_from_grid(self, grid, field):
- if not fieldInfo[field].particle_type:
+ if not self.pf.field_info[field].particle_type:
pointI = self._get_point_indices(grid)
if grid[field].size == 1: # dx, dy, dz, cellvolume
t = grid[field] * na.ones(grid.ActiveDimensions)
@@ -925,7 +923,7 @@
for field in fields + [self._weight]:
if field is None: continue
dls.append(just_one(grid['d%s' % axis_names[self.axis]]))
- convs.append(self.pf.units[fieldInfo[field].projection_conversion])
+ convs.append(self.pf.units[self.pf.field_info[field].projection_conversion])
return na.array(dls), na.array(convs)
def __project_level(self, level, fields):
@@ -963,25 +961,9 @@
field_data *= convs[...,na.newaxis]
mylog.info("Level %s done: %s final", \
level, coord_data.shape[1])
- dx = grids_to_project[0].dx# * na.ones(coord_data.shape[0], dtype='float64')
+ dx = grids_to_project[0]['dds'][self.axis] # this is our dl
return coord_data, dx, field_data
- def __cleanup_level(self, level):
- pass
- grids_to_project = self.source.select_grids(level)
- coord_data = []
- field_data = []
- for grid in grids_to_project:
- if self.__retval_coords[grid.id][0].size == 0: continue
- if self._weight is not None:
- weightedData = grid.coarseData[2] / grid.coarseData[4]
- else:
- weightedData = grid.coarseData[2]
- all_data.append([grid.coarseData[0], grid.coarseData[1],
- weightedData,
- na.ones(grid.coarseData[0].shape,dtype='float64')*grid.dx])
- return na.concatenate(all_data, axis=1)
-
def __combine_grids_on_level(self, level):
grids = self.source.select_grids(level)
grids_i = self.source.levelIndices[level]
@@ -1022,7 +1004,8 @@
args = []
args += self.__retval_coords[grid2.id] + [self.__retval_fields[grid2.id]]
args += self.__retval_coords[grid1.id] + [self.__retval_fields[grid1.id]]
- args.append(int(grid2.dx / grid1.dx))
+ # Refinement factor, which is same in all directions
+ args.append(int(grid2.dx / grid1.dx))
args.append(na.ones(args[0].shape, dtype='int64'))
kk = PointCombine.CombineGrids(*args)
goodI = args[-1].astype('bool')
@@ -1040,7 +1023,7 @@
#@time_execution
def get_data(self, fields = None):
- if fields is None: fields = self.fields
+ if fields is None: fields = ensure_list(self.fields)[:]
fields = ensure_list(fields)
coord_data = []
field_data = []
@@ -1219,11 +1202,11 @@
@restore_grid_state
def _get_data_from_grid(self, grid, field):
- if field in fieldInfo and fieldInfo[field].particle_type:
+ if field in self.pf.field_info and self.pf.field_info[field].particle_type:
if grid.NumberOfParticles == 0: return na.array([])
pointI = self._get_particle_indices(grid)
return grid[field][pointI].ravel()
- if field in fieldInfo and fieldInfo[field].vector_field:
+ if field in self.pf.field_info and self.pf.field_info[field].vector_field:
pointI = self._get_point_indices(grid)
f = grid[field]
return na.array([f[i,:][pointI] for i in range(3)])
@@ -1251,16 +1234,16 @@
i += np
def _generate_field(self, field):
- if fieldInfo.has_key(field):
+ if self.pf.field_info.has_key(field):
# First we check the validator
try:
- fieldInfo[field].check_available(self)
+ self.pf.field_info[field].check_available(self)
except NeedsGridType, ngt_exception:
# We leave this to be implementation-specific
self._generate_field_in_grids(field, ngt_exception.ghost_zones)
return False
else:
- self[field] = fieldInfo[field](self)
+ self[field] = self.pf.field_info[field](self)
return True
else: # Can't find the field, try as it might
raise exceptions.KeyError(field)
@@ -1777,7 +1760,7 @@
def _get_level_array(self, level, fields):
fields = ensure_list(fields)
# We assume refinement by a factor of two
- rf = float(2**(self.level - level))
+ rf = float(self.pf["RefineBy"]**(self.level - level))
dims = na.maximum(1,self.ActiveDimensions/rf) + 2
dx = (self.right_edge-self.left_edge)/(dims-2)
x,y,z = (na.mgrid[0:dims[0],0:dims[1],0:dims[2]].astype('float64')+0.5)\
@@ -1793,7 +1776,7 @@
for ax in 'xyz': self['cd%s'%ax] = fake_grid['d%s'%ax]
for field in fields:
# Generate the new grid field
- if field in fieldInfo and fieldInfo[field].take_log:
+ if field in self.pf.field_info and self.pf.field_info[field].take_log:
interpolator = TrilinearFieldInterpolator(
na.log10(self[field]), bounds, ['x','y','z'],
truncate = True)
Modified: trunk/yt/lagos/BaseGridType.py
==============================================================================
--- trunk/yt/lagos/BaseGridType.py (original)
+++ trunk/yt/lagos/BaseGridType.py Wed Nov 19 09:56:37 2008
@@ -44,7 +44,26 @@
if filename: self.set_filename(filename)
self.overlap_masks = [None, None, None]
self._overlap_grids = [None, None, None]
+ self._file_access_pooling = False
self.pf = self.hierarchy.parameter_file # weakref already
+
+ def _generate_field(self, field):
+ if self.pf.field_info.has_key(field):
+ # First we check the validator
+ try:
+ self.pf.field_info[field].check_available(self)
+ except NeedsGridType, ngt_exception:
+ # This is only going to be raised if n_gz > 0
+ n_gz = ngt_exception.ghost_zones
+ f_gz = ngt_exception.fields
+ gz_grid = self.retrieve_ghost_zones(n_gz, f_gz)
+ temp_array = self.pf.field_info[field](gz_grid)
+ sl = [slice(n_gz,-n_gz)] * 3
+ self[field] = temp_array[sl]
+ else:
+ self[field] = self.pf.field_info[field](self)
+ else: # Can't find the field, try as it might
+ raise exceptions.KeyError, field
def get_data(self, field):
"""
@@ -53,8 +72,8 @@
if not self.data.has_key(field):
if field in self.hierarchy.field_list:
conv_factor = 1.0
- if fieldInfo.has_key(field):
- conv_factor = fieldInfo[field]._convert_function(self)
+ if self.pf.field_info.has_key(field):
+ conv_factor = self.pf.field_info[field]._convert_function(self)
try:
if hasattr(self.hierarchy, 'queue'):
temp = self.hierarchy.queue.pop(self, field)
@@ -62,34 +81,16 @@
temp = self.readDataFast(field)
self[field] = temp * conv_factor
except self._read_exception, exc:
- if field in fieldInfo:
- if fieldInfo[field].particle_type:
+ if field in self.pf.field_info:
+ if self.pf.field_info[field].particle_type:
self[field] = na.array([],dtype='int64')
- if fieldInfo[field].not_in_all:
+ if self.pf.field_info[field].not_in_all:
self[field] = na.zeros(self.ActiveDimensions, dtype='float64')
else: raise
else:
self._generate_field(field)
return self.data[field]
- def _generate_field(self, field):
- if fieldInfo.has_key(field):
- # First we check the validator
- try:
- fieldInfo[field].check_available(self)
- except NeedsGridType, ngt_exception:
- # This is only going to be raised if n_gz > 0
- n_gz = ngt_exception.ghost_zones
- f_gz = ngt_exception.fields
- gz_grid = self.retrieve_ghost_zones(n_gz, f_gz)
- temp_array = fieldInfo[field](gz_grid)
- sl = [slice(n_gz,-n_gz)] * 3
- self[field] = temp_array[sl]
- else:
- self[field] = fieldInfo[field](self)
- else: # Can't find the field, try as it might
- raise exceptions.KeyError, field
-
def _setup_dx(self):
# So first we figure out what the index is. We don't assume
# that dx=dy=dz , at least here. We probably do elsewhere.
@@ -100,6 +101,7 @@
self.data['dx'] = self.dx
self.data['dy'] = self.dy
self.data['dz'] = self.dz
+ self.data['dds'] = na.array([self.dx, self.dy, self.dz])
self._corners = self.hierarchy.gridCorners[:,:,id]
def _generate_overlap_masks(self, axis, LE, RE):
@@ -162,19 +164,20 @@
# Now we give it pointers to all of its attributes
# Note that to keep in line with Enzo, we have broken PEP-8
h = self.hierarchy # cache it
- self.Dimensions = h.gridDimensions[self.id-1]
- self.StartIndices = h.gridStartIndices[self.id-1]
- self.EndIndices = h.gridEndIndices[self.id-1]
- self.LeftEdge = h.gridLeftEdge[self.id-1]
- self.RightEdge = h.gridRightEdge[self.id-1]
- self.Level = h.gridLevels[self.id-1,0]
- self.Time = h.gridTimes[self.id-1,0]
- self.NumberOfParticles = h.gridNumberOfParticles[self.id-1,0]
+ my_ind = self.id - self._id_offset
+ self.Dimensions = h.gridDimensions[my_ind]
+ self.StartIndices = h.gridStartIndices[my_ind]
+ self.EndIndices = h.gridEndIndices[my_ind]
+ self.LeftEdge = h.gridLeftEdge[my_ind]
+ self.RightEdge = h.gridRightEdge[my_ind]
+ self.Level = h.gridLevels[my_ind,0]
+ self.Time = h.gridTimes[my_ind,0]
+ self.NumberOfParticles = h.gridNumberOfParticles[my_ind,0]
self.ActiveDimensions = (self.EndIndices - self.StartIndices + 1)
- self.Children = h.gridTree[self.id-1]
- pID = h.gridReverseTree[self.id-1]
+ self.Children = h.gridTree[my_ind]
+ pID = h.gridReverseTree[my_ind]
if pID != None and pID != -1:
- self.Parent = weakref.proxy(h.grids[pID - 1])
+ self.Parent = weakref.proxy(h.grids[pID - self._id_offset])
else:
self.Parent = None
@@ -395,21 +398,23 @@
space-filling tiling of grids, possibly due to the finite accuracy in a
standard Enzo hierarchy file.
"""
+ rf = self.pf["RefineBy"]
+ my_ind = self.id - self._id_offset
le = self.LeftEdge
- self.dx = self.Parent.dx/2.0
- self.dy = self.Parent.dy/2.0
- self.dz = self.Parent.dz/2.0
+ self.dx = self.Parent.dx/rf
+ self.dy = self.Parent.dy/rf
+ self.dz = self.Parent.dz/rf
ParentLeftIndex = na.rint((self.LeftEdge-self.Parent.LeftEdge)/self.Parent.dx)
- self.start_index = 2*(ParentLeftIndex + self.Parent.get_global_startindex()).astype('int64')
+ self.start_index = rf*(ParentLeftIndex + self.Parent.get_global_startindex()).astype('int64')
self.LeftEdge = self.Parent.LeftEdge + self.Parent.dx * ParentLeftIndex
self.RightEdge = self.LeftEdge + \
self.ActiveDimensions*na.array([self.dx,self.dy,self.dz])
- self.hierarchy.gridDxs[self.id-1,0] = self.dx
- self.hierarchy.gridDys[self.id-1,0] = self.dy
- self.hierarchy.gridDzs[self.id-1,0] = self.dz
- self.hierarchy.gridLeftEdge[self.id-1,:] = self.LeftEdge
- self.hierarchy.gridRightEdge[self.id-1,:] = self.RightEdge
- self.hierarchy.gridCorners[:,:,self.id-1] = na.array([ # Unroll!
+ self.hierarchy.gridDxs[my_ind,0] = self.dx
+ self.hierarchy.gridDys[my_ind,0] = self.dy
+ self.hierarchy.gridDzs[my_ind,0] = self.dz
+ self.hierarchy.gridLeftEdge[my_ind,:] = self.LeftEdge
+ self.hierarchy.gridRightEdge[my_ind,:] = self.RightEdge
+ self.hierarchy.gridCorners[:,:,my_ind] = na.array([ # Unroll!
[self.LeftEdge[0], self.LeftEdge[1], self.LeftEdge[2]],
[self.RightEdge[0], self.LeftEdge[1], self.LeftEdge[2]],
[self.RightEdge[0], self.RightEdge[1], self.LeftEdge[2]],
@@ -437,7 +442,7 @@
pdx = na.array([self.Parent.dx, self.Parent.dy, self.Parent.dz]).ravel()
start_index = (self.Parent.get_global_startindex()) + \
na.rint((self.LeftEdge - self.Parent.LeftEdge)/pdx)
- self.start_index = (start_index*2).astype('int64').ravel()
+ self.start_index = (start_index*self.pf["RefineBy"]).astype('int64').ravel()
return self.start_index
def set_filename(self, filename):
Modified: trunk/yt/lagos/DataReadingFuncs.py
==============================================================================
--- trunk/yt/lagos/DataReadingFuncs.py (original)
+++ trunk/yt/lagos/DataReadingFuncs.py Wed Nov 19 09:56:37 2008
@@ -166,6 +166,36 @@
def getExceptionInMemory():
return KeyError
+def readDataSlicePacked2D(self, grid, field, axis, coord):
+ """
+ Reads a slice through the HDF5 data
+
+ @param grid: Grid to slice
+ @type grid: L{EnzoGrid<EnzoGrid>}
+ @param field: field to get
+ @type field: string
+ @param sl: region to get
+ @type sl: SliceType
+ """
+ t = HDF5LightReader.ReadData(grid.filename, "/Grid%08i/%s" %
+ (grid.id, field)).transpose()
+ return t
+
+def readDataSlicePacked1D(self, grid, field, axis, coord):
+ """
+ Reads a slice through the HDF5 data
+
+ @param grid: Grid to slice
+ @type grid: L{EnzoGrid<EnzoGrid>}
+ @param field: field to get
+ @type field: string
+ @param sl: region to get
+ @type sl: SliceType
+ """
+ t = HDF5LightReader.ReadData(grid.filename, "/Grid%08i/%s" %
+ (grid.id, field))
+ return t
+
class BaseDataQueue(object):
def __init__(self):
@@ -254,3 +284,19 @@
def preload(self, grids, sets):
pass
+
+class DataQueuePacked2D(BaseDataQueue):
+ def _read_set(self, grid, field):
+ return HDF5LightReader.ReadData(grid.filename,
+ "/Grid%08i/%s" % (grid.id, field)).transpose()[:,:,None]
+
+ def modify(self, field):
+ pass
+
+class DataQueuePacked1D(BaseDataQueue):
+ def _read_set(self, grid, field):
+ return HDF5LightReader.ReadData(grid.filename,
+ "/Grid%08i/%s" % (grid.id, field)).transpose()[:,None,None]
+
+ def modify(self, field):
+ pass
Modified: trunk/yt/lagos/DerivedFields.py
==============================================================================
--- trunk/yt/lagos/DerivedFields.py (original)
+++ trunk/yt/lagos/DerivedFields.py Wed Nov 19 09:56:37 2008
@@ -52,6 +52,7 @@
Msun2g = 1.989e33
MJ_constant = (((5*kboltz)/(G*mh))**(1.5)) * (3/(4*pi))**(0.5) / Msun2g
rho_crit_now = 1.8788e-29 # g cm^-3
+axis_names = 'xyz'
class FieldInfoContainer: # We are all Borg.
_shared_state = {}
@@ -449,9 +450,12 @@
bulk_velocity = data.get_field_parameter("bulk_velocity")
if bulk_velocity == None:
bulk_velocity = na.zeros(3)
- return ( (data["x-velocity"]-bulk_velocity[0])**2.0 + \
- (data["y-velocity"]-bulk_velocity[1])**2.0 + \
- (data["z-velocity"]-bulk_velocity[2])**2.0 )**(1.0/2.0)
+ tr = (data["x-velocity"]-bulk_velocity[0])**2.0
+ if data.pf["TopGridRank"] > 1:
+ tr += (data["y-velocity"]-bulk_velocity[1])**2.0
+ if data.pf["TopGridRank"] > 2:
+ tr += (data["z-velocity"]-bulk_velocity[2])**2.0
+ return na.sqrt(tr)
add_field("VelocityMagnitude", take_log=False, units=r"\rm{cm}/\rm{s}")
def _TangentialOverVelocityMagnitude(field, data):
Modified: trunk/yt/lagos/HierarchyType.py
==============================================================================
--- trunk/yt/lagos/HierarchyType.py (original)
+++ trunk/yt/lagos/HierarchyType.py Wed Nov 19 09:56:37 2008
@@ -34,6 +34,8 @@
5: (readDataHDF5, readAllDataHDF5, getFieldsHDF5, readDataSliceHDF5, getExceptionHDF5),
6: (readDataPacked, readAllDataPacked, getFieldsPacked, readDataSlicePacked, getExceptionHDF5),
8: (readDataInMemory, readAllDataInMemory, getFieldsInMemory, readDataSliceInMemory, getExceptionInMemory),
+ 'enzo_packed_2d': (readDataPacked, readAllDataPacked, getFieldsPacked, readDataSlicePacked2D, getExceptionHDF5),
+ 'enzo_packed_1d': (readDataPacked, readAllDataPacked, getFieldsPacked, readDataSlicePacked1D, getExceptionHDF5),
}
class AMRHierarchy:
@@ -47,7 +49,6 @@
# 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()
@@ -252,10 +253,6 @@
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 maxGrid, maxCoord, maxVal, pos
@@ -555,7 +552,7 @@
if line.startswith("Grid "):
self.num_grids = testGridID = int(line.split("=")[-1])
break
- self.__guess_data_style(testGrid, testGridID)
+ self.__guess_data_style(pf["TopGridRank"], testGrid, testGridID)
# For some reason, r8 seems to want Float64
if pf.has_key("CompilerPrecision") \
and pf["CompilerPrecision"] == "r4":
@@ -572,7 +569,7 @@
self.grid = classobj("EnzoGrid",(EnzoGridBase,), dd)
AMRHierarchy._setup_classes(self, dd)
- def __guess_data_style(self, testGrid, testGridID):
+ def __guess_data_style(self, rank, testGrid, testGridID):
if self.data_style: return
if testGrid[0] != os.path.sep:
testGrid = os.path.join(self.directory, testGrid)
@@ -588,14 +585,24 @@
self.queue = DataQueueHDF4()
except:
list_of_sets = HDF5LightReader.ReadListOfDatasets(testGrid, "/")
- if len(list_of_sets) == 0:
+ if len(list_of_sets) == 0 and rank == 3:
mylog.debug("Detected packed HDF5")
self.data_style = 6
self.queue = DataQueuePackedHDF5()
- else:
+ elif len(list_of_sets) > 0 and rank == 3:
mylog.debug("Detected unpacked HDF5")
self.data_style = 5
self.queue = DataQueueHDF5()
+ elif len(list_of_sets) == 0 and rank == 2:
+ mylog.debug("Detect packed 2D")
+ self.data_style = 'enzo_packed_2d'
+ self.queue = DataQueuePacked2D()
+ elif len(list_of_sets) == 0 and rank == 1:
+ mylog.debug("Detect packed 1D")
+ self.data_style = 'enzo_packed_1d'
+ self.queue = DataQueuePacked1D()
+ else:
+ raise TypeError
def __setup_filemap(self, grid):
if not self.data_style == 6:
@@ -828,7 +835,7 @@
def _setup_field_lists(self):
field_list = self.get_data("/", "DataFields")
- if field_list == None:
+ if field_list is None:
mylog.info("Gathering a field list (this may take a moment.)")
field_list = sets.Set()
random_sample = self._generate_random_grids()
@@ -845,16 +852,16 @@
self.save_data(list(field_list),"/","DataFields")
self.field_list = list(field_list)
for field in self.field_list:
- if field in fieldInfo: continue
+ if field in self.parameter_file.field_info: continue
mylog.info("Adding %s to list of fields", field)
cf = None
if self.parameter_file.has_key(field):
cf = lambda d: d.convert(field)
add_field(field, lambda a, b: None, convert_function=cf)
self.derived_field_list = []
- for field in fieldInfo:
+ for field in self.parameter_file.field_info:
try:
- fd = fieldInfo[field].get_dependencies(pf = self.parameter_file)
+ fd = self.parameter_file.field_info[field].get_dependencies(pf = self.parameter_file)
except:
continue
available = na.all([f in self.field_list for f in fd.requested])
@@ -957,6 +964,27 @@
random_sample = na.mgrid[0:max(len(gg)-1,1)].astype("int32")
return gg[(random_sample,)]
+class EnzoHierarchy1D(EnzoHierarchy):
+ def __init__(self, *args, **kwargs):
+ EnzoHierarchy.__init__(self, *args, **kwargs)
+ self.gridRightEdge[:,1:3] = 1.0
+ self.gridDimensions[:,1:3] = 1.0
+ self.gridDys[:,0] = 1.0
+ self.gridDzs[:,0] = 1.0
+ for g in self.grids:
+ g._prepare_grid()
+ g._setup_dx()
+
+class EnzoHierarchy2D(EnzoHierarchy):
+ def __init__(self, *args, **kwargs):
+ EnzoHierarchy.__init__(self, *args, **kwargs)
+ self.gridRightEdge[:,2] = 1.0
+ self.gridDimensions[:,2] = 1.0
+ self.gridDzs[:,0] = 1.0
+ for g in self.grids:
+ g._prepare_grid()
+ g._setup_dx()
+
scanf_regex = {}
scanf_regex['e'] = r"[-+]?\d+\.?\d*?|\.\d+[eE][-+]?\d+?"
scanf_regex['g'] = scanf_regex['e']
@@ -982,7 +1010,7 @@
# http://www.reddit.com/r/Python/comments/6hj75/reverse_file_iterator/c03vms4
# Credit goes to "Brian" on Reddit
-def rblocks(f, blocksize=4096*256):
+def rblocks(f, blocksize=4096):
"""Read file as series of blocks from end of file to start.
The data itself is in normal order, only the order of the blocks is reversed.
Modified: trunk/yt/lagos/OutputTypes.py
==============================================================================
--- trunk/yt/lagos/OutputTypes.py (original)
+++ trunk/yt/lagos/OutputTypes.py Wed Nov 19 09:56:37 2008
@@ -114,6 +114,7 @@
Enzo-specific output, set at a fixed time.
"""
_hierarchy_class = EnzoHierarchy
+ _fieldinfo_class = EnzoFieldContainer
def __init__(self, filename, data_style=None,
parameter_override = None,
conversion_override = None):
@@ -139,6 +140,29 @@
if os.path.exists(cp):
self.cool = EnzoTable(cp, cool_out_key)
+ # Now fixes for different types of Hierarchies
+ # This includes changing the fieldinfo class!
+ if self["TopGridRank"] == 1: self._setup_1d()
+ elif self["TopGridRank"] == 2: self._setup_2d()
+
+ self.field_info = self._fieldinfo_class()
+
+ def _setup_1d(self):
+ self._hierarchy_class = EnzoHierarchy1D
+ self._fieldinfo_class = Enzo1DFieldContainer
+ self.parameters["DomainLeftEdge"] = \
+ na.concatenate([self["DomainLeftEdge"], [0.0, 0.0]])
+ self.parameters["DomainRightEdge"] = \
+ na.concatenate([self["DomainRightEdge"], [1.0, 1.0]])
+
+ def _setup_2d(self):
+ self._hierarchy_class = EnzoHierarchy2D
+ self._fieldinfo_class = Enzo2DFieldContainer
+ self.parameters["DomainLeftEdge"] = \
+ na.concatenate([self["DomainLeftEdge"], [0.0]])
+ self.parameters["DomainRightEdge"] = \
+ na.concatenate([self["DomainRightEdge"], [1.0]])
+
def _parse_parameter_file(self):
"""
Parses the parameter file and establishes the various
Modified: trunk/yt/lagos/ParallelTools.py
==============================================================================
--- trunk/yt/lagos/ParallelTools.py (original)
+++ trunk/yt/lagos/ParallelTools.py Wed Nov 19 09:56:37 2008
@@ -278,6 +278,7 @@
def _get_dependencies(self, fields):
deps = []
+ fi = self.pf.field_info
for field in fields:
- deps += ensure_list(fieldInfo[field].get_dependencies().requested)
+ deps += ensure_list(fi[field].get_dependencies().requested)
return list(set(deps))
Modified: trunk/yt/lagos/Profiles.py
==============================================================================
--- trunk/yt/lagos/Profiles.py (original)
+++ trunk/yt/lagos/Profiles.py Wed Nov 19 09:56:37 2008
@@ -51,6 +51,7 @@
class BinnedProfile(ParallelAnalysisInterface):
def __init__(self, data_source, lazy_reader):
self._data_source = data_source
+ self.pf = data_source.pf
self._data = {}
self._pdata = {}
self._lazy_reader = lazy_reader
@@ -144,7 +145,8 @@
data = []
for field in _field_mapping.get(field, (field,)):
if check_cut:
- if field in fieldInfo and fieldInfo[field].particle_type:
+ if field in self.pf.field_info \
+ and self.pf.field_info[field].particle_type:
pointI = self._data_source._get_particle_indices(source)
else:
pointI = self._data_source._get_point_indices(source)
Copied: trunk/yt/lagos/UniversalFields.py (from r949, /branches/yt-non-3d/yt/lagos/UniversalFields.py)
==============================================================================
--- /branches/yt-non-3d/yt/lagos/UniversalFields.py (original)
+++ trunk/yt/lagos/UniversalFields.py Wed Nov 19 09:56:37 2008
@@ -30,17 +30,16 @@
import inspect
import copy
-# All our math stuff here:
-try:
- import scipy.signal
-except ImportError:
- pass
-
from math import pi
from yt.funcs import *
from FieldInfoContainer import *
+try:
+ import cic_deposit
+except ImportError:
+ pass
+
mh = 1.67e-24 # g
me = 9.11e-28 # g
sigma_thompson = 6.65e-25 # cm^2
@@ -318,6 +317,12 @@
function=_CellMassCode,
convert_function=_convertCellMassCode)
+def _TotalMass(field,data):
+ return (data["Density"]+data["particle_density"]) * data["CellVolume"]
+add_field("TotalMassMsun", units=r"M_{\odot}",
+ function=_TotalMass,
+ convert_function=_convertCellMassMsun)
+
def _CellVolume(field, data):
if data['dx'].size == 1:
try:
Modified: trunk/yt/lagos/__init__.py
==============================================================================
--- trunk/yt/lagos/__init__.py (original)
+++ trunk/yt/lagos/__init__.py Wed Nov 19 09:56:37 2008
@@ -77,8 +77,14 @@
import PointCombine
import HDF5LightReader
from EnzoDefs import *
-from DerivedFields import *
+
from ParallelTools import *
+# Now our fields
+#from DerivedFields import *
+from FieldInfoContainer import *
+from UniversalFields import *
+from EnzoFields import *
+
from DerivedQuantities import DerivedQuantityCollection, GridChildMaskWrapper
from DataReadingFuncs import *
from ClusterFiles import *
Modified: trunk/yt/mods.py
==============================================================================
--- trunk/yt/mods.py (original)
+++ trunk/yt/mods.py Wed Nov 19 09:56:37 2008
@@ -40,9 +40,14 @@
# Now individual component imports from lagos
from yt.lagos import EnzoStaticOutput, \
BinnedProfile1D, BinnedProfile2D, BinnedProfile3D, \
- add_field, fieldInfo, \
+ add_field, FieldInfo, EnzoFieldInfo, Enzo2DFieldInfo, OrionFieldInfo, \
Clump, write_clump_hierarchy, find_clumps, write_clumps
+# This is a temporary solution -- in the future, we will allow the user to
+# select this via ytcfg.
+
+fieldInfo = lagos.fieldInfo = EnzoFieldInfo
+
# Now individual component imports from raven
from yt.raven import PlotCollection, PlotCollectionInteractive, \
QuiverCallback, ParticleCallback, ContourCallback, \
Modified: trunk/yt/raven/Callbacks.py
==============================================================================
--- trunk/yt/raven/Callbacks.py (original)
+++ trunk/yt/raven/Callbacks.py Wed Nov 19 09:56:37 2008
@@ -199,6 +199,8 @@
z = plot.data[self.field][wI]
if self.take_log: z=na.log10(z)
zi = self.de.Triangulation(x,y).nn_interpolator(z)(xi,yi)
+ print z.min(), z.max(), na.nanmin(z), na.nanmax(z)
+ print zi.min(), zi.max(), na.nanmin(zi), na.nanmax(zi)
plot._axes.contour(xi,yi,zi,self.ncont, **self.plot_args)
plot._axes.set_xlim(xx0,xx1)
plot._axes.set_ylim(yy0,yy1)
Modified: trunk/yt/raven/PlotCollection.py
==============================================================================
--- trunk/yt/raven/PlotCollection.py (original)
+++ trunk/yt/raven/PlotCollection.py Wed Nov 19 09:56:37 2008
@@ -391,6 +391,14 @@
p["Axis"] = "na"
return p
+ def add_ortho_ray(self, axis, coords, field, axes = None,
+ figure = None, **kwargs):
+ data_source = self.pf.h.ortho_ray(axis, coords, field)
+ p = self._add_plot(PlotTypes.LineQueryPlot(data_source,
+ [axis_names[axis], field], self._get_new_id(),
+ figure, axes, plot_options=kwargs))
+ return p
+
def _get_new_id(self):
self.__id_counter += 1
return self.__id_counter-1
Modified: trunk/yt/raven/PlotTypes.py
==============================================================================
--- trunk/yt/raven/PlotTypes.py (original)
+++ trunk/yt/raven/PlotTypes.py Wed Nov 19 09:56:37 2008
@@ -53,6 +53,7 @@
self.set_autoscale(True)
self.im = defaultdict(lambda: "")
self["ParameterFile"] = "%s" % self.data.pf
+ self.pf = self.data.pf
self.axis_names = {}
if not figure:
self._figure = matplotlib.figure.Figure(size)
@@ -196,7 +197,7 @@
def __setup_from_field(self, field):
self.set_log_field(field in lagos.log_fields
- or lagos.fieldInfo[field].take_log)
+ or self.pf.field_info[field].take_log)
self.axis_names["Z"] = field
def set_log_field(self, val):
@@ -246,7 +247,7 @@
self.data['pdx'],
self.data['pdy'],
self[self.axis_names["Z"]],
- int(width), int(width),
+ int(height), int(width),
(x0, x1, y0, y1),aa,self._period).transpose()
return buff
@@ -266,9 +267,10 @@
newmax = na.nanmax(buff)
if self.do_autoscale:
self.norm.autoscale(na.array((newmin,newmax)))
+ aspect = (self.ylim[1]-self.ylim[0])/(self.xlim[1]-self.xlim[0])
self.image = \
self._axes.imshow(buff, interpolation='nearest', norm = self.norm,
- aspect=1.0, picker=True, origin='lower')
+ aspect=aspect, picker=True, origin='lower')
self._reset_image_parameters()
self._run_callbacks()
@@ -338,7 +340,7 @@
def switch_z(self, field):
self.set_log_field(field in lagos.log_fields
- or lagos.fieldInfo[field].take_log)
+ or self.pf.field_info[field].take_log)
self.axis_names["Z"] = field
self._redraw_image()
@@ -362,11 +364,11 @@
return
field_name = self.axis_names["Z"]
data_label = r"$\rm{%s}" % field_name.replace("_","\hspace{0.5}")
- if lagos.fieldInfo.has_key(field_name):
+ if self.pf.field_info.has_key(field_name):
if self._projected:
- data_label += r"\/\/ (%s)" % (lagos.fieldInfo[field_name].get_projected_units())
+ data_label += r"\/\/ (%s)" % (self.pf.field_info[field_name].get_projected_units())
else:
- data_label += r"\/\/ (%s)" % (lagos.fieldInfo[field_name].get_units())
+ data_label += r"\/\/ (%s)" % (self.pf.field_info[field_name].get_units())
data_label += r"$"
if self.colorbar != None: self.colorbar.set_label(str(data_label))
@@ -389,8 +391,8 @@
return
field_name = self.axis_names["Z"]
data_label = r"$\rm{%s}" % field_name.replace("_","\hspace{0.5}")
- if lagos.fieldInfo.has_key(field_name):
- data_label += r"\/\/ (%s)" % (lagos.fieldInfo[field_name].get_units())
+ if self.pf.field_info.has_key(field_name):
+ data_label += r"\/\/ (%s)" % (self.pf.field_info[field_name].get_units())
data_label += r"$"
if self.colorbar != None: self.colorbar.set_label(str(data_label))
@@ -436,8 +438,8 @@
return
field_name = self.axis_names["Z"]
data_label = r"$\rm{%s}" % field_name.replace("_","\hspace{0.5}")
- if lagos.fieldInfo.has_key(field_name):
- data_label += r"\/\/ (%s)" % (lagos.fieldInfo[field_name].get_projected_units())
+ if self.pf.field_info.has_key(field_name):
+ data_label += r"\/\/ (%s)" % (self.pf.field_info[field_name].get_projected_units())
data_label += r"$"
if self.colorbar != None: self.colorbar.set_label(str(data_label))
@@ -553,7 +555,7 @@
class ProfilePlot(RavenPlot):
def setup_bins(self, field, func=None):
- if field in lagos.fieldInfo and lagos.fieldInfo[field].take_log:
+ if field in self.pf.field_info and self.pf.field_info[field].take_log:
log_field = True
if func: func('log')
else:
@@ -564,8 +566,8 @@
def autoset_label(self, field, func):
dataLabel = r"$\rm{%s}" % (field.replace("_","\hspace{0.5}"))
- if field in lagos.fieldInfo:
- dataLabel += r" (%s)" % (lagos.fieldInfo[field].get_units())
+ if field in self.pf.field_info:
+ dataLabel += r" (%s)" % (self.pf.field_info[field].get_units())
dataLabel += r"$"
func(str(dataLabel))
@@ -752,14 +754,86 @@
func(str(self.datalabel))
return
data_label = r"$\rm{%s}" % field_name.replace("_"," ")
- if field_name in lagos.fieldInfo:
- data_label += r"\/\/ (%s)" % (lagos.fieldInfo[field_name].get_units())
+ if field_name in self.pf.field_info:
+ data_label += r"\/\/ (%s)" % (self.pf.field_info[field_name].get_units())
data_label += r"$"
func(str(data_label))
def set_width(self, width, unit):
mylog.warning("Choosing not to change the width of a phase plot instance")
+class LineQueryPlot(RavenPlot):
+ _type_name = "LineQueryPlot"
+
+ def __init__(self, data, fields, id, ticker=None,
+ figure=None, axes=None, plot_options=None):
+ self._semi_unique_id = id
+ RavenPlot.__init__(self, data, fields, figure, axes)
+
+ self.axis_names["X"] = fields[0]
+ self.axis_names["Y"] = fields[1]
+
+ self._log_x = False
+ if fields[1] in self.pf.field_info and \
+ self.pf.field_info[fields[1]].take_log:
+ self._log_y = True
+ else:
+ self._log_y = False
+
+ if plot_options is None: plot_options = {}
+ self.plot_options = plot_options
+
+ def _generate_prefix(self, prefix):
+ self.prefix = "_".join([prefix, self._type_name,
+ str(self._semi_unique_id),
+ self.axis_names['X'], self.axis_names['Y']])
+ self["Field1"] = self.axis_names["X"]
+ self["Field2"] = self.axis_names["Y"]
+
+ def _redraw_image(self):
+ self._axes.clear()
+ if not self._log_x and not self._log_y:
+ func = self._axes.plot
+ elif self._log_x and not self._log_y:
+ func = self._axes.semilogx
+ elif not self._log_x and self._log_y:
+ func = self._axes.semilogy
+ elif self._log_x and self._log_y:
+ func = self._axes.loglog
+ indices = na.argsort(self.data[self.fields[0]])
+ func(self.data[self.fields[0]][indices],
+ self.data[self.fields[1]][indices],
+ **self.plot_options)
+ self.autoset_label(self.fields[0], self._axes.set_xlabel)
+ self.autoset_label(self.fields[1], self._axes.set_ylabel)
+ self._run_callbacks()
+
+ def set_log_field(self, val):
+ if val:
+ self._log_y = True
+ else:
+ self._log_y = False
+
+ def switch_x(self, field, weight="CellMassMsun", accumulation=False):
+ self.fields[0] = field
+ self.axis_names["X"] = field
+
+ def switch_z(self, field, weight="CellMassMsun", accumulation=False):
+ self.fields[1] = field
+ self.axis_names["Y"] = field
+
+ def autoset_label(self, field_name, func):
+ if self.datalabel != None:
+ func(str(self.datalabel))
+ return
+ data_label = r"$\rm{%s}" % field_name.replace("_"," ")
+ if field_name in self.pf.field_info:
+ data_label += r"\/\/ (%s)" % (self.pf.field_info[field_name].get_units())
+ data_label += r"$"
+ func(str(data_label))
+
+
+ switch_y = switch_z # Compatibility...
# Now we provide some convenience functions to get information about plots.
# With Matplotlib 0.98.x, the 'transforms' branch broke backwards
Modified: trunk/yt/raven/_MPL.c
==============================================================================
--- trunk/yt/raven/_MPL.c (original)
+++ trunk/yt/raven/_MPL.c Wed Nov 19 09:56:37 2008
@@ -127,7 +127,7 @@
//npy_float64 *gridded = (npy_float64 *) my_array->data;
for(i=0;i<rows;i++)for(j=0;j<cols;j++)
- *(npy_float64*) PyArray_GETPTR2(my_array, j, i) = 0.0;
+ *(npy_float64*) PyArray_GETPTR2(my_array, i, j) = 0.0;
for(p=0;p<nx;p++)
{
// these are cell-centered
Modified: trunk/yt/raven/__init__.py
==============================================================================
--- trunk/yt/raven/__init__.py (original)
+++ trunk/yt/raven/__init__.py Wed Nov 19 09:56:37 2008
@@ -37,6 +37,9 @@
except:
mylog.warning("Deliverator import failed; all deliverator actions will fail!")
+import matplotlib
+matplotlib.rc('contour', negative_linestyle='solid')
+
import matplotlib.image
import matplotlib.ticker
import matplotlib.axes
Modified: trunk/yt/reason/Functions.py
==============================================================================
--- trunk/yt/reason/Functions.py (original)
+++ trunk/yt/reason/Functions.py Wed Nov 19 09:56:37 2008
@@ -37,10 +37,11 @@
def QueryFields(outputfile, only_display_fields = False):
fields = []
for f in outputfile.hierarchy.derived_field_list:
- if f in lagos.fieldInfo and lagos.fieldInfo[f].particle_type: continue
+ if f in outputfile.field_info and
+ outputfile.field_info[f].particle_type: continue
if only_display_fields and \
- f in lagos.fieldInfo and \
- not lagos.fieldInfo[f].display_field: continue
+ f in outputfile.field_info[f] and \
+ not outputfile.field_info[f].display_field: continue
fields.append(f)
return sorted(fields)
More information about the yt-svn
mailing list