|
ExaFMM 1
Fast-multipole Method for exascale systems
|
include/bottomup.h
Go to the documentation of this file.
00001 /* 00002 Copyright (C) 2011 by Rio Yokota, Simon Layton, Lorena Barba 00003 00004 Permission is hereby granted, free of charge, to any person obtaining a copy 00005 of this software and associated documentation files (the "Software"), to deal 00006 in the Software without restriction, including without limitation the rights 00007 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 00008 copies of the Software, and to permit persons to whom the Software is 00009 furnished to do so, subject to the following conditions: 00010 00011 The above copyright notice and this permission notice shall be included in 00012 all copies or substantial portions of the Software. 00013 00014 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 00015 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 00016 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 00017 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 00018 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 00019 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 00020 THE SOFTWARE. 00021 */ 00022 #ifndef bottomup_h 00023 #define bottomup_h 00024 #include "topdown.h" 00025 00026 //! Bottomup tree constructor 00027 template<Equation equation> 00028 class BottomUp : public TopDown<equation> { 00029 public: 00030 using Kernel<equation>::printNow; //!< Switch to print timings 00031 using Kernel<equation>::startTimer; //!< Start timer for given event 00032 using Kernel<equation>::stopTimer; //!< Stop timer for given event 00033 using Kernel<equation>::sortBodies; //!< Sort bodies according to cell index 00034 using Kernel<equation>::X0; //!< Center of root cell 00035 using Kernel<equation>::R0; //!< Radius of root cell 00036 using TreeStructure<equation>::buffer; //!< Buffer for MPI communication & sorting 00037 using TreeStructure<equation>::getLevel; //!< Get level from cell index 00038 00039 protected: 00040 //! Max level for bottom up tree build 00041 int getMaxLevel(Bodies &bodies) { 00042 const long N = bodies.size() * MPISIZE; // Number of bodies 00043 int level; // Max level 00044 level = N >= NCRIT ? 1 + int(log(N / NCRIT)/M_LN2/3) : 0; // Decide max level from N/Ncrit 00045 int MPIlevel = int(log(MPISIZE-1) / M_LN2 / 3) + 1; // Level of local root cell 00046 if( MPISIZE == 1 ) MPIlevel = 0; // For serial execution local root cell is root cell 00047 if( MPIlevel > level ) { // If process hierarchy is deeper than tree 00048 // std::cout << "Process hierarchy is deeper than tree @ rank" << MPIRANK << std::endl; 00049 level = MPIlevel; 00050 } 00051 return level; // Return max level 00052 } 00053 00054 public: 00055 //! Constructor 00056 BottomUp() : TopDown<equation>() {} 00057 //! Destructor 00058 ~BottomUp() {} 00059 00060 //! Set cell index of all bodies 00061 void setIndex(Bodies &bodies, int level=-1, int begin=0, int end=0, bool update=false) { 00062 startTimer("Set index "); // Start timer 00063 bigint i; // Levelwise cell index 00064 if( level == -1 ) level = getMaxLevel(bodies); // Decide max level 00065 bigint off = ((1 << 3*level) - 1) / 7; // Offset for each level 00066 real r = R0 / (1 << (level-1)); // Radius at finest level 00067 vec<3,int> nx; // 3-D cell index 00068 if( end == 0 ) end = bodies.size(); // Default size is all bodies 00069 for( int b=begin; b!=end; ++b ) { // Loop over bodies 00070 for( int d=0; d!=3; ++d ) { // Loop over dimension 00071 nx[d] = int( ( bodies[b].X[d] - (X0[d]-R0) ) / r ); // 3-D cell index 00072 } // End loop over dimension 00073 i = 0; // Initialize cell index 00074 for( int l=0; l!=level; ++l ) { // Loop over levels of tree 00075 for( int d=0; d!=3; ++d ) { // Loop over dimension 00076 i += nx[d] % 2 << (3 * l + d); // Accumulate cell index 00077 nx[d] >>= 1; // Bitshift 3-D cell index 00078 } // End loop over dimension 00079 } // End loop over levels 00080 if( !update ) { // If this not an update 00081 bodies[b].ICELL = i+off; // Store index in bodies 00082 } else if( i+off > bodies[b].ICELL ) { // If the new cell index is larger 00083 bodies[b].ICELL = i+off; // Store index in bodies 00084 } // Endif for update 00085 } // End loop over bodies 00086 stopTimer("Set index ",printNow); // Stop timer 00087 } 00088 00089 //! Prune tree by merging cells 00090 void prune(Bodies &bodies) { 00091 startTimer("Prune tree "); // Start timer 00092 int maxLevel = getMaxLevel(bodies); // Max level for bottom up tree build 00093 for( int l=maxLevel; l>0; --l ) { // Loop upwards from bottom level 00094 int level = getLevel(bodies[0].ICELL); // Current level 00095 bigint cOff = ((1 << 3 * level) - 1) / 7; // Current ce;; index offset 00096 bigint pOff = ((1 << 3 * (l-1)) - 1) / 7; // Parent cell index offset 00097 bigint index = ((bodies[0].ICELL-cOff) >> 3*(level-l+1)) + pOff;// Current cell index 00098 int begin = 0; // Begin cell index for bodies in cell 00099 int size = 0; // Number of bodies in cell 00100 int b = 0; // Current body index 00101 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B,++b ) {// Loop over bodies 00102 level = getLevel(B->ICELL); // Level of twig 00103 cOff = ((1 << 3*level) - 1) / 7; // Offset of twig 00104 bigint p = ((B->ICELL-cOff) >> 3*(level-l+1)) + pOff; // Cell index of parent cell 00105 if( p != index ) { // If it's a new parent cell 00106 if( size < NCRIT ) { // If parent cell has few enough bodies 00107 for( int i=begin; i!=begin+size; ++i ) { // Loop over bodies in that cell 00108 bodies[i].ICELL = index; // Renumber cell index to parent cell 00109 } // End loop over bodies in cell 00110 } // Endif for merging 00111 begin = b; // Set new begin index 00112 size = 0; // Reset number of bodies 00113 index = p; // Set new parent cell 00114 } // Endif for new cell 00115 size++; // Increment body counter 00116 } // End loop over bodies 00117 if( size < NCRIT ) { // If last parent cell has few enough bodies 00118 for( int i=begin; i!=begin+size; ++i ) { // Loop over bodies in that cell 00119 bodies[i].ICELL = index; // Renumber cell index to parent cell 00120 } // End loop over bodies in cell 00121 } // Endif for merging 00122 } // End loop over levels 00123 stopTimer("Prune tree ",printNow); // Stop timer 00124 } 00125 00126 //! Grow tree by splitting cells 00127 void grow(Bodies &bodies, int level=0, int begin=0, int end=0) { 00128 bigint index = bodies[begin].ICELL; // Initialize cell index 00129 int off=begin, size=0; // Initialize offset, and size 00130 if( level == 0 ) level = getMaxLevel(bodies); // Max level for bottom up tree build 00131 if( end == 0 ) end = bodies.size(); // Default size is all bodies 00132 for( int b=begin; b!=end; ++b ) { // Loop over bodies under consideration 00133 if( bodies[b].ICELL != index ) { // If it's a new cell 00134 if( size >= NCRIT ) { // If the cell has too many bodies 00135 level++; // Increment level 00136 setIndex(bodies,level,off,off+size); // Set new cell index considering new level 00137 sortBodies(bodies,buffer,false,off,off+size); // Sort new cell index 00138 grow(bodies,level,off,off+size); // Recursively grow tree 00139 level--; // Go back to previous level 00140 } // Endif for splitting 00141 off = b; // Set new offset 00142 size = 0; // Reset number of bodies 00143 index = bodies[b].ICELL; // Set new cell 00144 } // Endif for new cell 00145 size++; // Increment body counter 00146 } // End loop over bodies 00147 if( size >= NCRIT ) { // If last cell has too many bodies 00148 level++; // Increment level 00149 setIndex(bodies,level,off,off+size); // Set new cell index considering new level 00150 sortBodies(bodies,buffer,false,off,off+size); // Sort new cell index 00151 grow(bodies,level,off,off+size); // Recursively grow tree 00152 level--; // Go back to previous level 00153 } // Endif for splitting 00154 } 00155 }; 00156 00157 #endif
