[yt-svn] commit/yt-3.0: 4 new changesets
Bitbucket
commits-noreply at bitbucket.org
Tue Jan 15 12:03:16 PST 2013
4 new commits in yt-3.0:
https://bitbucket.org/yt_analysis/yt-3.0/commits/10fb5dc22ecb/
changeset: 10fb5dc22ecb
branch: yt-3.0
user: MatthewTurk
date: 2013-01-11 14:32:23
summary: Adding a read_record helper to fortran_utils.
affected #: 1 file
diff -r 24cd112f134bde038875e3bdd2fb9c270b8f126c -r 10fb5dc22ecb8b14ad062761d220a8d65226a118 yt/utilities/fortran_utils.py
--- a/yt/utilities/fortran_utils.py
+++ b/yt/utilities/fortran_utils.py
@@ -33,7 +33,9 @@
Fortran unformatted files provide total bytesize at the beginning and end
of a record. By correlating the components of that record with attribute
- names, we construct a dictionary that gets returned.
+ names, we construct a dictionary that gets returned. Note that this
+ function is used for reading sequentially-written records. If you have
+ many written that were written simultaneously, see read_record.
Parameters
----------
@@ -133,3 +135,47 @@
ss = struct.unpack(fmt, f.read(struct.calcsize(fmt)))[0]
f.seek(ss + struct.calcsize("=I"), os.SEEK_CUR)
+def read_record(f, rspec):
+ r"""This function accepts a file pointer and reads from that file pointer
+ a single "record" with different components.
+
+ Fortran unformatted files provide total bytesize at the beginning and end
+ of a record. By correlating the components of that record with attribute
+ names, we construct a dictionary that gets returned.
+
+ Parameters
+ ----------
+ f : File object
+ An open file object. Should have been opened in mode rb.
+ rspec : iterable of iterables
+ This object should be an iterable of the format [ (attr_name, count,
+ struct type), ... ].
+
+ Returns
+ -------
+ values : dict
+ This will return a dict of iterables of the components of the values in
+ the file.
+
+ Examples
+ --------
+
+ >>> header = [ ("ncpu", 1, "i"), ("nfiles", 2, "i") ]
+ >>> f = open("fort.3", "rb")
+ >>> rv = read_record(f, header)
+ """
+ vv = {}
+ net_format = "=I" + "".join(["%s%s" % (n, t) for a, n, t in rspec]) + "I"
+ size = struct.calcsize(net_format)
+ vals = list(struct.unpack(net_format, f.read(size)))
+ vvv = vals[:]
+ s1, s2 = vals.pop(0), vals.pop(-1)
+ if s1 != s2:
+ print "S1 = %s ; S2 = %s ; SIZE = %s"
+ raise RuntimeError
+ pos = 0
+ for a, n, t in rspec:
+ vv[a] = vals[pos:pos+n]
+ pos += n
+ return vv, vvv
+
https://bitbucket.org/yt_analysis/yt-3.0/commits/0ef03c5db632/
changeset: 0ef03c5db632
branch: yt-3.0
user: MatthewTurk
date: 2013-01-11 19:50:18
summary: Initial import of Binary Gadget reader. Field IO still TBD.
affected #: 5 files
diff -r 10fb5dc22ecb8b14ad062761d220a8d65226a118 -r 0ef03c5db632882274f68dca0f259ca95a21e393 yt/frontends/sph/api.py
--- a/yt/frontends/sph/api.py
+++ b/yt/frontends/sph/api.py
@@ -29,7 +29,9 @@
"""
from .data_structures import \
- OWLSStaticOutput
+ OWLSStaticOutput, \
+ GadgetStaticOutput
from .io import \
- IOHandlerOWLS
+ IOHandlerOWLS, \
+ IOHandlerGadgetBinary
diff -r 10fb5dc22ecb8b14ad062761d220a8d65226a118 -r 0ef03c5db632882274f68dca0f259ca95a21e393 yt/frontends/sph/data_structures.py
--- a/yt/frontends/sph/data_structures.py
+++ b/yt/frontends/sph/data_structures.py
@@ -29,6 +29,7 @@
import weakref
from itertools import izip
+from yt.utilities.fortran_utils import read_record
from yt.funcs import *
from yt.geometry.oct_geometry_handler import \
OctreeGeometryHandler
@@ -42,7 +43,9 @@
mpc_conversion, sec_conversion
from .fields import \
OWLSFieldInfo, \
- KnownOWLSFields
+ KnownOWLSFields, \
+ GadgetFieldInfo, \
+ KnownGadgetFields
from yt.data_objects.field_info_container import \
FieldInfoContainer, NullFunc
@@ -53,7 +56,7 @@
self.io = weakref.proxy(io)
self.domain_filename = domain_filename
self.domain_number = domain_number
- self.total_particles = self.io._count_particles(domain_filename)
+ self.total_particles = self.io._count_particles(self)
def select(self, selector):
pass
@@ -108,9 +111,11 @@
self._setup_data_io()
template = self.parameter_file.domain_template
ndoms = self.parameter_file.domain_count
- self.domains = [ParticleDomainFile(self.parameter_file, self.io, template % i, i)
+ cls = self.parameter_file._domain_class
+ self.domains = [cls(self.parameter_file, self.io, template % {'num':i}, i)
for i in range(ndoms)]
- total_particles = sum(d.total_particles for d in self.domains)
+ total_particles = sum(sum(d.total_particles.values())
+ for d in self.domains)
self.oct_handler = ParticleOctreeContainer(
self.parameter_file.domain_dimensions,
self.parameter_file.domain_left_edge,
@@ -126,7 +131,7 @@
# TODO: Add additional fields
pfl = []
for dom in self.domains:
- fl = self.io._identify_fields(dom.domain_filename)
+ fl = self.io._identify_fields(dom)
for f in fl:
if f not in pfl: pfl.append(f)
self.field_list = pfl
@@ -163,6 +168,7 @@
class OWLSStaticOutput(StaticOutput):
_hierarchy_class = ParticleGeometryHandler
+ _domain_class = ParticleDomainFile
_fieldinfo_fallback = OWLSFieldInfo
_fieldinfo_known = KnownOWLSFields
@@ -204,7 +210,7 @@
prefix = self.parameter_filename.split(".", 1)[0]
suffix = self.parameter_filename.rsplit(".", 1)[-1]
- self.domain_template = "%s.%%i.%s" % (prefix, suffix)
+ self.domain_template = "%s.%%(num)i.%s" % (prefix, suffix)
self.domain_count = hvals["NumFilesPerSnapshot"]
handle.close()
@@ -221,3 +227,99 @@
except:
pass
return False
+
+class GadgetBinaryDomainFile(ParticleDomainFile):
+ def __init__(self, pf, io, domain_filename, domain_number):
+ with open(domain_filename, "rb") as f:
+ self.header = read_record(f, pf._header_spec)
+
+ super(GadgetBinaryDomainFile, self).__init__(pf, io,
+ domain_filename, domain_number)
+
+class GadgetStaticOutput(StaticOutput):
+ _hierarchy_class = ParticleGeometryHandler
+ _domain_class = GadgetBinaryDomainFile
+ _fieldinfo_fallback = GadgetFieldInfo
+ _fieldinfo_known = KnownGadgetFields
+ _header_spec = (('Npart', 6, 'i'),
+ ('Massarr', 6, 'd'),
+ ('Time', 1, 'd'),
+ ('Redshift', 1, 'd'),
+ ('FlagSfr', 1, 'i'),
+ ('FlagFeedback', 1, 'i'),
+ ('Nall', 6, 'i'),
+ ('FlagCooling', 1, 'i'),
+ ('NumFiles', 1, 'i'),
+ ('BoxSize', 1, 'd'),
+ ('Omega0', 1, 'd'),
+ ('OmegaLambda', 1, 'd'),
+ ('HubbleParam', 1, 'd'),
+ ('FlagAge', 1, 'i'),
+ ('FlagMEtals', 1, 'i'),
+ ('NallHW', 6, 'i'),
+ ('unused', 16, 'i') )
+
+ def __init__(self, filename, data_style="gadget_binary",
+ additional_fields = (), root_dimensions = 64):
+ self._root_dimensions = root_dimensions
+ # Set up the template for domain files
+ self.storage_filename = None
+ super(GadgetStaticOutput, self).__init__(filename, data_style)
+
+ def __repr__(self):
+ return os.path.basename(self.parameter_filename).split(".")[0]
+
+ def _set_units(self):
+ self.units = {}
+ self.time_units = {}
+ self.conversion_factors = {}
+
+ def _parse_parameter_file(self):
+
+ # The entries in this header are capitalized and named to match Table 4
+ # in the GADGET-2 user guide.
+
+ f = open(self.parameter_filename)
+ hvals = read_record(f, self._header_spec)
+ for i in hvals:
+ if len(hvals[i]) == 1:
+ hvals[i] = hvals[i][0]
+
+ self.dimensionality = 3
+ self.refine_by = 2
+ self.parameters["HydroMethod"] = "sph"
+ self.unique_identifier = \
+ int(os.stat(self.parameter_filename)[stat.ST_CTIME])
+ # Set standard values
+
+ # This may not be correct.
+ self.current_time = hvals["Time"] * sec_conversion["Gyr"]
+
+ self.domain_left_edge = np.zeros(3, "float64")
+ self.domain_right_edge = np.ones(3, "float64") * hvals["BoxSize"]
+ self.domain_dimensions = np.ones(3, "int32") * self._root_dimensions
+
+ self.cosmological_simulation = 1
+
+ self.current_redshift = hvals["Redshift"]
+ self.omega_lambda = hvals["OmegaLambda"]
+ self.omega_matter = hvals["Omega0"]
+ self.hubble_constant = hvals["HubbleParam"]
+ self.parameters = hvals
+
+ prefix = self.parameter_filename.split(".", 1)[0]
+ suffix = self.parameter_filename.rsplit(".", 1)[-1]
+
+ if hvals["NumFiles"] > 1:
+ self.domain_template = "%s.%%(num)s" % (prefix)
+ else:
+ self.domain_template = self.parameter_filename
+
+ self.domain_count = hvals["NumFiles"]
+
+ f.close()
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ # We do not allow load() of these files.
+ return False
diff -r 10fb5dc22ecb8b14ad062761d220a8d65226a118 -r 0ef03c5db632882274f68dca0f259ca95a21e393 yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -42,3 +42,9 @@
KnownOWLSFields = FieldInfoContainer()
add_OWLS_field = KnownOWLSFields.add_field
+GadgetFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo)
+add_Gadget_field = GadgetFieldInfo.add_field
+
+KnownGadgetFields = FieldInfoContainer()
+add_Gadget_field = KnownGadgetFields.add_field
+
diff -r 10fb5dc22ecb8b14ad062761d220a8d65226a118 -r 0ef03c5db632882274f68dca0f259ca95a21e393 yt/frontends/sph/io.py
--- a/yt/frontends/sph/io.py
+++ b/yt/frontends/sph/io.py
@@ -31,6 +31,8 @@
from yt.utilities.io_handler import \
BaseIOHandler
+from yt.utilities.fortran_utils import read_record
+
_vector_fields = ("Coordinates", "Velocity")
class IOHandlerOWLS(BaseIOHandler):
@@ -49,13 +51,13 @@
ptf[ftype].append(fname)
for chunk in chunks: # Will be OWLS domains
for subset in chunk.objs:
+ f = h5py.File(subset.domain.domain_filename, "r")
for ptype, field_list in sorted(ptf.items()):
- f = h5py.File(subset.domain.domain_filename, "r")
coords = f["/%s/Coordinates" % ptype][:].astype("float64")
psize[ptype] += selector.count_points(
coords[:,0], coords[:,1], coords[:,2])
del coords
- f.close()
+ f.close()
# Now we have all the sizes, and we can allocate
ind = {}
for field in fields:
@@ -68,8 +70,8 @@
ind[field] = 0
for chunk in chunks: # Will be OWLS domains
for subset in chunk.objs:
+ f = h5py.File(subset.domain.domain_filename, "r")
for ptype, field_list in sorted(ptf.items()):
- f = h5py.File(subset.domain.domain_filename, "r")
g = f["/%s" % ptype]
coords = g["Coordinates"][:].astype("float64")
mask = selector.select_points(
@@ -83,7 +85,7 @@
my_ind, my_ind+data.shape[0], field)
rv[ptype, field][my_ind:my_ind + data.shape[0],...] = data
ind[ptype, field] += data.shape[0]
- f.close()
+ f.close()
return rv
def _initialize_octree(self, domain, octree):
@@ -94,14 +96,15 @@
octree.add(pos, domain.domain_number)
f.close()
- def _count_particles(self, domain_filename):
- f = h5py.File(domain_filename, "r")
- npart = f["/Header"].attrs["NumPart_ThisFile"].sum()
+ def _count_particles(self, domain):
+ f = h5py.File(domain.domain_filename, "r")
+ np = f["/Header"].attrs["NumPart_ThisFile"][:]
f.close()
+ npart = dict(("PartType%s" % (i), v) for i, v in enumerate(np))
return npart
- def _identify_fields(self, domain_filename):
- f = h5py.File(domain_filename, "r")
+ def _identify_fields(self, domain):
+ f = h5py.File(domain.domain_filename, "r")
fields = []
for key in f.keys():
if not key.startswith("PartType"): continue
@@ -114,3 +117,82 @@
fields.append((ptype, str(k)))
f.close()
return fields
+
+
+ZeroMass = object()
+
+class IOHandlerGadgetBinary(BaseIOHandler):
+ _data_style = "gadget_binary"
+
+ # Particle types (Table 3 in GADGET-2 user guide)
+ _ptypes = ( "Gas",
+ "Halo",
+ "Disk",
+ "Bulge",
+ "Stars",
+ "Bndry" )
+ #
+ # Blocks in the file:
+ # HEAD
+ # POS
+ # VEL
+ # ID
+ # MASS (variable mass only)
+ # U (gas only)
+ # RHO (gas only)
+ # HSML (gas only)
+ # POT (only if enabled in makefile)
+ # ACCE (only if enabled in makefile)
+ # ENDT (only if enabled in makefile)
+ # TSTP (only if enabled in makefile)
+
+ _fields = ( "Coordinates",
+ "Velocities",
+ "ParticleIDs",
+ ("Masses", ZeroMass),
+ ("InternalEnergy", "Gas"),
+ ("Density", "Gas"),
+ ("SmoothingLength", "Gas"),
+ )
+
+ def _read_fluid_selection(self, chunks, selector, fields, size):
+ raise NotImplementedError
+
+ def _read_particle_selection(self, chunks, selector, fields):
+ raise NotImplementedError
+
+ def _initialize_octree(self, domain, octree):
+ count = sum(domain.total_particles.values())
+ dt = [("px", "float32"), ("py", "float32"), ("pz", "float32")]
+ with open(domain.domain_filename, "rb") as f:
+ # The first total_particles * 3 values are positions
+ pp = np.fromfile(f, dtype = dt, count = count)
+ pos = np.empty((count, 3), dtype="float64")
+ pos[:,0] = pp['px']
+ pos[:,1] = pp['py']
+ pos[:,2] = pp['pz']
+ del pp
+ octree.add(pos, domain.domain_number)
+
+ def _count_particles(self, domain):
+ npart = dict((self._ptypes[i], v)
+ for i, v in enumerate(domain.header["Npart"]))
+ return npart
+
+ def _identify_fields(self, domain):
+ # We can just look at the particle counts.
+ field_list = []
+ tp = domain.total_particles
+ for i, ptype in enumerate(self._ptypes):
+ count = tp[ptype]
+ if count == 0: continue
+ m = domain.header["Massarr"][i]
+ for field in self._fields:
+ if isinstance(field, types.TupleType):
+ field, req = field
+ if req is ZeroMass:
+ if m > 0.0 : continue
+ elif req != field:
+ continue
+ field_list.append((ptype, field))
+ return field_list
diff -r 10fb5dc22ecb8b14ad062761d220a8d65226a118 -r 0ef03c5db632882274f68dca0f259ca95a21e393 yt/utilities/fortran_utils.py
--- a/yt/utilities/fortran_utils.py
+++ b/yt/utilities/fortran_utils.py
@@ -177,5 +177,5 @@
for a, n, t in rspec:
vv[a] = vals[pos:pos+n]
pos += n
- return vv, vvv
+ return vv
https://bitbucket.org/yt_analysis/yt-3.0/commits/c5ad47a08dde/
changeset: c5ad47a08dde
branch: yt-3.0
user: MatthewTurk
date: 2013-01-11 22:02:33
summary: First pass at Gadget Binary data reading.
affected #: 2 files
diff -r 0ef03c5db632882274f68dca0f259ca95a21e393 -r c5ad47a08dde5fc1f2923f838cbb843e392fb1dc yt/frontends/sph/data_structures.py
--- a/yt/frontends/sph/data_structures.py
+++ b/yt/frontends/sph/data_structures.py
@@ -64,6 +64,9 @@
def count(self, selector):
pass
+ def _calculate_offsets(self, fields):
+ pass
+
class ParticleDomainSubset(object):
def __init__(self, domain, mask, count):
self.domain = domain
@@ -132,6 +135,7 @@
pfl = []
for dom in self.domains:
fl = self.io._identify_fields(dom)
+ dom._calculate_offsets(fl)
for f in fl:
if f not in pfl: pfl.append(f)
self.field_list = pfl
@@ -232,10 +236,15 @@
def __init__(self, pf, io, domain_filename, domain_number):
with open(domain_filename, "rb") as f:
self.header = read_record(f, pf._header_spec)
+ self._position_offset = f.tell()
super(GadgetBinaryDomainFile, self).__init__(pf, io,
domain_filename, domain_number)
+ def _calculate_offsets(self, field_list):
+ self.field_offsets = self.io._calculate_field_offsets(
+ field_list, self.total_particles)
+
class GadgetStaticOutput(StaticOutput):
_hierarchy_class = ParticleGeometryHandler
_domain_class = GadgetBinaryDomainFile
diff -r 0ef03c5db632882274f68dca0f259ca95a21e393 -r c5ad47a08dde5fc1f2923f838cbb843e392fb1dc yt/frontends/sph/io.py
--- a/yt/frontends/sph/io.py
+++ b/yt/frontends/sph/io.py
@@ -33,7 +33,7 @@
from yt.utilities.fortran_utils import read_record
-_vector_fields = ("Coordinates", "Velocity")
+_vector_fields = ("Coordinates", "Velocity", "Velocities")
class IOHandlerOWLS(BaseIOHandler):
_data_style = "OWLS"
@@ -52,6 +52,7 @@
for chunk in chunks: # Will be OWLS domains
for subset in chunk.objs:
f = h5py.File(subset.domain.domain_filename, "r")
+ # This double-reads
for ptype, field_list in sorted(ptf.items()):
coords = f["/%s/Coordinates" % ptype][:].astype("float64")
psize[ptype] += selector.count_points(
@@ -159,7 +160,76 @@
raise NotImplementedError
def _read_particle_selection(self, chunks, selector, fields):
- raise NotImplementedError
+ rv = {}
+ # We first need a set of masks for each particle type
+ ptf = defaultdict(list)
+ psize = defaultdict(lambda: 0)
+ chunks = list(chunks)
+ ptypes = set()
+ for ftype, fname in fields:
+ ptf[ftype].append(fname)
+ ptypes.add(ftype)
+ ptypes = list(ptypes)
+ ptypes.sort(key = lambda a: self._ptypes.index(a))
+ for chunk in chunks:
+ for subset in chunk.objs:
+ poff = subset.domain.field_offsets
+ tp = subset.domain.total_particles
+ f = open(subset.domain.domain_filename, "rb")
+ for ptype in ptypes:
+ f.seek(poff[ptype, "Coordinates"], os.SEEK_SET)
+ pos = self._read_field_from_file(f,
+ tp[ptype], "Coordinates")
+ psize[ptype] += selector.count_points(
+ pos[:,0], pos[:,1], pos[:,2])
+ del pos
+ f.close()
+ ind = {}
+ for field in fields:
+ mylog.debug("Allocating %s values for %s", psize[field[0]], field)
+ if field[1] in _vector_fields:
+ shape = (psize[field[0]], 3)
+ else:
+ shape = psize[field[0]]
+ rv[field] = np.empty(shape, dtype="float64")
+ ind[field] = 0
+ for chunk in chunks:
+ for subset in chunk.objs:
+ poff = subset.domain.field_offsets
+ tp = subset.domain.total_particles
+ f = open(subset.domain.domain_filename, "rb")
+ for ptype, field_list in sorted(ptf.items()):
+ f.seek(poff[ptype, "Coordinates"], os.SEEK_SET)
+ pos = self._read_field_from_file(f,
+ tp[ptype], "Coordinates")
+ mask = selector.select_points(
+ pos[:,0], pos[:,1], pos[:,2])
+ del pos
+ if mask is None: continue
+ for field in field_list:
+ f.seek(poff[ptype, field], os.SEEK_SET)
+ data = self._read_field_from_file(f, tp[ptype], field)
+ data = data[mask]
+ my_ind = ind[ptype, field]
+ mylog.debug("Filling from %s to %s with %s",
+ my_ind, my_ind+data.shape[0], field)
+ rv[ptype, field][my_ind:my_ind + data.shape[0],...] = data
+ ind[ptype, field] += data.shape[0]
+ f.close()
+ return rv
+
+ def _read_field_from_file(self, f, count, name):
+ if count == 0: return
+ if name == "ParticleIDs":
+ dt = "int32"
+ else:
+ dt = "float32"
+ if name in _vector_fields:
+ count *= 3
+ arr = np.fromfile(f, dtype=dt, count = count)
+ if name in _vector_fields:
+ arr = arr.reshape((count/3, 3), order="C")
+ return arr.astype("float64")
def _initialize_octree(self, domain, octree):
count = sum(domain.total_particles.values())
@@ -179,6 +249,26 @@
for i, v in enumerate(domain.header["Npart"]))
return npart
+ _header_offset = 256
+ _field_size = 4
+ def _calculate_field_offsets(self, field_list, pcount):
+ # field_list is (ftype, fname) but the blocks are ordered
+ # (fname, ftype) in the file.
+ pos = self._header_offset # 256 bytes for the header
+ fs = self._field_size
+ offsets = {}
+ for field in self._fields:
+ if not isinstance(field, types.StringTypes):
+ field = field[0]
+ for ptype in self._ptypes:
+ if (ptype, field) not in field_list: continue
+ offsets[(ptype, field)] = pos
+ if field in _vector_fields:
+ pos += 3 * pcount[ptype] * fs
+ else:
+ pos += pcount[ptype] * fs
+ return offsets
+
def _identify_fields(self, domain):
# We can just look at the particle counts.
field_list = []
https://bitbucket.org/yt_analysis/yt-3.0/commits/fdc02b6318b1/
changeset: fdc02b6318b1
branch: yt-3.0
user: MatthewTurk
date: 2013-01-12 00:06:50
summary: Fixes for ParticleOctreeContainer.
* Reverting domain_number => domain_id
* Fixing masking-by-domain of DomainSubsets.
* Adding max_level attribute for ParticleGeometry.
affected #: 3 files
diff -r c5ad47a08dde5fc1f2923f838cbb843e392fb1dc -r fdc02b6318b16e2cc3b32677c4c611c51479e99d yt/frontends/sph/data_structures.py
--- a/yt/frontends/sph/data_structures.py
+++ b/yt/frontends/sph/data_structures.py
@@ -51,11 +51,11 @@
FieldInfoContainer, NullFunc
class ParticleDomainFile(object):
- def __init__(self, pf, io, domain_filename, domain_number):
+ def __init__(self, pf, io, domain_filename, domain_id):
self.pf = pf
self.io = weakref.proxy(io)
self.domain_filename = domain_filename
- self.domain_number = domain_number
+ self.domain_id = domain_id
self.total_particles = self.io._count_particles(self)
def select(self, selector):
@@ -73,16 +73,19 @@
self.mask = mask
self.cell_count = count
self.oct_handler = domain.pf.h.oct_handler
+ level_counts = self.oct_handler.count_levels(
+ 99, self.domain.domain_id, mask)
+ level_counts[1:] = level_counts[:-1]
+ level_counts[0] = 0
+ self.level_counts = np.add.accumulate(level_counts)
def icoords(self, dobj):
return self.oct_handler.icoords(self.domain.domain_id, self.mask,
- self.cell_count,
- self.level_counts.copy())
+ self.cell_count)
def fcoords(self, dobj):
return self.oct_handler.fcoords(self.domain.domain_id, self.mask,
- self.cell_count,
- self.level_counts.copy())
+ self.cell_count)
def fwidth(self, dobj):
# Recall domain_dimensions is the number of cells, not octs
@@ -95,8 +98,7 @@
def ires(self, dobj):
return self.oct_handler.ires(self.domain.domain_id, self.mask,
- self.cell_count,
- self.level_counts.copy())
+ self.cell_count)
class ParticleGeometryHandler(OctreeGeometryHandler):
@@ -127,9 +129,11 @@
for dom in self.domains:
self.io._initialize_octree(dom, self.oct_handler)
self.oct_handler.finalize()
+ self.max_level = self.oct_handler.max_level
tot = self.oct_handler.linearly_count()
mylog.info("Identified %0.3e octs", tot)
+
def _detect_fields(self):
# TODO: Add additional fields
pfl = []
@@ -233,13 +237,13 @@
return False
class GadgetBinaryDomainFile(ParticleDomainFile):
- def __init__(self, pf, io, domain_filename, domain_number):
+ def __init__(self, pf, io, domain_filename, domain_id):
with open(domain_filename, "rb") as f:
self.header = read_record(f, pf._header_spec)
self._position_offset = f.tell()
super(GadgetBinaryDomainFile, self).__init__(pf, io,
- domain_filename, domain_number)
+ domain_filename, domain_id)
def _calculate_offsets(self, field_list):
self.field_offsets = self.io._calculate_field_offsets(
diff -r c5ad47a08dde5fc1f2923f838cbb843e392fb1dc -r fdc02b6318b16e2cc3b32677c4c611c51479e99d yt/frontends/sph/io.py
--- a/yt/frontends/sph/io.py
+++ b/yt/frontends/sph/io.py
@@ -94,7 +94,7 @@
for key in f.keys():
if not key.startswith("PartType"): continue
pos = f[key]["Coordinates"][:].astype("float64")
- octree.add(pos, domain.domain_number)
+ octree.add(pos, domain.domain_id)
f.close()
def _count_particles(self, domain):
@@ -242,7 +242,7 @@
pos[:,1] = pp['py']
pos[:,2] = pp['pz']
del pp
- octree.add(pos, domain.domain_number)
+ octree.add(pos, domain.domain_id)
def _count_particles(self, domain):
npart = dict((self._ptypes[i], v)
diff -r c5ad47a08dde5fc1f2923f838cbb843e392fb1dc -r fdc02b6318b16e2cc3b32677c4c611c51479e99d yt/geometry/oct_container.pyx
--- a/yt/geometry/oct_container.pyx
+++ b/yt/geometry/oct_container.pyx
@@ -558,6 +558,7 @@
cdef ParticleArrays *last_sd
cdef Oct** oct_list
cdef np.int64_t *dom_offsets
+ cdef public int max_level
def allocate_root(self):
cdef int i, j, k
@@ -599,7 +600,8 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- def icoords(self, np.ndarray[np.uint8_t, ndim=2, cast=True] mask,
+ def icoords(self, int domain_id,
+ np.ndarray[np.uint8_t, ndim=2, cast=True] mask,
np.int64_t cell_count):
cdef np.ndarray[np.int64_t, ndim=2] coords
coords = np.empty((cell_count, 3), dtype="int64")
@@ -607,6 +609,7 @@
ci = 0
for oi in range(self.nocts):
o = self.oct_list[oi]
+ if o.domain != domain_id: continue
for i in range(2):
for j in range(2):
for k in range(2):
@@ -621,7 +624,8 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- def ires(self, np.ndarray[np.uint8_t, ndim=2, cast=True] mask,
+ def ires(self, int domain_id,
+ np.ndarray[np.uint8_t, ndim=2, cast=True] mask,
np.int64_t cell_count):
cdef np.ndarray[np.int64_t, ndim=1] res
res = np.empty(cell_count, dtype="int64")
@@ -629,6 +633,7 @@
ci = 0
for oi in range(self.nocts):
o = self.oct_list[oi]
+ if o.domain != domain_id: continue
for i in range(8):
if mask[oi, i] == 1:
res[ci] = o.level
@@ -638,7 +643,8 @@
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
- def fcoords(self, np.ndarray[np.uint8_t, ndim=2, cast=True] mask,
+ def fcoords(self, int domain_id,
+ np.ndarray[np.uint8_t, ndim=2, cast=True] mask,
np.int64_t cell_count):
cdef np.ndarray[np.float64_t, ndim=2] coords
coords = np.empty((cell_count, 3), dtype="float64")
@@ -659,6 +665,7 @@
break
if proc == 0: continue
o = self.oct_list[oi]
+ if o.domain != domain_id: continue
for i in range(3):
# This gives the *grid* width for this level
dx[i] = base_dx[i] / (1 << o.level)
@@ -681,11 +688,13 @@
pass
def finalize(self):
+ cdef int max_level = 0
self.oct_list = <Oct**> malloc(sizeof(Oct*)*self.nocts)
cdef np.int64_t i = 0
cdef ParticleArrays *c = self.first_sd
while c != NULL:
self.oct_list[i] = c.oct
+ max_level = imax(max_level, c.oct.level)
if c.np >= 0:
for j in range(3):
free(c.pos[j])
@@ -693,6 +702,7 @@
c.pos = NULL
c = c.next
i += 1
+ self.max_level = max_level
qsort(self.oct_list, self.nocts, sizeof(Oct*), &compare_octs)
cdef int cur_dom = -1
# We always need at least 2, and if max_domain is 0, we need 3.
Repository URL: https://bitbucket.org/yt_analysis/yt-3.0/
--
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