Changeset b0225a for pcp

Timestamp:

Apr 22, 2008, 10:49:59 AM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

0da6d5

Parents:

2399875

Message:

added LinearInterpolationBetweenGrid(), LinearPullDownFromUpperLevel() and FirstDiscreteDerivative() for future Gradient Correction methods

File:

• pcp/src/perturbed.c (modified) (1 diff)

pcp/src/perturbed.c

@@ -1711 +1711 @@
 }
 
+/** Linear interpolation for coordinate \a R that lies between grid nodes of \a *grid.
+ * \param *P Problem at hand
+ * \param *Lat Lattice structure for grid axis
+ * \param *Lev LatticeLevel structure for grid axis node counts
+ * \param R[] coordinate vector
+ * \param *grid grid with fixed nodes
+ * \return linearly interpolated value of \a *grid for position \a R[NDIM]
+ */
+double LinearInterpolationBetweenGrid(struct Problem *P, struct Lattice *Lat, struct LatticeLevel *Lev, double R[NDIM], fftw_real *grid) 
+{
+  double x[2][NDIM];
+  const int myPE = P->Par.me_comm_ST_Psi;
+  int N[NDIM];
+  const int N0 = Lev->Plan0.plan->local_nx;
+  N[0] = Lev->Plan0.plan->N[0];
+  N[1] = Lev->Plan0.plan->N[1];
+  N[2] = Lev->Plan0.plan->N[2];
+  int g;
+  double n[NDIM];
+  int k[2][NDIM];
+  double sigma;
+  
+  RMat33Vec3(n, Lat->ReciBasis, &R[0]); // transform real coordinates to [0,1]^3 vector
+  for (g=0;g<NDIM;g++) {
+    // k[i] are right and left nearest neighbour node to true position
+    k[0][g] = floor(n[g]/(2.*PI)*(double)N[g]);  // n[2] is floor grid
+    k[1][g] = ceil(n[g]/(2.*PI)*(double)N[g]); // n[1] is ceil grid
+    // x[i] give weights of left and right neighbours (the nearer the true point is to one, the closer its weight to 1) 
+    x[0][g] = (k[1][g] - n[g]/(2.*PI)*(double)N[g]);
+    x[1][g] = 1. - x[0][g];
+    //fprintf(stderr,"(%i) n = %lg, n_floor[%i] = %i\tn_ceil[%i] = %i --- x_floor[%i] = %e\tx_ceil[%i] = %e\n",P->Par.me, n[g], g,k[0][g], g,k[1][g], g,x[0][g], g,x[1][g]);
+  }
+  sigma = 0.;
+  for (g=0;g<2;g++) { // interpolate linearly between adjacent grid points per axis
+    if ((k[g][0] >= N0*myPE) && (k[g][0] < N0*(myPE+1))) {
+      //fprintf(stderr,"(%i) grid[%i]: sigma = %e\n", P->Par.me,  k[g][2]+N[2]*(k[g][1]+N[1]*(k[g][0]-N0*myPE)), sigma);
+      sigma += (x[g][0]*x[0][1]*x[0][2])*grid[k[0][2]+N[2]*(k[0][1]+N[1]*(k[g][0]-N0*myPE))]*mu0; // if it's local and factor from inverse fft
+      //fprintf(stderr,"(%i) grid[%i]: sigma += %e * %e \n", P->Par.me,  k[g][2]+N[2]*(k[g][1]+N[1]*(k[g][0]-N0*myPE)), (x[g][0]*x[0][1]*x[0][2]), grid[k[0][2]+N[2]*(k[0][1]+N[1]*(k[g][0]-N0*myPE))]*mu0);
+      sigma += (x[g][0]*x[0][1]*x[1][2])*grid[k[1][2]+N[2]*(k[0][1]+N[1]*(k[g][0]-N0*myPE))]*mu0; // if it's local and factor from inverse fft
+      //fprintf(stderr,"(%i) grid[%i]: sigma += %e * %e \n", P->Par.me,  k[g][2]+N[2]*(k[g][1]+N[1]*(k[g][0]-N0*myPE)), (x[g][0]*x[0][1]*x[1][2]), grid[k[1][2]+N[2]*(k[0][1]+N[1]*(k[g][0]-N0*myPE))]*mu0);
+      sigma += (x[g][0]*x[1][1]*x[0][2])*grid[k[0][2]+N[2]*(k[1][1]+N[1]*(k[g][0]-N0*myPE))]*mu0; // if it's local and factor from inverse fft
+      //fprintf(stderr,"(%i) grid[%i]: sigma += %e * %e \n", P->Par.me,  k[g][2]+N[2]*(k[g][1]+N[1]*(k[g][0]-N0*myPE)), (x[g][0]*x[1][1]*x[0][2]), grid[k[0][2]+N[2]*(k[1][1]+N[1]*(k[g][0]-N0*myPE))]*mu0);
+      sigma += (x[g][0]*x[1][1]*x[1][2])*grid[k[1][2]+N[2]*(k[1][1]+N[1]*(k[g][0]-N0*myPE))]*mu0; // if it's local and factor from inverse fft
+      //fprintf(stderr,"(%i) grid[%i]: sigma += %e * %e \n", P->Par.me,  k[g][2]+N[2]*(k[g][1]+N[1]*(k[g][0]-N0*myPE)), (x[g][0]*x[1][1]*x[1][2]), grid[k[1][2]+N[2]*(k[1][1]+N[1]*(k[g][0]-N0*myPE))]*mu0);
+    }
+  }
+  return sigma;
+}
+
+/** Linear Interpolation from all eight corners of the box that singles down to a point on the lower level.
+ * \param *P Problem at hand
+ * \param *Lev LatticeLevel structure for node numbers
+ * \param upperNode Node around which to interpolate
+ * \param *upperGrid array of grid points
+ * \return summed up and then averaged octant around \a upperNode 
+ */  
+double LinearPullDownFromUpperLevel(struct Problem *P, struct LatticeLevel *Lev, int upperNode, fftw_real *upperGrid) 
+{
+  const int N0 = Lev->Plan0.plan->local_nx;
+  const int N1 = Lev->Plan0.plan->N[1];
+  const int N2 = Lev->Plan0.plan->N[2];
+  double lowerGrid = 0.;
+  int nr=1;
+  lowerGrid += upperGrid[upperNode];
+  if (upperNode % N0 != N0-1) {
+    lowerGrid += upperGrid[upperNode+1];
+    nr++;
+    if (upperNode % N1 != N1-1) {
+      lowerGrid += upperGrid[upperNode + 0 + N2*(1 + N1*1)];
+      nr++;
+      if (upperNode % N2 != N2-1) {
+       lowerGrid += upperGrid[upperNode + 1 + N2*(1 + N1*1)];
+        nr++;
+      }
+    }
+    if (upperNode % N2 != N2-1) {
+      lowerGrid += upperGrid[upperNode + 1 + N2*(0 + N1*1)];
+      nr++;
+    }
+  }
+  if (upperNode % N1 != N1-1) {
+    lowerGrid += upperGrid[upperNode + 0 + N2*(1 + N1*0)];
+    nr++;
+    if (upperNode % N2 != N2-1) {
+      lowerGrid += upperGrid[upperNode + 1 + N2*(1 + N1*0)];
+      nr++;
+    }
+  }
+  if (upperNode % N2 != N2-1) {
+    lowerGrid += upperGrid[upperNode + 1 + N2*(0 + N1*0)];
+    nr++;
+  }
+  return (lowerGrid/(double)nr);
+} 
+
+/** Evaluates the 1-stern in order to evaluate the first derivative on the grid.
+ * \param *P Problem at hand
+ * \param *Lev Level to interpret the \a *density on
+ * \param *density array with gridded values
+ * \param *n 3 vector with indices on the grid
+ * \param axis axis along which is derived
+ * \param myPE number of processes who share the density
+ * \return  [+1/2 -1/2] of \a *n
+ */
+double FirstDiscreteDerivative(struct Problem *P, struct LatticeLevel *Lev, fftw_real *density, int *n, int axis, int myPE)
+{
+  int *N = Lev->Plan0.plan->N;  // maximum nodes per axis
+  const int N0 = Lev->Plan0.plan->local_nx; // special local number due to parallel split up
+  double ret[NDIM],  Ret[NDIM];  // return value local/global
+  int i;
+
+  for (i=0;i<NDIM;i++) {
+    ret[i] = Ret[i] = 0.;
+  }
+  if (((n[0]+1)%N[0] >= N0*myPE) && ((n[0]+1)%N[0] < N0*(myPE+1)))   // next cell belongs to this process 
+    ret[0] += 1./2. * (density[n[2]+N[2]*(n[1]+N[1]*(n[0]+1-N0*myPE))]);
+  if (((n[0]-1)%N[0] >= N0*myPE) && ((n[0]-1)%N[0] < N0*(myPE+1))) // previous cell belongs to this process 
+    ret[0] -= 1./2. * (density[n[2]+N[2]*(n[1]+N[1]*(n[0]-1-N0*myPE))]);
+  if ((n[0] >= N0*myPE) && (n[0] < N0*(myPE+1))) { 
+    ret[1] += 1./2. * (density[n[2]+N[2]*((n[1]+1)%N[1] + N[1]*(n[0]%N0))]);
+    ret[1] -= 1./2. * (density[n[2]+N[2]*((n[1]-1)%N[1] + N[1]*(n[0]%N0))]);
+  }     
+  if ((n[0] >= N0*myPE) && (n[0] < N0*(myPE+1))) {
+    ret[2] += 1./2. * (density[(n[2]+1)%N[2] + N[2]*(n[1]+N[1]*(n[0]%N0))]);
+    ret[2] -= 1./2. * (density[(n[2]-1)%N[2] + N[2]*(n[1]+N[1]*(n[0]%N0))]);
+  }
+  
+  if (MPI_Allreduce(ret, Ret, 3, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi) != MPI_SUCCESS)
+    Error(SomeError, "FirstDiscreteDerivative: MPI_Allreduce failure!");
+   
+  for (i=0;i<NDIM;i++)  // transform from node count to [0,1]^3
+    Ret[i] *= N[i]; 
+  RMat33Vec3(ret, P->Lat.ReciBasis, Ret); // this actually divides it by mesh length in real coordinates
+  //fprintf(stderr, "(%i) sum at (%i,%i,%i) : %lg\n",P->Par.me, n[0],n[1],n[2], ret[axis]);  
+  return ret[axis]; ///(P->Lat.RealBasisQ[axis]/N[axis]);
+} 
 
 /** Fouriertransforms given \a source.