[yt-svn] commit/yt: 6 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Fri May 9 08:38:10 PDT 2014
6 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/b2cff9622c3e/
Changeset: b2cff9622c3e
Branch: yt-3.0
User: brittonsmith
Date: 2014-01-22 12:55:03
Summary: Merging.
Affected #: 20 files
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/analysis_modules/halo_analysis/finding_methods.py
--- /dev/null
+++ b/yt/analysis_modules/halo_analysis/finding_methods.py
@@ -0,0 +1,20 @@
+"""
+Halo Finding methods
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+from .operator_registry import \
+ hf_registry
+
+class HaloFindingMethod(object):
+ pass
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/analysis_modules/halo_analysis/halo_callbacks.py
--- /dev/null
+++ b/yt/analysis_modules/halo_analysis/halo_callbacks.py
@@ -0,0 +1,250 @@
+"""
+Halo callback object
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import h5py
+import numpy as np
+
+from yt.data_objects.profiles import \
+ Profile1D
+from yt.data_objects.yt_array import \
+ YTArray, YTQuantity
+from yt.funcs import \
+ ensure_list
+from yt.utilities.cosmology import \
+ Cosmology
+from yt.utilities.exceptions import YTUnitConversionError
+from yt.utilities.logger import ytLogger as mylog
+
+from .operator_registry import \
+ callback_registry
+
+def add_callback(name, function):
+ callback_registry[name] = HaloCallback(function)
+
+class HaloCallback(object):
+ def __init__(self, function, args=None, kwargs=None):
+ self.function = function
+ self.args = args
+ if self.args is None: self.args = []
+ self.kwargs = kwargs
+ if self.kwargs is None: self.kwargs = {}
+
+ def __call__(self, halo):
+ self.function(halo, *self.args, **self.kwargs)
+ return True
+
+def halo_sphere(halo, radius_field="virial_radius", factor=1.0):
+ r"""
+ Create a sphere data container to associate with a halo.
+
+ Parameters
+ ----------
+ halo : Halo object
+ The Halo object to be provided by the HaloCatalog.
+ radius_field : string
+ Field to be retrieved from the quantities dictionary as
+ the basis of the halo radius.
+ Default: "virial_radius".
+ factor : float
+ Factor to be multiplied by the base radius for defining
+ the radius of the sphere.
+ Defautl: 1.0.
+
+ """
+
+ dpf = halo.halo_catalog.data_pf
+ hpf = halo.halo_catalog.halos_pf
+ center = dpf.arr([halo.quantities["particle_position_%s" % axis] \
+ for axis in "xyz"]) / dpf.length_unit
+ radius = factor * halo.quantities[radius_field] / dpf.length_unit
+ sphere = dpf.h.sphere(center, (radius, "code_length"))
+ setattr(halo, "data_object", sphere)
+
+add_callback("sphere", halo_sphere)
+
+def profile(halo, x_field, y_fields, x_bins=32, x_range=None, x_log=True,
+ weight_field="cell_mass", accumulation=False, storage="profiles"):
+ r"""
+ Create 1d profiles.
+
+ Store profile data in a dictionary associated with the halo object.
+
+ Parameters
+ ----------
+ halo : Halo object
+ The Halo object to be provided by the HaloCatalog.
+ x_field : string
+ The binning field for the profile.
+ y_fields : string or list of strings
+ The fields to be profiled.
+ x_bins : int
+ The number of bins in the profile.
+ Default: 32
+ x_range : (float, float)
+ The range of the x_field. If None, the extrema are used.
+ Default: None
+ x_log : bool
+ Flag for logarithmmic binning.
+ Default: True
+ weight_field : string
+ Weight field for profiling.
+ Default : "cell_mass"
+ accumulation : bool
+ If True, profile data is a cumulative sum.
+ Default : False
+ storage : string
+ Name of the dictionary to store profiles.
+ Default: "profiles"
+
+ """
+
+ mylog.info("Calculating 1D profile for halo %d." %
+ halo.quantities["particle_identifier"])
+
+ dpf = halo.halo_catalog.data_pf
+
+ if x_range is None:
+ x_range = halo.data_object.quantities["Extrema"](x_field)[0]
+ # temporary check until derived quantities are fixed
+ # right now derived quantities do not work with units, so once they do, let us know
+ try:
+ x_range[0]._unit_repr_check_same("cm")
+ raise RuntimeError("Looks like derived quantities have been fixed. Fix this code!")
+ except YTUnitConversionError:
+ # for now, Extrema return dimensionless, but assume it is code_length
+ x_range = [dpf.arr(x.to_ndarray(), "cm") for x in x_range]
+
+ my_profile = Profile1D(halo.data_object, x_field, x_bins,
+ x_range[0], x_range[1], x_log,
+ weight_field=weight_field)
+ my_profile.add_fields(ensure_list(y_fields))
+
+ # temporary fix since profiles do not have units at the moment
+ for field in my_profile.field_data:
+ my_profile.field_data[field] = dpf.arr(my_profile[field],
+ dpf.field_info[field].units)
+
+ # accumulate, if necessary
+ if accumulation:
+ used = my_profile.used
+ for field in my_profile.field_data:
+ if weight_field is None:
+ my_profile.field_data[field][used] = \
+ np.cumsum(my_profile.field_data[field][used])
+ else:
+ my_weight = my_profile.weight[:, 0]
+ my_profile.field_data[field][used] = \
+ np.cumsum(my_profile.field_data[field][used] * my_weight[used]) / \
+ np.cumsum(my_weight[used])
+
+ # create storage dictionary
+ prof_store = dict([(field, my_profile[field]) \
+ for field in my_profile.field_data])
+ prof_store[my_profile.x_field] = my_profile.x
+ if "used" in prof_store:
+ prof_store["used"] &= my_profile.used
+ else:
+ prof_store["used"] = my_profile.used
+
+ setattr(halo, storage, prof_store)
+
+add_callback("profile", profile)
+
+def virial_quantities(halo, fields, critical_overdensity=200,
+ profile_storage="profiles"):
+ r"""
+ Calculate the value of the given fields at the virial radius defined at
+ the given critical density by interpolating from radial profiles.
+
+ Parameters
+ ----------
+ halo : Halo object
+ The Halo object to be provided by the HaloCatalog.
+ fields : (field, units) tuple or list of tuples
+ The fields whose virial values are to be calculated.
+ critical_density : float
+ The value of the overdensity at which to evaulate the virial quantities.
+ Overdensity is with respect to the critical density.
+ Default: 200
+ profile_storage : string
+ Name of the halo attribute that holds the profiles to be used.
+ Default: "profiles"
+
+ """
+
+ mylog.info("Calculating virial quantities for halo %d." %
+ halo.quantities["particle_identifier"])
+
+ fields = ensure_list(fields)
+ for field in fields:
+ q_tuple = ("%s_%d" % (field[0], critical_overdensity), "callback")
+ if q_tuple not in halo.halo_catalog.quantities:
+ halo.halo_catalog.quantities.append(q_tuple)
+
+ dpf = halo.halo_catalog.data_pf
+ co = Cosmology(hubble_constant=dpf.hubble_constant,
+ omega_matter=dpf.omega_matter,
+ omega_lambda=dpf.omega_lambda,
+ unit_registry=dpf.unit_registry)
+ profile_data = getattr(halo, profile_storage)
+
+ if ("index", "cell_volume") not in profile_data or \
+ ("gas", "matter_mass") not in profile_data:
+ raise RuntimeError('virial_quantities callback requires profiles of ("index", "cell_volume") and ("gas", "matter_mass").')
+
+ total_volume = profile_data[("index", "cell_volume")]
+ total_mass = profile_data[("gas", "matter_mass")]
+ overdensity = total_mass / total_volume / co.critical_density(dpf.current_redshift)
+ dfilter = np.isfinite(overdensity) & profile_data["used"] & (overdensity > 0)
+ print overdensity
+ print overdensity[dfilter]
+
+ vquantities = dict([("%s_%d" % (field[0], critical_overdensity), 0) \
+ for field in fields])
+
+ if dfilter.sum() < 2:
+ halo.quantities.update(vquantities)
+ return
+
+ # find interpolation index
+ # require a negative slope, but not monotonicity
+ vod = overdensity[dfilter]
+ if (vod > critical_overdensity).all():
+ if vod[-1] < vod[-2]:
+ index = -2
+ else:
+ halo.quantities.update(vquantities)
+ return
+ elif (vod < critical_overdensity).all():
+ if vod[0] > vod[1]:
+ index = 0
+ else:
+ halo.quantities.update(vquantities)
+ return
+ else:
+ # take first instance of downward intersection with critical value
+ index = np.where((vod[:-1] >= critical_overdensity) &
+ (vod[1:] < critical_overdensity))[0][0]
+
+ for field in fields:
+ v_prof = profile_data[field[0]]
+ slope = (v_prof[index + 1] - v_prof[index]) / (vod[index + 1] - vod[index])
+ vquantities["%s_%d" % (field[0], critical_overdensity)] = \
+ (slope * (critical_overdensity - vod[index]) + v_prof[index]).in_units(field[1])
+ print field[0], vquantities["%s_%d" % (field[0], critical_overdensity)]
+
+ halo.quantities.update(vquantities)
+
+add_callback("virial_quantities", virial_quantities)
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/analysis_modules/halo_analysis/halo_catalog.py
--- /dev/null
+++ b/yt/analysis_modules/halo_analysis/halo_catalog.py
@@ -0,0 +1,147 @@
+"""
+HaloCatalog object
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import h5py
+import numpy as np
+
+from yt.funcs import \
+ mylog
+from yt.utilities.parallel_tools.parallel_analysis_interface import \
+ ParallelAnalysisInterface, \
+ parallel_blocking_call, \
+ parallel_root_only, \
+ parallel_objects
+
+from .halo_object import \
+ Halo
+from .operator_registry import \
+ callback_registry, \
+ quantity_registry, \
+ hf_registry
+
+class HaloCatalog(ParallelAnalysisInterface):
+ r"""Create a HaloCatalog: an object that allows for the creation and association
+ of data with a set of halo objects.
+
+ Parameters
+ ----------
+ halos_pf : str
+ Dataset created by a halo finder. If None, a halo finder should be
+ provided with the finder_method keyword.
+ data_pf : str
+ Dataset created by a simulation.
+ data_source : data container
+ Data container associated with either the halos_pf or the data_pf.
+ finder_method : str
+ Halo finder to be used if no halos_pf is given.
+
+ Examples
+ --------
+
+ """
+
+ def __init__(self, halos_pf=None, data_pf=None,
+ data_source=None, finder_method=None):
+ ParallelAnalysisInterface.__init__(self)
+ self.halos_pf = halos_pf
+ self.data_pf = data_pf
+ self.finder_method = finder_method
+
+ if halos_pf is None:
+ if data_pf is None:
+ raise RuntimeError("Must specify a halos_pf, data_pf, or both.")
+ if finder_method is None:
+ raise RuntimeError("Must specify a halos_pf or a finder_method.")
+ if data_source is None:
+ if halos_pf is not None:
+ data_source = halos_pf.h.all_data()
+ else:
+ data_source = data_pf.h.all_data()
+ self.data_source = data_source
+
+ self.add_default_quantities()
+ self.callbacks = []
+
+ def add_callback(self, callback, *args, **kwargs):
+ callback = callback_registry.find(callback, *args, **kwargs)
+ self.callbacks.append(callback)
+
+ def add_quantity(self, key, field_type=None, *args, **kwargs):
+ if field_type is None:
+ quantity = quantity_registry.find(key, *args, **kwargs)
+ elif (field_type, key) in self.halos_pf.field_info:
+ quantity = (field_type, key)
+ else:
+ raise RuntimeError("HaloCatalog quantity must be a registered function or a field of a known type.")
+ self.quantities.append((key, quantity))
+
+ def add_filter(self, filter, *args, **kwargs):
+ filter = callback_registry.find(filter, *args, **kwargs)
+ self.callbacks.append(filter)
+
+ def run(self, njobs=-1, filename=None):
+ self.halo_list = []
+
+ if self.halos_pf is None:
+ # this is where we would do halo finding and assign halos_pf to
+ # the dataset that we have just created.
+ raise NotImplementedError
+
+ self.n_halos = self.data_source["particle_identifier"].size
+ for i in parallel_objects(xrange(self.n_halos), njobs=njobs):
+ new_halo = Halo(self)
+ for key, quantity in self.quantities:
+ if quantity in self.halos_pf.field_info:
+ new_halo.quantities[key] = self.data_source[quantity][i]
+ elif callable(quantity):
+ new_halo.quantities[key] = quantity(new_halo)
+
+ for callback in self.callbacks:
+ callback(new_halo)
+
+ self.halo_list.append(new_halo.quantities)
+
+ del new_halo
+
+ self.halo_list = self.comm.par_combine_object(self.halo_list,
+ datatype="list",
+ op="cat")
+ self.halo_list.sort(key=lambda a:a['particle_identifier'].to_ndarray())
+
+ if filename is not None:
+ self.save_catalog(filename)
+
+ @parallel_root_only
+ def save_catalog(self, filename):
+ "Write out hdf5 file with all halo quantities."
+
+ mylog.info("Saving halo catalog to %s." % filename)
+ out_file = h5py.File(filename, 'w')
+ field_data = np.empty(self.n_halos)
+ for key, quantity in self.quantities:
+ for i in xrange(self.n_halos):
+ field_data[i] = self.halo_list[i][key]
+ out_file.create_dataset(key, data=field_data)
+ out_file.close()
+
+ def add_default_quantities(self):
+ self.quantities = []
+ self.add_quantity("particle_identifier", field_type="halos")
+ self.add_quantity("particle_mass", field_type="halos")
+ self.add_quantity("particle_position_x", field_type="halos")
+ self.add_quantity("particle_position_y", field_type="halos")
+ self.add_quantity("particle_position_z", field_type="halos")
+ self.add_quantity("virial_radius", field_type="halos")
+
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/analysis_modules/halo_analysis/halo_filters.py
--- /dev/null
+++ b/yt/analysis_modules/halo_analysis/halo_filters.py
@@ -0,0 +1,23 @@
+"""
+Halo filter object
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+from .halo_callbacks import HaloCallback
+
+class HaloFilter(HaloCallback):
+ def __init__(self, function, args, kwargs):
+ HaloCallback.__init__(self, function, args, kwargs)
+
+ def __call__(self, halo_catalog, halo):
+ return self.function(halo_catalog, halo, *self.args, **self.kwargs)
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/analysis_modules/halo_analysis/halo_object.py
--- /dev/null
+++ b/yt/analysis_modules/halo_analysis/halo_object.py
@@ -0,0 +1,20 @@
+"""
+Halo object.
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+class Halo(object):
+ particles = None
+ def __init__(self, halo_catalog):
+ self.halo_catalog = halo_catalog
+ self.quantities = {}
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/analysis_modules/halo_analysis/halo_quantities.py
--- /dev/null
+++ b/yt/analysis_modules/halo_analysis/halo_quantities.py
@@ -0,0 +1,51 @@
+"""
+Halo quantity object
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+
+from .halo_callbacks import HaloCallback
+from .operator_registry import quantity_registry
+
+def add_quantity(name, function):
+ quantity_registry[name] = HaloQuantity(function)
+
+class HaloQuantity(HaloCallback):
+ def __init__(self, function, *args, **kwargs):
+ HaloCallback.__init__(self, function, args, kwargs)
+
+ def __call__(self, halo):
+ return self.function(halo, *self.args, **self.kwargs)
+
+def center_of_mass(halo):
+ if halo.particles is None:
+ raise RuntimeError("Center of mass requires halo to have particle data.")
+ return (halo.particles['particle_mass'] *
+ np.array([halo.particles['particle_position_x'],
+ halo.particles['particle_position_y'],
+ halo.particles['particle_position_z']])).sum(axis=1) / \
+ halo.particles['particle_mass'].sum()
+
+add_quantity('center_of_mass', center_of_mass)
+
+def bulk_velocity(halo):
+ if halo.particles is None:
+ raise RuntimeError("Bulk velocity requires halo to have particle data.")
+ return (halo.particles['particle_mass'] *
+ np.array([halo.particles['particle_velocity_x'],
+ halo.particles['particle_velocity_y'],
+ halo.particles['particle_velocity_z']])).sum(axis=1) / \
+ halo.particles['particle_mass'].sum()
+
+add_quantity('bulk_velocity', bulk_velocity)
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/analysis_modules/halo_analysis/operator_registry.py
--- /dev/null
+++ b/yt/analysis_modules/halo_analysis/operator_registry.py
@@ -0,0 +1,29 @@
+"""
+Operation registry class
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import types
+
+class OperatorRegistry(dict):
+ def find(self, op, *args, **kwargs):
+ if isinstance(op, types.StringTypes):
+ # Lookup, assuming string or hashable object
+ op = self[op]
+ op.args = args
+ op.kwargs = kwargs
+ return op
+
+callback_registry = OperatorRegistry()
+hf_registry = OperatorRegistry()
+quantity_registry = OperatorRegistry()
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/analysis_modules/halo_finding/rockstar/rockstar.py
--- a/yt/analysis_modules/halo_finding/rockstar/rockstar.py
+++ b/yt/analysis_modules/halo_finding/rockstar/rockstar.py
@@ -143,7 +143,7 @@
the width of the smallest grid element in the simulation from the
last data snapshot (i.e. the one where time has evolved the
longest) in the time series:
- ``pf_last.h.get_smallest_dx() * pf_last['mpch']``.
+ ``pf_last.h.get_smallest_dx().in_units("Mpc/h")``.
total_particles : int
If supplied, this is a pre-calculated total number of particles present
in the simulation. For example, this is useful when analyzing a series
@@ -203,7 +203,7 @@
self.outbase = outbase
if force_res is None:
tpf = ts[-1] # Cache a reference
- self.force_res = tpf.h.get_smallest_dx() * tpf['mpch']
+ self.force_res = tpf.h.get_smallest_dx().in_units("Mpc/h")
# We have to delete now to wipe the hierarchy
del tpf
else:
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/analysis_modules/halo_finding/rockstar/rockstar_interface.pyx
--- a/yt/analysis_modules/halo_finding/rockstar/rockstar_interface.pyx
+++ b/yt/analysis_modules/halo_finding/rockstar/rockstar_interface.pyx
@@ -198,7 +198,7 @@
p[0] = <particle *> malloc(sizeof(particle) * local_parts)
- conv[0] = conv[1] = conv[2] = pf["mpchcm"]
+ conv[0] = conv[1] = conv[2] = pf.length_unit.in_units("Mpccm/h")
conv[3] = conv[4] = conv[5] = 1e-5
left_edge[0] = pf.domain_left_edge[0]
left_edge[1] = pf.domain_left_edge[1]
@@ -301,7 +301,7 @@
PARTICLE_MASS = particle_mass
PERIODIC = periodic
BOX_SIZE = (tpf.domain_right_edge[0] -
- tpf.domain_left_edge[0]) * tpf['mpchcm']
+ tpf.domain_left_edge[0]).in_units("Mpccm/h")
setup_config()
rh = self
cdef LPG func = rh_read_particles
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/data_objects/profiles.py
--- a/yt/data_objects/profiles.py
+++ b/yt/data_objects/profiles.py
@@ -19,6 +19,8 @@
from yt.funcs import *
from yt.units.yt_array import uconcatenate
+from yt.data_objects.yt_array import \
+ YTArray, YTQuantity
from yt.data_objects.data_containers import YTFieldData
from yt.utilities.lib import bin_profile1d, bin_profile2d, bin_profile3d
from yt.utilities.lib import new_bin_profile1d, new_bin_profile2d, \
@@ -771,9 +773,12 @@
# We use our main comm here
# This also will fill _field_data
# FIXME: Add parallelism and combining std stuff
+ blank = ~temp_storage.used
+ self.used = temp_storage.used
if self.weight_field is not None:
temp_storage.values /= temp_storage.weight_values[...,None]
- blank = ~temp_storage.used
+ self.weight = temp_storage.weight_values[...,None]
+ self.weight[blank] = 0.0
for i, field in enumerate(fields):
self.field_data[field] = temp_storage.values[...,i]
self.field_data[field][blank] = 0.0
@@ -827,7 +832,9 @@
super(Profile1D, self).__init__(data_source, weight_field)
self.x_field = x_field
self.x_log = x_log
- self.x_bins = self._get_bins(x_min, x_max, x_n, x_log)
+ self.x_bins = data_source.pf.arr(self._get_bins(x_min.to_ndarray(),
+ x_max.to_ndarray(), x_n, x_log),
+ x_min.units)
self.size = (self.x_bins.size - 1,)
self.bin_fields = (self.x_field,)
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -212,10 +212,10 @@
max_nu = 1e30
good_u = None
for unit in ['mpc', 'kpc', 'pc', 'au', 'rsun', 'km', 'cm']:
- vv = v*self[unit]
+ vv = v * self.length_unit.in_units(unit)
if vv < max_nu and vv > 1.0:
good_u = unit
- max_nu = v*self[unit]
+ max_nu = v * self.length_unit.in_units(unit)
if good_u is None : good_u = 'cm'
return good_u
@@ -428,7 +428,6 @@
self.set_code_units()
-
def get_unit_from_registry(self, unit_str):
"""
Creates a unit object matching the string expression, using this
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/fields/api.py
--- a/yt/fields/api.py
+++ b/yt/fields/api.py
@@ -17,6 +17,7 @@
field_plugins
from . import angular_momentum
+from . import cosmology_fields
from . import fluid_fields
from . import magnetic_field
from . import geometric_fields
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/fields/cosmology_fields.py
--- /dev/null
+++ b/yt/fields/cosmology_fields.py
@@ -0,0 +1,128 @@
+"""
+Cosmology related fields.
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2014, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+
+from .field_plugin_registry import \
+ register_field_plugin
+
+from yt.utilities.cosmology import \
+ Cosmology
+from yt.utilities.physical_constants import \
+ speed_of_light_cgs
+
+ at register_field_plugin
+def setup_cosmology_fields(registry, ftype = "gas", slice_info = None):
+ # slice_info would be the left, the right, and the factor.
+ # For example, with the old Enzo-ZEUS fields, this would be:
+ # slice(None, -2, None)
+ # slice(1, -1, None)
+ # 1.0
+ # Otherwise, we default to a centered difference.
+ if slice_info is None:
+ sl_left = slice(None, -2, None)
+ sl_right = slice(2, None, None)
+ div_fac = 2.0
+ else:
+ sl_left, sl_right, div_face = slice_info
+
+ def _matter_density(field, data):
+ return data[ftype, "density"] + \
+ data[ftype, "dark_matter_density"]
+
+ registry.add_field((ftype, "matter_density"),
+ function=_matter_density,
+ units="g/cm**3")
+
+ def _matter_mass(field, data):
+ return data[ftype, "matter_density"] * data["index", "cell_volume"]
+
+ registry.add_field((ftype, "matter_mass"),
+ function=_matter_mass,
+ units="g")
+
+ # rho_total / rho_cr(z).
+ def _overdensity(field, data):
+ # consider moving cosmology object to pf attribute
+ co = Cosmology(hubble_constant=data.pf.hubble_constant,
+ omega_matter=data.pf.omega_matter,
+ omega_lambda=data.pf.omega_lambda)
+ return data["matter_density"] / co.critical_density(data.pf.current_redshift)
+
+ registry.add_field((ftype, "overdensity"),
+ function=_overdensity,
+ units="")
+
+ # rho_baryon / <rho_baryon>
+ def _baryon_overdensity(field, data):
+ omega_baryon = data.get_field_parameter("omega_baryon")
+ # consider moving cosmology object to pf attribute
+ co = Cosmology(hubble_constant=data.pf.hubble_constant,
+ omega_matter=data.pf.omega_matter,
+ omega_lambda=data.pf.omega_lambda,
+ unit_registry=data.pf.unit_registry)
+ # critical_density(z) ~ omega_lambda + omega_matter * (1 + z)^3
+ # mean density(z) ~ omega_matter * (1 + z)^3
+ return data["density"] / omega_baryon_now / co.critical_density(0.0) / \
+ (1.0 + data.pf.hubble_constant)**3
+
+ registry.add_field("baryon_overdensity",
+ function=_baryon_overdensity,
+ units="")
+
+ # rho_matter / <rho_matter>
+ def _matter_overdensity(field, data):
+ # consider moving cosmology object to pf attribute
+ co = Cosmology(hubble_constant=data.pf.hubble_constant,
+ omega_matter=data.pf.omega_matter,
+ omega_lambda=data.pf.omega_lambda,
+ unit_registry=data.pf.unit_registry)
+ # critical_density(z) ~ omega_lambda + omega_matter * (1 + z)^3
+ # mean density(z) ~ omega_matter * (1 + z)^3
+ return data["density"] / data.pf.omega_matter / co.critical_density(0.0) / \
+ (1.0 + data.pf.hubble_constant)**3
+
+ registry.add_field("matter_overdensity",
+ function=_matter_overdensity,
+ units="")
+
+ # Weak lensing convergence.
+ # Eqn 4 of Metzler, White, & Loken (2001, ApJ, 547, 560).
+ # This needs to be checked for accuracy.
+ def _weak_lensing_convergence(field, data):
+ # consider moving cosmology object to pf attribute
+ co = Cosmology(hubble_constant=data.pf.hubble_constant,
+ omega_matter=data.pf.omega_matter,
+ omega_lambda=data.pf.omega_lambda,
+ unit_registry=data.pf.unit_registry)
+ observer_redshift = data.get_field_parameter('observer_redshift')
+ source_redshift = data.get_field_parameter('source_redshift')
+
+ # observer to lens
+ dl = co.angular_diameter_distance(observer_redshift, data.pf.current_redshift)
+ # observer to source
+ ds = co.angular_diameter_distance(observer_redshift, source_redshift)
+ # lens to source
+ dls = co.angular_diameter_distance(data.pf.current_redshift, source_redshift)
+
+ # removed the factor of 1 / a to account for the fact that we are projecting
+ # with a proper distance.
+ return 1.5 * (co.hubble_constant / speed_of_light_cgs)**2 * (dl * dls / ds) * \
+ data["matter_overdensity"]
+
+ registry.add_field("weak_lensing_convergence",
+ function=_weak_lensing_convergence,
+ units="1/cm")
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/fields/universal_fields.py
--- a/yt/fields/universal_fields.py
+++ b/yt/fields/universal_fields.py
@@ -82,63 +82,6 @@
add_field("jeans_mass", function=jeans_mass, units="g")
-# This is rho_total / rho_cr(z).
-def _overdensity(field, data):
- return data["matter_density"] / (rho_crit_now * data.pf.hubble_constant**2 *
- (1+data.pf.current_redshift)**3)
-add_field("overdensity", function=_overdensity)
-
-# This is rho_matter / <rho_matter> - 1.0
-def _density_perturbation(field, data):
- omega_m = data.pf.omega_matter
- h = data.pf.hubble_constant
- z = data.pf.current_redshift
- rho_m = rho_crit_now * h**2 * omega_m * (1.0 + z)**3
- return data["matter_density"] / rho_m - 1.0
-
-add_field("density_perturbation", function=_overdensity)
-
-# This is rho_baryon / <rho_baryon>
-def _baryon_overdensity(field, data):
- # @todo: should we provide this field if the dataset doesn't have omega_b?
- if data.pf.has_key('omega_baryon_now'):
- omega_baryon_now = data.pf['omega_baryon_now']
- else:
- omega_baryon_now = 0.0441
-
- # These are enzo parameters and should be changed.
- return data["density"] / (omega_baryon_now * rho_crit_now *
- (data.pf["CosmologyHubbleConstantNow"]**2) *
- ((1.0 + data.pf["CosmologyCurrentRedshift"])**3))
-
-add_field("baryon_overdensity", function=_baryon_overdensity)
-
-#FIXME
-
-# Weak lensing convergence.
-# Eqn 4 of Metzler, White, & Loken (2001, ApJ, 547, 560).
-#def _convertConvergence(data):
-# if not data.pf.parameters.has_key('cosmology_calculator'):
-# data.pf.parameters['cosmology_calculator'] = Cosmology(
-# HubbleConstantNow=(100.*data.pf.hubble_constant),
-# OmegaMatterNow=data.pf.omega_matter, OmegaLambdaNow=data.pf.omega_lambda)
-# # observer to lens
-# DL = data.pf.parameters['cosmology_calculator'].AngularDiameterDistance(
-# data.pf.parameters['observer_redshift'], data.pf.current_redshift)
-# # observer to source
-# DS = data.pf.parameters['cosmology_calculator'].AngularDiameterDistance(
-# data.pf.parameters['observer_redshift'], data.pf.parameters['lensing_source_redshift'])
-# # lens to source
-# DLS = data.pf.parameters['cosmology_calculator'].AngularDiameterDistance(
-# data.pf.current_redshift, data.pf.parameters['lensing_source_redshift'])
-# # TODO: convert 1.5e14 to constants
-# return (((DL * DLS) / DS) * (1.5e14 * data.pf.omega_matter *
-# (data.pf.hubble_constant / speed_of_light_cgs)**2 *
-# (1 + data.pf.current_redshift)))
-#add_field("weak_lensing_convergence", function=_overdensity,
-# convert_function=_convertConvergence,
-# projection_conversion='mpccm')
-
### Begin block that should probably be in an analysis module ###
def _chandra_emissivity(field, data):
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -62,7 +62,7 @@
("RadAccel1", (ra_units, ["radiation_acceleration_x"], None)),
("RadAccel2", (ra_units, ["radiation_acceleration_y"], None)),
("RadAccel3", (ra_units, ["radiation_acceleration_z"], None)),
- ("Dark_Matter_Mass", (rho_units, ["dark_matter_mass"], None)),
+ ("Dark_Matter_Density", (rho_units, ["dark_matter_density"], None)),
("Temperature", ("K", ["temperature"], None)),
("Dust_Temperature", ("K", ["dust_temperature"], None)),
("x-velocity", (vel_units, ["velocity_x"], None)),
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/geometry/selection_routines.pyx
--- a/yt/geometry/selection_routines.pyx
+++ b/yt/geometry/selection_routines.pyx
@@ -1595,3 +1595,18 @@
return ("always", 1,)
always_selector = AlwaysSelector
+
+cdef class HaloParticlesSelector(SelectorObject):
+ cdef public object base_source
+ cdef SelectorObject base_selector
+ cdef object pind
+ cdef public np.int64_t halo_id
+ def __init__(self, dobj):
+ self.base_source = dobj.base_source
+ self.base_selector = self.base_source.selector
+ self.pind = dobj.particle_indices
+
+ def _hash_vals(self):
+ return ("halo_particles", self.halo_id)
+
+halo_particles_selector = HaloParticlesSelector
diff -r c0389c170da8570d5b1d1d97025977542e0b3316 -r b2cff9622c3e4a8f5d205e7bb84f07cb68033f50 yt/units/unit_object.py
--- a/yt/units/unit_object.py
+++ b/yt/units/unit_object.py
@@ -431,7 +431,7 @@
"""
if symbol_str in unit_symbol_lut:
# lookup successful, return the tuple directly
- return unit_symbol_lut[symbol_str]
+ return unit_symbol_lut[symbol_str][0:2]
# could still be a known symbol with a prefix
possible_prefix = symbol_str[0]
https://bitbucket.org/yt_analysis/yt/commits/824549a00c0e/
Changeset: 824549a00c0e
Branch: yt-3.0
User: brittonsmith
Date: 2014-03-17 15:01:47
Summary: Adding gravitational potential and fixing energy units.
Affected #: 1 file
diff -r 732be20f7ed12269689d3e5040c7e9e2acbce121 -r 824549a00c0e9c0078940d2f183950f4cf5f9899 yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -83,6 +83,7 @@
("z-velocity", (vel_units, ["velocity_z"], None)),
("RaySegments", ("", ["ray_segments"], None)),
("PhotoGamma", (ra_units, ["photo_gamma"], None)),
+ ("PotentialField", ("code_velocity**2", ["gravitational_potential"], None)),
("Density", (rho_units, ["density"], None)),
("Metal_Density", (rho_units, ["metal_density"], None)),
("SN_Colour", (rho_units, [], None)),
@@ -173,13 +174,13 @@
if self.pf.parameters["HydroMethod"] == 2:
self.add_output_field(("enzo", te_name),
- units="code_length**2/code_time**2")
+ units="code_velocity**2")
self.alias(("gas", "thermal_energy"), ("enzo", te_name))
elif self.pf.parameters["DualEnergyFormalism"] == 1:
self.add_output_field(
("enzo", ge_name),
- units="code_length**2/code_time**2")
+ units="code_velocity**2")
self.alias(
("gas", "thermal_energy"),
("enzo", ge_name),
@@ -187,7 +188,7 @@
elif self.pf.parameters["HydroMethod"] in (4, 6):
self.add_output_field(
("enzo", te_name),
- units="code_length**2/code_time**2")
+ units="code_velocity**2")
# Subtract off B-field energy
def _sub_b(field, data):
return data[te_name] - 0.5*(
@@ -201,7 +202,7 @@
else: # Otherwise, we assume TotalEnergy is kinetic+thermal
self.add_output_field(
("enzo", te_name),
- units = "code_length**2/code_time**2")
+ units = "code_velocity**2")
def _tot_minus_kin(field, data):
return data[te_name] - 0.5*(
data["x-velocity"]**2.0
https://bitbucket.org/yt_analysis/yt/commits/e80c4b8ba92f/
Changeset: e80c4b8ba92f
Branch: yt-3.0
User: brittonsmith
Date: 2014-05-05 17:57:28
Summary: Merging.
Affected #: 1 file
https://bitbucket.org/yt_analysis/yt/commits/1b91abfc9f84/
Changeset: 1b91abfc9f84
Branch: yt-3.0
User: brittonsmith
Date: 2014-05-05 20:36:11
Summary: No-op merging old head.
Affected #: 3 files
https://bitbucket.org/yt_analysis/yt/commits/d90d4a360d16/
Changeset: d90d4a360d16
Branch: yt-3.0
User: brittonsmith
Date: 2014-05-09 17:29:49
Summary: Merging.
Affected #: 4 files
https://bitbucket.org/yt_analysis/yt/commits/d6681af78917/
Changeset: d6681af78917
Branch: yt-3.0
User: brittonsmith
Date: 2014-05-09 17:35:43
Summary: Adding subfind frontend to setup.
Affected #: 1 file
diff -r d90d4a360d164adc4ebe3383152f21e7d559455a -r d6681af789175edcbcf4360dd03b1a8ec726ea25 yt/frontends/halo_catalogs/setup.py
--- a/yt/frontends/halo_catalogs/setup.py
+++ b/yt/frontends/halo_catalogs/setup.py
@@ -5,6 +5,7 @@
def configuration(parent_package='', top_path=None):
config = Configuration('halo_catalogs', parent_package, top_path)
config.add_subpackage("halo_catalog")
+ config.add_subpackage("owls_subfind")
config.add_subpackage("rockstar")
config.make_config_py()
return config
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