ExaFMM 1
Fast-multipole Method for exascale systems
|
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 #include "parallelfmm.h" 00023 #ifdef VTK 00024 #include "vtk.h" 00025 #endif 00026 00027 int main() { 00028 int numBodies = 10000; // Number of bodies 00029 int numTarget = 100; // Number of target points to be used for error eval. 00030 IMAGES = 0; // Level of periodic image tree (0 for non-periodic FMM) 00031 THETA = 1 / sqrtf(4); // Multipole acceptace criteria 00032 Bodies bodies, jbodies; // Define vector of bodies 00033 Cells cells, jcells; // Define vector of cells 00034 ParallelFMM<Laplace> FMM; // Instantiate ParallelFMM class 00035 FMM.initialize(); // Initialize FMM 00036 bool printNow = MPIRANK == 0; // Print only if MPIRANK == 0 00037 00038 for( int it=0; it!=25; ++it ) { // Loop from N = 10^4 to 10^7 00039 numBodies = int(pow(10,(it+32)/8.0)); // Exponentially increase N 00040 if(printNow) std::cout << "N : " << numBodies << std::endl;// Print N 00041 bodies.resize(numBodies); // Resize bodies vector 00042 FMM.random(bodies,MPIRANK+1); // Initialize bodies with random coordinates 00043 FMM.startTimer("FMM "); // Start timer 00044 FMM.setGlobDomain(bodies); // Set global domain size of FMM 00045 FMM.octsection(bodies); // Partition domain and redistribute bodies 00046 cells.clear(); // Make sure cells vector is empty 00047 #ifdef TOPDOWN 00048 FMM.topdown(bodies,cells); // Tree construction (top down) & upward sweep 00049 #else 00050 FMM.bottomup(bodies,cells); // Tree construction (bottom up) & upward sweep 00051 #endif 00052 FMM.commBodies(cells); // Send bodies (not receiving yet) 00053 jbodies = bodies; // Vector of source bodies 00054 jcells = cells; // Vector of source cells 00055 FMM.commCells(jbodies,jcells); // Communicate cells (receive bodies here) 00056 00057 FMM.downward(cells,jcells); // Downward sweep 00058 FMM.stopTimer("FMM ",printNow); // Stop timer 00059 FMM.eraseTimer("FMM "); // Erase entry from timer to avoid timer overlap 00060 00061 FMM.startTimer("Direct sum "); // Start timer 00062 Bodies bodies2 = bodies; // Define new bodies vector for direct sum 00063 #if 1 00064 FMM.initTarget(bodies2); // Reset target values to 0 00065 if( IMAGES != 0 ) { // For periodic boundary condition 00066 jbodies = FMM.periodicBodies(bodies2); // Copy source bodies for all periodic images 00067 } else { // For free field boundary condition 00068 jbodies = bodies2; // Source bodies = target bodies 00069 } // End if for periodic boundary condition 00070 bodies2.resize(numTarget); // Shrink target bodies vector to save time 00071 for( int i=0; i!=MPISIZE; ++i ) { // Loop over all MPI processes 00072 FMM.shiftBodies(jbodies); // Communicate bodies round-robin 00073 FMM.evalP2P(bodies2,jbodies); // Direct summation between bodies2 and jbodies 00074 if(FMM.printNow) std::cout << "Direct loop : " << i+1 << "/" << MPISIZE << std::endl;// Print loop counter 00075 } // End loop over all MPI processes 00076 FMM.writeTarget(bodies2); // Write direct summation results to file 00077 #else 00078 FMM.readTarget(bodies2); // Read direct summation results from file 00079 #endif 00080 FMM.stopTimer("Direct sum ",printNow); // Stop timer 00081 FMM.eraseTimer("Direct sum "); // Erase entry from timer to avoid timer overlap 00082 if(printNow) FMM.writeTime(); // Write timings of all events to file 00083 if(printNow) FMM.resetTimer(); // Erase all events in timer 00084 00085 real diff1 = 0, norm1 = 0, diff2 = 0, norm2 = 0, diff3 = 0, norm3 = 0, diff4 = 0, norm4 = 0; 00086 bodies.resize(numTarget); // Shrink bodies to match bodies2 00087 FMM.evalError(bodies,bodies2,diff1,norm1,diff2,norm2); // Evaluate error on the reduced set of bodies 00088 MPI_Datatype MPI_TYPE = FMM.getType(diff1); // Get MPI datatype 00089 MPI_Reduce(&diff1,&diff3,1,MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD);// Reduce difference in potential 00090 MPI_Reduce(&norm1,&norm3,1,MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD);// Reduce norm of potential 00091 MPI_Reduce(&diff2,&diff4,1,MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD);// Reduce difference in force 00092 MPI_Reduce(&norm2,&norm4,1,MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD);// Recude norm of force 00093 if(printNow) FMM.printError(diff3,norm3,diff4,norm4); // Print the L2 norm error of potential & force 00094 } // End loop over N 00095 FMM.finalize(); // Finalize FMM 00096 }