[Yt-svn] yt-commit r1293 - trunk/yt/lagos/hop
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Tue May 5 22:58:18 PDT 2009
Author: mturk
Date: Tue May 5 22:58:17 2009
New Revision: 1293
URL: http://yt.spacepope.org/changeset/1293
Log:
HOP now works in place rather than via copies.
Modified:
trunk/yt/lagos/hop/EnzoHop.c
trunk/yt/lagos/hop/hop_hop.c
trunk/yt/lagos/hop/hop_kd.c
trunk/yt/lagos/hop/hop_smooth.c
trunk/yt/lagos/hop/kd.h
Modified: trunk/yt/lagos/hop/EnzoHop.c
==============================================================================
--- trunk/yt/lagos/hop/EnzoHop.c (original)
+++ trunk/yt/lagos/hop/EnzoHop.c Tue May 5 22:58:17 2009
@@ -30,6 +30,7 @@
#include <ctype.h>
#include "kd.h"
#include "hop.h"
+#include "hop_numpy.h"
#include "numpy/ndarrayobject.h"
@@ -114,14 +115,18 @@
}
/* Copy positions into kd structure. */
+ PyArrayObject *particle_density = (PyArrayObject *)
+ PyArray_SimpleNewFromDescr(1, PyArray_DIMS(xpos),
+ PyArray_DescrFromType(NPY_FLOAT64));
- fprintf(stdout, "Filling in %d particles\n", num_particles);
- for (i = 0; i < num_particles; i++) {
- kd->p[i].r[0] = (float)(*(npy_float64*) PyArray_GETPTR1(xpos, i));
- kd->p[i].r[1] = (float)(*(npy_float64*) PyArray_GETPTR1(ypos, i));
- kd->p[i].r[2] = (float)(*(npy_float64*) PyArray_GETPTR1(zpos, i));
- kd->p[i].fMass = (float)(*(npy_float64*) PyArray_GETPTR1(mass, i)/totalmass);
- }
+ fprintf(stdout, "Copying arrays for %d particles\n", num_particles);
+ kd->np_masses = mass;
+ kd->np_pos[0] = xpos;
+ kd->np_pos[1] = ypos;
+ kd->np_pos[2] = zpos;
+ kd->np_densities = particle_density;
+ kd->totalmass = totalmass;
+ for (i = 0; i < num_particles; i++) kd->p[i].np_index = i;
HC my_comm;
my_comm.s = newslice();
@@ -145,17 +150,11 @@
// All we need to do is provide density and group information.
// Tags (as per writetagsf77) are in gl.s->ntag+1 and there are gl.s->numlist of them.
- PyArrayObject *particle_density = (PyArrayObject *)
- PyArray_SimpleNewFromDescr(1, PyArray_DIMS(xpos),
- PyArray_DescrFromType(NPY_FLOAT64));
PyArrayObject *particle_group_id = (PyArrayObject *)
PyArray_SimpleNewFromDescr(1, PyArray_DIMS(xpos),
PyArray_DescrFromType(NPY_INT32));
for (i = 0; i < num_particles; i++) {
- // Density is in kd->p[i].fDensity
- *(npy_float64*)(PyArray_GETPTR1(particle_density, i)) =
- (npy_float64) kd->p[i].fDensity;
// tag is in gl.s->ntag[i+1]
*(npy_int32*)(PyArray_GETPTR1(particle_group_id, i)) =
(npy_int32) my_comm.s->ntag[i+1];
Modified: trunk/yt/lagos/hop/hop_hop.c
==============================================================================
--- trunk/yt/lagos/hop/hop_hop.c (original)
+++ trunk/yt/lagos/hop/hop_hop.c Tue May 5 22:58:17 2009
@@ -22,6 +22,7 @@
#include "kd.h"
#include "hop.h"
#include "smooth.h"
+#include "hop_numpy.h"
//#include "macros_and_parameters.h"
#define ISYM "d"
@@ -169,13 +170,13 @@
int j;
float totalmass;
for (j=0,totalmass=0.0; j<nSmooth; j++)
- totalmass += smx->kd->p[pList[j]].fMass;
- smx->kd->p[pi].fDensity = totalmass*0.75*M_1_PI/
+ totalmass += NP_MASS(smx->kd, pList[j]);
+ NP_DENS(smx->kd, pi) = totalmass*0.75*M_1_PI/
smx->pfBall2[pi]/sqrt(smx->pfBall2[pi]);
#else
/* This case is simple: the total mass is nSmooth times the mass
per particle */
- smx->kd->p[pi].fDensity = (nSmooth)*smx->kd->fMass*0.75*M_1_PI/
+ NP_DENS(smx->kd, pi) = (nSmooth)*smx->kd->fMass*0.75*M_1_PI/
smx->pfBall2[pi]/sqrt(smx->pfBall2[pi]);
#endif
return;
@@ -198,7 +199,7 @@
void ssort(float X[], int Y[], int N, int KFLAG);
/* If the density is less than the threshold requirement, then assign 0 */
- if (smx->kd->p[pi].fDensity<smx->fDensThresh) {
+ if (NP_DENS(smx->kd, pi)<smx->fDensThresh) {
smx->kd->p[pi].iHop = 0;
return;
}
@@ -218,9 +219,9 @@
max = 0;
maxden = 0.0;
for (i=0;i<search;++i) {
- if (smx->kd->p[pList[i]].fDensity>maxden) {
+ if (NP_DENS(smx->kd, pList[i])>maxden) {
max = i;
- maxden = smx->kd->p[pList[i]].fDensity;
+ maxden = NP_DENS(smx->kd, pList[i]);
}
}
smx->kd->p[pi].iHop = -1-pList[max];
@@ -409,8 +410,8 @@
/* It's in a different group; we need to connect the two */
if (group<g2) g1=group;
else {g1=g2; g2=group;}
- averdensity = 0.5*(smx->kd->p[pi].fDensity +
- smx->kd->p[pList[j]].fDensity);
+ averdensity = 0.5*(NP_DENS(smx->kd, pi) +
+ NP_DENS(smx->kd, pList[j]));
hashpoint = (g1+1)*g2; /* Avoid multiplying by 0 */
hashpoint = hashpoint % smx->nHashLength;
hp = smx->hash+hashpoint;
@@ -492,15 +493,15 @@
** Calculate Bounds.
*/
for (j=0;j<3;++j) {
- bnd.fMin[j] = kd->p[0].r[j];
- bnd.fMax[j] = kd->p[0].r[j];
+ bnd.fMin[j] = NP_POS(kd, 0, j);
+ bnd.fMax[j] = NP_POS(kd, 0, j);
}
for (i=1;i<kd->nActive;++i) {
for (j=0;j<3;++j) {
- if (bnd.fMin[j] > kd->p[i].r[j])
- bnd.fMin[j] = kd->p[i].r[j];
- else if (bnd.fMax[j] < kd->p[i].r[j])
- bnd.fMax[j] = kd->p[i].r[j];
+ if (bnd.fMin[j] > NP_POS(kd, i, j))
+ bnd.fMin[j] = NP_POS(kd, i, j);
+ else if (bnd.fMax[j] < NP_POS(kd, i, j))
+ bnd.fMax[j] = NP_POS(kd, i, j);
}
}
kd->bnd = bnd;
@@ -522,7 +523,7 @@
//s->ID = ivector(1,s->numlist);
for (j=0;j<smx->kd->nActive;j++) {
//s->ID[1+j] = smx->kd->p[j].iID; /* S Skory's addition */
- if (smx->kd->p[j].fDensity < densthres) s->ntag[j+1] = -1;
+ if (NP_DENS(smx->kd,j) < densthres) s->ntag[j+1] = -1;
else s->ntag[j+1] = smx->kd->p[j].iHop;
}
@@ -545,7 +546,7 @@
my_comm->gdensity = vector(0,smx->nGroups-1);
for (j=0;j<smx->nGroups;j++) {
den = smx->densestingroup[j];
- my_comm->gdensity[j]=smx->kd->p[den].fDensity;
+ my_comm->gdensity[j]=NP_DENS(smx->kd, den);
}
int nb = 0;
for (j=0, hp=smx->hash;j<smx->nHashLength; j++,hp++)
Modified: trunk/yt/lagos/hop/hop_kd.c
==============================================================================
--- trunk/yt/lagos/hop/hop_kd.c (original)
+++ trunk/yt/lagos/hop/hop_kd.c Tue May 5 22:58:17 2009
@@ -16,6 +16,7 @@
#include <sys/resource.h>
#include <assert.h>
#include "kd.h"
+#include "hop_numpy.h"
//#include "macros_and_parameters.h"
/* #include "tipsydefs.h" */ /* Don't need this, since I removed kdReadTipsy()*/
@@ -55,7 +56,7 @@
*/
int kdMedianJst(KD kd,int d,int l,int u)
{
- float fm;
+ npy_float64 fm;
int i,k,m;
PARTICLE *p,t;
@@ -64,15 +65,15 @@
m = k;
while (l < u) {
m = (l+u)/2;
- fm = p[m].r[d];
+ fm = NP_POS(kd, m, d);
t = p[m];
p[m] = p[u];
p[u] = t;
i = u-1;
m = l;
- while (p[m].r[d] < fm) ++m;
+ while (NP_POS(kd, m, d) < fm) ++m;
while (m < i) {
- while (p[i].r[d] >= fm) if (--i == m) break;
+ while (NP_POS(kd, i, d) >= fm) if (--i == m) break;
/*
** Swap
*/
@@ -80,7 +81,7 @@
p[m] = p[i];
p[i] = t;
--i;
- while (p[m].r[d] < fm) ++m;
+ while (NP_POS(kd, m, d) < fm) ++m;
}
t = p[m];
p[m] = p[u];
@@ -129,15 +130,15 @@
l = c[iCell].pLower;
u = c[iCell].pUpper;
for (j=0;j<3;++j) {
- c[iCell].bnd.fMin[j] = kd->p[u].r[j];
- c[iCell].bnd.fMax[j] = kd->p[u].r[j];
+ c[iCell].bnd.fMin[j] = NP_POS(kd, u, j);
+ c[iCell].bnd.fMax[j] = NP_POS(kd, u, j);
}
for (pj=l;pj<u;++pj) {
for (j=0;j<3;++j) {
- if (kd->p[pj].r[j] < c[iCell].bnd.fMin[j])
- c[iCell].bnd.fMin[j] = kd->p[pj].r[j];
- if (kd->p[pj].r[j] > c[iCell].bnd.fMax[j])
- c[iCell].bnd.fMax[j] = kd->p[pj].r[j];
+ if (NP_POS(kd, pj, j) < c[iCell].bnd.fMin[j])
+ c[iCell].bnd.fMin[j] = NP_POS(kd, pj, j);
+ if (NP_POS(kd, pj, j) > c[iCell].bnd.fMax[j])
+ c[iCell].bnd.fMax[j] = NP_POS(kd, pj, j);
}
}
}
@@ -179,7 +180,7 @@
}
c[i].iDim = d;
m = kdMedianJst(kd,d,c[i].pLower,c[i].pUpper);
- c[i].fSplit = kd->p[m].r[d];
+ c[i].fSplit = NP_POS(kd, m, d);
c[LOWER(i)].bnd = c[i].bnd;
c[LOWER(i)].bnd.fMax[d] = c[i].fSplit;
c[LOWER(i)].pLower = c[i].pLower;
Modified: trunk/yt/lagos/hop/hop_smooth.c
==============================================================================
--- trunk/yt/lagos/hop/hop_smooth.c (original)
+++ trunk/yt/lagos/hop/hop_smooth.c Tue May 5 22:58:17 2009
@@ -19,6 +19,7 @@
#include <assert.h>
#include "smooth.h"
#include "kd.h"
+#include "hop_numpy.h"
#define ISYM "d"
#define GSYM "g"
@@ -60,7 +61,7 @@
** Initialize arrays for calculated quantities.--DJE
*/
for (pi=0;pi<smx->kd->nActive;++pi) {
- smx->kd->p[pi].fDensity = 0.0;
+ NP_DENS(smx->kd, pi) = 0.0;
smx->kd->p[pi].iHop = 0;
}
*psmx = smx;
@@ -108,9 +109,9 @@
** Now start the search from the bucket given by cell!
*/
for (pj=c[cell].pLower;pj<=c[cell].pUpper;++pj) {
- dx = x - p[pj].r[0];
- dy = y - p[pj].r[1];
- dz = z - p[pj].r[2];
+ dx = x - NP_POS(smx->kd, pj, 0);
+ dy = y - NP_POS(smx->kd, pj, 1);
+ dz = z - NP_POS(smx->kd, pj, 2);
fDist2 = dx*dx + dy*dy + dz*dz;
if (fDist2 < fBall2) {
if (smx->iMark[pj]) continue;
@@ -140,9 +141,9 @@
}
else {
for (pj=c[cp].pLower;pj<=c[cp].pUpper;++pj) {
- dx = sx - p[pj].r[0];
- dy = sy - p[pj].r[1];
- dz = sz - p[pj].r[2];
+ dx = sx - NP_POS(smx->kd, pj, 0);
+ dy = sy - NP_POS(smx->kd, pj, 1);
+ dz = sz - NP_POS(smx->kd, pj, 2);
fDist2 = dx*dx + dy*dy + dz*dz;
if (fDist2 < fBall2) {
if (smx->iMark[pj]) continue;
@@ -197,9 +198,9 @@
}
else {
for (pj=c[cp].pLower;pj<=c[cp].pUpper;++pj) {
- dx = sx - p[pj].r[0];
- dy = sy - p[pj].r[1];
- dz = sz - p[pj].r[2];
+ dx = sx - NP_POS(smx->kd, pj, 0);
+ dy = sy - NP_POS(smx->kd, pj, 1);
+ dz = sz - NP_POS(smx->kd, pj, 2);
fDist2 = dx*dx + dy*dy + dz*dz;
if (fDist2 < fBall2) {
smx->fList[nCnt] = fDist2;
@@ -229,6 +230,7 @@
int cell;
int pi,pin,pj,pNext,nCnt,nSmooth;
float dx,dy,dz,x,y,z,h2,ax,ay,az;
+ float temp_ri[3];
for (pi=0;pi<smx->kd->nActive;++pi) {
@@ -271,9 +273,9 @@
if (pNext == smx->kd->nActive) break;
pi = pNext;
++pNext;
- x = p[pi].r[0];
- y = p[pi].r[1];
- z = p[pi].r[2];
+ x = NP_POS(smx->kd, pi, 0);
+ y = NP_POS(smx->kd, pi, 1);
+ z = NP_POS(smx->kd, pi, 2);
/* printf("%"ISYM": %"GSYM" %"GSYM" %"GSYM"\n", pi, x, y, z); */
/*
** First find the "local" Bucket.
@@ -281,7 +283,7 @@
*/
cell = ROOT;
while (cell < smx->kd->nSplit) {
- if (p[pi].r[c[cell].iDim] < c[cell].fSplit)
+ if (NP_POS(smx->kd, pi, c[cell].iDim) < c[cell].fSplit)
cell = LOWER(cell);
else
cell = UPPER(cell);
@@ -299,9 +301,9 @@
pj = smx->kd->nActive - nSmooth;
for (pq=smx->pq;pq<=pqLast;++pq) {
smx->iMark[pj] = 1;
- dx = x - p[pj].r[0];
- dy = y - p[pj].r[1];
- dz = z - p[pj].r[2];
+ dx = x - NP_POS(smx->kd, pj, 0);
+ dy = y - NP_POS(smx->kd, pj, 1);
+ dz = z - NP_POS(smx->kd, pj, 2);
pq->fKey = dx*dx + dy*dy + dz*dz;
pq->p = pj++;
pq->ax = 0.0;
@@ -315,17 +317,17 @@
** Calculate the priority queue using the previous particles!
*/
pi = pin;
- x = p[pi].r[0];
- y = p[pi].r[1];
- z = p[pi].r[2];
+ x = NP_POS(smx->kd, pi, 0);
+ y = NP_POS(smx->kd, pi, 1);
+ z = NP_POS(smx->kd, pi, 2);
smx->pqHead = NULL;
for (pq=smx->pq;pq<=pqLast;++pq) {
pq->ax -= ax;
pq->ay -= ay;
pq->az -= az;
- dx = x + pq->ax - p[pq->p].r[0];
- dy = y + pq->ay - p[pq->p].r[1];
- dz = z + pq->az - p[pq->p].r[2];
+ dx = x + pq->ax - NP_POS(smx->kd, pq->p, 0);
+ dy = y + pq->ay - NP_POS(smx->kd, pq->p, 1);
+ dz = z + pq->az - NP_POS(smx->kd, pq->p, 2);
pq->fKey = dx*dx + dy*dy + dz*dz;
}
PQ_BUILD(smx->pq,nSmooth,smx->pqHead);
@@ -333,7 +335,10 @@
ay = 0.0;
az = 0.0;
}
- smBallSearch(smx,smx->pqHead->fKey,p[pi].r);
+ temp_ri[0] = NP_POS(smx->kd, pi, 0);
+ temp_ri[1] = NP_POS(smx->kd, pi, 1);
+ temp_ri[2] = NP_POS(smx->kd, pi, 2);
+ smBallSearch(smx,smx->pqHead->fKey,temp_ri);
smx->pfBall2[pi] = smx->pqHead->fKey;
/*
** Pick next particle, 'pin'.
@@ -364,6 +369,7 @@
{
PARTICLE *p;
int pi,nSmooth;
+ float temp_ri[3];
p = smx->kd->p;
for (pi=0;pi<smx->kd->nActive;++pi) {
@@ -372,7 +378,10 @@
** Do a Ball Gather at the radius of the most distant particle
** which is smDensity sets in smx->pBall[pi].
*/
- nSmooth = smBallGather(smx,smx->pfBall2[pi],p[pi].r);
+ temp_ri[0] = NP_POS(smx->kd, pi, 0);
+ temp_ri[1] = NP_POS(smx->kd, pi, 1);
+ temp_ri[2] = NP_POS(smx->kd, pi, 2);
+ nSmooth = smBallGather(smx,smx->pfBall2[pi],temp_ri);
(*fncSmooth)(smx,pi,nSmooth,smx->pList,smx->fList);
}
}
@@ -392,12 +401,12 @@
if (r2 < 1.0) rs = (1.0 - 0.75*rs*r2);
else rs = 0.25*rs*rs*rs;
#ifdef DIFFERENT_MASSES
- fDensity += rs*smx->kd->p[pj].fMass;
+ fDensity += rs*NP_MASS(smx->kd, pj);
#else
fDensity += rs*smx->kd->fMass;
#endif
}
- smx->kd->p[pi].fDensity = M_1_PI*sqrt(ih2)*ih2*fDensity;
+ NP_DENS(smx->kd, pi) = M_1_PI*sqrt(ih2)*ih2*fDensity;
}
@@ -416,8 +425,8 @@
else rs = 0.25*rs*rs*rs;
rs *= fNorm;
#ifdef DIFFERENT_MASSES
- smx->kd->p[pi].fDensity += rs*smx->kd->p[pj].fMass;
- smx->kd->p[pj].fDensity += rs*smx->kd->p[pi].fMass;
+ NP_DENS(smx->kd, pi) += rs*NP_MASS(smx->kd, pj);
+ NP_DENS(smx->kd, pj) += rs*NP_MASS(smx->kd, pi);
#else
smx->kd->p[pi].fDensity += rs*smx->kd->fMass;
smx->kd->p[pj].fDensity += rs*smx->kd->fMass;
@@ -438,7 +447,7 @@
for (i=0;i<smx->kd->nGas;++i) {
if (smx->kd->bGas) {
if (IMARK)
- fprintf(fp,"%.8"GSYM"\n",smx->kd->p[iCnt].fDensity);
+ fprintf(fp,"%.8"GSYM"\n",NP_DENS(smx->kd, iCnt));
else fprintf(fp,"0\n");
++iCnt;
}
@@ -447,7 +456,7 @@
for (i=0;i<smx->kd->nDark;++i) {
if (smx->kd->bDark) {
if (IMARK)
- fprintf(fp,"%.8"GSYM"\n",smx->kd->p[iCnt].fDensity);
+ fprintf(fp,"%.8"GSYM"\n",NP_DENS(smx->kd, iCnt));
else fprintf(fp,"0\n");
++iCnt;
}
@@ -456,7 +465,7 @@
for (i=0;i<smx->kd->nStar;++i) {
if (smx->kd->bStar) {
if (IMARK)
- fprintf(fp,"%.8"GSYM"\n",smx->kd->p[iCnt].fDensity);
+ fprintf(fp,"%.8"GSYM"\n",NP_DENS(smx->kd, iCnt));
else fprintf(fp,"0\n");
++iCnt;
}
Modified: trunk/yt/lagos/hop/kd.h
==============================================================================
--- trunk/yt/lagos/hop/kd.h (original)
+++ trunk/yt/lagos/hop/kd.h Tue May 5 22:58:17 2009
@@ -18,6 +18,9 @@
#ifndef KD_HINCLUDED
#define KD_HINCLUDED
+#include "Python.h"
+#include "numpy/ndarrayobject.h"
+
#define ROOT 1
#define LOWER(i) (i<<1)
#define UPPER(i) ((i<<1)+1)
@@ -34,8 +37,11 @@
#define STAR 4
typedef struct Particle {
- float r[3];
+ int np_index;
+ int iHop;
int iOrder;
+#if 0
+ float r[3];
float fDensity;
// int iID; /* the real ID of the particle S. Skory */
int iHop; /* DJE: The number of the highest-density neighbor;
@@ -49,6 +55,7 @@
/* int iMark; */
/* float vMean[3]; */
/* float fVelDisp2; */
+#endif
} PARTICLE;
typedef struct bndBound {
@@ -84,6 +91,10 @@
KDN *kdNodes;
int uSecond;
int uMicro;
+ PyArrayObject *np_densities;
+ PyArrayObject *np_pos[3];
+ PyArrayObject *np_masses;
+ float totalmass;
} * KD;
More information about the yt-svn
mailing list