Source

exafmm / unit_test / ewald_fmm.cxx

Full commit
#include "serialfmm.h"

int main() {
  const int numBodies = 3000;                                   // Number of bodies - DEF=1000.
  const real delBodies = pow(numBodies, 1.0/3.0)/10.0;         // change in numBodies
  const real xmax = 3.10736;                                      // Size of domain - DEF=100.0
  const real delXmax = 100.0/xmax;                             // Change in xmax
  const real ksize = 44.0*delBodies;                           // Ewald wave number
  const real alpha = 1.0/.289108; //0.2*delXmax*delBodies;                    // Ewald alpha value
  const real sigma = .25 / M_PI;                               // Ewald sigma value
  IMAGES = 8;                                                  // Level of periodic image tree (0 for non-periodic)
  THETA = 1 / sqrt(4);                                         // Multipole acceptance criteria
  
  std::cout << "numBodies: " << numBodies << ", delBodies: " << delBodies << std::endl << "xmax: " << xmax << ", delXmax: " << delXmax << std::endl << "ksize: " << ksize << ", alpha: " << alpha << std::endl << "sigma: " << sigma <<", IMAGES: "  << IMAGES << ", THETA: " << THETA << std::endl;

  Bodies bodies(numBodies);                                     // Define vector of bodies
  Bodies jbodies;                                               // Define vector of source bodies
  Cells cells, jcells;                                          // Define vector of cells
  SerialFMM<Laplace> FMM;                                       // Instantiate SerialFMM class
  FMM.initialize();                                             // Initialize FMM
  FMM.printNow = true;                                          // Print timings

  FMM.startTimer("Set bodies");                                 // Start timer
  srand48(0);                                                   // Seed for random number generator
  real average = 0;                                             // Initialize average charge
  for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) {        // Loop over bodies
    for( int d=0; d!=3; ++d ) {                                 //  Loop over dimensions
      B->X[d] = drand48() * xmax;                               //   Initialize positions
    }                                                           //  End loop over dimensions
    B->SRC = drand48() / numBodies;                             //  Set charges
    average += B->SRC;                                          //  Accumulate charges
    B->TRG = 0;                                                 //  Initialize target values
    B->IBODY = B-bodies.begin();
  }                                                             // End loop over bodies
  average /= numBodies;                                         // Divide by total to get average
  for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) {        // Loop over bodies
    B->SRC -= average;                                          //  Make average charge 0
  }                                                             // End loop over bodies
  FMM.stopTimer("Set bodies",FMM.printNow);                     // Stop timer
  FMM.eraseTimer("Set bodies");                                 // Erase entry from timer to avoid timer overlap

  FMM.startTimer("Set domain");                                 // Start timer
  FMM.setDomain(bodies,xmax/2,xmax/2);                          // Set domain size of FMM
  FMM.stopTimer("Set domain",FMM.printNow);                     // Stop timer
  FMM.eraseTimer("Set domain");                                 // Erase entry from timer to avoid timer overlap

  FMM.bottomup(bodies,cells);                                   // Tree construction (bottom up) & upward sweep
  FMM.setEwald(ksize,alpha,sigma);                              // Set Ewald method paramters
  jcells = cells;                                               // Copy cells to jcells
  FMM.Ewald(bodies,cells,jcells);                               // Ewald summation

  Bodies bodies2 = bodies;                                      // Define new bodies vector for direct sum
  for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) B->TRG = 0;
  FMM.startTimer("Downward");                                   // Start timer
  FMM.downward(cells,jcells);                                   // Downward sweep
  FMM.stopTimer("Downward",FMM.printNow);                       // Stop timer
  FMM.eraseTimer("Downward");                                   // Erase entry from timer to avoid timer overlap

  real diff1 = 0, norm1 = 0, diff2 = 0, norm2 = 0;              // Initialize accumulators
  FMM.evalError(bodies,bodies2,diff1,norm1,diff2,norm2,true);   // Evaluate error
  FMM.printError(diff1,norm1,diff2,norm2);                      // Print the L2 norm error
  FMM.finalize();                                               // Finalize FMM
}