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++)
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)
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)
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;
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;
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++)
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)
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)
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;
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;
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