Homec4science
Diffusion pNbody (master) Browse

pNbody/src/streelib/treelib.c
6fbb9ed07f94master

Last Change
Raw File

treelib.c
View Options

#include <Python.h>
#include <numpy/arrayobject.h>
#include "structmember.h"
#define FLOAT float
typedef long long peanokey;
#define MAXTOPNODES 200000
#define MAX_REAL_NUMBER 1e37
#define BITS_PER_DIMENSION 18
#define PEANOCELLS (((peanokey)1)<<(3*BITS_PER_DIMENSION))
#ifndef TWODIMS
#define NUMDIMS 3 /*!< For 3D-normalized kernel */
#define KERNEL_COEFF_1 2.546479089470 /*!< Coefficients for SPH spline kernel and its derivative */
#define KERNEL_COEFF_2 15.278874536822
#define KERNEL_COEFF_3 45.836623610466
#define KERNEL_COEFF_4 30.557749073644
#define KERNEL_COEFF_5 5.092958178941
#define KERNEL_COEFF_6 (-15.278874536822)
#define NORM_COEFF 4.188790204786 /*!< Coefficient for kernel normalization. Note: 4.0/3 * PI = 4.188790204786 */
#else
#define NUMDIMS 2 /*!< For 2D-normalized kernel */
#define KERNEL_COEFF_1 (5.0/7*2.546479089470) /*!< Coefficients for SPH spline kernel and its derivative */
#define KERNEL_COEFF_2 (5.0/7*15.278874536822)
#define KERNEL_COEFF_3 (5.0/7*45.836623610466)
#define KERNEL_COEFF_4 (5.0/7*30.557749073644)
#define KERNEL_COEFF_5 (5.0/7*5.092958178941)
#define KERNEL_COEFF_6 (5.0/7*(-15.278874536822))
#define NORM_COEFF M_PI /*!< Coefficient for kernel normalization. */
#endif
/******************************************************************************
SYSTEM
*******************************************************************************/
/*! returns the maximum of two double
*/
double dmax(double x, double y)
{
if(x > y)
return x;
else
return y;
}
/******************************************************************************
TREE STRUCTURE
*******************************************************************************/
struct global_data_all_processes
{
long long TotNumPart;
long long TotN_gas;
int MaxPart;
int MaxPartSph;
double SofteningTable[6];
double ForceSoftening[6];
double PartAllocFactor;
double TreeAllocFactor;
double ErrTolTheta;
double DesNumNgb;
double MaxNumNgbDeviation;
double MinGasHsmlFractional;
double MinGasHsml;
};
struct particle_data
{
FLOAT Pos[3]; /*!< particle position at its current time */
FLOAT Mass; /*!< particle mass */
FLOAT Vel[3]; /*!< particle velocity at its current time */
FLOAT GravAccel[3]; /*!< particle acceleration due to gravity */
FLOAT Potential; /*!< gravitational potential */
FLOAT OldAcc; /*!< magnitude of old gravitational force. Used in relative opening criterion */
int Type; /*!< flags particle type. 0=gas, 1=halo, 2=disk, 3=bulge, 4=stars, 5=bndry */
int Ti_endstep; /*!< marks start of current timestep of particle on integer timeline */
int Ti_begstep; /*!< marks end of current timestep of particle on integer timeline */
FLOAT Density;
FLOAT Observable;
int next; /* yr : indice of the next particule in the node list */
};
struct topnode_exchange
{
peanokey Startkey;
int Count;
};
struct topnode_data
{
int Daughter; /*!< index of first daughter cell (out of 8) of top-level node */
int Pstart; /*!< for the present top-level node, this gives the index of the first node in the concatenated list of topnodes collected from all processors */
int Blocks; /*!< for the present top-level node, this gives the number of corresponding nodes in the concatenated list of topnodes collected from all processors */
int Leaf; /*!< if the node is a leaf, this gives its number when all leaves are traversed in Peano-Hilbert order */
peanokey Size; /*!< number of Peano-Hilbert mesh-cells represented by top-level node */
peanokey StartKey; /*!< first Peano-Hilbert key in top-level node */
long long Count; /*!< counts the number of particles in this top-level node */
};
typedef struct {
PyObject_HEAD
PyObject *first; /* first name */
PyObject *list;
int number;
/* allvars */
int Numnodestree;
int MaxNodes;
int NumPart;
int N_gas;
long long Ntype[6];
int ThisTask;
int NTask;
struct NODE *Nodes_base;
struct NODE *Nodes;
struct topnode_data *TopNodes;
peanokey *Key;
peanokey *KeySorted;
int *DomainNodeIndex;
struct global_data_all_processes All;
struct particle_data *P;
double DomainCorner[3];
double DomainCenter[3];
double DomainLen;
double DomainFac;
int NTopnodes;
int *Nextnode;
int *Father;
int NTopleaves;
/* allvars.c */
int NtypeLocal[6];
/* domain */
long long maxload, maxloadsph;
//int *list_NumPart;
//int *list_N_gas;
/* force */
int last;
/* ngb */
int *Ngblist;
} Tree;
/******************************************************************************
ENDRUN
*******************************************************************************/
/*! This function aborts the simulations. If a single processors wants an
* immediate termination, the function needs to be called with ierr>0. A
* bunch of MPI-error messages may also appear in this case. For ierr=0,
* MPI is gracefully cleaned up, but this requires that all processors
* call endrun().
*/
void endrun(Tree *self,int ierr)
{
if(ierr)
{
printf("task %d: endrun called with an error level of %d\n\n\n", self->ThisTask, ierr);
fflush(stdout);
//#ifdef DEBUG
// terminate_processes();
// raise(SIGABRT);
// sleep(60);
//#else
// MPI_Abort(MPI_COMM_WORLD, ierr);
//#endif
exit(0);
}
// MPI_Finalize();
exit(0);
};
/******************************************************************************
PEANO THINGS
*******************************************************************************/
static int quadrants[24][2][2][2] = {
/* rotx=0, roty=0-3 */
{{{0, 7}, {1, 6}}, {{3, 4}, {2, 5}}},
{{{7, 4}, {6, 5}}, {{0, 3}, {1, 2}}},
{{{4, 3}, {5, 2}}, {{7, 0}, {6, 1}}},
{{{3, 0}, {2, 1}}, {{4, 7}, {5, 6}}},
/* rotx=1, roty=0-3 */
{{{1, 0}, {6, 7}}, {{2, 3}, {5, 4}}},
{{{0, 3}, {7, 4}}, {{1, 2}, {6, 5}}},
{{{3, 2}, {4, 5}}, {{0, 1}, {7, 6}}},
{{{2, 1}, {5, 6}}, {{3, 0}, {4, 7}}},
/* rotx=2, roty=0-3 */
{{{6, 1}, {7, 0}}, {{5, 2}, {4, 3}}},
{{{1, 2}, {0, 3}}, {{6, 5}, {7, 4}}},
{{{2, 5}, {3, 4}}, {{1, 6}, {0, 7}}},
{{{5, 6}, {4, 7}}, {{2, 1}, {3, 0}}},
/* rotx=3, roty=0-3 */
{{{7, 6}, {0, 1}}, {{4, 5}, {3, 2}}},
{{{6, 5}, {1, 2}}, {{7, 4}, {0, 3}}},
{{{5, 4}, {2, 3}}, {{6, 7}, {1, 0}}},
{{{4, 7}, {3, 0}}, {{5, 6}, {2, 1}}},
/* rotx=4, roty=0-3 */
{{{6, 7}, {5, 4}}, {{1, 0}, {2, 3}}},
{{{7, 0}, {4, 3}}, {{6, 1}, {5, 2}}},
{{{0, 1}, {3, 2}}, {{7, 6}, {4, 5}}},
{{{1, 6}, {2, 5}}, {{0, 7}, {3, 4}}},
/* rotx=5, roty=0-3 */
{{{2, 3}, {1, 0}}, {{5, 4}, {6, 7}}},
{{{3, 4}, {0, 7}}, {{2, 5}, {1, 6}}},
{{{4, 5}, {7, 6}}, {{3, 2}, {0, 1}}},
{{{5, 2}, {6, 1}}, {{4, 3}, {7, 0}}}
};
static int rotxmap_table[24] = { 4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14, 15, 0, 1, 2, 3, 17, 18, 19, 16, 23, 20, 21, 22
};
static int rotymap_table[24] = { 1, 2, 3, 0, 16, 17, 18, 19,
11, 8, 9, 10, 22, 23, 20, 21, 14, 15, 12, 13, 4, 5, 6, 7
};
static int rotx_table[8] = { 3, 0, 0, 2, 2, 0, 0, 1 };
static int roty_table[8] = { 0, 1, 1, 2, 2, 3, 3, 0 };
static int sense_table[8] = { -1, -1, -1, +1, +1, -1, -1, -1 };
static int flag_quadrants_inverse = 1;
static char quadrants_inverse_x[24][8];
static char quadrants_inverse_y[24][8];
static char quadrants_inverse_z[24][8];
/*! This function computes a Peano-Hilbert key for an integer triplet (x,y,z),
* with x,y,z in the range between 0 and 2^bits-1.
*/
peanokey peano_hilbert_key(int x, int y, int z, int bits)
{
int i, quad, bitx, bity, bitz;
int mask, rotation, rotx, roty, sense;
peanokey key;
mask = 1 << (bits - 1);
key = 0;
rotation = 0;
sense = 1;
for(i = 0; i < bits; i++, mask >>= 1)
{
bitx = (x & mask) ? 1 : 0;
bity = (y & mask) ? 1 : 0;
bitz = (z & mask) ? 1 : 0;
quad = quadrants[rotation][bitx][bity][bitz];
key <<= 3;
key += (sense == 1) ? (quad) : (7 - quad);
rotx = rotx_table[quad];
roty = roty_table[quad];
sense *= sense_table[quad];
while(rotx > 0)
{
rotation = rotxmap_table[rotation];
rotx--;
}
while(roty > 0)
{
rotation = rotymap_table[rotation];
roty--;
}
}
return key;
}
/******************************************************************************
DOMAIN THINGS
*******************************************************************************/
/*! This routine finds the extent of the global domain grid.
*/
void domain_findExtent(Tree *self)
{
int i, j;
double len, xmin[3], xmax[3], xmin_glob[3], xmax_glob[3];
/* determine local extension */
for(j = 0; j < 3; j++)
{
xmin[j] = MAX_REAL_NUMBER;
xmax[j] = -MAX_REAL_NUMBER;
}
for(i = 0; i < self->NumPart; i++)
{
for(j = 0; j < 3; j++)
{
if(xmin[j] > self->P[i].Pos[j])
xmin[j] = self->P[i].Pos[j];
if(xmax[j] < self->P[i].Pos[j])
xmax[j] = self->P[i].Pos[j];
}
}
for(j = 0; j < 3; j++)
{
xmin_glob[j] = xmin[j];
xmax_glob[j] = xmax[j];
}
len = 0;
for(j = 0; j < 3; j++)
if(xmax_glob[j] - xmin_glob[j] > len)
len = xmax_glob[j] - xmin_glob[j];
len *= 1.001;
for(j = 0; j < 3; j++)
{
self->DomainCenter[j] = 0.5 * (xmin_glob[j] + xmax_glob[j]);
self->DomainCorner[j] = 0.5 * (xmin_glob[j] + xmax_glob[j]) - 0.5 * len;
}
self->DomainLen = len;
self->DomainFac = 1.0 / len * (((peanokey) 1) << (BITS_PER_DIMENSION));
}
/******************************************************************************
FORCE THINGS
*******************************************************************************/
struct NODE
{
FLOAT len; /*!< sidelength of treenode */
FLOAT center[3]; /*!< geometrical center of node */
int count[8]; /* number of particles contained by each daughter node*/
int first[8]; /* index of first particle inside the daughter node*/
int last[8]; /* index of last particle inside the daughter node*/
union
{
int suns[8]; /*!< temporary pointers to daughter nodes */
struct
{
FLOAT s[3]; /*!< center of mass of node */
FLOAT mass; /*!< mass of node */
int bitflags; /*!< a bit-field with various information on the node */
int sibling; /*!< this gives the next node in the walk in case the current node can be used */
int nextnode; /*!< this gives the next node in case the current node needs to be opened */
int father; /*!< this gives the parent node of each node (or -1 if we have the root node) */
}
d;
}
u;
};
int nNodes=0;
PyArrayObject *NodesLen,*NodesCenter;
/*
walk the tree and gather info : center, size
this function uses
u.suns[j]
so, it must be used before force_update_node_recursive
*/
void tree_walk_recursive(Tree *self,int no, int sib, int father)
{
int j,p;
int nextsib;
if(no >= self->All.MaxPart && no < self->All.MaxPart + self->MaxNodes)
{
nNodes++;
/* here, we can gather info */
//printf("%g %g %g %g \n",self->Nodes[no].len,self->Nodes[no].center[0],self->Nodes[no].center[1],self->Nodes[no].center[2]);
self->Nodes[no].center[0];
/* loop over */
for(j = 0; j < 8; j++)
{
if((p = self->Nodes[no].u.suns[j]) >= 0)
{
/* useless ... ??? */
nextsib = sib;
/* go to next node */
tree_walk_recursive(self,p,nextsib,no);
}
}
}
else /* single particle or pseudo particle */
{
//printf("found single particle or pseudo particle\n");
}
}
/*
walk the tree and gather info : center, size
this function uses
Nextnode and self->Nodes[no].u.d.nextnode
so, it must be used after force_update_node_recursive
*/
void tree_walk(Tree *self,int no, int sib, int father)
{
int j,p;
int nextsib;
while(no >= 0)
{
if(no >= self->All.MaxPart && no < self->All.MaxPart + self->MaxNodes)
{
nNodes++;
//printf("%g %g %g %g \n",self->Nodes[no].len,self->Nodes[no].center[0],self->Nodes[no].center[1],self->Nodes[no].center[2]);
no = self->Nodes[no].u.d.nextnode;
}
else /* single particle or pseudo particle */
{
//printf("found single particle or pseudo particle\n");
no = self->Nextnode[no];
}
}
}
/*
This function simply walk the tree and gather info, like
the size of nodes and their position.
The output is printed on the standard io.
*/
void tree_walk_and_gather(Tree *self,int mode)
{
int no;
int sib=-1;
int father=-1;
nNodes=0;
no = self->All.MaxPart; /* root node */
printf("start tree_walk_and_gather...\n");
if (mode==0)
tree_walk_recursive(self,no,sib,father);
else
tree_walk(self,no,sib,father);
printf("number of nodes : %d\n",nNodes);
printf("tree_walk_and_gather done.\n");
}
/*! This routine computes the gravitational force for a given local
* particle, or for a particle in the communication buffer. Depending on
* the value of TypeOfOpeningCriterion, either the geometrical BH
* cell-opening criterion, or the `relative' opening criterion is used.
*/
int force_treeevaluate_local(Tree* self,double pos_x, double pos_y, double pos_z, double h, double *acc_x, double *acc_y, double *acc_z)
{
struct NODE *nop = 0;
int no, ninteractions, ptype;
double r2, dx, dy, dz, mass, r, fac, u, h_inv, h3_inv;
#if defined(UNEQUALSOFTENINGS) && !defined(ADAPTIVE_GRAVSOFT_FORGAS)
int maxsofttype;
#endif
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
double soft = 0;
#endif
#ifdef PERIODIC
double boxsize, boxhalf;
boxsize = All.BoxSize;
boxhalf = 0.5 * All.BoxSize;
#endif
*acc_x = 0;
*acc_y = 0;
*acc_z = 0;
ninteractions = 0;
ptype = 0; /* ????? */
h_inv = 1.0 / h;
h3_inv = h_inv * h_inv * h_inv;
no = self->All.MaxPart; /* root node */
while(no >= 0)
{
if(no < self->All.MaxPart) /* single particle */
{
/* the index of the node is the index of the particle */
/* observe the sign */
dx = self->P[no].Pos[0] - pos_x;
dy = self->P[no].Pos[1] - pos_y;
dz = self->P[no].Pos[2] - pos_z;
mass = self->P[no].Mass;
}
else /* inernal node */
{
nop = &self->Nodes[no];
dx = nop->u.d.s[0] - pos_x;
dy = nop->u.d.s[1] - pos_y;
dz = nop->u.d.s[2] - pos_z;
mass = nop->u.d.mass;
}
r2 = dx * dx + dy * dy + dz * dz;
if(no < self->All.MaxPart) /* single particle */
{
no = self->Nextnode[no];
}
else /* inernal node */
{
if(nop->len * nop->len > r2 * self->All.ErrTolTheta * self->All.ErrTolTheta)
{
/* open cell */
no = nop->u.d.nextnode;
continue;
}
no = nop->u.d.sibling; /* ok, node can be used */
}
/* here we use a plummer softening */
if (r2>0)
{
fac = mass/pow(r2+h*h ,3.0/2.0);
*acc_x += dx * fac;
*acc_y += dy * fac;
*acc_z += dz * fac;
}
ninteractions++;
}
return ninteractions;
}
/*! This routine computes the gravitational potential by walking the
* tree. The same opening criteria is used as for the gravitational force
* walk.
*/
double force_treeevaluate_local_potential(Tree* self,double pos_x, double pos_y, double pos_z, double h)
{
struct NODE *nop = 0;
int no, ptype;
double r2, dx, dy, dz, mass, r, u, h_inv, wp;
double pot;
#if defined(UNEQUALSOFTENINGS) && !defined(ADAPTIVE_GRAVSOFT_FORGAS)
int maxsofttype;
#endif
//#ifdef ADAPTIVE_GRAVSOFT_FORGAS
// double soft = 0;
//#endif
#ifdef PERIODIC
double boxsize, boxhalf;
boxsize = All.BoxSize;
boxhalf = 0.5 * All.BoxSize;
#endif
pot = 0;
ptype = 0; /* ????? */
h_inv = 1.0 / h;
no = self->All.MaxPart;
while(no >= 0)
{
if(no < self->All.MaxPart) /* single particle */
{
/* the index of the node is the index of the particle */
/* observe the sign */
dx = self->P[no].Pos[0] - pos_x;
dy = self->P[no].Pos[1] - pos_y;
dz = self->P[no].Pos[2] - pos_z;
mass = self->P[no].Mass;
}
else /* inernal node */
{
nop = &self->Nodes[no];
dx = nop->u.d.s[0] - pos_x;
dy = nop->u.d.s[1] - pos_y;
dz = nop->u.d.s[2] - pos_z;
mass = nop->u.d.mass;
}
r2 = dx * dx + dy * dy + dz * dz;
if(no < self->All.MaxPart) /* single particle */
{
no = self->Nextnode[no];
}
else /* inernal node */
{
if(self->All.ErrTolTheta) /* check Barnes-Hut opening criterion */
{
if(nop->len * nop->len > r2 * self->All.ErrTolTheta * self->All.ErrTolTheta)
{
/* open cell */
no = nop->u.d.nextnode;
continue;
}
}
no = nop->u.d.sibling; /* node can be used */
}
/* here we use a plummer softening */
if (r2>0)
pot += -mass / sqrt(r2+h*h);
}
/* return result */
return pot;
}
/*! this routine determines the multipole moments for a given internal node
* and all its subnodes using a recursive computation. The result is
* stored in the Nodes[] structure in the sequence of this tree-walk.
*
* Note that the bitflags-variable for each node is used to store in the
* lowest bits some special information: Bit 0 flags whether the node
* belongs to the top-level tree corresponding to the domain
* decomposition, while Bit 1 signals whether the top-level node is
* dependent on local mass.
*
* If UNEQUALSOFTENINGS is set, bits 2-4 give the particle type with
* the maximum softening among the particles in the node, and bit 5
* flags whether the node contains any particles with lower softening
* than that.
*/
void force_update_node_recursive(Tree *self,int no, int sib, int father)
{
int j, jj, p, pp, nextsib, suns[8];
FLOAT hmax;
//#ifdef UNEQUALSOFTENINGS
//#ifndef ADAPTIVE_GRAVSOFT_FORGAS
int maxsofttype, diffsoftflag;
//#else
//FLOAT maxsoft;
//#endif
//#endif
struct particle_data *pa;
double s[3], vs[3], mass;
if(no >= self->All.MaxPart && no < self->All.MaxPart + self->MaxNodes) /* internal node */
{
for(j = 0; j < 8; j++)
suns[j] = self->Nodes[no].u.suns[j]; /* this "backup" is necessary because the nextnode entry will
overwrite one element (union!) */
if(self->last >= 0)
{
if(self->last >= self->All.MaxPart) /* the previous (last) particle exists and is a node */
{
self->Nodes[self->last].u.d.nextnode = no; /* make a link to the current */
}
else /* the pervious (last) particle exists and is a particle */
self->Nextnode[self->last] = no; /* make a link to the current */
}
self->last = no;
mass = 0;
s[0] = 0;
s[1] = 0;
s[2] = 0;
vs[0] = 0;
vs[1] = 0;
vs[2] = 0;
hmax = 0;
//#ifdef UNEQUALSOFTENINGS
//#ifndef ADAPTIVE_GRAVSOFT_FORGAS
maxsofttype = 7;
diffsoftflag = 0;
//#else
// maxsoft = 0;
//#endif
//#endif
for(j = 0; j < 8; j++) /* loop over daughters */
{
if((p = suns[j]) >= 0) /* the node is not empty */
{
/* check if we have a sibling on the same level */
for(jj = j + 1; jj < 8; jj++)
if((pp = suns[jj]) >= 0)
break;
if(jj < 8) /* yes, we do */
nextsib = pp;
else
nextsib = sib;
force_update_node_recursive(self,p, nextsib, no);
if(p >= self->All.MaxPart) /* an internal node or pseudo particle */
{
mass += self->Nodes[p].u.d.mass;
s[0] += self->Nodes[p].u.d.mass * self->Nodes[p].u.d.s[0];
s[1] += self->Nodes[p].u.d.mass * self->Nodes[p].u.d.s[1];
s[2] += self->Nodes[p].u.d.mass * self->Nodes[p].u.d.s[2];
diffsoftflag |= (self->Nodes[p].u.d.bitflags >> 5) & 1;
if(maxsofttype == 7)
{
maxsofttype = (self->Nodes[p].u.d.bitflags >> 2) & 7;
}
else
{
if(((self->Nodes[p].u.d.bitflags >> 2) & 7) != 7)
{
if(self->All.ForceSoftening[((self->Nodes[p].u.d.bitflags >> 2) & 7)] >
self->All.ForceSoftening[maxsofttype])
{
maxsofttype = ((self->Nodes[p].u.d.bitflags >> 2) & 7);
diffsoftflag = 1;
}
else
{
if(self->All.ForceSoftening[((self->Nodes[p].u.d.bitflags >> 2) & 7)] <
self->All.ForceSoftening[maxsofttype])
diffsoftflag = 1;
}
}
}
}
else /* a particle */
{
pa = &self->P[p];
mass += pa->Mass;
s[0] += pa->Mass * pa->Pos[0];
s[1] += pa->Mass * pa->Pos[1];
s[2] += pa->Mass * pa->Pos[2];
vs[0] += pa->Mass * pa->Vel[0];
vs[1] += pa->Mass * pa->Vel[1];
vs[2] += pa->Mass * pa->Vel[2];
if(maxsofttype == 7)
{
maxsofttype = pa->Type;
}
else
{
if(self->All.ForceSoftening[pa->Type] > self->All.ForceSoftening[maxsofttype])
{
maxsofttype = pa->Type;
diffsoftflag = 1;
}
else
{
if(self->All.ForceSoftening[pa->Type] < self->All.ForceSoftening[maxsofttype])
diffsoftflag = 1;
}
}
}
}
}
if(mass)
{
s[0] /= mass;
s[1] /= mass;
s[2] /= mass;
vs[0] /= mass;
vs[1] /= mass;
vs[2] /= mass;
}
else
{
s[0] = self->Nodes[no].center[0];
s[1] = self->Nodes[no].center[1];
s[2] = self->Nodes[no].center[2];
}
self->Nodes[no].u.d.s[0] = s[0];
self->Nodes[no].u.d.s[1] = s[1];
self->Nodes[no].u.d.s[2] = s[2];
self->Nodes[no].u.d.mass = mass;
self->Nodes[no].u.d.bitflags = 4 * maxsofttype + 32 * diffsoftflag;
self->Nodes[no].u.d.sibling = sib;
self->Nodes[no].u.d.father = father;
}
else /* single particle */
{
if(self->last >= 0)
{
if(self->last >= self->All.MaxPart)
{
self->Nodes[self->last].u.d.nextnode = no;
}
else
self->Nextnode[self->last] = no;
}
self->last = no;
if(no < self->All.MaxPart) /* only set it for single particles */
self->Father[no] = father;
}
}
/*! Constructs the gravitational oct-tree.
*
* The index convention for accessing tree nodes is the following: the
* indices 0...NumPart-1 reference single particles, the indices
* All.MaxPart.... All.MaxPart+nodes-1 reference tree nodes. `Nodes_base'
* points to the first tree node, while `nodes' is shifted such that
* nodes[All.MaxPart] gives the first tree node. Finally, node indices
* with values 'All.MaxPart + MaxNodes' and larger indicate "pseudo
* particles", i.e. multipole moments of top-level nodes that lie on
* different CPUs. If such a node needs to be opened, the corresponding
* particle must be exported to that CPU. The 'Extnodes' structure
* parallels that of 'Nodes'. Its information is only needed for the SPH
* part of the computation. (The data is split onto these two structures
* as a tuning measure. If it is merged into 'Nodes' a somewhat bigger
* size of the nodes also for gravity would result, which would reduce
* cache utilization slightly.
*/
int force_treebuild_single(Tree *self,int npart)
{
int i, j, subnode = 0, parent, numnodes;
int nfree, th, nn, no;
struct NODE *nfreep;
double lenhalf, epsilon;
peanokey key;
/* create an empty root node */
nfree = self->All.MaxPart; /* index of first free node */
nfreep = &self->Nodes[nfree]; /* select first node */
nfreep->len = self->DomainLen;
for(j = 0; j < 3; j++)
nfreep->center[j] = self->DomainCenter[j];
for(j = 0; j < 8; j++)
nfreep->u.suns[j] = -1;
/* (yr) init list */
for(j = 0; j < 8; j++)
{
nfreep->count[j]=0;
nfreep->first[j]=-1;
nfreep->last[j]=-1;
}
numnodes = 1;
nfreep++;
nfree++;
/* if a high-resolution region in a global tree is used, we need to generate
* an additional set empty nodes to make sure that we have a complete
* top-level tree for the high-resolution inset
*/
nfreep = &self->Nodes[nfree];
parent = -1; /* note: will not be used below before it is changed */
/* now we insert all particles */
for(i = 0; i < npart; i++)
{
//printf("do i=%d\n",i);
/* (yr) init particule */
self->P[i].next = -1;
/* the softening is only used to check whether particles are so close
* that the tree needs not to be refined further
*/
epsilon = self->All.ForceSoftening[self->P[i].Type];
/* start from root */
th = self->All.MaxPart;
while(1)
{
if(th >= self->All.MaxPart) /* we are dealing with an internal node */
{
subnode = 0;
if(self->P[i].Pos[0] > self->Nodes[th].center[0])
subnode += 1;
if(self->P[i].Pos[1] > self->Nodes[th].center[1])
subnode += 2;
if(self->P[i].Pos[2] > self->Nodes[th].center[2])
subnode += 4;
nn = self->Nodes[th].u.suns[subnode];
if(nn >= 0) /* ok, something is in the daughter slot already, need to continue */
{
//printf(" meet a node : i=%d th=%d parent=%d next=%d\n",i,th,parent,nn);
parent = th;
th = nn;
}
else
{
/* here we have found an empty slot where we can attach
* the new particle as a leaf.
*/
//printf(" empty a node : i=%d th=%d parent=%d\n",i,th,parent);
self->Nodes[th].u.suns[subnode] = i;
/* (yr) set parent properties */
if (self->Nodes[th].count[subnode]!=0)
{
printf("we are in trouble here !\n");
endrun(self,9876);
}
//printf(" ++ add : parent node=%d\n",th);
self->Nodes[th].count[subnode] = 1;
self->Nodes[th].first[subnode] = i;
self->Nodes[th].last[subnode] = i;
break; /* done for this particle */
}
}
else
{
/* We try to insert into a leaf with a single particle.*/
/* (yr) check the number of particles in the node */
//printf(" --> a particle i=%d th=%d parent=%d count=%d\n",i,th,parent,self->Nodes[parent].count[subnode]);
if (self->Nodes[parent].count[subnode]>=1)
{
/* the node is full */
/* Need to generate a new internal node at this point.*/
self->Nodes[parent].u.suns[subnode] = nfree;
nfreep->len = 0.5 * self->Nodes[parent].len;
lenhalf = 0.25 * self->Nodes[parent].len;
if(subnode & 1)
nfreep->center[0] = self->Nodes[parent].center[0] + lenhalf;
else
nfreep->center[0] = self->Nodes[parent].center[0] - lenhalf;
if(subnode & 2)
nfreep->center[1] = self->Nodes[parent].center[1] + lenhalf;
else
nfreep->center[1] = self->Nodes[parent].center[1] - lenhalf;
if(subnode & 4)
nfreep->center[2] = self->Nodes[parent].center[2] + lenhalf;
else
nfreep->center[2] = self->Nodes[parent].center[2] - lenhalf;
nfreep->u.suns[0] = -1;
nfreep->u.suns[1] = -1;
nfreep->u.suns[2] = -1;
nfreep->u.suns[3] = -1;
nfreep->u.suns[4] = -1;
nfreep->u.suns[5] = -1;
nfreep->u.suns[6] = -1;
nfreep->u.suns[7] = -1;
/* (yr) init list */
for(j = 0; j < 8; j++)
{
nfreep->count[j]=0;
nfreep->first[j]=-1;
nfreep->last[j]=-1;
}
/* start with the first in the list */
th = self->Nodes[parent].first[subnode];
/* (yr) here, we need to list over the list */
for(j=0;j<self->Nodes[parent].count[subnode];j++) /* here, we should use a while... */
{
subnode = 0;
if(self->P[th].Pos[0] > nfreep->center[0])
subnode += 1;
if(self->P[th].Pos[1] > nfreep->center[1])
subnode += 2;
if(self->P[th].Pos[2] > nfreep->center[2])
subnode += 4;
nfreep->u.suns[subnode] = th;
/* (yr) update parent */
if (nfreep->count[subnode]==0)
nfreep->first[subnode]=th; /* if the particle comes first, ad it on top of the list */
else
self->P[nfreep->first[subnode]].next=th; /* add it to the list */
nfreep->last[subnode]=th; /* put it as last */
nfreep->count[subnode]++; /* this is not a problem as we cannot exceed max part per node */
//printf(" update parent node %d\n",parent);
/* next in the list */
th = self->P[th].next;
}
th = nfree; /* resume trying to insert the new particle at
* the newly created internal node
*/
numnodes++;
nfree++;
nfreep++;
if((numnodes) >= self->MaxNodes)
{
printf("task %d: maximum number %d of tree-nodes reached.\n", self->ThisTask, self->MaxNodes);
printf("for particle %d\n", i);
//dump_particles();
endrun(self,1);
}
}
else
{
/* there is still room for a particle, add it to the list */
/* add it to the list */
self->P[self->Nodes[parent].last[subnode]].next=i;
/* update parent */
self->Nodes[parent].count[subnode]++;
self->Nodes[parent].last[subnode]=i;
}
}
}
}
printf("!!!particles are now inserted numnodes=%d\n",numnodes);
tree_walk_and_gather(self,0);
/* now compute the multipole moments recursively */
self->last = -1;
force_update_node_recursive(self,self->All.MaxPart, -1, -1);
if(self->last >= self->All.MaxPart)
{
if(self->last >= self->All.MaxPart + self->MaxNodes) /* a pseudo-particle */
self->Nextnode[self->last - self->MaxNodes] = -1;
else
self->Nodes[self->last].u.d.nextnode = -1;
}
else
self->Nextnode[self->last] = -1;
tree_walk_and_gather(self,1);
return numnodes;
}
/*! This function is a driver routine for constructing the gravitational
* oct-tree, which is done by calling a small number of other functions.
*/
int force_treebuild(Tree *self,int npart)
{
self->Numnodestree = force_treebuild_single(self,npart);
return self->Numnodestree;
}
/*! This function allocates the memory used for storage of the tree and of
* auxiliary arrays needed for tree-walk and link-lists. Usually,
* maxnodes approximately equal to 0.7*maxpart is sufficient to store the
* tree for up to maxpart particles.
*/
void force_treeallocate(Tree *self,int maxnodes, int maxpart)
{
int i;
size_t bytes;
double allbytes = 0;
double u;
self->MaxNodes = maxnodes;
if(!(self->Nodes_base = malloc(bytes = (self->MaxNodes + 1) * sizeof(struct NODE))))
{
printf("failed to allocate memory for %d tree-nodes (%g MB).\n", self->MaxNodes, bytes / (1024.0 * 1024.0));
endrun(self,3);
}
allbytes += bytes;
self->Nodes = self->Nodes_base - self->All.MaxPart;
if(!(self->Nextnode = malloc(bytes = (maxpart + MAXTOPNODES) * sizeof(int))))
{
printf("Failed to allocate %d spaces for 'Nextnode' array (%g MB)\n", maxpart + MAXTOPNODES,
bytes / (1024.0 * 1024.0));
exit(0);
}
allbytes += bytes;
if(!(self->Father = malloc(bytes = (maxpart) * sizeof(int))))
{
printf("Failed to allocate %d spaces for 'Father' array (%g MB)\n", maxpart, bytes / (1024.0 * 1024.0));
exit(0);
}
allbytes += bytes;
}
/******************************************************************************
TREE OBJECT
*******************************************************************************/
static void
Tree_dealloc(Tree* self)
{
free(self->P);
free(self->TopNodes);
free(self->DomainNodeIndex);
free(self->Nodes_base);
free(self->Nextnode);
free(self->Father);
free(self->Ngblist);
Py_XDECREF(self->first);
Py_XDECREF(self->list);
self->ob_type->tp_free((PyObject*)self);
}
static PyObject *
Tree_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
{
Tree *self;
self = (Tree *)type->tp_alloc(type, 0);
if (self != NULL) {
self->first = PyString_FromString("");
if (self->first == NULL)
{
Py_DECREF(self);
return NULL;
}
self->list = PyList_New(0);
if (self->list == NULL)
{
Py_DECREF(self);
return NULL;
}
self->number = 7;
}
return (PyObject *)self;
}
static int
Tree_init(Tree *self, PyObject *args, PyObject *kwds)
{
import_array();
int i;
PyObject *first=NULL, *tmp, *list=NULL;
PyArrayObject *ntype,*pos,*vel,*mass;
double ErrTolTheta;
static char *kwlist[] = {"first", "npart", "pos","vel", "mass", "ErrTolTheta", NULL};
if (! PyArg_ParseTupleAndKeywords(args, kwds, "|OOOOOd",kwlist,&first,&ntype,&pos,&vel,&mass,&ErrTolTheta))
return -1;
if (first) {
tmp = self->first;
Py_INCREF(first);
self->first = first;
Py_XDECREF(tmp);
}
/* variables related to nbody */
/* here, we should pass direcly the pointer */
if ((ntype->nd != 1) || (ntype->dimensions[0]!=6))
{
PyErr_SetString(PyExc_ValueError,"Tree_init, npart must be an array of dimension 1x6.");
return -1;
}
self->NtypeLocal[0] = *(int*) (ntype->data + 0*(ntype->strides[0]));
self->NtypeLocal[1] = *(int*) (ntype->data + 1*(ntype->strides[0]));
self->NtypeLocal[2] = *(int*) (ntype->data + 2*(ntype->strides[0]));
self->NtypeLocal[3] = *(int*) (ntype->data + 3*(ntype->strides[0]));
self->NtypeLocal[4] = *(int*) (ntype->data + 4*(ntype->strides[0]));
self->NtypeLocal[5] = *(int*) (ntype->data + 5*(ntype->strides[0]));
self->NumPart = 0;
self->N_gas = self->NtypeLocal[0];
for (i = 0; i < 6; i++)
self->NumPart += self->NtypeLocal[i];
self->All.TotNumPart = self->NumPart;
self->All.TotN_gas = self->N_gas;
/* global variables */
self->ThisTask = 0;
self->NTask = 1;
/* All vars */
for (i = 0; i < 6; i++)
self->All.SofteningTable[i] = 0.1; /* a changer !!!! */
for (i = 0; i < 6; i++)
self->All.ForceSoftening[i] = 0.1; /* a changer !!!! */
self->All.PartAllocFactor = 1.5;
self->All.TreeAllocFactor = 4.0;
self->All.ErrTolTheta = ErrTolTheta;
self->All.MaxPart = self->All.PartAllocFactor * (self->All.TotNumPart / self->NTask);
self->All.MaxPartSph = self->All.PartAllocFactor * (self->All.TotN_gas / self->NTask);
self->All.DesNumNgb = 33;
self->All.MaxNumNgbDeviation = 3;
self->All.MinGasHsmlFractional = 0.25;
self->All.MinGasHsml = self->All.MinGasHsmlFractional * self->All.ForceSoftening[0];
/* create P */
size_t bytes;
if(!(self->P = malloc(bytes = self->All.MaxPart * sizeof(struct particle_data))))
{
printf("failed to allocate memory for `P' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(self,1);
}
for (i = 0; i < pos->dimensions[0]; i++)
{
self->P[i].Pos[0] = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
self->P[i].Pos[1] = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
self->P[i].Pos[2] = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
self->P[i].Vel[0] = *(float *) (vel->data + i*(vel->strides[0]) + 0*vel->strides[1]);
self->P[i].Vel[1] = *(float *) (vel->data + i*(vel->strides[0]) + 1*vel->strides[1]);
self->P[i].Vel[2] = *(float *) (vel->data + i*(vel->strides[0]) + 2*vel->strides[1]);
self->P[i].Mass = *(float *) (mass->data + i*(mass->strides[0]));
self->P[i].Type = 0; /* this should be changed... */
}
/***************************************
* find domain extent *
/***************************************/
domain_findExtent(self);
/************************
/* tree construction *
/************************/
force_treeallocate(self,self->All.TreeAllocFactor * self->All.MaxPart, self->All.MaxPart);
force_treebuild(self,self->NumPart);
return 0;
}
static PyMemberDef Tree_members[] = {
{"first", T_OBJECT_EX, offsetof(Tree, first), 0,
"first name"},
{"list", T_OBJECT_EX, offsetof(Tree, list), 0,
"list of"},
{"number", T_INT, offsetof(Tree, number), 0,
"Tree number"},
{NULL} /* Sentinel */
};
static PyObject *
Tree_name(Tree* self)
{
static PyObject *format = NULL;
PyObject *args, *result;
if (format == NULL) {
format = PyString_FromString("%s %s");
if (format == NULL)
return NULL;
}
if (self->first == NULL) {
PyErr_SetString(PyExc_AttributeError, "first");
return NULL;
}
result = PyString_Format(format, args);
Py_DECREF(args);
return result;
}
static PyObject *
Tree_info(Tree* self)
{
//static PyObject *format = NULL;
//PyObject *args, *result;
printf("NumPart = %d\n",self->NumPart);
printf("N_gas = %d\n",self->N_gas);
printf("DomainLen = %g\n",self->DomainLen);
printf("DomainCenter x = %g\n",self->DomainCenter[0]);
printf("DomainCenter y = %g\n",self->DomainCenter[1]);
printf("DomainCenter z = %g\n",self->DomainCenter[2]);
printf("DomainCorner x = %g\n",self->DomainCorner[0]);
printf("DomainCorner y = %g\n",self->DomainCorner[1]);
printf("DomainCorner z = %g\n",self->DomainCorner[2]);
return Py_BuildValue("i",1);
}
static PyObject *
Tree_Acceleration(Tree* self, PyObject *args)
{
PyArrayObject *pos;
float eps;
if (! PyArg_ParseTuple(args, "Of",&pos,&eps))
return PyString_FromString("error");
PyArrayObject *acc;
int i;
double x,y,z,ax,ay,az;
int input_dimension;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
if (pos->descr->type_num != PyArray_FLOAT)
PyErr_SetString(PyExc_ValueError,"argument 1 must be of type Float32");
/* create a NumPy object */
npy_intp ld[2];
ld[0]=pos->dimensions[0];
ld[1]=pos->dimensions[1];
/* there is a kind of bug here ! I cannt replace ld by pos->dimensions */
//acc = (PyArrayObject *) PyArray_FromDims(pos->nd,ld,pos->descr->type_num);
acc = (PyArrayObject *) PyArray_SimpleNew(pos->nd,ld,pos->descr->type_num);
for (i = 0; i < pos->dimensions[0]; i++)
{
x = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
y = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
z = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
force_treeevaluate_local(self,x,y,z,(double)eps,&ax,&ay,&az);
*(float *)(acc->data + i*(acc->strides[0]) + 0*acc->strides[1]) = ax;
*(float *)(acc->data + i*(acc->strides[0]) + 1*acc->strides[1]) = ay;
*(float *)(acc->data + i*(acc->strides[0]) + 2*acc->strides[1]) = az;
}
return PyArray_Return(acc);
}
static PyObject *
Tree_Potential(Tree* self, PyObject *args)
{
PyArrayObject *pos;
float eps;
if (! PyArg_ParseTuple(args, "Of",&pos,&eps))
return PyString_FromString("error");
PyArrayObject *pot;
int i;
double x,y,z,lpot;
npy_intp ld[1];
int input_dimension;
input_dimension =pos->nd;
if (input_dimension != 2)
PyErr_SetString(PyExc_ValueError,"dimension of first argument must be 2");
if (pos->descr->type_num != PyArray_FLOAT)
PyErr_SetString(PyExc_ValueError,"argument 1 must be of type Float32");
/* create a NumPy object */
ld[0]=pos->dimensions[0];
//pot = (PyArrayObject *) PyArray_FromDims(1,ld,PyArray_FLOAT);
pot = (PyArrayObject *) PyArray_SimpleNew(1,ld,NPY_FLOAT);
for (i = 0; i < pos->dimensions[0]; i++)
{
x = *(float *) (pos->data + i*(pos->strides[0]) + 0*pos->strides[1]);
y = *(float *) (pos->data + i*(pos->strides[0]) + 1*pos->strides[1]);
z = *(float *) (pos->data + i*(pos->strides[0]) + 2*pos->strides[1]);
lpot = force_treeevaluate_local_potential(self,x,y,z,(double)eps);
*(float *)(pot->data + i*(pot->strides[0])) = lpot;
}
return PyArray_Return(pot);
}
static PyMethodDef Tree_methods[] = {
{"name", (PyCFunction)Tree_name, METH_NOARGS,
"Return the name, combining the first and last name"
},
{"info", (PyCFunction)Tree_info, METH_NOARGS,
"Return info"
},
{"Potential", (PyCFunction)Tree_Potential, METH_VARARGS,
"This function computes the potential at a given position using the tree"
},
{"Acceleration", (PyCFunction)Tree_Acceleration, METH_VARARGS,
"This function computes the acceleration at a given position using the tree"
},
{NULL} /* Sentinel */
};
static PyTypeObject TreeType = {
PyObject_HEAD_INIT(NULL)
0, /*ob_size*/
"tree.Tree", /*tp_name*/
sizeof(Tree), /*tp_basicsize*/
0, /*tp_itemsize*/
(destructor)Tree_dealloc, /*tp_dealloc*/
0, /*tp_print*/
0, /*tp_getattr*/
0, /*tp_setattr*/
0, /*tp_compare*/
0, /*tp_repr*/
0, /*tp_as_number*/
0, /*tp_as_sequence*/
0, /*tp_as_mapping*/
0, /*tp_hash */
0, /*tp_call*/
0, /*tp_str*/
0, /*tp_getattro*/
0, /*tp_setattro*/
0, /*tp_as_buffer*/
Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
"Tree objects", /* tp_doc */
0, /* tp_traverse */
0, /* tp_clear */
0, /* tp_richcompare */
0, /* tp_weaklistoffset */
0, /* tp_iter */
0, /* tp_iternext */
Tree_methods, /* tp_methods */
Tree_members, /* tp_members */
0, /* tp_getset */
0, /* tp_base */
0, /* tp_dict */
0, /* tp_descr_get */
0, /* tp_descr_set */
0, /* tp_dictoffset */
(initproc)Tree_init, /* tp_init */
0, /* tp_alloc */
Tree_new, /* tp_new */
};
static PyMethodDef module_methods[] = {
{NULL} /* Sentinel */
};
#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */
#define PyMODINIT_FUNC void
#endif
PyMODINIT_FUNC
inittreelib(void)
{
PyObject* m;
if (PyType_Ready(&TreeType) < 0)
return;
m = Py_InitModule3("treelib", module_methods,
"Example module that creates an extension type.");
if (m == NULL)
return;
Py_INCREF(&TreeType);
PyModule_AddObject(m, "Tree", (PyObject *)&TreeType);
}
c4science · Help