Commits

Rio Yokota committed 61b34eb

Created Box structure for passing around bounding box.

Comments (0)

Files changed (6)

examples/parallel.cxx

     FMM.startTimer("Total FMM");
 
     FMM.partition(bodies);
-    FMM.setBounds(bodies);
-    FMM.buildTree(bodies,cells);
+    Box box = FMM.setBounds(bodies);
+    FMM.buildTree(bodies, cells, box);
     FMM.upwardPass(cells);
     FMM.startPAPI();
 
     FMM.setLET(cells);
     FMM.commBodies();
     FMM.commCells();
-    FMM.dualTreeTraversal(cells, cells);
+    FMM.dualTreeTraversal(cells, cells, FMM.periodicCycle, ARGS.mutual);
     jbodies = bodies;
     for( int irank=1; irank<FMM.MPISIZE; irank++ ) {
       FMM.getLET(jcells,(FMM.MPIRANK+irank)%FMM.MPISIZE);
       }
       assert( ic == int(icells.size()) );
 #endif
-      FMM.dualTreeTraversal(cells, jcells);
+      FMM.dualTreeTraversal(cells, jcells, FMM.periodicCycle, ARGS.mutual);
     }
 #else
     jbodies = bodies;
       FMM.setBounds(jbodies);
       FMM.buildTree(jbodies, jcells);
       FMM.upwardPass(jcells);
-      FMM.dualTreeTraversal(cells, jcells);
+      FMM.dualTreeTraversal(cells, jcells, FMM.periodicCycle, ARGS.mutual);
     }
 #endif
 

examples/serial.cxx

       << "Num bodies           : " << numBodies << std::endl;
     bodies.resize(numBodies);
     DATA.initBodies(bodies, ARGS.distribution);
-    FMM.setBounds(bodies);
-    FMM.buildTree(bodies, cells);                               // TODO : make it work without this
+    Box box = FMM.setBounds(bodies);
+    FMM.buildTree(bodies, cells, box);                          //TODO : make it work without this
     FMM.resetTimer();
     FMM.startTimer("Total FMM");
-    FMM.buildTree(bodies, cells);
+    FMM.buildTree(bodies, cells, box);
     FMM.upwardPass(cells);
     FMM.startPAPI();
-    FMM.dualTreeTraversal(cells, cells, ARGS.mutual);
+    FMM.dualTreeTraversal(cells, cells, FMM.periodicCycle, ARGS.mutual);
     FMM.stopPAPI();
     FMM.downwardPass(cells);
     std::cout << "----------------------------------" << std::endl;

include/serialfmm.h

  private:
   typedef std::pair<vec3,vec3> vec3Pair;                        //!< Pair of vec3
 
+ public:
+  real_t periodicCycle;                                         //!< Length of the global domain
+  real_t localRadius;                                           //!< Radius of local root cell
+  vec3   localCenter;                                           //!< Center of local root cell
+  fvec3  localXmin;                                             //!< Local Xmin for a given rank
+  fvec3  localXmax;                                             //!< Local Xmax for a given rank
+
  private:
 //! Error optimization of Rcrit
   void setRcrit(C_iter C, C_iter C0, real_t c) {
 
  public:
 //! Set center and size of root cell
-  void setBounds(Bodies &bodies) {
+  Box setBounds(Bodies &bodies) {
     startTimer("Set bounds");                                   // Start timer
     vec3Pair bounds = getBounds(bodies.begin(), bodies.end());  // Get Xmin (first) and Xmax (second) of domain
     for (int d=0; d<3; d++) localXmin[d] = bounds.first[d];     // Set local Xmin
     } else {                                                    // If periodic boundary condition
       periodicCycle = 2 * M_PI;                                 //  Set global radius
     }                                                           // End if for periodic boundary condition
+    Box box;
+    box.X = localCenter;
+    box.R = localRadius;
     stopTimer("Set bounds",printNow);
+    return box;
   }
 
 //! Upward pass (P2M, M2M)

include/traversal.h

   int IMAGES;                                                   //!< Number of periodic image sublevels
   float THETA;                                                  //!< Multipole acceptance criteria 
 
-  real_t localRadius;                                           //!< Radius of local root cell
-  vec3   localCenter;                                           //!< Center of local root cell
-  fvec3  localXmin;                                             //!< Local Xmin for a given rank
-  fvec3  localXmax;                                             //!< Local Xmax for a given rank
-
-  real_t periodicCycle;                                         //!< Diameter of global root cell
-
  private:
 //! Calculate Bmax
   real_t getBmax(vec3 const &X, C_iter C) const {
 
 
 //! Tree traversal of periodic cells
-  void traversePeriodic(real_t Length) {
+  void traversePeriodic(real_t length) {
     startTimer("Traverse periodic");                            // Start timer
     Xperiodic = 0;                                              // Periodic coordinate offset
     Cells pcells(27);                                           // Create cells
               for (int cx=-1; cx<=1; cx++) {                    //      Loop over x periodic direction (child)
                 for (int cy=-1; cy<=1; cy++) {                  //       Loop over y periodic direction (child)
                   for (int cz=-1; cz<=1; cz++) {                //        Loop over z periodic direction (child)
-                    Xperiodic[0] = (ix * 3 + cx) * Length;      //         Coordinate offset for x periodic direction
-                    Xperiodic[1] = (iy * 3 + cy) * Length;      //         Coordinate offset for y periodic direction
-                    Xperiodic[2] = (iz * 3 + cz) * Length;      //         Coordinate offset for z periodic direction
+                    Xperiodic[0] = (ix * 3 + cx) * length;      //         Coordinate offset for x periodic direction
+                    Xperiodic[1] = (iy * 3 + cy) * length;      //         Coordinate offset for y periodic direction
+                    Xperiodic[2] = (iz * 3 + cz) * length;      //         Coordinate offset for z periodic direction
                     M2L(Ci0,Ci,false);                          //         Perform M2L kernel
                   }                                             //        End loop over z periodic direction (child)
                 }                                               //       End loop over y periodic direction (child)
         for (int iy=-1; iy<=1; iy++) {                          //   Loop over y periodic direction
           for (int iz=-1; iz<=1; iz++, Cj++) {                  //    Loop over z periodic direction
             if( ix != 0 || iy != 0 || iz != 0 ) {               //     If periodic cell is not at center
-              Cj->X[0] = Ci->X[0] + ix * Length;                //      Set new x coordinate for periodic image
-              Cj->X[1] = Ci->X[0] + iy * Length;                //      Set new y cooridnate for periodic image
-              Cj->X[2] = Ci->X[0] + iz * Length;                //      Set new z coordinate for periodic image
+              Cj->X[0] = Ci->X[0] + ix * length;                //      Set new x coordinate for periodic image
+              Cj->X[1] = Ci->X[0] + iy * length;                //      Set new y cooridnate for periodic image
+              Cj->X[2] = Ci->X[0] + iz * length;                //      Set new z coordinate for periodic image
               Cj->M    = Ci->M;                                 //      Copy multipoles to new periodic image
             }                                                   //     Endif for periodic center cell
           }                                                     //    End loop over z periodic direction
       Ci->M = 0;                                                //  Reset multipoles of periodic parent
       setCenter(Ci,Cj0);                                        //  Set center of mass for periodic parent
       M2M(Ci,Cj0);                                              //  Evaluate periodic M2M kernels for this sublevel
-      Length *= 3;                                              //  Increase center cell size three times
+      length *= 3;                                              //  Increase center cell size three times
       Cj0 = C0;                                                 //  Reset Cj0 back
     }                                                           // End loop over sublevels of tree
     stopTimer("Traverse periodic",printNow);                    // Stop timer
   ~Traversal() {}
 
 //! Evaluate P2P and M2L using dual tree traversal
-  void dualTreeTraversal(Cells &icells, Cells &jcells, bool mutual=false) {
+  void dualTreeTraversal(Cells &icells, Cells &jcells, real_t length, bool mutual=false) {
     Ci0 = icells.begin();                                       // Set iterator of target root cell
     Cj0 = jcells.begin();                                       // Set iterator of source root cell
     startTimer("Traverse");                                     // Start timer
       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 * periodicCycle;                  //     Coordinate shift for x periodic direction
-            Xperiodic[1] = iy * periodicCycle;                  //     Coordinate shift for y periodic direction
-            Xperiodic[2] = iz * periodicCycle;                  //     Coordinate shift for z periodic direction
+            Xperiodic[0] = ix * length;                         //     Coordinate shift for x periodic direction
+            Xperiodic[1] = iy * length;                         //     Coordinate shift for y periodic direction
+            Xperiodic[2] = iz * length;                         //     Coordinate shift for z periodic direction
             traverse(Ci0,Cj0,false);                            //     Traverse the tree for this periodic image
           }                                                     //    End loop over z periodic direction
         }                                                       //   End loop over y periodic direction
       }                                                         //  End loop over x periodic direction
-      traversePeriodic(periodicCycle);                          //  Traverse tree for periodic images
+      traversePeriodic(length);                                 //  Traverse tree for periodic images
     }                                                           // End if for periodic boundary condition
     stopTimer("Traverse",printNow);                             // Stop timer
   }

include/treebuilder.h

 
 //! Build nodes of octree adaptively using a top-down approach based on recursion (uses task based thread parallelism)
   OctreeNode * buildNodes(Bodies& bodies, Bodies& buffer, int begin, int end,
-                          BinaryTreeNode * binNode, vec3 X, int level=0, bool direction=false) {
+                          BinaryTreeNode * binNode, vec3 X, real_t R0, int level=0, bool direction=false) {
     assert(getMaxBinNode(end - begin) <= binNode->END - binNode->BEGIN);
     if (begin == end) return NULL;                              // If no bodies left, return null pointer
     if (end - begin <= NCRIT) {                                 // If number of bodies is less than threshold
 	assert(binNodeOffset + maxBinNode <= binNode->END);
 	spawn_task2(buffer, bodies, {                           //  Spawn new task using Intel TBB or MTHREAD
 	    vec3 Xchild = X;                                    //   Initialize center position of child node
-	    real_t r = localRadius / (1 << (level + 1));        //   Radius of cells for child's level
+	    real_t r = R0 / (1 << (level + 1));                 //   Radius of cells for child's level
 	    for (int d=0; d<3; d++) {                           //   Loop over dimensions
 	      Xchild[d] += r * (((i & 1 << d) >> d) * 2 - 1);   //    Shift center position to that of child node
 	    }                                                   //   End loop over dimensions
 	    binNodeChild->BEGIN = binNodeOffset;                //   Assign first memory address from offset
 	    binNodeChild->END = binNodeOffset + maxBinNode;     //   Keep track of last memory address
 	    octNode->CHILD[i] = buildNodes(buffer, bodies,      //   Recursive call for each child
-	       octantOffset[i], octantOffset[i] + binNode->NBODY[i],//   Range of bodies is calcuated from octant offset
-	       binNodeChild, Xchild, level+1, !direction);      //   Alternate copy direction bodies <-> buffer
+	      octantOffset[i], octantOffset[i] + binNode->NBODY[i],//   Range of bodies is calcuated from octant offset
+	      binNodeChild, Xchild, R0, level+1, !direction);   //   Alternate copy direction bodies <-> buffer
 	  });                                                   //  Close lambda expression
 	binNodeOffset += maxBinNode;                            //  Increment offset for binNode memory address
       }                                                         // End loop over children
   }
 
 //! Create cell data structure from nodes
-  void nodes2cells(OctreeNode * octNode, C_iter C, C_iter C0, C_iter CN, int level=0, int iparent=0) {
+  void nodes2cells(OctreeNode * octNode, C_iter C, C_iter C0, C_iter CN, real_t R0, int level=0, int iparent=0) {
     C->PARENT = iparent;                                        // Index of parent cell
-    C->R      = localRadius / (1 << level);                     // Cell radius
+    C->R      = R0 / (1 << level);                              // Cell radius
     C->X      = octNode->X;                                     // Cell center
     C->NDBODY = octNode->NBODY;                                 // Number of decendant bodies
     C->BODY   = B0 + octNode->BODY;                             // Iterator of first body in cell
 	for (int i=0; i<nchild; i++) {                          //  Loop over children
 	  int octant = octants[i];                              //   Get octant from child index
 	  spawn_task0_if(octNode->NNODE > 1000,                 //   Spawn task if number of sub-nodes is large
-			 nodes2cells(octNode->CHILD[octant], Ci, C0, CN, level+1, C-C0));// Recursive call for each child
+	    nodes2cells(octNode->CHILD[octant], Ci, C0, CN, R0, level+1, C-C0));// Recursive call for each child
 	  Ci++;                                                 //   Increment cell iterator
 	  CN += octNode->CHILD[octant]->NNODE - 1;              //   Increment next free memory address
 	}                                                       //  End loop over children
   }
 
 //! Grow tree structure top down
-  void growTree(Bodies &bodies) {
+  void growTree(Bodies &bodies, vec3 X0, real_t R0) {
     Bodies buffer = bodies;                                     // Copy bodies to buffer
     startTimer("Grow tree");                                    // Start timer
     B0 = bodies.begin();                                        // Bodies iterator
     int maxBinNode = getMaxBinNode(bodies.size());              // Get maximum size of binary tree
     binNode->BEGIN = new BinaryTreeNode[maxBinNode];            // Allocate array for binary tree nodes
     binNode->END = binNode->BEGIN + maxBinNode;                 // Set end pointer
-    N0 = buildNodes(bodies, buffer, 0, bodies.size(), binNode, localCenter);// Build tree recursively
+    N0 = buildNodes(bodies, buffer, 0, bodies.size(), binNode, X0, R0);// Build tree recursively
     delete[] binNode->BEGIN;                                    // Deallocate binary tree array
     stopTimer("Grow tree",printNow);                            // Stop timer
   }
 
 //! Link tree structure
-  void linkTree(Cells &cells) {
+  void linkTree(Cells &cells, real_t R0) {
     startTimer("Link tree");                                    // Start timer
     cells.resize(N0->NNODE);                                    // Allocate cells array
     C_iter C0 = cells.begin();                                  // Cell begin iterator
-    nodes2cells(N0, C0, C0, C0+1);                              // Convert nodes to cells recursively
+    nodes2cells(N0, C0, C0, C0+1, R0);                          // Convert nodes to cells recursively
     delete N0;                                                  // Deallocate nodes
     stopTimer("Link tree",printNow);                            // Stop timer
   }
   ~TreeBuilder() {}
 
 //! Build tree structure top down
-  void buildTree(Bodies &bodies, Cells &cells) {
-    growTree(bodies);                                           // Grow tree from root
-    linkTree(cells);                                            // Form parent-child links in tree
+  void buildTree(Bodies &bodies, Cells &cells, Box box) {
+    growTree(bodies,box.X,box.R);                               // Grow tree from root
+    linkTree(cells,box.R);                                      // Form parent-child links in tree
   }
 
 //! Print tree structure statistics
 typedef vec<NTERM,complex_t> vecL;                              //!< Local coefficient type for spherical
 #endif
 
+//! Structure for defining bounding box
+struct Box {
+  vec3   X;                                                     //!< Box center
+  real_t R;                                                     //!< Box radius
+};
+
 //! Structure of aligned source for SIMD
 struct Source {
   vec3   X;                                                     //!< Position