Commits

Rio Yokota committed 9fc5d0d

Merged wrapper/ewald.cxx into charmm.cxx.

Comments (0)

Files changed (5)

examples/Makefile

 # Ewald vs. periodic FMM
 ewald: ewald.o $(OBJECTS)
 	$(CXX) $? $(LFLAGS)
-	mpirun -np 2 ./a.out --numBodies 1000 --ncrit 32 --images 3
+	mpirun -np 2 ./a.out --numBodies 500 --ncrit 32 --images 3

examples/ewald.cxx

   logger.startTimer("Total FMM");
   logger.startPAPI();
   Bodies bodies = data.initBodies(args.numBodies, args.distribution, LET.mpirank, LET.mpisize, 3-LET.mpisize);
+  //data.writeSources(bodies, LET.mpirank);
   Bounds localBounds = boundbox.getBounds(bodies);
   Bounds globalBounds = LET.allreduceBounds(localBounds);
   localBounds = LET.partition(bodies, globalBounds);
 #if 1
   Bodies jbodies = bodies;
   for (int i=0; i<LET.mpisize; i++) {
+    if (args.verbose) std::cout << "Ewald loop           : " << i+1 << "/" << LET.mpisize << std::endl;
     LET.shiftBodies(jbodies);
-    if (args.verbose) std::cout << "Ewald loop           : " << i+1 << "/" << LET.mpisize << std::endl;
     localBounds = boundbox.getBounds(jbodies);
     Cells jcells = tree.buildTree(jbodies, localBounds);
     ewald.wavePart(bodies, jbodies);

include/dataset.h

     return bodies;                                              // Return bodies
   }
 
+//! Read source values from file
+  void readSources(Bodies &bodies, int mpirank) {
+    std::stringstream name;                                     // File name
+    name << "source" << std::setfill('0') << std::setw(4)       // Set format
+         << mpirank << ".dat";                                  // Create file name
+    std::ifstream file(name.str().c_str(),std::ios::in);        // Open file
+    file.seekg(filePosition);                                   // Set position in file
+    for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {       // Loop over bodies
+      file >> B->X[0];                                          //  Read data for x coordinates
+      file >> B->X[1];                                          //  Read data for y coordinates
+      file >> B->X[2];                                          //  Read data for z coordinates
+      file >> B->SRC;                                           //  Read data for charge
+    }                                                           // End loop over bodies
+    filePosition = file.tellg();                                // Get position in file
+    file.close();                                               // Close file
+  }
+
+//! Write source values to file
+  void writeSources(Bodies &bodies, int mpirank) {
+    std::stringstream name;                                     // File name
+    name << "source" << std::setfill('0') << std::setw(4)       // Set format
+         << mpirank << ".dat";                                  // Create file name
+    std::ofstream file(name.str().c_str(),std::ios::out);       // Open file
+    for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {       // Loop over bodies
+      file << B->X[0] << std::endl;                             //  Write data for x coordinates
+      file << B->X[1] << std::endl;                             //  Write data for y coordinates
+      file << B->X[2] << std::endl;                             //  Write data for z coordinates
+      file << B->SRC  << std::endl;                             //  Write data for charge
+    }                                                           // End loop over bodies
+    file.close();                                               // Close file
+  }
+
 //! Read target values from file
-  void readTarget(Bodies &bodies, int mpirank) {
+  void readTargets(Bodies &bodies, int mpirank) {
     std::stringstream name;                                     // File name
-    name << "direct" << std::setfill('0') << std::setw(4)       // Set format
-         << mpirank << ".dat";                                  // Create file name for saving direct calculation values
-    std::ifstream file(name.str().c_str(),std::ios::in | std::ios::binary);// Open file
+    name << "target" << std::setfill('0') << std::setw(4)       // Set format
+         << mpirank << ".dat";                                  // Create file name
+    std::ifstream file(name.str().c_str(),std::ios::in);        // Open file
     file.seekg(filePosition);                                   // Set position in file
     for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {       // Loop over bodies
       file >> B->TRG[0];                                        //  Read data for potential
   }
 
 //! Write target values to file
-  void writeTarget(Bodies &bodies, int mpirank) {
+  void writeTargets(Bodies &bodies, int mpirank) {
     std::stringstream name;                                     // File name
-    name << "direct" << std::setfill('0') << std::setw(4)       // Set format
-         << mpirank << ".dat";                                  // Create file name for saving direct calculation values
-    std::ofstream file(name.str().c_str(),std::ios::out | std::ios::app | std::ios::binary);// Open file
+    name << "target" << std::setfill('0') << std::setw(4)       // Set format
+         << mpirank << ".dat";                                  // Create file name
+    std::ofstream file(name.str().c_str(),std::ios::out);       // Open file
     for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {       // Loop over bodies
       file << B->TRG[0] << std::endl;                           //  Write data for potential
       file << B->TRG[1] << std::endl;                           //  Write data for x acceleration

wrappers/charmm.cxx

 #include "args.h"
 #include "boundbox.h"
 #include "buildtree.h"
+#include "ewald.h"
 #include "sort.h"
 #include "traversal.h"
 #include "updownpass.h"
 LocalEssentialTree *LET;
 
 extern "C" void fmm_init_(int * images) {
-  const int ncrit = 16;
+  const int ncrit = 32;
   const int nspawn = 1000;
   const real_t theta = 0.4;
   args = new Args;
   bodies = LET->commBodies(bodies);
   Cells cells = tree->buildTree(bodies, localBounds);
   pass->upwardPass(cells);
-
   for (int i=0; i<*nglobal; i++) {
     icpumap[i] = 0;
   }
     f[3*i+2] = B->TRG[3];
   }
 }
+
+extern "C" void ewald_(int * nglobal, int * icpumap, double * x, double * q, double * p, double * f,
+                       int * ksize, double * alpha, double * sigma, double * cutoff, double * cycle) {
+  Ewald * ewald = new Ewald(*ksize, *alpha, *sigma, *cutoff, *cycle);
+  if (args->verbose) ewald->verbose = true;
+  int nlocal = 0;
+  for (int i=0; i<*nglobal; i++) {
+    if (icpumap[i] == 1) nlocal++;
+  }
+  args->numBodies = nlocal;
+  logger->printTitle("Ewald Parameters");
+  args->print(logger->stringLength, P, LET->mpirank);
+  ewald->print(logger->stringLength);
+#if _OPENMP
+#pragma omp parallel
+#pragma omp master
+#endif
+  logger->printTitle("Ewald Profiling");
+  logger->startTimer("Total Ewald");
+  logger->startPAPI();
+  Bodies bodies(nlocal);
+  B_iter B = bodies.begin();
+  for (int i=0; i<*nglobal; i++) {
+    if (icpumap[i] == 1) {
+      B->X[0] = x[3*i+0];
+      B->X[1] = x[3*i+1];
+      B->X[2] = x[3*i+2];
+      B->SRC = q[i];
+      B->TRG[0] = p[i];
+      B->TRG[1] = f[3*i+0];
+      B->TRG[2] = f[3*i+1];
+      B->TRG[3] = f[3*i+2];
+      B->IBODY = i;
+      B++;
+    }
+  }
+  for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {
+    if( B->X[0] < -*cycle/2 ) B->X[0] += *cycle;
+    if( B->X[1] < -*cycle/2 ) B->X[1] += *cycle;
+    if( B->X[2] < -*cycle/2 ) B->X[2] += *cycle;
+    if( B->X[0] >  *cycle/2 ) B->X[0] -= *cycle;
+    if( B->X[1] >  *cycle/2 ) B->X[1] -= *cycle;
+    if( B->X[2] >  *cycle/2 ) B->X[2] -= *cycle;
+  }
+  Cells cells = tree->buildTree(bodies, localBounds);
+  Bodies jbodies = bodies;
+  for (int i=0; i<LET->mpisize; i++) {
+    if (args->verbose) std::cout << "Ewald loop           : " << i+1 << "/" << LET->mpisize << std::endl;
+    LET->shiftBodies(jbodies);
+    localBounds = boundbox->getBounds(jbodies);
+    Cells jcells = tree->buildTree(jbodies, localBounds);
+    ewald->wavePart(bodies, jbodies);
+    ewald->realPart(cells, jcells);
+  }
+  ewald->selfTerm(bodies);
+  logger->stopPAPI();
+  logger->stopTimer("Total Ewald");
+  logger->printTitle("Total runtime");
+  logger->printTime("Total Ewald");
+
+  for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {
+    int i = B->IBODY;
+    p[i]     = B->TRG[0];
+    f[3*i+0] = B->TRG[1];
+    f[3*i+1] = B->TRG[2];
+    f[3*i+2] = B->TRG[3];
+  }
+}

wrappers/test_charmm.f90

       implicit real*8(a-h,o-z)
       include 'mpif.h'
       parameter(nmax = 1000000, pi = 3.14159265358979312d0)
+      character(128) filename
       integer prange, d
       real(8) norm, norm1, norm2, norm3, norm4
       type(c_ptr) ctx
       call mpi_comm_size(mpi_comm_world, mpisize, ierr);
       call mpi_comm_rank(mpi_comm_world, mpirank, ierr);
       nglobal = 1000;
-      images = 2
+      images = 3
       ksize = 11
       pcycle = 2 * pi
       alpha = 10 / pcycle
+      sigma = .25 / pi;
+      cutoff = pcycle * alpha / 3;
       allocate( x(3*nmax),  q(nmax),  p(nmax),  f(3*nmax), icpumap(nmax) )
       allocate( x2(3*nmax), q2(nmax), p2(nmax), f2(3*nmax) )
-
+#if 0
       do i = 1, 128
         iseed(i) = 0
       end do
       do i = 1, nglobal
         q(i) = q(i) - average
       end do
+#else
+      write (filename, '("source", i4.4, ".dat")') 0
+      open(10, file=filename)
+      do i = 1, nglobal
+        read(10,*) x(3*i-2), x(3*i-1), x(3*i-0), q(i)
+        p(i) = 0
+        f(3*i-2) = 0
+        f(3*i-1) = 0
+        f(3*i-0) = 0
+        icpumap(i) = 0
+      end do
+#endif
       ista = 1
       iend = nglobal
       call split_range(ista,iend,mpirank,mpisize)
         f2(3*i-1) = 0
         f2(3*i-0) = 0
       end do
-#if 0
-      call fmm_ewald(nglobal, x2, q2, p2, f2, ksize, alpha, pcycle)
+#if 1
+      call ewald(nglobal, icpumap, x2, q2, p2, f2, ksize, alpha, sigma, cutoff, pcycle)
 #else
       prange = 0
       do i = 0, images-1
                   dz = x(3*i-0) - x2(3*j-0) - xperiodic(3)
                   R2 = dx * dx + dy * dy + dz * dz
                   Rinv = 1 / sqrt(R2)
-                  if(R2.eq.0) Rinv = 0
+                  if (R2.eq.0) Rinv = 0
                   R3inv = q2(j) * Rinv * Rinv * Rinv
                   pp = pp + q2(j) * Rinv
                   fx = fx + dx * R3inv