diff --git a/src/streelib/Readme b/src/streelib/Readme index a240546..be665bf 100644 --- a/src/streelib/Readme +++ b/src/streelib/Readme @@ -1,294 +1,323 @@ treelib.000.c : version avec peut-etre encore un peu de TopNodes treelib.001.c : version épurée, plus de topnodes plus de pseudo particles pas d'unequal softenning pas de periodic box pas de gravite adaptatif treelib.c : on essaye d'ajouter le nombre de particules permises par noeuds see end of the file... but : avoir un tree minimal domain_Decomposition domain_findExtent ok domain_determineTopTree self->TopNodes[0] domain_topsplit_local domain_walktoptree force_treeallocate self->Nodes force_treebuild force_treebuild_single *force_create_empty_nodes (no longer necessary) * -- insert particles -- self->last = -1; /* pointe sur le dernier noeud fait */ *force_update_node_recursive TOUT OK implementation de la suite: -------------------------- 1) chaque cellule enregistre ca quantité de particules --> on peut crée une routine pour vérifier... la suite : -------- - regarder un peu force_tree - stoquer plusieurs particules par noeuds... - construction --> ok, on peut utiliser Nextnode pour cela a) on pointe sur Node => fin de la liste b) on pointe sur P => suite de la liste -force_tree --> si un noeud a plusieurs particules => il n'a pas de descendants => on calcul la somme de toutes les particules présentes --> ok, facile... - revoir la recherche des voisins - enlever les branches vides... ??? - sur GPU : --> on definit des cellules contenant des particules --> avec GPU : pour chacune de ces cellules, c'est facile de calculer l'accel de toutes sur toutes... >>> chaque particules a les forces correspondant a ces voisins de cellules a) --> on calcul par la suite les autres contribution avec un arbre normal b) --> lorsque une cellule doit etre ouverte (et ne contiend que des particules) l'interaction peut se faire sur gpu également ##################################### how to remove topNodes ? ---> ok,done :-) ....> j'ai un treecode minmal... ##################################### self->Father : non utilisé : utique que lorsque le tree est updated dynamically self->Nextnode ???? utilité ? * initialllisé dans : force_update_node_recursive * utilisé lors du tree walking, que pour les particules et psudo-particules ###################################################### # ###################################################### -> arbre sur les particules triée avec peano -> pas d'info sur la position... (si ce n'est peanokey StartKey) TopNodes[no].Daughter == -1 --> pas cde cellules filles 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 */ }; si une cellule contiend plus que self->All.TotNumPart / (TOPNODEFACTOR * self->NTask * self->NTask) (dans topsplit local) particules, on l'ouvre. ###################################################### extern int *DomainNodeIndex; /*!< this table gives for each leaf of the top-level tree the corresponding node of the gravitational tree */ ???? quand est que DomainNodeIndex est remplis ? ###################################################### # count the number of part in a node ###################################################### struct 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*/ !!! this requires a lot of memory !!! !!! should do this only if a daugher contains a list of particle !!! !!! .first[subnode] should be equal to .u.suns[subnode] !!! count is not necessary, can check if first=last struct particle_data int next; actuellement: ------------ th est un index si th pointe sur un noeud on regarde ou la particule va -> subnode -> nn = self->Nodes[th].u.suns[subnode]; nn>=0 le noeud es occupé parent=th th=nn --> on continue nn = -1 le noeud est vide --> on y attache la particule self->Nodes[th].u.suns[subnode] sinon la cellule contient une particule il faut generer des cellules filles 1) on lie la cellule parent a la cellule fille=next free on pousse la particule qui réside déjà dans une des nouvelle cellule th = nfree : on essaye d'incerer la particule dans la nouvelle cellule Thu Oct 28 23:31:53 CEST 2010 need to test if (self->Nodes[parent].count[subnode]>=1) --> if (self->Nodes[parent].count[subnode]>=2) Fri Oct 29 21:30:32 CEST 2010 ok it seems to work -> now, we need a tool to explore the tree in python * draw the nodes * draw the particles do it in C : tree_walk_recursive Sat Oct 30 21:18:54 CEST 2010 -> ok, la construction du tree est ok -> need to flag the cells add : ListOfPartLists (size=maxpart + MAXTOPNODES) !!! bien trops grand... do not forget : self->ListOfPartLists[-1] = -1 pour stoper self->last = -1 (au debut) no : node sib : ? father : noeud d'en haut Sat Oct 30 22:33:06 CEST 2010 -> force_update_node_recursive, ok (pas 100% testé) +Sun Oct 31 21:28:42 CET 2010 + + +suite: modifier le calcul des potentiel, + --> somme sur la liste + + +force_treeevaluate_local_potential + +ok, semble marcher + +...> mais le temps de calcul se déteriore... + +!!! il faut du calcul directe sur les particules dans la meme cellule !!! + + + + + + + + + + + + + + + diff --git a/src/streelib/test01.py b/src/streelib/test01.py index 7d4509c..9cb6159 100755 --- a/src/streelib/test01.py +++ b/src/streelib/test01.py @@ -1,40 +1,42 @@ #!/usr/bin/env python from pNbody import * from pNbody import ic import time import treelib import Ptools as pt ErrTolTheta = 0.8 eps = 0.1 # open a model #nb = Nbody("../../examples/snap.dat",ftype='gadget') #nb = nb.selectc(nb.rxyz()<2) random.seed(1) -nb = ic.plummer(1000000,1,1,1,5,20) +nb = ic.plummer(500000,1,1,1,5,20) -MaxPartPerNode=1024 +MaxPartPerNode=4 # create a tree object print "start tree" MyTree = treelib.Tree(npart=array(nb.npart),pos=nb.pos,vel=nb.vel,mass=nb.mass,ErrTolTheta=ErrTolTheta,MaxPartPerNode=MaxPartPerNode) R = arange(-20,20,0.1) print len(R) pos = concatenate((R,zeros(len(R)),zeros(len(R)))) pos.shape = (3,len(R)) pos = transpose(pos) pos = pos.astype(float32) -pot = MyTree.Potential(pos,eps) +#pot = MyTree.Potential(pos,eps) +pot = MyTree.Potential(nb.pos,eps) +#print pot -pt.plot(R,pot) -pt.show() +#pt.plot(R,pot) +#pt.show() diff --git a/src/streelib/treelib.c b/src/streelib/treelib.c index a762a74..be8d218 100644 --- a/src/streelib/treelib.c +++ b/src/streelib/treelib.c @@ -1,1898 +1,1993 @@ #include #include #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; int MaxPartPerNode; 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; int *ListOfPartLists; /* 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,nPart; 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 i,j,p; int nextsib; int n; if(no >= self->All.MaxPart && no < self->All.MaxPart + self->MaxNodes) { /* a node */ nNodes++; /* here, we can gather info */ //printf("(%d) %g %g %g %g \n",no,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 the daughters */ for(j = 0; j < 8; j++) { if( (p = self->Nodes[no].u.suns[j]) >= 0 ) /* the node is not emtpy */ { /* useless ... ??? */ nextsib = sib; /* go to next node */ //printf(" (%d.%d)-> goto %d\n",no,j,p); tree_walk_recursive(self,p,nextsib,no); } } } else /* single particle or pseudo particle */ { //printf("found a single particle\n"); /* loop over particles in the list */ i = no; while(1) { nPart++; //printf(" particle_number %d (next=%d) (no=%d)\n",i,self->P[i].next,no); i = self->P[i].next; if (i<0) break; } } } /* 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 j,p,pp; 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"); + + + pp = no; + while(1) + { + nPart++; + pp = self->P[pp].next; + + if (pp<0) + break; + + } + 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; nPart=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",nPart); + printf("number of part : %d\n",nPart); 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; + int no, ptype, pp; 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 */ + + + + 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 */ + + pp = no; + + while(1) /* loop over particles in the list */ + { + + dx = self->P[pp].Pos[0] - pos_x; + dy = self->P[pp].Pos[1] - pos_y; + dz = self->P[pp].Pos[2] - pos_z; + mass = self->P[pp].Mass; + + r2 = dx * dx + dy * dy + dz * dz; + + + /* here we use a plummer softening */ + if (r2>0) + pot += -mass / sqrt(r2+h*h); + + + if (self->P[pp].next==-1) + break; + + pp = self->P[pp].next; + + + + } + + + no = self->Nextnode[no]; + } + 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; + 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); + /* here we use a plummer softening */ + if (r2>0) + pot += -mass / sqrt(r2+h*h); + + + } + + + + + + +// 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 */ { /* (yr) here, we need to loop over particles... */ int pp; pp = p; while(1) { nPart++; pa = &self->P[pp]; 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 (self->P[pp].next==-1) break; pp = self->P[pp].next; } } } } 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 */ { /* yr : here, we need to deal with the list */ /* I think, there is nothing to do here */ if(self->last >= 0) { if(self->last >= self->All.MaxPart) /* last is a node or pseudo */ { 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,next; 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 */ { //printf("(%d)++> a node %d\n",i,th); 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(" node not empty : subnode=%d next=%d\n",subnode,nn); parent = th; th = nn; } else { /* here we have found an empty slot where we can attach * the new particle as a leaf. */ //printf(" node empty : subnode=%d\n",subnode); //printf("add particle %d to the empty branch %d.%d\n",i,th,subnode); 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 %d (n in list=%d)\n",th,self->Nodes[parent].count[subnode]); if (self->Nodes[parent].count[subnode]>=self->All.MaxPartPerNode) { //printf(" (%d) need to distibute the particles in the list of particle %d \n",i,th); /* 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]; //printf("count=%d, for %d.%d\n",self->Nodes[parent].count[subnode],parent,subnode); /* (yr) here, we need to list over the list */ while(1) { //printf(" realocate particule %d\n",th); 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; nn = nfreep->u.suns[subnode]; if(nn >= 0) /* ok, something is in the daughter slot add to the list */ { if (nn>self->All.MaxPart) { printf("at this point, the node should contains only particles... stop !\n"); endrun(self,9878); } /* add to the list */ //printf(" add particle %d to the list of %d\n",th,nfreep->first[subnode]); /* add it to the list */ self->P[nfreep->last[subnode]].next=th; /* update parent */ nfreep->count[subnode]++; nfreep->last[subnode]=th; //printf(" count=%d, for %d.%d\n",self->Nodes[parent].count[subnode],parent,subnode); } else { /* here we have found an empty slot where we can attach * the new particle as a leaf. */ //printf(" add particle %d to the empty branch %d.%d\n",i,th,subnode); nfreep->u.suns[subnode] = th; /* (yr) set parent properties */ if (nfreep->count[subnode]!=0) { printf("we are in trouble here !\n"); endrun(self,9877); } nfreep->count[subnode] = 1; nfreep->first[subnode] = th; nfreep->last[subnode] = th; } /* next in the list */ next = self->P[th].next; self->P[th].next=-1; /* remove the link */ th = next; if (th<0) break; } 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 */ //printf("add particle %d to the list of %d\n",i,self->Nodes[parent].first[subnode]); /* 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("count=%d, for %d.%d\n",self->Nodes[parent].count[subnode],parent,subnode); break; } } } } - //exit(-1); printf("!!!particles are now inserted numnodes=%d\n\n\n",numnodes); tree_walk_and_gather(self,0); /* now compute the multipole moments recursively */ self->last = -1; nPart=0; force_update_node_recursive(self,self->All.MaxPart, -1, -1); printf("out of force_update_node_recursive : nPart = %d\n",nPart); 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; /* yr */ if(!(self->ListOfPartLists = malloc(bytes = (maxpart) * sizeof(int)))) { printf("Failed to allocate %d spaces for 'ListOfPartLists' 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); /* yr */ free(self->ListOfPartLists); 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; int MaxPartPerNode; static char *kwlist[] = {"first", "npart", "pos","vel", "mass", "ErrTolTheta", "MaxPartPerNode", NULL}; if (! PyArg_ParseTupleAndKeywords(args, kwds, "|OOOOOdi",kwlist,&first,&ntype,&pos,&vel,&mass,&ErrTolTheta,&MaxPartPerNode)) 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.MaxPartPerNode = MaxPartPerNode; 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); }