sandbox/Antoonvh/structure_function.h

    Second order structure function

    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,\ \mathbf{x}_1+\mathbf{l}). Generally speaking, the vectors \mathbf{u} at these points are different. 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,\mathbf{l})= \left( \mathbf{v}(\mathbf{x}_1)-\mathbf{v}(\mathbf{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,\mathbf{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.

    There are many sub-functions so that maybe later, it may be modified/extended to e.g. 2D and 1D.

    First, a function is defined that selects j cells for analysis and stores the relevant data in an array.

    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
    }

    Since this function is not necessarily cheap, each thread does part of the calculations. This function sets the variables that determine what point combinations need to be assessed by the specific thread.

    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;
          }
        }
      }
    }

    This function actually calculates the the structure function for each thread.

    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;
    }

    This function can be used to print the result to a file

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

    This function comminucates the points that each thread found in the obtain_a_cache function. Note that this requires all-to-all communication, a feature that is warranted by the non-local nature of the statistic.

    #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

    User interface function

    The function requires a file pointer to write the result, the vector field, the number of cells that will be used (n, so that there are \mathcal{O}(n^2) queries), the maximum separation distance, and the number of bins for the discretized represenation of l.

    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; 
    }

    Usage