Homec4science
Diffusion pNbody (master) Browse

pNbody/src/streelib/treelib.000.c
ea88594ae72emaster

Last Change
Raw File

treelib.000.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;
};
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
*******************************************************************************/
#define REDUC_FAC 0.98
#define TOPNODEFACTOR 20.0
/*! This is a comparison kernel used in a sort routine.
*/
int domain_compare_toplist(const void *a, const void *b)
{
if(((struct topnode_exchange *) a)->Startkey < (((struct topnode_exchange *) b)->Startkey))
return -1;
if(((struct topnode_exchange *) a)->Startkey > (((struct topnode_exchange *) b)->Startkey))
return +1;
return 0;
}
/*! This is a comparison kernel used in a sort routine.
*/
int domain_compare_key(const void *a, const void *b)
{
if(*(peanokey *) a < *(peanokey *) b)
return -1;
if(*(peanokey *) a > *(peanokey *) b)
return +1;
return 0;
}
/*! This function walks the global top tree in order to establish the
* number of leaves it has. These leaves are distributed to different
* processors.
*/
void domain_walktoptree(Tree *self,int no)
{
int i;
if(self->TopNodes[no].Daughter == -1)
{
self->TopNodes[no].Leaf = self->NTopleaves;
self->NTopleaves++;
}
else
{
for(i = 0; i < 8; i++)
domain_walktoptree(self,self->TopNodes[no].Daughter + i);
}
}
/*! This function is responsible for constructing the local top-level
* Peano-Hilbert segments. A segment is cut into 8 pieces recursively
* until the number of particles in the segment has fallen below
* All.TotNumPart / (TOPNODEFACTOR * NTask * NTask).
*/
void domain_topsplit_local(Tree *self,int node, peanokey startkey)
{
int i, p, sub, bin;
if(self->TopNodes[node].Size >= 8)
{
self->TopNodes[node].Daughter = self->NTopnodes;
for(i = 0; i < 8; i++)
{
if(self->NTopnodes < MAXTOPNODES)
{
sub = self->TopNodes[node].Daughter + i;
self->TopNodes[sub].Size = self->TopNodes[node].Size / 8;
self->TopNodes[sub].Count = 0;
self->TopNodes[sub].Daughter = -1;
self->TopNodes[sub].StartKey = startkey + i * self->TopNodes[sub].Size;
self->TopNodes[sub].Pstart = self->TopNodes[node].Pstart;
self->NTopnodes++;
}
else
{
printf("task=%d: We are out of Topnodes. Increasing the constant MAXTOPNODES might help.\n",self->ThisTask);
fflush(stdout);
endrun(self,13213);
}
}
for(p = self->TopNodes[node].Pstart; p < self->TopNodes[node].Pstart + self->TopNodes[node].Count; p++)
{
bin = (self->KeySorted[p] - startkey) / (self->TopNodes[node].Size / 8);
if(bin < 0 || bin > 7)
{
printf("task=%d: something odd has happened here. bin=%d\n", self->ThisTask, bin);
fflush(stdout);
endrun(self,13123123);
}
sub = self->TopNodes[node].Daughter + bin;
if(self->TopNodes[sub].Count == 0)
self->TopNodes[sub].Pstart = p;
self->TopNodes[sub].Count++;
}
for(i = 0; i < 8; i++)
{
sub = self->TopNodes[node].Daughter + i;
if(self->TopNodes[sub].Count > self->All.TotNumPart / (TOPNODEFACTOR * self->NTask * self->NTask))
domain_topsplit_local(self,sub, self->TopNodes[sub].StartKey);
}
}
}
/*! This function constructs the global top-level tree node that is used
* for the domain decomposition. This is done by considering the string of
* Peano-Hilbert keys for all particles, which is recursively chopped off
* in pieces of eight segments until each segment holds at most a certain
* number of particles.
*/
void domain_determineTopTree(Tree *self)
{
int i, ntop_local, ntop;
int *ntopnodelist, *ntopoffset;
for(i = 0; i < self->NumPart; i++)
{
self->KeySorted[i] = self->Key[i] = peano_hilbert_key((self->P[i].Pos[0] - self->DomainCorner[0]) * self->DomainFac,
(self->P[i].Pos[1] - self->DomainCorner[1]) * self->DomainFac,
(self->P[i].Pos[2] - self->DomainCorner[2]) * self->DomainFac,
BITS_PER_DIMENSION);
}
qsort(self->KeySorted, self->NumPart, sizeof(peanokey), domain_compare_key);
self->NTopnodes = 1;
self->TopNodes[0].Daughter = -1;
self->TopNodes[0].Size = PEANOCELLS;
self->TopNodes[0].StartKey = 0;
self->TopNodes[0].Count = self->NumPart;
self->TopNodes[0].Pstart = 0;
domain_topsplit_local(self,0, 0);
/* the rest of the function is only used for parallelisme */
}
/*! 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));
}
/*! This is the main routine for the domain decomposition. It acts as a
* driver routine that allocates various temporary buffers, maps the
* particles back onto the periodic box if needed, and then does the
* domain decomposition, and a final Peano-Hilbert order of all particles
* as a tuning measure.
*/
void domain_Decomposition(Tree *self)
{
double t0, t1;
self->Key = malloc(sizeof(peanokey) * self->All.MaxPart);
self->KeySorted = malloc(sizeof(peanokey) * self->All.MaxPart);
domain_findExtent(self);
domain_determineTopTree(self);
free(self->KeySorted);
free(self->Key);
}
void domain_domainallocate(Tree *self)
{
self->TopNodes = malloc(MAXTOPNODES * sizeof(struct topnode_data));
self->DomainNodeIndex = malloc(MAXTOPNODES * sizeof(int));
}
/******************************************************************************
FORCE THINGS
*******************************************************************************/
struct NODE
{
FLOAT len; /*!< sidelength of treenode */
FLOAT center[3]; /*!< geometrical center of 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
*/
void tree_walk_recursive(Tree *self,int no, int sib, int father)
{
int j,p;
int nextsib;
nNodes++;
if(no >= self->All.MaxPart && no < self->All.MaxPart + self->MaxNodes)
{
/* 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");
}
}
/*
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 no;
int sib=-1;
int father=-1;
no = self->All.MaxPart; /* root node */
printf("start tree_walk_and_gather...\n");
tree_walk_recursive(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;
//#ifndef UNEQUALSOFTENINGS
// h = All.ForceSoftening[ptype];
// h_inv = 1.0 / h;
// h3_inv = h_inv * h_inv * h_inv;
//#endif
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
{
if(no >= self->All.MaxPart + self->MaxNodes) /* pseudo particle */
{
//if(mode == 0)
//{
// Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
//}
//no = Nextnode[no - MaxNodes];
//continue;
}
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;
}
#ifdef PERIODIC
dx = NEAREST(dx);
dy = NEAREST(dy);
dz = NEAREST(dz);
#endif
r2 = dx * dx + dy * dy + dz * dz;
if(no < self->All.MaxPart)
{
//#ifdef UNEQUALSOFTENINGS
//#ifdef ADAPTIVE_GRAVSOFT_FORGAS
// if(ptype == 0)
// h = soft;
// else
// h = All.ForceSoftening[ptype];
//
// if(P[no].Type == 0)
// {
// if(h < SphP[no].Hsml)
// h = SphP[no].Hsml;
// }
// else
// {
// if(h < All.ForceSoftening[P[no].Type])
// h = All.ForceSoftening[P[no].Type];
// }
//#else
// h = All.ForceSoftening[ptype];
// if(h < All.ForceSoftening[P[no].Type])
// h = All.ForceSoftening[P[no].Type];
//#endif
//#endif
no = self->Nextnode[no];
}
else /* we have an internal node. Need to check opening criterion */
{
// if(mode == 1)
// {
// if((nop->u.d.bitflags & 3) == 1) /* if it's a top-level node
// * which does not contain
// * local particles we can
// * continue to do a short-cut */
// {
// no = nop->u.d.sibling;
// continue;
// }
// }
// 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;
}
}
// else /* check relative opening criterion */
// {
// if(mass * nop->len * nop->len > r2 * r2 * aold)
// {
// /* open cell */
// no = nop->u.d.nextnode;
// continue;
// }
//
// /* check in addition whether we lie inside the cell */
//
// if(fabs(nop->center[0] - pos_x) < 0.60 * nop->len)
// {
// if(fabs(nop->center[1] - pos_y) < 0.60 * nop->len)
// {
// if(fabs(nop->center[2] - pos_z) < 0.60 * nop->len)
// {
// no = nop->u.d.nextnode;
// continue;
// }
// }
// }
// }
//#ifdef UNEQUALSOFTENINGS
//#ifndef ADAPTIVE_GRAVSOFT_FORGAS
// h = All.ForceSoftening[ptype];
// maxsofttype = (nop->u.d.bitflags >> 2) & 7;
// if(maxsofttype == 7) /* may only occur for zero mass top-level nodes */
// {
// if(mass > 0)
// endrun(986);
// no = nop->u.d.nextnode;
// continue;
// }
// else
// {
// if(h < All.ForceSoftening[maxsofttype])
// {
// h = All.ForceSoftening[maxsofttype];
// if(r2 < h * h)
// {
// if(((nop->u.d.bitflags >> 5) & 1)) /* bit-5 signals that there are particles of different softening in the node */
// {
// no = nop->u.d.nextnode;
// continue;
// }
// }
// }
// }
//#else
// if(ptype == 0)
// h = soft;
// else
// h = All.ForceSoftening[ptype];
//
// if(h < nop->maxsoft)
// {
// h = nop->maxsoft;
// if(r2 < h * h)
// {
// no = nop->u.d.nextnode;
// continue;
// }
// }
//#endif
//#endif
no = nop->u.d.sibling; /* ok, node can be used */
// if(mode == 1)
// {
// if(((nop->u.d.bitflags) & 1)) /* Bit 0 signals that this node belongs to top-level tree */
// continue;
// }
}
// r = sqrt(r2);
//
// if(r >= h)
// fac = mass / (r2 * r);
// else
// {
////#ifdef UNEQUALSOFTENINGS
//// h_inv = 1.0 / h;
//// h3_inv = h_inv * h_inv * h_inv;
////#endif
// u = r * h_inv;
// if(u < 0.5)
// fac = mass * h3_inv * (10.666666666667 + u * u * (32.0 * u - 38.4));
// else
// fac =
// mass * h3_inv * (21.333333333333 - 48.0 * u +
// 38.4 * u * u - 10.666666666667 * u * u * u - 0.066666666667 / (u * u * u));
// }
/* 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++;
}
#ifdef PERIODIC
*ewaldcountsum += force_treeevaluate_ewald_correction(target, mode, pos_x, pos_y, pos_z, aold);
#endif
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
{
if(no >= self->All.MaxPart + self->MaxNodes) /* pseudo particle */
{
//if(mode == 0)
//{
// Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
//}
//no = Nextnode[no - MaxNodes];
//continue;
printf("force_treeevaluate_local_potential : pseudo particle !\n");
endrun(self,1234);
}
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;
}
#ifdef PERIODIC
dx = NEAREST(dx);
dy = NEAREST(dy);
dz = NEAREST(dz);
#endif
r2 = dx * dx + dy * dy + dz * dz;
if(no < self->All.MaxPart) /* single particle */
{
//#ifdef UNEQUALSOFTENINGS
//#ifdef ADAPTIVE_GRAVSOFT_FORGAS
// if(ptype == 0)
// h = soft;
// else
// h = All.ForceSoftening[ptype];
//
// if(P[no].Type == 0)
// {
// if(h < SphP[no].Hsml)
// h = SphP[no].Hsml;
// }
// else
// {
// if(h < All.ForceSoftening[P[no].Type])
// h = All.ForceSoftening[P[no].Type];
// }
//#else
// h = All.ForceSoftening[ptype];
// if(h < All.ForceSoftening[P[no].Type])
// h = All.ForceSoftening[P[no].Type];
//#endif
//#endif
no = self->Nextnode[no];
}
else /* we have an internal node. Need to check opening criterion */
{
// if(mode == 1)
// {
// if((nop->u.d.bitflags & 3) == 1) /* if it's a top-level node
// * which does not contain
// * local particles we can make
// * a short-cut
// */
// {
// no = nop->u.d.sibling;
// continue;
// }
// }
// 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;
}
}
// else /* check relative opening criterion */
// {
// if(mass * nop->len * nop->len > r2 * r2 * aold)
// {
// /* open cell */
// no = nop->u.d.nextnode;
// continue;
// }
//
// if(fabs(nop->center[0] - pos_x) < 0.60 * nop->len)
// {
// if(fabs(nop->center[1] - pos_y) < 0.60 * nop->len)
// {
// if(fabs(nop->center[2] - pos_z) < 0.60 * nop->len)
// {
// no = nop->u.d.nextnode;
// continue;
// }
// }
// }
// }
//#ifdef UNEQUALSOFTENINGS
//#ifndef ADAPTIVE_GRAVSOFT_FORGAS
// h = All.ForceSoftening[ptype];
// maxsofttype = (nop->u.d.bitflags >> 2) & 7;
// if(maxsofttype == 7) /* may only occur for zero mass top-level nodes */
// {
// if(mass > 0)
// endrun(988);
// no = nop->u.d.nextnode;
// continue;
// }
// else
// {
// if(h < All.ForceSoftening[maxsofttype])
// {
// h = All.ForceSoftening[maxsofttype];
// if(r2 < h * h)
// {
// if(((nop->u.d.bitflags >> 5) & 1)) /* bit-5 signals that there are particles of different softening in the node */
// {
// no = nop->u.d.nextnode;
// continue;
// }
// }
// }
// }
//#else
// if(ptype == 0)
// h = soft;
// else
// h = All.ForceSoftening[ptype];
//
// if(h < nop->maxsoft)
// {
// h = nop->maxsoft;
// if(r2 < h * h)
// {
// no = nop->u.d.nextnode;
// continue;
// }
// }
//#endif
//#endif
no = nop->u.d.sibling; /* node can be used */
// if(mode == 1)
// {
// if(((nop->u.d.bitflags) & 1)) /* Bit 0 signals that this node belongs to top-level tree */
// continue;
// }
}
// r = sqrt(r2);
//
// if(r >= h)
// pot -= mass / r;
// else
// {
////#ifdef UNEQUALSOFTENINGS
//// h_inv = 1.0 / h;
////#endif
// u = r * h_inv;
//
// if(u < 0.5)
// wp = -2.8 + u * u * (5.333333333333 + u * u * (6.4 * u - 9.6));
// else
// wp =
// -3.2 + 0.066666666667 / u + u * u * (10.666666666667 +
// u * (-16.0 + u * (9.6 - 2.133333333333 * u)));
//
// pot += mass * h_inv * wp;
//
// }
/* here we use a plummer softening */
if (r2>0)
pot += -mass / sqrt(r2+h*h);
#ifdef PERIODIC
pot += mass * ewald_pot_corr(dx, dy, dz);
#endif
}
/* return result */
return pot;
}
/*! This function flags nodes in the top-level tree that are dependent on
* local particle data.
*/
void force_flag_localnodes(Tree *self)
{
int no, i;
/* mark all top-level nodes */
for(i = 0; i < self->NTopleaves; i++)
{
no = self->DomainNodeIndex[i];
while(no >= 0)
{
if((self->Nodes[no].u.d.bitflags & 1))
break;
self->Nodes[no].u.d.bitflags |= 1;
no = self->Nodes[no].u.d.father;
}
}
/* mark top-level nodes that contain local particles */
//for(i = DomainMyStart; i <= DomainMyLast; i++) /* !!! warning !!! */
for(i = 0; i < self->NTopleaves; i++)
{
/*
if(DomainMoment[i].mass > 0)
*/
{
no = self->DomainNodeIndex[i];
while(no >= 0)
{
if((self->Nodes[no].u.d.bitflags & 2))
break;
self->Nodes[no].u.d.bitflags |= 2;
no = self->Nodes[no].u.d.father;
}
}
}
}
/*! 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)
{
if(self->last >= self->All.MaxPart + self->MaxNodes) /* a pseudo-particle */
self->Nextnode[self->last - self->MaxNodes] = no;
else
self->Nodes[self->last].u.d.nextnode = no;
}
else
self->Nextnode[self->last] = no;
}
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++)
{
if((p = suns[j]) >= 0)
{
/* 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 */
{
if(p >= self->All.MaxPart + self->MaxNodes) /* a pseudo particle */
{
/* nothing to be done here because the mass of the
* pseudo-particle is still zero. This will be changed
* later.
*/
}
else
{
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];
//vs[0] += self->Nodes[p].u.d.mass * self->Extnodes[p].vs[0];
//vs[1] += self->Nodes[p].u.d.mass * self->Extnodes[p].vs[1];
//vs[2] += self->Nodes[p].u.d.mass * self->Extnodes[p].vs[2];
//if(Extnodes[p].hmax > hmax)
// hmax = Extnodes[p].hmax;
//#ifdef UNEQUALSOFTENINGS
//#ifndef ADAPTIVE_GRAVSOFT_FORGAS
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
// if(self->Nodes[p].maxsoft > maxsoft)
// maxsoft = Nodes[p].maxsoft;
//#endif
//#endif
}
}
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];
//#ifdef UNEQUALSOFTENINGS
//#ifndef ADAPTIVE_GRAVSOFT_FORGAS
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;
}
}
//#else
// if(pa->Type == 0)
// {
// if(SphP[p].Hsml > maxsoft)
// maxsoft = SphP[p].Hsml;
// }
// else
// {
// if(All.ForceSoftening[pa->Type] > maxsoft)
// maxsoft = All.ForceSoftening[pa->Type];
// }
//#endif
//#endif
// if(pa->Type == 0)
// if(SphP[p].Hsml > hmax)
// hmax = SphP[p].Hsml;
}
}
}
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;
//#ifdef UNEQUALSOFTENINGS
//#ifndef ADAPTIVE_GRAVSOFT_FORGAS
self->Nodes[no].u.d.bitflags = 4 * maxsofttype + 32 * diffsoftflag;
//#else
// self->Nodes[no].u.d.bitflags = 0;
// self->Nodes[no].maxsoft = maxsoft;
//#endif
//#else
// Nodes[no].u.d.bitflags = 0;
//#endif
//Extnodes[no].vs[0] = vs[0];
//Extnodes[no].vs[1] = vs[1];
//Extnodes[no].vs[2] = vs[2];
//Extnodes[no].hmax = hmax;
self->Nodes[no].u.d.sibling = sib;
self->Nodes[no].u.d.father = father;
}
else /* single particle or pseudo particle */
{
if(self->last >= 0)
{
if(self->last >= self->All.MaxPart)
{
if(self->last >= self->All.MaxPart + self->MaxNodes) /* a pseudo-particle */
self->Nextnode[self->last - self->MaxNodes] = no;
else
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;
}
}
/*! This function recursively creates a set of empty tree nodes which
* corresponds to the top-level tree for the domain grid. This is done to
* ensure that this top-level tree is always "complete" so that we can
* easily associate the pseudo-particles of other CPUs with tree-nodes at
* a given level in the tree, even when the particle population is so
* sparse that some of these nodes are actually empty.
*/
void force_create_empty_nodes(Tree *self,int no, int topnode, int bits, int x, int y, int z, int *nodecount,
int *nextfree)
{
int i, j, k, n, sub, count;
if(self->TopNodes[topnode].Daughter >= 0)
{
for(i = 0; i < 2; i++)
for(j = 0; j < 2; j++)
for(k = 0; k < 2; k++)
{
/* gives values betwen 0 and 7 : index of daughter TopNodes */
/* self->TopNodes[topnode].Daughter + sub : index to daughter TopNode */
sub = 7 & peano_hilbert_key((x << 1) + i, (y << 1) + j, (z << 1) + k, bits);
count = i + 2 * j + 4 * k;
self->Nodes[no].u.suns[count] = *nextfree;
/* set size and center */
self->Nodes[*nextfree].len = 0.5 * self->Nodes[no].len;
self->Nodes[*nextfree].center[0] = self->Nodes[no].center[0] + (2 * i - 1) * 0.25 * self->Nodes[no].len;
self->Nodes[*nextfree].center[1] = self->Nodes[no].center[1] + (2 * j - 1) * 0.25 * self->Nodes[no].len;
self->Nodes[*nextfree].center[2] = self->Nodes[no].center[2] + (2 * k - 1) * 0.25 * self->Nodes[no].len;
/* all daughters are leaves (empty) */
for(n = 0; n < 8; n++)
self->Nodes[*nextfree].u.suns[n] = -1;
/* set the correspondance between the Leaves of TopNodes and Nodes of the tree */
if(self->TopNodes[self->TopNodes[topnode].Daughter + sub].Daughter == -1)
self->DomainNodeIndex[self->TopNodes[self->TopNodes[topnode].Daughter + sub].Leaf] = *nextfree;
*nextfree = *nextfree + 1;
*nodecount = *nodecount + 1;
if((*nodecount) >= self->MaxNodes)
{
printf("task %d: maximum number %d of tree-nodes reached.\n", self->ThisTask, self->MaxNodes);
printf("in create empty nodes\n");
//dump_particles();
endrun(self,11);
}
force_create_empty_nodes(self,*nextfree - 1, self->TopNodes[topnode].Daughter + sub,
bits + 1, 2 * x + i, 2 * y + j, 2 * z + k, nodecount, nextfree);
}
}
}
/*! 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;
numnodes = 1;
nfreep++;
nfree++;
/* create a set of empty nodes corresponding to the top-level domain
* grid. We need to generate these nodes first to make sure that we have a
* complete top-level tree which allows the easy insertion of the
* pseudo-particles at the right place
*/
//force_create_empty_nodes(self,self->All.MaxPart, 0, 1, 0, 0, 0, &numnodes, &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++)
{
/* 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];
/* si on utilise TopTree
key = peano_hilbert_key((self->P[i].Pos[0] - self->DomainCorner[0]) * self->DomainFac,
(self->P[i].Pos[1] - self->DomainCorner[1]) * self->DomainFac,
(self->P[i].Pos[2] - self->DomainCorner[2]) * self->DomainFac, BITS_PER_DIMENSION);
no = 0;
while(self->TopNodes[no].Daughter >= 0)
no = self->TopNodes[no].Daughter + (key - self->TopNodes[no].StartKey) / (self->TopNodes[no].Size / 8); // trouve la feuille auquel la particule apartiend...
no = self->TopNodes[no].Leaf;
th = self->DomainNodeIndex[no];
*/
/* we always start from the first cell */
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 */
{
parent = th;
th = nn;
}
else
{
/* here we have found an empty slot where we can attach
* the new particle as a leaf.
*/
self->Nodes[th].u.suns[subnode] = i;
break; /* done for this particle */
}
}
else
{
/* We try to insert into a leaf with a single particle. 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;
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;
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);
}
}
}
}
/* insert the pseudo particles that represent the mass distribution of other domains */
//force_insert_pseudo_particles();
/* 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;
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);
force_flag_localnodes(self);
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... */
}
/***************************************
* domain decomposition construction *
/***************************************/
domain_domainallocate(self);
domain_Decomposition(self);
self->NTopleaves = 0; /* normally in domain_sumCost */
domain_walktoptree(self,0); /* normally in domain_sumCost */
printf("NTopnodes = %d\n",self->NTopnodes);
printf("NTopleaves = %d\n",self->NTopleaves);
/* count the number of particles per topnodes */
/*
int npart;
peanokey key;
float epsilon;
int no;
npart = self->NumPart;
for(i = 0; i < npart; i++)
{
epsilon = self->All.ForceSoftening[self->P[i].Type];
key = peano_hilbert_key((self->P[i].Pos[0] - self->DomainCorner[0]) * self->DomainFac,
(self->P[i].Pos[1] - self->DomainCorner[1]) * self->DomainFac,
(self->P[i].Pos[2] - self->DomainCorner[2]) * self->DomainFac, BITS_PER_DIMENSION);
no = 0;
while(self->TopNodes[no].Daughter >= 0)
//trouve la feuille auquel la particule apartiend...
no = self->TopNodes[no].Daughter + (key - self->TopNodes[no].StartKey) / (self->TopNodes[no].Size / 8);
printf("%d %d\n",i,no);
}
/************************
/* 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