Commits

Rio Yokota committed 89254c1

Added Ewald but not working yet.

Comments (0)

Files changed (15)

 help:
-	@echo "Help documentation will be available soon...\n"
+	@echo "Please read the README file in the root directory."
 clean:
 	find . -name "*.o" -o -name "*.out*" | xargs rm -f
 cleandat:
 OBJECTS	= $(SOURCES:.cxx=.o)
 endif
 
+MAKEFLAGS += --no-print-directory
+
 .cxx.o  :
 	$(CXX) -c $^ -o $@ $(LFLAGS)
 .cu.o   :
-1. Quick start guide
-There is no need to configure or CMake. Just do the following.
+--- 1. Quick start guide
+There is no need to configure or CMake.
+You may need to change the CXX command in Makefile.include.
+Then do the following.
 
 To run the serial FMM example:
 cd examples
 To get the list of command line options:
 ./a.out --help
 
-2. Where are the important files?
+
+--- 2. Where are the important files?
 Most of the source code is in the directory "include" in the form of header files.
 Different versions of P2M,M2M,M2L,L2L,L2P,P2P kernel implementations are in "kernels".
 File names and variable names are self-explanatory, and every line is commented.
 The best documentation is most likely the code itself.
 
 
-3. What's in the directory "gpu"?
+--- 3. What's in the directory "gpu"?
 "gpu" is a completely independent code at the moment.
 It will be integrated into the main code soon.
 To run it:
 make
 
 
-4. What's in Makefile.include?
+--- 4. What's in Makefile.include?
 This file is included from each Makefile in all subdirectories.
 Most options of exafmm can be controlled from this file.
 
 LFLAGS : Linker options and preprocessor options (see comments)
 
 
-5. Why is there a CMakeLists.txt?
+--- 5. Why is there a CMakeLists.txt?
 This project is not yet large enough to justify the use of CMake,
 but we are just testing a few features of CMake.
 You probably don't want to use this feature.

examples/Makefile

 
 #LFLAGS	+= -DIneJ # Different target and source points
 
+# Test for serial FMM
 serial: serial.o $(OBJECTS)
 	$(CXX) $? $(LFLAGS)
 	./a.out
 
+# Test for parallel FMM
 parallel: parallel.o $(OBJECTS)
 	$(CXX) $? $(LFLAGS)
 	mpirun -np 8 ./a.out
 
+# Checking O(N) complexity
 complexity: serial.o $(OBJECTS)
 	$(CXX) $? $(LFLAGS)
 	for N in 10000 20000 50000 100000 200000 500000 1000000; do \
         ./a.out --numBodies $$N --verbose 0; echo; \
 	done
 
+# Test for kernels only
 kernel: kernel.o $(OBJECTS)
 	$(CXX) $? $(LFLAGS)
 	./a.out
 
+# Checking O(theta^p) convergence (Turn on VTK to see a graph at the end)
 convergence: kernel.cxx $(SOURCES)
 	rm -f kernel.dat
 	for P in 3 4 5 6 7 8; do echo P = $$P compiling && \
 	$(CXX) $? $(LFLAGS) -DEXPANSION=$$P && \
-	echo P = $$P executing && ./a.out; done
+	echo P = $$P executing && ./a.out; done
+
+# Ewald vs. periodic FMM
+ewald: ewald.o $(OBJECTS)
+	$(CXX) $? $(LFLAGS)
+	./a.out --numBodies 10000

examples/ewald.cxx

+#include "args.h"
+#include "boundbox.h"
+#include "buildtree.h"
+#include "dataset.h"
+#include "ewald.h"
+#include "logger.h"
+#include "traversal.h"
+#include "updownpass.h"
+#if VTK
+#include "vtk.h"
+#endif
+
+int main(int argc, char ** argv) {
+  Args args(argc, argv);
+  Bodies bodies;
+  Cells cells;
+  Dataset data;
+  Logger logger;
+
+  const real_t ksize = 11.;
+  const real_t alpha = .2;
+  const real_t sigma = .25 / M_PI;
+  const real_t theta = .5;
+  const real_t cycle = 100.;
+  BoundBox boundbox(args.NSPAWN);
+  BuildTree tree(args.NCRIT,args.NSPAWN);
+  UpDownPass pass(args.THETA);
+  Traversal traversal(args.NSPAWN,args.IMAGES);
+  Ewald ewald(ksize,alpha,sigma,theta,cycle);
+  logger.verbose = true;
+  if (args.verbose) {
+    boundbox.verbose = true;
+    tree.verbose = true;
+    pass.verbose = true;
+    traversal.verbose = true;
+    logger.printTitle("Parameters");
+  }
+  args.print(logger.stringLength,P);
+#if AUTO
+  traversal.timeKernels();
+#endif
+#if _OPENMP
+#pragma omp parallel
+#pragma omp master
+#endif
+  if(args.verbose) logger.printTitle("Profiling");
+  logger.startTimer("Total FMM");
+  logger.startPAPI();
+  bodies.resize(args.numBodies);
+  data.initBodies(bodies, args.distribution);
+  Bounds bounds = ewald.rescale(bodies);
+  Box box = boundbox.bounds2box(bounds);
+  tree.buildTree(bodies, cells, box);
+  pass.upwardPass(cells);
+  traversal.dualTreeTraversal(cells, cells, cycle, args.mutual);
+  pass.downwardPass(cells);
+  if (args.verbose) logger.printTitle("Total runtime");
+  logger.stopPAPI();
+  logger.stopTimer("Total FMM", logger.verbose);
+  boundbox.writeTime();
+  tree.writeTime();
+  pass.writeTime();
+  traversal.writeTime();
+  boundbox.resetTimer();
+  tree.resetTimer();
+  pass.resetTimer();
+  traversal.resetTimer();
+  logger.resetTimer();
+#if 1
+  Bodies bodies2 = bodies;
+  data.initTarget(bodies2);
+  logger.startTimer("Total Ewald");
+  ewald.realPart(cells,cells);
+  ewald.wavePart(bodies2);
+  logger.stopTimer("Total Ewald", args.verbose);
+#else
+  Bodies jbodies = bodies;
+  if (int(bodies.size()) > args.numTarget) data.sampleBodies(bodies, args.numTarget);
+  Bodies bodies2 = bodies;
+  data.initTarget(bodies2);
+  logger.startTimer("Total Direct");
+  traversal.direct(bodies2, jbodies, cycle);
+  traversal.normalize(bodies2);
+  logger.stopTimer("Total Direct", args.verbose);
+#endif
+  double diff1 = 0, norm1 = 0, diff2 = 0, norm2 = 0;
+  data.evalError(bodies, bodies2, diff1, norm1, diff2, norm2);
+  if (args.verbose) {
+    logger.printTitle("FMM vs. Ewald");
+    logger.printError(diff1, norm1, diff2, norm2);
+    tree.printTreeData(cells);
+    traversal.printTraversalData();
+    logger.printPAPI();
+  }
+#if VTK
+  for (B_iter B=bodies.begin(); B!=bodies.end(); B++) B->ICELL = 0;
+  for (C_iter C=cells.begin(); C!=cells.end(); C++) {
+    Body body;
+    body.ICELL = 1;
+    body.X     = C->X;
+    body.SRC   = 0;
+    bodies.push_back(body);
+  }
+  vtk3DPlot vtk;
+  vtk.setBounds(M_PI,0);
+  vtk.setGroupOfPoints(bodies);
+  vtk.plot();
+#endif
+  return 0;
+}

examples/kernel.cxx

 #include <vector>
 #if VTK
 #include "vtk.h"
-#else
-#error Turn on VTK for this test
 #endif
 #if COMkernel
 #error Turn off COMkernel for this test

examples/parallel.cxx

   MPI_Reduce(&diff2, &diff4, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
   MPI_Reduce(&norm2, &norm4, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
   if(args.verbose) {
+    logger.printTitle("FMM vs. direct");
     logger.printError(diff3, norm3, diff4, norm4);
     tree.printTreeData(cells);
     traversal.printTraversalData();

examples/serial.cxx

   double diff1 = 0, norm1 = 0, diff2 = 0, norm2 = 0;
   data.evalError(bodies, bodies2, diff1, norm1, diff2, norm2);
   if (args.verbose) {
+    logger.printTitle("FMM vs. direct");
     logger.printError(diff1, norm1, diff2, norm2);
     tree.printTreeData(cells);
     traversal.printTraversalData();

include/boundbox.h

 
  public:
   BoundBox(int nspawn) : NSPAWN(nspawn) {}
-  ~BoundBox() {}
 
   // ! Get Xmin and Xmax of domain
   Bounds getBounds(Bodies bodies) {

include/buildtree.h

 
  public:
   BuildTree(int ncrit, int nspawn) : NCRIT(ncrit), NSPAWN(nspawn), MAXLEVEL(0) {}
-  ~BuildTree() {}
 
 //! Build tree structure top down
   void buildTree(Bodies &bodies, Cells &cells, Box box) {

include/dataset.h

  public:
 //! Constructor
   Dataset() : filePosition(0) {}
-//! Destructor
-  ~Dataset() {}
 
 //! Initialize target values
   void initTarget(Bodies &bodies) {
+#ifndef ewald_h
+#define ewald_h
+#include <stack>
+#include "types.h"
+
+class Ewald {
+//! Wave structure for Ewald summation
+  struct Wave {
+    vec3   K;                                                   //!< 3-D wave number vector
+    real_t REAL;                                                //!< real part of wave
+    real_t IMAG;                                                //!< imaginary part of wave
+  };
+  typedef std::vector<Wave>           Waves;                    //!< Vector of Wave types
+  typedef std::vector<Wave>::iterator W_iter;                   //!< Iterator for Wave types
+
+ private:
+  real_t KSIZE;                                                 //!< Number of waves in Ewald summation
+  real_t ALPHA;                                                 //!< Scaling parameter for Ewald summation
+  real_t SIGMA;                                                 //!< Scaling parameter for Ewald summation
+  real_t THETA;                                                 //!< Neighbor acceptance criteria
+  real_t CYCLE;                                                 //!< Periodic cycle
+  vec3   Xperiodic;                                             //!< Coordinate offset for periodic B.C.
+
+ private:
+//! Forward DFT
+  void dft(Waves &waves, Bodies &bodies) {
+    real_t scale = 2 * M_PI / CYCLE;                            // Scale conversion
+    for (W_iter W=waves.begin(); W!=waves.end(); W++) {         // Loop over waves
+      W->REAL = W->IMAG = 0;                                    //  Initialize waves
+      for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {     //  Loop over bodies
+	real_t th = 0;                                          //   Initialize phase
+	for (int d=0; d<3; d++) th += W->K[d] * B->X[d] * scale;//   Determine phase
+	W->REAL += B->SRC * cos(th);                            //   Accumulate real component
+	W->IMAG += B->SRC * sin(th);                            //   Accumulate imaginary component
+      }                                                         //  End loop over bodies
+    }                                                           // End loop over waves
+  }
+
+//! Inverse DFT
+  void idft(Waves &waves, Bodies &bodies) {
+    real_t scale = 2 * M_PI / CYCLE;                            // Scale conversion
+    for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {       // Loop over bodies
+      kvec4 TRG = 0;                                            //  Initialize target values
+      for (W_iter W=waves.begin(); W!=waves.end(); W++) {       //   Loop over waves
+	real_t th = 0;                                          //    Initialzie phase
+	for (int d=0; d<3; d++) th += W->K[d] * B->X[d] * scale;//    Determine phase
+	real_t dtmp = W->REAL * sin(th) - W->IMAG * cos(th);    //    Temporary value
+	TRG[0]     += W->REAL * cos(th) + W->IMAG * sin(th);    //    Accumulate potential
+	for (int d=0; d<3; d++) TRG[d+1] -= dtmp * W->K[d];     //    Accumulate force
+      }                                                         //   End loop over waves
+      B->TRG = TRG;                                             //  Copy results to bodies
+    }                                                           // End loop over bodies
+  }
+
+//! Initialize wave vector
+  Waves initWaves() {
+    Waves waves;                                                // Initialzie wave vector
+    real_t kmaxsq = KSIZE * KSIZE;                              // kmax squared
+    int kmax = int(KSIZE);                                      // kmax as integer
+    for (int l=0; l<=kmax; l++) {                               // Loop over x component
+      int mmin = -kmax;                                         //  Determine minimum y component
+      if (l==0) mmin = 0;                                       //  Exception for minimum y component
+      for (int m=mmin; m<=kmax; m++) {                          //  Loop over y component
+	int nmin = -kmax;                                       //   Determine minimum z component
+	if (l==0 && m==0) nmin=1;                               //   Exception for minimum z component
+	for (int n=nmin; n<=kmax; n++) {                        //   Loop over z component
+	  real_t ksq = l * l + m * m + n * n;                   //    Wave number squared
+	  if (ksq <= kmaxsq) {                                  //    If wave number is below kmax
+	    Wave wave;                                          //     Initialzie wave structure
+	    wave.K[0] = l;                                      //     x component of k
+	    wave.K[1] = m;                                      //     y component of k
+	    wave.K[2] = n;                                      //     z component of k
+	    wave.REAL = wave.IMAG = 0;                          //     Initialize amplitude
+	    waves.push_back(wave);                              //     Push wave to vector
+	  }                                                     //    End if for wave number
+	}                                                       //   End loop over z component
+      }                                                         //  End loop over y component
+    }                                                           // End loop over x component
+    return waves;                                               // Return wave vector
+  }
+//! Ewald real part P2P kernel
+  void P2P(C_iter Ci, C_iter Cj) {
+    for (B_iter Bi=Ci->BODY; Bi!=Ci->BODY+Ci->NDBODY; Bi++) {   // Loop over target bodies
+      for (B_iter Bj=Cj->BODY; Bj!=Cj->BODY+Cj->NDBODY; Bj++) { //  Loop over source bodies
+	vec3 dist = Bi->X - Bj->X - Xperiodic;                  //   Distance vector from source to target
+	for (int d=0; d<3; d++) {                               //   Loop over dimensions
+	  if (dist[d] < -CYCLE/2) dist[d] += CYCLE;             //    Wrap domain so that target is always at
+	  if (dist[d] >= CYCLE/2) dist[d] -= CYCLE;             //     the center of a [-CYCLE/2,CYCLE/2]^3 source cube
+	}                                                       //   End loop over dimensions
+	real_t R2 = norm(dist);                                 //   R^2
+	if (R2 != 0) {                                          //   Exclude self interaction
+	  real_t R2s = R2 * ALPHA * ALPHA;                      //    (R * alpha)^2
+	  real_t Rs = std::sqrt(R2s);                           //    R * alpha
+	  real_t invRs = 1 / Rs;                                //    1 / (R * alpha)
+	  real_t invR2s = invRs * invRs;                        //    1 / (R * alpha)^2
+	  real_t invR3s = invR2s * invRs;                       //    1 / (R * alpha)^3
+	  real_t dtmp = Bj->SRC * (M_2_SQRTPI * exp(-R2s) * invR2s + erfc(Rs) * invR3s);
+	  dtmp *= ALPHA * ALPHA * ALPHA;                        //    Scale temporary value
+	  Bi->TRG[0] += Bj->SRC * erfc(Rs) * invRs * ALPHA;     //    Ewald real potential
+	  Bi->TRG[1] -= dist[0] * dtmp;                         //    x component of Ewald real force
+	  Bi->TRG[2] -= dist[1] * dtmp;                         //    y component of Ewald real force
+	  Bi->TRG[3] -= dist[2] * dtmp;                         //    z component of Ewald real force
+	}                                                       //   End if for self interaction
+      }                                                         //  End loop over source bodies
+    }                                                           // End loop over target bodies
+  }
+
+//! Traverse tree to find neighbors
+  void traverse(C_iter Ci, C_iter C, C_iter C0) {
+    std::stack<C_iter> cellStack;                               // Stack of cell iterators
+    cellStack.push(C);                                          // Push pair to queue
+    while (!cellStack.empty()) {                                // While traversal stack is not empty
+      C = cellStack.top();                                      //  Get cell from top of stack
+      cellStack.pop();                                          //  Pop traversal stack
+      for (C_iter Cj=C0+C->CHILD; Cj!=C0+C->CHILD+C->NCHILD; Cj++) {// Loop over cell's children
+        vec3 dX = Ci->X - Cj->X - Xperiodic;                    //   Distance vector from source to target
+        real_t Rq = std::sqrt(norm(dX));                        //   Scalar distance
+        if (Rq * THETA < Ci->R + Cj->R && Cj->NCHILD == 0) {    //   If twigs are close
+          P2P(Ci,Cj);                                           //    Ewald real part
+        } else if (Cj->NCHILD != 0) {                           //   If cells are not twigs
+          cellStack.push(Cj);                                   //    Push source cell to stack
+        }                                                       //   End if for twig cells
+      }                                                         //  End loop over cell's children
+    }                                                           // End while loop for traversal stack
+  }
+
+ public:
+//! Constructor
+ Ewald(real_t ksize, real_t alpha, real_t sigma, real_t theta, real_t cycle) :
+  KSIZE(ksize), ALPHA(alpha), SIGMA(sigma), THETA(theta), CYCLE(cycle), Xperiodic(0) {}
+
+//! Rescale positions and charges
+  Bounds rescale(Bodies &bodies) {
+    const int numBodies = bodies.size();                        // Number of bodies
+    real_t average = 0;                                         // Initialize average charge
+    for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {       // Loop over bodies
+      B->X *= CYCLE / (2 * M_PI);                               //  Rescale positions
+      B->SRC = drand48() / numBodies;                           //  Randomize charges
+      average += B->SRC;                                        //  Accumulate charges
+    }                                                           // 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
+    Bounds bounds;                                              // Initialize bounds
+    bounds.Xmin = -CYCLE/2;                                     // Force Xmin to be -CYCLE/2
+    bounds.Xmax = CYCLE/2;                                      // Force Xmax to be CYCLE/2
+    return bounds;                                              // Return bounds
+  }
+
+//! Ewald real part
+  void realPart(Cells &cells, Cells &jcells) {
+    C_iter Cj = jcells.begin();                                 // Set begin iterator for source cells
+    for (C_iter Ci=cells.begin(); Ci!=cells.end(); Ci++) {      // Loop over target cells
+      if (Ci->NCHILD == 0) {                                    //  If cell is a twig
+        for (int ix=-1; ix<=1; ix++) {                          //   Loop over x periodic direction
+          for (int iy=-1; iy<=1; iy++) {                        //    Loop over y periodic direction
+            for (int iz=-1; iz<=1; iz++) {                      //     Loop over z periodic direction
+              Xperiodic[0] = ix * CYCLE;                        //      Coordinate offset for x periodic direction
+              Xperiodic[1] = iy * CYCLE;                        //      Coordinate offset for y periodic direction
+              Xperiodic[2] = iz * CYCLE;                        //      Coordinate offset for z periodic direction
+              traverse(Ci,Cj,Cj);                               //      Traverse the source tree
+            }                                                   //     End loop over z periodic direction
+          }                                                     //    End loop over y periodic direction
+        }                                                       //   End loop over x periodic direction
+        for (B_iter B=Ci->BODY; B!=Ci->BODY+Ci->NCBODY; B++) {  //   Loop over all bodies in cell
+          B->TRG[0] -= M_2_SQRTPI * B->SRC * ALPHA;             //    Self term of Ewald real part
+        }                                                       //   End loop over all bodies in cell
+      }                                                         //  End if for twig cells
+    }                                                           // End loop over target cells
+  }
+
+//! Ewald wave part
+  void wavePart(Bodies &bodies) {
+    Waves waves = initWaves();                                  // Initialize wave vector
+    dft(waves,bodies);                                          // Apply DFT to bodies to get waves
+    real_t scale = 2 * M_PI / CYCLE;                            // Scale conversion
+    real_t coef = .5 / M_PI / M_PI / SIGMA / CYCLE;             // First constant
+    real_t coef2 = scale * scale / (4 * ALPHA * ALPHA);         // Second constant
+    for (W_iter W=waves.begin(); W!=waves.end(); W++) {         // Loop over waves
+      real_t K2 = norm(W->K);                                   //  Wave number squared
+      real_t factor = coef * exp(-K2 * coef2) / K2;             //  Wave factor
+      W->REAL *= factor;                                        //  Apply wave factor to real part
+      W->IMAG *= factor;                                        //  Apply wave factor to imaginary part
+    }                                                           // End loop over waves
+    idft(waves,bodies);                                         // Inverse DFT
+    for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {       // Loop over bodies
+      for (int d=0; d<3; d++) B->TRG[d+1] *= scale;             //  Scale forces
+    }                                                           // End loop over bodies
+
+    vec3 dipole = 0;                                            // Initialize dipole correction
+    for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {       // Loop over bodies
+      dipole += (B->X - CYCLE/2) * B->SRC;                      //  Calcuate dipole of the whole system
+    }                                                           // End loop over bodies
+    coef = 4 * M_PI / (3 * CYCLE * CYCLE * CYCLE);              // Precalcualte constant
+    for (B_iter B=bodies.begin(); B!=bodies.end(); B++) {       // Loop over bodies
+      B->TRG[0] += coef * norm(dipole) / bodies.size() / B->SRC;//  Dipole correction for potential
+      for (int d=0; d!=3; d++) {                                //  Loop over dimensions
+	B->TRG[d+1] += coef * dipole[d];                        //   Dipole correction for forces
+      }                                                         //  End loop over dimensions
+    }                                                           // End loop over bodies
+  }
+};
+#endif
 
   //! Print relative L2 norm error
   void printError(double diff1, double norm1, double diff2, double norm2) {
-    printTitle("FMM vs. direct");                               // Print title
     std::cout << std::setw(stringLength) << std::left           // Set format
               << "Rel. L2 Error (pot)" << " : " << std::sqrt(diff1/norm1) << std::endl;// Print potential error
     if( std::abs(diff2) > 0 ) {                                 // If acceleration was calculated

include/traversal.h

     stopTimer("Traverse periodic",verbose);                     // Stop timer
   }
 
- protected:
  public:
   Traversal(int nspawn, int images) : NSPAWN(nspawn), IMAGES(images), NP2P(0), NM2L(0) {}
-  ~Traversal() {}
 
 //! Evaluate P2P and M2L using dual tree traversal
   void dualTreeTraversal(Cells &icells, Cells &jcells, real_t cycle, bool mutual=false) {

include/updownpass.h

   }
 
 //! Recursive call for upward pass
-  void upwardRecursion(C_iter C, C_iter C0) {
+  void postOrderTraversal(C_iter C, C_iter C0) {
     spawn_tasks {                                               // Initialize tasks
       for (C_iter CC=C0+C->CHILD; CC!=C0+C->CHILD+C->NCHILD; CC++) {// Loop over child cells
-	spawn_task0(upwardRecursion(CC, C0));                   //  Recursive call with new task
+	spawn_task0(postOrderTraversal(CC, C0));                //  Recursive call with new task
       }                                                         // End loop over child cells
       sync_tasks;                                               // Synchronize tasks
     }
   }
 
 //! Recursive call for downward pass
-  void downwardRecursion(C_iter C, C_iter C0) const {
+  void preOrderTraversal(C_iter C, C_iter C0) const {
     L2L(C,C0);                                                  // L2L kernel
     L2P(C);                                                     // L2P kernel
     spawn_tasks {                                               // Initialize tasks
       for (C_iter CC=C0+C->CHILD; CC!=C0+C->CHILD+C->NCHILD; CC++) {// Loop over child cells
-	spawn_task0(downwardRecursion(CC, C0));                 //  Recursive call with new task
+	spawn_task0(preOrderTraversal(CC, C0));                 //  Recursive call with new task
       }                                                         // End loop over chlid cells
       sync_tasks;                                               // Synchronize tasks
     }
 
  public:
   UpDownPass(real_t theta) : THETA(theta) {}
-  ~UpDownPass() {}
 
 //! Upward pass (P2M, M2M)
   void upwardPass(Cells &cells) {
     startTimer("Upward pass");                                  // Start timer
     C_iter C0 = cells.begin();                                  // Set iterator of target root cell
-    upwardRecursion(C0, C0);                                    // Recursive call for upward pass
+    postOrderTraversal(C0, C0);                                 // Recursive call for upward pass
     real_t c = (1 - THETA) * (1 - THETA) / std::pow(THETA,P+2) / std::pow(std::abs(C0->M[0]),1.0/3); // Root coefficient
     setRcrit(C0, C0, c);                                        // Error optimization of Rcrit
     if( cells.size() > 9 ) {                                    // If tree has more than 2 levels
     L2P(C0);                                                    // If root is the only cell do L2P
     spawn_tasks {                                               // Initialize tasks
       for (C_iter CC=C0+C0->CHILD; CC!=C0+C0->CHILD+C0->NCHILD; CC++) {// Loop over child cells
-	spawn_task0(downwardRecursion(CC, C0));                 //  Recursive call for downward pass
+	spawn_task0(preOrderTraversal(CC, C0));                 //  Recursive call for downward pass
       }                                                         // End loop over child cells
       sync_tasks;                                               // Synchronize tasks
     }