[yt-svn] commit/yt: brittonsmith: Merged in ngoldbaum/yt (pull request #1787)

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Thu Oct 22 02:34:49 PDT 2015


1 new commit in yt:

https://bitbucket.org/yt_analysis/yt/commits/a5d76f647fbb/
Changeset:   a5d76f647fbb
Branch:      yt
User:        brittonsmith
Date:        2015-10-22 09:34:37+00:00
Summary:     Merged in ngoldbaum/yt (pull request #1787)

Tipsy cosmology unit fixes
Affected #:  3 files

diff -r d77e881944f88d7d4fce3c90b4fec115cf0180ce -r a5d76f647fbb9a740e21cf5d8269660aea6736cc yt/frontends/tipsy/data_structures.py
--- a/yt/frontends/tipsy/data_structures.py
+++ b/yt/frontends/tipsy/data_structures.py
@@ -23,6 +23,7 @@
 
 from yt.frontends.sph.data_structures import \
     ParticleDataset
+from yt.funcs import deprecate
 from yt.geometry.particle_geometry_handler import \
     ParticleIndex
 from yt.data_objects.static_output import \
@@ -167,12 +168,14 @@
         periodic = self.parameters.get('bPeriodic', True)
         period = self.parameters.get('dPeriod', None)
         self.periodicity = (periodic, periodic, periodic)
-        self.comoving = self.parameters.get('bComove', False)
-        if self.comoving and period is None:
+        self.cosmological_simulation = float(self.parameters.get(
+            'bComove', self._cosmology_parameters is not None))
+        if self.cosmological_simulation and period is None:
             period = 1.0
         if self.bounding_box is None:
             if periodic and period is not None:
-                # If we are periodic, that sets our domain width to either 1 or dPeriod.
+                # If we are periodic, that sets our domain width to
+                # either 1 or dPeriod.
                 self.domain_left_edge = np.zeros(3, "float64") - 0.5*period
                 self.domain_right_edge = np.zeros(3, "float64") + 0.5*period
             else:
@@ -187,19 +190,22 @@
 
         # If the cosmology parameters dictionary got set when data is
         # loaded, we can assume it's a cosmological data set
-        if self.comoving or self._cosmology_parameters is not None:
+        if self.cosmological_simulation == 1.0:
             cosm = self._cosmology_parameters or {}
-            self.scale_factor = hvals["time"]#In comoving simulations, time stores the scale factor a
-            self.cosmological_simulation = 1
-            dcosm = dict(current_redshift=(1.0/self.scale_factor)-1.0,
-                         omega_lambda=self.parameters.get('dLambda', cosm.get('omega_lambda',0.0)),
-                         omega_matter=self.parameters.get('dOmega0', cosm.get('omega_matter',0.0)),
-                         hubble_constant=self.parameters.get('dHubble0', cosm.get('hubble_constant',1.0)))
+            # In comoving simulations, time stores the scale factor a
+            self.scale_factor = hvals["time"]
+            dcosm = dict(
+                current_redshift=(1.0/self.scale_factor)-1.0,
+                omega_lambda=self.parameters.get(
+                    'dLambda', cosm.get('omega_lambda',0.0)),
+                omega_matter=self.parameters.get(
+                    'dOmega0', cosm.get('omega_matter',0.0)),
+                hubble_constant=self.parameters.get(
+                    'dHubble0', cosm.get('hubble_constant',1.0)))
             for param in dcosm.keys():
                 pval = dcosm[param]
                 setattr(self, param, pval)
         else:
-            self.cosmological_simulation = 0.0
             kpc_unit = self.parameters.get('dKpcUnit', 1.0)
             self._unit_base['cm'] = 1.0 / (kpc_unit * cm_per_kpc)
 
@@ -218,51 +224,57 @@
         super(TipsyDataset, self)._set_derived_attrs()
 
     def _set_code_unit_attributes(self):
+        # First try to set units based on parameter file
         if self.cosmological_simulation:
             mu = self.parameters.get('dMsolUnit', 1.)
+            self.mass_unit = self.quan(mu, 'Msun')
             lu = self.parameters.get('dKpcUnit', 1000.)
             # In cosmological runs, lengths are stored as length*scale_factor
             self.length_unit = self.quan(lu, 'kpc')*self.scale_factor
-            self.mass_unit = self.quan(mu, 'Msun')
-            density_unit = self.mass_unit/ (self.length_unit/self.scale_factor)**3
-
-            # If self.comoving is set, we know this is a gasoline data set,
-            # and we do the conversion on the hubble constant.
-            if self.comoving:
-                # Gasoline's hubble constant, dHubble0, is stored units of
-                # proper code time.
-                self.hubble_constant *= np.sqrt(G.in_units(
-                    'kpc**3*Msun**-1*s**-2') * density_unit).value / (
-                    3.2407793e-18)
-            cosmo = Cosmology(self.hubble_constant,
-                              self.omega_matter, self.omega_lambda)
-            self.current_time = cosmo.hubble_time(self.current_redshift)
+            density_unit = self.mass_unit / (self.length_unit / self.scale_factor)**3
+            if 'dHubble0' in self.parameters:
+                # Gasoline's internal hubble constant, dHubble0, is stored in
+                # units of proper code time
+                self.hubble_constant *= np.sqrt(G * density_unit)
+                # Finally, we scale the hubble constant by 100 km/s/Mpc
+                self.hubble_constant /= self.quan(100, 'km/s/Mpc')
+                # If we leave it as a YTQuantity, the cosmology object
+                # used below will add units back on.
+                self.hubble_constant = self.hubble_constant.in_units("").d
         else:
             mu = self.parameters.get('dMsolUnit', 1.0)
             self.mass_unit = self.quan(mu, 'Msun')
             lu = self.parameters.get('dKpcUnit', 1.0)
             self.length_unit = self.quan(lu, 'kpc')
-            density_unit = self.mass_unit / self.length_unit**3
-        self.time_unit = 1.0 / np.sqrt(G * density_unit)
 
         # If unit base is defined by the user, override all relevant units
         if self._unit_base is not None:
-            length = self._unit_base.get('length', self.length_unit)
-            length = self.quan(*length) if isinstance(length, tuple) else self.quan(length)
-            self.length_unit = length
+            for my_unit in ["length", "mass", "time"]:
+                if my_unit in self._unit_base:
+                    my_val = self._unit_base[my_unit]
+                    my_val = \
+                      self.quan(*my_val) if isinstance(my_val, tuple) \
+                      else self.quan(my_val)
+                    setattr(self, "%s_unit" % my_unit, my_val)
 
-            mass = self._unit_base.get('mass', self.mass_unit)
-            mass = self.quan(*mass) if isinstance(mass, tuple) else self.quan(mass)
-            self.mass_unit = mass
+        # Finally, set the dependent units
+        if self.cosmological_simulation:
+            cosmo = Cosmology(self.hubble_constant,
+                              self.omega_matter, self.omega_lambda)
+            self.current_time = cosmo.hubble_time(self.current_redshift)
+            # mass units are rho_crit(z=0) * domain volume
+            mu = cosmo.critical_density(0.0) * \
+              (1 + self.current_redshift)**3 * self.length_unit**3
+            self.mass_unit = self.quan(mu.in_units("Msun"), "Msun")
+            density_unit = self.mass_unit / (self.length_unit / self.scale_factor)**3
+            # need to do this again because we've modified the hubble constant
+            self.unit_registry.modify("h", self.hubble_constant)
+        else:
+            density_unit = self.mass_unit / self.length_unit**3
 
-            density_unit = self.mass_unit / self.length_unit**3
+        if not hasattr(self, "time_unit"):
             self.time_unit = 1.0 / np.sqrt(G * density_unit)
 
-            time = self._unit_base.get('time', self.time_unit)
-            time = self.quan(*time) if isinstance(time, tuple) else self.quan(time)
-            self.time_unit = time
-
-
     @staticmethod
     def _validate_header(filename):
         '''
@@ -304,3 +316,8 @@
     @classmethod
     def _is_valid(self, *args, **kwargs):
         return TipsyDataset._validate_header(args[0])[0]
+
+    @property
+    @deprecate(replacement='cosmological_simulation')
+    def comoving(self):
+        return self.cosmological_simulation == 1.0

diff -r d77e881944f88d7d4fce3c90b4fec115cf0180ce -r a5d76f647fbb9a740e21cf5d8269660aea6736cc yt/frontends/tipsy/tests/test_outputs.py
--- a/yt/frontends/tipsy/tests/test_outputs.py
+++ b/yt/frontends/tipsy/tests/test_outputs.py
@@ -42,7 +42,7 @@
                                 hubble_constant = 0.702)
     kwargs = dict(field_dtypes = {"Coordinates": "d"},
                   cosmology_parameters = cosmology_parameters,
-                  unit_base = {'length': (1.0/60.0, "Mpccm/h")},
+                  unit_base = {'length': (60.0, "Mpccm/h")},
                   n_ref = 64)
     ds = data_dir_load(pkdgrav, TipsyDataset, (), kwargs)
     yield assert_equal, str(ds), "halo1e11_run1.00400"
@@ -73,7 +73,7 @@
                                 omega_matter = 0.272,
                                 hubble_constant = 0.702)
     kwargs = dict(cosmology_parameters = cosmology_parameters,
-                  unit_base = {'length': (1.0/60.0, "Mpccm/h")},
+                  unit_base = {'length': (60.0, "Mpccm/h")},
                   n_ref = 64)
     ds = data_dir_load(gasoline_dmonly, TipsyDataset, (), kwargs)
     yield assert_equal, str(ds), "agora_1e11.00400"

diff -r d77e881944f88d7d4fce3c90b4fec115cf0180ce -r a5d76f647fbb9a740e21cf5d8269660aea6736cc yt/funcs.py
--- a/yt/funcs.py
+++ b/yt/funcs.py
@@ -221,24 +221,29 @@
     if ytcfg.getint("yt", "__topcomm_parallel_rank") > 0: return
     mylog.info(*args)
 
-def deprecate(func):
-    """
-    This decorator issues a deprecation warning.
+def deprecate(replacement):
+    def real_deprecate(func):
+        """
+        This decorator issues a deprecation warning.
 
-    This can be used like so:
+        This can be used like so:
 
-    .. code-block:: python
+        .. code-block:: python
 
-       @deprecate
-       def some_really_old_function(...):
+        @deprecate("new_function")
+        def some_really_old_function(...):
 
-    """
-    @wraps(func)
-    def run_func(*args, **kwargs):
-        warnings.warn("%s has been deprecated and may be removed without notice!" \
-                % func.__name__, DeprecationWarning, stacklevel=2)
-        func(*args, **kwargs)
-    return run_func
+        """
+        @wraps(func)
+        def run_func(*args, **kwargs):
+            message = "%s has been deprecated and may be removed without notice!"
+            if replacement is not None:
+                message += " Use %s instead." % replacement
+            warnings.warn(message % func.__name__, DeprecationWarning,
+                          stacklevel=2)
+            func(*args, **kwargs)
+        return run_func
+    return real_deprecate
 
 def pdb_run(func):
     """

Repository URL: https://bitbucket.org/yt_analysis/yt/

--

This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.



More information about the yt-svn mailing list