sandbox/Antoonvh/structurefunction.h
The functions within this file aim to calculate a second-order longitudional structure function from a vector field. We can define a distance vector \mathbf{l} separating two points in space (\mathbf{x_1,x_1+l}). Generally speaking, the vectors \mathbf{u} at these points are differt. The difference in the magnitude of their projections on the \mathbf{l} vector may have some interesting properties. We define:
\displaystyle \delta v_{\|}(\mathbf{x_1,l})= \left( \mathbf{v(x_1+l)-v(x_1)}\right) \cdot \frac{\mathbf{l}}{l},
where l is \|\mathbf{l}\|. By definition, from a statistical perspective, within a homogeneous and isotropic region, one would expect \|\delta v_{\|}(\mathbf{x_1,l})\| only to be dependend on l, i.e. the length of the separation vector. Therefore we can define a second order statistic related to a vector field (\mathbf{u}) according to:
\displaystyle S_2(l)=\langle \delta v_{\|}(\mathbf{x_1,l})^2 \rangle
where the variable \langle x \rangle represents a region-averaged value of a dummy variable x. The code belows aims to evaluate some discretized version of S_2(l) from a vector field on a three-dimensional-tree grid.
int obtain_a_cache(int j,double d[][6],int lev[],vector u){
#if _MPI
MPI_Barrier(MPI_COMM_WORLD);
#endif
srand(time(NULL));// In order to obtain a smooth-ish pair distances distribution we add a stogastic component;
int index[j];
int n=0;
foreach()
n++;
double meaninc=(double)(n)/(double)(j); //The mean increment to obtain approx j cells.
index[0]=0;
int g=1;
while (g<j){
double randsteprel =((sign(noise())*pow((rand()/(double)RAND_MAX),0.2))+1.0);
int inc = (int)((1.0*randsteprel*meaninc)+1.0);
index[g]=index[g-1]+inc;
g++;
}
n=0;
int jj=0;
foreach(){
if(n==index[jj]){
d[jj][0]=x; d[jj][1]=y; d[jj][2]=z;
d[jj][3]=u.x[]; d[jj][4]=u.y[]; d[jj][5]=u.z[];
lev[jj]=level;
jj++;
}
n++;
}
printf("there are %d in cache and %d requested\n",jj,j);
return jj; //Return jj (<=j) the number of cells in the cache
}
void set_g_and_gg(int ggg[2], int j){
int total = j*(j-1)/2;
int start = (pid()*total/npe())+1;
int m = 0;
for (int g=0;g<j;g++){
for (int gg=g+1;gg<j;gg++){
m++;
if (m==start){
ggg[0]=g;
ggg[1]=gg;
printf("%d %d %d %d %d %d\n",ggg[0],ggg[1],total,pid(),start,j);
return;
}
}
}
}
int structure2(double cach[][6],int lev[],int j,double L,int Len,double dvarr[][4]){
int g=0;
int gg=1;
#if _MPI
int ggg[2];
set_g_and_gg(ggg,j);
g=ggg[0];
gg=ggg[1];
#endif
int iter=1;
int todo=((j*(j-1))/2)/npe();
double Del = L/(double)Len;
double l;
while (iter<=todo){
l=pow(sq(cach[g][0]-cach[gg][0])+
sq(cach[g][1]-cach[gg][1])+
sq(cach[g][2]-cach[gg][2]),0.5);
if (l>0.&&l<L0){ //distance within range
double weight = (double)(1<<(3*(2*depth()-lev[gg]-lev[g]))); //weigh with cells' volumes
double dv = (sq((cach[g][3]-cach[gg][3])*(cach[g][0]-cach[gg][0]))+
sq((cach[g][4]-cach[gg][4])*(cach[g][1]-cach[gg][1]))+
sq((cach[g][5]-cach[gg][5])*(cach[g][2]-cach[gg][2])))/(sq(l));
int index = (int)((l/Del)+0.5);
if (index<Len){
dvarr[index][1]+=dv*weight;
dvarr[index][2]+=weight;
dvarr[index][3]+=1.;
}
}
gg++;
if (gg==j){
g++;
gg=g+1;
}
iter++;
}
return 1;
}
void printdvarr(FILE * fp, double dvarr[][4], int Len){
int it = 0;
while (it<Len){
for (int gg=0;gg<4;gg++)
fprintf(fp,"%g\t",dvarr[it][gg]);
fprintf(fp,"\n");
it++;
}
}
#if _MPI
int combine_MPI_caches(double cach[][6],int lev[],int j){
int tot=0;
MPI_Allreduce(&j,&tot,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
double cachsend[j*6];
double * cachrecv= (double *) malloc(sizeof(double)* tot * 6);
int levrecv[tot];
for (int g=0;g<j;g++){
for (int gg=0;gg<6;gg++)
cachsend[6*g+gg]=cach[g][gg];
}
MPI_Allgather(&cachsend[0],6*j,MPI_DOUBLE,cachrecv,6*j,MPI_DOUBLE,MPI_COMM_WORLD);
MPI_Allgather(&lev[0],j,MPI_INT,&levrecv[0],j,MPI_INT,MPI_COMM_WORLD);
for (int g=0;g<tot;g++){
lev[g]=levrecv[g];
for (int gg=0;gg<6;gg++)
cach[g][gg]=cachrecv[6*g+gg];
}
return tot;
}
#endif
//The user interface function;
int long2structure(FILE * fp,vector u,int n,double L,int Len){
double Del = L/(double)Len;
double dvarr[Len][4];
int it=0;
while (it<Len){
for (int gg=0;gg<4;gg++)
dvarr[it][gg]=0.;
it++;
}
double cach[n][6];
int lev[n];
int locpnts = ((n/npe())+0.5);
int j=obtain_a_cache(locpnts,cach,lev,u);
#if _MPI
int d=j;
j=combine_MPI_caches(cach,lev,d);
#endif
structure2(cach,lev,j,L,Len,dvarr);
it=0;
#if _MPI
double dvarrrecv[Len][4]; // An array to store the reduced version of *dvarr*, according to MPI_SUM
MPI_Reduce(&dvarr[0][0],&dvarrrecv[0][0],Len*4,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
if (pid()==0){
while (it<Len){
for (int gg=0;gg<4;gg++)
dvarr[it][gg]=dvarrrecv[it][gg];
it++;
}
it=0;
#endif
double aa;
while (it<Len){
if (dvarr[it][2]>0.){
aa = dvarr[it][1]/dvarr[it][2];
dvarr[it][1] = aa;
}
dvarr[it][0]=it*Del;
it++;
}
printdvarr(fp,dvarr,Len);
#if _MPI
}
#endif
return 1;
}