GADGET-2: ngb.c Source File

00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <string.h>
00004 #include <math.h>
00005 #include <time.h>
00006 #include <mpi.h>
00007 
00008 #include "allvars.h"
00009 #include "proto.h"
00010 
00011 
00019 #ifdef PERIODIC
00020 static double boxSize, boxHalf;
00021 
00022 #ifdef LONG_X
00023 static double boxSize_X, boxHalf_X;
00024 #else
00025 #define boxSize_X boxSize
00026 #define boxHalf_X boxHalf
00027 #endif
00028 #ifdef LONG_Y
00029 static double boxSize_Y, boxHalf_Y;
00030 #else
00031 #define boxSize_Y boxSize
00032 #define boxHalf_Y boxHalf
00033 #endif
00034 #ifdef LONG_Z
00035 static double boxSize_Z, boxHalf_Z;
00036 #else
00037 #define boxSize_Z boxSize
00038 #define boxHalf_Z boxHalf
00039 #endif
00040 #endif
00041 
00042 
00047 #define NGB_PERIODIC_X(x) (xtmp=(x),(xtmp>boxHalf_X)?(xtmp-boxSize_X):((xtmp<-boxHalf_X)?(xtmp+boxSize_X):xtmp))
00048 #define NGB_PERIODIC_Y(x) (xtmp=(x),(xtmp>boxHalf_Y)?(xtmp-boxSize_Y):((xtmp<-boxHalf_Y)?(xtmp+boxSize_Y):xtmp))
00049 #define NGB_PERIODIC_Z(x) (xtmp=(x),(xtmp>boxHalf_Z)?(xtmp-boxSize_Z):((xtmp<-boxHalf_Z)?(xtmp+boxSize_Z):xtmp))
00050 
00051 
00052 
00064 int ngb_treefind_pairs(FLOAT searchcenter[3], FLOAT hsml, int *startnode)
00065 {
00066   int k, no, p, numngb;
00067   FLOAT hdiff;
00068   FLOAT searchmin[3], searchmax[3];
00069   struct NODE *this;
00070 
00071 #ifdef PERIODIC
00072   double xtmp;
00073 #endif
00074 
00075   for(k = 0; k < 3; k++)        /* cube-box window */
00076     {
00077       searchmin[k] = searchcenter[k] - hsml;
00078       searchmax[k] = searchcenter[k] + hsml;
00079     }
00080 
00081   numngb = 0;
00082   no = *startnode;
00083 
00084   while(no >= 0)
00085     {
00086       if(no < All.MaxPart)      /* single particle */
00087         {
00088           p = no;
00089           no = Nextnode[no];
00090 
00091           if(P[p].Type > 0)
00092             continue;
00093 
00094           hdiff = SphP[p].Hsml - hsml;
00095           if(hdiff < 0)
00096             hdiff = 0;
00097 
00098 #ifdef PERIODIC
00099           if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) < (-hsml - hdiff))
00100             continue;
00101           if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) > (hsml + hdiff))
00102             continue;
00103           if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) < (-hsml - hdiff))
00104             continue;
00105           if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) > (hsml + hdiff))
00106             continue;
00107           if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) < (-hsml - hdiff))
00108             continue;
00109           if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) > (hsml + hdiff))
00110             continue;
00111 #else
00112           if(P[p].Pos[0] < (searchmin[0] - hdiff))
00113             continue;
00114           if(P[p].Pos[0] > (searchmax[0] + hdiff))
00115             continue;
00116           if(P[p].Pos[1] < (searchmin[1] - hdiff))
00117             continue;
00118           if(P[p].Pos[1] > (searchmax[1] + hdiff))
00119             continue;
00120           if(P[p].Pos[2] < (searchmin[2] - hdiff))
00121             continue;
00122           if(P[p].Pos[2] > (searchmax[2] + hdiff))
00123             continue;
00124 #endif
00125           Ngblist[numngb++] = p;
00126 
00127           if(numngb == MAX_NGB)
00128             {
00129               printf
00130                 ("ThisTask=%d: Need to do a second neighbour loop in hydro-force for (%g|%g|%g) hsml=%g no=%d\n",
00131                  ThisTask, searchcenter[0], searchcenter[1], searchcenter[2], hsml, no);
00132               *startnode = no;
00133               return numngb;
00134             }
00135         }
00136       else
00137         {
00138           if(no >= All.MaxPart + MaxNodes)      /* pseudo particle */
00139             {
00140               Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
00141               no = Nextnode[no - MaxNodes];
00142               continue;
00143             }
00144 
00145           this = &Nodes[no];
00146           hdiff = Extnodes[no].hmax - hsml;
00147           if(hdiff < 0)
00148             hdiff = 0;
00149 
00150           no = this->u.d.sibling;       /* in case the node can be discarded */
00151 
00152 #ifdef PERIODIC
00153           if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) + 0.5 * this->len) < (-hsml - hdiff))
00154             continue;
00155           if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) - 0.5 * this->len) > (hsml + hdiff))
00156             continue;
00157           if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) + 0.5 * this->len) < (-hsml - hdiff))
00158             continue;
00159           if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) - 0.5 * this->len) > (hsml + hdiff))
00160             continue;
00161           if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) + 0.5 * this->len) < (-hsml - hdiff))
00162             continue;
00163           if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) - 0.5 * this->len) > (hsml + hdiff))
00164             continue;
00165 #else
00166           if((this->center[0] + 0.5 * this->len) < (searchmin[0] - hdiff))
00167             continue;
00168           if((this->center[0] - 0.5 * this->len) > (searchmax[0] + hdiff))
00169             continue;
00170           if((this->center[1] + 0.5 * this->len) < (searchmin[1] - hdiff))
00171             continue;
00172           if((this->center[1] - 0.5 * this->len) > (searchmax[1] + hdiff))
00173             continue;
00174           if((this->center[2] + 0.5 * this->len) < (searchmin[2] - hdiff))
00175             continue;
00176           if((this->center[2] - 0.5 * this->len) > (searchmax[2] + hdiff))
00177             continue;
00178 #endif
00179           no = this->u.d.nextnode;      /* ok, we need to open the node */
00180         }
00181     }
00182 
00183   *startnode = -1;
00184   return numngb;
00185 }
00186 
00187 
00188 
00194 int ngb_treefind_variable(FLOAT searchcenter[3], FLOAT hsml, int *startnode)
00195 {
00196   int k, numngb;
00197   int no, p;
00198   struct NODE *this;
00199   FLOAT searchmin[3], searchmax[3];
00200 
00201 #ifdef PERIODIC
00202   double xtmp;
00203 #endif
00204 
00205   for(k = 0; k < 3; k++)        /* cube-box window */
00206     {
00207       searchmin[k] = searchcenter[k] - hsml;
00208       searchmax[k] = searchcenter[k] + hsml;
00209     }
00210 
00211   numngb = 0;
00212   no = *startnode;
00213 
00214   while(no >= 0)
00215     {
00216       if(no < All.MaxPart)      /* single particle */
00217         {
00218           p = no;
00219           no = Nextnode[no];
00220 
00221           if(P[p].Type > 0)
00222             continue;
00223 
00224 #ifdef PERIODIC
00225           if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) < -hsml)
00226             continue;
00227           if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) > hsml)
00228             continue;
00229           if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) < -hsml)
00230             continue;
00231           if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) > hsml)
00232             continue;
00233           if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) < -hsml)
00234             continue;
00235           if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) > hsml)
00236             continue;
00237 #else
00238           if(P[p].Pos[0] < searchmin[0])
00239             continue;
00240           if(P[p].Pos[0] > searchmax[0])
00241             continue;
00242           if(P[p].Pos[1] < searchmin[1])
00243             continue;
00244           if(P[p].Pos[1] > searchmax[1])
00245             continue;
00246           if(P[p].Pos[2] < searchmin[2])
00247             continue;
00248           if(P[p].Pos[2] > searchmax[2])
00249             continue;
00250 #endif
00251           Ngblist[numngb++] = p;
00252 
00253           if(numngb == MAX_NGB)
00254             {
00255               numngb = ngb_clear_buf(searchcenter, hsml, numngb);
00256               if(numngb == MAX_NGB)
00257                 {
00258                   printf("ThisTask=%d: Need to do a second neighbour loop for (%g|%g|%g) hsml=%g no=%d\n",
00259                          ThisTask, searchcenter[0], searchcenter[1], searchcenter[2], hsml, no);
00260                   *startnode = no;
00261                   return numngb;
00262                 }
00263             }
00264         }
00265       else
00266         {
00267           if(no >= All.MaxPart + MaxNodes)      /* pseudo particle */
00268             {
00269               Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
00270               no = Nextnode[no - MaxNodes];
00271               continue;
00272             }
00273 
00274           this = &Nodes[no];
00275 
00276           no = this->u.d.sibling;       /* in case the node can be discarded */
00277 #ifdef PERIODIC
00278           if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) + 0.5 * this->len) < -hsml)
00279             continue;
00280           if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) - 0.5 * this->len) > hsml)
00281             continue;
00282           if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) + 0.5 * this->len) < -hsml)
00283             continue;
00284           if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) - 0.5 * this->len) > hsml)
00285             continue;
00286           if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) + 0.5 * this->len) < -hsml)
00287             continue;
00288           if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) - 0.5 * this->len) > hsml)
00289             continue;
00290 #else
00291           if((this->center[0] + 0.5 * this->len) < (searchmin[0]))
00292             continue;
00293           if((this->center[0] - 0.5 * this->len) > (searchmax[0]))
00294             continue;
00295           if((this->center[1] + 0.5 * this->len) < (searchmin[1]))
00296             continue;
00297           if((this->center[1] - 0.5 * this->len) > (searchmax[1]))
00298             continue;
00299           if((this->center[2] + 0.5 * this->len) < (searchmin[2]))
00300             continue;
00301           if((this->center[2] - 0.5 * this->len) > (searchmax[2]))
00302             continue;
00303 #endif
00304           no = this->u.d.nextnode;      /* ok, we need to open the node */
00305         }
00306     }
00307 
00308   *startnode = -1;
00309   return numngb;
00310 }
00311 
00312 
00313 
00314 
00320 int ngb_clear_buf(FLOAT searchcenter[3], FLOAT hsml, int numngb)
00321 {
00322   int i, p;
00323   FLOAT dx, dy, dz, r2;
00324 
00325 #ifdef PERIODIC
00326   double xtmp;
00327 #endif
00328 
00329   for(i = 0; i < numngb; i++)
00330     {
00331       p = Ngblist[i];
00332 #ifdef PERIODIC
00333       dx = NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]);
00334       dy = NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]);
00335       dz = NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]);
00336 #else
00337       dx = P[p].Pos[0] - searchcenter[0];
00338       dy = P[p].Pos[1] - searchcenter[1];
00339       dz = P[p].Pos[2] - searchcenter[2];
00340 #endif
00341       r2 = dx * dx + dy * dy + dz * dz;
00342 
00343       if(r2 > hsml * hsml)
00344         {
00345           Ngblist[i] = Ngblist[numngb - 1];
00346           i--;
00347           numngb--;
00348         }
00349     }
00350 
00351   return numngb;
00352 }
00353 
00354 
00355 
00358 void ngb_treeallocate(int npart)
00359 {
00360   double totbytes = 0;
00361   size_t bytes;
00362 
00363 #ifdef PERIODIC
00364   boxSize = All.BoxSize;
00365   boxHalf = 0.5 * All.BoxSize;
00366 #ifdef LONG_X
00367   boxHalf_X = boxHalf * LONG_X;
00368   boxSize_X = boxSize * LONG_X;
00369 #endif
00370 #ifdef LONG_Y
00371   boxHalf_Y = boxHalf * LONG_Y;
00372   boxSize_Y = boxSize * LONG_Y;
00373 #endif
00374 #ifdef LONG_Z
00375   boxHalf_Z = boxHalf * LONG_Z;
00376   boxSize_Z = boxSize * LONG_Z;
00377 #endif
00378 #endif
00379 
00380   if(!(Ngblist = malloc(bytes = npart * (long) sizeof(int))))
00381     {
00382       printf("Failed to allocate %g MB for ngblist array\n", bytes / (1024.0 * 1024.0));
00383       endrun(78);
00384     }
00385   totbytes += bytes;
00386 
00387   if(ThisTask == 0)
00388     printf("allocated %g Mbyte for ngb search.\n", totbytes / (1024.0 * 1024.0));
00389 }
00390 
00391 
00394 void ngb_treefree(void)
00395 {
00396   free(Ngblist);
00397 }
00398 
00403 void ngb_treebuild(void)
00404 {
00405   if(ThisTask == 0)
00406     printf("Begin Ngb-tree construction.\n");
00407 
00408   force_treebuild(N_gas);
00409 
00410   if(ThisTask == 0)
00411     printf("Ngb-Tree contruction finished \n");
00412 }
00413