[yt-svn] commit/yt: 12 new changesets
Bitbucket
commits-noreply at bitbucket.org
Tue Sep 4 19:07:01 PDT 2012
12 new commits in yt:
https://bitbucket.org/yt_analysis/yt/changeset/37540bd65173/
changeset: 37540bd65173
branch: yt
user: samskillman
date: 2012-08-31 21:21:51
summary: Adding preliminary Athena support. From what I can tell all basic functionality works. Is fragile to how you layout our data.
affected #: 11 files
diff -r 514f5a6a60c893b2a05f929150f352b153d68c65 -r 37540bd65173d499ac8d6397496e716753af2145 yt/frontends/athena/api.py
--- /dev/null
+++ b/yt/frontends/athena/api.py
@@ -0,0 +1,42 @@
+"""
+API for yt.frontends.athena
+
+Author: Samuel W. Skillman <samskillman at gmail.com>
+Affiliation: University of Colorado at Boulder
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: J.S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: MSU
+License:
+ Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+from .data_structures import \
+ AthenaGrid, \
+ AthenaHierarchy, \
+ AthenaStaticOutput
+
+from .fields import \
+ AthenaFieldInfo, \
+ KnownAthenaFields, \
+ add_athena_field
+
+from .io import \
+ IOHandlerAthena
diff -r 514f5a6a60c893b2a05f929150f352b153d68c65 -r 37540bd65173d499ac8d6397496e716753af2145 yt/frontends/athena/data_structures.py
--- /dev/null
+++ b/yt/frontends/athena/data_structures.py
@@ -0,0 +1,396 @@
+"""
+Data structures for Athena.
+
+Author: Samuel W. Skillman <samskillman at gmail.com>
+Affiliation: University of Colorado at Boulder
+Author: Matthew Turk <matthewturk at gmail.com>
+Author: J. S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2008-2011 Samuel W. Skillman, Matthew Turk, J. S. Oishi.
+ All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import h5py
+import numpy as na
+import weakref
+from yt.funcs import *
+from yt.data_objects.grid_patch import \
+ AMRGridPatch
+from yt.data_objects.hierarchy import \
+ AMRHierarchy
+from yt.data_objects.static_output import \
+ StaticOutput
+from yt.utilities.definitions import \
+ mpc_conversion, sec_conversion
+
+from .fields import AthenaFieldInfo, KnownAthenaFields
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, NullFunc
+import pdb
+
+def _get_convert(fname):
+ def _conv(data):
+ return data.convert(fname)
+ return _conv
+
+class AthenaGrid(AMRGridPatch):
+ _id_offset = 0
+ def __init__(self, id, hierarchy, level, start, dimensions):
+ df = hierarchy.storage_filename
+ if id == 0:
+ gname = 'id0/' + df + '.vtk'
+ else:
+ gname = 'id%i/' % id + df[:-5] + '-id%i'%id + df[-5:] + '.vtk'
+ AMRGridPatch.__init__(self, id, filename = gname,
+ hierarchy = hierarchy)
+ self.filename = gname
+ self.Parent = []
+ self.Children = []
+ self.Level = level
+ self.start_index = start.copy()
+ self.stop_index = self.start_index + dimensions
+ self.ActiveDimensions = dimensions.copy()
+
+ 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.
+ id = self.id - self._id_offset
+ if len(self.Parent) > 0:
+ self.dds = self.Parent[0].dds / self.pf.refine_by
+ else:
+ LE, RE = self.hierarchy.grid_left_edge[id,:], \
+ self.hierarchy.grid_right_edge[id,:]
+ self.dds = na.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.field_data['dx'], self.field_data['dy'], self.field_data['dz'] = self.dds
+
+def parse_line(line, grid):
+ # grid is a dictionary
+ splitup = line.strip().split()
+ if "vtk" in splitup:
+ grid['vtk_version'] = splitup[-1]
+ elif "Really" in splitup:
+ grid['time'] = splitup[-1]
+ elif 'PRIMITIVE' in splitup:
+ grid['time'] = float(splitup[4].rstrip(','))
+ grid['level'] = int(splitup[6].rstrip(','))
+ grid['domain'] = int(splitup[8].rstrip(','))
+ elif "DIMENSIONS" in splitup:
+ grid['dimensions'] = na.array(splitup[-3:]).astype('int')
+ elif "ORIGIN" in splitup:
+ grid['left_edge'] = na.array(splitup[-3:]).astype('float64')
+ elif "SPACING" in splitup:
+ grid['dds'] = na.array(splitup[-3:]).astype('float64')
+ elif "CELL_DATA" in splitup:
+ grid["ncells"] = int(splitup[-1])
+ elif "SCALARS" in splitup:
+ field = splitup[1]
+ grid['read_field'] = field
+ grid['read_type'] = 'scalar'
+ elif "VECTORS" in splitup:
+ field = splitup[1]
+ grid['read_field'] = field
+ grid['read_type'] = 'vector'
+
+
+
+class AthenaHierarchy(AMRHierarchy):
+
+ grid = AthenaGrid
+ _data_style='athena'
+
+ def __init__(self, pf, data_style='athena'):
+ self.parameter_file = weakref.proxy(pf)
+ self.data_style = data_style
+ # for now, the hierarchy file is the parameter file!
+ self.storage_filename = self.parameter_file.storage_filename
+ self.hierarchy_filename = self.parameter_file.filename
+ #self.directory = os.path.dirname(self.hierarchy_filename)
+ self._fhandle = file(self.hierarchy_filename,'rb')
+ AMRHierarchy.__init__(self,pf,data_style)
+
+ self._fhandle.close()
+
+ def _initialize_data_storage(self):
+ pass
+
+ def _detect_fields(self):
+ field_map = {}
+ f = open(self.hierarchy_filename,'rb')
+ line = f.readline()
+ while line != '':
+ splitup = line.strip().split()
+ if "DIMENSIONS" in splitup:
+ grid_dims = na.array(splitup[-3:]).astype('int')
+ line = f.readline()
+ continue
+ elif "CELL_DATA" in splitup:
+ grid_ncells = int(splitup[-1])
+ line = f.readline()
+ if na.prod(grid_dims) != grid_ncells:
+ grid_dims -= 1
+ grid_dims[grid_dims==0]=1
+ if na.prod(grid_dims) != grid_ncells:
+ mylog.error('product of dimensions %i not equal to number of cells %i' %
+ (na.prod(grid_dims), grid_ncells))
+ raise TypeError
+ break
+ else:
+ del line
+ line = f.readline()
+ read_table = False
+ read_table_offset = f.tell()
+ while line != '':
+ if len(line) == 0: break
+ splitup = line.strip().split()
+ if 'SCALARS' in splitup:
+ field = splitup[1]
+ if not read_table:
+ line = f.readline() # Read the lookup table line
+ read_table = True
+ field_map[field] = 'scalar',f.tell() - read_table_offset
+ read_table=False
+
+ elif 'VECTORS' in splitup:
+ field = splitup[1]
+ vfield = field+'_x'
+ field_map[vfield] = 'vector',f.tell() - read_table_offset
+ vfield = field+'_y'
+ field_map[vfield] = 'vector',f.tell() - read_table_offset
+ vfield = field+'_z'
+ field_map[vfield] = 'vector',f.tell() - read_table_offset
+ del line
+ line = f.readline()
+
+ f.close()
+ del f
+
+ self.field_list = field_map.keys()
+ self._field_map = field_map
+
+ def _setup_classes(self):
+ dd = self._get_data_reader_dict()
+ AMRHierarchy._setup_classes(self, dd)
+ self.object_types.sort()
+
+ def _count_grids(self):
+ self.num_grids = self.parameter_file.nvtk
+
+ def _parse_hierarchy(self):
+ f = open(self.hierarchy_filename,'rb')
+ grid = {}
+ grid['read_field'] = None
+ grid['read_type'] = None
+ table_read=False
+ line = f.readline()
+ while grid['read_field'] is None:
+ parse_line(line, grid)
+ if "SCALAR" in line.strip().split():
+ break
+ if "VECTOR" in line.strip().split():
+ break
+ if 'TABLE' in line.strip().split():
+ break
+ if len(line) == 0: break
+ del line
+ line = f.readline()
+ f.close()
+ del f
+
+ if na.prod(grid['dimensions']) != grid['ncells']:
+ grid['dimensions'] -= 1
+ grid['dimensions'][grid['dimensions']==0]=1
+ if na.prod(grid['dimensions']) != grid['ncells']:
+ mylog.error('product of dimensions %i not equal to number of cells %i' %
+ (na.prod(grid['dimensions']), grid['ncells']))
+ raise TypeError
+
+ dxs=[]
+ self.grids = na.empty(self.num_grids, dtype='object')
+ levels = na.zeros(self.num_grids, dtype='int32')
+ single_grid_width = grid['dds']*grid['dimensions']
+ grids_per_dim = (self.parameter_file.domain_width/single_grid_width).astype('int32')
+ glis = na.empty((self.num_grids,3), dtype='int64')
+ for i in range(self.num_grids):
+ procz = i/(grids_per_dim[0]*grids_per_dim[1])
+ procy = (i - procz*(grids_per_dim[0]*grids_per_dim[1]))/grids_per_dim[0]
+ procx = i - procz*(grids_per_dim[0]*grids_per_dim[1]) - procy*grids_per_dim[1]
+ glis[i, 0] = procx*grid['dimensions'][0]
+ glis[i, 1] = procy*grid['dimensions'][1]
+ glis[i, 2] = procz*grid['dimensions'][2]
+ gdims = na.ones_like(glis)
+ gdims[:] = grid['dimensions']
+ for i in range(levels.shape[0]):
+ self.grids[i] = self.grid(i, self, levels[i],
+ glis[i],
+ gdims[i])
+ self.grids[i]._level_id = levels[i]
+
+ dx = (self.parameter_file.domain_right_edge-
+ self.parameter_file.domain_left_edge)/self.parameter_file.domain_dimensions
+ dx = dx/self.parameter_file.refine_by**(levels[i])
+ dxs.append(grid['dds'])
+ dx = na.array(dxs)
+ self.grid_left_edge = self.parameter_file.domain_left_edge + dx*glis
+ self.grid_dimensions = gdims.astype("int32")
+ self.grid_right_edge = self.grid_left_edge + dx*self.grid_dimensions
+ self.grid_particle_count = na.zeros([self.num_grids, 1], dtype='int64')
+ del levels, glis, gdims
+
+ def _populate_grid_objects(self):
+ for g in self.grids:
+ g._prepare_grid()
+ g._setup_dx()
+
+ for g in self.grids:
+ g.Children = self._get_grid_children(g)
+ for g1 in g.Children:
+ g1.Parent.append(g)
+ self.max_level = self.grid_levels.max()
+
+ def _setup_derived_fields(self):
+ self.derived_field_list = []
+
+ def _get_grid_children(self, grid):
+ mask = na.zeros(self.num_grids, dtype='bool')
+ grids, grid_ind = self.get_box_grids(grid.LeftEdge, grid.RightEdge)
+ mask[grid_ind] = True
+ return [g for g in self.grids[mask] if g.Level == grid.Level + 1]
+
+class AthenaStaticOutput(StaticOutput):
+ _hierarchy_class = AthenaHierarchy
+ _fieldinfo_fallback = AthenaFieldInfo
+ _fieldinfo_known = KnownAthenaFields
+ _data_style = "athena"
+
+ def __init__(self, filename, data_style='athena',
+ storage_filename = None, parameters = {}):
+ StaticOutput.__init__(self, filename, data_style)
+ self.filename = filename
+ self.storage_filename = filename[4:-4]
+ self.specified_parameters = parameters
+
+ def _set_units(self):
+ """
+ Generates the conversion to various physical _units based on the parameter file
+ """
+ self.units = {}
+ self.time_units = {}
+ if len(self.parameters) == 0:
+ self._parse_parameter_file()
+ self.time_units['1'] = 1
+ self.units['1'] = 1.0
+ self.units['cm'] = 1.0
+ self.units['unitary'] = 1.0 / (self.domain_right_edge - self.domain_left_edge).max()
+ for unit in mpc_conversion.keys():
+ self.units[unit] = 1.0 * mpc_conversion[unit] / mpc_conversion["cm"]
+ for unit in sec_conversion.keys():
+ self.time_units[unit] = 1.0 / sec_conversion[unit]
+
+ # Here should read through and add fields.
+
+ #default_fields=['density']
+ # for field in self.field_list:
+ # self.units[field] = 1.0
+ # self._fieldinfo_known.add_field(field, function=NullFunc, take_log=False,
+ # units="", projected_units="",
+ # convert_function=None)
+
+ # This should be improved.
+ # self._handle = h5py.File(self.parameter_filename, "r")
+ # for field_name in self._handle["/field_types"]:
+ # current_field = self._handle["/field_types/%s" % field_name]
+ # try:
+ # self.units[field_name] = current_field.attrs['field_to_cgs']
+ # except:
+ # self.units[field_name] = 1.0
+ # try:
+ # current_fields_unit = current_field.attrs['field_units'][0]
+ # except:
+ # current_fields_unit = ""
+ # self._fieldinfo_known.add_field(field_name, function=NullFunc, take_log=False,
+ # units=current_fields_unit, projected_units="",
+ # convert_function=_get_convert(field_name))
+
+ # self._handle.close()
+ # del self._handle
+
+ def _parse_parameter_file(self):
+ self._handle = open(self.parameter_filename, "rb")
+ # Read the start of a grid to get simulation parameters.
+ grid = {}
+ grid['read_field'] = None
+ line = self._handle.readline()
+ while grid['read_field'] is None:
+ parse_line(line, grid)
+ if "SCALAR" in line.strip().split():
+ break
+ if "VECTOR" in line.strip().split():
+ break
+ if 'TABLE' in line.strip().split():
+ break
+ if len(line) == 0: break
+ del line
+ line = self._handle.readline()
+
+ self.domain_left_edge = grid['left_edge']
+ self.domain_right_edge = -grid['left_edge']
+ self.domain_width = self.domain_right_edge-self.domain_left_edge
+ self.domain_dimensions = self.domain_width/grid['dds']
+ refine_by = None
+ if refine_by is None: refine_by = 2
+ self.refine_by = refine_by
+ self.dimensionality = 3
+ self.current_time = grid["time"]
+ self.unique_identifier = None
+ self.cosmological_simulation = False
+ self.num_ghost_zones = 0
+ self.field_ordering = 'fortran'
+ self.boundary_conditions = [1]*6
+
+ self.nvtk = int(na.product(self.domain_dimensions/(grid['dimensions']-1)))
+
+ # if self.cosmological_simulation:
+ # self.current_redshift = sp["current_redshift"]
+ # self.omega_lambda = sp["omega_lambda"]
+ # self.omega_matter = sp["omega_matter"]
+ # self.hubble_constant = sp["hubble_constant"]
+ # else:
+ self.current_redshift = self.omega_lambda = self.omega_matter = \
+ self.hubble_constant = self.cosmological_simulation = 0.0
+ self.parameters['Time'] = self.current_time # Hardcode time conversion for now.
+ self.parameters["HydroMethod"] = 0 # Hardcode for now until field staggering is supported.
+ self._handle.close()
+ del self._handle
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ try:
+ fileh = file(args[0],'rb')
+ if "gridded_data_format" in fileh:
+ return True
+ except:
+ pass
+ return False
+
+ def __repr__(self):
+ return self.basename.rsplit(".", 1)[0]
+
diff -r 514f5a6a60c893b2a05f929150f352b153d68c65 -r 37540bd65173d499ac8d6397496e716753af2145 yt/frontends/athena/definitions.py
--- /dev/null
+++ b/yt/frontends/athena/definitions.py
@@ -0,0 +1,25 @@
+"""
+Various definitions for various other modules and routines
+
+Author: J. S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2008-2011 J.S. Oishi. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
diff -r 514f5a6a60c893b2a05f929150f352b153d68c65 -r 37540bd65173d499ac8d6397496e716753af2145 yt/frontends/athena/fields.py
--- /dev/null
+++ b/yt/frontends/athena/fields.py
@@ -0,0 +1,91 @@
+"""
+Athena-specific fields
+
+Author: Samuel W. Skillman <samskillman at gmail.com>
+Affiliation: University of Colorado at Boulder
+Author: J. S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2008-2011 Samuel W. Skillman, Matthew Turk, J. S. Oishi.
+ All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, \
+ FieldInfo, \
+ ValidateParameter, \
+ ValidateDataField, \
+ ValidateProperty, \
+ ValidateSpatial, \
+ ValidateGridType, \
+ NullFunc, \
+ TranslationFunc
+import yt.data_objects.universal_fields
+
+log_translation_dict = {"Density": "density",
+ "Pressure": "pressure"}
+
+translation_dict = {"x-velocity": "velocity_x",
+ "y-velocity": "velocity_y",
+ "z-velocity": "velocity_z"}
+
+# translation_dict = {"mag_field_x": "cell_centered_B_x ",
+# "mag_field_y": "cell_centered_B_y ",
+# "mag_field_z": "cell_centered_B_z "}
+
+AthenaFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
+add_field = AthenaFieldInfo.add_field
+
+KnownAthenaFields = FieldInfoContainer()
+add_athena_field = KnownAthenaFields.add_field
+
+add_athena_field("density", function=NullFunc, take_log=True,
+ units=r"\rm{g}/\rm{cm}^3",
+ projected_units =r"\rm{g}/\rm{cm}^2")
+
+add_athena_field("specific_energy", function=NullFunc, take_log=True,
+ units=r"\rm{erg}/\rm{g}")
+
+add_athena_field("pressure", function=NullFunc, take_log=True,
+ units=r"\rm{erg}/\rm{g}")
+
+add_athena_field("velocity_x", function=NullFunc, take_log=False,
+ units=r"\rm{cm}/\rm{s}")
+
+add_athena_field("velocity_y", function=NullFunc, take_log=False,
+ units=r"\rm{cm}/\rm{s}")
+
+add_athena_field("velocity_z", function=NullFunc, take_log=False,
+ units=r"\rm{cm}/\rm{s}")
+
+add_athena_field("mag_field_x", function=NullFunc, take_log=False,
+ units=r"\rm{cm}/\rm{s}")
+
+add_athena_field("mag_field_y", function=NullFunc, take_log=False,
+ units=r"\rm{cm}/\rm{s}")
+
+add_athena_field("mag_field_z", function=NullFunc, take_log=False,
+ units=r"\rm{cm}/\rm{s}")
+
+for f,v in log_translation_dict.items():
+ add_field(f, TranslationFunc(v), take_log=True)
+
+for f,v in translation_dict.items():
+ add_field(f, TranslationFunc(v), take_log=False)
+
diff -r 514f5a6a60c893b2a05f929150f352b153d68c65 -r 37540bd65173d499ac8d6397496e716753af2145 yt/frontends/athena/io.py
--- /dev/null
+++ b/yt/frontends/athena/io.py
@@ -0,0 +1,116 @@
+"""
+The data-file handling functions
+
+Author: Samuel W. Skillman <samskillman at gmail.com>
+Affiliation: University of Colorado at Boulder
+Author: Matthew Turk <matthewturk at gmail.com>
+Author: J. S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2007-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+from yt.utilities.io_handler import \
+ BaseIOHandler
+import numpy as na
+
+class IOHandlerAthena(BaseIOHandler):
+ _data_style = "athena"
+ _offset_string = 'data:offsets=0'
+ _data_string = 'data:datatype=0'
+ _read_table_offset = None
+
+ def _field_dict(self,fhandle):
+ keys = fhandle['field_types'].keys()
+ val = fhandle['field_types'].keys()
+ return dict(zip(keys,val))
+
+ def _read_field_names(self,grid):
+ pass
+
+ def _read_data_set(self,grid,field):
+ f = file(grid.filename, 'rb')
+ dtype, offset = grid.hierarchy._field_map[field]
+ grid_ncells = na.prod(grid.ActiveDimensions)
+ grid_dims = grid.ActiveDimensions
+ line = f.readline()
+ while True:
+ splitup = line.strip().split()
+ if 'CELL_DATA' in splitup:
+ f.readline()
+ read_table_offset = f.tell()
+ del line
+ break
+ del line; line = f.readline()
+
+
+ f.seek(read_table_offset+offset)
+ if dtype == 'scalar':
+ data = na.fromfile(f, dtype='>f4', count=grid_ncells).reshape(grid_dims,order='F').copy()
+ if dtype == 'vector':
+ data = na.fromfile(f, dtype='>f4', count=3*grid_ncells)
+ if '_x' in field:
+ data = data[0::3].reshape(grid_dims,order='F').copy()
+ elif '_y' in field:
+ data = data[1::3].reshape(grid_dims,order='F').copy()
+ elif '_z' in field:
+ data = data[2::3].reshape(grid_dims,order='F').copy()
+ f.close()
+ if grid.pf.field_ordering == 1:
+ return data.T
+ else:
+ return data
+
+ def _read_data_slice(self, grid, field, axis, coord):
+ sl = [slice(None), slice(None), slice(None)]
+ sl[axis] = slice(coord, coord + 1)
+ if grid.pf.field_ordering == 1:
+ sl.reverse()
+
+ f = file(grid.filename, 'rb')
+ dtype, offset = grid.hierarchy._field_map[field]
+ grid_ncells = na.prod(grid.ActiveDimensions)
+
+ line = f.readline()
+ while True:
+ splitup = line.strip().split()
+ if 'CELL_DATA' in splitup:
+ f.readline()
+ read_table_offset = f.tell()
+ del line
+ break
+ del line; line = f.readline()
+
+ f.seek(read_table_offset+offset)
+ if dtype == 'scalar':
+ data = na.fromfile(f, dtype='>f4', count=grid_ncells).reshape(grid.ActiveDimensions,order='F')[sl].copy()
+ if dtype == 'vector':
+ data = na.fromfile(f, dtype='>f4', count=3*grid_ncells)
+ if '_x' in field:
+ data = data[0::3].reshape(grid.ActiveDimensions,order='F')[sl].copy()
+ elif '_y' in field:
+ data = data[1::3].reshape(grid.ActiveDimensions,order='F')[sl].copy()
+ elif '_z' in field:
+ data = data[2::3].reshape(grid.ActiveDimensions,order='F')[sl].copy()
+
+ f.close()
+ if grid.pf.field_ordering == 1:
+ return data.T
+ else:
+ return data
+
diff -r 514f5a6a60c893b2a05f929150f352b153d68c65 -r 37540bd65173d499ac8d6397496e716753af2145 yt/frontends/athena/setup.py
--- /dev/null
+++ b/yt/frontends/athena/setup.py
@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+import setuptools
+import os
+import sys
+import os.path
+
+
+def configuration(parent_package='', top_path=None):
+ from numpy.distutils.misc_util import Configuration
+ config = Configuration('athena', parent_package, top_path)
+ config.make_config_py() # installs __config__.py
+ #config.make_svn_version_py()
+ return config
diff -r 514f5a6a60c893b2a05f929150f352b153d68c65 -r 37540bd65173d499ac8d6397496e716753af2145 yt/frontends/setup.py
--- a/yt/frontends/setup.py
+++ b/yt/frontends/setup.py
@@ -7,6 +7,7 @@
config = Configuration('frontends', parent_package, top_path)
config.make_config_py() # installs __config__.py
#config.make_svn_version_py()
+ config.add_subpackage("athena")
config.add_subpackage("gdf")
config.add_subpackage("chombo")
config.add_subpackage("enzo")
diff -r 514f5a6a60c893b2a05f929150f352b153d68c65 -r 37540bd65173d499ac8d6397496e716753af2145 yt/mods.py
--- a/yt/mods.py
+++ b/yt/mods.py
@@ -95,6 +95,9 @@
from yt.frontends.gdf.api import \
GDFStaticOutput, GDFFieldInfo, add_gdf_field
+from yt.frontends.athena.api import \
+ AthenaStaticOutput, AthenaFieldInfo, add_athena_field
+
from yt.frontends.art.api import \
ARTStaticOutput, ARTFieldInfo, add_art_field
diff -r 514f5a6a60c893b2a05f929150f352b153d68c65 -r 37540bd65173d499ac8d6397496e716753af2145 yt/utilities/io_handler.py
--- a/yt/utilities/io_handler.py
+++ b/yt/utilities/io_handler.py
@@ -42,6 +42,7 @@
def __init__(cls, name, b, d):
type.__init__(cls, name, b, d)
if hasattr(cls, "_data_style"):
+ print 'Registering Class ', cls ,' with datastyle ', cls._data_style
io_registry[cls._data_style] = cls
def __init__(self):
https://bitbucket.org/yt_analysis/yt/changeset/0128c2a41436/
changeset: 0128c2a41436
branch: yt
user: samskillman
date: 2012-08-31 22:48:44
summary: yt load, load() now work. Also now works for serial runs. Assumes all .vtk files are AthenaStaticOutput.
affected #: 2 files
diff -r 37540bd65173d499ac8d6397496e716753af2145 -r 0128c2a41436595a037b92ad572a2a539b0e77f3 yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -54,10 +54,13 @@
_id_offset = 0
def __init__(self, id, hierarchy, level, start, dimensions):
df = hierarchy.storage_filename
- if id == 0:
- gname = 'id0/' + df + '.vtk'
+ if 'id0' not in hierarchy.parameter_file.filename:
+ gname = hierarchy.parameter_file.filename
else:
- gname = 'id%i/' % id + df[:-5] + '-id%i'%id + df[-5:] + '.vtk'
+ if id == 0:
+ gname = 'id0/' + df + '.vtk'
+ else:
+ gname = 'id%i/' % id + df[:-5] + '-id%i'%id + df[-5:] + '.vtk'
AMRGridPatch.__init__(self, id, filename = gname,
hierarchy = hierarchy)
self.filename = gname
@@ -384,8 +387,7 @@
@classmethod
def _is_valid(self, *args, **kwargs):
try:
- fileh = file(args[0],'rb')
- if "gridded_data_format" in fileh:
+ if 'vtk' in args[0]:
return True
except:
pass
diff -r 37540bd65173d499ac8d6397496e716753af2145 -r 0128c2a41436595a037b92ad572a2a539b0e77f3 yt/utilities/io_handler.py
--- a/yt/utilities/io_handler.py
+++ b/yt/utilities/io_handler.py
@@ -42,7 +42,6 @@
def __init__(cls, name, b, d):
type.__init__(cls, name, b, d)
if hasattr(cls, "_data_style"):
- print 'Registering Class ', cls ,' with datastyle ', cls._data_style
io_registry[cls._data_style] = cls
def __init__(self):
https://bitbucket.org/yt_analysis/yt/changeset/c80c81adee02/
changeset: c80c81adee02
branch: yt
user: samskillman
date: 2012-09-03 19:59:56
summary: Allow over-riding of domain_right edge in parameters = {} dictionary. i.e. load('id0/athena.vtk', parameters={'domain_right_edge':[1.0,1.0,1.0]}
affected #: 1 file
diff -r 0128c2a41436595a037b92ad572a2a539b0e77f3 -r c80c81adee02a586cbda215b892e40ce65fa6bde yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -355,7 +355,11 @@
line = self._handle.readline()
self.domain_left_edge = grid['left_edge']
- self.domain_right_edge = -grid['left_edge']
+ try:
+ self.domain_right_edge = na.array(self.specified_parameters['domain_right_edge'])
+ except:
+ mylog.info("Please set 'domain_right_edge' in parameters dictionary argument " +
+ "if it is not equal to -domain_left_edge.")
self.domain_width = self.domain_right_edge-self.domain_left_edge
self.domain_dimensions = self.domain_width/grid['dds']
refine_by = None
https://bitbucket.org/yt_analysis/yt/changeset/0470c0e4eecf/
changeset: 0470c0e4eecf
branch: yt
user: samskillman
date: 2012-09-04 19:00:19
summary: na->np., remove pdb.
affected #: 2 files
diff -r c80c81adee02a586cbda215b892e40ce65fa6bde -r 0470c0e4eecff2bb942057822a6d2c1a6ee76ed0 yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -28,7 +28,7 @@
"""
import h5py
-import numpy as na
+import numpy as np
import weakref
from yt.funcs import *
from yt.data_objects.grid_patch import \
@@ -43,7 +43,6 @@
from .fields import AthenaFieldInfo, KnownAthenaFields
from yt.data_objects.field_info_container import \
FieldInfoContainer, NullFunc
-import pdb
def _get_convert(fname):
def _conv(data):
@@ -80,7 +79,7 @@
else:
LE, RE = self.hierarchy.grid_left_edge[id,:], \
self.hierarchy.grid_right_edge[id,:]
- self.dds = na.array((RE-LE)/self.ActiveDimensions)
+ 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.field_data['dx'], self.field_data['dy'], self.field_data['dz'] = self.dds
@@ -97,11 +96,11 @@
grid['level'] = int(splitup[6].rstrip(','))
grid['domain'] = int(splitup[8].rstrip(','))
elif "DIMENSIONS" in splitup:
- grid['dimensions'] = na.array(splitup[-3:]).astype('int')
+ grid['dimensions'] = np.array(splitup[-3:]).astype('int')
elif "ORIGIN" in splitup:
- grid['left_edge'] = na.array(splitup[-3:]).astype('float64')
+ grid['left_edge'] = np.array(splitup[-3:]).astype('float64')
elif "SPACING" in splitup:
- grid['dds'] = na.array(splitup[-3:]).astype('float64')
+ grid['dds'] = np.array(splitup[-3:]).astype('float64')
elif "CELL_DATA" in splitup:
grid["ncells"] = int(splitup[-1])
elif "SCALARS" in splitup:
@@ -128,7 +127,7 @@
self.hierarchy_filename = self.parameter_file.filename
#self.directory = os.path.dirname(self.hierarchy_filename)
self._fhandle = file(self.hierarchy_filename,'rb')
- AMRHierarchy.__init__(self,pf,data_style)
+ AMRHierarchy.__init__(self, pf, data_style)
self._fhandle.close()
@@ -142,18 +141,18 @@
while line != '':
splitup = line.strip().split()
if "DIMENSIONS" in splitup:
- grid_dims = na.array(splitup[-3:]).astype('int')
+ grid_dims = np.array(splitup[-3:]).astype('int')
line = f.readline()
continue
elif "CELL_DATA" in splitup:
grid_ncells = int(splitup[-1])
line = f.readline()
- if na.prod(grid_dims) != grid_ncells:
+ if np.prod(grid_dims) != grid_ncells:
grid_dims -= 1
grid_dims[grid_dims==0]=1
- if na.prod(grid_dims) != grid_ncells:
+ if np.prod(grid_dims) != grid_ncells:
mylog.error('product of dimensions %i not equal to number of cells %i' %
- (na.prod(grid_dims), grid_ncells))
+ (np.prod(grid_dims), grid_ncells))
raise TypeError
break
else:
@@ -218,20 +217,20 @@
f.close()
del f
- if na.prod(grid['dimensions']) != grid['ncells']:
+ if np.prod(grid['dimensions']) != grid['ncells']:
grid['dimensions'] -= 1
grid['dimensions'][grid['dimensions']==0]=1
- if na.prod(grid['dimensions']) != grid['ncells']:
+ if np.prod(grid['dimensions']) != grid['ncells']:
mylog.error('product of dimensions %i not equal to number of cells %i' %
- (na.prod(grid['dimensions']), grid['ncells']))
+ (np.prod(grid['dimensions']), grid['ncells']))
raise TypeError
dxs=[]
- self.grids = na.empty(self.num_grids, dtype='object')
- levels = na.zeros(self.num_grids, dtype='int32')
+ self.grids = np.empty(self.num_grids, dtype='object')
+ levels = np.zeros(self.num_grids, dtype='int32')
single_grid_width = grid['dds']*grid['dimensions']
grids_per_dim = (self.parameter_file.domain_width/single_grid_width).astype('int32')
- glis = na.empty((self.num_grids,3), dtype='int64')
+ glis = np.empty((self.num_grids,3), dtype='int64')
for i in range(self.num_grids):
procz = i/(grids_per_dim[0]*grids_per_dim[1])
procy = (i - procz*(grids_per_dim[0]*grids_per_dim[1]))/grids_per_dim[0]
@@ -239,7 +238,7 @@
glis[i, 0] = procx*grid['dimensions'][0]
glis[i, 1] = procy*grid['dimensions'][1]
glis[i, 2] = procz*grid['dimensions'][2]
- gdims = na.ones_like(glis)
+ gdims = np.ones_like(glis)
gdims[:] = grid['dimensions']
for i in range(levels.shape[0]):
self.grids[i] = self.grid(i, self, levels[i],
@@ -251,11 +250,11 @@
self.parameter_file.domain_left_edge)/self.parameter_file.domain_dimensions
dx = dx/self.parameter_file.refine_by**(levels[i])
dxs.append(grid['dds'])
- dx = na.array(dxs)
+ dx = np.array(dxs)
self.grid_left_edge = self.parameter_file.domain_left_edge + dx*glis
self.grid_dimensions = gdims.astype("int32")
self.grid_right_edge = self.grid_left_edge + dx*self.grid_dimensions
- self.grid_particle_count = na.zeros([self.num_grids, 1], dtype='int64')
+ self.grid_particle_count = np.zeros([self.num_grids, 1], dtype='int64')
del levels, glis, gdims
def _populate_grid_objects(self):
@@ -273,7 +272,7 @@
self.derived_field_list = []
def _get_grid_children(self, grid):
- mask = na.zeros(self.num_grids, dtype='bool')
+ mask = np.zeros(self.num_grids, dtype='bool')
grids, grid_ind = self.get_box_grids(grid.LeftEdge, grid.RightEdge)
mask[grid_ind] = True
return [g for g in self.grids[mask] if g.Level == grid.Level + 1]
@@ -356,7 +355,7 @@
self.domain_left_edge = grid['left_edge']
try:
- self.domain_right_edge = na.array(self.specified_parameters['domain_right_edge'])
+ self.domain_right_edge = np.array(self.specified_parameters['domain_right_edge'])
except:
mylog.info("Please set 'domain_right_edge' in parameters dictionary argument " +
"if it is not equal to -domain_left_edge.")
@@ -373,7 +372,7 @@
self.field_ordering = 'fortran'
self.boundary_conditions = [1]*6
- self.nvtk = int(na.product(self.domain_dimensions/(grid['dimensions']-1)))
+ self.nvtk = int(np.product(self.domain_dimensions/(grid['dimensions']-1)))
# if self.cosmological_simulation:
# self.current_redshift = sp["current_redshift"]
diff -r c80c81adee02a586cbda215b892e40ce65fa6bde -r 0470c0e4eecff2bb942057822a6d2c1a6ee76ed0 yt/frontends/athena/io.py
--- a/yt/frontends/athena/io.py
+++ b/yt/frontends/athena/io.py
@@ -27,7 +27,7 @@
"""
from yt.utilities.io_handler import \
BaseIOHandler
-import numpy as na
+import numpy as np
class IOHandlerAthena(BaseIOHandler):
_data_style = "athena"
@@ -46,7 +46,7 @@
def _read_data_set(self,grid,field):
f = file(grid.filename, 'rb')
dtype, offset = grid.hierarchy._field_map[field]
- grid_ncells = na.prod(grid.ActiveDimensions)
+ grid_ncells = np.prod(grid.ActiveDimensions)
grid_dims = grid.ActiveDimensions
line = f.readline()
while True:
@@ -61,9 +61,9 @@
f.seek(read_table_offset+offset)
if dtype == 'scalar':
- data = na.fromfile(f, dtype='>f4', count=grid_ncells).reshape(grid_dims,order='F').copy()
+ data = np.fromfile(f, dtype='>f4', count=grid_ncells).reshape(grid_dims,order='F').copy()
if dtype == 'vector':
- data = na.fromfile(f, dtype='>f4', count=3*grid_ncells)
+ data = np.fromfile(f, dtype='>f4', count=3*grid_ncells)
if '_x' in field:
data = data[0::3].reshape(grid_dims,order='F').copy()
elif '_y' in field:
@@ -84,7 +84,7 @@
f = file(grid.filename, 'rb')
dtype, offset = grid.hierarchy._field_map[field]
- grid_ncells = na.prod(grid.ActiveDimensions)
+ grid_ncells = np.prod(grid.ActiveDimensions)
line = f.readline()
while True:
@@ -98,9 +98,9 @@
f.seek(read_table_offset+offset)
if dtype == 'scalar':
- data = na.fromfile(f, dtype='>f4', count=grid_ncells).reshape(grid.ActiveDimensions,order='F')[sl].copy()
+ data = np.fromfile(f, dtype='>f4', count=grid_ncells).reshape(grid.ActiveDimensions,order='F')[sl].copy()
if dtype == 'vector':
- data = na.fromfile(f, dtype='>f4', count=3*grid_ncells)
+ data = np.fromfile(f, dtype='>f4', count=3*grid_ncells)
if '_x' in field:
data = data[0::3].reshape(grid.ActiveDimensions,order='F')[sl].copy()
elif '_y' in field:
https://bitbucket.org/yt_analysis/yt/changeset/c883c2104519/
changeset: c883c2104519
branch: yt
user: samskillman
date: 2012-09-04 19:23:36
summary: First pass on some of Matts suggestions on teh PR.
affected #: 1 file
diff -r 0470c0e4eecff2bb942057822a6d2c1a6ee76ed0 -r c883c21045190ad9e628d28542cfeb7460d4e581 yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -57,9 +57,9 @@
gname = hierarchy.parameter_file.filename
else:
if id == 0:
- gname = 'id0/' + df + '.vtk'
+ gname = 'id0/%s.vtk' % df
else:
- gname = 'id%i/' % id + df[:-5] + '-id%i'%id + df[-5:] + '.vtk'
+ gname = 'id%i/%s-id%i%s.vtk' % (id, df[:-5], id, df[-5:] )
AMRGridPatch.__init__(self, id, filename = gname,
hierarchy = hierarchy)
self.filename = gname
@@ -143,7 +143,6 @@
if "DIMENSIONS" in splitup:
grid_dims = np.array(splitup[-3:]).astype('int')
line = f.readline()
- continue
elif "CELL_DATA" in splitup:
grid_ncells = int(splitup[-1])
line = f.readline()
@@ -151,7 +150,7 @@
grid_dims -= 1
grid_dims[grid_dims==0]=1
if np.prod(grid_dims) != grid_ncells:
- mylog.error('product of dimensions %i not equal to number of cells %i' %
+ mylog.error('product of dimensions %i not equal to number of cells %i' %
(np.prod(grid_dims), grid_ncells))
raise TypeError
break
@@ -168,17 +167,14 @@
if not read_table:
line = f.readline() # Read the lookup table line
read_table = True
- field_map[field] = 'scalar',f.tell() - read_table_offset
+ field_map[field] = ('scalar', f.tell() - read_table_offset)
read_table=False
elif 'VECTORS' in splitup:
field = splitup[1]
- vfield = field+'_x'
- field_map[vfield] = 'vector',f.tell() - read_table_offset
- vfield = field+'_y'
- field_map[vfield] = 'vector',f.tell() - read_table_offset
- vfield = field+'_z'
- field_map[vfield] = 'vector',f.tell() - read_table_offset
+ for ax in 'xyz':
+ field_map["%s_%s" % (field, ax)] =\
+ ('vector', f.tell() - read_table_offset)
del line
line = f.readline()
@@ -217,6 +213,8 @@
f.close()
del f
+ # It seems some datasets have a mismatch between ncells and
+ # the actual grid dimensions.
if np.prod(grid['dimensions']) != grid['ncells']:
grid['dimensions'] -= 1
grid['dimensions'][grid['dimensions']==0]=1
@@ -244,7 +242,6 @@
self.grids[i] = self.grid(i, self, levels[i],
glis[i],
gdims[i])
- self.grids[i]._level_id = levels[i]
dx = (self.parameter_file.domain_right_edge-
self.parameter_file.domain_left_edge)/self.parameter_file.domain_dimensions
@@ -268,8 +265,8 @@
g1.Parent.append(g)
self.max_level = self.grid_levels.max()
- def _setup_derived_fields(self):
- self.derived_field_list = []
+# def _setup_derived_fields(self):
+# self.derived_field_list = []
def _get_grid_children(self, grid):
mask = np.zeros(self.num_grids, dtype='bool')
@@ -293,6 +290,7 @@
def _set_units(self):
"""
Generates the conversion to various physical _units based on the parameter file
+ This is a stub for future development. Currently sets arbitrary.
"""
self.units = {}
self.time_units = {}
@@ -307,34 +305,6 @@
for unit in sec_conversion.keys():
self.time_units[unit] = 1.0 / sec_conversion[unit]
- # Here should read through and add fields.
-
- #default_fields=['density']
- # for field in self.field_list:
- # self.units[field] = 1.0
- # self._fieldinfo_known.add_field(field, function=NullFunc, take_log=False,
- # units="", projected_units="",
- # convert_function=None)
-
- # This should be improved.
- # self._handle = h5py.File(self.parameter_filename, "r")
- # for field_name in self._handle["/field_types"]:
- # current_field = self._handle["/field_types/%s" % field_name]
- # try:
- # self.units[field_name] = current_field.attrs['field_to_cgs']
- # except:
- # self.units[field_name] = 1.0
- # try:
- # current_fields_unit = current_field.attrs['field_units'][0]
- # except:
- # current_fields_unit = ""
- # self._fieldinfo_known.add_field(field_name, function=NullFunc, take_log=False,
- # units=current_fields_unit, projected_units="",
- # convert_function=_get_convert(field_name))
-
- # self._handle.close()
- # del self._handle
-
def _parse_parameter_file(self):
self._handle = open(self.parameter_filename, "rb")
# Read the start of a grid to get simulation parameters.
@@ -366,7 +336,7 @@
self.refine_by = refine_by
self.dimensionality = 3
self.current_time = grid["time"]
- self.unique_identifier = None
+ self.unique_identifier = self._handle.__hash__()
self.cosmological_simulation = False
self.num_ghost_zones = 0
self.field_ordering = 'fortran'
@@ -374,18 +344,11 @@
self.nvtk = int(np.product(self.domain_dimensions/(grid['dimensions']-1)))
- # if self.cosmological_simulation:
- # self.current_redshift = sp["current_redshift"]
- # self.omega_lambda = sp["omega_lambda"]
- # self.omega_matter = sp["omega_matter"]
- # self.hubble_constant = sp["hubble_constant"]
- # else:
self.current_redshift = self.omega_lambda = self.omega_matter = \
self.hubble_constant = self.cosmological_simulation = 0.0
self.parameters['Time'] = self.current_time # Hardcode time conversion for now.
self.parameters["HydroMethod"] = 0 # Hardcode for now until field staggering is supported.
self._handle.close()
- del self._handle
@classmethod
def _is_valid(self, *args, **kwargs):
https://bitbucket.org/yt_analysis/yt/changeset/3abd0449f7a4/
changeset: 3abd0449f7a4
branch: yt
user: samskillman
date: 2012-09-04 19:25:14
summary: Fix for domain_right_edge.
affected #: 1 file
diff -r c883c21045190ad9e628d28542cfeb7460d4e581 -r 3abd0449f7a4add4a2a379aebe0647b47101a52e yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -329,6 +329,7 @@
except:
mylog.info("Please set 'domain_right_edge' in parameters dictionary argument " +
"if it is not equal to -domain_left_edge.")
+ self.domain_right_edge = -self.domain_left_edge
self.domain_width = self.domain_right_edge-self.domain_left_edge
self.domain_dimensions = self.domain_width/grid['dds']
refine_by = None
https://bitbucket.org/yt_analysis/yt/changeset/e2ba49fe07b7/
changeset: e2ba49fe07b7
branch: yt
user: samskillman
date: 2012-09-04 20:03:17
summary: More cleanup on aisle Athena.
affected #: 3 files
diff -r 3abd0449f7a4add4a2a379aebe0647b47101a52e -r e2ba49fe07b7680e7d536207a4a82853c8146eda yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -155,12 +155,10 @@
raise TypeError
break
else:
- del line
line = f.readline()
read_table = False
read_table_offset = f.tell()
while line != '':
- if len(line) == 0: break
splitup = line.strip().split()
if 'SCALARS' in splitup:
field = splitup[1]
@@ -175,11 +173,9 @@
for ax in 'xyz':
field_map["%s_%s" % (field, ax)] =\
('vector', f.tell() - read_table_offset)
- del line
line = f.readline()
f.close()
- del f
self.field_list = field_map.keys()
self._field_map = field_map
@@ -208,10 +204,8 @@
if 'TABLE' in line.strip().split():
break
if len(line) == 0: break
- del line
line = f.readline()
f.close()
- del f
# It seems some datasets have a mismatch between ncells and
# the actual grid dimensions.
@@ -252,7 +246,6 @@
self.grid_dimensions = gdims.astype("int32")
self.grid_right_edge = self.grid_left_edge + dx*self.grid_dimensions
self.grid_particle_count = np.zeros([self.num_grids, 1], dtype='int64')
- del levels, glis, gdims
def _populate_grid_objects(self):
for g in self.grids:
@@ -282,10 +275,10 @@
def __init__(self, filename, data_style='athena',
storage_filename = None, parameters = {}):
+ self.specified_parameters = parameters
StaticOutput.__init__(self, filename, data_style)
self.filename = filename
self.storage_filename = filename[4:-4]
- self.specified_parameters = parameters
def _set_units(self):
"""
@@ -298,7 +291,6 @@
self._parse_parameter_file()
self.time_units['1'] = 1
self.units['1'] = 1.0
- self.units['cm'] = 1.0
self.units['unitary'] = 1.0 / (self.domain_right_edge - self.domain_left_edge).max()
for unit in mpc_conversion.keys():
self.units[unit] = 1.0 * mpc_conversion[unit] / mpc_conversion["cm"]
@@ -320,13 +312,12 @@
if 'TABLE' in line.strip().split():
break
if len(line) == 0: break
- del line
line = self._handle.readline()
self.domain_left_edge = grid['left_edge']
- try:
+ if 'domain_right_edge' in self.specified_parameters:
self.domain_right_edge = np.array(self.specified_parameters['domain_right_edge'])
- except:
+ else:
mylog.info("Please set 'domain_right_edge' in parameters dictionary argument " +
"if it is not equal to -domain_left_edge.")
self.domain_right_edge = -self.domain_left_edge
diff -r 3abd0449f7a4add4a2a379aebe0647b47101a52e -r e2ba49fe07b7680e7d536207a4a82853c8146eda yt/frontends/athena/fields.py
--- a/yt/frontends/athena/fields.py
+++ b/yt/frontends/athena/fields.py
@@ -44,10 +44,10 @@
translation_dict = {"x-velocity": "velocity_x",
"y-velocity": "velocity_y",
"z-velocity": "velocity_z"}
-
-# translation_dict = {"mag_field_x": "cell_centered_B_x ",
-# "mag_field_y": "cell_centered_B_y ",
-# "mag_field_z": "cell_centered_B_z "}
+
+translation_dict = {"mag_field_x": "cell_centered_B_x ",
+ "mag_field_y": "cell_centered_B_y ",
+ "mag_field_z": "cell_centered_B_z "}
AthenaFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
add_field = AthenaFieldInfo.add_field
@@ -55,36 +55,33 @@
KnownAthenaFields = FieldInfoContainer()
add_athena_field = KnownAthenaFields.add_field
-add_athena_field("density", function=NullFunc, take_log=True,
- units=r"\rm{g}/\rm{cm}^3",
- projected_units =r"\rm{g}/\rm{cm}^2")
+add_athena_field("density", function=NullFunc, take_log=False,
+ units=r"",
+ projected_units =r"")
-add_athena_field("specific_energy", function=NullFunc, take_log=True,
- units=r"\rm{erg}/\rm{g}")
-
-add_athena_field("pressure", function=NullFunc, take_log=True,
- units=r"\rm{erg}/\rm{g}")
+add_athena_field("pressure", function=NullFunc, take_log=False,
+ units=r"")
add_athena_field("velocity_x", function=NullFunc, take_log=False,
- units=r"\rm{cm}/\rm{s}")
+ units=r"")
add_athena_field("velocity_y", function=NullFunc, take_log=False,
- units=r"\rm{cm}/\rm{s}")
+ units=r"")
add_athena_field("velocity_z", function=NullFunc, take_log=False,
- units=r"\rm{cm}/\rm{s}")
+ units=r"")
-add_athena_field("mag_field_x", function=NullFunc, take_log=False,
- units=r"\rm{cm}/\rm{s}")
+add_athena_field("cell_centered_B_x", function=NullFunc, take_log=False,
+ units=r"")
-add_athena_field("mag_field_y", function=NullFunc, take_log=False,
- units=r"\rm{cm}/\rm{s}")
+add_athena_field("cell_centered_B_y", function=NullFunc, take_log=False,
+ units=r"")
-add_athena_field("mag_field_z", function=NullFunc, take_log=False,
- units=r"\rm{cm}/\rm{s}")
+add_athena_field("cell_centered_B_z", function=NullFunc, take_log=False,
+ units=r"")
for f,v in log_translation_dict.items():
- add_field(f, TranslationFunc(v), take_log=True)
+ add_field(f, TranslationFunc(v), take_log=False)
for f,v in translation_dict.items():
add_field(f, TranslationFunc(v), take_log=False)
diff -r 3abd0449f7a4add4a2a379aebe0647b47101a52e -r e2ba49fe07b7680e7d536207a4a82853c8146eda yt/frontends/athena/io.py
--- a/yt/frontends/athena/io.py
+++ b/yt/frontends/athena/io.py
@@ -48,20 +48,11 @@
dtype, offset = grid.hierarchy._field_map[field]
grid_ncells = np.prod(grid.ActiveDimensions)
grid_dims = grid.ActiveDimensions
- line = f.readline()
- while True:
- splitup = line.strip().split()
- if 'CELL_DATA' in splitup:
- f.readline()
- read_table_offset = f.tell()
- del line
- break
- del line; line = f.readline()
-
-
+ read_table_offset = get_read_table_offset(f)
f.seek(read_table_offset+offset)
if dtype == 'scalar':
- data = np.fromfile(f, dtype='>f4', count=grid_ncells).reshape(grid_dims,order='F').copy()
+ data = np.fromfile(f, dtype='>f4',
+ count=grid_ncells).reshape(grid_dims,order='F').copy()
if dtype == 'vector':
data = np.fromfile(f, dtype='>f4', count=3*grid_ncells)
if '_x' in field:
@@ -86,19 +77,11 @@
dtype, offset = grid.hierarchy._field_map[field]
grid_ncells = np.prod(grid.ActiveDimensions)
- line = f.readline()
- while True:
- splitup = line.strip().split()
- if 'CELL_DATA' in splitup:
- f.readline()
- read_table_offset = f.tell()
- del line
- break
- del line; line = f.readline()
-
+ read_table_offset = get_read_table_offset(f)
f.seek(read_table_offset+offset)
if dtype == 'scalar':
- data = np.fromfile(f, dtype='>f4', count=grid_ncells).reshape(grid.ActiveDimensions,order='F')[sl].copy()
+ data = np.fromfile(f, dtype='>f4',
+ count=grid_ncells).reshape(grid.ActiveDimensions,order='F')[sl].copy()
if dtype == 'vector':
data = np.fromfile(f, dtype='>f4', count=3*grid_ncells)
if '_x' in field:
@@ -107,10 +90,18 @@
data = data[1::3].reshape(grid.ActiveDimensions,order='F')[sl].copy()
elif '_z' in field:
data = data[2::3].reshape(grid.ActiveDimensions,order='F')[sl].copy()
+ f.close()
+ return data
- f.close()
- if grid.pf.field_ordering == 1:
- return data.T
- else:
- return data
+def get_read_table_offset(f):
+ line = f.readline()
+ while True:
+ splitup = line.strip().split()
+ if 'CELL_DATA' in splitup:
+ f.readline()
+ read_table_offset = f.tell()
+ break
+ line = f.readline()
+ return read_table_offset
+
https://bitbucket.org/yt_analysis/yt/changeset/87f9eddc4718/
changeset: 87f9eddc4718
branch: yt
user: samskillman
date: 2012-09-04 20:25:47
summary: Some units modifications for Athena.
affected #: 1 file
diff -r e2ba49fe07b7680e7d536207a4a82853c8146eda -r 87f9eddc471822783826a87abd3714d0209c8ec9 yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -283,19 +283,21 @@
def _set_units(self):
"""
Generates the conversion to various physical _units based on the parameter file
- This is a stub for future development. Currently sets arbitrary.
"""
self.units = {}
self.time_units = {}
if len(self.parameters) == 0:
self._parse_parameter_file()
+ self._setup_nounits_units()
+ self.conversion_factors = defaultdict(lambda: 1.0)
self.time_units['1'] = 1
self.units['1'] = 1.0
self.units['unitary'] = 1.0 / (self.domain_right_edge - self.domain_left_edge).max()
+
+ def _setup_nounits_units(self):
+ self.conversion_factors["Time"] = 1.0
for unit in mpc_conversion.keys():
- self.units[unit] = 1.0 * mpc_conversion[unit] / mpc_conversion["cm"]
- for unit in sec_conversion.keys():
- self.time_units[unit] = 1.0 / sec_conversion[unit]
+ self.units[unit] = mpc_conversion[unit] / mpc_conversion["cm"]
def _parse_parameter_file(self):
self._handle = open(self.parameter_filename, "rb")
https://bitbucket.org/yt_analysis/yt/changeset/28de3d40037a/
changeset: 28de3d40037a
branch: yt
user: samskillman
date: 2012-09-04 21:08:42
summary: Fixing some athena fields. Set everything to linear space, not log, by default.
affected #: 1 file
diff -r 87f9eddc471822783826a87abd3714d0209c8ec9 -r 28de3d40037a7a938b026e828fa1999b3cd4605d yt/frontends/athena/fields.py
--- a/yt/frontends/athena/fields.py
+++ b/yt/frontends/athena/fields.py
@@ -38,14 +38,14 @@
TranslationFunc
import yt.data_objects.universal_fields
-log_translation_dict = {"Density": "density",
- "Pressure": "pressure"}
+log_translation_dict = {}
-translation_dict = {"x-velocity": "velocity_x",
+translation_dict = {"Density": "density",
+ "Pressure": "pressure",
+ "x-velocity": "velocity_x",
"y-velocity": "velocity_y",
- "z-velocity": "velocity_z"}
-
-translation_dict = {"mag_field_x": "cell_centered_B_x ",
+ "z-velocity": "velocity_z",
+ "mag_field_x": "cell_centered_B_x ",
"mag_field_y": "cell_centered_B_y ",
"mag_field_z": "cell_centered_B_z "}
@@ -81,7 +81,7 @@
units=r"")
for f,v in log_translation_dict.items():
- add_field(f, TranslationFunc(v), take_log=False)
+ add_field(f, TranslationFunc(v), take_log=True)
for f,v in translation_dict.items():
add_field(f, TranslationFunc(v), take_log=False)
https://bitbucket.org/yt_analysis/yt/changeset/74fd72ecdba4/
changeset: 74fd72ecdba4
branch: yt
user: samskillman
date: 2012-09-04 22:28:12
summary: Adding johns fix for the conserved flag.
affected #: 1 file
diff -r 28de3d40037a7a938b026e828fa1999b3cd4605d -r 74fd72ecdba41721ab8dce549dee968772db47b9 yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -91,7 +91,7 @@
grid['vtk_version'] = splitup[-1]
elif "Really" in splitup:
grid['time'] = splitup[-1]
- elif 'PRIMITIVE' in splitup:
+ elif any(x in ['PRIMITIVE','CONSERVED'] for x in splitup):
grid['time'] = float(splitup[4].rstrip(','))
grid['level'] = int(splitup[6].rstrip(','))
grid['domain'] = int(splitup[8].rstrip(','))
https://bitbucket.org/yt_analysis/yt/changeset/aaae6f8d3286/
changeset: aaae6f8d3286
branch: yt
user: samskillman
date: 2012-09-05 04:01:35
summary: Wrong dimension in index calculation, now parallel with non-trivial processor counts is working. Fixes issue pointed out by John Zuhone.
affected #: 1 file
diff -r 74fd72ecdba41721ab8dce549dee968772db47b9 -r aaae6f8d3286bfa838c1c5b7c1d385c2eb777d3e yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -112,8 +112,6 @@
grid['read_field'] = field
grid['read_type'] = 'vector'
-
-
class AthenaHierarchy(AMRHierarchy):
grid = AthenaGrid
@@ -226,7 +224,7 @@
for i in range(self.num_grids):
procz = i/(grids_per_dim[0]*grids_per_dim[1])
procy = (i - procz*(grids_per_dim[0]*grids_per_dim[1]))/grids_per_dim[0]
- procx = i - procz*(grids_per_dim[0]*grids_per_dim[1]) - procy*grids_per_dim[1]
+ procx = i - procz*(grids_per_dim[0]*grids_per_dim[1]) - procy*grids_per_dim[0]
glis[i, 0] = procx*grid['dimensions'][0]
glis[i, 1] = procy*grid['dimensions'][1]
glis[i, 2] = procz*grid['dimensions'][2]
https://bitbucket.org/yt_analysis/yt/changeset/0c844cb7df71/
changeset: 0c844cb7df71
branch: yt
user: jzuhone
date: 2012-09-05 04:06:59
summary: Merged in samskillman/yt (pull request #259)
affected #: 11 files
diff -r 3e35b1adb541c3f1412d3381d57bb5c0f83e3bf2 -r 0c844cb7df71949206ff497c57ef437d6429d58e yt/frontends/athena/api.py
--- /dev/null
+++ b/yt/frontends/athena/api.py
@@ -0,0 +1,42 @@
+"""
+API for yt.frontends.athena
+
+Author: Samuel W. Skillman <samskillman at gmail.com>
+Affiliation: University of Colorado at Boulder
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: J.S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Author: Britton Smith <brittonsmith at gmail.com>
+Affiliation: MSU
+License:
+ Copyright (C) 2010-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+from .data_structures import \
+ AthenaGrid, \
+ AthenaHierarchy, \
+ AthenaStaticOutput
+
+from .fields import \
+ AthenaFieldInfo, \
+ KnownAthenaFields, \
+ add_athena_field
+
+from .io import \
+ IOHandlerAthena
diff -r 3e35b1adb541c3f1412d3381d57bb5c0f83e3bf2 -r 0c844cb7df71949206ff497c57ef437d6429d58e yt/frontends/athena/data_structures.py
--- /dev/null
+++ b/yt/frontends/athena/data_structures.py
@@ -0,0 +1,356 @@
+"""
+Data structures for Athena.
+
+Author: Samuel W. Skillman <samskillman at gmail.com>
+Affiliation: University of Colorado at Boulder
+Author: Matthew Turk <matthewturk at gmail.com>
+Author: J. S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2008-2011 Samuel W. Skillman, Matthew Turk, J. S. Oishi.
+ All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+import h5py
+import numpy as np
+import weakref
+from yt.funcs import *
+from yt.data_objects.grid_patch import \
+ AMRGridPatch
+from yt.data_objects.hierarchy import \
+ AMRHierarchy
+from yt.data_objects.static_output import \
+ StaticOutput
+from yt.utilities.definitions import \
+ mpc_conversion, sec_conversion
+
+from .fields import AthenaFieldInfo, KnownAthenaFields
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, NullFunc
+
+def _get_convert(fname):
+ def _conv(data):
+ return data.convert(fname)
+ return _conv
+
+class AthenaGrid(AMRGridPatch):
+ _id_offset = 0
+ def __init__(self, id, hierarchy, level, start, dimensions):
+ df = hierarchy.storage_filename
+ if 'id0' not in hierarchy.parameter_file.filename:
+ gname = hierarchy.parameter_file.filename
+ else:
+ if id == 0:
+ gname = 'id0/%s.vtk' % df
+ else:
+ gname = 'id%i/%s-id%i%s.vtk' % (id, df[:-5], id, df[-5:] )
+ AMRGridPatch.__init__(self, id, filename = gname,
+ hierarchy = hierarchy)
+ self.filename = gname
+ self.Parent = []
+ self.Children = []
+ self.Level = level
+ self.start_index = start.copy()
+ self.stop_index = self.start_index + dimensions
+ self.ActiveDimensions = dimensions.copy()
+
+ 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.
+ id = self.id - self._id_offset
+ if len(self.Parent) > 0:
+ self.dds = self.Parent[0].dds / self.pf.refine_by
+ else:
+ LE, RE = self.hierarchy.grid_left_edge[id,:], \
+ self.hierarchy.grid_right_edge[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.field_data['dx'], self.field_data['dy'], self.field_data['dz'] = self.dds
+
+def parse_line(line, grid):
+ # grid is a dictionary
+ splitup = line.strip().split()
+ if "vtk" in splitup:
+ grid['vtk_version'] = splitup[-1]
+ elif "Really" in splitup:
+ grid['time'] = splitup[-1]
+ elif any(x in ['PRIMITIVE','CONSERVED'] for x in splitup):
+ grid['time'] = float(splitup[4].rstrip(','))
+ grid['level'] = int(splitup[6].rstrip(','))
+ grid['domain'] = int(splitup[8].rstrip(','))
+ elif "DIMENSIONS" in splitup:
+ grid['dimensions'] = np.array(splitup[-3:]).astype('int')
+ elif "ORIGIN" in splitup:
+ grid['left_edge'] = np.array(splitup[-3:]).astype('float64')
+ elif "SPACING" in splitup:
+ grid['dds'] = np.array(splitup[-3:]).astype('float64')
+ elif "CELL_DATA" in splitup:
+ grid["ncells"] = int(splitup[-1])
+ elif "SCALARS" in splitup:
+ field = splitup[1]
+ grid['read_field'] = field
+ grid['read_type'] = 'scalar'
+ elif "VECTORS" in splitup:
+ field = splitup[1]
+ grid['read_field'] = field
+ grid['read_type'] = 'vector'
+
+class AthenaHierarchy(AMRHierarchy):
+
+ grid = AthenaGrid
+ _data_style='athena'
+
+ def __init__(self, pf, data_style='athena'):
+ self.parameter_file = weakref.proxy(pf)
+ self.data_style = data_style
+ # for now, the hierarchy file is the parameter file!
+ self.storage_filename = self.parameter_file.storage_filename
+ self.hierarchy_filename = self.parameter_file.filename
+ #self.directory = os.path.dirname(self.hierarchy_filename)
+ self._fhandle = file(self.hierarchy_filename,'rb')
+ AMRHierarchy.__init__(self, pf, data_style)
+
+ self._fhandle.close()
+
+ def _initialize_data_storage(self):
+ pass
+
+ def _detect_fields(self):
+ field_map = {}
+ f = open(self.hierarchy_filename,'rb')
+ line = f.readline()
+ while line != '':
+ splitup = line.strip().split()
+ if "DIMENSIONS" in splitup:
+ grid_dims = np.array(splitup[-3:]).astype('int')
+ line = f.readline()
+ elif "CELL_DATA" in splitup:
+ grid_ncells = int(splitup[-1])
+ line = f.readline()
+ if np.prod(grid_dims) != grid_ncells:
+ grid_dims -= 1
+ grid_dims[grid_dims==0]=1
+ if np.prod(grid_dims) != grid_ncells:
+ mylog.error('product of dimensions %i not equal to number of cells %i' %
+ (np.prod(grid_dims), grid_ncells))
+ raise TypeError
+ break
+ else:
+ line = f.readline()
+ read_table = False
+ read_table_offset = f.tell()
+ while line != '':
+ splitup = line.strip().split()
+ if 'SCALARS' in splitup:
+ field = splitup[1]
+ if not read_table:
+ line = f.readline() # Read the lookup table line
+ read_table = True
+ field_map[field] = ('scalar', f.tell() - read_table_offset)
+ read_table=False
+
+ elif 'VECTORS' in splitup:
+ field = splitup[1]
+ for ax in 'xyz':
+ field_map["%s_%s" % (field, ax)] =\
+ ('vector', f.tell() - read_table_offset)
+ line = f.readline()
+
+ f.close()
+
+ self.field_list = field_map.keys()
+ self._field_map = field_map
+
+ def _setup_classes(self):
+ dd = self._get_data_reader_dict()
+ AMRHierarchy._setup_classes(self, dd)
+ self.object_types.sort()
+
+ def _count_grids(self):
+ self.num_grids = self.parameter_file.nvtk
+
+ def _parse_hierarchy(self):
+ f = open(self.hierarchy_filename,'rb')
+ grid = {}
+ grid['read_field'] = None
+ grid['read_type'] = None
+ table_read=False
+ line = f.readline()
+ while grid['read_field'] is None:
+ parse_line(line, grid)
+ if "SCALAR" in line.strip().split():
+ break
+ if "VECTOR" in line.strip().split():
+ break
+ if 'TABLE' in line.strip().split():
+ break
+ if len(line) == 0: break
+ line = f.readline()
+ f.close()
+
+ # It seems some datasets have a mismatch between ncells and
+ # the actual grid dimensions.
+ if np.prod(grid['dimensions']) != grid['ncells']:
+ grid['dimensions'] -= 1
+ grid['dimensions'][grid['dimensions']==0]=1
+ if np.prod(grid['dimensions']) != grid['ncells']:
+ mylog.error('product of dimensions %i not equal to number of cells %i' %
+ (np.prod(grid['dimensions']), grid['ncells']))
+ raise TypeError
+
+ dxs=[]
+ self.grids = np.empty(self.num_grids, dtype='object')
+ levels = np.zeros(self.num_grids, dtype='int32')
+ single_grid_width = grid['dds']*grid['dimensions']
+ grids_per_dim = (self.parameter_file.domain_width/single_grid_width).astype('int32')
+ glis = np.empty((self.num_grids,3), dtype='int64')
+ for i in range(self.num_grids):
+ procz = i/(grids_per_dim[0]*grids_per_dim[1])
+ procy = (i - procz*(grids_per_dim[0]*grids_per_dim[1]))/grids_per_dim[0]
+ procx = i - procz*(grids_per_dim[0]*grids_per_dim[1]) - procy*grids_per_dim[0]
+ glis[i, 0] = procx*grid['dimensions'][0]
+ glis[i, 1] = procy*grid['dimensions'][1]
+ glis[i, 2] = procz*grid['dimensions'][2]
+ gdims = np.ones_like(glis)
+ gdims[:] = grid['dimensions']
+ for i in range(levels.shape[0]):
+ self.grids[i] = self.grid(i, self, levels[i],
+ glis[i],
+ gdims[i])
+
+ dx = (self.parameter_file.domain_right_edge-
+ self.parameter_file.domain_left_edge)/self.parameter_file.domain_dimensions
+ dx = dx/self.parameter_file.refine_by**(levels[i])
+ dxs.append(grid['dds'])
+ dx = np.array(dxs)
+ self.grid_left_edge = self.parameter_file.domain_left_edge + dx*glis
+ self.grid_dimensions = gdims.astype("int32")
+ self.grid_right_edge = self.grid_left_edge + dx*self.grid_dimensions
+ self.grid_particle_count = np.zeros([self.num_grids, 1], dtype='int64')
+
+ def _populate_grid_objects(self):
+ for g in self.grids:
+ g._prepare_grid()
+ g._setup_dx()
+
+ for g in self.grids:
+ g.Children = self._get_grid_children(g)
+ for g1 in g.Children:
+ g1.Parent.append(g)
+ self.max_level = self.grid_levels.max()
+
+# def _setup_derived_fields(self):
+# self.derived_field_list = []
+
+ def _get_grid_children(self, grid):
+ mask = np.zeros(self.num_grids, dtype='bool')
+ grids, grid_ind = self.get_box_grids(grid.LeftEdge, grid.RightEdge)
+ mask[grid_ind] = True
+ return [g for g in self.grids[mask] if g.Level == grid.Level + 1]
+
+class AthenaStaticOutput(StaticOutput):
+ _hierarchy_class = AthenaHierarchy
+ _fieldinfo_fallback = AthenaFieldInfo
+ _fieldinfo_known = KnownAthenaFields
+ _data_style = "athena"
+
+ def __init__(self, filename, data_style='athena',
+ storage_filename = None, parameters = {}):
+ self.specified_parameters = parameters
+ StaticOutput.__init__(self, filename, data_style)
+ self.filename = filename
+ self.storage_filename = filename[4:-4]
+
+ def _set_units(self):
+ """
+ Generates the conversion to various physical _units based on the parameter file
+ """
+ self.units = {}
+ self.time_units = {}
+ if len(self.parameters) == 0:
+ self._parse_parameter_file()
+ self._setup_nounits_units()
+ self.conversion_factors = defaultdict(lambda: 1.0)
+ self.time_units['1'] = 1
+ self.units['1'] = 1.0
+ self.units['unitary'] = 1.0 / (self.domain_right_edge - self.domain_left_edge).max()
+
+ def _setup_nounits_units(self):
+ self.conversion_factors["Time"] = 1.0
+ for unit in mpc_conversion.keys():
+ self.units[unit] = mpc_conversion[unit] / mpc_conversion["cm"]
+
+ def _parse_parameter_file(self):
+ self._handle = open(self.parameter_filename, "rb")
+ # Read the start of a grid to get simulation parameters.
+ grid = {}
+ grid['read_field'] = None
+ line = self._handle.readline()
+ while grid['read_field'] is None:
+ parse_line(line, grid)
+ if "SCALAR" in line.strip().split():
+ break
+ if "VECTOR" in line.strip().split():
+ break
+ if 'TABLE' in line.strip().split():
+ break
+ if len(line) == 0: break
+ line = self._handle.readline()
+
+ self.domain_left_edge = grid['left_edge']
+ if 'domain_right_edge' in self.specified_parameters:
+ self.domain_right_edge = np.array(self.specified_parameters['domain_right_edge'])
+ else:
+ mylog.info("Please set 'domain_right_edge' in parameters dictionary argument " +
+ "if it is not equal to -domain_left_edge.")
+ self.domain_right_edge = -self.domain_left_edge
+ self.domain_width = self.domain_right_edge-self.domain_left_edge
+ self.domain_dimensions = self.domain_width/grid['dds']
+ refine_by = None
+ if refine_by is None: refine_by = 2
+ self.refine_by = refine_by
+ self.dimensionality = 3
+ self.current_time = grid["time"]
+ self.unique_identifier = self._handle.__hash__()
+ self.cosmological_simulation = False
+ self.num_ghost_zones = 0
+ self.field_ordering = 'fortran'
+ self.boundary_conditions = [1]*6
+
+ self.nvtk = int(np.product(self.domain_dimensions/(grid['dimensions']-1)))
+
+ self.current_redshift = self.omega_lambda = self.omega_matter = \
+ self.hubble_constant = self.cosmological_simulation = 0.0
+ self.parameters['Time'] = self.current_time # Hardcode time conversion for now.
+ self.parameters["HydroMethod"] = 0 # Hardcode for now until field staggering is supported.
+ self._handle.close()
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ try:
+ if 'vtk' in args[0]:
+ return True
+ except:
+ pass
+ return False
+
+ def __repr__(self):
+ return self.basename.rsplit(".", 1)[0]
+
diff -r 3e35b1adb541c3f1412d3381d57bb5c0f83e3bf2 -r 0c844cb7df71949206ff497c57ef437d6429d58e yt/frontends/athena/definitions.py
--- /dev/null
+++ b/yt/frontends/athena/definitions.py
@@ -0,0 +1,25 @@
+"""
+Various definitions for various other modules and routines
+
+Author: J. S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2008-2011 J.S. Oishi. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
diff -r 3e35b1adb541c3f1412d3381d57bb5c0f83e3bf2 -r 0c844cb7df71949206ff497c57ef437d6429d58e yt/frontends/athena/fields.py
--- /dev/null
+++ b/yt/frontends/athena/fields.py
@@ -0,0 +1,88 @@
+"""
+Athena-specific fields
+
+Author: Samuel W. Skillman <samskillman at gmail.com>
+Affiliation: University of Colorado at Boulder
+Author: J. S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2008-2011 Samuel W. Skillman, Matthew Turk, J. S. Oishi.
+ All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from yt.data_objects.field_info_container import \
+ FieldInfoContainer, \
+ FieldInfo, \
+ ValidateParameter, \
+ ValidateDataField, \
+ ValidateProperty, \
+ ValidateSpatial, \
+ ValidateGridType, \
+ NullFunc, \
+ TranslationFunc
+import yt.data_objects.universal_fields
+
+log_translation_dict = {}
+
+translation_dict = {"Density": "density",
+ "Pressure": "pressure",
+ "x-velocity": "velocity_x",
+ "y-velocity": "velocity_y",
+ "z-velocity": "velocity_z",
+ "mag_field_x": "cell_centered_B_x ",
+ "mag_field_y": "cell_centered_B_y ",
+ "mag_field_z": "cell_centered_B_z "}
+
+AthenaFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
+add_field = AthenaFieldInfo.add_field
+
+KnownAthenaFields = FieldInfoContainer()
+add_athena_field = KnownAthenaFields.add_field
+
+add_athena_field("density", function=NullFunc, take_log=False,
+ units=r"",
+ projected_units =r"")
+
+add_athena_field("pressure", function=NullFunc, take_log=False,
+ units=r"")
+
+add_athena_field("velocity_x", function=NullFunc, take_log=False,
+ units=r"")
+
+add_athena_field("velocity_y", function=NullFunc, take_log=False,
+ units=r"")
+
+add_athena_field("velocity_z", function=NullFunc, take_log=False,
+ units=r"")
+
+add_athena_field("cell_centered_B_x", function=NullFunc, take_log=False,
+ units=r"")
+
+add_athena_field("cell_centered_B_y", function=NullFunc, take_log=False,
+ units=r"")
+
+add_athena_field("cell_centered_B_z", function=NullFunc, take_log=False,
+ units=r"")
+
+for f,v in log_translation_dict.items():
+ add_field(f, TranslationFunc(v), take_log=True)
+
+for f,v in translation_dict.items():
+ add_field(f, TranslationFunc(v), take_log=False)
+
diff -r 3e35b1adb541c3f1412d3381d57bb5c0f83e3bf2 -r 0c844cb7df71949206ff497c57ef437d6429d58e yt/frontends/athena/io.py
--- /dev/null
+++ b/yt/frontends/athena/io.py
@@ -0,0 +1,107 @@
+"""
+The data-file handling functions
+
+Author: Samuel W. Skillman <samskillman at gmail.com>
+Affiliation: University of Colorado at Boulder
+Author: Matthew Turk <matthewturk at gmail.com>
+Author: J. S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt-project.org/
+License:
+ Copyright (C) 2007-2011 Matthew Turk. All Rights Reserved.
+
+ This file is part of yt.
+
+ yt is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+"""
+from yt.utilities.io_handler import \
+ BaseIOHandler
+import numpy as np
+
+class IOHandlerAthena(BaseIOHandler):
+ _data_style = "athena"
+ _offset_string = 'data:offsets=0'
+ _data_string = 'data:datatype=0'
+ _read_table_offset = None
+
+ def _field_dict(self,fhandle):
+ keys = fhandle['field_types'].keys()
+ val = fhandle['field_types'].keys()
+ return dict(zip(keys,val))
+
+ def _read_field_names(self,grid):
+ pass
+
+ def _read_data_set(self,grid,field):
+ f = file(grid.filename, 'rb')
+ dtype, offset = grid.hierarchy._field_map[field]
+ grid_ncells = np.prod(grid.ActiveDimensions)
+ grid_dims = grid.ActiveDimensions
+ read_table_offset = get_read_table_offset(f)
+ f.seek(read_table_offset+offset)
+ if dtype == 'scalar':
+ data = np.fromfile(f, dtype='>f4',
+ count=grid_ncells).reshape(grid_dims,order='F').copy()
+ if dtype == 'vector':
+ data = np.fromfile(f, dtype='>f4', count=3*grid_ncells)
+ if '_x' in field:
+ data = data[0::3].reshape(grid_dims,order='F').copy()
+ elif '_y' in field:
+ data = data[1::3].reshape(grid_dims,order='F').copy()
+ elif '_z' in field:
+ data = data[2::3].reshape(grid_dims,order='F').copy()
+ f.close()
+ if grid.pf.field_ordering == 1:
+ return data.T
+ else:
+ return data
+
+ def _read_data_slice(self, grid, field, axis, coord):
+ sl = [slice(None), slice(None), slice(None)]
+ sl[axis] = slice(coord, coord + 1)
+ if grid.pf.field_ordering == 1:
+ sl.reverse()
+
+ f = file(grid.filename, 'rb')
+ dtype, offset = grid.hierarchy._field_map[field]
+ grid_ncells = np.prod(grid.ActiveDimensions)
+
+ read_table_offset = get_read_table_offset(f)
+ f.seek(read_table_offset+offset)
+ if dtype == 'scalar':
+ data = np.fromfile(f, dtype='>f4',
+ count=grid_ncells).reshape(grid.ActiveDimensions,order='F')[sl].copy()
+ if dtype == 'vector':
+ data = np.fromfile(f, dtype='>f4', count=3*grid_ncells)
+ if '_x' in field:
+ data = data[0::3].reshape(grid.ActiveDimensions,order='F')[sl].copy()
+ elif '_y' in field:
+ data = data[1::3].reshape(grid.ActiveDimensions,order='F')[sl].copy()
+ elif '_z' in field:
+ data = data[2::3].reshape(grid.ActiveDimensions,order='F')[sl].copy()
+ f.close()
+ return data
+
+def get_read_table_offset(f):
+ line = f.readline()
+ while True:
+ splitup = line.strip().split()
+ if 'CELL_DATA' in splitup:
+ f.readline()
+ read_table_offset = f.tell()
+ break
+ line = f.readline()
+ return read_table_offset
+
+
diff -r 3e35b1adb541c3f1412d3381d57bb5c0f83e3bf2 -r 0c844cb7df71949206ff497c57ef437d6429d58e yt/frontends/athena/setup.py
--- /dev/null
+++ b/yt/frontends/athena/setup.py
@@ -0,0 +1,13 @@
+#!/usr/bin/env python
+import setuptools
+import os
+import sys
+import os.path
+
+
+def configuration(parent_package='', top_path=None):
+ from numpy.distutils.misc_util import Configuration
+ config = Configuration('athena', parent_package, top_path)
+ config.make_config_py() # installs __config__.py
+ #config.make_svn_version_py()
+ return config
diff -r 3e35b1adb541c3f1412d3381d57bb5c0f83e3bf2 -r 0c844cb7df71949206ff497c57ef437d6429d58e yt/frontends/setup.py
--- a/yt/frontends/setup.py
+++ b/yt/frontends/setup.py
@@ -7,6 +7,7 @@
config = Configuration('frontends', parent_package, top_path)
config.make_config_py() # installs __config__.py
#config.make_svn_version_py()
+ config.add_subpackage("athena")
config.add_subpackage("gdf")
config.add_subpackage("chombo")
config.add_subpackage("enzo")
diff -r 3e35b1adb541c3f1412d3381d57bb5c0f83e3bf2 -r 0c844cb7df71949206ff497c57ef437d6429d58e yt/mods.py
--- a/yt/mods.py
+++ b/yt/mods.py
@@ -95,6 +95,9 @@
from yt.frontends.gdf.api import \
GDFStaticOutput, GDFFieldInfo, add_gdf_field
+from yt.frontends.athena.api import \
+ AthenaStaticOutput, AthenaFieldInfo, add_athena_field
+
from yt.frontends.art.api import \
ARTStaticOutput, ARTFieldInfo, add_art_field
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list