[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