[yt-svn] commit/yt: 11 new changesets
Bitbucket
commits-noreply at bitbucket.org
Mon Mar 12 22:52:13 PDT 2012
11 new commits in yt:
https://bitbucket.org/yt_analysis/yt/changeset/b8a7ceffd281/
changeset: b8a7ceffd281
branch: yt
user: Andrew Myers
date: 2011-07-06 08:49:40
summary: read in radiation-energy-density
affected #: 3 files
diff -r 64cbb6e520d36d2b147dc1677efd5c69426aa39f -r b8a7ceffd281cf614294677383c96e34c8f2a566 yt/frontends/chombo/data_structures.py
--- a/yt/frontends/chombo/data_structures.py
+++ b/yt/frontends/chombo/data_structures.py
@@ -55,7 +55,9 @@
from yt.utilities.parallel_tools.parallel_analysis_interface import \
parallel_root_only
-from .fields import ChomboFieldContainer
+from .fields import \
+ ChomboFieldContainer, \
+ add_field
class ChomboGrid(AMRGridPatch):
_id_offset = 0
@@ -100,11 +102,9 @@
self.hierarchy = os.path.abspath(self.hierarchy_filename)
self.directory = os.path.dirname(self.hierarchy_filename)
self._fhandle = h5py.File(self.hierarchy_filename)
-
self.float_type = self._fhandle['/level_0']['data:datatype=0'].dtype.name
self._levels = self._fhandle.listnames()[1:]
AMRHierarchy.__init__(self,pf,data_style)
-
self._fhandle.close()
def _initialize_data_storage(self):
@@ -163,7 +163,7 @@
def _setup_unknown_fields(self):
pass
-
+
def _setup_derived_fields(self):
self.derived_field_list = []
@@ -179,8 +179,8 @@
def __init__(self, filename, data_style='chombo_hdf5',
storage_filename = None, ini_filename = None):
- # hardcoded for now
- self.current_time = 0.0
+ fileh = h5py.File(filename,'r')
+ self.current_time = fileh.attrs['time']
self.ini_filename = ini_filename
StaticOutput.__init__(self,filename,data_style)
self.storage_filename = storage_filename
@@ -307,7 +307,6 @@
pass
return False
-
@parallel_root_only
def print_key_parameters(self):
for a in ["current_time", "domain_dimensions", "domain_left_edge",
@@ -317,3 +316,13 @@
continue
v = getattr(self, a)
mylog.info("Parameters: %-25s = %s", a, v)
+ if hasattr(self, "cosmological_simulation") and \
+ getattr(self, "cosmological_simulation"):
+ for a in ["current_redshift", "omega_lambda", "omega_matter",
+ "hubble_constant"]:
+ if not hasattr(self, a):
+ mylog.error("Missing %s in parameter file definition!", a)
+ continue
+ v = getattr(self, a)
+ mylog.info("Parameters: %-25s = %s", a, v)
+
diff -r 64cbb6e520d36d2b147dc1677efd5c69426aa39f -r b8a7ceffd281cf614294677383c96e34c8f2a566 yt/frontends/chombo/fields.py
--- a/yt/frontends/chombo/fields.py
+++ b/yt/frontends/chombo/fields.py
@@ -43,7 +43,6 @@
add_field("density", function=lambda a,b: None, take_log=True,
validators = [ValidateDataField("density")],
units=r"\rm{g}/\rm{cm}^3")
-
ChomboFieldInfo["density"]._projected_units =r"\rm{g}/\rm{cm}^2"
add_field("X-momentum", function=lambda a,b: None, take_log=False,
@@ -76,6 +75,16 @@
units=r"",display_name=r"B_z")
ChomboFieldInfo["Z-magnfield"]._projected_units=r""
+add_field("energy-density", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("energy-density")],
+ units=r"\rm{erg}/\rm{cm}^3")
+ChomboFieldInfo["energy-density"]._projected_units =r""
+
+add_field("radiation-energy-density", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("radiation-energy-density")],
+ units=r"\rm{erg}/\rm{cm}^3")
+ChomboFieldInfo["radiation-energy-density"]._projected_units =r""
+
def _MagneticEnergy(field,data):
return (data["X-magnfield"]**2 +
data["Y-magnfield"]**2 +
@@ -96,9 +105,6 @@
"""generate y-velocity from y-momentum and density
"""
- #try:
- # return data["xvel"]
- #except KeyError:
return data["Y-momentum"]/data["density"]
add_field("y-velocity",function=_yVelocity, take_log=False,
units=r'\rm{cm}/\rm{s}')
diff -r 64cbb6e520d36d2b147dc1677efd5c69426aa39f -r b8a7ceffd281cf614294677383c96e34c8f2a566 yt/frontends/chombo/io.py
--- a/yt/frontends/chombo/io.py
+++ b/yt/frontends/chombo/io.py
@@ -45,7 +45,7 @@
fhandle = h5py.File(grid.filename,'r')
ncomp = int(fhandle['/'].attrs['num_components'])
- return [c[1] for c in f['/'].attrs.listitems()[-ncomp:]]
+ return [c[1] for c in f['/'].attrs.listitems()[-ncomp-1:-1]]
def _read_data_set(self,grid,field):
fhandle = h5py.File(grid.hierarchy.hierarchy_filename,'r')
https://bitbucket.org/yt_analysis/yt/changeset/0d88c527ebc8/
changeset: 0d88c527ebc8
branch: yt
user: Andrew Myers
date: 2011-10-10 20:48:06
summary: a fix from Andy, so that the Chombo reader does not assume that the domain starts at [0, 0, 0]
affected #: 1 file
diff -r b8a7ceffd281cf614294677383c96e34c8f2a566 -r 0d88c527ebc8423b31cb454b9eb038d9461c56fc yt/frontends/chombo/data_structures.py
--- a/yt/frontends/chombo/data_structures.py
+++ b/yt/frontends/chombo/data_structures.py
@@ -1,328 +1,346 @@
-"""
-Data structures for Chombo.
-
-Author: Matthew Turk <matthewturk at gmail.com>
-Author: J. S. Oishi <jsoishi at gmail.com>
-Affiliation: KIPAC/SLAC/Stanford
-Homepage: http://yt.enzotools.org/
-License:
- Copyright (C) 2008-2010 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 re
-import os
-import weakref
-import numpy as na
-
-from collections import \
- defaultdict
-from string import \
- strip, \
- rstrip
-from stat import \
- ST_CTIME
-
-from .definitions import \
- pluto2enzoDict, \
- yt2plutoFieldsDict, \
- parameterDict \
-
-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
-from yt.utilities.parallel_tools.parallel_analysis_interface import \
- parallel_root_only
-
-from .fields import \
- ChomboFieldContainer, \
- add_field
-
-class ChomboGrid(AMRGridPatch):
- _id_offset = 0
- __slots__ = ["_level_id", "stop_index"]
- def __init__(self, id, hierarchy, level, start, stop):
- AMRGridPatch.__init__(self, id, filename = hierarchy.hierarchy_filename,
- hierarchy = hierarchy)
- self.Parent = []
- self.Children = []
- self.Level = level
- self.start_index = start.copy()#.transpose()
- self.stop_index = stop.copy()#.transpose()
- self.ActiveDimensions = stop - start + 1
-
- 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.data['dx'], self.data['dy'], self.data['dz'] = self.dds
-
-class ChomboHierarchy(AMRHierarchy):
-
- grid = ChomboGrid
-
- def __init__(self,pf,data_style='chombo_hdf5'):
- self.domain_left_edge = pf.domain_left_edge # need these to determine absolute grid locations
- self.domain_right_edge = pf.domain_right_edge # need these to determine absolute grid locations
- self.data_style = data_style
- self.field_info = ChomboFieldContainer()
- self.field_indexes = {}
- self.parameter_file = weakref.proxy(pf)
- # for now, the hierarchy file is the parameter file!
- self.hierarchy_filename = self.parameter_file.parameter_filename
- self.hierarchy = os.path.abspath(self.hierarchy_filename)
- self.directory = os.path.dirname(self.hierarchy_filename)
- self._fhandle = h5py.File(self.hierarchy_filename)
- self.float_type = self._fhandle['/level_0']['data:datatype=0'].dtype.name
- self._levels = self._fhandle.listnames()[1:]
- AMRHierarchy.__init__(self,pf,data_style)
- self._fhandle.close()
-
- def _initialize_data_storage(self):
- pass
-
- def _detect_fields(self):
- ncomp = int(self._fhandle['/'].attrs['num_components'])
- self.field_list = [c[1] for c in self._fhandle['/'].attrs.listitems()[-ncomp:]]
-
- 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 = 0
- for lev in self._levels:
- self.num_grids += self._fhandle[lev]['Processors'].len()
-
- def _parse_hierarchy(self):
- f = self._fhandle # shortcut
-
- # this relies on the first Group in the H5 file being
- # 'Chombo_global'
- levels = f.listnames()[1:]
- self.grids = []
- i = 0
- for lev in levels:
- level_number = int(re.match('level_(\d+)',lev).groups()[0])
- boxes = f[lev]['boxes'].value
- dx = f[lev].attrs['dx']
- for level_id, box in enumerate(boxes):
- si = na.array([box['lo_%s' % ax] for ax in 'ijk'])
- ei = na.array([box['hi_%s' % ax] for ax in 'ijk'])
- pg = self.grid(len(self.grids),self,level=level_number,
- start = si, stop = ei)
- self.grids.append(pg)
- self.grids[-1]._level_id = level_id
- self.grid_left_edge[i] = dx*si.astype(self.float_type) + self.domain_left_edge
- self.grid_right_edge[i] = dx*(ei.astype(self.float_type)+1) + self.domain_left_edge
- self.grid_particle_count[i] = 0
- self.grid_dimensions[i] = ei - si + 1
- i += 1
- self.grids = na.array(self.grids, dtype='object')
-
- 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_unknown_fields(self):
- pass
-
- 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 ChomboStaticOutput(StaticOutput):
- _hierarchy_class = ChomboHierarchy
- _fieldinfo_class = ChomboFieldContainer
-
- def __init__(self, filename, data_style='chombo_hdf5',
- storage_filename = None, ini_filename = None):
- fileh = h5py.File(filename,'r')
- self.current_time = fileh.attrs['time']
- self.ini_filename = ini_filename
- StaticOutput.__init__(self,filename,data_style)
- self.storage_filename = storage_filename
- self.field_info = self._fieldinfo_class()
-
- 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()
- seconds = 1 #self["Time"]
- self.time_units['years'] = seconds / (365*3600*24.0)
- self.time_units['days'] = seconds / (3600*24.0)
- for key in yt2plutoFieldsDict:
- self.conversion_factors[key] = 1.0
-
- def _setup_nounits_units(self):
- z = 0
- mylog.warning("Setting 1.0 in code units to be 1.0 cm")
- if not self.has_key("TimeUnits"):
- mylog.warning("No time units. Setting 1.0 = 1 second.")
- self.conversion_factors["Time"] = 1.0
- for unit in mpc_conversion.keys():
- self.units[unit] = mpc_conversion[unit] / mpc_conversion["cm"]
-
-
- def _localize(self, f, default):
- if f is None:
- return os.path.join(self.directory, default)
- return f
-
- def _parse_parameter_file(self):
- """
- Check to see whether a 'pluto.ini' or 'orion2.ini' file
- exists in the plot file directory. If one does, attempt to parse it.
- Otherwise, assume the left edge starts at 0 and get the right edge
- from the hdf5 file.
- """
- if os.path.isfile('pluto.ini'):
- self._parse_pluto_file('pluto.ini')
- elif os.path.isfile('orion2.ini'):
- self._parse_pluto_file('orion2.ini')
- else:
- self.unique_identifier = \
- int(os.stat(self.parameter_filename)[ST_CTIME])
- self.domain_left_edge = na.array([0.,0.,0.])
- self.domain_right_edge = self.__calc_right_edge()
- self.dimensionality = 3
- self.refine_by = 2
-
- def _parse_pluto_file(self, ini_filename):
- """
- Reads in an inputs file in the 'pluto.ini' format. Probably not
- especially robust at the moment.
- """
- self.fullplotdir = os.path.abspath(self.parameter_filename)
- self.ini_filename = self._localize( \
- self.ini_filename, ini_filename)
- self.unique_identifier = \
- int(os.stat(self.parameter_filename)[ST_CTIME])
- lines = open(self.ini_filename).readlines()
- # read the file line by line, storing important parameters
- for lineI, line in enumerate(lines):
- try:
- param, sep, vals = map(rstrip,line.partition(' '))
- except ValueError:
- mylog.error("ValueError: '%s'", line)
- if pluto2enzoDict.has_key(param):
- paramName = pluto2enzoDict[param]
- t = map(parameterDict[paramName], vals.split())
- if len(t) == 1:
- self.parameters[paramName] = t[0]
- else:
- if paramName == "RefineBy":
- self.parameters[paramName] = t[0]
- else:
- self.parameters[paramName] = t
-
- # assumes 3D for now
- elif param.startswith("X1-grid"):
- t = vals.split()
- low1 = float(t[1])
- high1 = float(t[4])
- N1 = int(t[2])
- elif param.startswith("X2-grid"):
- t = vals.split()
- low2 = float(t[1])
- high2 = float(t[4])
- N2 = int(t[2])
- elif param.startswith("X3-grid"):
- t = vals.split()
- low3 = float(t[1])
- high3 = float(t[4])
- N3 = int(t[2])
-
- self.dimensionality = 3
- self.domain_left_edge = na.array([low1,low2,low3])
- self.domain_right_edge = na.array([high1,high2,high3])
- self.domain_dimensions = na.array([N1,N2,N3])
- self.refine_by = self.parameters["RefineBy"]
-
- def __calc_right_edge(self):
- fileh = h5py.File(self.parameter_filename,'r')
- dx0 = fileh['/level_0'].attrs['dx']
- RE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain']))[3:] + 1)
- fileh.close()
- return RE
-
- @classmethod
- def _is_valid(self, *args, **kwargs):
- try:
- fileh = h5py.File(args[0],'r')
- if (fileh.listnames())[0] == 'Chombo_global':
- return True
- except:
- pass
- return False
-
- @parallel_root_only
- def print_key_parameters(self):
- for a in ["current_time", "domain_dimensions", "domain_left_edge",
- "domain_right_edge"]:
- if not hasattr(self, a):
- mylog.error("Missing %s in parameter file definition!", a)
- continue
- v = getattr(self, a)
- mylog.info("Parameters: %-25s = %s", a, v)
- if hasattr(self, "cosmological_simulation") and \
- getattr(self, "cosmological_simulation"):
- for a in ["current_redshift", "omega_lambda", "omega_matter",
- "hubble_constant"]:
- if not hasattr(self, a):
- mylog.error("Missing %s in parameter file definition!", a)
- continue
- v = getattr(self, a)
- mylog.info("Parameters: %-25s = %s", a, v)
-
+"""
+Data structures for Chombo.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Author: J. S. Oishi <jsoishi at gmail.com>
+Affiliation: KIPAC/SLAC/Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2008-2010 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 re
+import os
+import weakref
+import numpy as na
+
+from collections import \
+ defaultdict
+from string import \
+ strip, \
+ rstrip
+from stat import \
+ ST_CTIME
+
+from .definitions import \
+ pluto2enzoDict, \
+ yt2plutoFieldsDict, \
+ parameterDict \
+
+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
+from yt.utilities.parallel_tools.parallel_analysis_interface import \
+ parallel_root_only
+
+from .fields import \
+ ChomboFieldContainer, \
+ add_field
+
+class ChomboGrid(AMRGridPatch):
+ _id_offset = 0
+ __slots__ = ["_level_id", "stop_index"]
+ def __init__(self, id, hierarchy, level, start, stop):
+ AMRGridPatch.__init__(self, id, filename = hierarchy.hierarchy_filename,
+ hierarchy = hierarchy)
+ self.Parent = []
+ self.Children = []
+ self.Level = level
+ self.start_index = start.copy()#.transpose()
+ self.stop_index = stop.copy()#.transpose()
+ self.ActiveDimensions = stop - start + 1
+
+ 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.data['dx'], self.data['dy'], self.data['dz'] = self.dds
+
+class ChomboHierarchy(AMRHierarchy):
+
+ grid = ChomboGrid
+
+ def __init__(self,pf,data_style='chombo_hdf5'):
+ self.domain_left_edge = pf.domain_left_edge # need these to determine absolute grid locations
+ self.domain_right_edge = pf.domain_right_edge # need these to determine absolute grid locations
+ self.data_style = data_style
+ self.field_info = ChomboFieldContainer()
+ self.field_indexes = {}
+ self.parameter_file = weakref.proxy(pf)
+ # for now, the hierarchy file is the parameter file!
+ self.hierarchy_filename = self.parameter_file.parameter_filename
+ self.hierarchy = os.path.abspath(self.hierarchy_filename)
+ self.directory = os.path.dirname(self.hierarchy_filename)
+ self._fhandle = h5py.File(self.hierarchy_filename)
+ self.float_type = self._fhandle['/level_0']['data:datatype=0'].dtype.name
+ self._levels = self._fhandle.listnames()[1:]
+ AMRHierarchy.__init__(self,pf,data_style)
+ self._fhandle.close()
+
+ def _initialize_data_storage(self):
+ pass
+
+ def _detect_fields(self):
+ ncomp = int(self._fhandle['/'].attrs['num_components'])
+ self.field_list = [c[1] for c in self._fhandle['/'].attrs.listitems()[-ncomp:]]
+
+ 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 = 0
+ for lev in self._levels:
+ self.num_grids += self._fhandle[lev]['Processors'].len()
+
+ def _parse_hierarchy(self):
+ f = self._fhandle # shortcut
+
+ # this relies on the first Group in the H5 file being
+ # 'Chombo_global'
+ levels = f.listnames()[1:]
+ self.grids = []
+ i = 0
+ for lev in levels:
+ level_number = int(re.match('level_(\d+)',lev).groups()[0])
+ boxes = f[lev]['boxes'].value
+ dx = f[lev].attrs['dx']
+ for level_id, box in enumerate(boxes):
+ si = na.array([box['lo_%s' % ax] for ax in 'ijk'])
+ ei = na.array([box['hi_%s' % ax] for ax in 'ijk'])
+ pg = self.grid(len(self.grids),self,level=level_number,
+ start = si, stop = ei)
+ self.grids.append(pg)
+ self.grids[-1]._level_id = level_id
+ self.grid_left_edge[i] = dx*si.astype(self.float_type) #AJC + self.domain_left_edge
+ self.grid_right_edge[i] = dx*(ei.astype(self.float_type)+1) #AJC + self.domain_left_edge
+ self.grid_particle_count[i] = 0
+ self.grid_dimensions[i] = ei - si + 1
+ i += 1
+ self.grids = na.array(self.grids, dtype='object')
+
+ 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_unknown_fields(self):
+ pass
+
+ 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 ChomboStaticOutput(StaticOutput):
+ _hierarchy_class = ChomboHierarchy
+ _fieldinfo_class = ChomboFieldContainer
+
+ def __init__(self, filename, data_style='chombo_hdf5',
+ storage_filename = None, ini_filename = None):
+ fileh = h5py.File(filename,'r')
+ self.current_time = fileh.attrs['time']
+ self.ini_filename = ini_filename
+ StaticOutput.__init__(self,filename,data_style)
+ self.storage_filename = storage_filename
+ self.field_info = self._fieldinfo_class()
+
+ 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()
+ seconds = 1 #self["Time"]
+ self.time_units['years'] = seconds / (365*3600*24.0)
+ self.time_units['days'] = seconds / (3600*24.0)
+ for key in yt2plutoFieldsDict:
+ self.conversion_factors[key] = 1.0
+
+ def _setup_nounits_units(self):
+ z = 0
+ mylog.warning("Setting 1.0 in code units to be 1.0 cm")
+ if not self.has_key("TimeUnits"):
+ mylog.warning("No time units. Setting 1.0 = 1 second.")
+ self.conversion_factors["Time"] = 1.0
+ for unit in mpc_conversion.keys():
+ self.units[unit] = mpc_conversion[unit] / mpc_conversion["cm"]
+
+
+ def _localize(self, f, default):
+ if f is None:
+ return os.path.join(self.directory, default)
+ return f
+
+ def _parse_parameter_file(self):
+ """
+ Check to see whether a 'pluto.ini' or 'orion2.ini' file
+ exists in the plot file directory. If one does, attempt to parse it.
+ Otherwise, assume the left edge starts at 0 and get the right edge
+ from the hdf5 file.
+ """
+ if os.path.isfile('pluto.ini'):
+ self._parse_pluto_file('pluto.ini')
+ elif os.path.isfile('orion2.ini'):
+ # AJC - this shouldn't be nesecary anymore, now
+ # that orion2 has fixed the pluto domian-offset glitch
+ self._parse_pluto_file('orion2.ini')
+ else:
+ self.unique_identifier = \
+ int(os.stat(self.parameter_filename)[ST_CTIME])
+ self.domain_left_edge = self.__calc_left_edge()
+ self.domain_right_edge = self.__calc_right_edge()
+ self.dimensionality = 3
+ # AJC - make refine_by a list because different levels
+ # can have different refinment ratios.
+ self.refine_by = []
+ fileh = h5py.File(self.parameter_filename,'r')
+ for level in range(0,fileh.attrs['max_level']):
+ self.refine_by.append(fileh['/level_'+str(level)].attrs['ref_ratio'])
+
+ def _parse_pluto_file(self, ini_filename):
+ """
+ Reads in an inputs file in the 'pluto.ini' format. Probably not
+ especially robust at the moment.
+ """
+ self.fullplotdir = os.path.abspath(self.parameter_filename)
+ self.ini_filename = self._localize( \
+ self.ini_filename, ini_filename)
+ self.unique_identifier = \
+ int(os.stat(self.parameter_filename)[ST_CTIME])
+ lines = open(self.ini_filename).readlines()
+ # read the file line by line, storing important parameters
+ for lineI, line in enumerate(lines):
+ try:
+ param, sep, vals = map(rstrip,line.partition(' '))
+ except ValueError:
+ mylog.error("ValueError: '%s'", line)
+ if pluto2enzoDict.has_key(param):
+ paramName = pluto2enzoDict[param]
+ t = map(parameterDict[paramName], vals.split())
+ if len(t) == 1:
+ self.parameters[paramName] = t[0]
+ else:
+ if paramName == "RefineBy":
+ self.parameters[paramName] = t[0]
+ else:
+ self.parameters[paramName] = t
+
+ # assumes 3D for now
+ elif param.startswith("X1-grid"):
+ t = vals.split()
+ low1 = float(t[1])
+ high1 = float(t[4])
+ N1 = int(t[2])
+ elif param.startswith("X2-grid"):
+ t = vals.split()
+ low2 = float(t[1])
+ high2 = float(t[4])
+ N2 = int(t[2])
+ elif param.startswith("X3-grid"):
+ t = vals.split()
+ low3 = float(t[1])
+ high3 = float(t[4])
+ N3 = int(t[2])
+
+ self.dimensionality = 3
+ self.domain_left_edge = na.array([low1,low2,low3])
+ self.domain_right_edge = na.array([high1,high2,high3])
+ self.domain_dimensions = na.array([N1,N2,N3])
+ self.refine_by = self.parameters["RefineBy"]
+
+ # AJC - added this routine
+ # Note that vanilla Chombo files (independent of pluto or orion2)
+ # do not always have left_edge = [0,0,0]
+ def __calc_left_edge(self):
+ fileh = h5py.File(self.parameter_filename,'r')
+ dx0 = fileh['/level_0'].attrs['dx']
+ RE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain']))[0:3])
+ fileh.close()
+ return RE
+
+
+ def __calc_right_edge(self):
+ fileh = h5py.File(self.parameter_filename,'r')
+ dx0 = fileh['/level_0'].attrs['dx']
+ RE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain']))[3:] + 1)
+ fileh.close()
+ return RE
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ try:
+ fileh = h5py.File(args[0],'r')
+ if (fileh.listnames())[0] == 'Chombo_global':
+ return True
+ except:
+ pass
+ return False
+
+ @parallel_root_only
+ def print_key_parameters(self):
+ for a in ["current_time", "domain_dimensions", "domain_left_edge",
+ "domain_right_edge"]:
+ if not hasattr(self, a):
+ mylog.error("Missing %s in parameter file definition!", a)
+ continue
+ v = getattr(self, a)
+ mylog.info("Parameters: %-25s = %s", a, v)
+ if hasattr(self, "cosmological_simulation") and \
+ getattr(self, "cosmological_simulation"):
+ for a in ["current_redshift", "omega_lambda", "omega_matter",
+ "hubble_constant"]:
+ if not hasattr(self, a):
+ mylog.error("Missing %s in parameter file definition!", a)
+ continue
+ v = getattr(self, a)
+ mylog.info("Parameters: %-25s = %s", a, v)
+
https://bitbucket.org/yt_analysis/yt/changeset/2922825feafe/
changeset: 2922825feafe
branch: yt
user: Andrew Myers
date: 2011-10-10 21:07:23
summary: a fix to the Ramses reader by Andy Cunningham. Need for the code to work on his machine.
affected #: 1 file
diff -r 0d88c527ebc8423b31cb454b9eb038d9461c56fc -r 2922825feafeb209651d127da1fc8ee94f39f907 yt/frontends/ramses/_ramses_reader.pyx
--- a/yt/frontends/ramses/_ramses_reader.pyx
+++ b/yt/frontends/ramses/_ramses_reader.pyx
@@ -1,952 +1,952 @@
-"""
-Wrapping code for Oliver Hahn's RamsesRead++
-
-Author: Matthew Turk <matthewturk at gmail.com>
-Affiliation: UCSD
-Author: Oliver Hahn <ohahn at stanford.edu>
-Affiliation: KIPAC / Stanford
-Homepage: http://yt.enzotools.org/
-License:
- Copyright (C) 2010 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/>.
-"""
-
-# Cython wrapping code for Oliver Hahn's RAMSES reader
-from cython.operator cimport dereference as deref, preincrement as inc
-from libc.stdlib cimport malloc, free, abs, calloc, labs
-
-import numpy as np
-cimport numpy as np
-cimport cython
-
-cdef extern from "math.h":
- double ceil(double x)
- double log2(double x)
-
-cdef inline np.int64_t i64max(np.int64_t i0, np.int64_t i1):
- if i0 > i1: return i0
- return i1
-
-cdef inline np.int64_t i64min(np.int64_t i0, np.int64_t i1):
- if i0 < i1: return i0
- return i1
-
-cdef extern from "<vector>" namespace "std":
- cdef cppclass vector[T]:
- cppclass iterator:
- T operator*()
- iterator operator++()
- bint operator==(iterator)
- bint operator!=(iterator)
- vector()
- void push_back(T&)
- T& operator[](int)
- T& at(int)
- iterator begin()
- iterator end()
-
-cdef extern from "<map>" namespace "std":
- cdef cppclass map[A,B]:
- pass
-
-cdef extern from "string" namespace "std":
- cdef cppclass string:
- string(char *cstr)
- char *c_str()
- string operator*()
-
-cdef extern from "RAMSES_typedefs.h":
- pass
-
-cdef extern from "RAMSES_info.hh" namespace "RAMSES":
- enum codeversion:
- version1
- version2
- version3
-
- cdef cppclass snapshot:
- string m_filename
- codeversion m_version
- cppclass info_data:
- unsigned ncpu
- unsigned ndim
- unsigned levelmin
- unsigned levelmax
- unsigned ngridmax
- unsigned nstep_coarse
- double boxlen
- double time
- double aexp
- double H0
- double omega_m
- double omega_l
- double omega_k
- double omega_b
- double unit_l
- double unit_d
- double unit_t
-
- info_data m_header
- vector[double] ind_min
- vector[double] ind_max
-
- snapshot(string info_filename, codeversion ver)
-
- unsigned get_snapshot_num()
- unsigned getdomain_bykey( double key )
-
- #void hilbert3d(vector[double] &x, vector[double] &y, vector[double] &z,
- # vector[double] &order, unsigned bit_length)
-
-cdef extern from "RAMSES_amr_data.hh" namespace "RAMSES::AMR":
- cdef cppclass vec[real_t]:
- real_t x, y, z
- vec( real_t x_, real_t y_, real_t z_)
- vec ( vec& v)
- vec ( )
-
- # This class definition is out of date. I have unrolled the template
- # below.
- cdef cppclass cell_locally_essential[id_t, real_t]:
- id_t m_neighbor[6]
- id_t m_father
- id_t m_son[8]
- real_t m_xg[3]
- id_t m_cpu
-
- char m_pos
-
- cell_locally_essential()
-
- bint is_refined(int ison)
-
- cdef cppclass RAMSES_cell:
- unsigned m_neighbour[6]
- unsigned m_father
- unsigned m_son[8]
- float m_xg[3]
- unsigned m_cpu
-
- char m_pos
-
- RAMSES_cell()
-
- bint is_refined(int ison)
-
- cdef cppclass cell_simple[id_t, real_t]:
- id_t m_son[1]
- real_t m_xg[3]
- id_t m_cpu
-
- char m_pos
- cell_simple()
- bint is_refined(int ison)
-
- # AMR level implementation
-
- # This class definition is out of date. I have unrolled the template
- # below.
- cdef cppclass level[Cell_]:
- unsigned m_ilevel
- vector[Cell_] m_level_cells
- double m_xc[8]
- double m_yc[8]
- double m_zc[8]
-
- # I skipped the typedefs here
-
- double m_dx
- unsigned m_nx
- level (unsigned ilevel)
-
- void register_cell( Cell_ cell )
- vector[Cell_].iterator begin()
- vector[Cell_].iterator end()
-
- Cell_& operator[]( unsigned i)
- unsigned size()
-
- cdef cppclass RAMSES_level:
- unsigned m_ilevel
- vector[RAMSES_cell] m_level_cells
- double m_xc[8]
- double m_yc[8]
- double m_zc[8]
-
- # I skipped the typedefs here
-
- double m_dx
- unsigned m_nx
- RAMSES_level (unsigned ilevel)
-
- void register_cell( RAMSES_cell cell )
- vector[RAMSES_cell].iterator begin()
- vector[RAMSES_cell].iterator end()
-
- RAMSES_cell& operator[]( unsigned i)
- unsigned size()
-
- # This class definition is out of date. I have unrolled the template
- # below.
- cdef cppclass tree[Cell_, Level_]:
- cppclass header:
- vector[int] nx
- vector[int] nout
- vector[int] nsteps
- int ncpu
- int ndim
- int nlevelmax
- int ngridmax
- int nboundary
- int ngrid_current
- double boxlen
- vector[double] tout
- vector[double] aout
- vector[double] dtold
- vector[double] dtnew
- vector[double] cosm
- vector[double] timing
- double t
- double mass_sph
-
- vector[Level_] m_AMR_levels
- vector[unsigned] mheadl, m_numbl, m_taill
-
- int m_cpu
- int m_minlevel
- int m_maxlevel
- string m_fname
- unsigned m_ncoarse
- header m_header
-
- # This is from later on in the .hh file ... I don't think we need them
- # typedef proto_iterator<const tree*> const_iterator;
- # typedef proto_iterator<tree *> iterator;
-
- tree (snapshot &snap, int cpu, int maxlevel, int minlevel = 1)
-
- cppclass tree_iterator "RAMSES::AMR::RAMSES_tree::iterator":
- tree_iterator operator*()
- RAMSES_cell operator*()
- tree_iterator begin()
- tree_iterator end()
- tree_iterator to_parent()
- tree_iterator get_parent()
- void next()
- bint operator!=(tree_iterator other)
- unsigned get_cell_father()
- bint is_finest(int ison)
- int get_absolute_position()
- int get_domain()
-
- cdef cppclass RAMSES_tree:
- # This is, strictly speaking, not a header. But, I believe it is
- # going to work alright.
- cppclass header:
- vector[int] nx
- vector[int] nout
- vector[int] nsteps
- int ncpu
- int ndim
- int nlevelmax
- int ngridmax
- int nboundary
- int ngrid_current
- double boxlen
- vector[double] tout
- vector[double] aout
- vector[double] dtold
- vector[double] dtnew
- vector[double] cosm
- vector[double] timing
- double t
- double mass_sph
-
- vector[RAMSES_level] m_AMR_levels
- vector[unsigned] mheadl, m_numbl, m_taill
-
- int m_cpu
- int m_minlevel
- int m_maxlevel
- string m_fname
- unsigned m_ncoarse
- header m_header
-
- unsigned size()
-
- # This is from later on in the .hh file ... I don't think we need them
- # typedef proto_iterator<const tree*> const_iterator;
- # typedef proto_iterator<tree *> iterator;
-
- RAMSES_tree(snapshot &snap, int cpu, int maxlevel, int minlevel)
- void read()
-
- tree_iterator begin(int ilevel)
- tree_iterator end(int ilevel)
-
- tree_iterator begin()
- tree_iterator end()
-
- vec[double] cell_pos_double "cell_pos<double>" (tree_iterator it, unsigned ind)
- vec[double] grid_pos_double "grid_pos<double>" (tree_iterator it)
- vec[float] cell_pos_float "cell_pos<float>" (tree_iterator it, unsigned ind)
- vec[float] grid_pos_float "grid_pos<float>" (tree_iterator it)
-
-cdef extern from "RAMSES_amr_data.hh" namespace "RAMSES::HYDRO":
- enum hydro_var:
- density
- velocity_x
- velocity_y
- velocity_z
- pressure
- metallicit
-
- char ramses_hydro_variables[][64]
-
- # There are a number of classes here that we could wrap and utilize.
- # However, I will only wrap the methods I know we need.
-
- # I have no idea if this will work.
- cdef cppclass TreeTypeIterator[TreeType_]:
- pass
-
- # This class definition is out of date. I have unrolled the template
- # below.
- cdef cppclass data[TreeType_, Real_]:
- cppclass header:
- unsigned ncpu
- unsigned nvar
- unsigned ndim
- unsigned nlevelmax
- unsigned nboundary
- double gamma
- string m_fname
- header m_header
-
- # I don't want to implement proto_data, so we put this here
- Real_& cell_value(TreeTypeIterator[TreeType_] &it, int ind)
-
- unsigned m_nvars
- vector[string] m_varnames
- map[string, unsigned] m_var_name_map
-
- data(TreeType_ &AMRtree)
- #_OutputIterator get_var_names[_OutputIterator](_OutputIterator names)
- void read(string varname)
-
- cdef cppclass RAMSES_hydro_data:
- cppclass header:
- unsigned ncpu
- unsigned nvar
- unsigned ndim
- unsigned nlevelmax
- unsigned nboundary
- double gamma
- string m_fname
- header m_header
-
- # I don't want to implement proto_data, so we put this here
- double cell_value (tree_iterator &it, int ind)
-
- unsigned m_nvars
- vector[string] m_varnames
- map[string, unsigned] m_var_name_map
-
- RAMSES_hydro_data(RAMSES_tree &AMRtree)
- #_OutputIterator get_var_names[_OutputIterator](_OutputIterator names)
- void read(string varname)
-
- vector[vector[double]] m_var_array
-
-cdef class RAMSES_tree_proxy:
- cdef string *snapshot_name
- cdef snapshot *rsnap
-
- cdef RAMSES_tree** trees
- cdef RAMSES_hydro_data*** hydro_datas
-
- cdef int **loaded
-
- cdef public object field_ind
- cdef public object field_names
-
- # We will store this here so that we have a record, independent of the
- # header, of how many things we have allocated
- cdef int ndomains, nfields
-
- def __cinit__(self, char *fn):
- cdef int idomain, ifield, ii
- cdef RAMSES_tree *local_tree
- cdef RAMSES_hydro_data *local_hydro_data
- self.snapshot_name = new string(fn)
- self.rsnap = new snapshot(deref(self.snapshot_name), version3)
- # We now have to get our field names to fill our array
- self.trees = <RAMSES_tree**>\
- malloc(sizeof(RAMSES_tree*) * self.rsnap.m_header.ncpu)
- self.hydro_datas = <RAMSES_hydro_data ***>\
- malloc(sizeof(RAMSES_hydro_data**) * self.rsnap.m_header.ncpu)
- self.ndomains = self.rsnap.m_header.ncpu
- #for ii in range(self.ndomains): self.trees[ii] = NULL
- # Note we don't do ncpu + 1
- for idomain in range(self.rsnap.m_header.ncpu):
- # we don't delete local_tree
- local_tree = new RAMSES_tree(deref(self.rsnap), idomain + 1,
- self.rsnap.m_header.levelmax, 0)
- local_tree.read()
- local_hydro_data = new RAMSES_hydro_data(deref(local_tree))
- self.hydro_datas[idomain] = <RAMSES_hydro_data **>\
- malloc(sizeof(RAMSES_hydro_data*) * local_hydro_data.m_nvars)
- for ii in range(local_hydro_data.m_nvars):
- self.hydro_datas[idomain][ii] = \
- new RAMSES_hydro_data(deref(local_tree))
- self.trees[idomain] = local_tree
- # We do not delete the final snapshot, which we'll use later
- if idomain + 1 < self.rsnap.m_header.ncpu:
- del local_hydro_data
- # Only once, we read all the field names
- self.nfields = local_hydro_data.m_nvars
- cdef string *field_name
- self.field_names = []
- self.field_ind = {}
- self.loaded = <int **> malloc(sizeof(int) * local_hydro_data.m_nvars)
- for idomain in range(self.ndomains):
- self.loaded[idomain] = <int *> malloc(
- sizeof(int) * local_hydro_data.m_nvars)
- for ifield in range(local_hydro_data.m_nvars):
- self.loaded[idomain][ifield] = 0
- for ifield in range(local_hydro_data.m_nvars):
- field_name = &(local_hydro_data.m_varnames[ifield])
- # Does this leak?
- self.field_names.append(field_name.c_str())
- self.field_ind[self.field_names[-1]] = ifield
- # This all needs to be cleaned up in the deallocator
- del local_hydro_data
-
- def __dealloc__(self):
- import traceback; traceback.print_stack()
- cdef int idomain, ifield
- # To ensure that 'delete' is used, not 'free',
- # we allocate temporary variables.
- cdef RAMSES_tree *temp_tree
- cdef RAMSES_hydro_data *temp_hdata
- for idomain in range(self.ndomains):
- for ifield in range(self.nfields):
- temp_hdata = self.hydro_datas[idomain][ifield]
- del temp_hdata
- temp_tree = self.trees[idomain]
- del temp_tree
- free(self.loaded[idomain])
- free(self.hydro_datas[idomain])
- free(self.trees)
- free(self.hydro_datas)
- free(self.loaded)
- if self.snapshot_name != NULL: del self.snapshot_name
- if self.rsnap != NULL: del self.rsnap
-
- def count_zones(self):
- # We need to do simulation domains here
-
- cdef unsigned idomain, ilevel
- cdef RAMSES_tree *local_tree
- cdef RAMSES_hydro_data *local_hydro_data
- cdef RAMSES_level *local_level
- cdef tree_iterator grid_it, grid_end
-
- # All the loop-local pointers must be declared up here
-
- cell_count = []
- cdef int local_count = 0
- for ilevel in range(self.rsnap.m_header.levelmax + 1):
- cell_count.append(0)
- for idomain in range(1, self.rsnap.m_header.ncpu + 1):
- local_tree = new RAMSES_tree(deref(self.rsnap), idomain,
- self.rsnap.m_header.levelmax, 0)
- local_tree.read()
- local_hydro_data = new RAMSES_hydro_data(deref(local_tree))
- for ilevel in range(local_tree.m_maxlevel + 1):
- local_count = 0
- local_level = &local_tree.m_AMR_levels[ilevel]
- grid_it = local_tree.begin(ilevel)
- grid_end = local_tree.end(ilevel)
- while grid_it != grid_end:
- local_count += (grid_it.get_domain() == idomain)
- grid_it.next()
- cell_count[ilevel] += local_count
- del local_tree, local_hydro_data
-
- return cell_count
-
- def ensure_loaded(self, char *varname, int domain_index):
- # this domain_index must be zero-indexed
- cdef int varindex = self.field_ind[varname]
- cdef string *field_name = new string(varname)
- if self.loaded[domain_index][varindex] == 1:
- return
- print "READING FROM DISK", varname, domain_index, varindex
- self.hydro_datas[domain_index][varindex].read(deref(field_name))
- self.loaded[domain_index][varindex] = 1
- del field_name
-
- def clear_tree(self, char *varname, int domain_index):
- # this domain_index must be zero-indexed
- # We delete and re-create
- cdef int varindex = self.field_ind[varname]
- cdef string *field_name = new string(varname)
- if self.loaded[domain_index][varindex] == 0: return
- cdef RAMSES_hydro_data *temp_hdata = self.hydro_datas[domain_index][varindex]
- del temp_hdata
- self.hydro_datas[domain_index - 1][varindex] = \
- new RAMSES_hydro_data(deref(self.trees[domain_index]))
- self.loaded[domain_index][varindex] = 0
- del field_name
-
- def get_file_info(self):
- header_info = {}
- header_info["ncpu"] = self.rsnap.m_header.ncpu
- header_info["ndim"] = self.rsnap.m_header.ndim
- header_info["levelmin"] = self.rsnap.m_header.levelmin
- header_info["levelmax"] = self.rsnap.m_header.levelmax
- header_info["ngridmax"] = self.rsnap.m_header.ngridmax
- header_info["nstep_coarse"] = self.rsnap.m_header.nstep_coarse
- header_info["boxlen"] = self.rsnap.m_header.boxlen
- header_info["time"] = self.rsnap.m_header.time
- header_info["aexp"] = self.rsnap.m_header.aexp
- header_info["H0"] = self.rsnap.m_header.H0
- header_info["omega_m"] = self.rsnap.m_header.omega_m
- header_info["omega_l"] = self.rsnap.m_header.omega_l
- header_info["omega_k"] = self.rsnap.m_header.omega_k
- header_info["omega_b"] = self.rsnap.m_header.omega_b
- header_info["unit_l"] = self.rsnap.m_header.unit_l
- header_info["unit_d"] = self.rsnap.m_header.unit_d
- header_info["unit_t"] = self.rsnap.m_header.unit_t
-
- # Now we grab some from the trees
- cdef np.ndarray[np.int32_t, ndim=1] top_grid_dims = np.zeros(3, "int32")
- cdef int i
- for i in range(3):
- # Note that nx is really boundary conditions. We always have 2.
- top_grid_dims[i] = self.trees[0].m_header.nx[i]
- top_grid_dims[i] = 2
- header_info["nx"] = top_grid_dims
-
- return header_info
-
- def fill_hierarchy_arrays(self,
- np.ndarray[np.int32_t, ndim=1] top_grid_dims,
- np.ndarray[np.float64_t, ndim=2] left_edges,
- np.ndarray[np.float64_t, ndim=2] right_edges,
- np.ndarray[np.int32_t, ndim=2] grid_levels,
- np.ndarray[np.int64_t, ndim=2] grid_file_locations,
- np.ndarray[np.uint64_t, ndim=1] grid_hilbert_indices,
- np.ndarray[np.int32_t, ndim=2] child_mask,
- np.ndarray[np.float64_t, ndim=1] domain_left,
- np.ndarray[np.float64_t, ndim=1] domain_right):
- # We need to do simulation domains here
-
- cdef unsigned idomain, ilevel
-
- # All the loop-local pointers must be declared up here
- cdef RAMSES_tree *local_tree
- cdef RAMSES_hydro_data *local_hydro_data
- cdef unsigned father
-
- cdef tree_iterator grid_it, grid_end, father_it
- cdef vec[double] gvec
- cdef int grid_ind = 0, grid_aind = 0
- cdef unsigned parent_ind
- cdef bint ci
- cdef double pos[3]
- cdef double grid_half_width
- cdef unsigned long rv
-
- cdef np.int32_t rr
- cdef int i
- cell_count = []
- level_cell_counts = {}
- for idomain in range(1, self.rsnap.m_header.ncpu + 1):
- local_tree = new RAMSES_tree(deref(self.rsnap), idomain,
- self.rsnap.m_header.levelmax, 0)
- local_tree.read()
- local_hydro_data = new RAMSES_hydro_data(deref(local_tree))
- for ilevel in range(local_tree.m_maxlevel + 1):
- # this gets overwritten for every domain, which is okay
- level_cell_counts[ilevel] = grid_ind
- #grid_half_width = self.rsnap.m_header.boxlen / \
- grid_half_width = 1.0 / \
- (2**(ilevel) * top_grid_dims[0])
- grid_it = local_tree.begin(ilevel)
- grid_end = local_tree.end(ilevel)
- while grid_it != grid_end:
- if grid_it.get_domain() != idomain:
- grid_ind += 1
- grid_it.next()
- continue
- gvec = local_tree.grid_pos_double(grid_it)
- left_edges[grid_aind, 0] = pos[0] = gvec.x - grid_half_width
- left_edges[grid_aind, 1] = pos[1] = gvec.y - grid_half_width
- left_edges[grid_aind, 2] = pos[2] = gvec.z - grid_half_width
- for i in range(3):
- pos[i] = (pos[i] - domain_left[i]) / (domain_right[i] - domain_left[i])
- right_edges[grid_aind, 0] = gvec.x + grid_half_width
- right_edges[grid_aind, 1] = gvec.y + grid_half_width
- right_edges[grid_aind, 2] = gvec.z + grid_half_width
- grid_levels[grid_aind, 0] = ilevel
- # Now the harder part
- father_it = grid_it.get_parent()
- grid_file_locations[grid_aind, 0] = <np.int64_t> idomain
- grid_file_locations[grid_aind, 1] = grid_ind - level_cell_counts[ilevel]
- parent_ind = father_it.get_absolute_position()
- if ilevel > 0:
- # We calculate the REAL parent index
- grid_file_locations[grid_aind, 2] = \
- level_cell_counts[ilevel - 1] + parent_ind
- else:
- grid_file_locations[grid_aind, 2] = -1
- for ci in range(8):
- rr = <np.int32_t> grid_it.is_finest(ci)
- child_mask[grid_aind, ci] = rr
- grid_ind += 1
- grid_aind += 1
- grid_it.next()
- del local_tree, local_hydro_data
-
- def read_oct_grid(self, char *field, int level, int domain, int grid_id):
-
- self.ensure_loaded(field, domain - 1)
- cdef int varindex = self.field_ind[field]
- cdef int i
-
- cdef np.ndarray[np.float64_t, ndim=3] tr = np.empty((2,2,2), dtype='float64',
- order='F')
- cdef tree_iterator grid_it, grid_end
- cdef double* data = <double*> tr.data
-
- cdef RAMSES_tree *local_tree = self.trees[domain - 1]
- cdef RAMSES_hydro_data *local_hydro_data = self.hydro_datas[domain - 1][varindex]
-
- #inline ValueType_& cell_value( const typename TreeType_::iterator& it,
- # int ind )
- #{
- # unsigned ipos = it.get_absolute_position();
- # unsigned ilevel = it.get_level();//-m_minlevel;
- # return (m_var_array[ilevel])[m_twotondim*ipos+ind];
- #}
-
- for i in range(8):
- data[i] = local_hydro_data.m_var_array[level][8*grid_id+i]
- return tr
-
- def read_grid(self, char *field,
- np.ndarray[np.int64_t, ndim=1] start_index,
- np.ndarray[np.int32_t, ndim=1] grid_dims,
- np.ndarray[np.float64_t, ndim=3] data,
- np.ndarray[np.int32_t, ndim=3] filled,
- int level, int ref_factor,
- component_grid_info):
- cdef int varindex = self.field_ind[field]
- cdef RAMSES_tree *local_tree = NULL
- cdef RAMSES_hydro_data *local_hydro_data = NULL
-
- cdef int gi, i, j, k, domain, offset
- cdef int ir, jr, kr
- cdef int offi, offj, offk, odind
- cdef np.int64_t di, dj, dk
- cdef np.ndarray[np.int64_t, ndim=1] ogrid_info
- cdef np.ndarray[np.int64_t, ndim=1] og_start_index
- cdef np.float64_t temp_data
- cdef np.int64_t end_index[3]
- cdef int to_fill = 0
- # Note that indexing into a cell is:
- # (k*2 + j)*2 + i
- for i in range(3):
- end_index[i] = start_index[i] + grid_dims[i]
- for gi in range(len(component_grid_info)):
- ogrid_info = component_grid_info[gi]
- domain = ogrid_info[0]
- #print "Loading", domain, ogrid_info
- self.ensure_loaded(field, domain - 1)
- local_tree = self.trees[domain - 1]
- local_hydro_data = self.hydro_datas[domain - 1][varindex]
- offset = ogrid_info[1]
- og_start_index = ogrid_info[3:]
- for i in range(2*ref_factor):
- di = i + og_start_index[0] * ref_factor
- if di < start_index[0] or di >= end_index[0]: continue
- ir = <int> (i / ref_factor)
- for j in range(2 * ref_factor):
- dj = j + og_start_index[1] * ref_factor
- if dj < start_index[1] or dj >= end_index[1]: continue
- jr = <int> (j / ref_factor)
- for k in range(2 * ref_factor):
- dk = k + og_start_index[2] * ref_factor
- if dk < start_index[2] or dk >= end_index[2]: continue
- kr = <int> (k / ref_factor)
- offi = di - start_index[0]
- offj = dj - start_index[1]
- offk = dk - start_index[2]
- #print offi, filled.shape[0],
- #print offj, filled.shape[1],
- #print offk, filled.shape[2]
- if filled[offi, offj, offk] == 1: continue
-
- odind = (kr*2 + jr)*2 + ir
- temp_data = local_hydro_data.m_var_array[
- level][8*offset + odind]
- data[offi, offj, offk] = temp_data
- filled[offi, offj, offk] = 1
- to_fill += 1
- return to_fill
-
-cdef class ProtoSubgrid:
- cdef np.int64_t *signature[3]
- cdef np.int64_t left_edge[3]
- cdef np.int64_t right_edge[3]
- cdef np.int64_t dimensions[3]
- cdef public np.float64_t efficiency
- cdef public object sigs
- cdef public object grid_file_locations
- cdef public object dd
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- def __cinit__(self,
- np.ndarray[np.int64_t, ndim=1] left_index,
- np.ndarray[np.int64_t, ndim=1] dimensions,
- np.ndarray[np.int64_t, ndim=2] left_edges,
- np.ndarray[np.int64_t, ndim=2] right_edges,
- np.ndarray[np.int64_t, ndim=2] grid_dimensions,
- np.ndarray[np.int64_t, ndim=2] grid_file_locations):
- # This also includes the shrinking step.
- cdef int i, ci, ng = left_edges.shape[0]
- cdef np.ndarray temp_arr
- cdef int l0, r0, l1, r1, l2, r2, i0, i1, i2
- cdef np.int64_t temp_l[3], temp_r[3], ncells
- cdef np.float64_t efficiency
- self.sigs = []
- for i in range(3):
- temp_l[i] = left_index[i] + dimensions[i]
- temp_r[i] = left_index[i]
- self.signature[i] = NULL
- for gi in range(ng):
- if left_edges[gi,0] > left_index[0]+dimensions[0] or \
- right_edges[gi,0] < left_index[0] or \
- left_edges[gi,1] > left_index[1]+dimensions[1] or \
- right_edges[gi,1] < left_index[1] or \
- left_edges[gi,2] > left_index[2]+dimensions[2] or \
- right_edges[gi,2] < left_index[2]:
- #print "Skipping grid", gi, "which lies outside out box"
- continue
- for i in range(3):
- temp_l[i] = i64min(left_edges[gi,i], temp_l[i])
- temp_r[i] = i64max(right_edges[gi,i], temp_r[i])
- for i in range(3):
- self.left_edge[i] = i64max(temp_l[i], left_index[i])
- self.right_edge[i] = i64min(temp_r[i], left_index[i] + dimensions[i])
- self.dimensions[i] = self.right_edge[i] - self.left_edge[i]
- if self.dimensions[i] <= 0:
- self.efficiency = -1.0
- return
- self.sigs.append(np.zeros(self.dimensions[i], 'int64'))
- #print self.sigs[0].size, self.sigs[1].size, self.sigs[2].size
-
- # My guess is that this whole loop could be done more efficiently.
- # However, this is clear and straightforward, so it is a good first
- # pass.
- cdef np.ndarray[np.int64_t, ndim=1] sig0, sig1, sig2
- sig0 = self.sigs[0]
- sig1 = self.sigs[1]
- sig2 = self.sigs[2]
- efficiency = 0.0
- cdef int used
- self.grid_file_locations = []
- for gi in range(ng):
- used = 0
- nnn = 0
- for l0 in range(grid_dimensions[gi, 0]):
- i0 = left_edges[gi, 0] + l0
- if i0 < self.left_edge[0]: continue
- if i0 >= self.right_edge[0]: break
- for l1 in range(grid_dimensions[gi, 1]):
- i1 = left_edges[gi, 1] + l1
- if i1 < self.left_edge[1]: continue
- if i1 >= self.right_edge[1]: break
- for l2 in range(grid_dimensions[gi, 2]):
- i2 = left_edges[gi, 2] + l2
- if i2 < self.left_edge[2]: continue
- if i2 >= self.right_edge[2]: break
- i = i0 - self.left_edge[0]
- sig0[i] += 1
- i = i1 - self.left_edge[1]
- sig1[i] += 1
- i = i2 - self.left_edge[2]
- sig2[i] += 1
- efficiency += 1
- used = 1
- if used == 1:
- grid_file_locations[gi,3] = left_edges[gi, 0]
- grid_file_locations[gi,4] = left_edges[gi, 1]
- grid_file_locations[gi,5] = left_edges[gi, 2]
- self.grid_file_locations.append(grid_file_locations[gi,:])
-
- self.dd = np.ones(3, dtype='int64')
- for i in range(3):
- efficiency /= self.dimensions[i]
- self.dd[i] = self.dimensions[i]
- #print "Efficiency is %0.3e" % (efficiency)
- self.efficiency = efficiency
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- def find_split(self):
- # First look for zeros
- cdef int i, center, ax
- cdef np.ndarray[ndim=1, dtype=np.int64_t] axes
- cdef np.int64_t strength, zcstrength, zcp
- axes = np.argsort(self.dd)[::-1]
- cdef np.ndarray[np.int64_t] sig
- for axi in range(3):
- ax = axes[axi]
- center = self.dimensions[ax] / 2
- sig = self.sigs[ax]
- for i in range(self.dimensions[ax]):
- if sig[i] == 0 and i > 0 and i < self.dimensions[ax] - 1:
- #print "zero: %s (%s)" % (i, self.dimensions[ax])
- return 0, ax, i
- zcstrength = 0
- zcp = 0
- zca = -1
- cdef int temp
- cdef np.int64_t *sig2d
- for axi in range(3):
- ax = axes[axi]
- sig = self.sigs[ax]
- sig2d = <np.int64_t *> malloc(sizeof(np.int64_t) * self.dimensions[ax])
- sig2d[0] = sig2d[self.dimensions[ax]-1] = 0
- for i in range(1, self.dimensions[ax] - 1):
- sig2d[i] = sig[i-1] - 2*sig[i] + sig[i+1]
- for i in range(1, self.dimensions[ax] - 1):
- if sig2d[i] * sig2d[i+1] <= 0:
- strength = labs(sig2d[i] - sig2d[i+1])
- if (strength > zcstrength) or \
- (strength == zcstrength and (abs(center - i) <
- abs(center - zcp))):
- zcstrength = strength
- zcp = i
- zca = ax
- free(sig2d)
- #print "zcp: %s (%s)" % (zcp, self.dimensions[ax])
- return 1, ax, zcp
-
- @cython.boundscheck(False)
- @cython.wraparound(False)
- def get_properties(self):
- cdef np.ndarray[np.int64_t, ndim=2] tr = np.empty((3,3), dtype='int64')
- cdef int i
- for i in range(3):
- tr[0,i] = self.left_edge[i]
- tr[1,i] = self.right_edge[i]
- tr[2,i] = self.dimensions[i]
- return tr
-
- at cython.cdivision
-cdef np.int64_t graycode(np.int64_t x):
- return x^(x>>1)
-
- at cython.cdivision
-cdef np.int64_t igraycode(np.int64_t x):
- cdef np.int64_t i, j
- if x == 0:
- return x
- m = <np.int64_t> ceil(log2(x)) + 1
- i, j = x, 1
- while j < m:
- i = i ^ (x>>j)
- j += 1
- return i
-
- at cython.cdivision
-cdef np.int64_t direction(np.int64_t x, np.int64_t n):
- #assert x < 2**n
- if x == 0:
- return 0
- elif x%2 == 0:
- return tsb(x-1, n)%n
- else:
- return tsb(x, n)%n
-
- at cython.cdivision
-cdef np.int64_t tsb(np.int64_t x, np.int64_t width):
- #assert x < 2**width
- cdef np.int64_t i = 0
- while x&1 and i <= width:
- x = x >> 1
- i += 1
- return i
-
- at cython.cdivision
-cdef np.int64_t bitrange(np.int64_t x, np.int64_t width,
- np.int64_t start, np.int64_t end):
- return x >> (width-end) & ((2**(end-start))-1)
-
- at cython.cdivision
-cdef np.int64_t rrot(np.int64_t x, np.int64_t i, np.int64_t width):
- i = i%width
- x = (x>>i) | (x<<width-i)
- return x&(2**width-1)
-
- at cython.cdivision
-cdef np.int64_t lrot(np.int64_t x, np.int64_t i, np.int64_t width):
- i = i%width
- x = (x<<i) | (x>>width-i)
- return x&(2**width-1)
-
- at cython.cdivision
-cdef np.int64_t transform(np.int64_t entry, np.int64_t direction,
- np.int64_t width, np.int64_t x):
- return rrot((x^entry), direction + 1, width)
-
- at cython.cdivision
-cdef np.int64_t entry(np.int64_t x):
- if x == 0: return 0
- return graycode(2*((x-1)/2))
-
- at cython.cdivision
-def get_hilbert_indices(int order, np.ndarray[np.int64_t, ndim=2] left_index):
- # This is inspired by the scurve package by user cortesi on GH.
- cdef int o, i
- cdef np.int64_t x, w, h, e, d, l, b
- cdef np.int64_t p[3]
- cdef np.ndarray[np.int64_t, ndim=1] hilbert_indices
- hilbert_indices = np.zeros(left_index.shape[0], 'int64')
- for o in range(left_index.shape[0]):
- p[0] = left_index[o, 0]
- p[1] = left_index[o, 1]
- p[2] = left_index[o, 2]
- h = e = d = 0
- for i in range(order):
- l = 0
- for x in range(3):
- b = bitrange(p[3-x-1], order, i, i+1)
- l |= (b<<x)
- l = transform(e, d, 3, l)
- w = igraycode(l)
- e = e ^ lrot(entry(w), d+1, 3)
- d = (d + direction(w, 3) + 1)%3
- h = (h<<3)|w
- hilbert_indices[o] = h
- return hilbert_indices
-
+"""
+Wrapping code for Oliver Hahn's RamsesRead++
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: UCSD
+Author: Oliver Hahn <ohahn at stanford.edu>
+Affiliation: KIPAC / Stanford
+Homepage: http://yt.enzotools.org/
+License:
+ Copyright (C) 2010 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/>.
+"""
+
+# Cython wrapping code for Oliver Hahn's RAMSES reader
+from cython.operator cimport dereference as deref, preincrement as inc
+from libc.stdlib cimport malloc, free, abs, calloc, labs
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+cdef extern from "math.h":
+ double ceil(double x)
+ double log2(double x)
+
+cdef inline np.int64_t i64max(np.int64_t i0, np.int64_t i1):
+ if i0 > i1: return i0
+ return i1
+
+cdef inline np.int64_t i64min(np.int64_t i0, np.int64_t i1):
+ if i0 < i1: return i0
+ return i1
+
+cdef extern from "<vector>" namespace "std":
+ cdef cppclass vector[T]:
+ cppclass iterator:
+ T operator*()
+ iterator operator++()
+ bint operator==(iterator)
+ bint operator!=(iterator)
+ vector()
+ void push_back(T&)
+ T& operator[](int)
+ T& at(int)
+ iterator begin()
+ iterator end()
+
+cdef extern from "<map>" namespace "std":
+ cdef cppclass map[A,B]:
+ pass
+
+cdef extern from "string" namespace "std":
+ cdef cppclass string:
+ string(char *cstr)
+ char *c_str()
+ string operator*()
+
+cdef extern from "RAMSES_typedefs.h":
+ pass
+
+cdef extern from "RAMSES_info.hh" namespace "RAMSES":
+ enum codeversion:
+ version1
+ version2
+ version3
+
+ cdef cppclass snapshot:
+ string m_filename
+ codeversion m_version
+ cppclass info_data:
+ unsigned ncpu
+ unsigned ndim
+ unsigned levelmin
+ unsigned levelmax
+ unsigned ngridmax
+ unsigned nstep_coarse
+ double boxlen
+ double time
+ double aexp
+ double H0
+ double omega_m
+ double omega_l
+ double omega_k
+ double omega_b
+ double unit_l
+ double unit_d
+ double unit_t
+
+ info_data m_header
+ vector[double] ind_min
+ vector[double] ind_max
+
+ snapshot(string info_filename, codeversion ver)
+
+ unsigned get_snapshot_num()
+ unsigned getdomain_bykey( double key )
+
+ #void hilbert3d(vector[double] &x, vector[double] &y, vector[double] &z,
+ # vector[double] &order, unsigned bit_length)
+
+cdef extern from "RAMSES_amr_data.hh" namespace "RAMSES::AMR":
+ cdef cppclass vec[real_t]:
+ real_t x, y, z
+ vec( real_t x_, real_t y_, real_t z_)
+ vec ( vec& v)
+ vec ( )
+
+ # This class definition is out of date. I have unrolled the template
+ # below.
+ cdef cppclass cell_locally_essential[id_t, real_t]:
+ id_t m_neighbor[6]
+ id_t m_father
+ id_t m_son[8]
+ real_t m_xg[3]
+ id_t m_cpu
+
+ char m_pos
+
+ cell_locally_essential()
+
+ bint is_refined(int ison)
+
+ cdef cppclass RAMSES_cell:
+ unsigned m_neighbour[6]
+ unsigned m_father
+ unsigned m_son[8]
+ float m_xg[3]
+ unsigned m_cpu
+
+ char m_pos
+
+ RAMSES_cell()
+
+ bint is_refined(int ison)
+
+ cdef cppclass cell_simple[id_t, real_t]:
+ id_t m_son[1]
+ real_t m_xg[3]
+ id_t m_cpu
+
+ char m_pos
+ cell_simple()
+ bint is_refined(int ison)
+
+ # AMR level implementation
+
+ # This class definition is out of date. I have unrolled the template
+ # below.
+ cdef cppclass level[Cell_]:
+ unsigned m_ilevel
+ vector[Cell_] m_level_cells
+ double m_xc[8]
+ double m_yc[8]
+ double m_zc[8]
+
+ # I skipped the typedefs here
+
+ double m_dx
+ unsigned m_nx
+ level (unsigned ilevel)
+
+ void register_cell( Cell_ cell )
+ vector[Cell_].iterator begin()
+ vector[Cell_].iterator end()
+
+ Cell_& operator[]( unsigned i)
+ unsigned size()
+
+ cdef cppclass RAMSES_level:
+ unsigned m_ilevel
+ vector[RAMSES_cell] m_level_cells
+ double m_xc[8]
+ double m_yc[8]
+ double m_zc[8]
+
+ # I skipped the typedefs here
+
+ double m_dx
+ unsigned m_nx
+ RAMSES_level (unsigned ilevel)
+
+ void register_cell( RAMSES_cell cell )
+ vector[RAMSES_cell].iterator begin()
+ vector[RAMSES_cell].iterator end()
+
+ RAMSES_cell& operator[]( unsigned i)
+ unsigned size()
+
+ # This class definition is out of date. I have unrolled the template
+ # below.
+ cdef cppclass tree[Cell_, Level_]:
+ cppclass header:
+ vector[int] nx
+ vector[int] nout
+ vector[int] nsteps
+ int ncpu
+ int ndim
+ int nlevelmax
+ int ngridmax
+ int nboundary
+ int ngrid_current
+ double boxlen
+ vector[double] tout
+ vector[double] aout
+ vector[double] dtold
+ vector[double] dtnew
+ vector[double] cosm
+ vector[double] timing
+ double t
+ double mass_sph
+
+ vector[Level_] m_AMR_levels
+ vector[unsigned] mheadl, m_numbl, m_taill
+
+ int m_cpu
+ int m_minlevel
+ int m_maxlevel
+ string m_fname
+ unsigned m_ncoarse
+ header m_header
+
+ # This is from later on in the .hh file ... I don't think we need them
+ # typedef proto_iterator<const tree*> const_iterator;
+ # typedef proto_iterator<tree *> iterator;
+
+ tree (snapshot &snap, int cpu, int maxlevel, int minlevel = 1)
+
+ cppclass tree_iterator "RAMSES::AMR::RAMSES_tree::iterator":
+ tree_iterator operator*()
+ RAMSES_cell operator*()
+ tree_iterator begin()
+ tree_iterator end()
+ tree_iterator to_parent()
+ tree_iterator get_parent()
+ void next()
+ bint operator!=(tree_iterator other)
+ unsigned get_cell_father()
+ bint is_finest(int ison)
+ int get_absolute_position()
+ int get_domain()
+
+ cdef cppclass RAMSES_tree:
+ # This is, strictly speaking, not a header. But, I believe it is
+ # going to work alright.
+ cppclass header:
+ vector[int] nx
+ vector[int] nout
+ vector[int] nsteps
+ int ncpu
+ int ndim
+ int nlevelmax
+ int ngridmax
+ int nboundary
+ int ngrid_current
+ double boxlen
+ vector[double] tout
+ vector[double] aout
+ vector[double] dtold
+ vector[double] dtnew
+ vector[double] cosm
+ vector[double] timing
+ double t
+ double mass_sph
+
+ vector[RAMSES_level] m_AMR_levels
+ vector[unsigned] mheadl, m_numbl, m_taill
+
+ int m_cpu
+ int m_minlevel
+ int m_maxlevel
+ string m_fname
+ unsigned m_ncoarse
+ header m_header
+
+ unsigned size()
+
+ # This is from later on in the .hh file ... I don't think we need them
+ # typedef proto_iterator<const tree*> const_iterator;
+ # typedef proto_iterator<tree *> iterator;
+
+ RAMSES_tree(snapshot &snap, int cpu, int maxlevel, int minlevel)
+ void read()
+
+ tree_iterator begin(int ilevel)
+ tree_iterator end(int ilevel)
+
+ tree_iterator begin()
+ tree_iterator end()
+
+ vec[double] cell_pos_double "cell_pos<double>" (tree_iterator it, unsigned ind)
+ vec[double] grid_pos_double "grid_pos<double>" (tree_iterator it)
+ vec[float] cell_pos_float "cell_pos<float>" (tree_iterator it, unsigned ind)
+ vec[float] grid_pos_float "grid_pos<float>" (tree_iterator it)
+
+cdef extern from "RAMSES_amr_data.hh" namespace "RAMSES::HYDRO":
+ enum hydro_var:
+ density
+ velocity_x
+ velocity_y
+ velocity_z
+ pressure
+ metallicit
+
+ char ramses_hydro_variables[][64]
+
+ # There are a number of classes here that we could wrap and utilize.
+ # However, I will only wrap the methods I know we need.
+
+ # I have no idea if this will work.
+ cdef cppclass TreeTypeIterator[TreeType_]:
+ pass
+
+ # This class definition is out of date. I have unrolled the template
+ # below.
+ cdef cppclass data[TreeType_, Real_]:
+ cppclass header:
+ unsigned ncpu
+ unsigned nvar
+ unsigned ndim
+ unsigned nlevelmax
+ unsigned nboundary
+ double gamma
+ string m_fname
+ header m_header
+
+ # I don't want to implement proto_data, so we put this here
+ Real_& cell_value(TreeTypeIterator[TreeType_] &it, int ind)
+
+ unsigned m_nvars
+ vector[string] m_varnames
+ map[string, unsigned] m_var_name_map
+
+ data(TreeType_ &AMRtree)
+ #_OutputIterator get_var_names[_OutputIterator](_OutputIterator names)
+ void read(string varname)
+
+ cdef cppclass RAMSES_hydro_data:
+ cppclass header:
+ unsigned ncpu
+ unsigned nvar
+ unsigned ndim
+ unsigned nlevelmax
+ unsigned nboundary
+ double gamma
+ string m_fname
+ header m_header
+
+ # I don't want to implement proto_data, so we put this here
+ double cell_value (tree_iterator &it, int ind)
+
+ unsigned m_nvars
+ vector[string] m_varnames
+ map[string, unsigned] m_var_name_map
+
+ RAMSES_hydro_data(RAMSES_tree &AMRtree)
+ #_OutputIterator get_var_names[_OutputIterator](_OutputIterator names)
+ void read(string varname)
+
+ vector[vector[double]] m_var_array
+
+cdef class RAMSES_tree_proxy:
+ cdef string *snapshot_name
+ cdef snapshot *rsnap
+
+ cdef RAMSES_tree** trees
+ cdef RAMSES_hydro_data*** hydro_datas
+
+ cdef int **loaded
+
+ cdef public object field_ind
+ cdef public object field_names
+
+ # We will store this here so that we have a record, independent of the
+ # header, of how many things we have allocated
+ cdef int ndomains, nfields
+
+ def __cinit__(self, char *fn):
+ cdef int idomain, ifield, ii
+ cdef RAMSES_tree *local_tree
+ cdef RAMSES_hydro_data *local_hydro_data
+ self.snapshot_name = new string(fn)
+ self.rsnap = new snapshot(deref(self.snapshot_name), version3)
+ # We now have to get our field names to fill our array
+ self.trees = <RAMSES_tree**>\
+ malloc(sizeof(RAMSES_tree*) * self.rsnap.m_header.ncpu)
+ self.hydro_datas = <RAMSES_hydro_data ***>\
+ malloc(sizeof(RAMSES_hydro_data**) * self.rsnap.m_header.ncpu)
+ self.ndomains = self.rsnap.m_header.ncpu
+ #for ii in range(self.ndomains): self.trees[ii] = NULL
+ # Note we don't do ncpu + 1
+ for idomain in range(self.rsnap.m_header.ncpu):
+ # we don't delete local_tree
+ local_tree = new RAMSES_tree(deref(self.rsnap), idomain + 1,
+ self.rsnap.m_header.levelmax, 0)
+ local_tree.read()
+ local_hydro_data = new RAMSES_hydro_data(deref(local_tree))
+ self.hydro_datas[idomain] = <RAMSES_hydro_data **>\
+ malloc(sizeof(RAMSES_hydro_data*) * local_hydro_data.m_nvars)
+ for ii in range(local_hydro_data.m_nvars):
+ self.hydro_datas[idomain][ii] = \
+ new RAMSES_hydro_data(deref(local_tree))
+ self.trees[idomain] = local_tree
+ # We do not delete the final snapshot, which we'll use later
+ if idomain + 1 < self.rsnap.m_header.ncpu:
+ del local_hydro_data
+ # Only once, we read all the field names
+ self.nfields = local_hydro_data.m_nvars
+ cdef string *field_name
+ self.field_names = []
+ self.field_ind = {}
+ self.loaded = <int **> malloc(sizeof(int) * local_hydro_data.m_nvars)
+ for idomain in range(self.ndomains):
+ self.loaded[idomain] = <int *> malloc(
+ sizeof(int) * local_hydro_data.m_nvars)
+ for ifield in range(local_hydro_data.m_nvars):
+ self.loaded[idomain][ifield] = 0
+ for ifield in range(local_hydro_data.m_nvars):
+ field_name = &(local_hydro_data.m_varnames[ifield])
+ # Does this leak?
+ self.field_names.append(field_name.c_str())
+ self.field_ind[self.field_names[-1]] = ifield
+ # This all needs to be cleaned up in the deallocator
+ del local_hydro_data
+
+ def __dealloc__(self):
+ import traceback; traceback.print_stack()
+ cdef int idomain, ifield
+ # To ensure that 'delete' is used, not 'free',
+ # we allocate temporary variables.
+ cdef RAMSES_tree *temp_tree
+ cdef RAMSES_hydro_data *temp_hdata
+ for idomain in range(self.ndomains):
+ for ifield in range(self.nfields):
+ temp_hdata = self.hydro_datas[idomain][ifield]
+ del temp_hdata
+ temp_tree = self.trees[idomain]
+ del temp_tree
+ free(self.loaded[idomain])
+ free(self.hydro_datas[idomain])
+ free(self.trees)
+ free(self.hydro_datas)
+ free(self.loaded)
+ if self.snapshot_name != NULL: del self.snapshot_name
+ if self.rsnap != NULL: del self.rsnap
+
+ def count_zones(self):
+ # We need to do simulation domains here
+
+ cdef unsigned idomain, ilevel
+ cdef RAMSES_tree *local_tree
+ cdef RAMSES_hydro_data *local_hydro_data
+ cdef RAMSES_level *local_level
+ cdef tree_iterator grid_it, grid_end
+
+ # All the loop-local pointers must be declared up here
+
+ cell_count = []
+ cdef int local_count = 0
+ for ilevel in range(self.rsnap.m_header.levelmax + 1):
+ cell_count.append(0)
+ for idomain in range(1, self.rsnap.m_header.ncpu + 1):
+ local_tree = new RAMSES_tree(deref(self.rsnap), idomain,
+ self.rsnap.m_header.levelmax, 0)
+ local_tree.read()
+ local_hydro_data = new RAMSES_hydro_data(deref(local_tree))
+ for ilevel in range(local_tree.m_maxlevel + 1):
+ local_count = 0
+ local_level = &local_tree.m_AMR_levels[ilevel]
+ grid_it = local_tree.begin(ilevel)
+ grid_end = local_tree.end(ilevel)
+ while grid_it != grid_end:
+ local_count += (grid_it.get_domain() == idomain)
+ grid_it.next()
+ cell_count[ilevel] += local_count
+ del local_tree, local_hydro_data
+
+ return cell_count
+
+ def ensure_loaded(self, char *varname, int domain_index):
+ # this domain_index must be zero-indexed
+ cdef int varindex = self.field_ind[varname]
+ cdef string *field_name = new string(varname)
+ if self.loaded[domain_index][varindex] == 1:
+ return
+ print "READING FROM DISK", varname, domain_index, varindex
+ self.hydro_datas[domain_index][varindex].read(deref(field_name))
+ self.loaded[domain_index][varindex] = 1
+ del field_name
+
+ def clear_tree(self, char *varname, int domain_index):
+ # this domain_index must be zero-indexed
+ # We delete and re-create
+ cdef int varindex = self.field_ind[varname]
+ cdef string *field_name = new string(varname)
+ if self.loaded[domain_index][varindex] == 0: return
+ cdef RAMSES_hydro_data *temp_hdata = self.hydro_datas[domain_index][varindex]
+ del temp_hdata
+ self.hydro_datas[domain_index - 1][varindex] = \
+ new RAMSES_hydro_data(deref(self.trees[domain_index]))
+ self.loaded[domain_index][varindex] = 0
+ del field_name
+
+ def get_file_info(self):
+ header_info = {}
+ header_info["ncpu"] = self.rsnap.m_header.ncpu
+ header_info["ndim"] = self.rsnap.m_header.ndim
+ header_info["levelmin"] = self.rsnap.m_header.levelmin
+ header_info["levelmax"] = self.rsnap.m_header.levelmax
+ header_info["ngridmax"] = self.rsnap.m_header.ngridmax
+ header_info["nstep_coarse"] = self.rsnap.m_header.nstep_coarse
+ header_info["boxlen"] = self.rsnap.m_header.boxlen
+ header_info["time"] = self.rsnap.m_header.time
+ header_info["aexp"] = self.rsnap.m_header.aexp
+ header_info["H0"] = self.rsnap.m_header.H0
+ header_info["omega_m"] = self.rsnap.m_header.omega_m
+ header_info["omega_l"] = self.rsnap.m_header.omega_l
+ header_info["omega_k"] = self.rsnap.m_header.omega_k
+ header_info["omega_b"] = self.rsnap.m_header.omega_b
+ header_info["unit_l"] = self.rsnap.m_header.unit_l
+ header_info["unit_d"] = self.rsnap.m_header.unit_d
+ header_info["unit_t"] = self.rsnap.m_header.unit_t
+
+ # Now we grab some from the trees
+ cdef np.ndarray[np.int32_t, ndim=1] top_grid_dims = np.zeros(3, "int32")
+ cdef int i
+ for i in range(3):
+ # Note that nx is really boundary conditions. We always have 2.
+ top_grid_dims[i] = self.trees[0].m_header.nx[i]
+ top_grid_dims[i] = 2
+ header_info["nx"] = top_grid_dims
+
+ return header_info
+
+ def fill_hierarchy_arrays(self,
+ np.ndarray[np.int32_t, ndim=1] top_grid_dims,
+ np.ndarray[np.float64_t, ndim=2] left_edges,
+ np.ndarray[np.float64_t, ndim=2] right_edges,
+ np.ndarray[np.int32_t, ndim=2] grid_levels,
+ np.ndarray[np.int64_t, ndim=2] grid_file_locations,
+ np.ndarray[np.uint64_t, ndim=1] grid_hilbert_indices,
+ np.ndarray[np.int32_t, ndim=2] child_mask,
+ np.ndarray[np.float64_t, ndim=1] domain_left,
+ np.ndarray[np.float64_t, ndim=1] domain_right):
+ # We need to do simulation domains here
+
+ cdef unsigned idomain, ilevel
+
+ # All the loop-local pointers must be declared up here
+ cdef RAMSES_tree *local_tree
+ cdef RAMSES_hydro_data *local_hydro_data
+ cdef unsigned father
+
+ cdef tree_iterator grid_it, grid_end, father_it
+ cdef vec[double] gvec
+ cdef int grid_ind = 0, grid_aind = 0
+ cdef unsigned parent_ind
+ cdef bint ci
+ cdef double pos[3]
+ cdef double grid_half_width
+ cdef unsigned long rv
+
+ cdef np.int32_t rr
+ cdef int i
+ cell_count = []
+ level_cell_counts = {}
+ for idomain in range(1, self.rsnap.m_header.ncpu + 1):
+ local_tree = new RAMSES_tree(deref(self.rsnap), idomain,
+ self.rsnap.m_header.levelmax, 0)
+ local_tree.read()
+ local_hydro_data = new RAMSES_hydro_data(deref(local_tree))
+ for ilevel in range(local_tree.m_maxlevel + 1):
+ # this gets overwritten for every domain, which is okay
+ level_cell_counts[ilevel] = grid_ind
+ #grid_half_width = self.rsnap.m_header.boxlen / \
+ grid_half_width = 1.0 / \
+ (2**(ilevel) * top_grid_dims[0])
+ grid_it = local_tree.begin(ilevel)
+ grid_end = local_tree.end(ilevel)
+ while grid_it != grid_end:
+ if grid_it.get_domain() != idomain:
+ grid_ind += 1
+ grid_it.next()
+ continue
+ gvec = local_tree.grid_pos_double(grid_it)
+ left_edges[grid_aind, 0] = pos[0] = gvec.x - grid_half_width
+ left_edges[grid_aind, 1] = pos[1] = gvec.y - grid_half_width
+ left_edges[grid_aind, 2] = pos[2] = gvec.z - grid_half_width
+ for i in range(3):
+ pos[i] = (pos[i] - domain_left[i]) / (domain_right[i] - domain_left[i])
+ right_edges[grid_aind, 0] = gvec.x + grid_half_width
+ right_edges[grid_aind, 1] = gvec.y + grid_half_width
+ right_edges[grid_aind, 2] = gvec.z + grid_half_width
+ grid_levels[grid_aind, 0] = ilevel
+ # Now the harder part
+ father_it = grid_it.get_parent()
+ grid_file_locations[grid_aind, 0] = <np.int64_t> idomain
+ grid_file_locations[grid_aind, 1] = grid_ind - level_cell_counts[ilevel]
+ parent_ind = father_it.get_absolute_position()
+ if ilevel > 0:
+ # We calculate the REAL parent index
+ grid_file_locations[grid_aind, 2] = \
+ level_cell_counts[ilevel - 1] + parent_ind
+ else:
+ grid_file_locations[grid_aind, 2] = -1
+ for ci in range(8):
+ rr = <np.int32_t> grid_it.is_finest(ci)
+ child_mask[grid_aind, ci] = rr
+ grid_ind += 1
+ grid_aind += 1
+ grid_it.next()
+ del local_tree, local_hydro_data
+
+ def read_oct_grid(self, char *field, int level, int domain, int grid_id):
+
+ self.ensure_loaded(field, domain - 1)
+ cdef int varindex = self.field_ind[field]
+ cdef int i
+
+ cdef np.ndarray[np.float64_t, ndim=3] tr = np.empty((2,2,2), dtype='float64',
+ order='F')
+ cdef tree_iterator grid_it, grid_end
+ cdef double* data = <double*> tr.data
+
+ cdef RAMSES_tree *local_tree = self.trees[domain - 1]
+ cdef RAMSES_hydro_data *local_hydro_data = self.hydro_datas[domain - 1][varindex]
+
+ #inline ValueType_& cell_value( const typename TreeType_::iterator& it,
+ # int ind )
+ #{
+ # unsigned ipos = it.get_absolute_position();
+ # unsigned ilevel = it.get_level();//-m_minlevel;
+ # return (m_var_array[ilevel])[m_twotondim*ipos+ind];
+ #}
+
+ for i in range(8):
+ data[i] = local_hydro_data.m_var_array[level][8*grid_id+i]
+ return tr
+
+ def read_grid(self, char *field,
+ np.ndarray[np.int64_t, ndim=1] start_index,
+ np.ndarray[np.int32_t, ndim=1] grid_dims,
+ np.ndarray[np.float64_t, ndim=3] data,
+ np.ndarray[np.int32_t, ndim=3] filled,
+ int level, int ref_factor,
+ component_grid_info):
+ cdef int varindex = self.field_ind[field]
+ cdef RAMSES_tree *local_tree = NULL
+ cdef RAMSES_hydro_data *local_hydro_data = NULL
+
+ cdef int gi, i, j, k, domain, offset
+ cdef int ir, jr, kr
+ cdef int offi, offj, offk, odind
+ cdef np.int64_t di, dj, dk
+ cdef np.ndarray[np.int64_t, ndim=1] ogrid_info
+ cdef np.ndarray[np.int64_t, ndim=1] og_start_index
+ cdef np.float64_t temp_data
+ cdef np.int64_t end_index[3]
+ cdef int to_fill = 0
+ # Note that indexing into a cell is:
+ # (k*2 + j)*2 + i
+ for i in range(3):
+ end_index[i] = start_index[i] + grid_dims[i]
+ for gi in range(len(component_grid_info)):
+ ogrid_info = component_grid_info[gi]
+ domain = ogrid_info[0]
+ #print "Loading", domain, ogrid_info
+ self.ensure_loaded(field, domain - 1)
+ local_tree = self.trees[domain - 1]
+ local_hydro_data = self.hydro_datas[domain - 1][varindex]
+ offset = ogrid_info[1]
+ og_start_index = ogrid_info[3:]
+ for i in range(2*ref_factor):
+ di = i + og_start_index[0] * ref_factor
+ if di < start_index[0] or di >= end_index[0]: continue
+ ir = <int> (i / ref_factor)
+ for j in range(2 * ref_factor):
+ dj = j + og_start_index[1] * ref_factor
+ if dj < start_index[1] or dj >= end_index[1]: continue
+ jr = <int> (j / ref_factor)
+ for k in range(2 * ref_factor):
+ dk = k + og_start_index[2] * ref_factor
+ if dk < start_index[2] or dk >= end_index[2]: continue
+ kr = <int> (k / ref_factor)
+ offi = di - start_index[0]
+ offj = dj - start_index[1]
+ offk = dk - start_index[2]
+ #print offi, filled.shape[0],
+ #print offj, filled.shape[1],
+ #print offk, filled.shape[2]
+ if filled[offi, offj, offk] == 1: continue
+
+ odind = (kr*2 + jr)*2 + ir
+ temp_data = local_hydro_data.m_var_array[
+ level][8*offset + odind]
+ data[offi, offj, offk] = temp_data
+ filled[offi, offj, offk] = 1
+ to_fill += 1
+ return to_fill
+
+cdef class ProtoSubgrid:
+ cdef np.int64_t *signature[3]
+ cdef np.int64_t left_edge[3]
+ cdef np.int64_t right_edge[3]
+ cdef np.int64_t dimensions[3]
+ cdef public np.float64_t efficiency
+ cdef public object sigs
+ cdef public object grid_file_locations
+ cdef public object dd
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def __cinit__(self,
+ np.ndarray[np.int64_t, ndim=1] left_index,
+ np.ndarray[np.int64_t, ndim=1] dimensions,
+ np.ndarray[np.int64_t, ndim=2] left_edges,
+ np.ndarray[np.int64_t, ndim=2] right_edges,
+ np.ndarray[np.int64_t, ndim=2] grid_dimensions,
+ np.ndarray[np.int64_t, ndim=2] grid_file_locations):
+ # This also includes the shrinking step.
+ cdef int i, ci, ng = left_edges.shape[0]
+ cdef np.ndarray temp_arr
+ cdef int l0, r0, l1, r1, l2, r2, i0, i1, i2
+ cdef np.int64_t temp_l[3], temp_r[3], ncells
+ cdef np.float64_t efficiency
+ self.sigs = []
+ for i in range(3):
+ temp_l[i] = left_index[i] + dimensions[i]
+ temp_r[i] = left_index[i]
+ self.signature[i] = NULL
+ for gi in range(ng):
+ if left_edges[gi,0] > left_index[0]+dimensions[0] or \
+ right_edges[gi,0] < left_index[0] or \
+ left_edges[gi,1] > left_index[1]+dimensions[1] or \
+ right_edges[gi,1] < left_index[1] or \
+ left_edges[gi,2] > left_index[2]+dimensions[2] or \
+ right_edges[gi,2] < left_index[2]:
+ #print "Skipping grid", gi, "which lies outside out box"
+ continue
+ for i in range(3):
+ temp_l[i] = i64min(left_edges[gi,i], temp_l[i])
+ temp_r[i] = i64max(right_edges[gi,i], temp_r[i])
+ for i in range(3):
+ self.left_edge[i] = i64max(temp_l[i], left_index[i])
+ self.right_edge[i] = i64min(temp_r[i], left_index[i] + dimensions[i])
+ self.dimensions[i] = self.right_edge[i] - self.left_edge[i]
+ if self.dimensions[i] <= 0:
+ self.efficiency = -1.0
+ return
+ self.sigs.append(np.zeros(self.dimensions[i], 'int64'))
+ #print self.sigs[0].size, self.sigs[1].size, self.sigs[2].size
+
+ # My guess is that this whole loop could be done more efficiently.
+ # However, this is clear and straightforward, so it is a good first
+ # pass.
+ cdef np.ndarray[np.int64_t, ndim=1] sig0, sig1, sig2
+ sig0 = self.sigs[0]
+ sig1 = self.sigs[1]
+ sig2 = self.sigs[2]
+ efficiency = 0.0
+ cdef int used
+ self.grid_file_locations = []
+ for gi in range(ng):
+ used = 0
+ nnn = 0
+ for l0 in range(grid_dimensions[gi, 0]):
+ i0 = left_edges[gi, 0] + l0
+ if i0 < self.left_edge[0]: continue
+ if i0 >= self.right_edge[0]: break
+ for l1 in range(grid_dimensions[gi, 1]):
+ i1 = left_edges[gi, 1] + l1
+ if i1 < self.left_edge[1]: continue
+ if i1 >= self.right_edge[1]: break
+ for l2 in range(grid_dimensions[gi, 2]):
+ i2 = left_edges[gi, 2] + l2
+ if i2 < self.left_edge[2]: continue
+ if i2 >= self.right_edge[2]: break
+ i = i0 - self.left_edge[0]
+ sig0[i] += 1
+ i = i1 - self.left_edge[1]
+ sig1[i] += 1
+ i = i2 - self.left_edge[2]
+ sig2[i] += 1
+ efficiency += 1
+ used = 1
+ if used == 1:
+ grid_file_locations[gi,3] = left_edges[gi, 0]
+ grid_file_locations[gi,4] = left_edges[gi, 1]
+ grid_file_locations[gi,5] = left_edges[gi, 2]
+ self.grid_file_locations.append(grid_file_locations[gi,:])
+
+ self.dd = np.ones(3, dtype='int64')
+ for i in range(3):
+ efficiency /= self.dimensions[i]
+ self.dd[i] = self.dimensions[i]
+ #print "Efficiency is %0.3e" % (efficiency)
+ self.efficiency = efficiency
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def find_split(self):
+ # First look for zeros
+ cdef int i, center, ax
+ cdef np.ndarray[ndim=1, dtype=np.int64_t] axes
+ cdef np.int64_t strength, zcstrength, zcp
+ axes = np.argsort(self.dd)[::-1]
+ cdef np.ndarray[np.int64_t] sig
+ for axi in range(3):
+ ax = axes[axi]
+ center = self.dimensions[ax] / 2
+ sig = self.sigs[ax]
+ for i in range(self.dimensions[ax]):
+ if sig[i] == 0 and i > 0 and i < self.dimensions[ax] - 1:
+ #print "zero: %s (%s)" % (i, self.dimensions[ax])
+ return 0, ax, i
+ zcstrength = 0
+ zcp = 0
+ zca = -1
+ cdef int temp
+ cdef np.int64_t *sig2d
+ for axi in range(3):
+ ax = axes[axi]
+ sig = self.sigs[ax]
+ sig2d = <np.int64_t *> malloc(sizeof(np.int64_t) * self.dimensions[ax])
+ sig2d[0] = sig2d[self.dimensions[ax]-1] = 0
+ for i in range(1, self.dimensions[ax] - 1):
+ sig2d[i] = sig[i-1] - 2*sig[i] + sig[i+1]
+ for i in range(1, self.dimensions[ax] - 1):
+ if sig2d[i] * sig2d[i+1] <= 0:
+ strength = labs(sig2d[i] - sig2d[i+1])
+ if (strength > zcstrength) or \
+ (strength == zcstrength and (abs(center - i) <
+ abs(center - zcp))):
+ zcstrength = strength
+ zcp = i
+ zca = ax
+ free(sig2d)
+ #print "zcp: %s (%s)" % (zcp, self.dimensions[ax])
+ return 1, ax, zcp
+
+ @cython.boundscheck(False)
+ @cython.wraparound(False)
+ def get_properties(self):
+ cdef np.ndarray[np.int64_t, ndim=2] tr = np.empty((3,3), dtype='int64')
+ cdef int i
+ for i in range(3):
+ tr[0,i] = self.left_edge[i]
+ tr[1,i] = self.right_edge[i]
+ tr[2,i] = self.dimensions[i]
+ return tr
+
+ at cython.cdivision(True)
+cdef np.int64_t graycode(np.int64_t x):
+ return x^(x>>1)
+
+ at cython.cdivision(True)
+cdef np.int64_t igraycode(np.int64_t x):
+ cdef np.int64_t i, j
+ if x == 0:
+ return x
+ m = <np.int64_t> ceil(log2(x)) + 1
+ i, j = x, 1
+ while j < m:
+ i = i ^ (x>>j)
+ j += 1
+ return i
+
+ at cython.cdivision(True)
+cdef np.int64_t direction(np.int64_t x, np.int64_t n):
+ #assert x < 2**n
+ if x == 0:
+ return 0
+ elif x%2 == 0:
+ return tsb(x-1, n)%n
+ else:
+ return tsb(x, n)%n
+
+ at cython.cdivision(True)
+cdef np.int64_t tsb(np.int64_t x, np.int64_t width):
+ #assert x < 2**width
+ cdef np.int64_t i = 0
+ while x&1 and i <= width:
+ x = x >> 1
+ i += 1
+ return i
+
+ at cython.cdivision(True)
+cdef np.int64_t bitrange(np.int64_t x, np.int64_t width,
+ np.int64_t start, np.int64_t end):
+ return x >> (width-end) & ((2**(end-start))-1)
+
+ at cython.cdivision(True)
+cdef np.int64_t rrot(np.int64_t x, np.int64_t i, np.int64_t width):
+ i = i%width
+ x = (x>>i) | (x<<width-i)
+ return x&(2**width-1)
+
+ at cython.cdivision(True)
+cdef np.int64_t lrot(np.int64_t x, np.int64_t i, np.int64_t width):
+ i = i%width
+ x = (x<<i) | (x>>width-i)
+ return x&(2**width-1)
+
+ at cython.cdivision(True)
+cdef np.int64_t transform(np.int64_t entry, np.int64_t direction,
+ np.int64_t width, np.int64_t x):
+ return rrot((x^entry), direction + 1, width)
+
+ at cython.cdivision(True)
+cdef np.int64_t entry(np.int64_t x):
+ if x == 0: return 0
+ return graycode(2*((x-1)/2))
+
+ at cython.cdivision(True)
+def get_hilbert_indices(int order, np.ndarray[np.int64_t, ndim=2] left_index):
+ # This is inspired by the scurve package by user cortesi on GH.
+ cdef int o, i
+ cdef np.int64_t x, w, h, e, d, l, b
+ cdef np.int64_t p[3]
+ cdef np.ndarray[np.int64_t, ndim=1] hilbert_indices
+ hilbert_indices = np.zeros(left_index.shape[0], 'int64')
+ for o in range(left_index.shape[0]):
+ p[0] = left_index[o, 0]
+ p[1] = left_index[o, 1]
+ p[2] = left_index[o, 2]
+ h = e = d = 0
+ for i in range(order):
+ l = 0
+ for x in range(3):
+ b = bitrange(p[3-x-1], order, i, i+1)
+ l |= (b<<x)
+ l = transform(e, d, 3, l)
+ w = igraycode(l)
+ e = e ^ lrot(entry(w), d+1, 3)
+ d = (d + direction(w, 3) + 1)%3
+ h = (h<<3)|w
+ hilbert_indices[o] = h
+ return hilbert_indices
+
https://bitbucket.org/yt_analysis/yt/changeset/6cbf80a388ce/
changeset: 6cbf80a388ce
branch: yt
user: Andrew Myers
date: 2011-11-12 03:18:54
summary: fixing the a grid offset problem in the Chombo reader. Also adding a redundant 'Density' field, since PlotCollection tries to use this if you don't pass in a center as an argument.
affected #: 2 files
diff -r 2922825feafeb209651d127da1fc8ee94f39f907 -r 6cbf80a388ceae4a41da960362d7b2ceeecfb6bf yt/frontends/chombo/data_structures.py
--- a/yt/frontends/chombo/data_structures.py
+++ b/yt/frontends/chombo/data_structures.py
@@ -68,9 +68,25 @@
self.Parent = []
self.Children = []
self.Level = level
- self.start_index = start.copy()#.transpose()
- self.stop_index = stop.copy()#.transpose()
self.ActiveDimensions = stop - start + 1
+
+ def get_global_startindex(self):
+ """
+ Return the integer starting index for each dimension at the current
+ level.
+
+ """
+ if self.start_index != None:
+ return self.start_index
+ if self.Parent == []:
+ iLE = self.LeftEdge - self.pf.domain_left_edge
+ start_index = iLE / self.dds
+ return na.rint(start_index).astype('int64').ravel()
+ pdx = self.Parent.dds
+ start_index = (self.Parent.get_global_startindex()) + \
+ na.rint((self.LeftEdge - self.Parent.LeftEdge)/pdx)
+ self.start_index = (start_index*self.pf.refine_by).astype('int64').ravel()
+ return self.start_index
def _setup_dx(self):
# So first we figure out what the index is. We don't assume
@@ -91,8 +107,8 @@
grid = ChomboGrid
def __init__(self,pf,data_style='chombo_hdf5'):
- self.domain_left_edge = pf.domain_left_edge # need these to determine absolute grid locations
- self.domain_right_edge = pf.domain_right_edge # need these to determine absolute grid locations
+ self.domain_left_edge = pf.domain_left_edge
+ self.domain_right_edge = pf.domain_right_edge
self.data_style = data_style
self.field_info = ChomboFieldContainer()
self.field_indexes = {}
@@ -143,8 +159,8 @@
start = si, stop = ei)
self.grids.append(pg)
self.grids[-1]._level_id = level_id
- self.grid_left_edge[i] = dx*si.astype(self.float_type) #AJC + self.domain_left_edge
- self.grid_right_edge[i] = dx*(ei.astype(self.float_type)+1) #AJC + self.domain_left_edge
+ self.grid_left_edge[i] = dx*si.astype(self.float_type)
+ self.grid_right_edge[i] = dx*(ei.astype(self.float_type)+1)
self.grid_particle_count[i] = 0
self.grid_dimensions[i] = ei - si + 1
i += 1
@@ -229,84 +245,77 @@
"""
if os.path.isfile('pluto.ini'):
self._parse_pluto_file('pluto.ini')
- elif os.path.isfile('orion2.ini'):
- # AJC - this shouldn't be nesecary anymore, now
- # that orion2 has fixed the pluto domian-offset glitch
- self._parse_pluto_file('orion2.ini')
+# elif os.path.isfile('orion2.ini'):
+# self._parse_pluto_file('orion2.ini')
else:
self.unique_identifier = \
int(os.stat(self.parameter_filename)[ST_CTIME])
self.domain_left_edge = self.__calc_left_edge()
self.domain_right_edge = self.__calc_right_edge()
+ self.domain_dimensions = self.__calc_domain_dimensions()
self.dimensionality = 3
- # AJC - make refine_by a list because different levels
- # can have different refinment ratios.
self.refine_by = []
fileh = h5py.File(self.parameter_filename,'r')
for level in range(0,fileh.attrs['max_level']):
self.refine_by.append(fileh['/level_'+str(level)].attrs['ref_ratio'])
- def _parse_pluto_file(self, ini_filename):
- """
- Reads in an inputs file in the 'pluto.ini' format. Probably not
- especially robust at the moment.
- """
- self.fullplotdir = os.path.abspath(self.parameter_filename)
- self.ini_filename = self._localize( \
- self.ini_filename, ini_filename)
- self.unique_identifier = \
- int(os.stat(self.parameter_filename)[ST_CTIME])
- lines = open(self.ini_filename).readlines()
- # read the file line by line, storing important parameters
- for lineI, line in enumerate(lines):
- try:
- param, sep, vals = map(rstrip,line.partition(' '))
- except ValueError:
- mylog.error("ValueError: '%s'", line)
- if pluto2enzoDict.has_key(param):
- paramName = pluto2enzoDict[param]
- t = map(parameterDict[paramName], vals.split())
- if len(t) == 1:
- self.parameters[paramName] = t[0]
- else:
- if paramName == "RefineBy":
- self.parameters[paramName] = t[0]
- else:
- self.parameters[paramName] = t
+ # def _parse_pluto_file(self, ini_filename):
+ # """
+ # Reads in an inputs file in the 'pluto.ini' format. Probably not
+ # especially robust at the moment.
+ # """
+ # self.fullplotdir = os.path.abspath(self.parameter_filename)
+ # self.ini_filename = self._localize( \
+ # self.ini_filename, ini_filename)
+ # self.unique_identifier = \
+ # int(os.stat(self.parameter_filename)[ST_CTIME])
+ # lines = open(self.ini_filename).readlines()
+ # # read the file line by line, storing important parameters
+ # for lineI, line in enumerate(lines):
+ # try:
+ # param, sep, vals = map(rstrip,line.partition(' '))
+ # except ValueError:
+ # mylog.error("ValueError: '%s'", line)
+ # if pluto2enzoDict.has_key(param):
+ # paramName = pluto2enzoDict[param]
+ # t = map(parameterDict[paramName], vals.split())
+ # if len(t) == 1:
+ # self.parameters[paramName] = t[0]
+ # else:
+ # if paramName == "RefineBy":
+ # self.parameters[paramName] = t[0]
+ # else:
+ # self.parameters[paramName] = t
- # assumes 3D for now
- elif param.startswith("X1-grid"):
- t = vals.split()
- low1 = float(t[1])
- high1 = float(t[4])
- N1 = int(t[2])
- elif param.startswith("X2-grid"):
- t = vals.split()
- low2 = float(t[1])
- high2 = float(t[4])
- N2 = int(t[2])
- elif param.startswith("X3-grid"):
- t = vals.split()
- low3 = float(t[1])
- high3 = float(t[4])
- N3 = int(t[2])
+ # # assumes 3D for now
+ # elif param.startswith("X1-grid"):
+ # t = vals.split()
+ # low1 = float(t[1])
+ # high1 = float(t[4])
+ # N1 = int(t[2])
+ # elif param.startswith("X2-grid"):
+ # t = vals.split()
+ # low2 = float(t[1])
+ # high2 = float(t[4])
+ # N2 = int(t[2])
+ # elif param.startswith("X3-grid"):
+ # t = vals.split()
+ # low3 = float(t[1])
+ # high3 = float(t[4])
+ # N3 = int(t[2])
- self.dimensionality = 3
- self.domain_left_edge = na.array([low1,low2,low3])
- self.domain_right_edge = na.array([high1,high2,high3])
- self.domain_dimensions = na.array([N1,N2,N3])
- self.refine_by = self.parameters["RefineBy"]
+ # self.dimensionality = 3
+ # self.domain_left_edge = na.array([low1,low2,low3])
+ # self.domain_right_edge = na.array([high1,high2,high3])
+ # self.domain_dimensions = na.array([N1,N2,N3])
+ # self.refine_by = self.parameters["RefineBy"]
- # AJC - added this routine
- # Note that vanilla Chombo files (independent of pluto or orion2)
- # do not always have left_edge = [0,0,0]
def __calc_left_edge(self):
fileh = h5py.File(self.parameter_filename,'r')
dx0 = fileh['/level_0'].attrs['dx']
- RE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain']))[0:3])
+ LE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain']))[0:3])
fileh.close()
- return RE
-
+ return LE
def __calc_right_edge(self):
fileh = h5py.File(self.parameter_filename,'r')
@@ -314,7 +323,13 @@
RE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain']))[3:] + 1)
fileh.close()
return RE
-
+
+ def __calc_domain_dimensions(self):
+ fileh = h5py.File(self.parameter_filename,'r')
+ L_index = ((na.array(fileh['/level_0'].attrs['prob_domain']))[0:3])
+ R_index = ((na.array(fileh['/level_0'].attrs['prob_domain']))[3:] + 1)
+ return R_index - L_index
+
@classmethod
def _is_valid(self, *args, **kwargs):
try:
diff -r 2922825feafeb209651d127da1fc8ee94f39f907 -r 6cbf80a388ceae4a41da960362d7b2ceeecfb6bf yt/frontends/chombo/fields.py
--- a/yt/frontends/chombo/fields.py
+++ b/yt/frontends/chombo/fields.py
@@ -116,4 +116,15 @@
return data["Z-momentum"]/data["density"]
add_field("z-velocity",function=_zVelocity, take_log=False,
units=r'\rm{cm}/\rm{s}')
+
+def _Density(field,data):
+ """A duplicate of the density field. This is needed because when you try
+ to instantiate a PlotCollection without passing in a center, the code
+ will try to generate one for you using the "Density" field, which gives an error
+ if it isn't defined.
+
+ """
+ return data["density"]
+add_field("Density",function=_Density, take_log=True,
+ units=r'\rm{g}/\rm{cm^3}')
https://bitbucket.org/yt_analysis/yt/changeset/971ae8ee85ba/
changeset: 971ae8ee85ba
branch: yt
user: Andrew Myers
date: 2011-11-12 21:33:46
summary: previous checkin contained a bug that only showed itself on multilevel datasets. Thanks to Andy Cunningham.
affected #: 1 file
diff -r 6cbf80a388ceae4a41da960362d7b2ceeecfb6bf -r 971ae8ee85ba951eb420d075ee847cf24aaaecdb yt/frontends/chombo/data_structures.py
--- a/yt/frontends/chombo/data_structures.py
+++ b/yt/frontends/chombo/data_structures.py
@@ -82,9 +82,9 @@
iLE = self.LeftEdge - self.pf.domain_left_edge
start_index = iLE / self.dds
return na.rint(start_index).astype('int64').ravel()
- pdx = self.Parent.dds
- start_index = (self.Parent.get_global_startindex()) + \
- na.rint((self.LeftEdge - self.Parent.LeftEdge)/pdx)
+ pdx = self.Parent[0].dds
+ start_index = (self.Parent[0].get_global_startindex()) + \
+ na.rint((self.LeftEdge - self.Parent[0].LeftEdge)/pdx)
self.start_index = (start_index*self.pf.refine_by).astype('int64').ravel()
return self.start_index
https://bitbucket.org/yt_analysis/yt/changeset/474656466409/
changeset: 474656466409
branch: yt
user: acunning
date: 2012-02-01 21:42:21
summary: Resolved some issues that appeared with the most recent version of h5py that I had to fix to compile on kraken.
affected #: 2 files
diff -r 971ae8ee85ba951eb420d075ee847cf24aaaecdb -r 47465646640936c1f42803ffd3447f0478319a20 yt/frontends/chombo/data_structures.py
--- a/yt/frontends/chombo/data_structures.py
+++ b/yt/frontends/chombo/data_structures.py
@@ -146,6 +146,7 @@
# this relies on the first Group in the H5 file being
# 'Chombo_global'
levels = f.listnames()[1:]
+ #levels = [fn for fn in f if fn != "Chombo_global"]
self.grids = []
i = 0
for lev in levels:
@@ -313,28 +314,28 @@
def __calc_left_edge(self):
fileh = h5py.File(self.parameter_filename,'r')
dx0 = fileh['/level_0'].attrs['dx']
- LE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain']))[0:3])
+ LE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain'].tolist()))[0:3])
fileh.close()
return LE
def __calc_right_edge(self):
fileh = h5py.File(self.parameter_filename,'r')
dx0 = fileh['/level_0'].attrs['dx']
- RE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain']))[3:] + 1)
+ RE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain'].tolist()))[3:] + 1)
fileh.close()
return RE
def __calc_domain_dimensions(self):
fileh = h5py.File(self.parameter_filename,'r')
- L_index = ((na.array(fileh['/level_0'].attrs['prob_domain']))[0:3])
- R_index = ((na.array(fileh['/level_0'].attrs['prob_domain']))[3:] + 1)
+ L_index = ((na.array(fileh['/level_0'].attrs['prob_domain'].tolist()))[0:3])
+ R_index = ((na.array(fileh['/level_0'].attrs['prob_domain'].tolist()))[3:] + 1)
return R_index - L_index
@classmethod
def _is_valid(self, *args, **kwargs):
try:
fileh = h5py.File(args[0],'r')
- if (fileh.listnames())[0] == 'Chombo_global':
+ if (fileh.keys())[0] == 'Chombo_global':
return True
except:
pass
diff -r 971ae8ee85ba951eb420d075ee847cf24aaaecdb -r 47465646640936c1f42803ffd3447f0478319a20 yt/frontends/chombo/fields.py
--- a/yt/frontends/chombo/fields.py
+++ b/yt/frontends/chombo/fields.py
@@ -47,17 +47,17 @@
add_field("X-momentum", function=lambda a,b: None, take_log=False,
validators = [ValidateDataField("X-Momentum")],
- units=r"",display_name=r"B_x")
+ units=r"",display_name=r"M_x")
ChomboFieldInfo["X-momentum"]._projected_units=r""
add_field("Y-momentum", function=lambda a,b: None, take_log=False,
validators = [ValidateDataField("Y-Momentum")],
- units=r"",display_name=r"B_y")
+ units=r"",display_name=r"M_y")
ChomboFieldInfo["Y-momentum"]._projected_units=r""
add_field("Z-momentum", function=lambda a,b: None, take_log=False,
validators = [ValidateDataField("Z-Momentum")],
- units=r"",display_name=r"B_z")
+ units=r"",display_name=r"M_z")
ChomboFieldInfo["Z-momentum"]._projected_units=r""
add_field("X-magnfield", function=lambda a,b: None, take_log=False,
@@ -128,3 +128,10 @@
add_field("Density",function=_Density, take_log=True,
units=r'\rm{g}/\rm{cm^3}')
+def _pressure(field,data):
+ return (data["X-magnfield"]**2 +
+ data["Y-magnfield"]**2 +
+ data["Z-magnfield"]**2)/2.
+add_field("pressure", function=_MagneticEnergy, take_log=True,
+ units=r"",display_name=r"\rm{erg}/\rm{cm}^3")
+ChomboFieldInfo["MagneticEnergy"]._projected_units=r""
https://bitbucket.org/yt_analysis/yt/changeset/9d5048db0b22/
changeset: 9d5048db0b22
branch: yt
user: acunning
date: 2012-03-02 23:51:27
summary: Fixes so that gamma gets parsed out of orion2.ini and left/right edge parsing
works both with older and newer versions of h5py.
affected #: 1 file
diff -r 47465646640936c1f42803ffd3447f0478319a20 -r 9d5048db0b22e3a10b92af120e691e3c193dd548 yt/frontends/chombo/data_structures.py
--- a/yt/frontends/chombo/data_structures.py
+++ b/yt/frontends/chombo/data_structures.py
@@ -246,9 +246,8 @@
"""
if os.path.isfile('pluto.ini'):
self._parse_pluto_file('pluto.ini')
-# elif os.path.isfile('orion2.ini'):
-# self._parse_pluto_file('orion2.ini')
else:
+ if os.path.isfile('orion2.ini'): self._parse_pluto_file('orion2.ini')
self.unique_identifier = \
int(os.stat(self.parameter_filename)[ST_CTIME])
self.domain_left_edge = self.__calc_left_edge()
@@ -260,75 +259,52 @@
for level in range(0,fileh.attrs['max_level']):
self.refine_by.append(fileh['/level_'+str(level)].attrs['ref_ratio'])
- # def _parse_pluto_file(self, ini_filename):
- # """
- # Reads in an inputs file in the 'pluto.ini' format. Probably not
- # especially robust at the moment.
- # """
- # self.fullplotdir = os.path.abspath(self.parameter_filename)
- # self.ini_filename = self._localize( \
- # self.ini_filename, ini_filename)
- # self.unique_identifier = \
- # int(os.stat(self.parameter_filename)[ST_CTIME])
- # lines = open(self.ini_filename).readlines()
- # # read the file line by line, storing important parameters
- # for lineI, line in enumerate(lines):
- # try:
- # param, sep, vals = map(rstrip,line.partition(' '))
- # except ValueError:
- # mylog.error("ValueError: '%s'", line)
- # if pluto2enzoDict.has_key(param):
- # paramName = pluto2enzoDict[param]
- # t = map(parameterDict[paramName], vals.split())
- # if len(t) == 1:
- # self.parameters[paramName] = t[0]
- # else:
- # if paramName == "RefineBy":
- # self.parameters[paramName] = t[0]
- # else:
- # self.parameters[paramName] = t
-
- # # assumes 3D for now
- # elif param.startswith("X1-grid"):
- # t = vals.split()
- # low1 = float(t[1])
- # high1 = float(t[4])
- # N1 = int(t[2])
- # elif param.startswith("X2-grid"):
- # t = vals.split()
- # low2 = float(t[1])
- # high2 = float(t[4])
- # N2 = int(t[2])
- # elif param.startswith("X3-grid"):
- # t = vals.split()
- # low3 = float(t[1])
- # high3 = float(t[4])
- # N3 = int(t[2])
-
- # self.dimensionality = 3
- # self.domain_left_edge = na.array([low1,low2,low3])
- # self.domain_right_edge = na.array([high1,high2,high3])
- # self.domain_dimensions = na.array([N1,N2,N3])
- # self.refine_by = self.parameters["RefineBy"]
+ def _parse_pluto_file(self, ini_filename):
+ """
+ Reads in an inputs file in the 'pluto.ini' format. Probably not
+ especially robust at the moment.
+ """
+ self.fullplotdir = os.path.abspath(self.parameter_filename)
+ self.ini_filename = self._localize( \
+ self.ini_filename, ini_filename)
+ self.unique_identifier = \
+ int(os.stat(self.parameter_filename)[ST_CTIME])
+ lines = open(self.ini_filename).readlines()
+ # read the file line by line, storing important parameters
+ for lineI, line in enumerate(lines):
+ try:
+ param, sep, vals = map(rstrip,line.partition(' '))
+ except ValueError:
+ mylog.error("ValueError: '%s'", line)
+ if pluto2enzoDict.has_key(param):
+ paramName = pluto2enzoDict[param]
+ t = map(parameterDict[paramName], vals.split())
+ if len(t) == 1:
+ self.parameters[paramName] = t[0]
+ else:
+ if paramName == "RefineBy":
+ self.parameters[paramName] = t[0]
+ else:
+ self.parameters[paramName] = t
def __calc_left_edge(self):
fileh = h5py.File(self.parameter_filename,'r')
dx0 = fileh['/level_0'].attrs['dx']
- LE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain'].tolist()))[0:3])
+ LE = dx0*((na.array(list(fileh['/level_0'].attrs['prob_domain'])))[0:3])
fileh.close()
return LE
def __calc_right_edge(self):
fileh = h5py.File(self.parameter_filename,'r')
dx0 = fileh['/level_0'].attrs['dx']
- RE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain'].tolist()))[3:] + 1)
+ RE = dx0*((na.array(list(fileh['/level_0'].attrs['prob_domain'])))[3:] + 1)
fileh.close()
return RE
def __calc_domain_dimensions(self):
fileh = h5py.File(self.parameter_filename,'r')
- L_index = ((na.array(fileh['/level_0'].attrs['prob_domain'].tolist()))[0:3])
- R_index = ((na.array(fileh['/level_0'].attrs['prob_domain'].tolist()))[3:] + 1)
+ L_index = ((na.array(list(fileh['/level_0'].attrs['prob_domain'])))[0:3])
+ R_index = ((na.array(list(fileh['/level_0'].attrs['prob_domain'])))[3:] + 1)
return R_index - L_index
@classmethod
https://bitbucket.org/yt_analysis/yt/changeset/96b906cea1ae/
changeset: 96b906cea1ae
branch: yt
user: acunning
date: 2012-03-12 10:43:29
summary: a fix for non-constant refinment ratios
affected #: 2 files
diff -r 9d5048db0b22e3a10b92af120e691e3c193dd548 -r 96b906cea1ae2e5dc57327610ac3c581746966a6 yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -364,7 +364,7 @@
#@time_execution
def __fill_child_mask(self, child, mask, tofill):
- rf = self.pf.refine_by
+ rf = self.pf.refine_by[child.Level-1]
gi, cgi = self.get_global_startindex(), child.get_global_startindex()
startIndex = na.maximum(0, cgi/rf - gi)
endIndex = na.minimum( (cgi+child.ActiveDimensions)/rf - gi,
diff -r 9d5048db0b22e3a10b92af120e691e3c193dd548 -r 96b906cea1ae2e5dc57327610ac3c581746966a6 yt/frontends/chombo/data_structures.py
--- a/yt/frontends/chombo/data_structures.py
+++ b/yt/frontends/chombo/data_structures.py
@@ -85,7 +85,8 @@
pdx = self.Parent[0].dds
start_index = (self.Parent[0].get_global_startindex()) + \
na.rint((self.LeftEdge - self.Parent[0].LeftEdge)/pdx)
- self.start_index = (start_index*self.pf.refine_by).astype('int64').ravel()
+ self.start_index = (start_index*self.pf.refine_by[self.Level-1]).astype('int64').ravel()
+ print self.start_index
return self.start_index
def _setup_dx(self):
https://bitbucket.org/yt_analysis/yt/changeset/e64d2b1a04ae/
changeset: e64d2b1a04ae
branch: yt
user: Andrew Myers
date: 2012-03-12 19:27:09
summary: merging
affected #: 4 files
diff -r afd860cd46d8dfca26bf7aa7a33049b654cd79c5 -r e64d2b1a04aea4dc2aef24622a916cfcdf6ab6bf yt/frontends/chombo/data_structures.py
--- a/yt/frontends/chombo/data_structures.py
+++ b/yt/frontends/chombo/data_structures.py
@@ -68,10 +68,26 @@
self.Parent = []
self.Children = []
self.Level = level
- self.start_index = start.copy()#.transpose()
- self.stop_index = stop.copy()#.transpose()
self.ActiveDimensions = stop - start + 1
+ def get_global_startindex(self):
+ """
+ Return the integer starting index for each dimension at the current
+ level.
+
+ """
+ if self.start_index != None:
+ return self.start_index
+ if self.Parent == []:
+ iLE = self.LeftEdge - self.pf.domain_left_edge
+ start_index = iLE / self.dds
+ return na.rint(start_index).astype('int64').ravel()
+ pdx = self.Parent[0].dds
+ start_index = (self.Parent[0].get_global_startindex()) + \
+ na.rint((self.LeftEdge - self.Parent[0].LeftEdge)/pdx)
+ self.start_index = (start_index*self.pf.refine_by).astype('int64').ravel()
+ return self.start_index
+
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.
@@ -91,8 +107,8 @@
grid = ChomboGrid
def __init__(self,pf,data_style='chombo_hdf5'):
- self.domain_left_edge = pf.domain_left_edge # need these to determine absolute grid locations
- self.domain_right_edge = pf.domain_right_edge # need these to determine absolute grid locations
+ self.domain_left_edge = pf.domain_left_edge
+ self.domain_right_edge = pf.domain_right_edge
self.data_style = data_style
self.field_indexes = {}
self.parameter_file = weakref.proxy(pf)
@@ -100,12 +116,11 @@
self.hierarchy_filename = self.parameter_file.parameter_filename
self.hierarchy = os.path.abspath(self.hierarchy_filename)
self.directory = os.path.dirname(self.hierarchy_filename)
- self._fhandle = h5py.File(self.hierarchy_filename)
+ self._fhandle = h5py.File(self.hierarchy_filename, 'r')
self.float_type = self._fhandle['/level_0']['data:datatype=0'].dtype.name
- self._levels = [fn for fn in self._fhandle if fn != "Chombo_global"]
+ self._levels = self._fhandle.keys()[1:]
AMRHierarchy.__init__(self,pf,data_style)
-
self._fhandle.close()
def _initialize_data_storage(self):
@@ -113,7 +128,7 @@
def _detect_fields(self):
ncomp = int(self._fhandle['/'].attrs['num_components'])
- self.field_list = [c[1] for c in self._fhandle['/'].attrs.listitems()[-ncomp:]]
+ self.field_list = [c[1] for c in self._fhandle['/'].attrs.items()[-ncomp:]]
def _setup_classes(self):
dd = self._get_data_reader_dict()
@@ -130,8 +145,8 @@
# this relies on the first Group in the H5 file being
# 'Chombo_global'
- levels = [fn for fn in f if fn != "Chombo_global"]
- self.grids = []
+ levels = f.keys()[1:]
+ grids = []
i = 0
for lev in levels:
level_number = int(re.match('level_(\d+)',lev).groups()[0])
@@ -140,17 +155,18 @@
for level_id, box in enumerate(boxes):
si = na.array([box['lo_%s' % ax] for ax in 'ijk'])
ei = na.array([box['hi_%s' % ax] for ax in 'ijk'])
- pg = self.grid(len(self.grids),self,level=level_number,
+ pg = self.grid(len(grids),self,level=level_number,
start = si, stop = ei)
- self.grids.append(pg)
- self.grids[-1]._level_id = level_id
- self.grid_left_edge[i] = dx*si.astype(self.float_type) + self.domain_left_edge
- self.grid_right_edge[i] = dx*(ei.astype(self.float_type)+1) + self.domain_left_edge
+ grids.append(pg)
+ grids[-1]._level_id = level_id
+ self.grid_left_edge[i] = dx*si.astype(self.float_type)
+ self.grid_right_edge[i] = dx*(ei.astype(self.float_type)+1)
self.grid_particle_count[i] = 0
self.grid_dimensions[i] = ei - si + 1
i += 1
self.grids = na.empty(len(grids), dtype='object')
for gi, g in enumerate(grids): self.grids[gi] = g
+# self.grids = na.array(self.grids, dtype='object')
def _populate_grid_objects(self):
for g in self.grids:
@@ -179,8 +195,8 @@
def __init__(self, filename, data_style='chombo_hdf5',
storage_filename = None, ini_filename = None):
- # hardcoded for now
- self.current_time = 0.0
+ fileh = h5py.File(filename,'r')
+ self.current_time = fileh.attrs['time']
self.ini_filename = ini_filename
StaticOutput.__init__(self,filename,data_style)
self.storage_filename = storage_filename
@@ -230,15 +246,19 @@
"""
if os.path.isfile('pluto.ini'):
self._parse_pluto_file('pluto.ini')
- elif os.path.isfile('orion2.ini'):
- self._parse_pluto_file('orion2.ini')
else:
+ if os.path.isfile('orion2.ini'): self._parse_pluto_file('orion2.ini')
self.unique_identifier = \
- int(os.stat(self.parameter_filename)[ST_CTIME])
- self.domain_left_edge = na.array([0.,0.,0.])
+ int(os.stat(self.parameter_filename)[ST_CTIME])
+ self.domain_left_edge = self.__calc_left_edge()
self.domain_right_edge = self.__calc_right_edge()
+ self.domain_dimensions = self.__calc_domain_dimensions()
self.dimensionality = 3
- self.refine_by = 2
+# self.refine_by = []
+ fileh = h5py.File(self.parameter_filename,'r')
+# for level in range(0,fileh.attrs['max_level']):
+# self.refine_by.append(fileh['/level_'+str(level)].attrs['ref_ratio'])
+ self.refine_by = fileh['/level_0'].attrs['ref_ratio']
def _parse_pluto_file(self, ini_filename):
"""
@@ -268,36 +288,26 @@
else:
self.parameters[paramName] = t
- # assumes 3D for now
- elif param.startswith("X1-grid"):
- t = vals.split()
- low1 = float(t[1])
- high1 = float(t[4])
- N1 = int(t[2])
- elif param.startswith("X2-grid"):
- t = vals.split()
- low2 = float(t[1])
- high2 = float(t[4])
- N2 = int(t[2])
- elif param.startswith("X3-grid"):
- t = vals.split()
- low3 = float(t[1])
- high3 = float(t[4])
- N3 = int(t[2])
+ def __calc_left_edge(self):
+ fileh = h5py.File(self.parameter_filename,'r')
+ dx0 = fileh['/level_0'].attrs['dx']
+ LE = dx0*((na.array(list(fileh['/level_0'].attrs['prob_domain'])))[0:3])
+ fileh.close()
+ return LE
- self.dimensionality = 3
- self.domain_left_edge = na.array([low1,low2,low3])
- self.domain_right_edge = na.array([high1,high2,high3])
- self.domain_dimensions = na.array([N1,N2,N3])
- self.refine_by = self.parameters["RefineBy"]
-
def __calc_right_edge(self):
fileh = h5py.File(self.parameter_filename,'r')
dx0 = fileh['/level_0'].attrs['dx']
- RE = dx0*((na.array(fileh['/level_0'].attrs['prob_domain']))[3:] + 1)
+ RE = dx0*((na.array(list(fileh['/level_0'].attrs['prob_domain'])))[3:] + 1)
fileh.close()
return RE
-
+
+ def __calc_domain_dimensions(self):
+ fileh = h5py.File(self.parameter_filename,'r')
+ L_index = ((na.array(list(fileh['/level_0'].attrs['prob_domain'])))[0:3])
+ R_index = ((na.array(list(fileh['/level_0'].attrs['prob_domain'])))[3:] + 1)
+ return R_index - L_index
+
@classmethod
def _is_valid(self, *args, **kwargs):
try:
@@ -309,7 +319,6 @@
pass
return False
-
@parallel_root_only
def print_key_parameters(self):
for a in ["current_time", "domain_dimensions", "domain_left_edge",
diff -r afd860cd46d8dfca26bf7aa7a33049b654cd79c5 -r e64d2b1a04aea4dc2aef24622a916cfcdf6ab6bf yt/frontends/chombo/fields.py
--- a/yt/frontends/chombo/fields.py
+++ b/yt/frontends/chombo/fields.py
@@ -38,45 +38,64 @@
add_chombo_field = KnownChomboFields.add_field
ChomboFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
-add_chombo_field = ChomboFieldInfo.add_field
+add_field = ChomboFieldInfo.add_field
-add_field = add_chombo_field
+add_chombo_field("density", function=NullFunc, take_log=True,
+ validators = [ValidateDataField("density")],
+ units=r"\rm{g}/\rm{cm}^3")
-add_field("density", function=NullFunc, take_log=True,
- validators = [ValidateDataField("density")],
- units=r"\rm{g}/\rm{cm}^3")
+KnownChomboFields["density"]._projected_units =r"\rm{g}/\rm{cm}^2"
-ChomboFieldInfo["density"]._projected_units =r"\rm{g}/\rm{cm}^2"
+add_chombo_field("X-momentum", function=NullFunc, take_log=False,
+ validators = [ValidateDataField("X-Momentum")],
+ units=r"",display_name=r"M_x")
+KnownChomboFields["X-momentum"]._projected_units=r""
-add_field("X-momentum", function=NullFunc, take_log=False,
- validators = [ValidateDataField("X-Momentum")],
- units=r"",display_name=r"B_x")
-ChomboFieldInfo["X-momentum"]._projected_units=r""
+add_chombo_field("Y-momentum", function=NullFunc, take_log=False,
+ validators = [ValidateDataField("Y-Momentum")],
+ units=r"",display_name=r"M_y")
+KnownChomboFields["Y-momentum"]._projected_units=r""
-add_field("Y-momentum", function=NullFunc, take_log=False,
- validators = [ValidateDataField("Y-Momentum")],
- units=r"",display_name=r"B_y")
-ChomboFieldInfo["Y-momentum"]._projected_units=r""
+add_chombo_field("Z-momentum", function=NullFunc, take_log=False,
+ validators = [ValidateDataField("Z-Momentum")],
+ units=r"",display_name=r"M_z")
+KnownChomboFields["Z-momentum"]._projected_units=r""
-add_field("Z-momentum", function=NullFunc, take_log=False,
- validators = [ValidateDataField("Z-Momentum")],
- units=r"",display_name=r"B_z")
-ChomboFieldInfo["Z-momentum"]._projected_units=r""
+add_chombo_field("X-magnfield", function=NullFunc, take_log=False,
+ validators = [ValidateDataField("X-Magnfield")],
+ units=r"",display_name=r"B_x")
+KnownChomboFields["X-magnfield"]._projected_units=r""
-add_field("X-magnfield", function=NullFunc, take_log=False,
- validators = [ValidateDataField("X-Magnfield")],
- units=r"",display_name=r"B_x")
-ChomboFieldInfo["X-magnfield"]._projected_units=r""
+add_chombo_field("Y-magnfield", function=NullFunc, take_log=False,
+ validators = [ValidateDataField("Y-Magnfield")],
+ units=r"",display_name=r"B_y")
+KnownChomboFields["Y-magnfield"]._projected_units=r""
-add_field("Y-magnfield", function=NullFunc, take_log=False,
- validators = [ValidateDataField("Y-Magnfield")],
- units=r"",display_name=r"B_y")
-ChomboFieldInfo["Y-magnfield"]._projected_units=r""
+add_chombo_field("Z-magnfield", function=NullFunc, take_log=False,
+ validators = [ValidateDataField("Z-Magnfield")],
+ units=r"",display_name=r"B_z")
+KnownChomboFields["Z-magnfield"]._projected_units=r""
-add_field("Z-magnfield", function=NullFunc, take_log=False,
- validators = [ValidateDataField("Z-Magnfield")],
- units=r"",display_name=r"B_z")
-ChomboFieldInfo["Z-magnfield"]._projected_units=r""
+add_chombo_field("energy-density", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("energy-density")],
+ units=r"\rm{erg}/\rm{cm}^3")
+KnownChomboFields["energy-density"]._projected_units =r""
+
+add_chombo_field("radiation-energy-density", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("radiation-energy-density")],
+ units=r"\rm{erg}/\rm{cm}^3")
+KnownChomboFields["radiation-energy-density"]._projected_units =r""
+
+def _Density(field,data):
+ """A duplicate of the density field. This is needed because when you try
+ to instantiate a PlotCollection without passing in a center, the code
+ will try to generate one for you using the "Density" field, which gives an error
+ if it isn't defined.
+
+ """
+ return data["density"]
+add_field("Density",function=_Density, take_log=True,
+ units=r'\rm{g}/\rm{cm^3}')
def _MagneticEnergy(field,data):
return (data["X-magnfield"]**2 +
diff -r afd860cd46d8dfca26bf7aa7a33049b654cd79c5 -r e64d2b1a04aea4dc2aef24622a916cfcdf6ab6bf yt/frontends/chombo/io.py
--- a/yt/frontends/chombo/io.py
+++ b/yt/frontends/chombo/io.py
@@ -36,7 +36,7 @@
def _field_dict(self,fhandle):
ncomp = int(fhandle['/'].attrs['num_components'])
- temp = fhandle['/'].attrs.listitems()[-ncomp:]
+ temp = fhandle['/'].attrs.items()[-ncomp:]
val, keys = zip(*temp)
val = [int(re.match('component_(\d+)',v).groups()[0]) for v in val]
return dict(zip(keys,val))
@@ -45,7 +45,7 @@
fhandle = h5py.File(grid.filename,'r')
ncomp = int(fhandle['/'].attrs['num_components'])
- fns = [c[1] for c in f['/'].attrs.listitems()[-ncomp:]]
+ fns = [c[1] for c in f['/'].attrs.items()[-ncomp-1:-1]]
fhandle.close()
def _read_data_set(self,grid,field):
@@ -64,7 +64,6 @@
fhandle.close()
return data.reshape(dims, order='F')
-
def _read_data_slice(self, grid, field, axis, coord):
sl = [slice(None), slice(None), slice(None)]
diff -r afd860cd46d8dfca26bf7aa7a33049b654cd79c5 -r e64d2b1a04aea4dc2aef24622a916cfcdf6ab6bf yt/frontends/orion/fields.py
--- a/yt/frontends/orion/fields.py
+++ b/yt/frontends/orion/fields.py
@@ -44,26 +44,26 @@
OrionFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
add_field = OrionFieldInfo.add_field
-add_field("density", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("density")],
- units=r"\rm{g}/\rm{cm}^3")
-OrionFieldInfo["density"]._projected_units =r"\rm{g}/\rm{cm}^2"
+add_orion_field("density", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("density")],
+ units=r"\rm{g}/\rm{cm}^3")
+KnownOrionFields["density"]._projected_units =r"\rm{g}/\rm{cm}^2"
-add_field("eden", function=lambda a,b: None, take_log=True,
- validators = [ValidateDataField("eden")],
- units=r"\rm{erg}/\rm{cm}^3")
+add_orion_field("eden", function=lambda a,b: None, take_log=True,
+ validators = [ValidateDataField("eden")],
+ units=r"\rm{erg}/\rm{cm}^3")
-add_field("xmom", function=lambda a,b: None, take_log=False,
- validators = [ValidateDataField("xmom")],
- units=r"\rm{g}/\rm{cm^2\ s}")
+add_orion_field("xmom", function=lambda a,b: None, take_log=False,
+ validators = [ValidateDataField("xmom")],
+ units=r"\rm{g}/\rm{cm^2\ s}")
-add_field("ymom", function=lambda a,b: None, take_log=False,
- validators = [ValidateDataField("ymom")],
- units=r"\rm{gm}/\rm{cm^2\ s}")
+add_orion_field("ymom", function=lambda a,b: None, take_log=False,
+ validators = [ValidateDataField("ymom")],
+ units=r"\rm{gm}/\rm{cm^2\ s}")
-add_field("zmom", function=lambda a,b: None, take_log=False,
- validators = [ValidateDataField("zmom")],
- units=r"\rm{g}/\rm{cm^2\ s}")
+add_orion_field("zmom", function=lambda a,b: None, take_log=False,
+ validators = [ValidateDataField("zmom")],
+ units=r"\rm{g}/\rm{cm^2\ s}")
translation_dict = {"x-velocity": "xvel",
"y-velocity": "yvel",
@@ -88,11 +88,11 @@
def _xVelocity(field, data):
"""generate x-velocity from x-momentum and density
-
+
"""
return data["xmom"]/data["density"]
-add_field("x-velocity",function=_xVelocity, take_log=False,
- units=r'\rm{cm}/\rm{s}')
+add_orion_field("x-velocity",function=_xVelocity, take_log=False,
+ units=r'\rm{cm}/\rm{s}')
def _yVelocity(field,data):
"""generate y-velocity from y-momentum and density
@@ -102,16 +102,16 @@
# return data["xvel"]
#except KeyError:
return data["ymom"]/data["density"]
-add_field("y-velocity",function=_yVelocity, take_log=False,
- units=r'\rm{cm}/\rm{s}')
+add_orion_field("y-velocity",function=_yVelocity, take_log=False,
+ units=r'\rm{cm}/\rm{s}')
def _zVelocity(field,data):
"""generate z-velocity from z-momentum and density
-
+
"""
return data["zmom"]/data["density"]
-add_field("z-velocity",function=_zVelocity, take_log=False,
- units=r'\rm{cm}/\rm{s}')
+add_orion_field("z-velocity",function=_zVelocity, take_log=False,
+ units=r'\rm{cm}/\rm{s}')
def _ThermalEnergy(field, data):
"""generate thermal (gas energy). Dual Energy Formalism was
@@ -125,19 +125,19 @@
data["x-velocity"]**2.0
+ data["y-velocity"]**2.0
+ data["z-velocity"]**2.0 )
-add_field("ThermalEnergy", function=_ThermalEnergy,
- units=r"\rm{ergs}/\rm{cm^3}")
+add_orion_field("ThermalEnergy", function=_ThermalEnergy,
+ units=r"\rm{ergs}/\rm{cm^3}")
def _Pressure(field,data):
"""M{(Gamma-1.0)*e, where e is thermal energy density
NB: this will need to be modified for radiation
"""
return (data.pf["Gamma"] - 1.0)*data["ThermalEnergy"]
-add_field("Pressure", function=_Pressure, units=r"\rm{dyne}/\rm{cm}^{2}")
+add_orion_field("Pressure", function=_Pressure, units=r"\rm{dyne}/\rm{cm}^{2}")
def _Temperature(field,data):
return (data.pf["Gamma"]-1.0)*data.pf["mu"]*mh*data["ThermalEnergy"]/(kboltz*data["Density"])
-add_field("Temperature",function=_Temperature,units=r"\rm{Kelvin}",take_log=False)
+add_orion_field("Temperature",function=_Temperature,units=r"\rm{Kelvin}",take_log=False)
# particle fields
@@ -170,6 +170,6 @@
for pf in _particle_field_list:
pfunc = particle_func("particle_%s" % (pf))
- add_field("particle_%s" % pf, function=pfunc,
- validators = [ValidateSpatial(0)],
- particle_type=True)
+ add_orion_field("particle_%s" % pf, function=pfunc,
+ validators = [ValidateSpatial(0)],
+ particle_type=True)
https://bitbucket.org/yt_analysis/yt/changeset/9610fc389e61/
changeset: 9610fc389e61
branch: yt
user: Andrew Myers
date: 2012-03-12 19:41:48
summary: merging
affected #: 5 files
diff -r e64d2b1a04aea4dc2aef24622a916cfcdf6ab6bf -r 9610fc389e6188cb6e9c620febca1e2dc18c711f yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -367,7 +367,7 @@
#@time_execution
def __fill_child_mask(self, child, mask, tofill):
- rf = self.pf.refine_by
+ rf = self.pf.refine_by[child.Level-1]
gi, cgi = self.get_global_startindex(), child.get_global_startindex()
startIndex = na.maximum(0, cgi / rf - gi)
endIndex = na.minimum((cgi + child.ActiveDimensions) / rf - gi,
diff -r e64d2b1a04aea4dc2aef24622a916cfcdf6ab6bf -r 9610fc389e6188cb6e9c620febca1e2dc18c711f yt/frontends/chombo/data_structures.py
--- a/yt/frontends/chombo/data_structures.py
+++ b/yt/frontends/chombo/data_structures.py
@@ -85,7 +85,7 @@
pdx = self.Parent[0].dds
start_index = (self.Parent[0].get_global_startindex()) + \
na.rint((self.LeftEdge - self.Parent[0].LeftEdge)/pdx)
- self.start_index = (start_index*self.pf.refine_by).astype('int64').ravel()
+ self.start_index = (start_index*self.pf.refine_by[self.Level-1]).astype('int64').ravel()
return self.start_index
def _setup_dx(self):
@@ -93,7 +93,7 @@
# 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
+ self.dds = self.Parent[0].dds / self.pf.refine_by[self.Level-1]
else:
LE, RE = self.hierarchy.grid_left_edge[id,:], \
self.hierarchy.grid_right_edge[id,:]
@@ -254,11 +254,10 @@
self.domain_right_edge = self.__calc_right_edge()
self.domain_dimensions = self.__calc_domain_dimensions()
self.dimensionality = 3
-# self.refine_by = []
+ self.refine_by = []
fileh = h5py.File(self.parameter_filename,'r')
-# for level in range(0,fileh.attrs['max_level']):
-# self.refine_by.append(fileh['/level_'+str(level)].attrs['ref_ratio'])
- self.refine_by = fileh['/level_0'].attrs['ref_ratio']
+ for level in range(0,fileh.attrs['max_level']):
+ self.refine_by.append(fileh['/level_'+str(level)].attrs['ref_ratio'])
def _parse_pluto_file(self, ini_filename):
"""
https://bitbucket.org/yt_analysis/yt/changeset/75acdcb27757/
changeset: 75acdcb27757
branch: yt
user: Andrew Myers
date: 2012-03-12 21:19:23
summary: one more fix to the chombo reader, for datasets where the max level is not present in the data dump
affected #: 2 files
diff -r 9610fc389e6188cb6e9c620febca1e2dc18c711f -r 75acdcb27757dcecbe11759cc5c409dfa5781de2 yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -367,7 +367,10 @@
#@time_execution
def __fill_child_mask(self, child, mask, tofill):
- rf = self.pf.refine_by[child.Level-1]
+ try:
+ rf = self.pf.refine_by[child.Level-1]
+ except TypeError:
+ rf = self.pf.refine_by
gi, cgi = self.get_global_startindex(), child.get_global_startindex()
startIndex = na.maximum(0, cgi / rf - gi)
endIndex = na.minimum((cgi + child.ActiveDimensions) / rf - gi,
diff -r 9610fc389e6188cb6e9c620febca1e2dc18c711f -r 75acdcb27757dcecbe11759cc5c409dfa5781de2 yt/frontends/chombo/data_structures.py
--- a/yt/frontends/chombo/data_structures.py
+++ b/yt/frontends/chombo/data_structures.py
@@ -256,7 +256,7 @@
self.dimensionality = 3
self.refine_by = []
fileh = h5py.File(self.parameter_filename,'r')
- for level in range(0,fileh.attrs['max_level']):
+ for level in range(0,fileh.attrs['num_levels']):
self.refine_by.append(fileh['/level_'+str(level)].attrs['ref_ratio'])
def _parse_pluto_file(self, ini_filename):
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