Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65689760
ngb.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Jun 5, 13:41
Size
32 KB
Mime Type
text/x-c
Expires
Fri, Jun 7, 13:41 (2 d)
Engine
blob
Format
Raw Data
Handle
18112808
Attached To
rGEAR Gear
ngb.c
View Options
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
/*! \file ngb.c
* \brief neighbour search by means of the tree
*
* This file contains routines for neighbour finding. We use the
* gravity-tree and a range-searching technique to find neighbours.
*/
#ifdef PERIODIC
static double boxSize, boxHalf;
#ifdef LONG_X
static double boxSize_X, boxHalf_X;
#else
#define boxSize_X boxSize
#define boxHalf_X boxHalf
#endif
#ifdef LONG_Y
static double boxSize_Y, boxHalf_Y;
#else
#define boxSize_Y boxSize
#define boxHalf_Y boxHalf
#endif
#ifdef LONG_Z
static double boxSize_Z, boxHalf_Z;
#else
#define boxSize_Z boxSize
#define boxHalf_Z boxHalf
#endif
#endif
/*! these macros maps a coordinate difference to the nearest periodic
* image
*/
#define NGB_PERIODIC_X(x) (xtmp=(x),(xtmp>boxHalf_X)?(xtmp-boxSize_X):((xtmp<-boxHalf_X)?(xtmp+boxSize_X):xtmp))
#define NGB_PERIODIC_Y(x) (xtmp=(x),(xtmp>boxHalf_Y)?(xtmp-boxSize_Y):((xtmp<-boxHalf_Y)?(xtmp+boxSize_Y):xtmp))
#define NGB_PERIODIC_Z(x) (xtmp=(x),(xtmp>boxHalf_Z)?(xtmp-boxSize_Z):((xtmp<-boxHalf_Z)?(xtmp+boxSize_Z):xtmp))
/*! This routine finds all neighbours `j' that can interact with the
* particle `i' in the communication buffer.
*
* Note that an interaction can take place if
* \f$ r_{ij} < h_i \f$ OR if \f$ r_{ij} < h_j \f$.
*
* In the range-search this is taken into account, i.e. it is guaranteed that
* all particles are found that fulfil this condition, including the (more
* difficult) second part of it. For this purpose, each node knows the
* maximum h occuring among the particles it represents.
*/
int ngb_treefind_pairs(FLOAT searchcenter[3], FLOAT hsml, int phase, int *startnode)
{
int k, no, p, numngb;
FLOAT hdiff;
FLOAT searchmin[3], searchmax[3];
struct NODE *this;
#ifdef PERIODIC
double xtmp;
#endif
for(k = 0; k < 3; k++) /* cube-box window */
{
searchmin[k] = searchcenter[k] - hsml;
searchmax[k] = searchcenter[k] + hsml;
}
numngb = 0;
no = *startnode;
while(no >= 0)
{
if(no < All.MaxPart) /* single particle */
{
p = no;
no = Nextnode[no];
if(P[p].Type > 0)
continue;
#ifdef MULTIPHASE
/* WARNING - WARNING - WARNING - WARNING */
/* keep only particles of same phase */
//if(SphP[p].Phase != phase)
// continue;
/* WARNING - WARNING - WARNING - WARNING */
/* this may be problematic for the feedback */
/* sticky sees only sticky */
//if ((phase==GAS_STICKY) && SphP[p].Phase!=GAS_STICKY)
// continue;
#endif
hdiff = SphP[p].Hsml - hsml;
if(hdiff < 0)
hdiff = 0;
#ifdef PERIODIC
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) < (-hsml - hdiff))
continue;
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) > (hsml + hdiff))
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) < (-hsml - hdiff))
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) > (hsml + hdiff))
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) < (-hsml - hdiff))
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) > (hsml + hdiff))
continue;
#else
if(P[p].Pos[0] < (searchmin[0] - hdiff))
continue;
if(P[p].Pos[0] > (searchmax[0] + hdiff))
continue;
if(P[p].Pos[1] < (searchmin[1] - hdiff))
continue;
if(P[p].Pos[1] > (searchmax[1] + hdiff))
continue;
if(P[p].Pos[2] < (searchmin[2] - hdiff))
continue;
if(P[p].Pos[2] > (searchmax[2] + hdiff))
continue;
#endif
Ngblist[numngb++] = p;
if(numngb == MAX_NGB)
{
printf
("ThisTask=%d: Need to do a second neighbour loop in hydro-force for (%g|%g|%g) hsml=%g no=%d\n",
ThisTask, searchcenter[0], searchcenter[1], searchcenter[2], hsml, no);
*startnode = no;
return numngb;
}
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
no = Nextnode[no - MaxNodes];
continue;
}
this = &Nodes[no];
hdiff = Extnodes[no].hmax - hsml;
if(hdiff < 0)
hdiff = 0;
no = this->u.d.sibling; /* in case the node can be discarded */
#ifdef PERIODIC
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) + 0.5 * this->len) < (-hsml - hdiff))
continue;
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) - 0.5 * this->len) > (hsml + hdiff))
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) + 0.5 * this->len) < (-hsml - hdiff))
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) - 0.5 * this->len) > (hsml + hdiff))
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) + 0.5 * this->len) < (-hsml - hdiff))
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) - 0.5 * this->len) > (hsml + hdiff))
continue;
#else
if((this->center[0] + 0.5 * this->len) < (searchmin[0] - hdiff))
continue;
if((this->center[0] - 0.5 * this->len) > (searchmax[0] + hdiff))
continue;
if((this->center[1] + 0.5 * this->len) < (searchmin[1] - hdiff))
continue;
if((this->center[1] - 0.5 * this->len) > (searchmax[1] + hdiff))
continue;
if((this->center[2] + 0.5 * this->len) < (searchmin[2] - hdiff))
continue;
if((this->center[2] - 0.5 * this->len) > (searchmax[2] + hdiff))
continue;
#endif
no = this->u.d.nextnode; /* ok, we need to open the node */
}
}
*startnode = -1;
return numngb;
}
/*! This routine finds all neighbours `j' that can interact with the
* particle `i' in the communication buffer.
*
* Note that an interaction can take place if
* \f$ r_{ij} < h_i \f$ OR if \f$ r_{ij} < h_j \f$.
*
* In the range-search this is taken into account, i.e. it is guaranteed that
* all particles are found that fulfil this condition, including the (more
* difficult) second part of it. For this purpose, each node knows the
* maximum h occuring among the particles it represents.
*
* This function is similar to ngb_treefind_pairs but returns only
* pairs of same phase
*/
#ifdef MULTIPHASE
int ngb_treefind_phase_pairs(FLOAT searchcenter[3], FLOAT hsml, int phase, int *startnode)
{
int k, no, p, numngb;
FLOAT hdiff;
FLOAT searchmin[3], searchmax[3];
struct NODE *this;
#ifdef PERIODIC
double xtmp;
#endif
for(k = 0; k < 3; k++) /* cube-box window */
{
searchmin[k] = searchcenter[k] - hsml;
searchmax[k] = searchcenter[k] + hsml;
}
numngb = 0;
no = *startnode;
while(no >= 0)
{
if(no < All.MaxPart) /* single particle */
{
p = no;
no = Nextnode[no];
if(P[p].Type > 0)
continue;
/* here, we select all gas particles */
#ifdef MULTIPHASE
/* keep only particles of same phase */
if(SphP[p].Phase != phase)
continue;
#endif
hdiff = SphP[p].Hsml - hsml;
if(hdiff < 0)
hdiff = 0;
#ifdef PERIODIC
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) < (-hsml - hdiff))
continue;
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) > (hsml + hdiff))
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) < (-hsml - hdiff))
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) > (hsml + hdiff))
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) < (-hsml - hdiff))
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) > (hsml + hdiff))
continue;
#else
if(P[p].Pos[0] < (searchmin[0] - hdiff))
continue;
if(P[p].Pos[0] > (searchmax[0] + hdiff))
continue;
if(P[p].Pos[1] < (searchmin[1] - hdiff))
continue;
if(P[p].Pos[1] > (searchmax[1] + hdiff))
continue;
if(P[p].Pos[2] < (searchmin[2] - hdiff))
continue;
if(P[p].Pos[2] > (searchmax[2] + hdiff))
continue;
#endif
Ngblist[numngb++] = p;
if(numngb == MAX_NGB)
{
printf
("ThisTask=%d: Need to do a second neighbour loop in hydro-force for (%g|%g|%g) hsml=%g no=%d\n",
ThisTask, searchcenter[0], searchcenter[1], searchcenter[2], hsml, no);
*startnode = no;
return numngb;
}
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
no = Nextnode[no - MaxNodes];
continue;
}
this = &Nodes[no];
hdiff = Extnodes[no].hmax - hsml;
if(hdiff < 0)
hdiff = 0;
no = this->u.d.sibling; /* in case the node can be discarded */
#ifdef PERIODIC
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) + 0.5 * this->len) < (-hsml - hdiff))
continue;
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) - 0.5 * this->len) > (hsml + hdiff))
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) + 0.5 * this->len) < (-hsml - hdiff))
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) - 0.5 * this->len) > (hsml + hdiff))
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) + 0.5 * this->len) < (-hsml - hdiff))
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) - 0.5 * this->len) > (hsml + hdiff))
continue;
#else
if((this->center[0] + 0.5 * this->len) < (searchmin[0] - hdiff))
continue;
if((this->center[0] - 0.5 * this->len) > (searchmax[0] + hdiff))
continue;
if((this->center[1] + 0.5 * this->len) < (searchmin[1] - hdiff))
continue;
if((this->center[1] - 0.5 * this->len) > (searchmax[1] + hdiff))
continue;
if((this->center[2] + 0.5 * this->len) < (searchmin[2] - hdiff))
continue;
if((this->center[2] - 0.5 * this->len) > (searchmax[2] + hdiff))
continue;
#endif
no = this->u.d.nextnode; /* ok, we need to open the node */
}
}
*startnode = -1;
return numngb;
}
#endif
/*! This function returns neighbours with distance <= hsml and returns them in
* Ngblist. Actually, particles in a box of half side length hsml are
* returned, i.e. the reduction to a sphere still needs to be done in the
* calling routine.
*/
int ngb_treefind_variable(FLOAT searchcenter[3], FLOAT hsml, int phase, int *startnode)
{
int k, numngb;
int no, p;
struct NODE *this;
FLOAT searchmin[3], searchmax[3];
#ifdef PERIODIC
double xtmp;
#endif
for(k = 0; k < 3; k++) /* cube-box window */
{
searchmin[k] = searchcenter[k] - hsml;
searchmax[k] = searchcenter[k] + hsml;
}
numngb = 0;
no = *startnode;
while(no >= 0)
{
if(no < All.MaxPart) /* single particle */
{
p = no;
no = Nextnode[no];
if(P[p].Type > 0)
continue;
#ifdef MULTIPHASE
/* GAS_SPH sees only GAS_SPH */
if((phase == GAS_SPH) && (SphP[p].Phase != GAS_SPH))
continue;
/* GAS_SPH sees only GAS_STICKY (better for sticky interaction) */
if((phase == GAS_STICKY) && (SphP[p].Phase != GAS_STICKY))
continue;
/* GAS_DARK sees everybody */
#endif
#ifdef PERIODIC
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) < -hsml)
continue;
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) > hsml)
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) < -hsml)
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) > hsml)
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) < -hsml)
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) > hsml)
continue;
#else
if(P[p].Pos[0] < searchmin[0])
continue;
if(P[p].Pos[0] > searchmax[0])
continue;
if(P[p].Pos[1] < searchmin[1])
continue;
if(P[p].Pos[1] > searchmax[1])
continue;
if(P[p].Pos[2] < searchmin[2])
continue;
if(P[p].Pos[2] > searchmax[2])
continue;
#endif
Ngblist[numngb++] = p;
if(numngb == MAX_NGB)
{
numngb = ngb_clear_buf(searchcenter, hsml, numngb);
if(numngb == MAX_NGB)
{
printf("ThisTask=%d: Need to do a second neighbour loop for (%g|%g|%g) hsml=%g no=%d\n",
ThisTask, searchcenter[0], searchcenter[1], searchcenter[2], hsml, no);
*startnode = no;
return numngb;
}
}
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
no = Nextnode[no - MaxNodes];
continue;
}
this = &Nodes[no];
no = this->u.d.sibling; /* in case the node can be discarded */
#ifdef PERIODIC
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) - 0.5 * this->len) > hsml)
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) - 0.5 * this->len) > hsml)
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) - 0.5 * this->len) > hsml)
continue;
#else
if((this->center[0] + 0.5 * this->len) < (searchmin[0]))
continue;
if((this->center[0] - 0.5 * this->len) > (searchmax[0]))
continue;
if((this->center[1] + 0.5 * this->len) < (searchmin[1]))
continue;
if((this->center[1] - 0.5 * this->len) > (searchmax[1]))
continue;
if((this->center[2] + 0.5 * this->len) < (searchmin[2]))
continue;
if((this->center[2] - 0.5 * this->len) > (searchmax[2]))
continue;
#endif
no = this->u.d.nextnode; /* ok, we need to open the node */
}
}
*startnode = -1;
return numngb;
}
#ifdef CHIMIE
/*! This function returns neighbours with distance <= hsml and returns them in
* Ngblist. Actually, particles in a box of half side length hsml are
* returned, i.e. the reduction to a sphere still needs to be done in the
* calling routine.
*/
int ngb_treefind_variable_for_chimie(FLOAT searchcenter[3], FLOAT hsml, int *startnode)
{
int k, numngb;
int no, p;
struct NODE *this;
FLOAT searchmin[3], searchmax[3];
#ifdef PERIODIC
double xtmp;
#endif
for(k = 0; k < 3; k++) /* cube-box window */
{
searchmin[k] = searchcenter[k] - hsml;
searchmax[k] = searchcenter[k] + hsml;
}
numngb = 0;
no = *startnode;
while(no >= 0)
{
if(no < All.MaxPart) /* single particle */
{
p = no;
no = Nextnode[no];
if(P[p].Type > 0)
continue;
#ifdef CHIMIE_KINETIC_FEEDBACK /* no, we need to keep all particles, especially when the wind is low */
// if (SphP[p].WindTime >= (All.Time-All.ChimieWindTime))
// continue;
#endif
#ifdef PERIODIC
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) < -hsml)
continue;
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) > hsml)
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) < -hsml)
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) > hsml)
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) < -hsml)
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) > hsml)
continue;
#else
if(P[p].Pos[0] < searchmin[0])
continue;
if(P[p].Pos[0] > searchmax[0])
continue;
if(P[p].Pos[1] < searchmin[1])
continue;
if(P[p].Pos[1] > searchmax[1])
continue;
if(P[p].Pos[2] < searchmin[2])
continue;
if(P[p].Pos[2] > searchmax[2])
continue;
#endif
Ngblist[numngb++] = p;
if(numngb == MAX_NGB)
{
numngb = ngb_clear_buf(searchcenter, hsml, numngb);
if(numngb == MAX_NGB)
{
printf("ThisTask=%d: Need to do a second neighbour loop for (%g|%g|%g) hsml=%g no=%d\n",
ThisTask, searchcenter[0], searchcenter[1], searchcenter[2], hsml, no);
*startnode = no;
return numngb;
}
}
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
no = Nextnode[no - MaxNodes];
continue;
}
this = &Nodes[no];
no = this->u.d.sibling; /* in case the node can be discarded */
#ifdef PERIODIC
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) - 0.5 * this->len) > hsml)
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) - 0.5 * this->len) > hsml)
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) - 0.5 * this->len) > hsml)
continue;
#else
if((this->center[0] + 0.5 * this->len) < (searchmin[0]))
continue;
if((this->center[0] - 0.5 * this->len) > (searchmax[0]))
continue;
if((this->center[1] + 0.5 * this->len) < (searchmin[1]))
continue;
if((this->center[1] - 0.5 * this->len) > (searchmax[1]))
continue;
if((this->center[2] + 0.5 * this->len) < (searchmin[2]))
continue;
if((this->center[2] - 0.5 * this->len) > (searchmax[2]))
continue;
#endif
no = this->u.d.nextnode; /* ok, we need to open the node */
}
}
*startnode = -1;
return numngb;
}
#endif
/*! This function returns the first particle that may sticky-interact with
* the given particle, in the sph kernel.
*/
#ifdef MULTIPHASE
int ngb_treefind_sticky_collisions(FLOAT searchcenter[3], FLOAT hsml, int phase, int *startnode)
{
int k, numngb;
int no, p;
struct NODE *this;
FLOAT searchmin[3], searchmax[3];
double dx,dy,dz,dr2;
#ifdef PERIODIC
double xtmp;
#endif
for(k = 0; k < 3; k++) /* cube-box window */
{
searchmin[k] = searchcenter[k] - hsml;
searchmax[k] = searchcenter[k] + hsml;
}
numngb = 0;
no = *startnode;
while(no >= 0)
{
if(no < All.MaxPart) /* single particle */
{
p = no;
no = Nextnode[no];
if((P[p].Type > 0) || (SphP[p].Phase!=GAS_STICKY)) /* keep only GAS_STICKY */
continue;
if(SphP[p].StickyFlag==0) /* the particle may not collide */
continue;
#ifdef PERIODIC
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) < -hsml)
continue;
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) > hsml)
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) < -hsml)
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) > hsml)
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) < -hsml)
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) > hsml)
continue;
#else
if(P[p].Pos[0] < searchmin[0])
continue;
if(P[p].Pos[0] > searchmax[0])
continue;
if(P[p].Pos[1] < searchmin[1])
continue;
if(P[p].Pos[1] > searchmax[1])
continue;
if(P[p].Pos[2] < searchmin[2])
continue;
if(P[p].Pos[2] > searchmax[2])
continue;
#endif
dx = searchcenter[0] - P[p].Pos[0];
dy = searchcenter[1] - P[p].Pos[1];
dz = searchcenter[2] - P[p].Pos[2];
dr2= dx * dx + dy * dy + dz * dz;
if (dr2>hsml*hsml)
continue;
if (dr2==0)
continue;
return p;
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
/* here, we do nothing */
no = Nextnode[no - MaxNodes];
continue;
}
this = &Nodes[no];
no = this->u.d.sibling; /* in case the node can be discarded */
#ifdef PERIODIC
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) - 0.5 * this->len) > hsml)
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) - 0.5 * this->len) > hsml)
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) - 0.5 * this->len) > hsml)
continue;
#else
if((this->center[0] + 0.5 * this->len) < (searchmin[0]))
continue;
if((this->center[0] - 0.5 * this->len) > (searchmax[0]))
continue;
if((this->center[1] + 0.5 * this->len) < (searchmin[1]))
continue;
if((this->center[1] - 0.5 * this->len) > (searchmax[1]))
continue;
if((this->center[2] + 0.5 * this->len) < (searchmin[2]))
continue;
if((this->center[2] - 0.5 * this->len) > (searchmax[2]))
continue;
#endif
no = this->u.d.nextnode; /* ok, we need to open the node */
}
}
return -1;
}
#endif
#ifdef TESSEL
#define NGB_PERIODIC_FLAG_X(x) (xtmp=(x),(xtmp>boxHalf_X)?(1):((xtmp<-boxHalf_X)?(1):0))
#define NGB_PERIODIC_FLAG_Y(x) (xtmp=(x),(xtmp>boxHalf_Y)?(1):((xtmp<-boxHalf_Y)?(1):0))
#define NGB_PERIODIC_FLAG_Z(x) (xtmp=(x),(xtmp>boxHalf_Z)?(1):((xtmp<-boxHalf_Z)?(1):0))
/*! This function returns neighbours with distance <= hsml and returns them in
* Ngblist. Actually, particles in a box of half side length hsml are
* returned, i.e. the reduction to a sphere still needs to be done in the
* calling routine.
*
*
*/
int ngb_treefind_variable_for_tessel(FLOAT searchcenter[3], FLOAT hsml, int phase, int *startnode)
{
int k, numngb;
int no, p;
struct NODE *this;
FLOAT searchmin[3], searchmax[3];
int flag;
#ifdef PERIODIC
double xtmp;
#endif
for(k = 0; k < 3; k++) /* cube-box window */
{
searchmin[k] = searchcenter[k] - hsml;
searchmax[k] = searchcenter[k] + hsml;
}
numngb = 0;
no = *startnode;
while(no >= 0)
{
if(no < All.MaxPart) /* single particle */
{
p = no;
no = Nextnode[no];
if(P[p].Type > 0)
continue;
//if(P[p].IsAdded) /* the point has already be added by another point, skip it */
// continue;
#ifdef MULTIPHASE
/* GAS_SPH sees only GAS_SPH */
if((phase == GAS_SPH) && (SphP[p].Phase != GAS_SPH))
continue;
/* GAS_SPH sees only GAS_STICKY (better for sticky interaction) */
if((phase == GAS_STICKY) && (SphP[p].Phase != GAS_STICKY))
continue;
/* GAS_DARK sees everybody */
#endif
#ifdef PERIODIC
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) < -hsml)
continue;
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) > hsml)
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) < -hsml)
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) > hsml)
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) < -hsml)
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) > hsml)
continue;
#else
if(P[p].Pos[0] < searchmin[0])
continue;
if(P[p].Pos[0] > searchmax[0])
continue;
if(P[p].Pos[1] < searchmin[1])
continue;
if(P[p].Pos[1] > searchmax[1])
continue;
if(P[p].Pos[2] < searchmin[2])
continue;
if(P[p].Pos[2] > searchmax[2])
continue;
#endif
/* the particle is inside the search radius */
#ifdef PERIODIC
/* keep it only if the particle is in the radius due to periodicity */
flag = 0; /* !!! here, we should stop once flag=1 */
flag += (NGB_PERIODIC_FLAG_X(P[p].Pos[0] - searchcenter[0]));
flag += (NGB_PERIODIC_FLAG_X(P[p].Pos[0] - searchcenter[0]));
flag += (NGB_PERIODIC_FLAG_Y(P[p].Pos[1] - searchcenter[1]));
flag += (NGB_PERIODIC_FLAG_Y(P[p].Pos[1] - searchcenter[1]));
flag += (NGB_PERIODIC_FLAG_Z(P[p].Pos[2] - searchcenter[2]));
flag += (NGB_PERIODIC_FLAG_Z(P[p].Pos[2] - searchcenter[2]));
if (flag==0)
continue;
#endif
Ngblist[numngb++] = p;
if(numngb == MAX_NGB)
{
numngb = ngb_clear_buf(searchcenter, hsml, numngb);
if(numngb == MAX_NGB)
{
printf("ThisTask=%d: Need to do a second neighbour loop for (%g|%g|%g) hsml=%g no=%d\n",
ThisTask, searchcenter[0], searchcenter[1], searchcenter[2], hsml, no);
*startnode = no;
return numngb;
}
}
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
no = Nextnode[no - MaxNodes];
continue;
}
this = &Nodes[no];
no = this->u.d.sibling; /* in case the node can be discarded */
#ifdef PERIODIC
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) - 0.5 * this->len) > hsml)
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) - 0.5 * this->len) > hsml)
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) - 0.5 * this->len) > hsml)
continue;
#else
if((this->center[0] + 0.5 * this->len) < (searchmin[0]))
continue;
if((this->center[0] - 0.5 * this->len) > (searchmax[0]))
continue;
if((this->center[1] + 0.5 * this->len) < (searchmin[1]))
continue;
if((this->center[1] - 0.5 * this->len) > (searchmax[1]))
continue;
if((this->center[2] + 0.5 * this->len) < (searchmin[2]))
continue;
if((this->center[2] - 0.5 * this->len) > (searchmax[2]))
continue;
#endif
no = this->u.d.nextnode; /* ok, we need to open the node */
}
}
*startnode = -1;
return numngb;
}
#endif /* TESSEL */
/*! The buffer for the neighbour list has a finite length MAX_NGB. For a large
* search region, this buffer can get full, in which case this routine can be
* called to eliminate some of the superfluous particles in the "corners" of
* the search box - only the ones in the inscribed sphere need to be kept.
*/
int ngb_clear_buf(FLOAT searchcenter[3], FLOAT hsml, int numngb)
{
int i, p;
FLOAT dx, dy, dz, r2;
#ifdef PERIODIC
double xtmp;
#endif
for(i = 0; i < numngb; i++)
{
p = Ngblist[i];
#ifdef PERIODIC
dx = NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]);
dy = NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]);
dz = NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]);
#else
dx = P[p].Pos[0] - searchcenter[0];
dy = P[p].Pos[1] - searchcenter[1];
dz = P[p].Pos[2] - searchcenter[2];
#endif
r2 = dx * dx + dy * dy + dz * dz;
if(r2 > hsml * hsml)
{
Ngblist[i] = Ngblist[numngb - 1];
i--;
numngb--;
}
}
return numngb;
}
/*! Allocates memory for the neighbour list buffer.
*/
void ngb_treeallocate(int npart)
{
double totbytes = 0;
size_t bytes;
#ifdef PERIODIC
boxSize = All.BoxSize;
boxHalf = 0.5 * All.BoxSize;
#ifdef LONG_X
boxHalf_X = boxHalf * LONG_X;
boxSize_X = boxSize * LONG_X;
#endif
#ifdef LONG_Y
boxHalf_Y = boxHalf * LONG_Y;
boxSize_Y = boxSize * LONG_Y;
#endif
#ifdef LONG_Z
boxHalf_Z = boxHalf * LONG_Z;
boxSize_Z = boxSize * LONG_Z;
#endif
#endif
if(!(Ngblist = malloc(bytes = npart * (long) sizeof(int))))
{
printf("Failed to allocate %g MB for ngblist array\n", bytes / (1024.0 * 1024.0));
endrun(78);
}
totbytes += bytes;
if(ThisTask == 0)
printf("allocated %g Mbyte for ngb search.\n", totbytes / (1024.0 * 1024.0));
}
/*! free memory allocated for neighbour list buffer.
*/
void ngb_treefree(void)
{
free(Ngblist);
}
/*! This function constructs the neighbour tree. To this end, we actually need
* to construct the gravitational tree, because we use it now for the
* neighbour search.
*/
void ngb_treebuild(void)
{
if(ThisTask == 0)
printf("Begin Ngb-tree construction.\n");
#if defined(STELLAR_PROP) && defined(CHIMIE) && !defined(CHIMIE_INPUT_ALL)
force_treebuild(N_gas+N_stars);
#else
force_treebuild(N_gas);
#endif
if(ThisTask == 0)
printf("Ngb-Tree contruction finished \n");
}
#ifdef SIDM_ON
/*! This is the same as ngb_treefind_variable except the particle
* type being looked for is SIDM_type not 0 (gas).
* This function returns neighbours with distance <= hsml and returns them in
* Ngblist. Actually, particles in a box of half side length hsml are
* returned, i.e. the reduction to a sphere still needs to be done in the
* calling routine.
*/
int ngb_treefind_variable_SIDM(FLOAT searchcenter[3], FLOAT hsml, int *startnode)
{
int k, numngb;
int no, p;
struct NODE *this;
FLOAT searchmin[3], searchmax[3];
#ifdef PERIODIC
double xtmp;
#endif
for(k = 0; k < 3; k++) /* cube-box window */
{
searchmin[k] = searchcenter[k] - hsml;
searchmax[k] = searchcenter[k] + hsml;
}
numngb = 0;
no = *startnode;
while(no >= 0)
{
if(no < All.MaxPart) /* single particle */
{
p = no;
no = Nextnode[no];
if(P[p].Type != SIDM_type)
continue;
#ifdef PERIODIC
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) < -hsml)
continue;
if(NGB_PERIODIC_X(P[p].Pos[0] - searchcenter[0]) > hsml)
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) < -hsml)
continue;
if(NGB_PERIODIC_Y(P[p].Pos[1] - searchcenter[1]) > hsml)
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) < -hsml)
continue;
if(NGB_PERIODIC_Z(P[p].Pos[2] - searchcenter[2]) > hsml)
continue;
#else
if(P[p].Pos[0] < searchmin[0])
continue;
if(P[p].Pos[0] > searchmax[0])
continue;
if(P[p].Pos[1] < searchmin[1])
continue;
if(P[p].Pos[1] > searchmax[1])
continue;
if(P[p].Pos[2] < searchmin[2])
continue;
if(P[p].Pos[2] > searchmax[2])
continue;
#endif
Ngblist[numngb++] = p;
if(numngb == MAX_NGB)
{
numngb = ngb_clear_buf(searchcenter, hsml, numngb);
if(numngb == MAX_NGB)
{
printf("ThisTask=%d: Need to do a second neighbour loop for (%g|%g|%g) hsml=%g no=%d\n",
ThisTask, searchcenter[0], searchcenter[1], searchcenter[2], hsml, no);
*startnode = no;
return numngb;
}
}
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
Exportflag[DomainTask[no - (All.MaxPart + MaxNodes)]] = 1;
no = Nextnode[no - MaxNodes];
continue;
}
this = &Nodes[no];
no = this->u.d.sibling; /* in case the node can be discarded */
#ifdef PERIODIC
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_X(this->center[0] - searchcenter[0]) - 0.5 * this->len) > hsml)
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_Y(this->center[1] - searchcenter[1]) - 0.5 * this->len) > hsml)
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) + 0.5 * this->len) < -hsml)
continue;
if((NGB_PERIODIC_Z(this->center[2] - searchcenter[2]) - 0.5 * this->len) > hsml)
continue;
#else
if((this->center[0] + 0.5 * this->len) < (searchmin[0]))
continue;
if((this->center[0] - 0.5 * this->len) > (searchmax[0]))
continue;
if((this->center[1] + 0.5 * this->len) < (searchmin[1]))
continue;
if((this->center[1] - 0.5 * this->len) > (searchmax[1]))
continue;
if((this->center[2] + 0.5 * this->len) < (searchmin[2]))
continue;
if((this->center[2] - 0.5 * this->len) > (searchmax[2]))
continue;
#endif
no = this->u.d.nextnode; /* ok, we need to open the node */
}
}
*startnode = -1;
return numngb;
}
#endif
Event Timeline
Log In to Comment