Commits

Anonymous committed 998dfe8

Add serial_coulomb.cxx

Comments (0)

Files changed (2)

wrapper/serial_coulomb.cxx

+#include "serialfmm.h"
+
+extern "C" void FMMcalccoulomb_ij(int ni, double* xi, double* qi, double* fi,
+  int nj, double* xj, double* qj, double, int tblno, double size, int periodicflag) {
+  std::cout << "tblno: " << tblno << std::endl;
+  IMAGES = ((periodicflag & 0x1) == 0) ? 0 : 3;
+  THETA = .5;
+  Bodies bodies(ni),jbodies(nj);
+  Cells cells,jcells;
+  SerialFMM<Laplace> FMM;
+
+  for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) {
+    int i = B-bodies.begin();
+    B->X[0] = xi[3*i+0];
+    B->X[1] = xi[3*i+1];
+    B->X[2] = xi[3*i+2];
+
+
+    if( B->X[0] < -size/2 ) B->X[0] += size;
+    if( B->X[1] < -size/2 ) B->X[1] += size;
+    if( B->X[2] < -size/2 ) B->X[2] += size;
+    if( B->X[0] >  size/2 ) B->X[0] -= size;
+    if( B->X[1] >  size/2 ) B->X[1] -= size;
+    if( B->X[2] >  size/2 ) B->X[2] -= size;
+    B->SRC  = qi[i];
+    switch (tblno) {
+    case 0 :
+      B->TRG[1] = -fi[3*i+0];
+      B->TRG[2] = -fi[3*i+1];
+      B->TRG[3] = -fi[3*i+2];
+      break;
+    case 1 :
+      B->TRG[0] = fi[3*i+0];
+      break;
+    }
+    B->IBODY = i;
+    B->IPROC = MPIRANK;
+  }
+//    printf("%f %f %f\n", xi[0], xi[1], xi[2]);
+//    printf("%f %f %f\n", xi[96], xi[97], xi[98]);
+//    printf("%f %f %f\n", xi[240], xi[241], xi[242]);
+
+  for( B_iter B=jbodies.begin(); B!=jbodies.end(); ++B ) {
+    int i = B-jbodies.begin();
+    B->X[0] = xj[3*i+0];
+    B->X[1] = xj[3*i+1];
+    B->X[2] = xj[3*i+2];
+    if( B->X[0] < -size/2 ) B->X[0] += size;
+    if( B->X[1] < -size/2 ) B->X[1] += size;
+    if( B->X[2] < -size/2 ) B->X[2] += size;
+    if( B->X[0] >  size/2 ) B->X[0] -= size;
+    if( B->X[1] >  size/2 ) B->X[1] -= size;
+    if( B->X[2] >  size/2 ) B->X[2] -= size;
+    B->SRC  = qj[i];
+  }
+
+  FMM.initialize();
+  FMM.setDomain(bodies,0,size/2);
+  FMM.bottomup(bodies,cells);
+  FMM.bottomup(jbodies,jcells);
+
+  FMM.downward(cells,jcells);
+  std::sort(bodies.begin(),bodies.end());
+  FMM.finalize();
+
+#if 1
+   std::cout<< "check 1" <<std::endl;
+  //for( B_iter B=bodies.begin(); B!=bodies.end()-ni+300; ++B ) {
+  for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) {
+    int i = B-bodies.begin();
+
+ //   std::cout << xi[3*i+0] << " " << B->X[0] << std::endl;
+    //assert(fabs(xi[3*i+0] - B->X[0])< 1e-6);
+    //xi[3*i+0] = B->X[0];
+   // xi[3*i+1] = B->X[1];
+    //xi[3*i+2] = B->X[2];
+
+    qi[i]     = B->SRC;
+    switch (tblno) {
+    case 0 :
+      fi[3*i+0] = -B->SRC * B->TRG[1];
+      fi[3*i+1] = -B->SRC * B->TRG[2];
+      fi[3*i+2] = -B->SRC * B->TRG[3];
+      break;
+    case 1 :
+      fi[3*i+0] = B->SRC * B->TRG[0];
+      break;
+    }
+  }
+
+//    printf("%f %f %f\n", xi[0], xi[1], xi[2]);
+//    printf("%f %f %f\n", xi[96], xi[97], xi[98]);
+//   printf("%f %f %f\n", xi[240], xi[241], xi[242]);
+
+  double fc[3];
+  for( int d=0; d!=3; ++d ) fc[d]=0;
+  for( int i=0; i!=ni; ++i ) { 
+    for( int d=0; d!=3; ++d ) { 
+      fc[d] += qi[i] * xi[3*i+d];
+    }   
+  }
+  if( (tblno % 2) == 1 ) { 
+    for( int i=0; i!=ni; ++i ) { 
+      fi[3*i+0] += 2.0 * M_PI / (3.0 * size * size * size)
+                * (fc[0] * fc[0] + fc[1] * fc[1] + fc[2] * fc[2]) / ni;
+    }   
+  } else {
+    for( int i=0; i!=ni; ++i ) { 
+      for( int d=0; d!=3; ++d ) {
+        fi[3*i+d] -= 4.0 * M_PI * qi[i] * fc[d] / (3.0 * size * size * size);
+      }
+    }   
+  }
+#else
+   std::cout<< "check not 1" <<std::endl;
+  switch (tblno) {
+  case 0 :
+    for( int i=0; i!=ni; ++i ) {
+      double Fx = 0, Fy = 0, Fz = 0;
+      for( int j=0; j!=nj; ++j ) {
+        double dx = xi[3*i+0] - xj[3*j+0];
+        double dy = xi[3*i+1] - xj[3*j+1];
+        double dz = xi[3*i+2] - xj[3*j+2];
+        double R2 = dx * dx + dy * dy + dz * dz;
+        double invR = 1 / std::sqrt(R2);
+        if( R2 == 0 ) invR = 0;
+        double invR3 = qj[j] * invR * invR * invR;
+        Fx += dx * invR3;
+        Fy += dy * invR3;
+        Fz += dz * invR3;
+      }
+      fi[3*i+0] += Fx;
+      fi[3*i+1] += Fy;
+      fi[3*i+2] += Fz;
+    }
+    break;
+  case 1:
+    for( int i=0; i!=ni; ++i ) {
+      double Po = 0;
+      for( int j=0; j!=nj; ++j ) {
+        double dx = xi[3*i+0] - xj[3*j+0];
+        double dy = xi[3*i+1] - xj[3*j+1];
+        double dz = xi[3*i+2] - xj[3*j+2];
+        double R2 = dx * dx + dy * dy + dz * dz;
+        double invR = 1 / std::sqrt(R2);
+        if( R2 == 0 ) invR = 0;
+        Po += qj[j] * invR;
+      } 
+      fi[3*i+0] += Po;
+   //   std::cout << i << " " << fi[3*i+0] << std::endl;
+    }
+    break;
+  } 
+#endif 
+}
+
+extern "C" void getaccfmm_(int *ni, double* xi, double* qi, double* fi,
+  int *nj, double* xj, double* qj, double *rscale, int *tblno, double *size, int *periodicflag) {
+  std::cout << "Starting FMM" << std::endl;
+  // std::cout<< "check" <<std::endl;
+  FMMcalccoulomb_ij(*ni,xi,qi,fi,*nj,xj,qj,*rscale,*tblno,*size,*periodicflag);
+}

wrapper/test_serial_coulomb.cxx

+#include <cmath>
+#include <cstdlib>
+#include <iostream>
+#include "mr3.h"
+
+const double R2MIN = 0.25;
+const double R2MAX = 64;
+
+extern "C" void FMMcalccoulomb_ij(int ni, double* xi, double* qi, double* fi,
+  int nj, double* xj, double* qj, double rscale, int tblno, double size, int periodicflag);
+
+int main() {
+  const int N = 10000;
+  const double size = 2;
+  double *xi     = new double [3*N];
+  double *qi     = new double [N];
+  double *pi     = new double [3*N];
+  double *fi     = new double [3*N];
+  double *pd     = new double [3*N];
+  double *fd     = new double [3*N];
+  double *xj     = new double [3*N];
+  double *qj     = new double [N];
+
+  srand48(0);
+  float average = 0;
+  for( int i=0; i!=N; ++i ) {
+    xi[3*i+0] = drand48() * size - size/2;
+    xi[3*i+1] = drand48() * size - size/2;
+    xi[3*i+2] = drand48() * size - size/2;
+    qi[i] = drand48()*2.0-1.0;
+    average += qi[i];
+    pi[3*i+0] = pi[3*i+1] = pi[3*i+2] = 0;
+    fi[3*i+0] = fi[3*i+1] = fi[3*i+2] = 0;
+    pd[3*i+0] = pd[3*i+1] = pd[3*i+2] = 0;
+    fd[3*i+0] = fd[3*i+1] = fd[3*i+2] = 0;
+  }
+  average /= N;
+  for( int i=0; i!=N; ++i ) {
+    qi[i] -= average;
+  }
+  average = 0;
+  for( int i=0; i!=N; ++i ) {
+    xj[3*i+0] = drand48() * size - size/2;
+    xj[3*i+1] = drand48() * size - size/2;
+    xj[3*i+2] = drand48() * size - size/2;
+    qj[i] = drand48()*2.0-1.0;
+  }
+  average /= N;
+  for( int i=0; i!=N; ++i ) {
+    qj[i] -= average;
+  }
+
+  FMMcalccoulomb_ij(N, xi, qi, pi, N, xj, qj, 0.0, 1, size, 0);
+  FMMcalccoulomb_ij(N, xi, qi, fi, N, xj, qj, 0.0, 0, size, 0);
+#if 0
+  MR3calccoulomb_ij(N, xi, qi, pd, N, xj, qj, 1.0, 1, size, 0);
+  MR3calccoulomb_ij(N, xi, qi, fd, N, xj, qj, 1.0, 0, size, 0);
+#else
+  for( int i=0; i!=N; ++i ) {
+    double P = 0, Fx = 0, Fy = 0, Fz = 0;
+    for( int j=0; j!=N; ++j ) {
+      double dx = xi[3*i+0] - xj[3*j+0];
+      double dy = xi[3*i+1] - xj[3*j+1];
+      double dz = xi[3*i+2] - xj[3*j+2];
+      double R2 = dx * dx + dy * dy + dz * dz;
+      double invR = 1 / std::sqrt(R2);
+      if( R2 == 0 ) invR = 0;
+      double invR3 = qj[j] * invR * invR * invR;
+      P += qj[j] * invR;
+      Fx += dx * invR3;
+      Fy += dy * invR3;
+      Fz += dz * invR3;
+    }
+    pd[3*i+0] += P;
+    fd[3*i+0] += Fx;
+    fd[3*i+1] += Fy;
+    fd[3*i+2] += Fz;
+  }
+#endif
+  double Pd = 0, Pn = 0, Fd = 0, Fn = 0;
+  float sums = 0;
+  for( int i=0; i!=N; ++i ) {
+    Pd += (pi[3*i+0] - pd[3*i+0]) * (pi[3*i+0] - pd[3*i+0]);
+    Pn += pd[3*i+0] * pd[3*i+0];
+    Fd += (fi[3*i+0] - fd[3*i+0]) * (fi[3*i+0] - fd[3*i+0])
+        + (fi[3*i+1] - fd[3*i+1]) * (fi[3*i+1] - fd[3*i+1])
+        + (fi[3*i+2] - fd[3*i+2]) * (fi[3*i+2] - fd[3*i+2]);
+    Fn += fd[3*i+0] * fd[3*i+0] + fd[3*i+1] * fd[3*i+1] + fd[3*i+2] * fd[3*i+2];
+    sums += qi[i]*fd[3*i+0];
+  }
+  std::cout << "Coulomb       potential : " << sqrtf(Pd/Pn) << std::endl;
+  std::cout << "Coulomb       force     : " << sqrtf(Fd/Fn) << std::endl;
+  std::cout << sums << std::endl;
+
+  delete[] xi;
+  delete[] qi;
+  delete[] pi;
+  delete[] fi;
+  delete[] pd;
+  delete[] fd;
+  delete[] xj;
+  delete[] qj;
+}