[Yt-svn] yt-commit r1641 - trunk/yt/lagos

sskory at wrangler.dreamhost.com sskory at wrangler.dreamhost.com
Thu Feb 25 12:43:26 PST 2010


Author: sskory
Date: Thu Feb 25 12:43:25 2010
New Revision: 1641
URL: http://yt.enzotools.org/changeset/1641

Log:
By request, adding a RMS velocity calculation to the halo finding module. It works for all kinds of halos, and in particular, in parallel with parallelHF halos.

Modified:
   trunk/yt/lagos/HaloFinding.py

Modified: trunk/yt/lagos/HaloFinding.py
==============================================================================
--- trunk/yt/lagos/HaloFinding.py	(original)
+++ trunk/yt/lagos/HaloFinding.py	Thu Feb 25 12:43:25 2010
@@ -58,7 +58,7 @@
 
     def __init__(self, halo_list, id, indices = None, size=None, CoM=None,
         max_dens_point=None, group_total_mass=None, max_radius=None, bulk_vel=None,
-        tasks=None):
+        tasks=None, rms_vel=None):
         self.halo_list = halo_list
         self.id = id
         self.data = halo_list._data_source
@@ -75,6 +75,7 @@
         self.max_radius = max_radius
         self.bulk_vel = bulk_vel
         self.tasks = tasks
+        self.rms_vel = rms_vel
         self.bin_count = None
         self.overdensity = None
 
@@ -121,6 +122,18 @@
         vz = (self["particle_velocity_z"] * pm).sum()
         return na.array([vx,vy,vz])/pm.sum()
 
+    def rms_velocity(self):
+        """
+        Returns the RMS velocity for the halo particles in code units.
+        """
+        bv = self.bulk_velocity()
+        vx = self["particle_velocity_x"] - bv[0]
+        vy = self["particle_velocity_y"] - bv[1]
+        vz = self["particle_velocity_z"] - bv[2]
+        s = vx**2 + vy**2 + vz**2
+        ms = na.mean(s)
+        return na.sqrt(ms)
+
     def maximum_radius(self, center_of_mass=True):
         """
         Returns the maximum radius in the halo for all particles,
@@ -271,7 +284,8 @@
     dont_wrap = ["maximum_density","maximum_density_location",
         "center_of_mass","total_mass","bulk_velocity","maximum_radius",
         "get_size","get_sphere", "write_particle_list","__getitem__", 
-        "virial_info", "virial_bin", "virial_mass", "virial_radius"]
+        "virial_info", "virial_bin", "virial_mass", "virial_radius",
+        "rms_velocity"]
 
     def maximum_density(self):
         """
@@ -363,6 +377,27 @@
         global_bv = self._mpi_allsum(bv)
         return global_bv[:3]/global_bv[3]
 
+    def rms_velocity(self):
+        """
+        Returns the RMS velocity for the halo particles in code units.
+        """
+        if self.rms_vel is not None:
+            return self.rms_vel
+        bv = self.bulk_velocity()
+        if self.indices is not None:
+            vx = self["particle_velocity_x"] - bv[0]
+            vy = self["particle_velocity_y"] - bv[1]
+            vz = self["particle_velocity_z"] - bv[2]
+            s = vx**2 + vy**2 + vz**2
+            s = na.sum(s)
+            size = vx.size
+            ss = na.array([s, float(size)])
+        else:
+            ss = na.array([0.,0.])
+        global_ss = self._mpi_allsum(ss)
+        ms = global_ss[0] / global_ss[1]
+        return na.sqrt(ms)
+
     def maximum_radius(self, center_of_mass=True):
         """
         Returns the maximum radius in the halo for all particles,
@@ -695,7 +730,7 @@
         f.write("\t".join(["# Group","Mass","# part","max dens"
                            "x","y","z", "center-of-mass",
                            "x","y","z",
-                           "vx","vy","vz","max_r","\n"]))
+                           "vx","vy","vz","max_r","rms_v","\n"]))
         for group in self:
             f.write("%10i\t" % group.id)
             f.write("%0.9e\t" % group.total_mass())
@@ -708,6 +743,7 @@
             f.write("\t".join(["%0.9e" % v for v in group.bulk_velocity()]))
             f.write("\t")
             f.write("%0.9e\t" % group.maximum_radius())
+            f.write("%0.9e\t" % group.rms_velocity())
             f.write("\n")
             f.flush()
         f.close()
@@ -850,12 +886,9 @@
         yt_counters("bulk vel. reading data")
         pm = self.particle_fields["ParticleMassMsun"]
         if ytcfg.getboolean("yt","inline") == False:
-            xv = self._data_source.particles["particle_velocity_x"][self._base_indices] * \
-                self._data_source.convert("x-velocity")
-            yv = self._data_source.particles["particle_velocity_y"][self._base_indices] * \
-                self._data_source.convert("y-velocity")
-            zv = self._data_source.particles["particle_velocity_z"][self._base_indices] * \
-                self._data_source.convert("z-velocity")
+            xv = self._data_source.particles["particle_velocity_x"][self._base_indices]
+            yv = self._data_source.particles["particle_velocity_y"][self._base_indices]
+            zv = self._data_source.particles["particle_velocity_z"][self._base_indices]
         else:
             xv = self._data_source["particle_velocity_x"][self._base_indices]
             yv = self._data_source["particle_velocity_y"][self._base_indices]
@@ -881,14 +914,39 @@
             marks = na.concatenate(([0], marks, [calc]))
             for i, u in enumerate(uniq_subchain):
                 self.bulk_vel[u] = na.sum(vel[marks[i]:marks[i+1]], axis=0)
-            del select, vel, ms, subchain, sort, sort_subchain, uniq_subchain
-            del diff_subchain, marks
+            del vel, ms, subchain, sort_subchain
+            del diff_subchain
         # Bring it together, and divide by the previously computed total mass
         # of each halo.
         self.bulk_vel = self._mpi_Allsum_double(self.bulk_vel)
         for groupID in xrange(self.group_count):
             self.bulk_vel[groupID] = self.bulk_vel[groupID] / self.Tot_M[groupID]
         yt_counters("bulk vel. computing")
+        # Now calculate the RMS velocity of the groups in parallel, very
+        # similarly to the bulk velocity and re-using some of the arrays.
+        yt_counters("rms vel computing")
+        rms_vel_temp = na.zeros((self.group_count,2), dtype='float64')
+        if calc:
+            vel = na.empty((calc, 3), dtype='float64')
+            vel[:,0] = xv[select]
+            vel[:,1] = yv[select]
+            vel[:,2] = zv[select]
+            vel = vel[sort]
+            vel = vel**2 # Squared
+            for i, u in enumerate(uniq_subchain):
+                # This finds the Mean
+                rms_vel_temp[u][0] = na.sum(vel[marks[i]:marks[i+1]] - self.bulk_vel[u])
+                rms_vel_temp[u][1] = marks[i+1] - marks[i]
+            del vel, marks, uniq_subchain
+        # Bring it together.
+        rms_vel_temp = self._mpi_Allsum_double(rms_vel_temp)
+        self.rms_vel = na.empty(self.group_count, dtype='float64')
+        for groupID in xrange(self.group_count):
+            # Here we do the Root.
+            self.rms_vel[groupID] = \
+                na.sqrt(rms_vel_temp[groupID][0] / rms_vel_temp[groupID][1])
+        del rms_vel_temp
+        yt_counters("rms vel computing")
         self.taskID = obj.mine
         self.halo_taskmap = obj.halo_taskmap # A defaultdict.
         del obj
@@ -923,7 +981,8 @@
                     size=self.group_sizes[index], CoM=self.CoM[index], \
                     max_dens_point=self.max_dens_point[index], \
                     group_total_mass=self.Tot_M[index], max_radius=self.max_radius[index],
-                    bulk_vel=self.bulk_vel[index], tasks=self.halo_taskmap[index])
+                    bulk_vel=self.bulk_vel[index], tasks=self.halo_taskmap[index],
+                    rms_vel=self.rms_vel[index])
                 # I don't own this halo
                 self._do_not_claim_object(self._groups[index])
                 self._max_dens[index] = [self.max_dens_point[index][0], self.max_dens_point[index][1], \
@@ -935,7 +994,8 @@
                 size=self.group_sizes[i], CoM=self.CoM[i], \
                 max_dens_point=self.max_dens_point[i], \
                 group_total_mass=self.Tot_M[i], max_radius=self.max_radius[i],
-                bulk_vel=self.bulk_vel[i], tasks=self.halo_taskmap[index])
+                bulk_vel=self.bulk_vel[i], tasks=self.halo_taskmap[index],
+                rms_vel=self.rms_vel[i])
             # This halo may be owned by many, including this task
             self._claim_object(self._groups[index])
             self._max_dens[index] = [self.max_dens_point[i][0], self.max_dens_point[i][1], \
@@ -948,14 +1008,15 @@
                 size=self.group_sizes[index], CoM=self.CoM[index], \
                 max_dens_point=self.max_dens_point[i], \
                 group_total_mass=self.Tot_M[index], max_radius=self.max_radius[index],
-                bulk_vel=self.bulk_vel[index], tasks=self.halo_taskmap[index])
+                bulk_vel=self.bulk_vel[index], tasks=self.halo_taskmap[index],
+                rms_vel=self.rms_vel[index])
             self._do_not_claim_object(self._groups[index])
             self._max_dens[index] = [self.max_dens_point[index][0], self.max_dens_point[index][1], \
                 self.max_dens_point[index][2], self.max_dens_point[index][3]]
             index += 1
         # Clean up
         del self.max_dens_point, self.max_radius, self.bulk_vel
-        del self.halo_taskmap, self.tags
+        del self.halo_taskmap, self.tags, self.rms_vel
 
     def __len__(self):
         return self.group_count



More information about the yt-svn mailing list