Homec4science
Diffusion Gear (master) Browse

Gear/src/init_sub.c
06d12d3a7f76master

Last Change
Raw File

init_sub.c
View Options

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
#ifdef PY_INTERFACE
#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 function is used to find an initial smoothing length for each SPH
* particle. It guarantees that the number of neighbours will be between
* desired_ngb-MAXDEV and desired_ngb+MAXDEV. For simplicity, a first guess
* of the smoothing length is provided to the function density(), which will
* then iterate if needed to find the right smoothing length.
*/
void setup_smoothinglengths_sub(void)
{
/* before, we can not use the tree, as Q particles do not belong to them... */
//density_sub();
int no, p;
int i,k,startnode;
struct NODE *this;
double old_len,old_mass;
FLOAT searchmin[3], searchmax[3],searchcenter[3];
FLOAT hsml;
#ifdef PERIODIC
double xtmp;
#endif
hsml = 0;
startnode = All.MaxPart;
for (i=0;i < NumPartQ;i++)
{
old_len = MAX_REAL_NUMBER; /* a big number */
old_mass = 0;
//if (i==859)
// printf("doing ... %d %g %g %g\n",i,Q[i].Pos[0],Q[i].Pos[1],Q[i].Pos[2]);
for(k = 0; k < 3; k++) /* cube-box window */
{
searchcenter[k] = Q[i].Pos[k];
searchmin[k] = Q[i].Pos[k] - hsml;
searchmax[k] = Q[i].Pos[k] + hsml;
}
no = startnode;
while(no >= 0)
{
//if (i==859)
// printf("\nno =%d\n",no);
if(no < All.MaxPart) /* single particle */
{
p = no;
//no = Nextnode[no];
/* get attribute from the nearest particle */
SphQ[i].Hsml = SphP[p].Hsml;
SphQ[i].Density = SphP[p].Density;
//if (i==859)
// printf("we find a particle i=%d: %g %g %g -> %g %g %g hsml=%g\n",i,Q[i].Pos[0],Q[i].Pos[1],Q[i].Pos[2],P[p].Pos[0],P[p].Pos[1],P[p].Pos[2],SphQ[i].Hsml);
break;
}
else
{
if(no >= All.MaxPart + MaxNodes) /* pseudo particle */
{
endrun(12321123); /* does not work in parallel */
/* to work in // simply wrap it as done for sph */
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 */
if (this->len > old_len) /* this should not occure !!! */
{
SphQ[i].Hsml = old_len;
SphQ[i].Density = old_mass / pow(old_len,3);
//if (i==859)
//{
//printf("this len %g old_len %g\n",this->len,old_len);
//printf("we find an empty node i=%d: %g %g %g hsml=%g\n",i,Q[i].Pos[0],Q[i].Pos[1],Q[i].Pos[2],SphQ[i].Hsml);
//printf("potential sibling no=%d\n",no);
//printf("current boxsize %g mass = %g\n",this->len,this->u.d.mass);
//printf("min=%g max=%g %g\n",this->center[0] - 0.5 * this->len,this->center[0] + 0.5 * this->len,searchmin[0]);
//printf("min=%g max=%g %g\n",this->center[1] - 0.5 * this->len,this->center[1] + 0.5 * this->len,searchmin[1]);
//printf("min=%g max=%g %g\n",this->center[2] - 0.5 * this->len,this->center[2] + 0.5 * this->len,searchmin[2]);
//}
break;
}
old_len = this->len;
old_mass = this->u.d.mass;
#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 */
//if (i==859)
// printf("try to open next node no=%d\n",no);
}
}
}
}
#endif /*PY_INTERFACE*/
c4science · Help