# HG changeset patch # User Matthew Turk # Date 1283553011 25200 # Branch yt # Node ID 9f8fb27e2fc129a992cf3ab779d7d3fbddd7fe3a # Parent a332a51c273f041fcd525127b56ad8f8bbc82f1d Fixes to Gadget data from the GDF format. Mostly stylistic, but also set the particle fields to be recognized as such, which means they can be accessed directly. This script should now work: import yt.frontends.gadget.data_structures as gds file='decay_100.gyt.hdf5' gf = gds.GadgetStaticOutput(file) dd = gf.h.all_data() px = dd["particle_position_x"] diff -r a332a51c273f041fcd525127b56ad8f8bbc82f1d -r 9f8fb27e2fc129a992cf3ab779d7d3fbddd7fe3a yt/frontends/gadget/data_structures.py --- a/yt/frontends/gadget/data_structures.py Thu Sep 02 21:42:44 2010 -0700 +++ b/yt/frontends/gadget/data_structures.py Fri Sep 03 15:30:11 2010 -0700 @@ -26,8 +26,8 @@ """ import h5py -import numpy as np -import ipdb +import numpy as na +from itertools import izip from yt.funcs import * from yt.data_objects.grid_patch import \ @@ -42,41 +42,25 @@ class GadgetGrid(AMRGridPatch): _id_offset = 0 def __init__(self, hierarchy, id,dimensions,left_index, - level,parent_id,particle_count,Parent=[],Children=[]): + level, parent_id, particle_count): AMRGridPatch.__init__(self, id, filename = hierarchy.filename, hierarchy = hierarchy) self.id = id - self.address = '/data/grid_%010i'%id - self.dimensions = dimensions - self.left_index = left_index - self.level = level + self.ActiveDimensions = dimensions + self.start_index = left_index self.Level = level - self.parent_id = parent_id - self.particle_count = particle_count + self._parent_id = parent_id + self.Parent = None # Only one parent per grid + self.Children = [] + self.NumberOfParticles = particle_count try: - padd =self.address+'/particles' + padd = '/data/grid_%010i/particles' % id self.particle_types = self.hierarchy._handle[padd].keys() except: self.particle_types = () self.filename = hierarchy.filename - self.Parent = Parent - self.Children = Children - - 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. - if len(self.Parent)>0: - self.dds = self.Parent[0].dds / self.pf.refine_by - else: - LE, RE = self.hierarchy.grid_left_edge[self.id,:], \ - self.hierarchy.grid_right_edge[self.id,:] - self.dds = np.array((RE-LE)/self.ActiveDimensions) - if self.pf.dimensionality < 2: self.dds[1] = 1.0 - if self.pf.dimensionality < 3: self.dds[2] = 1.0 - self.data['dx'], self.data['dy'], self.data['dz'] = self.dds - def __repr__(self): return 'GadgetGrid_%05i'%self.id @@ -116,30 +100,31 @@ def _parse_hierarchy(self): f = self._handle # shortcut - npa = np.array + npa = na.array - grid_levels = npa(f['/grid_level'],dtype='int') - self.grid_left_edge = npa(f['/grid_left_index'],dtype='int') - self.grid_dimensions = npa(f['/grid_dimensions'],dtype='int') - self.grid_right_edge = self.grid_left_edge + self.grid_dimensions - grid_particle_count = npa(f['/grid_particle_count'],dtype='int') - self.grid_parent_id = npa(f['/grid_parent_id'],dtype='int') - self.max_level = np.max(self.grid_levels) + self.grid_levels.flat[:] = f['/grid_level'][:].astype("int32") + self.grid_left_edge[:] = f['/grid_left_index'][:] + self.grid_dimensions[:] = f['/grid_dimensions'][:] + self.grid_right_edge[:] = self.grid_left_edge + self.grid_dimensions + self.grid_particle_count.flat[:] = f['/grid_particle_count'][:].astype("int32") + grid_parent_id = f['/grid_parent_id'][:] + self.max_level = na.max(self.grid_levels) - args = zip(range(self.num_grids),grid_levels,self.grid_parent_id, - self.grid_left_edge,self.grid_dimensions, grid_particle_count) - grids = [self.grid(self,j,d,le,lvl,p,n) for j,lvl,p, le, d, n in args] - self.grids = npa(grids,dtype='object') + args = izip(xrange(self.num_grids), self.grid_levels, + grid_parent_id, self.grid_left_edge, + self.grid_dimensions, self.grid_particle_count) + self.grids = na.array([self.grid(self,j,d,le,lvl,p,n) + for j,lvl,p, le, d, n in args], + dtype='object') def _populate_grid_objects(self): - ipdb.set_trace() for g in self.grids: - if g.parent_id>=0: - parent = self.grids[g.parent_id] - g.Parent = g.Parent + [parent,] - parent.Children = parent.Children + [g,] + if g._parent_id >= 0: + parent = self.grids[g._parent_id] + g.Parent = parent + parent.Children.append(g) for g in self.grids: g._prepare_grid() g._setup_dx() diff -r a332a51c273f041fcd525127b56ad8f8bbc82f1d -r 9f8fb27e2fc129a992cf3ab779d7d3fbddd7fe3a yt/frontends/gadget/fields.py --- a/yt/frontends/gadget/fields.py Thu Sep 02 21:42:44 2010 -0700 +++ b/yt/frontends/gadget/fields.py Fri Sep 03 15:30:11 2010 -0700 @@ -40,6 +40,26 @@ add_field = add_gadget_field +translation_dict = {"particle_position_x" : "position_x", + "particle_position_y" : "position_y", + "particle_position_z" : "position_z", + } + +def _generate_translation(mine, theirs): + pfield = mine.startswith("particle") + add_field(theirs, function=lambda a, b: b[mine], take_log=True, + particle_type = pfield) + +for f,v in translation_dict.items(): + if v not in GadgetFieldInfo: + # Note here that it's the yt field that we check for particle nature + pfield = f.startswith("particle") + add_field(v, function=lambda a,b: None, take_log=False, + validators = [ValidateDataField(v)], + particle_type = pfield) + print "Setting up translator from %s to %s" % (v, f) + _generate_translation(v, f) + #for f,v in translation_dict.items(): # add_field(f, function=lambda a,b: None, take_log=True, @@ -49,9 +69,21 @@ # validators = [ValidateDataField(v)], # units=r"\rm{cm}") + add_field("position_x", function=lambda a,b: None, take_log=True, validators = [ValidateDataField("position_x")], + particle_type = True, + units=r"\rm{cm}") + +add_field("position_y", function=lambda a,b: None, take_log=True, + validators = [ValidateDataField("position_y")], + particle_type = True, + units=r"\rm{cm}") + +add_field("position_z", function=lambda a,b: None, take_log=True, + validators = [ValidateDataField("position_z")], + particle_type = True, units=r"\rm{cm}") add_field("VEL", function=lambda a,b: None, take_log=True, diff -r a332a51c273f041fcd525127b56ad8f8bbc82f1d -r 9f8fb27e2fc129a992cf3ab779d7d3fbddd7fe3a yt/frontends/gadget/io.py --- a/yt/frontends/gadget/io.py Thu Sep 02 21:42:44 2010 -0700 +++ b/yt/frontends/gadget/io.py Fri Sep 03 15:30:11 2010 -0700 @@ -24,7 +24,7 @@ """ import h5py -import numpy as np +import numpy as na from yt.utilities.io_handler import \ BaseIOHandler @@ -32,15 +32,15 @@ class IOHandlerGadget(BaseIOHandler): _data_style = 'gadget_infrastructure' def _read_data_set(self, grid, field): - data = () + data = [] fh = h5py.File(grid.filename,mode='r') for ptype in grid.particle_types: - address = '%s/particles/%s/%s'%(grid.address,field) - dat += fh[address].read(), - if len(data)>0: - data = np.vstack(data) + address = '/data/grid_%010i/particles/%s/%s' % (grid.id, ptype, field) + data.append(fh[address][:]) + if len(data) > 0: + data = na.vstack(data) fh.close() - return np.array(data) + return na.array(data) def _read_field_names(self,grid): adr = grid.Address fh = h5py.File(grid.filename,mode='r')