[yt-svn] commit/yt: 5 new changesets
Bitbucket
commits-noreply at bitbucket.org
Thu Oct 18 12:21:08 PDT 2012
5 new commits in yt:
https://bitbucket.org/yt_analysis/yt/changeset/ae6617541438/
changeset: ae6617541438
branch: yt
user: MatthewTurk
date: 2012-10-18 19:13:32
summary: Adding a routine to test fields. However, the aspect of making the parameter
file "realistic" points out how some fields are probably misplaced as
"universal." I've also updated a few that rely on cosmological parameters to
use the attributes, not the dict-style access.
The test_fields.py routine should eventually be extended to handle particle
fields, but fake_random_pf does not yet do that.
Additionally, many of the tests currently fail because of buffer size
mismatches in the coordinate transformations.
affected #: 3 files
diff -r 06462f53a16bbbb175a88920ed0c219aacbdc4d9 -r ae66175414388a369d488d58fab308f84e80203b yt/data_objects/field_info_container.py
--- a/yt/data_objects/field_info_container.py
+++ b/yt/data_objects/field_info_container.py
@@ -160,7 +160,8 @@
# required attrs
pf = fake_parameter_file(lambda: 1)
pf.current_redshift = pf.omega_lambda = pf.omega_matter = \
- pf.hubble_constant = pf.cosmological_simulation = 0.0
+ pf.cosmological_simulation = 0.0
+ pf.hubble_constant = 0.7
pf.domain_left_edge = np.zeros(3, 'float64')
pf.domain_right_edge = np.ones(3, 'float64')
pf.dimensionality = 3
diff -r 06462f53a16bbbb175a88920ed0c219aacbdc4d9 -r ae66175414388a369d488d58fab308f84e80203b yt/data_objects/tests/test_fields.py
--- /dev/null
+++ b/yt/data_objects/tests/test_fields.py
@@ -0,0 +1,90 @@
+from yt.testing import *
+import numpy as np
+from yt.data_objects.field_info_container import \
+ FieldInfo
+import yt.data_objects.universal_fields
+from yt.utilities.definitions import \
+ mpc_conversion, sec_conversion
+
+def setup():
+ from yt.config import ytcfg
+ ytcfg["yt","__withintesting"] = "True"
+ np.seterr(all = 'ignore')
+
+_sample_parameters = dict(
+ axis = 0,
+ center = np.array((0.0, 0.0, 0.0)),
+ bulk_velocity = np.array((0.0, 0.0, 0.0)),
+ normal = np.array((0.0, 0.0, 1.0)),
+ cp_x_vec = np.array((1.0, 0.0, 0.0)),
+ cp_y_vec = np.array((0.0, 1.0, 0.0)),
+ cp_z_vec = np.array((0.0, 0.0, 1.0)),
+)
+
+_base_fields = ["Density", "x-velocity", "y-velocity", "z-velocity"]
+
+def realistic_pf(fields, nprocs):
+ pf = fake_random_pf(16, fields = fields, nprocs = nprocs)
+ pf.parameters["HydroMethod"] = "streaming"
+ pf.parameters["Gamma"] = 5.0/3.0
+ pf.parameters["EOSType"] = 1.0
+ pf.parameters["EOSSoundSpeed"] = 1.0
+ pf.conversion_factors["Time"] = 1.0
+ pf.conversion_factors.update( dict((f, 1.0) for f in fields) )
+ pf.current_redshift = 0.0001
+ pf.hubble_constant = 0.7
+ for unit in mpc_conversion:
+ pf.units[unit+'h'] = pf.units[unit]
+ pf.units[unit+'cm'] = pf.units[unit]
+ pf.units[unit+'hcm'] = pf.units[unit]
+ return pf
+
+class TestFieldAccess(object):
+ description = None
+
+ def __init__(self, field_name, nproc):
+ # Note this should be a field name
+ self.field_name = field_name
+ self.description = "Accessing_%s_%s" % (field_name, nproc)
+ self.nproc = nproc
+
+ def __call__(self):
+ field = FieldInfo[self.field_name]
+ deps = field.get_dependencies()
+ fields = deps.requested + _base_fields
+ skip_grids = False
+ needs_spatial = False
+ for v in field.validators:
+ f = getattr(v, "fields", None)
+ if f: fields += f
+ if getattr(v, "ghost_zones", 0) > 0:
+ skip_grids = True
+ if hasattr(v, "ghost_zones"):
+ needs_spatial = True
+ pf = realistic_pf(fields, self.nproc)
+ # This gives unequal sized grids as well as subgrids
+ dd1 = pf.h.all_data()
+ dd2 = pf.h.all_data()
+ dd1.field_parameters.update(_sample_parameters)
+ dd2.field_parameters.update(_sample_parameters)
+ v1 = dd1[self.field_name]
+ conv = field._convert_function(dd1) or 1.0
+ if not needs_spatial:
+ assert_equal(v1, conv*field._function(field, dd2))
+ if not skip_grids:
+ for g in pf.h.grids:
+ g.field_parameters.update(_sample_parameters)
+ conv = field._convert_function(g) or 1.0
+ v1 = g[self.field_name]
+ g.clear_data()
+ g.field_parameters.update(_sample_parameters)
+ assert_equal(v1, conv*field._function(field, g))
+
+def test_all_fields():
+ for field in FieldInfo:
+ if field.startswith("CuttingPlane"): continue
+ if field.startswith("particle"): continue
+ if field.startswith("CIC"): continue
+ if FieldInfo[field].particle_type: continue
+ for nproc in [1, 4, 8]:
+ yield TestFieldAccess(field, nproc)
diff -r 06462f53a16bbbb175a88920ed0c219aacbdc4d9 -r ae66175414388a369d488d58fab308f84e80203b yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -448,7 +448,7 @@
# This is rho_total / rho_cr(z).
def _Convert_Overdensity(data):
- return 1 / (rho_crit_now * data.pf.hubble_constant**2 *
+ return 1.0 / (rho_crit_now * data.pf.hubble_constant**2 *
(1+data.pf.current_redshift)**3)
add_field("Overdensity",function=_Matter_Density,
convert_function=_Convert_Overdensity, units=r"")
@@ -468,8 +468,8 @@
else:
omega_baryon_now = 0.0441
return data['Density'] / (omega_baryon_now * rho_crit_now *
- (data.pf['CosmologyHubbleConstantNow']**2) *
- ((1+data.pf['CosmologyCurrentRedshift'])**3))
+ (data.pf.hubble_constant**2) *
+ ((1+data.pf.current_redshift)**3))
add_field("Baryon_Overdensity", function=_Baryon_Overdensity,
units=r"")
https://bitbucket.org/yt_analysis/yt/changeset/a98f8afb7386/
changeset: a98f8afb7386
branch: yt
user: ngoldbaum
date: 2012-10-18 20:53:46
summary: Fixing field access for fields that depend on coordinate transformations.
affected #: 2 files
diff -r ae66175414388a369d488d58fab308f84e80203b -r a98f8afb73867ce458753dd52a54a8881e001c49 yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -224,9 +224,9 @@
def _sph_r(field, data):
center = data.get_field_parameter("center")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
- return get_sph_r(vectors, center)
+ return get_sph_r(coords)
def _Convert_sph_r_CGS(data):
return data.convert("cm")
@@ -241,7 +241,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
return get_sph_theta(coords, normal)
@@ -254,7 +254,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
return get_sph_phi(coords, normal)
@@ -266,7 +266,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
return get_cyl_r(coords, normal)
@@ -286,7 +286,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
return get_cyl_z(coords, normal)
@@ -303,7 +303,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
return get_cyl_theta(coords, normal)
@@ -344,9 +344,9 @@
def _cyl_RadialVelocity(field, data):
normal = data.get_field_parameter("normal")
- velocities = obtain_rv_vec(data).transpose()
+ velocities = obtain_rv_vec(data)
- theta = np.tile(data['cyl_theta'], (3, 1)).transpose()
+ theta = data['cyl_theta']
return get_cyl_r_component(velocities, theta, normal)
@@ -369,8 +369,8 @@
def _cyl_TangentialVelocity(field, data):
normal = data.get_field_parameter("normal")
- velocities = obtain_rv_vec(data).transpose()
- theta = np.tile(data['cyl_theta'], (3, 1)).transpose()
+ velocities = obtain_rv_vec(data)
+ theta = data['cyl_theta']
return get_cyl_theta_component(velocities, theta, normal)
@@ -876,9 +876,9 @@
def _RadialVelocity(field, data):
normal = data.get_field_parameter("normal")
- velocities = obtain_rv_vec(data).transpose()
- theta = np.tile(data['sph_theta'], (3, 1)).transpose()
- phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+ velocities = obtain_rv_vec(data)
+ theta = data['sph_theta']
+ phi = data['sph_phi']
return get_sph_r_component(velocities, theta, phi, normal)
@@ -1022,8 +1022,8 @@
Bfields = np.array([data['Bx'], data['By'], data['Bz']])
- theta = np.tile(data['sph_theta'], (3, 1)).transpose()
- phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+ theta = data['sph_theta']
+ phi = data['sph_phi']
return get_sph_theta_component(Bfields, theta, phi, normal)
@@ -1036,7 +1036,7 @@
Bfields = np.array([data['Bx'], data['By'], data['Bz']])
- phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+ phi = data['sph_phi']
return get_sph_phi_component(Bfields, phi, normal)
@@ -1049,8 +1049,8 @@
Bfields = np.array([data['Bx'], data['By'], data['Bz']])
- theta = np.tile(data['sph_theta'], (3, 1)).transpose()
- phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+ theta = data['sph_theta']
+ phi = data['sph_phi']
return get_sph_r_component(Bfields, theta, phi, normal)
diff -r ae66175414388a369d488d58fab308f84e80203b -r a98f8afb73867ce458753dd52a54a8881e001c49 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py
+++ b/yt/utilities/math_utils.py
@@ -686,20 +686,29 @@
# The spherical coordinates radius is simply the magnitude of the
# coordinate vector.
- return np.sqrt(np.sum(coords**2, axis=-1))
+ return np.sqrt(np.sum(coords**2, axis=0))
+def resize_vector(vector,vector_array):
+ if len(vector_array.shape) == 4:
+ res_vector = np.resize(vector,(3,1,1,1))
+ else:
+ res_vector = np.resize(vector,(3,1))
+ return res_vector
def get_sph_theta(coords, normal):
# The angle (theta) with respect to the normal (J), is the arccos
# of the dot product of the normal with the normalized coordinate
# vector.
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
+ res_normal = resize_vector(normal, coords)
- JdotCoords = np.sum(J*coords,axis=-1)
+ tile_shape = [1] + list(coords.shape)[1:]
- return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=-1)) )
+ J = np.tile(res_normal,tile_shape)
+
+ JdotCoords = np.sum(J*coords,axis=0)
+
+ return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=0)) )
def get_sph_phi(coords, normal):
# We have freedom with respect to what axis (xprime) to define
@@ -713,13 +722,16 @@
# vector.
(xprime, yprime, zprime) = get_ortho_basis(normal)
-
- tile_shape = list(coords.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
- Px = np.sum(Jx*coords,axis=-1)
- Py = np.sum(Jy*coords,axis=-1)
+ res_xprime = resize_vector(xprime, coords)
+ res_yprime = resize_vector(yprime, coords)
+
+ tile_shape = [1] + list(coords.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
+
+ Px = np.sum(Jx*coords,axis=0)
+ Py = np.sum(Jy*coords,axis=0)
return np.arctan2(Py,Px)
@@ -727,20 +739,24 @@
# The cross product of the normal (J) with a coordinate vector
# gives a vector of magnitude equal to the cylindrical radius.
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal, tile_shape)
+ res_normal = resize_vector(normal, coords)
+
+ tile_shape = [1] + list(coords.shape)[1:]
+ J = np.tile(res_normal, tile_shape)
- JcrossCoords = np.cross(J, coords)
- return np.sqrt(np.sum(JcrossCoords**2, axis=-1))
+ JcrossCoords = np.cross(J, coords, axisa=0, axisb=0, axisc=0)
+ return np.sqrt(np.sum(JcrossCoords**2, axis=0))
def get_cyl_z(coords, normal):
# The dot product of the normal (J) with the coordinate vector
# gives the cylindrical height.
+
+ res_normal = resize_vector(normal, coords)
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal, tile_shape)
+ tile_shape = [1] + list(coords.shape)[1:]
+ J = np.tile(res_normal, tile_shape)
- return np.sum(J*coords, axis=-1)
+ return np.sum(J*coords, axis=0)
def get_cyl_theta(coords, normal):
# This is identical to the spherical phi component
@@ -753,77 +769,96 @@
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
+ res_xprime = resize_vector(xprime, vectors)
+ res_yprime = resize_vector(yprime, vectors)
+
+ tile_shape = [1] + list(vectors.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
rhat = Jx*np.cos(theta) + Jy*np.sin(theta)
- return np.sum(vectors*rhat,axis=-1)
+ return np.sum(vectors*rhat,axis=0)
def get_cyl_theta_component(vectors, theta, normal):
# The theta component of a vector is the vector dotted with thetahat
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
+ res_xprime = resize_vector(xprime, vectors)
+ res_yprime = resize_vector(yprime, vectors)
+
+ tile_shape = [1] + list(vectors.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
thetahat = -Jx*np.sin(theta) + Jy*np.cos(theta)
- return np.sum(vectors*thetahat, axis=-1)
+ return np.sum(vectors*thetahat, axis=0)
def get_cyl_z_component(vectors, normal):
# The z component of a vector is the vector dotted with zhat
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- zhat = np.tile(zprime, tile_shape)
+ res_zprime = resize_vector(zprime, vectors)
- return np.sum(vectors*zhat, axis=-1)
+ tile_shape = [1] + list(vectors.shape)[1:]
+ zhat = np.tile(res_zprime, tile_shape)
+
+ return np.sum(vectors*zhat, axis=0)
def get_sph_r_component(vectors, theta, phi, normal):
# The r component of a vector is the vector dotted with rhat
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
- Jz = np.tile(zprime,tile_shape)
+ res_xprime = resize_vector(xprime, vectors)
+ res_yprime = resize_vector(yprime, vectors)
+ res_zprime = resize_vector(zprime, vectors)
+
+ tile_shape = [1] + list(vectors.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
+ Jz = np.tile(res_zprime,tile_shape)
rhat = Jx*np.sin(theta)*np.cos(phi) + \
Jy*np.sin(theta)*np.sin(phi) + \
Jz*np.cos(theta)
- return np.sum(vectors*rhat, axis=-1)
+ return np.sum(vectors*rhat, axis=0)
def get_sph_phi_component(vectors, phi, normal):
# The phi component of a vector is the vector dotted with phihat
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
+ res_xprime = resize_vector(xprime, vectors)
+ res_yprime = resize_vector(yprime, vectors)
+
+ tile_shape = [1] + list(vectors.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
phihat = -Jx*np.sin(phi) + Jy*np.cos(phi)
- return np.sum(vectors*phihat, axis=-1)
+ return np.sum(vectors*phihat, axis=0)
def get_sph_theta_component(vectors, theta, phi, normal):
# The theta component of a vector is the vector dotted with thetahat
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
- Jz = np.tile(zprime,tile_shape)
+ res_xprime = resize_vector(xprime, vectors)
+ res_yprime = resize_vector(yprime, vectors)
+ res_zprime = resize_vector(zprime, vectors)
+
+ tile_shape = [1] + list(vectors.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
+ Jz = np.tile(res_zprime,tile_shape)
thetahat = Jx*np.cos(theta)*np.cos(phi) + \
Jy*np.cos(theta)*np.sin(phi) - \
Jz*np.sin(theta)
- return np.sum(vectors*thetahat, axis=-1)
+ return np.sum(vectors*thetahat, axis=0)
https://bitbucket.org/yt_analysis/yt/changeset/ac14b13905d0/
changeset: ac14b13905d0
branch: yt
user: MatthewTurk
date: 2012-10-18 21:04:12
summary: Skipping WeakLensingConvergence.
affected #: 1 file
diff -r a98f8afb73867ce458753dd52a54a8881e001c49 -r ac14b13905d0563fd0c03d1ad2b92ab107e605d7 yt/data_objects/tests/test_fields.py
--- a/yt/data_objects/tests/test_fields.py
+++ b/yt/data_objects/tests/test_fields.py
@@ -85,6 +85,7 @@
if field.startswith("CuttingPlane"): continue
if field.startswith("particle"): continue
if field.startswith("CIC"): continue
+ if field.startswith("WeakLensingConvergence"): continue
if FieldInfo[field].particle_type: continue
for nproc in [1, 4, 8]:
yield TestFieldAccess(field, nproc)
https://bitbucket.org/yt_analysis/yt/changeset/4f48186e3644/
changeset: 4f48186e3644
branch: yt
user: ngoldbaum
date: 2012-10-18 21:18:49
summary: Updating the coordinate conversion tests.
affected #: 1 file
diff -r ac14b13905d0563fd0c03d1ad2b92ab107e605d7 -r 4f48186e3644fc5a4a4ce418272b210a277299a8 yt/utilities/tests/test_coordinate_conversions.py
--- a/yt/utilities/tests/test_coordinate_conversions.py
+++ b/yt/utilities/tests/test_coordinate_conversions.py
@@ -18,7 +18,7 @@
[ 0.73186508, -0.3109153 , 0.75728738],
[ 0.8757989 , -0.41475119, -0.57039201],
[ 0.58040762, 0.81969082, 0.46759728],
- [-0.89983356, -0.9853683 , -0.38355343]])
+ [-0.89983356, -0.9853683 , -0.38355343]]).T
def test_spherical_coordinate_conversion():
normal = [0, 0, 1]
@@ -85,44 +85,41 @@
def test_spherical_coordinate_projections():
normal = [0, 0, 1]
theta = get_sph_theta(coords, normal)
- theta_pass = np.tile(theta, (3, 1)).T
phi = get_sph_phi(coords, normal)
- phi_pass = np.tile(phi, (3, 1)).T
- zero = np.tile(0,coords.shape[0])
+ zero = np.tile(0,coords.shape[1])
# Purely radial field
- vecs = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)]).T
- assert_array_almost_equal(zero, get_sph_theta_component(vecs, theta_pass, phi_pass, normal))
- assert_array_almost_equal(zero, get_sph_phi_component(vecs, phi_pass, normal))
+ vecs = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)])
+ assert_array_almost_equal(zero, get_sph_theta_component(vecs, theta, phi, normal))
+ assert_array_almost_equal(zero, get_sph_phi_component(vecs, phi, normal))
# Purely toroidal field
- vecs = np.array([-np.sin(phi), np.cos(phi), zero]).T
- assert_array_almost_equal(zero, get_sph_theta_component(vecs, theta_pass, phi_pass, normal))
- assert_array_almost_equal(zero, get_sph_r_component(vecs, theta_pass, phi_pass, normal))
+ vecs = np.array([-np.sin(phi), np.cos(phi), zero])
+ assert_array_almost_equal(zero, get_sph_theta_component(vecs, theta, phi, normal))
+ assert_array_almost_equal(zero, get_sph_r_component(vecs, theta, phi, normal))
# Purely poloidal field
- vecs = np.array([np.cos(theta)*np.cos(phi), np.cos(theta)*np.sin(phi), -np.sin(theta)]).T
- assert_array_almost_equal(zero, get_sph_phi_component(vecs, phi_pass, normal))
- assert_array_almost_equal(zero, get_sph_r_component(vecs, theta_pass, phi_pass, normal))
+ vecs = np.array([np.cos(theta)*np.cos(phi), np.cos(theta)*np.sin(phi), -np.sin(theta)])
+ assert_array_almost_equal(zero, get_sph_phi_component(vecs, phi, normal))
+ assert_array_almost_equal(zero, get_sph_r_component(vecs, theta, phi, normal))
def test_cylindrical_coordinate_projections():
normal = [0, 0, 1]
theta = get_cyl_theta(coords, normal)
- theta_pass = np.tile(theta, (3, 1)).T
z = get_cyl_z(coords, normal)
- zero = np.tile(0, coords.shape[0])
+ zero = np.tile(0, coords.shape[1])
# Purely radial field
- vecs = np.array([np.cos(theta), np.sin(theta), zero]).T
- assert_array_almost_equal(zero, get_cyl_theta_component(vecs, theta_pass, normal))
+ vecs = np.array([np.cos(theta), np.sin(theta), zero])
+ assert_array_almost_equal(zero, get_cyl_theta_component(vecs, theta, normal))
assert_array_almost_equal(zero, get_cyl_z_component(vecs, normal))
# Purely toroidal field
- vecs = np.array([-np.sin(theta), np.cos(theta), zero]).T
+ vecs = np.array([-np.sin(theta), np.cos(theta), zero])
assert_array_almost_equal(zero, get_cyl_z_component(vecs, normal))
- assert_array_almost_equal(zero, get_cyl_r_component(vecs, theta_pass, normal))
+ assert_array_almost_equal(zero, get_cyl_r_component(vecs, theta, normal))
# Purely z field
- vecs = np.array([zero, zero, z]).T
- assert_array_almost_equal(zero, get_cyl_theta_component(vecs, theta_pass, normal))
- assert_array_almost_equal(zero, get_cyl_r_component(vecs, theta_pass, normal))
+ vecs = np.array([zero, zero, z])
+ assert_array_almost_equal(zero, get_cyl_theta_component(vecs, theta, normal))
+ assert_array_almost_equal(zero, get_cyl_r_component(vecs, theta, normal))
https://bitbucket.org/yt_analysis/yt/changeset/0e2a183fc784/
changeset: 0e2a183fc784
branch: yt
user: ngoldbaum
date: 2012-10-18 21:21:06
summary: Merged in MatthewTurk/yt (pull request #306)
affected #: 5 files
diff -r 21333cf31d213c5bac726e8b8f71f6c2f7520afa -r 0e2a183fc7848ef746d5608cfb9cc837fda95ea5 yt/data_objects/field_info_container.py
--- a/yt/data_objects/field_info_container.py
+++ b/yt/data_objects/field_info_container.py
@@ -160,7 +160,8 @@
# required attrs
pf = fake_parameter_file(lambda: 1)
pf.current_redshift = pf.omega_lambda = pf.omega_matter = \
- pf.hubble_constant = pf.cosmological_simulation = 0.0
+ pf.cosmological_simulation = 0.0
+ pf.hubble_constant = 0.7
pf.domain_left_edge = np.zeros(3, 'float64')
pf.domain_right_edge = np.ones(3, 'float64')
pf.dimensionality = 3
diff -r 21333cf31d213c5bac726e8b8f71f6c2f7520afa -r 0e2a183fc7848ef746d5608cfb9cc837fda95ea5 yt/data_objects/tests/test_fields.py
--- /dev/null
+++ b/yt/data_objects/tests/test_fields.py
@@ -0,0 +1,91 @@
+from yt.testing import *
+import numpy as np
+from yt.data_objects.field_info_container import \
+ FieldInfo
+import yt.data_objects.universal_fields
+from yt.utilities.definitions import \
+ mpc_conversion, sec_conversion
+
+def setup():
+ from yt.config import ytcfg
+ ytcfg["yt","__withintesting"] = "True"
+ np.seterr(all = 'ignore')
+
+_sample_parameters = dict(
+ axis = 0,
+ center = np.array((0.0, 0.0, 0.0)),
+ bulk_velocity = np.array((0.0, 0.0, 0.0)),
+ normal = np.array((0.0, 0.0, 1.0)),
+ cp_x_vec = np.array((1.0, 0.0, 0.0)),
+ cp_y_vec = np.array((0.0, 1.0, 0.0)),
+ cp_z_vec = np.array((0.0, 0.0, 1.0)),
+)
+
+_base_fields = ["Density", "x-velocity", "y-velocity", "z-velocity"]
+
+def realistic_pf(fields, nprocs):
+ pf = fake_random_pf(16, fields = fields, nprocs = nprocs)
+ pf.parameters["HydroMethod"] = "streaming"
+ pf.parameters["Gamma"] = 5.0/3.0
+ pf.parameters["EOSType"] = 1.0
+ pf.parameters["EOSSoundSpeed"] = 1.0
+ pf.conversion_factors["Time"] = 1.0
+ pf.conversion_factors.update( dict((f, 1.0) for f in fields) )
+ pf.current_redshift = 0.0001
+ pf.hubble_constant = 0.7
+ for unit in mpc_conversion:
+ pf.units[unit+'h'] = pf.units[unit]
+ pf.units[unit+'cm'] = pf.units[unit]
+ pf.units[unit+'hcm'] = pf.units[unit]
+ return pf
+
+class TestFieldAccess(object):
+ description = None
+
+ def __init__(self, field_name, nproc):
+ # Note this should be a field name
+ self.field_name = field_name
+ self.description = "Accessing_%s_%s" % (field_name, nproc)
+ self.nproc = nproc
+
+ def __call__(self):
+ field = FieldInfo[self.field_name]
+ deps = field.get_dependencies()
+ fields = deps.requested + _base_fields
+ skip_grids = False
+ needs_spatial = False
+ for v in field.validators:
+ f = getattr(v, "fields", None)
+ if f: fields += f
+ if getattr(v, "ghost_zones", 0) > 0:
+ skip_grids = True
+ if hasattr(v, "ghost_zones"):
+ needs_spatial = True
+ pf = realistic_pf(fields, self.nproc)
+ # This gives unequal sized grids as well as subgrids
+ dd1 = pf.h.all_data()
+ dd2 = pf.h.all_data()
+ dd1.field_parameters.update(_sample_parameters)
+ dd2.field_parameters.update(_sample_parameters)
+ v1 = dd1[self.field_name]
+ conv = field._convert_function(dd1) or 1.0
+ if not needs_spatial:
+ assert_equal(v1, conv*field._function(field, dd2))
+ if not skip_grids:
+ for g in pf.h.grids:
+ g.field_parameters.update(_sample_parameters)
+ conv = field._convert_function(g) or 1.0
+ v1 = g[self.field_name]
+ g.clear_data()
+ g.field_parameters.update(_sample_parameters)
+ assert_equal(v1, conv*field._function(field, g))
+
+def test_all_fields():
+ for field in FieldInfo:
+ if field.startswith("CuttingPlane"): continue
+ if field.startswith("particle"): continue
+ if field.startswith("CIC"): continue
+ if field.startswith("WeakLensingConvergence"): continue
+ if FieldInfo[field].particle_type: continue
+ for nproc in [1, 4, 8]:
+ yield TestFieldAccess(field, nproc)
diff -r 21333cf31d213c5bac726e8b8f71f6c2f7520afa -r 0e2a183fc7848ef746d5608cfb9cc837fda95ea5 yt/data_objects/universal_fields.py
--- a/yt/data_objects/universal_fields.py
+++ b/yt/data_objects/universal_fields.py
@@ -224,9 +224,9 @@
def _sph_r(field, data):
center = data.get_field_parameter("center")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
- return get_sph_r(vectors, center)
+ return get_sph_r(coords)
def _Convert_sph_r_CGS(data):
return data.convert("cm")
@@ -241,7 +241,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
return get_sph_theta(coords, normal)
@@ -254,7 +254,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
return get_sph_phi(coords, normal)
@@ -266,7 +266,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
return get_cyl_r(coords, normal)
@@ -286,7 +286,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
return get_cyl_z(coords, normal)
@@ -303,7 +303,7 @@
center = data.get_field_parameter("center")
normal = data.get_field_parameter("normal")
- coords = obtain_rvec(data).transpose()
+ coords = obtain_rvec(data)
return get_cyl_theta(coords, normal)
@@ -344,9 +344,9 @@
def _cyl_RadialVelocity(field, data):
normal = data.get_field_parameter("normal")
- velocities = obtain_rv_vec(data).transpose()
+ velocities = obtain_rv_vec(data)
- theta = np.tile(data['cyl_theta'], (3, 1)).transpose()
+ theta = data['cyl_theta']
return get_cyl_r_component(velocities, theta, normal)
@@ -369,8 +369,8 @@
def _cyl_TangentialVelocity(field, data):
normal = data.get_field_parameter("normal")
- velocities = obtain_rv_vec(data).transpose()
- theta = np.tile(data['cyl_theta'], (3, 1)).transpose()
+ velocities = obtain_rv_vec(data)
+ theta = data['cyl_theta']
return get_cyl_theta_component(velocities, theta, normal)
@@ -448,7 +448,7 @@
# This is rho_total / rho_cr(z).
def _Convert_Overdensity(data):
- return 1 / (rho_crit_now * data.pf.hubble_constant**2 *
+ return 1.0 / (rho_crit_now * data.pf.hubble_constant**2 *
(1+data.pf.current_redshift)**3)
add_field("Overdensity",function=_Matter_Density,
convert_function=_Convert_Overdensity, units=r"")
@@ -468,8 +468,8 @@
else:
omega_baryon_now = 0.0441
return data['Density'] / (omega_baryon_now * rho_crit_now *
- (data.pf['CosmologyHubbleConstantNow']**2) *
- ((1+data.pf['CosmologyCurrentRedshift'])**3))
+ (data.pf.hubble_constant**2) *
+ ((1+data.pf.current_redshift)**3))
add_field("Baryon_Overdensity", function=_Baryon_Overdensity,
units=r"")
@@ -876,9 +876,9 @@
def _RadialVelocity(field, data):
normal = data.get_field_parameter("normal")
- velocities = obtain_rv_vec(data).transpose()
- theta = np.tile(data['sph_theta'], (3, 1)).transpose()
- phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+ velocities = obtain_rv_vec(data)
+ theta = data['sph_theta']
+ phi = data['sph_phi']
return get_sph_r_component(velocities, theta, phi, normal)
@@ -1022,8 +1022,8 @@
Bfields = np.array([data['Bx'], data['By'], data['Bz']])
- theta = np.tile(data['sph_theta'], (3, 1)).transpose()
- phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+ theta = data['sph_theta']
+ phi = data['sph_phi']
return get_sph_theta_component(Bfields, theta, phi, normal)
@@ -1036,7 +1036,7 @@
Bfields = np.array([data['Bx'], data['By'], data['Bz']])
- phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+ phi = data['sph_phi']
return get_sph_phi_component(Bfields, phi, normal)
@@ -1049,8 +1049,8 @@
Bfields = np.array([data['Bx'], data['By'], data['Bz']])
- theta = np.tile(data['sph_theta'], (3, 1)).transpose()
- phi = np.tile(data['sph_phi'], (3, 1)).transpose()
+ theta = data['sph_theta']
+ phi = data['sph_phi']
return get_sph_r_component(Bfields, theta, phi, normal)
diff -r 21333cf31d213c5bac726e8b8f71f6c2f7520afa -r 0e2a183fc7848ef746d5608cfb9cc837fda95ea5 yt/utilities/math_utils.py
--- a/yt/utilities/math_utils.py
+++ b/yt/utilities/math_utils.py
@@ -686,20 +686,29 @@
# The spherical coordinates radius is simply the magnitude of the
# coordinate vector.
- return np.sqrt(np.sum(coords**2, axis=-1))
+ return np.sqrt(np.sum(coords**2, axis=0))
+def resize_vector(vector,vector_array):
+ if len(vector_array.shape) == 4:
+ res_vector = np.resize(vector,(3,1,1,1))
+ else:
+ res_vector = np.resize(vector,(3,1))
+ return res_vector
def get_sph_theta(coords, normal):
# The angle (theta) with respect to the normal (J), is the arccos
# of the dot product of the normal with the normalized coordinate
# vector.
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal,tile_shape)
+ res_normal = resize_vector(normal, coords)
- JdotCoords = np.sum(J*coords,axis=-1)
+ tile_shape = [1] + list(coords.shape)[1:]
- return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=-1)) )
+ J = np.tile(res_normal,tile_shape)
+
+ JdotCoords = np.sum(J*coords,axis=0)
+
+ return np.arccos( JdotCoords / np.sqrt(np.sum(coords**2,axis=0)) )
def get_sph_phi(coords, normal):
# We have freedom with respect to what axis (xprime) to define
@@ -713,13 +722,16 @@
# vector.
(xprime, yprime, zprime) = get_ortho_basis(normal)
-
- tile_shape = list(coords.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
- Px = np.sum(Jx*coords,axis=-1)
- Py = np.sum(Jy*coords,axis=-1)
+ res_xprime = resize_vector(xprime, coords)
+ res_yprime = resize_vector(yprime, coords)
+
+ tile_shape = [1] + list(coords.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
+
+ Px = np.sum(Jx*coords,axis=0)
+ Py = np.sum(Jy*coords,axis=0)
return np.arctan2(Py,Px)
@@ -727,20 +739,24 @@
# The cross product of the normal (J) with a coordinate vector
# gives a vector of magnitude equal to the cylindrical radius.
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal, tile_shape)
+ res_normal = resize_vector(normal, coords)
+
+ tile_shape = [1] + list(coords.shape)[1:]
+ J = np.tile(res_normal, tile_shape)
- JcrossCoords = np.cross(J, coords)
- return np.sqrt(np.sum(JcrossCoords**2, axis=-1))
+ JcrossCoords = np.cross(J, coords, axisa=0, axisb=0, axisc=0)
+ return np.sqrt(np.sum(JcrossCoords**2, axis=0))
def get_cyl_z(coords, normal):
# The dot product of the normal (J) with the coordinate vector
# gives the cylindrical height.
+
+ res_normal = resize_vector(normal, coords)
- tile_shape = list(coords.shape)[:-1] + [1]
- J = np.tile(normal, tile_shape)
+ tile_shape = [1] + list(coords.shape)[1:]
+ J = np.tile(res_normal, tile_shape)
- return np.sum(J*coords, axis=-1)
+ return np.sum(J*coords, axis=0)
def get_cyl_theta(coords, normal):
# This is identical to the spherical phi component
@@ -753,77 +769,96 @@
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
+ res_xprime = resize_vector(xprime, vectors)
+ res_yprime = resize_vector(yprime, vectors)
+
+ tile_shape = [1] + list(vectors.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
rhat = Jx*np.cos(theta) + Jy*np.sin(theta)
- return np.sum(vectors*rhat,axis=-1)
+ return np.sum(vectors*rhat,axis=0)
def get_cyl_theta_component(vectors, theta, normal):
# The theta component of a vector is the vector dotted with thetahat
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
+ res_xprime = resize_vector(xprime, vectors)
+ res_yprime = resize_vector(yprime, vectors)
+
+ tile_shape = [1] + list(vectors.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
thetahat = -Jx*np.sin(theta) + Jy*np.cos(theta)
- return np.sum(vectors*thetahat, axis=-1)
+ return np.sum(vectors*thetahat, axis=0)
def get_cyl_z_component(vectors, normal):
# The z component of a vector is the vector dotted with zhat
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- zhat = np.tile(zprime, tile_shape)
+ res_zprime = resize_vector(zprime, vectors)
- return np.sum(vectors*zhat, axis=-1)
+ tile_shape = [1] + list(vectors.shape)[1:]
+ zhat = np.tile(res_zprime, tile_shape)
+
+ return np.sum(vectors*zhat, axis=0)
def get_sph_r_component(vectors, theta, phi, normal):
# The r component of a vector is the vector dotted with rhat
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
- Jz = np.tile(zprime,tile_shape)
+ res_xprime = resize_vector(xprime, vectors)
+ res_yprime = resize_vector(yprime, vectors)
+ res_zprime = resize_vector(zprime, vectors)
+
+ tile_shape = [1] + list(vectors.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
+ Jz = np.tile(res_zprime,tile_shape)
rhat = Jx*np.sin(theta)*np.cos(phi) + \
Jy*np.sin(theta)*np.sin(phi) + \
Jz*np.cos(theta)
- return np.sum(vectors*rhat, axis=-1)
+ return np.sum(vectors*rhat, axis=0)
def get_sph_phi_component(vectors, phi, normal):
# The phi component of a vector is the vector dotted with phihat
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
+ res_xprime = resize_vector(xprime, vectors)
+ res_yprime = resize_vector(yprime, vectors)
+
+ tile_shape = [1] + list(vectors.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
phihat = -Jx*np.sin(phi) + Jy*np.cos(phi)
- return np.sum(vectors*phihat, axis=-1)
+ return np.sum(vectors*phihat, axis=0)
def get_sph_theta_component(vectors, theta, phi, normal):
# The theta component of a vector is the vector dotted with thetahat
(xprime, yprime, zprime) = get_ortho_basis(normal)
- tile_shape = list(vectors.shape)[:-1] + [1]
- Jx = np.tile(xprime,tile_shape)
- Jy = np.tile(yprime,tile_shape)
- Jz = np.tile(zprime,tile_shape)
+ res_xprime = resize_vector(xprime, vectors)
+ res_yprime = resize_vector(yprime, vectors)
+ res_zprime = resize_vector(zprime, vectors)
+
+ tile_shape = [1] + list(vectors.shape)[1:]
+ Jx = np.tile(res_xprime,tile_shape)
+ Jy = np.tile(res_yprime,tile_shape)
+ Jz = np.tile(res_zprime,tile_shape)
thetahat = Jx*np.cos(theta)*np.cos(phi) + \
Jy*np.cos(theta)*np.sin(phi) - \
Jz*np.sin(theta)
- return np.sum(vectors*thetahat, axis=-1)
+ return np.sum(vectors*thetahat, axis=0)
diff -r 21333cf31d213c5bac726e8b8f71f6c2f7520afa -r 0e2a183fc7848ef746d5608cfb9cc837fda95ea5 yt/utilities/tests/test_coordinate_conversions.py
--- a/yt/utilities/tests/test_coordinate_conversions.py
+++ b/yt/utilities/tests/test_coordinate_conversions.py
@@ -18,7 +18,7 @@
[ 0.73186508, -0.3109153 , 0.75728738],
[ 0.8757989 , -0.41475119, -0.57039201],
[ 0.58040762, 0.81969082, 0.46759728],
- [-0.89983356, -0.9853683 , -0.38355343]])
+ [-0.89983356, -0.9853683 , -0.38355343]]).T
def test_spherical_coordinate_conversion():
normal = [0, 0, 1]
@@ -85,44 +85,41 @@
def test_spherical_coordinate_projections():
normal = [0, 0, 1]
theta = get_sph_theta(coords, normal)
- theta_pass = np.tile(theta, (3, 1)).T
phi = get_sph_phi(coords, normal)
- phi_pass = np.tile(phi, (3, 1)).T
- zero = np.tile(0,coords.shape[0])
+ zero = np.tile(0,coords.shape[1])
# Purely radial field
- vecs = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)]).T
- assert_array_almost_equal(zero, get_sph_theta_component(vecs, theta_pass, phi_pass, normal))
- assert_array_almost_equal(zero, get_sph_phi_component(vecs, phi_pass, normal))
+ vecs = np.array([np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)])
+ assert_array_almost_equal(zero, get_sph_theta_component(vecs, theta, phi, normal))
+ assert_array_almost_equal(zero, get_sph_phi_component(vecs, phi, normal))
# Purely toroidal field
- vecs = np.array([-np.sin(phi), np.cos(phi), zero]).T
- assert_array_almost_equal(zero, get_sph_theta_component(vecs, theta_pass, phi_pass, normal))
- assert_array_almost_equal(zero, get_sph_r_component(vecs, theta_pass, phi_pass, normal))
+ vecs = np.array([-np.sin(phi), np.cos(phi), zero])
+ assert_array_almost_equal(zero, get_sph_theta_component(vecs, theta, phi, normal))
+ assert_array_almost_equal(zero, get_sph_r_component(vecs, theta, phi, normal))
# Purely poloidal field
- vecs = np.array([np.cos(theta)*np.cos(phi), np.cos(theta)*np.sin(phi), -np.sin(theta)]).T
- assert_array_almost_equal(zero, get_sph_phi_component(vecs, phi_pass, normal))
- assert_array_almost_equal(zero, get_sph_r_component(vecs, theta_pass, phi_pass, normal))
+ vecs = np.array([np.cos(theta)*np.cos(phi), np.cos(theta)*np.sin(phi), -np.sin(theta)])
+ assert_array_almost_equal(zero, get_sph_phi_component(vecs, phi, normal))
+ assert_array_almost_equal(zero, get_sph_r_component(vecs, theta, phi, normal))
def test_cylindrical_coordinate_projections():
normal = [0, 0, 1]
theta = get_cyl_theta(coords, normal)
- theta_pass = np.tile(theta, (3, 1)).T
z = get_cyl_z(coords, normal)
- zero = np.tile(0, coords.shape[0])
+ zero = np.tile(0, coords.shape[1])
# Purely radial field
- vecs = np.array([np.cos(theta), np.sin(theta), zero]).T
- assert_array_almost_equal(zero, get_cyl_theta_component(vecs, theta_pass, normal))
+ vecs = np.array([np.cos(theta), np.sin(theta), zero])
+ assert_array_almost_equal(zero, get_cyl_theta_component(vecs, theta, normal))
assert_array_almost_equal(zero, get_cyl_z_component(vecs, normal))
# Purely toroidal field
- vecs = np.array([-np.sin(theta), np.cos(theta), zero]).T
+ vecs = np.array([-np.sin(theta), np.cos(theta), zero])
assert_array_almost_equal(zero, get_cyl_z_component(vecs, normal))
- assert_array_almost_equal(zero, get_cyl_r_component(vecs, theta_pass, normal))
+ assert_array_almost_equal(zero, get_cyl_r_component(vecs, theta, normal))
# Purely z field
- vecs = np.array([zero, zero, z]).T
- assert_array_almost_equal(zero, get_cyl_theta_component(vecs, theta_pass, normal))
- assert_array_almost_equal(zero, get_cyl_r_component(vecs, theta_pass, normal))
+ vecs = np.array([zero, zero, z])
+ assert_array_almost_equal(zero, get_cyl_theta_component(vecs, theta, normal))
+ assert_array_almost_equal(zero, get_cyl_r_component(vecs, theta, normal))
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