[yt-svn] commit/yt-3.0: 9 new changesets

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Wed Jun 12 10:35:52 PDT 2013


9 new commits in yt-3.0:

https://bitbucket.org/yt_analysis/yt-3.0/commits/e9af6a7b7fa2/
Changeset:   e9af6a7b7fa2
Branch:      yt-3.0
User:        MatthewTurk
Date:        2013-06-11 15:39:59
Summary:     Fixing 1.0/mass_sun_cgs in Enzo frontend.
Affected #:  1 file

diff -r 7a79dd9c7ebefddba48dfe0440a47ca3fddd6b98 -r e9af6a7b7fa22cd4eb45dff093a56776640035de yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -37,7 +37,8 @@
     ValidateGridType
 import yt.data_objects.universal_fields
 from yt.utilities.physical_constants import \
-    mh
+    mh, \
+    mass_sun_cgs
 from yt.funcs import *
 
 import yt.utilities.lib as amr_utils
@@ -81,7 +82,7 @@
     return data[sp] / _speciesMass[species]
 
 def _convertCellMassMsun(data):
-    return 5.027854e-34 # g^-1
+    return 1.0/mass_sun_cgs # g^-1
 def _ConvertNumberDensity(data):
     return 1.0/mh
 


https://bitbucket.org/yt_analysis/yt-3.0/commits/2deb3327354a/
Changeset:   2deb3327354a
Branch:      yt-3.0
User:        MatthewTurk
Date:        2013-06-11 15:42:22
Summary:     Adding Enzo-style species fields to RAMSES
Affected #:  1 file

diff -r e9af6a7b7fa22cd4eb45dff093a56776640035de -r 2deb3327354a362a9692147299fc989b9891303d yt/frontends/ramses/fields.py
--- a/yt/frontends/ramses/fields.py
+++ b/yt/frontends/ramses/fields.py
@@ -178,3 +178,71 @@
          units = r"\mathrm{g}/\mathrm{cm}^{3}",
          projected_units = r"\mathrm{g}/\mathrm{cm}^{-2}",
          projection_conversion = 'cm')
+
+# We'll add a bunch of species fields here.  In the not too distant future,
+# we'll be moving all of these to a unified field location, so they can be
+# shared between various frontends.
+
+# NOTE: No Electron here because I don't know how RAMSES handles them, and if
+# they are handled differently than Enzo does (where they are scaled to mh)
+
+_speciesList = ["HI", "HII",
+                "HeI", "HeII", "HeIII",
+                "H2I", "H2II", "HM",
+                "DI", "DII", "HDI"]
+_speciesMass = {"HI": 1.0, "HII": 1.0,
+                "HeI": 4.0, "HeII": 4.0, "HeIII": 4.0,
+                "H2I": 2.0, "H2II": 2.0, "HM": 1.0,
+                "DI": 2.0, "DII": 2.0, "HDI": 3.0}
+
+def _SpeciesComovingDensity(field, data):
+    sp = field.name.split("_")[0] + "_Density"
+    ef = (1.0 + data.pf.current_redshift)**3.0
+    return data[sp] / ef
+
+def _SpeciesFraction(field, data):
+    sp = field.name.split("_")[0] + "_Density"
+    return data[sp] / data["Density"]
+
+def _SpeciesMass(field, data):
+    sp = field.name.split("_")[0] + "_Density"
+    return data[sp] * data["CellVolume"]
+
+def _SpeciesNumberDensity(field, data):
+    species = field.name.split("_")[0]
+    sp = field.name.split("_")[0] + "_Density"
+    return data[sp] / _speciesMass[species]
+
+def _convertCellMassMsun(data):
+    return 1.0/mass_sun_cgs # g^-1
+def _ConvertNumberDensity(data):
+    return 1.0/mh
+
+for species in _speciesList:
+    add_ramses_field("%s_Density" % species,
+             function = NullFunc,
+             display_name = "%s\/Density" % species,
+             units = r"\rm{g}/\rm{cm}^3",
+             projected_units = r"\rm{g}/\rm{cm}^2")
+    add_field("%s_Fraction" % species,
+             function=_SpeciesFraction,
+             validators=ValidateDataField("%s_Density" % species),
+             display_name="%s\/Fraction" % species)
+    add_field("Comoving_%s_Density" % species,
+             function=_SpeciesComovingDensity,
+             validators=ValidateDataField("%s_Density" % species),
+             display_name="Comoving\/%s\/Density" % species)
+    add_field("%s_Mass" % species, units=r"\rm{g}", 
+              function=_SpeciesMass, 
+              validators=ValidateDataField("%s_Density" % species),
+              display_name="%s\/Mass" % species)
+    add_field("%s_MassMsun" % species, units=r"M_{\odot}", 
+              function=_SpeciesMass, 
+              convert_function=_convertCellMassMsun,
+              validators=ValidateDataField("%s_Density" % species),
+              display_name="%s\/Mass" % species)
+    if _speciesMass.has_key(species):
+        add_field("%s_NumberDensity" % species,
+                  function=_SpeciesNumberDensity,
+                  convert_function=_ConvertNumberDensity,
+                  validators=ValidateDataField("%s_Density" % species))


https://bitbucket.org/yt_analysis/yt-3.0/commits/3565d03875fe/
Changeset:   3565d03875fe
Branch:      yt-3.0
User:        MatthewTurk
Date:        2013-06-11 15:44:02
Summary:     We need a convert function for the density fields in RAMSES.
Affected #:  1 file

diff -r 2deb3327354a362a9692147299fc989b9891303d -r 3565d03875fed7118299145e6c2b9480a7e90ecb yt/frontends/ramses/fields.py
--- a/yt/frontends/ramses/fields.py
+++ b/yt/frontends/ramses/fields.py
@@ -222,6 +222,7 @@
     add_ramses_field("%s_Density" % species,
              function = NullFunc,
              display_name = "%s\/Density" % species,
+             convert_function = _convertDensity,
              units = r"\rm{g}/\rm{cm}^3",
              projected_units = r"\rm{g}/\rm{cm}^2")
     add_field("%s_Fraction" % species,


https://bitbucket.org/yt_analysis/yt-3.0/commits/4fb9244f4260/
Changeset:   4fb9244f4260
Branch:      yt-3.0
User:        MatthewTurk
Date:        2013-06-11 15:50:05
Summary:     Merging with OWLS and CIC deposition changes.
Affected #:  9 files

diff -r 3565d03875fed7118299145e6c2b9480a7e90ecb -r 4fb9244f42602d769afa0d55bbd903a49e592ee5 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -53,6 +53,7 @@
     MinimalProjectionData
 from yt.utilities.parallel_tools.parallel_analysis_interface import \
     parallel_objects, parallel_root_only, ParallelAnalysisInterface
+import yt.geometry.particle_deposit as particle_deposit
 
 from .field_info_container import\
     NeedsGridType,\
@@ -433,20 +434,31 @@
         fields_to_get = [f for f in fields if f not in self.field_data]
         fields_to_get = self._identify_dependencies(fields_to_get)
         if len(fields_to_get) == 0: return
-        fill, gen = self._split_fields(fields_to_get)
+        fill, gen, part = self._split_fields(fields_to_get)
+        if len(part) > 0: self._fill_particles(part)
         if len(fill) > 0: self._fill_fields(fill)
         if len(gen) > 0: self._generate_fields(gen)
 
     def _split_fields(self, fields_to_get):
         fill, gen = self.pf.h._split_fields(fields_to_get)
+        particles = []
         for field in gen:
             finfo = self.pf._get_field_info(*field)
             try:
                 finfo.check_available(self)
             except NeedsOriginalGrid:
                 fill.append(field)
+        for field in fill:
+            finfo = self.pf._get_field_info(*field)
+            if finfo.particle_type:
+                particles.append(field)
         gen = [f for f in gen if f not in fill]
-        return fill, gen
+        fill = [f for f in fill if f not in particles]
+        return fill, gen, particles
+
+    def _fill_particles(self, part):
+        for p in part:
+            self[p] = self._data_source[p]
 
     def _fill_fields(self, fields):
         output_fields = [np.zeros(self.ActiveDimensions, dtype="float64")
@@ -485,6 +497,67 @@
             raise KeyError(field)
         return rv
 
+    @property
+    def LeftEdge(self):
+        return self.left_edge
+
+    def deposit(self, positions, fields = None, method = None):
+        cls = getattr(particle_deposit, "deposit_%s" % method, None)
+        if cls is None:
+            raise YTParticleDepositionNotImplemented(method)
+        op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs
+        op.initialize()
+        op.process_grid(self, positions, fields)
+        vals = op.finalize()
+        return vals.reshape(self.ActiveDimensions, order="F")
+
+class YTArbitraryGridBase(YTCoveringGridBase):
+    """A 3D region with arbitrary bounds and dimensions.
+
+    In contrast to the Covering Grid, this object accepts a left edge, a right
+    edge, and dimensions.  This allows it to be used for creating 3D particle
+    deposition fields that are independent of the underlying mesh, whether that
+    is yt-generated or from the simulation data.  For example, arbitrary boxes
+    around particles can be drawn and particle deposition fields can be
+    created.  This object will refuse to generate any fluid fields.
+    
+    Parameters
+    ----------
+    left_edge : array_like
+        The left edge of the region to be extracted
+    rigth_edge : array_like
+        The left edge of the region to be extracted
+    dims : array_like
+        Number of cells along each axis of resulting grid.
+
+    Examples
+    --------
+    >>> obj = pf.h.arbitrary_grid([0.0, 0.0, 0.0], [0.99, 0.99, 0.99],
+    ...                          dims=[128, 128, 128])
+    """
+    _spatial = True
+    _type_name = "arbitrary_grid"
+    _con_args = ('left_edge', 'right_edge', 'ActiveDimensions')
+    _container_fields = ("dx", "dy", "dz", "x", "y", "z")
+    def __init__(self, left_edge, right_edge, dims,
+                 pf = None, field_parameters = None):
+        if field_parameters is None:
+            center = None
+        else:
+            center = field_parameters.get("center", None)
+        YTSelectionContainer3D.__init__(self, center, pf, field_parameters)
+        self.left_edge = np.array(left_edge)
+        self.right_edge = np.array(right_edge)
+        self.ActiveDimensions = np.array(dims, dtype='int32')
+        if self.ActiveDimensions.size == 1:
+            self.ActiveDimensions = np.array([dims, dims, dims], dtype="int32")
+        self.dds = (self.right_edge - self.left_edge)/self.ActiveDimensions
+        self.level = 99
+        self._setup_data_source()
+
+    def _fill_fields(self, fields):
+        raise NotImplementedError
+
 class LevelState(object):
     current_dx = None
     current_dims = None

diff -r 3565d03875fed7118299145e6c2b9480a7e90ecb -r 4fb9244f42602d769afa0d55bbd903a49e592ee5 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -477,8 +477,23 @@
         if fields is None: return
         fields = self._determine_fields(fields)
         # Now we collect all our fields
-        fields_to_get = [f for f in fields if f not in self.field_data]
-        if len(fields_to_get) == 0:
+        # Here is where we need to perform a validation step, so that if we
+        # have a field requested that we actually *can't* yet get, we put it
+        # off until the end.  This prevents double-reading fields that will
+        # need to be used in spatial fields later on.
+        fields_to_get = []
+        # This will be pre-populated with spatial fields
+        fields_to_generate = [] 
+        for field in self._determine_fields(fields):
+            if field in self.field_data: continue
+            finfo = self.pf._get_field_info(*field)
+            try:
+                finfo.check_available(self)
+            except NeedsGridType:
+                fields_to_generate.append(field)
+                continue
+            fields_to_get.append(field)
+        if len(fields_to_get) == 0 and fields_to_generate == 0:
             return
         elif self._locked == True:
             raise GenerationInProgress(fields)
@@ -502,12 +517,18 @@
         read_particles, gen_particles = self.hierarchy._read_particle_fields(
                                         particles, self, self._current_chunk)
         self.field_data.update(read_particles)
-        fields_to_generate = gen_fluids + gen_particles
+        fields_to_generate += gen_fluids + gen_particles
         self._generate_fields(fields_to_generate)
 
     def _generate_fields(self, fields_to_generate):
         index = 0
         with self._field_lock():
+            # At this point, we assume that any fields that are necessary to
+            # *generate* a field are in fact already available to us.  Note
+            # that we do not make any assumption about whether or not the
+            # fields have a spatial requirement.  This will be checked inside
+            # _generate_field, at which point additional dependencies may
+            # actually be noted.
             while any(f not in self.field_data for f in fields_to_generate):
                 field = fields_to_generate[index % len(fields_to_generate)]
                 index += 1
@@ -519,11 +540,6 @@
                         if f not in fields_to_generate:
                             fields_to_generate.append(f)
 
-    def deposit(self, positions, fields, op):
-        assert(self._current_chunk.chunk_type == "spatial")
-        fields = ensure_list(fields)
-        self.hierarchy._deposit_particle_fields(self, positions, fields, op)
-
     @contextmanager
     def _field_lock(self):
         self._locked = True

diff -r 3565d03875fed7118299145e6c2b9480a7e90ecb -r 4fb9244f42602d769afa0d55bbd903a49e592ee5 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -894,7 +894,7 @@
     Child cells are not returned.
     """
     _type_name = "data_collection"
-    _con_args = ("obj_list",)
+    _con_args = ("_obj_list",)
     def __init__(self, center, obj_list, pf = None, field_parameters = None):
         YTSelectionContainer3D.__init__(self, center, pf, field_parameters)
         self._obj_ids = np.array([o.id - o._id_offset for o in obj_list],

diff -r 3565d03875fed7118299145e6c2b9480a7e90ecb -r 4fb9244f42602d769afa0d55bbd903a49e592ee5 yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -697,3 +697,15 @@
          units = r"\mathrm{g}/\mathrm{cm}^{3}",
          projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
          projection_conversion = 'cm')
+
+def particle_cic(field, data):
+    pos = np.column_stack([data["particle_position_%s" % ax] for ax in 'xyz'])
+    d = data.deposit(pos, [data["ParticleMass"]], method="cic")
+    d /= data["CellVolume"]
+    return d
+add_field(("deposit", "all_cic"), function=particle_cic,
+          validators=[ValidateSpatial()],
+          display_name = "\\mathrm{All CIC Density}",
+          units = r"\mathrm{g}/\mathrm{cm}^{3}",
+          projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
+          projection_conversion = 'cm')

diff -r 3565d03875fed7118299145e6c2b9480a7e90ecb -r 4fb9244f42602d769afa0d55bbd903a49e592ee5 yt/frontends/sph/data_structures.py
--- a/yt/frontends/sph/data_structures.py
+++ b/yt/frontends/sph/data_structures.py
@@ -163,71 +163,6 @@
         for subset in oobjs:
             yield YTDataChunk(dobj, "io", [subset], subset.cell_count)
 
-class OWLSStaticOutput(StaticOutput):
-    _hierarchy_class = ParticleGeometryHandler
-    _domain_class = ParticleDomainFile
-    _fieldinfo_fallback = OWLSFieldInfo
-    _fieldinfo_known = KnownOWLSFields
-
-    def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
-        self._root_dimensions = root_dimensions
-        # Set up the template for domain files
-        self.storage_filename = None
-        super(OWLSStaticOutput, 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 = {}
-        DW = self.domain_right_edge - self.domain_left_edge
-        self.units["unitary"] = 1.0 / DW.max()
-
-    def _parse_parameter_file(self):
-        handle = h5py.File(self.parameter_filename)
-        hvals = {}
-        hvals.update(handle["/Header"].attrs)
-
-        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
-        self.current_time = hvals["Time_GYR"] * 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.periodicity = (True, True, True)
-        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]
-        self.domain_template = "%s.%%(num)i.%s" % (prefix, suffix)
-        self.domain_count = hvals["NumFilesPerSnapshot"]
-
-        handle.close()
-
-    @classmethod
-    def _is_valid(self, *args, **kwargs):
-        try:
-            fileh = h5py.File(args[0],'r')
-            if "Constants" in fileh["/"].keys() and \
-               "Header" in fileh["/"].keys():
-                fileh.close()
-                return True
-            fileh.close()
-        except:
-            pass
-        return False
-
 class GadgetBinaryDomainFile(ParticleDomainFile):
     def __init__(self, pf, io, domain_filename, domain_id):
         with open(domain_filename, "rb") as f:
@@ -404,6 +339,74 @@
         # We do not allow load() of these files.
         return False
 
+class OWLSStaticOutput(GadgetStaticOutput):
+    _hierarchy_class = ParticleGeometryHandler
+    _domain_class = ParticleDomainFile
+    _fieldinfo_fallback = OWLSFieldInfo # For now we have separate from Gadget
+    _fieldinfo_known = KnownOWLSFields
+    _header_spec = None # Override so that there's no confusion
+
+    def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
+        self._root_dimensions = root_dimensions
+        # Set up the template for domain files
+        self.storage_filename = None
+        super(OWLSStaticOutput, self).__init__(filename, data_style,
+                                               root_dimensions,
+                                               unit_base = None)
+
+    def __repr__(self):
+        return os.path.basename(self.parameter_filename).split(".")[0]
+
+    def _parse_parameter_file(self):
+        handle = h5py.File(self.parameter_filename)
+        hvals = {}
+        hvals.update((str(k), v) for k, v in handle["/Header"].attrs.items())
+
+        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
+        self.current_time = hvals["Time_GYR"] * 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.periodicity = (True, True, True)
+        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]
+        self.domain_template = "%s.%%(num)i.%s" % (prefix, suffix)
+        self.domain_count = hvals["NumFilesPerSnapshot"]
+
+        # To avoid having to open files twice
+        self._unit_base = {}
+        self._unit_base.update((str(k), v) for k, v in handle["/Units"].attrs.items())
+        # Comoving cm is given in the Units
+        self._unit_base['cmcm'] = 1.0 / self._unit_base["UnitLength_in_cm"]
+
+        handle.close()
+
+    @classmethod
+    def _is_valid(self, *args, **kwargs):
+        try:
+            fileh = h5py.File(args[0],'r')
+            if "Constants" in fileh["/"].keys() and \
+               "Header" in fileh["/"].keys():
+                fileh.close()
+                return True
+            fileh.close()
+        except:
+            pass
+        return False
+
 class TipsyDomainFile(ParticleDomainFile):
 
     def _calculate_offsets(self, field_list):

diff -r 3565d03875fed7118299145e6c2b9480a7e90ecb -r 4fb9244f42602d769afa0d55bbd903a49e592ee5 yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -96,6 +96,20 @@
              projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
              projection_conversion = 'cm')
 
+    def particle_cic(field, data):
+        pos = data[ptype, coord_name]
+        d = data.deposit(pos, [data[ptype, mass_name]], method = "cic")
+        d /= data["CellVolume"]
+        return d
+
+    registry.add_field(("deposit", "%s_cic" % ptype),
+             function = particle_cic,
+             validators = [ValidateSpatial()],
+             display_name = "\\mathrm{%s CIC Density}" % ptype,
+             units = r"\mathrm{g}/\mathrm{cm}^{3}",
+             projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
+             projection_conversion = 'cm')
+
     # Now some translation functions.
 
     registry.add_field((ptype, "ParticleMass"),
@@ -247,3 +261,60 @@
         GadgetFieldInfo.add_field(("all", oname + ax), function=func,
                 particle_type = True)
 
+# OWLS
+# ====
+
+# I am optimistic that some day we will be able to get rid of much of this, and
+# make OWLS a subclass of Gadget fields.
+
+_owls_ptypes = ("PartType0", "PartType1", "PartType2", "PartType3",
+                "PartType4")
+
+for fname in ["Coordinates", "Velocities", "ParticleIDs",
+              # Note: Mass, not Masses
+              "Mass"]:
+    func = _field_concat(fname)
+    OWLSFieldInfo.add_field(("all", fname), function=func,
+            particle_type = True)
+
+def _owls_particle_fields(ptype):
+    def _Mass(field, data):
+        pind = _owls_ptypes.index(ptype)
+        if data.pf["MassTable"][pind] == 0.0:
+            raise RuntimeError
+        mass = np.ones(data[ptype, "ParticleIDs"].shape[0], dtype="float64")
+        # Note that this is an alias, which is why we need to apply conversion
+        # here.  Otherwise we'd have an asymmetry.
+        mass *= data.pf["MassTable"][pind] 
+        return mass
+    OWLSFieldInfo.add_field((ptype, "Mass"), function=_Mass,
+                            convert_function = _get_conv("mass"),
+                            particle_type = True)
+
+for ptype in _owls_ptypes:
+    # Note that this adds a "Known" Mass field and a "Derived" Mass field.
+    # This way the "Known" will get used, and if it's not there, it will use
+    # the derived.
+    KnownOWLSFields.add_field((ptype, "Mass"), function=NullFunc,
+        particle_type = True,
+        convert_function=_get_conv("mass"),
+        units = r"\mathrm{g}")
+    _owls_particle_fields(ptype)
+    KnownOWLSFields.add_field((ptype, "Velocities"), function=NullFunc,
+        particle_type = True,
+        convert_function=_get_conv("velocity"),
+        units = r"\mathrm{cm}/\mathrm{s}")
+    _particle_functions(ptype, "Coordinates", "Mass", OWLSFieldInfo)
+    KnownOWLSFields.add_field((ptype, "Coordinates"), function=NullFunc,
+        particle_type = True)
+_particle_functions("all", "Coordinates", "Mass", OWLSFieldInfo)
+
+# Now we have to manually apply the splits for "all", since we don't want to
+# use the splits defined above.
+
+for iname, oname in [("Coordinates", "particle_position_"),
+                     ("Velocities", "particle_velocity_")]:
+    for axi, ax in enumerate("xyz"):
+        func = _field_concat_slice(iname, axi)
+        OWLSFieldInfo.add_field(("all", oname + ax), function=func,
+                particle_type = True)

diff -r 3565d03875fed7118299145e6c2b9480a7e90ecb -r 4fb9244f42602d769afa0d55bbd903a49e592ee5 yt/geometry/oct_geometry_handler.py
--- a/yt/geometry/oct_geometry_handler.py
+++ b/yt/geometry/oct_geometry_handler.py
@@ -78,7 +78,7 @@
             source.quantities["MaxLocation"](field)
         mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f", 
               max_val, mx, my, mz)
-        self.pf.parameters["Max%sValue" % (field)] = max_val
-        self.pf.parameters["Max%sPos" % (field)] = "%s" % ((mx,my,mz),)
+        self.pf.parameters["Max%sValue" % (field,)] = max_val
+        self.pf.parameters["Max%sPos" % (field,)] = "%s" % ((mx,my,mz),)
         return max_val, np.array((mx,my,mz), dtype='float64')
 

diff -r 3565d03875fed7118299145e6c2b9480a7e90ecb -r 4fb9244f42602d769afa0d55bbd903a49e592ee5 yt/geometry/particle_deposit.pxd
--- a/yt/geometry/particle_deposit.pxd
+++ b/yt/geometry/particle_deposit.pxd
@@ -29,6 +29,7 @@
 import numpy as np
 from libc.stdlib cimport malloc, free
 cimport cython
+from libc.math cimport sqrt
 
 from fp_utils cimport *
 from oct_container cimport Oct, OctAllocationContainer, OctreeContainer
@@ -39,6 +40,21 @@
 cdef inline int gind(int i, int j, int k, int dims[3]):
     return ((k*dims[1])+j)*dims[0]+i
 
+
+####################################################
+# Standard SPH kernel for use with the Grid method #
+####################################################
+
+cdef inline np.float64_t sph_kernel(np.float64_t x) nogil:
+    cdef np.float64_t kernel
+    if x <= 0.5:
+        kernel = 1.-6.*x*x*(1.-x)
+    elif x>0.5 and x<=1.0:
+        kernel = 2.*(1.-x)*(1.-x)*(1.-x)
+    else:
+        kernel = 0.
+    return kernel
+
 cdef class ParticleDepositOperation:
     # We assume each will allocate and define their own temporary storage
     cdef np.int64_t nvals

diff -r 3565d03875fed7118299145e6c2b9480a7e90ecb -r 4fb9244f42602d769afa0d55bbd903a49e592ee5 yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -29,6 +29,7 @@
 import numpy as np
 from libc.stdlib cimport malloc, free
 cimport cython
+from libc.math cimport sqrt
 
 from fp_utils cimport *
 from oct_container cimport Oct, OctAllocationContainer, \
@@ -145,6 +146,68 @@
 
 deposit_count = CountParticles
 
+cdef class SimpleSmooth(ParticleDepositOperation):
+    # Note that this does nothing at the edges.  So it will give a poor
+    # estimate there, and since Octrees are mostly edges, this will be a very
+    # poor SPH kernel.
+    cdef np.float64_t *data
+    cdef public object odata
+    cdef np.float64_t *temp
+    cdef public object otemp
+
+    def initialize(self):
+        self.odata = np.zeros(self.nvals, dtype="float64")
+        cdef np.ndarray arr = self.odata
+        self.data = <np.float64_t*> arr.data
+        self.otemp = np.zeros(self.nvals, dtype="float64")
+        arr = self.otemp
+        self.temp = <np.float64_t*> arr.data
+
+    @cython.cdivision(True)
+    cdef void process(self, int dim[3],
+                      np.float64_t left_edge[3],
+                      np.float64_t dds[3],
+                      np.int64_t offset,
+                      np.float64_t ppos[3],
+                      np.float64_t *fields
+                      ):
+        cdef int ii[3], half_len, ib0[3], ib1[3]
+        cdef int i, j, k
+        cdef np.float64_t idist[3], kernel_sum, dist
+        # Smoothing length is fields[0]
+        kernel_sum = 0.0
+        for i in range(3):
+            ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
+            half_len = <int>(fields[0]/dds[i]) + 1
+            ib0[i] = ii[i] - half_len
+            ib1[i] = ii[i] + half_len
+            if ib0[i] >= dim[i] or ib1[i] <0:
+                return
+            ib0[i] = iclip(ib0[i], 0, dim[i] - 1)
+            ib1[i] = iclip(ib1[i], 0, dim[i] - 1)
+        for i from ib0[0] <= i <= ib1[0]:
+            idist[0] = (ii[0] - i) * (ii[0] - i) * dds[0]
+            for j from ib0[1] <= j <= ib1[1]:
+                idist[1] = (ii[1] - j) * (ii[1] - j) * dds[1] 
+                for k from ib0[2] <= k <= ib1[2]:
+                    idist[2] = (ii[2] - k) * (ii[2] - k) * dds[2]
+                    dist = idist[0] + idist[1] + idist[2]
+                    # Calculate distance in multiples of the smoothing length
+                    dist = sqrt(dist) / fields[0]
+                    self.temp[gind(i,j,k,dim) + offset] = sph_kernel(dist)
+                    kernel_sum += self.temp[gind(i,j,k,dim) + offset]
+        # Having found the kernel, deposit accordingly into gdata
+        for i from ib0[0] <= i <= ib1[0]:
+            for j from ib0[1] <= j <= ib1[1]:
+                for k from ib0[2] <= k <= ib1[2]:
+                    dist = self.temp[gind(i,j,k,dim) + offset] / kernel_sum
+                    self.data[gind(i,j,k,dim) + offset] += fields[1] * dist
+        
+    def finalize(self):
+        return self.odata
+
+deposit_simple_smooth = SimpleSmooth
+
 cdef class SumParticleField(ParticleDepositOperation):
     cdef np.float64_t *sum
     cdef public object osum
@@ -232,3 +295,44 @@
 
 deposit_std = StdParticleField
 
+cdef class CICDeposit(ParticleDepositOperation):
+    cdef np.float64_t *field
+    cdef public object ofield
+    def initialize(self):
+        self.ofield = np.zeros(self.nvals, dtype="float64")
+        cdef np.ndarray arr = self.ofield
+        self.field = <np.float64_t *> arr.data
+
+    cdef void process(self, int dim[3],
+                      np.float64_t left_edge[3],
+                      np.float64_t dds[3],
+                      np.int64_t offset, # offset into IO field
+                      np.float64_t ppos[3], # this particle's position
+                      np.float64_t *fields # any other fields we need
+                      ):
+        
+        cdef int i, j, k, ind[3], ii
+        cdef np.float64_t rpos[3], rdds[3][2]
+        cdef np.float64_t fact, edge0, edge1, edge2
+        cdef np.float64_t le0, le1, le2
+        cdef np.float64_t dx, dy, dz, dx2, dy2, dz2
+
+        # Compute the position of the central cell
+        for i in range(3):
+            rpos[i] = (ppos[i]-left_edge[i])/dds[i]
+            rpos[i] = fclip(rpos[i], 0.5001, dim[i]-0.5001)
+            ind[i] = <int> (rpos[i] + 0.5)
+            # Note these are 1, then 0
+            rdds[i][1] = (<np.float64_t> ind[i]) + 0.5 - rpos[i]
+            rdds[i][0] = 1.0 - rdds[i][1]
+
+        for i in range(2):
+            for j in range(2):
+                for k in range(2):
+                    ii = gind(ind[0] - i, ind[1] - j, ind[2] - k, dim) + offset
+                    self.field[ii] += fields[0]*rdds[0][i]*rdds[1][j]*rdds[2][k]
+
+    def finalize(self):
+        return self.ofield
+
+deposit_cic = CICDeposit


https://bitbucket.org/yt_analysis/yt-3.0/commits/a4bbc055e03e/
Changeset:   a4bbc055e03e
Branch:      yt-3.0
User:        MatthewTurk
Date:        2013-06-11 17:07:39
Summary:     Adding a particle fields registry.
Affected #:  1 file

diff -r 4fb9244f42602d769afa0d55bbd903a49e592ee5 -r a4bbc055e03edf2c8eccdf0988b7c463a27656e9 yt/data_objects/particle_fields.py
--- /dev/null
+++ b/yt/data_objects/particle_fields.py
@@ -0,0 +1,147 @@
+"""
+These are common particle deposition fields.
+
+Author: Matthew Turk <matthewturk at gmail.com>
+Affiliation: Columbia University
+Homepage: http://yt-project.org/
+License:
+  Copyright (C) 2013 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/>.
+"""
+
+import numpy as np
+
+from yt.funcs import *
+from yt.data_objects.field_info_container import \
+    FieldInfoContainer, \
+    FieldInfo, \
+    ValidateParameter, \
+    ValidateDataField, \
+    ValidateProperty, \
+    ValidateSpatial, \
+    ValidateGridType, \
+    NullFunc, \
+    TranslationFunc
+from yt.utilities.physical_constants import \
+    mass_hydrogen_cgs, \
+    mass_sun_cgs, \
+    mh
+
+def particle_deposition_functions(ptype, coord_name, mass_name, registry):
+    def particle_count(field, data):
+        pos = data[ptype, coord_name]
+        d = data.deposit(pos, method = "count")
+        return d
+    registry.add_field(("deposit", "%s_count" % ptype),
+             function = particle_count,
+             validators = [ValidateSpatial()],
+             display_name = "\\mathrm{%s Count}" % ptype,
+             projection_conversion = '1')
+
+    def particle_mass(field, data):
+        pos = data[ptype, coord_name]
+        d = data.deposit(pos, [data[ptype, mass_name]], method = "sum")
+        return d
+
+    registry.add_field(("deposit", "%s_mass" % ptype),
+             function = particle_mass,
+             validators = [ValidateSpatial()],
+             display_name = "\\mathrm{%s Mass}" % ptype,
+             units = r"\mathrm{g}",
+             projected_units = r"\mathrm{g}\/\mathrm{cm}",
+             projection_conversion = 'cm')
+
+    def particle_density(field, data):
+        pos = data[ptype, coord_name]
+        d = data.deposit(pos, [data[ptype, mass_name]], method = "sum")
+        d /= data["CellVolume"]
+        return d
+
+    registry.add_field(("deposit", "%s_density" % ptype),
+             function = particle_density,
+             validators = [ValidateSpatial()],
+             display_name = "\\mathrm{%s Density}" % ptype,
+             units = r"\mathrm{g}/\mathrm{cm}^{3}",
+             projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
+             projection_conversion = 'cm')
+
+    def particle_cic(field, data):
+        pos = data[ptype, coord_name]
+        d = data.deposit(pos, [data[ptype, mass_name]], method = "cic")
+        d /= data["CellVolume"]
+        return d
+
+    registry.add_field(("deposit", "%s_cic" % ptype),
+             function = particle_cic,
+             validators = [ValidateSpatial()],
+             display_name = "\\mathrm{%s CIC Density}" % ptype,
+             units = r"\mathrm{g}/\mathrm{cm}^{3}",
+             projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
+             projection_conversion = 'cm')
+
+    # Now some translation functions.
+
+    registry.add_field((ptype, "ParticleMass"),
+            function = TranslationFunc((ptype, mass_name)),
+            particle_type = True,
+            units = r"\mathrm{g}")
+
+    def _ParticleMassMsun(field, data):
+        return data[ptype, mass_name].copy()
+    def _conv_Msun(data):
+        return 1.0/mass_sun_cgs
+
+    registry.add_field((ptype, "ParticleMassMsun"),
+            function = _ParticleMassMsun,
+            convert_function = _conv_Msun,
+            particle_type = True,
+            units = r"\mathrm{M}_\odot")
+
+def particle_scalar_functions(ptype, coord_name, vel_name, registry):
+
+    # Now we have to set up the various velocity and coordinate things.  In the
+    # future, we'll actually invert this and use the 3-component items
+    # elsewhere, and stop using these.
+    
+    # Note that we pass in _ptype here so that it's defined inside the closure.
+    def _get_coord_funcs(axi, _ptype):
+        def _particle_velocity(field, data):
+            return data[_ptype, vel_name][:,axi]
+        def _particle_position(field, data):
+            return data[_ptype, coord_name][:,axi]
+        return _particle_velocity, _particle_position
+    for axi, ax in enumerate("xyz"):
+        v, p = _get_coord_funcs(axi, ptype)
+        registry.add_field((ptype, "particle_velocity_%s" % ax),
+            particle_type = True, function = v)
+        registry.add_field((ptype, "particle_position_%s" % ax),
+            particle_type = True, function = p)
+
+def particle_vector_functions(ptype, coord_names, vel_names, registry):
+
+    # This will column_stack a set of scalars to create vector fields.
+
+    def _get_vec_func(_ptype, names):
+        def particle_vectors(field, data):
+            return np.column_stack([data[_ptype, name] for name in names])
+        return particle_vectors
+    registry.add_field((ptype, "Coordinates"),
+                       function=_get_vec_func(ptype, coord_names),
+                       particle_type=True)
+    registry.add_field((ptype, "Velocities"),
+                       function=_get_vec_func(ptype, vel_names),
+                       particle_type=True)


https://bitbucket.org/yt_analysis/yt-3.0/commits/27f5b715b5cf/
Changeset:   27f5b715b5cf
Branch:      yt-3.0
User:        MatthewTurk
Date:        2013-06-11 17:09:29
Summary:     This converts Enzo to use the particle field registry.

Note that I also change a couple of the ways fields are detected, so that
fields which are found in the HDF5 files but are in fact particle fields will
(I believe) be correctly identified as such.  We now also explicitly add
("all", "whatever") fields.
Affected #:  2 files

diff -r a4bbc055e03edf2c8eccdf0988b7c463a27656e9 -r 27f5b715b5cf07575fe4c42bef402e9f9039f798 yt/frontends/enzo/data_structures.py
--- a/yt/frontends/enzo/data_structures.py
+++ b/yt/frontends/enzo/data_structures.py
@@ -491,7 +491,13 @@
         else:
             field_list = None
         field_list = self.comm.mpi_bcast(field_list)
-        self.field_list = list(field_list)
+        self.field_list = []
+        # Now we will, avoiding the problem of particle types not having names.
+        for field in field_list:
+            if ("all", field) in KnownEnzoFields:
+                self.field_list.append(("all", field))
+            else:
+                self.field_list.append(field)
 
     def _generate_random_grids(self):
         if self.num_grids > 40:

diff -r a4bbc055e03edf2c8eccdf0988b7c463a27656e9 -r 27f5b715b5cf07575fe4c42bef402e9f9039f798 yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -36,6 +36,9 @@
     ValidateSpatial, \
     ValidateGridType
 import yt.data_objects.universal_fields
+from yt.data_objects.particle_fields import \
+    particle_deposition_functions, \
+    particle_vector_functions
 from yt.utilities.physical_constants import \
     mh, \
     mass_sun_cgs
@@ -528,11 +531,11 @@
 
 for pf in ["type", "mass"] + \
           ["position_%s" % ax for ax in 'xyz']:
-    add_enzo_field("particle_%s" % pf, NullFunc, particle_type=True)
+    add_enzo_field(("all", "particle_%s" % pf), NullFunc, particle_type=True)
     
 def _convRetainInt(data):
     return 1
-add_enzo_field("particle_index", function=NullFunc,
+add_enzo_field(("all", "particle_index"), function=NullFunc,
           particle_type=True, convert_function=_convRetainInt)
 
 def _get_vel_convert(ax):
@@ -542,11 +545,11 @@
 for ax in 'xyz':
     pf = "particle_velocity_%s" % ax
     cfunc = _get_vel_convert(ax)
-    add_enzo_field(pf, function=NullFunc, convert_function=cfunc,
+    add_enzo_field(("all", pf), function=NullFunc, convert_function=cfunc,
               particle_type=True)
 
 for pf in ["creation_time", "dynamical_time", "metallicity_fraction"]:
-    add_enzo_field(pf, function=NullFunc,
+    add_enzo_field(("all", pf), function=NullFunc,
               validators = [ValidateDataField(pf)],
               particle_type=True)
 add_field("particle_mass", function=NullFunc, particle_type=True)
@@ -661,51 +664,8 @@
                function=TranslationFunc(("CenOstriker","position_%s" % ax)),
                particle_type = True)
 
-def particle_count(field, data):
-    pos = np.column_stack([data["particle_position_%s" % ax] for ax in 'xyz'])
-    d = data.deposit(pos, method = "count")
-    return d
-EnzoFieldInfo.add_field(("deposit", "%s_count" % "all"),
-         function = particle_count,
-         validators = [ValidateSpatial()],
-         display_name = "\\mathrm{%s Count}" % "all",
-         projection_conversion = '1')
-
-def particle_mass(field, data):
-    pos = np.column_stack([data["particle_position_%s" % ax] for ax in 'xyz'])
-    d = data.deposit(pos, [data["ParticleMass"]], method = "sum")
-    return d
-
-EnzoFieldInfo.add_field(("deposit", "%s_mass" % "all"),
-         function = particle_mass,
-         validators = [ValidateSpatial()],
-         display_name = "\\mathrm{%s Mass}" % "all",
-         units = r"\mathrm{g}",
-         projected_units = r"\mathrm{g}\/\mathrm{cm}",
-         projection_conversion = 'cm')
-
-def particle_density(field, data):
-    pos = np.column_stack([data["particle_position_%s" % ax] for ax in 'xyz'])
-    d = data.deposit(pos, [data["ParticleMass"]], method = "sum")
-    d /= data["CellVolume"]
-    return d
-
-EnzoFieldInfo.add_field(("deposit", "%s_density" % "all"),
-         function = particle_density,
-         validators = [ValidateSpatial()],
-         display_name = "\\mathrm{%s Density}" % "all",
-         units = r"\mathrm{g}/\mathrm{cm}^{3}",
-         projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
-         projection_conversion = 'cm')
-
-def particle_cic(field, data):
-    pos = np.column_stack([data["particle_position_%s" % ax] for ax in 'xyz'])
-    d = data.deposit(pos, [data["ParticleMass"]], method="cic")
-    d /= data["CellVolume"]
-    return d
-add_field(("deposit", "all_cic"), function=particle_cic,
-          validators=[ValidateSpatial()],
-          display_name = "\\mathrm{All CIC Density}",
-          units = r"\mathrm{g}/\mathrm{cm}^{3}",
-          projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
-          projection_conversion = 'cm')
+particle_vector_functions("all", ["particle_position_%s" % ax for ax in 'xyz'],
+                                 ["particle_velocity_%s" % ax for ax in 'xyz'],
+                          EnzoFieldInfo)
+particle_deposition_functions("all", "Coordinates", "ParticleMass",
+                               EnzoFieldInfo)


https://bitbucket.org/yt_analysis/yt-3.0/commits/358f9ca48941/
Changeset:   358f9ca48941
Branch:      yt-3.0
User:        MatthewTurk
Date:        2013-06-11 17:10:14
Summary:     Convert RAMSES to use ("all", "field_name") for particle fields and the
particle field registry.
Affected #:  3 files

diff -r 27f5b715b5cf07575fe4c42bef402e9f9039f798 -r 358f9ca489417c9a0800ec302e6729be9ee34c76 yt/frontends/ramses/data_structures.py
--- a/yt/frontends/ramses/data_structures.py
+++ b/yt/frontends/ramses/data_structures.py
@@ -152,8 +152,8 @@
         _pfields = {}
         for field, vtype in particle_fields:
             if f.tell() >= flen: break
-            field_offsets[field] = f.tell()
-            _pfields[field] = vtype
+            field_offsets["all", field] = f.tell()
+            _pfields["all", field] = vtype
             fpu.skip(f, 1)
         self.particle_field_offsets = field_offsets
         self.particle_field_types = _pfields

diff -r 27f5b715b5cf07575fe4c42bef402e9f9039f798 -r 358f9ca489417c9a0800ec302e6729be9ee34c76 yt/frontends/ramses/fields.py
--- a/yt/frontends/ramses/fields.py
+++ b/yt/frontends/ramses/fields.py
@@ -34,10 +34,14 @@
     ValidateSpatial, \
     ValidateGridType
 import yt.data_objects.universal_fields
+from yt.data_objects.particle_fields import \
+    particle_deposition_functions, \
+    particle_vector_functions
 from yt.utilities.physical_constants import \
     boltzmann_constant_cgs, \
     mass_hydrogen_cgs, \
-    mass_sun_cgs
+    mass_sun_cgs, \
+    mh
 import numpy as np
 
 RAMSESFieldInfo = FieldInfoContainer.create_with_fallback(FieldInfo, "RFI")
@@ -106,29 +110,19 @@
 ]
 
 for f in known_ramses_particle_fields:
-    if f not in KnownRAMSESFields:
-        add_ramses_field(f, function=NullFunc, take_log=True,
-                  validators = [ValidateDataField(f)],
-                  particle_type = True)
+    add_ramses_field(("all", f), function=NullFunc, take_log=True,
+              particle_type = True)
 
 for ax in 'xyz':
-    KnownRAMSESFields["particle_velocity_%s" % ax]._convert_function = \
+    KnownRAMSESFields["all", "particle_velocity_%s" % ax]._convert_function = \
         _convertVelocity
 
 def _convertParticleMass(data):
     return data.convert("mass")
 
-KnownRAMSESFields["particle_mass"]._convert_function = \
+KnownRAMSESFields["all", "particle_mass"]._convert_function = \
         _convertParticleMass
-KnownRAMSESFields["particle_mass"]._units = r"\mathrm{g}"
-
-def _convertParticleMassMsun(data):
-    return 1.0/mass_sun_cgs
-add_field("ParticleMass", function=TranslationFunc("particle_mass"), 
-          particle_type=True)
-add_field("ParticleMassMsun",
-          function=TranslationFunc("particle_mass"), 
-          particle_type=True, convert_function=_convertParticleMassMsun)
+KnownRAMSESFields["all", "particle_mass"]._units = r"\mathrm{g}"
 
 def _Temperature(field, data):
     rv = data["Pressure"]/data["Density"]
@@ -136,49 +130,6 @@
     return rv
 add_field("Temperature", function=_Temperature, units=r"\rm{K}")
 
-
-# We now set up a couple particle fields.  This should eventually be abstracted
-# into a single particle field function that adds them all on and is used
-# across frontends, but that will need to wait until moving to using
-# Coordinates, or vector fields.
-
-def particle_count(field, data):
-    pos = np.column_stack([data["particle_position_%s" % ax] for ax in 'xyz'])
-    d = data.deposit(pos, method = "count")
-    return d
-RAMSESFieldInfo.add_field(("deposit", "%s_count" % "all"),
-         function = particle_count,
-         validators = [ValidateSpatial()],
-         display_name = "\\mathrm{%s Count}" % "all",
-         projection_conversion = '1')
-
-def particle_mass(field, data):
-    pos = np.column_stack([data["particle_position_%s" % ax] for ax in 'xyz'])
-    d = data.deposit(pos, [data["ParticleMass"]], method = "sum")
-    return d
-
-RAMSESFieldInfo.add_field(("deposit", "%s_mass" % "all"),
-         function = particle_mass,
-         validators = [ValidateSpatial()],
-         display_name = "\\mathrm{%s Mass}" % "all",
-         units = r"\mathrm{g}",
-         projected_units = r"\mathrm{g}\/\mathrm{cm}",
-         projection_conversion = 'cm')
-
-def particle_density(field, data):
-    pos = np.column_stack([data["particle_position_%s" % ax] for ax in 'xyz'])
-    d = data.deposit(pos, [data["ParticleMass"]], method = "sum")
-    d /= data["CellVolume"]
-    return d
-
-RAMSESFieldInfo.add_field(("deposit", "%s_density" % "all"),
-         function = particle_density,
-         validators = [ValidateSpatial()],
-         display_name = "\\mathrm{%s Density}" % "all",
-         units = r"\mathrm{g}/\mathrm{cm}^{3}",
-         projected_units = r"\mathrm{g}/\mathrm{cm}^{-2}",
-         projection_conversion = 'cm')
-
 # We'll add a bunch of species fields here.  In the not too distant future,
 # we'll be moving all of these to a unified field location, so they can be
 # shared between various frontends.
@@ -247,3 +198,10 @@
                   function=_SpeciesNumberDensity,
                   convert_function=_ConvertNumberDensity,
                   validators=ValidateDataField("%s_Density" % species))
+
+# PARTICLE FIELDS
+particle_vector_functions("all", ["particle_position_%s" % ax for ax in 'xyz'],
+                                 ["particle_velocity_%s" % ax for ax in 'xyz'],
+                          RAMSESFieldInfo)
+particle_deposition_functions("all", "Coordinates", "particle_mass",
+                               RAMSESFieldInfo)

diff -r 27f5b715b5cf07575fe4c42bef402e9f9039f798 -r 358f9ca489417c9a0800ec302e6729be9ee34c76 yt/frontends/ramses/io.py
--- a/yt/frontends/ramses/io.py
+++ b/yt/frontends/ramses/io.py
@@ -86,7 +86,7 @@
                 for field in fields:
                     ti = selection.pop(field)[mask]
                     if field not in tr:
-                        dt = subset.domain.particle_field_types[field[1]]
+                        dt = subset.domain.particle_field_types[field]
                         tr[field] = np.empty(size, dt)
                     tr[field][pos:pos+count] = ti
                 pos += count
@@ -98,7 +98,7 @@
         tr = {}
         #for field in sorted(fields, key=lambda a:foffsets[a]):
         for field in fields:
-            f.seek(foffsets[field[1]])
-            dt = subset.domain.particle_field_types[field[1]]
+            f.seek(foffsets[field])
+            dt = subset.domain.particle_field_types[field]
             tr[field] = fpu.read_vector(f, dt)
         return tr


https://bitbucket.org/yt_analysis/yt-3.0/commits/00a361746464/
Changeset:   00a361746464
Branch:      yt-3.0
User:        MatthewTurk
Date:        2013-06-11 17:10:28
Summary:     Convert SPH to use the particle field registry.
Affected #:  1 file

diff -r 358f9ca489417c9a0800ec302e6729be9ee34c76 -r 00a361746464f7adf59927eb04fbbc49bb145d45 yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -37,6 +37,9 @@
     NullFunc, \
     TranslationFunc
 import yt.data_objects.universal_fields
+from yt.data_objects.particle_fields import \
+    particle_deposition_functions, \
+    particle_scalar_functions
 from yt.utilities.physical_constants import \
     mass_sun_cgs
 
@@ -58,98 +61,6 @@
 KnownTipsyFields = FieldInfoContainer()
 add_tipsy_field = KnownTipsyFields.add_field
 
-def _particle_functions(ptype, coord_name, mass_name, registry):
-    def particle_count(field, data):
-        pos = data[ptype, coord_name]
-        d = data.deposit(pos, method = "count")
-        return d
-    registry.add_field(("deposit", "%s_count" % ptype),
-             function = particle_count,
-             validators = [ValidateSpatial()],
-             display_name = "\\mathrm{%s Count}" % ptype,
-             projection_conversion = '1')
-
-    def particle_mass(field, data):
-        pos = data[ptype, coord_name]
-        d = data.deposit(pos, [data[ptype, mass_name]], method = "sum")
-        return d
-
-    registry.add_field(("deposit", "%s_mass" % ptype),
-             function = particle_mass,
-             validators = [ValidateSpatial()],
-             display_name = "\\mathrm{%s Mass}" % ptype,
-             units = r"\mathrm{g}",
-             projected_units = r"\mathrm{g}\/\mathrm{cm}",
-             projection_conversion = 'cm')
-
-    def particle_density(field, data):
-        pos = data[ptype, coord_name]
-        d = data.deposit(pos, [data[ptype, mass_name]], method = "sum")
-        d /= data["CellVolume"]
-        return d
-
-    registry.add_field(("deposit", "%s_density" % ptype),
-             function = particle_density,
-             validators = [ValidateSpatial()],
-             display_name = "\\mathrm{%s Density}" % ptype,
-             units = r"\mathrm{g}/\mathrm{cm}^{3}",
-             projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
-             projection_conversion = 'cm')
-
-    def particle_cic(field, data):
-        pos = data[ptype, coord_name]
-        d = data.deposit(pos, [data[ptype, mass_name]], method = "cic")
-        d /= data["CellVolume"]
-        return d
-
-    registry.add_field(("deposit", "%s_cic" % ptype),
-             function = particle_cic,
-             validators = [ValidateSpatial()],
-             display_name = "\\mathrm{%s CIC Density}" % ptype,
-             units = r"\mathrm{g}/\mathrm{cm}^{3}",
-             projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
-             projection_conversion = 'cm')
-
-    # Now some translation functions.
-
-    registry.add_field((ptype, "ParticleMass"),
-            function = TranslationFunc((ptype, mass_name)),
-            particle_type = True,
-            units = r"\mathrm{g}")
-
-    def _ParticleMassMsun(field, data):
-        return data[ptype, mass_name].copy()
-    def _conv_Msun(data):
-        return 1.0/mass_sun_cgs
-
-    registry.add_field((ptype, "ParticleMassMsun"),
-            function = _ParticleMassMsun,
-            convert_function = _conv_Msun,
-            particle_type = True,
-            units = r"\mathrm{M}_\odot")
-
-    # For 'all', which is a special field, we skip adding a few types.
-    
-    if ptype == "all": return
-
-    # Now we have to set up the various velocity and coordinate things.  In the
-    # future, we'll actually invert this and use the 3-component items
-    # elsewhere, and stop using these.
-    
-    # Note that we pass in _ptype here so that it's defined inside the closure.
-    def _get_coord_funcs(axi, _ptype):
-        def _particle_velocity(field, data):
-            return data[_ptype, "Velocities"][:,axi]
-        def _particle_position(field, data):
-            return data[_ptype, "Coordinates"][:,axi]
-        return _particle_velocity, _particle_position
-    for axi, ax in enumerate("xyz"):
-        v, p = _get_coord_funcs(axi, ptype)
-        registry.add_field((ptype, "particle_velocity_%s" % ax),
-            particle_type = True, function = v)
-        registry.add_field((ptype, "particle_position_%s" % ax),
-            particle_type = True, function = p)
-
 # Here are helper functions for things like vector fields and so on.
 
 def _get_conv(cf):
@@ -191,8 +102,11 @@
         units = r"\mathrm{cm}/\mathrm{s}")
     # Note that we have to do this last so that TranslationFunc operates
     # correctly.
-    _particle_functions(ptype, "Coordinates", "Mass", TipsyFieldInfo)
-_particle_functions("all", "Coordinates", "Mass", TipsyFieldInfo)
+    particle_deposition_functions(ptype, "Coordinates", "Mass",
+                                  TipsyFieldInfo)
+    particle_scalar_functions(ptype, "Coordinates", "Velocities",
+                              TipsyFieldInfo)
+particle_deposition_functions("all", "Coordinates", "Mass", TipsyFieldInfo)
 
 for fname in ["Coordinates", "Velocities", "ParticleIDs", "Mass",
               "Epsilon", "Phi"]:
@@ -246,10 +160,11 @@
         particle_type = True,
         convert_function=_get_conv("velocity"),
         units = r"\mathrm{cm}/\mathrm{s}")
-    _particle_functions(ptype, "Coordinates", "Mass", GadgetFieldInfo)
+    particle_deposition_functions(ptype, "Coordinates", "Mass", GadgetFieldInfo)
+    particle_scalar_functions(ptype, "Coordinates", "Velocities", GadgetFieldInfo)
     KnownGadgetFields.add_field((ptype, "Coordinates"), function=NullFunc,
         particle_type = True)
-_particle_functions("all", "Coordinates", "Mass", GadgetFieldInfo)
+particle_deposition_functions("all", "Coordinates", "Mass", GadgetFieldInfo)
 
 # Now we have to manually apply the splits for "all", since we don't want to
 # use the splits defined above.
@@ -304,10 +219,11 @@
         particle_type = True,
         convert_function=_get_conv("velocity"),
         units = r"\mathrm{cm}/\mathrm{s}")
-    _particle_functions(ptype, "Coordinates", "Mass", OWLSFieldInfo)
+    particle_deposition_functions(ptype, "Coordinates", "Mass", OWLSFieldInfo)
+    particle_scalar_functions(ptype, "Coordinates", "Velocities", OWLSFieldInfo)
     KnownOWLSFields.add_field((ptype, "Coordinates"), function=NullFunc,
         particle_type = True)
-_particle_functions("all", "Coordinates", "Mass", OWLSFieldInfo)
+particle_deposition_functions("all", "Coordinates", "Mass", OWLSFieldInfo)
 
 # Now we have to manually apply the splits for "all", since we don't want to
 # use the splits defined above.


https://bitbucket.org/yt_analysis/yt-3.0/commits/d07242c7c90b/
Changeset:   d07242c7c90b
Branch:      yt-3.0
User:        MatthewTurk
Date:        2013-06-12 00:17:16
Summary:     Adding a better comment and fixing a typo.
Affected #:  1 file

diff -r 00a361746464f7adf59927eb04fbbc49bb145d45 -r d07242c7c90beb454362e7779ce2ae0790f7a69b yt/frontends/enzo/data_structures.py
--- a/yt/frontends/enzo/data_structures.py
+++ b/yt/frontends/enzo/data_structures.py
@@ -492,7 +492,11 @@
             field_list = None
         field_list = self.comm.mpi_bcast(field_list)
         self.field_list = []
-        # Now we will, avoiding the problem of particle types not having names.
+        # Now we will iterate over all fields, trying to avoid the problem of
+        # particle types not having names.  This should convert all known
+        # particle fields that exist in Enzo outputs into the construction
+        # ("all", field) and should not otherwise affect ActiveParticle
+        # simulations.
         for field in field_list:
             if ("all", field) in KnownEnzoFields:
                 self.field_list.append(("all", field))

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