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

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Wed Jun 29 11:06:47 PDT 2016


1 new commit in yt:

https://bitbucket.org/yt_analysis/yt/commits/9b2d8e598c95/
Changeset:   9b2d8e598c95
Branch:      yt
User:        xarthisius
Date:        2016-06-29 18:06:34+00:00
Summary:     Merged in ngoldbaum/yt (pull request #2213)

MKS code unit system fixes. Closes #1229
Affected #:  11 files

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -2,7 +2,7 @@
   local_artio_000:
     - yt/frontends/artio/tests/test_outputs.py
 
-  local_athena_000:
+  local_athena_001:
     - yt/frontends/athena
 
   local_chombo_000:

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -43,6 +43,7 @@
     ParameterFileStore, \
     NoParameterShelf, \
     output_type_registry
+from yt.units.dimensions import current_mks
 from yt.units.unit_object import Unit, unit_system_registry
 from yt.units.unit_registry import UnitRegistry
 from yt.fields.derived_field import \
@@ -237,15 +238,19 @@
         self.no_cgs_equiv_length = False
 
         self._create_unit_registry()
-        self._parse_parameter_file()
-        self.set_units()
-        self._setup_coordinate_handler()
 
         create_code_unit_system(self)
         if unit_system == "code":
             unit_system = str(self)
+        else:
+            unit_system = str(unit_system).lower()
+
         self.unit_system = unit_system_registry[unit_system]
 
+        self._parse_parameter_file()
+        self.set_units()
+        self._setup_coordinate_handler()
+
         # Because we need an instantiated class to check the ds's existence in
         # the cache, we move that check to here from __new__.  This avoids
         # double-instantiation.
@@ -865,18 +870,30 @@
         self._set_code_unit_attributes()
         # here we override units, if overrides have been provided.
         self._override_code_units()
+
         self.unit_registry.modify("code_length", self.length_unit)
         self.unit_registry.modify("code_mass", self.mass_unit)
         self.unit_registry.modify("code_time", self.time_unit)
         if hasattr(self, 'magnetic_unit'):
-            # If we do not have this set, but some fields come in in
-            # "code_magnetic", this will allow them to remain in that unit.
+            # if the magnetic unit is in T, we need to recreate the code unit
+            # system as an MKS-like system
+            if current_mks in self.magnetic_unit.units.dimensions.free_symbols:
+
+                if self.unit_system == unit_system_registry[str(self)]:
+                    unit_system_registry.pop(str(self))
+                    create_code_unit_system(self, current_mks_unit='A')
+                    self.unit_system = unit_system_registry[str(self)]
+                elif str(self.unit_system) == 'mks':
+                    pass
+                else:
+                    self.magnetic_unit = \
+                        self.magnetic_unit.to_equivalent('gauss', 'CGS')
             self.unit_registry.modify("code_magnetic", self.magnetic_unit)
         vel_unit = getattr(
             self, "velocity_unit", self.length_unit / self.time_unit)
         pressure_unit = getattr(
             self, "pressure_unit",
-            self.mass_unit / (self.length_unit * self.time_unit)**2)
+            self.mass_unit / (self.length_unit * (self.time_unit)**2))
         temperature_unit = getattr(self, "temperature_unit", 1.0)
         density_unit = getattr(self, "density_unit", self.mass_unit / self.length_unit**3)
         self.unit_registry.modify("code_velocity", vel_unit)

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d yt/frontends/athena/data_structures.py
--- a/yt/frontends/athena/data_structures.py
+++ b/yt/frontends/athena/data_structures.py
@@ -327,11 +327,6 @@
         self.dataset.domain_dimensions = \
                 np.round(self.dataset.domain_width/gdds[0]).astype('int')
 
-        # Need to reset the units in the dataset based on the correct
-        # domain left/right/dimensions.
-        # DEV: Is this really necessary?
-        #self.dataset._set_code_unit_attributes()
-
         if self.dataset.dimensionality <= 2 :
             self.dataset.domain_dimensions[2] = np.int(1)
         if self.dataset.dimensionality == 1 :
@@ -483,19 +478,10 @@
             # We set these to cgs for now, but they may be overridden later.
             mylog.warning("Assuming 1.0 = 1.0 %s", cgs)
             setattr(self, "%s_unit" % unit, self.quan(1.0, cgs))
-
-    def set_code_units(self):
-        super(AthenaDataset, self).set_code_units()
-        mag_unit = getattr(self, "magnetic_unit", None)
-        vel_unit = getattr(self, "velocity_unit", None)
-        if mag_unit is None:
-            self.magnetic_unit = np.sqrt(4*np.pi * self.mass_unit /
-                                         (self.time_unit**2 * self.length_unit))
+        self.magnetic_unit = np.sqrt(4*np.pi * self.mass_unit /
+                                     (self.time_unit**2 * self.length_unit))
         self.magnetic_unit.convert_to_units("gauss")
-        self.unit_registry.modify("code_magnetic", self.magnetic_unit)
-        if vel_unit is None:
-            self.velocity_unit = self.length_unit/self.time_unit
-        self.unit_registry.modify("code_velocity", self.velocity_unit)
+        self.velocity_unit = self.length_unit / self.time_unit
 
     def _parse_parameter_file(self):
         self._handle = open(self.parameter_filename, "rb")

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d yt/frontends/athena/fields.py
--- a/yt/frontends/athena/fields.py
+++ b/yt/frontends/athena/fields.py
@@ -19,7 +19,7 @@
     kboltz, mh
 
 b_units = "code_magnetic"
-pres_units = "code_mass/(code_length*code_time**2)"
+pres_units = "code_pressure"
 erg_units = "code_mass * (code_length/code_time)**2"
 rho_units = "code_mass / code_length**3"
 

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d yt/frontends/flash/data_structures.py
--- a/yt/frontends/flash/data_structures.py
+++ b/yt/frontends/flash/data_structures.py
@@ -253,12 +253,9 @@
         self.time_unit = self.quan(1.0, "s")
         self.velocity_unit = self.quan(1.0, "cm/s")
         self.temperature_unit = self.quan(temperature_factor, "K")
-        self.unit_registry.modify("code_magnetic", self.magnetic_unit)
-        
+
     def set_code_units(self):
         super(FLASHDataset, self).set_code_units()
-        self.unit_registry.modify("code_temperature",
-                                  self.temperature_unit.value)
 
     def _find_parameter(self, ptype, pname, scalar = False):
         nn = "/%s %s" % (ptype,

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d yt/frontends/gdf/data_structures.py
--- a/yt/frontends/gdf/data_structures.py
+++ b/yt/frontends/gdf/data_structures.py
@@ -26,6 +26,10 @@
     GridIndex
 from yt.data_objects.static_output import \
     Dataset
+from yt.units.dimensions import \
+    dimensionless as sympy_one
+from yt.units.unit_object import \
+    Unit
 from yt.utilities.exceptions import \
     YTGDFUnknownGeometry
 from yt.utilities.lib.misc_utilities import \
@@ -211,12 +215,22 @@
                 current_unit = h5f["/dataset_units/%s" % unit_name]
                 value = current_unit.value
                 unit = current_unit.attrs["unit"]
+                # need to convert to a Unit object and check dimensions
+                # because unit can be things like
+                # 'dimensionless/dimensionless**3' so naive string
+                # comparisons are insufficient
+                unit = Unit(unit, registry=self.unit_registry)
+                if unit_name.endswith('_unit') and unit.dimensions is sympy_one:
+                    un = unit_name[:-5]
+                    un = un.replace('magnetic', 'magnetic_field', 1)
+                    unit = self.unit_system[un]
+                    setattr(self, unit_name, self.quan(value, unit))
                 setattr(self, unit_name, self.quan(value, unit))
                 if unit_name in h5f["/field_types"]:
                     if unit_name in self.field_units:
                         mylog.warning("'field_units' was overridden by 'dataset_units/%s'"
                                       % (unit_name))
-                    self.field_units[unit_name] = unit.decode('utf8')
+                    self.field_units[unit_name] = str(unit)
         else:
             self.length_unit = self.quan(1.0, "cm")
             self.mass_unit = self.quan(1.0, "g")

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d yt/units/tests/test_unit_systems.py
--- a/yt/units/tests/test_unit_systems.py
+++ b/yt/units/tests/test_unit_systems.py
@@ -15,7 +15,7 @@
 from yt.units.unit_systems import UnitSystem
 from yt.units import dimensions
 from yt.convenience import load
-from yt.testing import assert_almost_equal, requires_file
+from yt.testing import assert_almost_equal, assert_allclose, requires_file
 from yt.config import ytcfg
 
 def test_unit_systems():
@@ -46,7 +46,7 @@
                           "velocity_divergence": "1/Myr",
                           "density_gradient_x": "Msun/kpc**4"}
 test_units["code"] = {"density": "code_mass/code_length**3",
-                      "kinetic_energy": "code_mass/(code_length*code_time**2)",
+                      "kinetic_energy": "code_pressure",
                       "velocity_magnitude": "code_velocity",
                       "velocity_divergence": "code_velocity/code_length",
                       "density_gradient_x": "code_mass/code_length**4"}
@@ -96,4 +96,44 @@
                 v1 = dd_cgs[field].in_base(us)
             v2 = dd[field]
             assert_almost_equal(v1.v, v2.v)
-            assert str(v2.units) == test_units[us][field]
\ No newline at end of file
+            assert str(v2.units) == test_units[us][field]
+
+wdm = 'WDMerger_hdf5_chk_1000/WDMerger_hdf5_chk_1000.hdf5'
+ at requires_file(wdm)
+def test_tesla_magnetic_unit():
+    ytcfg["yt", "skip_dataset_cache"] = "True"
+    for us in ['cgs', 'mks', 'code']:
+        ds = load(wdm, unit_system=us,
+                  units_override={'magnetic_unit': (1.0, 'T')})
+        ad = ds.all_data()
+        dens = ad['density']
+        magx = ad['magx']
+        magnetic_field_x = ad['magnetic_field_x']
+
+        if us == 'cgs':
+            assert str(dens.units) == 'g/cm**3'
+            assert str(magx.units) == 'code_magnetic'
+            assert magx.uq == ds.quan(1e4, 'G')
+            assert str(magnetic_field_x.units) == 'gauss'
+            assert_allclose(magx.value, magnetic_field_x.value/1e4)
+            assert_allclose(
+                magnetic_field_x.to_equivalent('T', 'SI').value,
+                magnetic_field_x.value/1e4)
+
+        if us == 'mks':
+            assert str(dens.units) == 'kg/m**3'
+            assert str(magx.units) == 'code_magnetic'
+            assert magx.uq == ds.quan(1, 'T')
+            assert str(magnetic_field_x.units) == 'T'
+            assert_allclose(magx.value, magnetic_field_x.value)
+            assert_allclose(magnetic_field_x.to_equivalent('G', 'CGS').value,
+                            magnetic_field_x.value*1e4)
+
+        if us == 'code':
+            assert str(dens.units) == 'code_mass/code_length**3'
+            assert str(magx.units) == 'code_magnetic'
+            assert magx.uq == ds.quan(1, 'T')
+            assert str(magnetic_field_x.units) == 'code_magnetic'
+            assert_allclose(magx.value, magnetic_field_x.value)
+            assert_allclose(magnetic_field_x.to_equivalent('G', 'CGS').value,
+                            magnetic_field_x.value*1e4)

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d yt/units/unit_object.py
--- a/yt/units/unit_object.py
+++ b/yt/units/unit_object.py
@@ -26,6 +26,8 @@
     base_dimensions, temperature, \
     dimensionless, current_mks, \
     angle
+from yt.units.equivalencies import \
+    equivalence_registry
 from yt.units.unit_lookup_table import \
     unit_prefixes, prefixable_units, latex_prefixes, \
     default_base_units
@@ -432,6 +434,25 @@
                 return False
         return True
 
+    def list_equivalencies(self):
+        """
+        Lists the possible equivalencies associated with this unit object
+        """
+        for k, v in equivalence_registry.items():
+            if self.has_equivalent(k):
+                print(v())
+
+    def has_equivalent(self, equiv):
+        """
+        Check to see if this unit object as an equivalent unit in *equiv*.
+        """
+        try:
+            this_equiv = equivalence_registry[equiv]()
+        except KeyError:
+            raise KeyError("No such equivalence \"%s\"." % equiv)
+        old_dims = self.dimensions
+        return old_dims in this_equiv.dims
+
     def get_base_equivalent(self, unit_system="cgs"):
         """
         Create and return dimensionally-equivalent units in a specified base.
@@ -652,6 +673,8 @@
     my_dims = dimensions.expand()
     if my_dims is dimensionless:
         return ""
+    if my_dims in base_units:
+        return base_units[my_dims]
     for factor in my_dims.as_ordered_factors():
         dim = list(factor.free_symbols)[0]
         unit_string = str(base_units[dim])

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d yt/units/unit_registry.py
--- a/yt/units/unit_registry.py
+++ b/yt/units/unit_registry.py
@@ -93,11 +93,13 @@
                 "in this registry." % symbol)
 
         del self.lut[symbol]
+        if symbol in self.unit_objs:
+            del self.unit_objs[symbol]
 
     def modify(self, symbol, base_value):
         """
-        Change the base value of a dimension.  Useful for adjusting code units
-        after parsing parameters."
+        Change the base value of a unit symbol.  Useful for adjusting code units
+        after parsing parameters.
 
         """
         if symbol not in self.lut:
@@ -106,9 +108,17 @@
                 "in this registry." % symbol)
 
         if hasattr(base_value, "in_base"):
-            base_value = base_value.in_base().value
+            new_dimensions = base_value.units.dimensions
+            base_value = base_value.in_base('cgs-ampere')
+            base_value = base_value.value
+        else:
+            new_dimensions = self.lut[symbol][1]
 
-        self.lut[symbol] = (float(base_value), ) + self.lut[symbol][1:]
+        self.lut[symbol] = ((float(base_value), new_dimensions) +
+                            self.lut[symbol][2:])
+
+        if symbol in self.unit_objs:
+            del self.unit_objs[symbol]
 
     def keys(self):
         """

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d yt/units/unit_systems.py
--- a/yt/units/unit_systems.py
+++ b/yt/units/unit_systems.py
@@ -79,10 +79,9 @@
 
     def __getitem__(self, key):
         if isinstance(key, string_types):
-            if key not in self._dims:
-                self._dims.append(key)
             key = getattr(dimensions, key)
-        if key not in self.units_map:
+        um = self.units_map
+        if key not in um or um[key].dimensions is not key:
             units = _get_system_unit_string(key, self.units_map)
             self.units_map[key] = Unit(units, registry=self.registry)
         return self.units_map[key]
@@ -109,10 +108,16 @@
                 repr += "  %s: %s\n" % (key, self.units_map[dim])
         return repr
 
-def create_code_unit_system(ds):
-    code_unit_system = UnitSystem(str(ds), "code_length", "code_mass", "code_time",
-                                  "code_temperature", registry=ds.unit_registry)
+def create_code_unit_system(ds, current_mks_unit=None):
+    code_unit_system = UnitSystem(
+        str(ds), "code_length", "code_mass", "code_time", "code_temperature",
+        current_mks_unit=current_mks_unit, registry=ds.unit_registry)
     code_unit_system["velocity"] = "code_velocity"
+    if current_mks_unit:
+        code_unit_system["magnetic_field_mks"] = "code_magnetic"
+    else:
+        code_unit_system["magnetic_field_cgs"] = "code_magnetic"
+    code_unit_system["pressure"] = "code_pressure"
 
 cgs_unit_system = UnitSystem("cgs", "cm", "g", "s")
 cgs_unit_system["energy"] = "erg"
@@ -147,3 +152,14 @@
 planck_unit_system = UnitSystem("planck", "l_pl", "m_pl", "t_pl", temperature_unit="T_pl")
 planck_unit_system["energy"] = "E_pl"
 planck_unit_system["charge_cgs"] = "q_pl"
+
+
+cgs_ampere_unit_system = UnitSystem('cgs-ampere', 'cm', 'g', 's',
+                                    current_mks_unit='A')
+cgs_ampere_unit_system["energy"] = "erg"
+cgs_ampere_unit_system["specific_energy"] = "erg/g"
+cgs_ampere_unit_system["pressure"] = "dyne/cm**2"
+cgs_ampere_unit_system["force"] = "dyne"
+cgs_ampere_unit_system["magnetic_field_cgs"] = "gauss"
+cgs_ampere_unit_system["charge_cgs"] = "esu"
+cgs_ampere_unit_system["current_cgs"] = "statA"

diff -r 3215234b6e95186ebf3a1858937e0e476f2b260f -r 9b2d8e598c9587b5508dc0c9e5530c32fccc329d yt/units/yt_array.py
--- a/yt/units/yt_array.py
+++ b/yt/units/yt_array.py
@@ -596,18 +596,21 @@
             The unit that you wish to convert to.
         equiv : string
             The equivalence you wish to use. To see which equivalencies are
-            supported for this unitful quantity, try the :meth:`list_equivalencies`
-            method.
+            supported for this unitful quantity, try the
+            :meth:`list_equivalencies` method.
 
         Examples
         --------
         >>> a = yt.YTArray(1.0e7,"K")
         >>> a.to_equivalent("keV", "thermal")
         """
-        unit_quan = YTQuantity(1.0, unit, registry=self.units.registry)
+        conv_unit = Unit(unit, registry=self.units.registry)
         this_equiv = equivalence_registry[equiv]()
-        if self.has_equivalent(equiv) and (unit_quan.has_equivalent(equiv) or this_equiv._one_way):
-            new_arr = this_equiv.convert(self, unit_quan.units.dimensions, **kwargs)
+        oneway_or_equivalent = (
+            conv_unit.has_equivalent(equiv) or this_equiv._one_way)
+        if self.has_equivalent(equiv) and oneway_or_equivalent:
+            new_arr = this_equiv.convert(
+                self, conv_unit.dimensions, **kwargs)
             if isinstance(new_arr, tuple):
                 try:
                     return YTArray(new_arr[0], new_arr[1]).in_units(unit)
@@ -623,21 +626,14 @@
         Lists the possible equivalencies associated with this YTArray or
         YTQuantity.
         """
-        for k,v in equivalence_registry.items():
-            if self.has_equivalent(k):
-                print(v())
+        self.units.list_equivalencies()
 
     def has_equivalent(self, equiv):
         """
         Check to see if this YTArray or YTQuantity has an equivalent unit in
         *equiv*.
         """
-        try:
-            this_equiv = equivalence_registry[equiv]()
-        except KeyError:
-            raise KeyError("No such equivalence \"%s\"." % equiv)
-        old_dims = self.units.dimensions
-        return old_dims in this_equiv.dims
+        return self.units.has_equivalent(equiv)
 
     def ndarray_view(self):
         """

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