Changeset b924cd for pcp

Timestamp:

Apr 22, 2008, 12:07:53 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

986488

Parents:

60a9f9

Message:

CalculateChemicalShieldingByReciprocalCurrentDensity(): Linear interpolation to core position R shifted to function LinearInterpolationBetweenGrid()

new function LinearInterpolationBetweenGrid() which does the linear interpolation.

File:

• pcp/src/perturbed.c (modified) (2 diffs)

pcp/src/perturbed.c

@@ -3058 +3058 @@
   fftw_real *sigma_real[NDIM_NDIM];
   double sigma,Sigma;
-  double x[2][NDIM];
   int it, ion, in, dex, g, Index, i;
   int *N = Lev0->Plan0.plan->N;
   //const double FFTfactor = 1.;///Lev0->MaxN;
-  int n[2][NDIM];
-  double eta, delta_sigma, S, A, iso, tmp;
+  double eta, delta_sigma, S, A, iso;
   FILE *SigmaFile;
   char suffixsigma[255];
-  const int myPE = P->Par.me_comm_ST_Psi;
   int cross_lookup[4]; // cross lookup table
-  const int N0 = Lev0->Plan0.plan->local_nx;
   const double factorDC = R->FactorDensityC;
   gsl_matrix_complex *H = gsl_matrix_complex_calloc(NDIM,NDIM);
@@ -3144 +3140 @@
         for (ion=0; ion < I->I[it].Max_IonsOfType; ion++) { // and each ion of type
           // read transformed sigma at core position and MPI_Allreduce
-          for (g=0;g<NDIM;g++) {
-            n[0][g] = floor(I->I[it].R[NDIM*ion+g]/sqrt(Lat->RealBasisSQ[g])*(double)N[g]);
-            n[1][g] = ceil(I->I[it].R[NDIM*ion+g]/sqrt(Lat->RealBasisSQ[g])*(double)N[g]);
-            x[1][g] = I->I[it].R[NDIM*ion+g]/sqrt(Lat->RealBasisSQ[g])*(double)N[g] - (double)n[0][g];
-            x[0][g] = 1. - x[1][g];
-            //fprintf(stderr,"(%i) n_floor[%i] = %i\tn_ceil[%i] = %i --- x_floor[%i] = %e\tx_ceil[%i] = %e\n",P->Par.me, g,n[0][g], g,n[1][g], g,x[0][g], g,x[1][g]);
-          }
-          sigma = 0.;
-          for (i=0;i<2;i++) { // interpolate linearly between adjacent grid points per axis
-            if ((n[i][0] >= N0*myPE) && (n[i][0] < N0*(myPE+1))) {
-//              fprintf(stderr,"(%i) field[%i]: sigma = %e\n", P->Par.me,  n[i][2]+N[2]*(n[i][1]+N[1]*(n[i][0]-N0*myPE)), sigma);
-              sigma += (x[i][0]*x[0][1]*x[0][2])*sigma_real[in+dex*NDIM][n[0][2]+N[2]*(n[0][1]+N[1]*(n[i][0]-N0*myPE))]*mu0; // if it's local and factor from inverse fft
-              //fprintf(stderr,"(%i) field[%i]: sigma += %e * %e \n", P->Par.me,  n[i][2]+N[2]*(n[i][1]+N[1]*(n[i][0]-N0*myPE)), (x[i][0]*x[0][1]*x[0][2]), field[n[0][2]+N[2]*(n[0][1]+N[1]*(n[i][0]-N0*myPE))]*mu0);
-              sigma += (x[i][0]*x[0][1]*x[1][2])*sigma_real[in+dex*NDIM][n[1][2]+N[2]*(n[0][1]+N[1]*(n[i][0]-N0*myPE))]*mu0; // if it's local and factor from inverse fft
-              //fprintf(stderr,"(%i) field[%i]: sigma += %e * %e \n", P->Par.me,  n[i][2]+N[2]*(n[i][1]+N[1]*(n[i][0]-N0*myPE)), (x[i][0]*x[0][1]*x[1][2]), field[n[1][2]+N[2]*(n[0][1]+N[1]*(n[i][0]-N0*myPE))]*mu0);
-              sigma += (x[i][0]*x[1][1]*x[0][2])*sigma_real[in+dex*NDIM][n[0][2]+N[2]*(n[1][1]+N[1]*(n[i][0]-N0*myPE))]*mu0; // if it's local and factor from inverse fft
-              //fprintf(stderr,"(%i) field[%i]: sigma += %e * %e \n", P->Par.me,  n[i][2]+N[2]*(n[i][1]+N[1]*(n[i][0]-N0*myPE)), (x[i][0]*x[1][1]*x[0][2]), field[n[0][2]+N[2]*(n[1][1]+N[1]*(n[i][0]-N0*myPE))]*mu0);
-              sigma += (x[i][0]*x[1][1]*x[1][2])*sigma_real[in+dex*NDIM][n[1][2]+N[2]*(n[1][1]+N[1]*(n[i][0]-N0*myPE))]*mu0; // if it's local and factor from inverse fft
-              //fprintf(stderr,"(%i) field[%i]: sigma += %e * %e \n", P->Par.me,  n[i][2]+N[2]*(n[i][1]+N[1]*(n[i][0]-N0*myPE)), (x[i][0]*x[1][1]*x[1][2]), field[n[1][2]+N[2]*(n[1][1]+N[1]*(n[i][0]-N0*myPE))]*mu0);
-            }
-          }
-          sigma *= -R->FactorDensityR; // factor from inverse fft? (and its defined as negative proportionaly factor)
+          sigma = -LinearInterpolationBetweenGrid(P, Lat, Lev0, &I->I[it].R[NDIM*ion], sigma_real[in+dex*NDIM]) * R->FactorDensityR; // factor from inverse fft
+
           MPI_Allreduce ( &sigma, &Sigma, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);  // sum local to total
           I->I[it].sigma_rezi[ion][in+dex*NDIM] = Sigma;