[Yt-svn] yt-commit r466 - trunk/yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Tue May 13 11:09:58 PDT 2008
Author: mturk
Date: Tue May 13 11:09:57 2008
New Revision: 466
URL: http://yt.spacepope.org/changeset/466
Log:
Added a new function to do the binding energy in C. Right now it's optional,
with the use_c kwarg.
Modified:
trunk/yt/lagos/DerivedQuantities.py
trunk/yt/lagos/PointCombine.c
Modified: trunk/yt/lagos/DerivedQuantities.py
==============================================================================
--- trunk/yt/lagos/DerivedQuantities.py (original)
+++ trunk/yt/lagos/DerivedQuantities.py Tue May 13 11:09:57 2008
@@ -155,7 +155,7 @@
add_quantity("BaryonSpinParameter", function=_BaryonSpinParameter,
combine_function=_combBaryonSpinParameter, n_ret=4)
-def _IsBound(data, truncate = True, include_thermal_energy = False):
+def _IsBound(data, truncate = True, include_thermal_energy = False, use_c = False):
# Kinetic energy
bv_x,bv_y,bv_z = data.quantities["BulkVelocity"]()
kinetic = 0.5 * (data["CellMass"] * (
@@ -175,19 +175,30 @@
cells_done = 0
potential = 0.0
pb = get_pbar("Calculating gravitational potential ", (0.5*(total_cells**2 - total_cells)))
- for q in xrange(total_cells-1):
- pot = data['CellMass'][(q+1):total_cells] / \
- (((data['x'][(q+1):total_cells]-data['x'][q])**2 + \
- (data['y'][(q+1):total_cells]-data['y'][q])**2 + \
- (data['z'][(q+1):total_cells]-data['z'][q])**2)**(0.5))
- potential += 2 * G * data['CellMass'][q] * pot.sum()
- cells_done += (total_cells - q - 1)
- pb.update(cells_done)
- if truncate and (potential > kinetic):
- pb.finish()
- return [1.0]
+ import time
+ t1 = time.time()
+ if use_c:
+ pot = PointCombine.FindBindingEnergy(data["CellMass"],
+ data["x"], data["y"], data["z"],
+ truncate, kinetic/(2*G))
+ potential += 2*G*pot
+ else:
+ t1 = time.time()
+ for q in xrange(total_cells-1):
+ pot = data['CellMass'][(q+1):total_cells] / \
+ (((data['x'][(q+1):total_cells]-data['x'][q])**2 + \
+ (data['y'][(q+1):total_cells]-data['y'][q])**2 + \
+ (data['z'][(q+1):total_cells]-data['z'][q])**2)**(0.5))
+ potential += 2 * G * data['CellMass'][q] * pot.sum()
+ cells_done += (total_cells - q - 1)
+ pb.update(cells_done)
+ if truncate and (potential > kinetic):
+ pb.finish()
+ break
pb.finish()
- return [(kinetic / potential)]
+ print "Took %0.9e" % (time.time() - t1)
+ #if truncate: return [potential > kinetic]
+ return [(potential / kinetic)]
def _combIsBound(data,bound):
return bound
add_quantity("IsBound",function=_IsBound,combine_function=_combIsBound,n_ret=1,
Modified: trunk/yt/lagos/PointCombine.c
==============================================================================
--- trunk/yt/lagos/PointCombine.c (original)
+++ trunk/yt/lagos/PointCombine.c Tue May 13 11:09:57 2008
@@ -789,6 +789,105 @@
return NULL;
}
+static PyObject *_findBindingEnergyError;
+
+static PyObject *
+Py_FindBindingEnergy(PyObject *obj, PyObject *args)
+{
+ PyObject *omass, *ox, *oy, *oz;
+ PyArrayObject *mass, *x, *y, *z;
+ int truncate;
+ double kinetic_energy;
+
+ if (!PyArg_ParseTuple(args, "OOOOid",
+ &omass, &ox, &oy, &oz, &truncate, &kinetic_energy))
+ return PyErr_Format(_findBindingEnergyError,
+ "FindBindingEnergy: Invalid parameters.");
+
+ mass = (PyArrayObject *) PyArray_FromAny(omass,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((mass==NULL) || (mass->nd != 1)) {
+ PyErr_Format(_findBindingEnergyError,
+ "FindBindingEnergy: One dimension required for mass.");
+ goto _fail;
+ }
+
+ x = (PyArrayObject *) PyArray_FromAny(ox,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((x==NULL) || (x->nd != 1)
+ || (PyArray_SIZE(x) != PyArray_SIZE(mass))) {
+ PyErr_Format(_findBindingEnergyError,
+ "FindBindingEnergy: x must be same size as mass.");
+ goto _fail;
+ }
+
+ y = (PyArrayObject *) PyArray_FromAny(oy,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((y==NULL) || (y->nd != 1)
+ || (PyArray_SIZE(y) != PyArray_SIZE(mass))) {
+ PyErr_Format(_findBindingEnergyError,
+ "FindBindingEnergy: y must be same size as mass.");
+ goto _fail;
+ }
+
+ z = (PyArrayObject *) PyArray_FromAny(oz,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((z==NULL) || (z->nd != 1)
+ || (PyArray_SIZE(z) != PyArray_SIZE(mass))) {
+ PyErr_Format(_findBindingEnergyError,
+ "FindBindingEnergy: z must be same size as mass.");
+ goto _fail;
+ }
+
+ /* Do the work here. */
+ int q_outer, q_inner;
+ double this_potential, total_potential;
+ total_potential = 0;
+ npy_float64 mass_o, x_o, y_o, z_o;
+ npy_float64 mass_i, x_i, y_i, z_i;
+
+ for (q_outer = 0; q_outer < PyArray_SIZE(mass) ; q_outer++) {
+ this_potential = 0;
+ mass_o = *(npy_float64*) PyArray_GETPTR1(mass, q_outer);
+ x_o = *(npy_float64*) PyArray_GETPTR1(x, q_outer);
+ y_o = *(npy_float64*) PyArray_GETPTR1(y, q_outer);
+ z_o = *(npy_float64*) PyArray_GETPTR1(z, q_outer);
+ for (q_inner = q_outer+1; q_inner < PyArray_SIZE(mass); q_inner++) {
+ mass_i = *(npy_float64*) PyArray_GETPTR1(mass, q_inner);
+ x_i = *(npy_float64*) PyArray_GETPTR1(x, q_inner);
+ y_i = *(npy_float64*) PyArray_GETPTR1(y, q_inner);
+ z_i = *(npy_float64*) PyArray_GETPTR1(z, q_inner);
+ this_potential += mass_o * mass_i /
+ sqrtl( (x_i-x_o)*(x_i-x_o)
+ + (y_i-y_o)*(y_i-y_o)
+ + (z_i-z_o)*(z_i-z_o) );
+ }
+ total_potential += this_potential;
+ if ((truncate == 1) && (total_potential > kinetic_energy)){
+ break;
+ }
+ }
+
+ Py_DECREF(mass);
+ Py_DECREF(x);
+ Py_DECREF(y);
+ Py_DECREF(z);
+
+ PyObject *status = PyFloat_FromDouble(total_potential);
+ return status;
+
+ _fail:
+ Py_XDECREF(mass);
+ Py_XDECREF(x);
+ Py_XDECREF(y);
+ Py_XDECREF(z);
+ return NULL;
+}
+
static PyMethodDef _combineMethods[] = {
{"CombineGrids", Py_CombineGrids, METH_VARARGS},
{"Interpolate", Py_Interpolate, METH_VARARGS},
@@ -796,6 +895,7 @@
{"DataCubeReplace", Py_DataCubeReplace, METH_VARARGS},
{"Bin2DProfile", Py_Bin2DProfile, METH_VARARGS},
{"FindContours", Py_FindContours, METH_VARARGS},
+ {"FindBindingEnergy", Py_FindBindingEnergy, METH_VARARGS},
{NULL, NULL} /* Sentinel */
};
More information about the yt-svn
mailing list