[Yt-svn] yt-commit r477 - in trunk/yt: . lagos reason
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Fri May 16 00:17:51 PDT 2008
Author: mturk
Date: Fri May 16 00:17:50 2008
New Revision: 477
URL: http://yt.spacepope.org/changeset/477
Log:
* The FieldDetector in the DerivedFields now works. Moving toward closing #86.
* 2DProfiles had some issues working with objects with data outside the
boundaries. Now fixed.
* just_one is more intelligent now for iterables versus arrays.
* VelocityMagnitude now accepts a bulk_velocity parameter
* Added some fields -- DiskAngle, some Tangential fields that I might remove.
* NumberDensity now uses a default mu of 0.6 if mu is not provided.
* PointCombine has some =NULL initializers to retain working behavior in case
of error.
* Notebook redraws the plot as well as the canvase on calls to UpdateCanvas.
Modified:
trunk/yt/funcs.py
trunk/yt/lagos/DerivedFields.py
trunk/yt/lagos/PointCombine.c
trunk/yt/lagos/Profiles.py
trunk/yt/reason/Notebook.py
Modified: trunk/yt/funcs.py
==============================================================================
--- trunk/yt/funcs.py (original)
+++ trunk/yt/funcs.py Fri May 16 00:17:50 2008
@@ -79,7 +79,9 @@
self._pbar.Destroy()
def just_one(obj):
- if iterable(obj):
+ if hasattr(obj,'flat'):
+ return obj.flat[0]
+ elif iterable(obj):
return obj[0]
return obj
Modified: trunk/yt/lagos/DerivedFields.py
==============================================================================
--- trunk/yt/lagos/DerivedFields.py (original)
+++ trunk/yt/lagos/DerivedFields.py Fri May 16 00:17:50 2008
@@ -79,18 +79,50 @@
name, function, **kwargs)
class FieldDetector(defaultdict):
- pf = defaultdict(lambda: 1)
- def __init__(self):
+ Level = 1
+ NumberOfParticles = 0
+ def __init__(self, nd = 16, pf = None):
+ self.nd = nd
+ self.ActiveDimensions = [nd,nd,nd]
+ self.LeftEdge = [0.0,0.0,0.0]
+ self.RightEdge = [1.0,1.0,1.0]
+ self.dx = self.dy = self.dz = na.array([1.0])
+ self.fields = []
+ if pf is None:
+ pf = defaultdict(lambda: 1)
+ self.pf = pf
self.requested = []
- defaultdict.__init__(self, lambda: na.ones(10))
+ self.requested_parameters = []
+ defaultdict.__init__(self, lambda: na.ones((nd,nd,nd)))
def __missing__(self, item):
- if fieldInfo.has_key(item):
- fieldInfo[item](self)
- return self[item]
+ if fieldInfo.has_key(item) and \
+ fieldInfo[item]._function.func_name != '<lambda>':
+ try:
+ vv = fieldInfo[item](self)
+ except NeedsGridType, exc:
+ ngz = exc.ghost_zones
+ nfd = FieldDetector(self.nd+ngz*2)
+ vv = fieldInfo[item](nfd)[ngz:-ngz,ngz:-ngz,ngz:-ngz]
+ for i in vv.requested:
+ if i not in self.requested: self.requested.append(i)
+ for i in vv.requested_parameters:
+ if i not in self.requested_parameters: self.requested_parameters.append(i)
+ if vv is not None:
+ self[item] = vv
+ return self[item]
self.requested.append(item)
return defaultdict.__missing__(self, item)
- def convert(self, item):
- return 1
+ def get_field_parameter(self, param):
+ self.requested_parameters.append(param)
+ if param in ['bulk_velocity','center','height_vector']:
+ return na.array([0,0,0])
+ else:
+ return 0.0
+ _spatial = True
+ _num_ghost_zones = 0
+ id = 1
+ def has_field_parameter(self, param): return True
+ def convert(self, item): return 1
class DerivedField:
def __init__(self, name, function,
@@ -121,10 +153,13 @@
validator(data)
# If we don't get an exception, we're good to go
return True
- def get_dependencies(self):
- e = FieldDetector()
- self(e)
- return e.requested
+ def get_dependencies(self, *args, **kwargs):
+ e = FieldDetector(*args, **kwargs)
+ if self._function.func_name == '<lambda>':
+ e.requested.append(self.name)
+ else:
+ self(e)
+ return e
def get_units(self):
return self._units
def get_projected_units(self):
@@ -151,7 +186,7 @@
def __call__(self, data):
doesnt_have = []
for p in self.parameters:
- if not data.field_parameters.has_key(p):
+ if not data.has_field_parameter(p):
doesnt_have.append(p)
if len(doesnt_have) > 0:
raise NeedsParameter(doesnt_have)
@@ -163,6 +198,7 @@
self.fields = ensure_list(field)
def __call__(self, data):
doesnt_have = []
+ if isinstance(data, FieldDetector): return True
for f in self.fields:
if f not in data.hierarchy.field_list:
doesnt_have.append(f)
@@ -191,6 +227,7 @@
def __call__(self, data):
# When we say spatial information, we really mean
# that it has a three-dimensional data structure
+ if isinstance(data, FieldDetector): return True
if not data._spatial:
raise NeedsGridType(self.ghost_zones,self.fields)
if self.ghost_zones == data._num_ghost_zones:
@@ -344,11 +381,23 @@
def _VelocityMagnitude(field, data):
"""M{|v|}"""
- return ( data["x-velocity"]**2.0 + \
- data["y-velocity"]**2.0 + \
- data["z-velocity"]**2.0 )**(1.0/2.0)
+ bulk_velocity = data.get_field_parameter("bulk_velocity")
+ if bulk_velocity == None:
+ bulk_velocity = na.zeros(3)
+ return ( (data["x-velocity"]-bulk_velocity[0])**2.0 + \
+ (data["y-velocity"]-bulk_velocity[1])**2.0 + \
+ (data["z-velocity"]-bulk_velocity[2])**2.0 )**(1.0/2.0)
add_field("VelocityMagnitude", take_log=False, units=r"\rm{cm}/\rm{s}")
+def _TangentialOverVelocityMagnitude(field, data):
+ return na.abs(data["TangentialVelocity"])/na.abs(data["VelocityMagnitude"])
+add_field("TangentialOverVelocityMagnitude", take_log=False)
+
+def _TangentialVelocity(field, data):
+ return na.sqrt(data["VelocityMagnitude"]**2.0
+ - data["RadialVelocity"]**2.0)
+add_field("TangentialVelocity", take_log=False, units=r"\rm{cm}/\rm{s}")
+
def _Pressure(field, data):
"""M{(Gamma-1.0)*rho*E}"""
return (data.pf["Gamma"] - 1.0) * \
@@ -358,7 +407,7 @@
def _ThermalEnergy(field, data):
if data.pf["HydroMethod"] == 2:
return data["Total_Energy"]
- if data.pf["HydroMethod"] == 0:
+ if data.pf["HydroMethod"] in [0,1]:
if data.pf["DualEnergyFormalism"]:
return data["Gas_Energy"]
else:
@@ -389,8 +438,31 @@
return na.abs(height)
def _convertHeight(data):
return data.convert("cm")
+def _convertHeightAU(data):
+ return data.convert("au")
add_field("Height", convert_function=_convertHeight,
- validators=[ValidateParameter("height_vector")])
+ validators=[ValidateParameter("height_vector")],
+ units=r"cm")
+add_field("HeightAU", function=_Height,
+ convert_function=_convertHeightAU,
+ validators=[ValidateParameter("height_vector")],
+ units=r"AU")
+
+def _DiskAngle(field, data):
+ # We make both r_vec and h_vec into unit vectors
+ center = data.get_field_parameter("center")
+ r_vec = na.array([data["x"] - center[0],
+ data["y"] - center[1],
+ data["z"] - center[2]])
+ r_vec = r_vec/na.sqrt((r_vec**2.0).sum(axis=0))
+ h_vec = na.array(data.get_field_parameter("height_vector"))
+ dp = r_vec[0,:] * h_vec[0] \
+ + r_vec[1,:] * h_vec[1] \
+ + r_vec[2,:] * h_vec[2]
+ return na.arccos(dp)
+add_field("DiskAngle", take_log=False,
+ validators=[ValidateParameter("height_vector"),
+ ValidateParameter("center")])
def _DynamicalTime(field, data):
"""
@@ -414,7 +486,11 @@
fieldData = na.zeros(data["Density"].shape,
dtype = data["Density"].dtype)
if data.pf["MultiSpecies"] == 0:
- fieldData += data["Density"] * data.get_field_parameter("mu", 0.6)
+ if data.has_field_parameter("mu"):
+ mu = data.get_field_parameter("mu")
+ else:
+ mu = 0.6
+ fieldData += data["Density"] * mu
if data.pf["MultiSpecies"] > 0:
fieldData += data["HI_Density"] / 1.0
fieldData += data["HII_Density"] / 1.0
@@ -531,7 +607,7 @@
def _DivV(field, data):
# We need to set up stencils
- if data.pf["HydroMethod"] == 0:
+ if data.pf["HydroMethod"] in [0,1]:
sl_left = slice(None,-2,None)
sl_right = slice(2,None,None)
div_fac = 2.0
@@ -541,13 +617,13 @@
div_fac = 1.0
div_x = (data["x-velocity"][sl_right,1:-1,1:-1] -
data["x-velocity"][sl_left,1:-1,1:-1]) \
- / (div_fac*data["dx"][1:-1,1:-1,1:-1])
+ / (div_fac*data["dx"].flat[0])
div_y = (data["y-velocity"][1:-1,sl_right,1:-1] -
data["y-velocity"][1:-1,sl_left,1:-1]) \
- / (div_fac*data["dy"][1:-1,1:-1,1:-1])
+ / (div_fac*data["dy"].flat[0])
div_z = (data["z-velocity"][1:-1,1:-1,sl_right] -
data["z-velocity"][1:-1,1:-1,sl_left]) \
- / (div_fac*data["dz"][1:-1,1:-1,1:-1])
+ / (div_fac*data["dz"].flat[0])
new_field = na.zeros(data["x-velocity"].shape)
new_field[1:-1,1:-1,1:-1] = div_x+div_y+div_z
return na.abs(new_field)
@@ -675,7 +751,7 @@
MJ_constant = (((5*kboltz)/(G*mh))**(1.5)) * \
(3/(4*3.1415926535897931))**(0.5) / 1.989e33
- return (MJ_constant *
+ return (MJ_constant *
((data["Temperature"]/data["MeanMolecularWeight"])**(1.5)) *
(data["Density"]**(-0.5)))
add_field("JeansMassMsun",function=_JeansMassMsun,
Modified: trunk/yt/lagos/PointCombine.c
==============================================================================
--- trunk/yt/lagos/PointCombine.c (original)
+++ trunk/yt/lagos/PointCombine.c Fri May 16 00:17:50 2008
@@ -273,6 +273,7 @@
int i, j;
PyObject *obins_x, *obins_y, *owsource, *obsource, *owresult, *obresult, *oused;
PyArrayObject *bins_x, *bins_y, *wsource, *bsource, *wresult, *bresult, *used;
+ bins_x = bins_y = wsource = bsource = wresult = bresult = used = NULL;
if (!PyArg_ParseTuple(args, "OOOOOOO",
&obins_x, &obins_y, &owsource, &obsource,
Modified: trunk/yt/lagos/Profiles.py
==============================================================================
--- trunk/yt/lagos/Profiles.py (original)
+++ trunk/yt/lagos/Profiles.py Fri May 16 00:17:50 2008
@@ -127,14 +127,14 @@
@preserve_source_parameters
def _bin_field(self, source, field, weight, accumulation,
args, check_cut=False):
- inv_bin_indices = args
+ mi, inv_bin_indices = args
if check_cut:
cm = self._data_source._get_point_indices(source)
- source_data = source[field][cm].astype('float64')
- if weight: weight_data = source[weight][cm].astype('float64')
+ source_data = source[field][cm].astype('float64')[mi]
+ if weight: weight_data = source[weight][cm].astype('float64')[mi]
else:
- source_data = source[field].astype('float64')
- if weight: weight_data = source[weight].astype('float64')
+ source_data = source[field].astype('float64')[mi]
+ if weight: weight_data = source[weight].astype('float64')[mi]
binned_field = self._get_empty_field()
weight_field = self._get_empty_field()
used_field = na.ones(weight_field.shape, dtype='bool')
@@ -167,7 +167,7 @@
inv_bin_indices = {}
for bin in range(self[self.bin_field].size):
inv_bin_indices[bin] = na.where(bin_indices == bin)
- return inv_bin_indices
+ return (mi, inv_bin_indices)
class BinnedProfile2D(BinnedProfile):
@@ -196,6 +196,11 @@
else:
self[y_bin_field] = na.linspace(
y_lower_bound*0.99, y_upper_bound*1.01, y_n_bins)
+ if na.any(na.isnan(self[x_bin_field])) \
+ or na.any(na.isnan(self[y_bin_field])):
+ mylog.error("Your min/max values for x, y have given me a nan.")
+ mylog.error("Usually this means you are asking for log, with a zero bound.")
+ raise ValueError
if not lazy_reader:
self._args = self._get_bins(data_source)
@@ -219,8 +224,11 @@
binned_field = self._get_empty_field()
weight_field = self._get_empty_field()
used_field = self._get_empty_field()
- bin_indices_x = args[0].ravel().astype('int64')
- bin_indices_y = args[1].ravel().astype('int64')
+ mi = args[0]
+ bin_indices_x = args[1].ravel().astype('int64')
+ bin_indices_y = args[2].ravel().astype('int64')
+ source_data = source_data[mi]
+ weight_data = weight_data[mi]
self.total_cells += bin_indices_x.size
nx = bin_indices_x.size
mylog.debug("Binning %s / %s times", source_data.size, nx)
@@ -257,4 +265,4 @@
bin_indices_x = na.digitize(sd_x, self[self.x_bin_field])
bin_indices_y = na.digitize(sd_y, self[self.y_bin_field])
# Now we set up our inverse bin indices
- return (bin_indices_x, bin_indices_y)
+ return (mi, bin_indices_x, bin_indices_y)
Modified: trunk/yt/reason/Notebook.py
==============================================================================
--- trunk/yt/reason/Notebook.py (original)
+++ trunk/yt/reason/Notebook.py Fri May 16 00:17:50 2008
@@ -610,6 +610,7 @@
def UpdateCanvas(self, *args):
if self.IsShown():
+ self.plot._redraw_image()
self.figure_canvas.draw()
#else: print "Opting not to update canvas"
More information about the yt-svn
mailing list