|
ExaFMM 1
Fast-multipole Method for exascale systems
|
include/dataset.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 dataset_h 00023 #define dataset_h 00024 #include "kernel.h" 00025 00026 //! Contains all the different datasets 00027 template<Equation equation> 00028 class Dataset : public Kernel<equation> { 00029 private: 00030 long filePosition; //!< Position of file stream 00031 00032 public: 00033 //! Constructor 00034 Dataset() : filePosition(0) {} 00035 //! Destructor 00036 ~Dataset() {} 00037 00038 //! Initialize source values 00039 void initSource(Bodies &bodies) { 00040 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) { // Loop over bodies 00041 B->IBODY = B-bodies.begin(); // Tag body with initial index 00042 B->IPROC = MPIRANK; // Tag body with initial MPI rank 00043 B->SRC = 1. / bodies.size() / MPISIZE; // Initialize mass/charge 00044 } // End loop over bodies 00045 } 00046 00047 //! Initialize target values 00048 void initTarget(Bodies &bodies, bool IeqJ=true) { 00049 srand48(0); // Set seed for random number generator 00050 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) { // Loop over bodies 00051 B->IBODY = B-bodies.begin(); // Tag body with initial index 00052 B->IPROC = MPIRANK; // Tag body with initial MPI rank 00053 B->TRG = 0; // Clear previous target values (IeqJ is dummy) 00054 if( EPS2 != 0 ) B->TRG[0] = -B->SRC / std::sqrt(EPS2) * IeqJ;// Initialize potential (0 if I != J) 00055 } // End loop over bodies 00056 } 00057 00058 //! Read target values from file 00059 void readTarget(Bodies &bodies) { 00060 char fname[256]; // File name for saving direct calculation values 00061 sprintf(fname,"direct%4.4d",MPIRANK); // Set file name 00062 std::ifstream file(fname,std::ios::in | std::ios::binary); // Open file 00063 file.seekg(filePosition); // Set position in file 00064 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) { // Loop over bodies 00065 file >> B->TRG[0]; // Read data for potential 00066 file >> B->TRG[1]; // Read data for x acceleration 00067 file >> B->TRG[2]; // Read data for y acceleration 00068 file >> B->TRG[3]; // Read data for z acceleration 00069 } // End loop over bodies 00070 filePosition = file.tellg(); // Get position in file 00071 file.close(); // Close file 00072 } 00073 00074 //! Write target values to file 00075 void writeTarget(Bodies &bodies) { 00076 char fname[256]; // File name for saving direct calculation values 00077 sprintf(fname,"direct%4.4d",MPIRANK); // Set file name 00078 std::ofstream file(fname,std::ios::out | std::ios::app | std::ios::binary);// Open file 00079 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) { // Loop over bodies 00080 file << B->TRG[0] << std::endl; // Write data for potential 00081 file << B->TRG[1] << std::endl; // Write data for x acceleration 00082 file << B->TRG[2] << std::endl; // Write data for y acceleration 00083 file << B->TRG[3] << std::endl; // Write data for z acceleration 00084 } // End loop over bodies 00085 file.close(); // Close file 00086 } 00087 00088 //! Evaluate relative L2 norm error 00089 void evalError(Bodies &bodies, Bodies &bodies2, 00090 real &diff1, real &norm1, real &diff2, real &norm2, bool ewald=false) { 00091 real p = 0, p2 = 0; // Total energy 00092 B_iter B2 = bodies2.begin(); // Set iterator for bodies2 00093 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B, ++B2 ) {// Loop over bodies & bodies2 00094 #ifdef DEBUG 00095 std::cout << B->ICELL << " " << B->TRG[0] << " " << B2->TRG[0] << std::endl;// Compare every element 00096 #endif 00097 if( ewald ) { // If ewald method is used 00098 p += B->TRG[0] * B->SRC; // total energy 00099 p2 += B2->TRG[0] * B2->SRC; // total energy 00100 } else { // If normal Laplace kernel is used 00101 diff1 += (B->TRG[0] - B2->TRG[0]) * (B->TRG[0] - B2->TRG[0]);// Difference of potential 00102 norm1 += B2->TRG[0] * B2->TRG[0]; // Value of potential 00103 } // End if for Ewald method 00104 diff2 += (B->TRG[1] - B2->TRG[1]) * (B->TRG[1] - B2->TRG[1]);// Difference of x acceleration 00105 diff2 += (B->TRG[2] - B2->TRG[2]) * (B->TRG[2] - B2->TRG[2]);// Difference of y acceleration 00106 diff2 += (B->TRG[3] - B2->TRG[3]) * (B->TRG[3] - B2->TRG[3]);// Difference of z acceleration 00107 norm2 += B2->TRG[1] * B2->TRG[1]; // Value of x acceleration 00108 norm2 += B2->TRG[2] * B2->TRG[2]; // Value of y acceleration 00109 norm2 += B2->TRG[3] * B2->TRG[3]; // Value of z acceleration 00110 } // End loop over bodies & bodies2 00111 if( ewald ) { // If ewald method is used 00112 diff1 = (p - p2) * (p - p2); // Difference of total energy 00113 norm1 = p2 * p2; // Value of total energy 00114 } // End if for Ewald method 00115 } 00116 00117 //! Print relative L2 norm error 00118 void printError(real diff1, real norm1, real diff2, real norm2) { 00119 std::cout << "Error (pot) : " << std::sqrt(diff1/norm1) << std::endl; 00120 std::cout << "Error (acc) : " << std::sqrt(diff2/norm2) << std::endl; 00121 } 00122 }; 00123 00124 template<> 00125 class Dataset<VanDerWaals> : public Kernel<VanDerWaals> { 00126 private: 00127 long filePosition; //!< Position of file stream 00128 00129 public: 00130 //! Constructor 00131 Dataset() : filePosition(0) {} 00132 //! Destructor 00133 ~Dataset() {} 00134 00135 //! Initialize source values 00136 void initSource(Bodies &bodies) { 00137 THETA = .1; // Force opening angle to be small 00138 ATOMS = 16; // Set number of atoms 00139 RSCALE.resize(ATOMS*ATOMS); // Resize rscale vector 00140 GSCALE.resize(ATOMS*ATOMS); // Resize gscale vector 00141 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) { // Loop over bodies 00142 B->IBODY = B-bodies.begin(); // Tag body with initial index 00143 B->IPROC = MPIRANK; // Tag body with initial MPI rank 00144 B->SRC = drand48() * ATOMS; // Initialize mass/charge 00145 } // End loop over bodies 00146 for( int i=0; i!=ATOMS; ++i ) { // Loop over atoms 00147 GSCALE[i*ATOMS+i] = drand48(); // Set VanDerWaals post scaling factor 00148 RSCALE[i*ATOMS+i] = drand48(); // Set VanDerWaals pre scaling factor 00149 } // End loop over atoms 00150 for( int i=0; i!=ATOMS; ++i ) { // Loop over target atoms 00151 for( int j=0; j!=ATOMS; ++j ) { // Loop over source atoms 00152 if( i != j ) { // If target and source are different 00153 GSCALE[i*ATOMS+j] = std::sqrt(GSCALE[i*ATOMS+i] * GSCALE[j*ATOMS+j]);// Set post scaling factor 00154 RSCALE[i*ATOMS+j] = (std::sqrt(RSCALE[i*ATOMS+i]) + std::sqrt(RSCALE[j*ATOMS+j])) * 0.5; 00155 RSCALE[i*ATOMS+j] *= RSCALE[i*ATOMS+j]; // Set pre scaling factor 00156 } // End if for different target and source 00157 } // End loop over source atoms 00158 } // End loop over target atoms 00159 } 00160 00161 //! Initialize target values 00162 void initTarget(Bodies &bodies, bool IeqJ=true) { 00163 srand48(0); // Set seed for random number generator 00164 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) { // Loop over bodies 00165 B->IBODY = B-bodies.begin(); // Tag body with initial index 00166 B->IPROC = MPIRANK; // Tag body with initial MPI rank 00167 B->TRG = 0 * IeqJ; // Clear previous target values 00168 } // End loop over bodies 00169 } 00170 00171 //! Read target values from file 00172 void readTarget(Bodies &bodies) { 00173 char fname[256]; // File name for saving direct calculation values 00174 sprintf(fname,"direct%4.4d",MPIRANK); // Set file name 00175 std::ifstream file(fname,std::ios::in | std::ios::binary); // Open file 00176 file.seekg(filePosition); // Set position in file 00177 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) { // Loop over bodies 00178 file >> B->TRG[0]; // Read data for potential 00179 file >> B->TRG[1]; // Read data for x acceleration 00180 file >> B->TRG[2]; // Read data for y acceleration 00181 file >> B->TRG[3]; // Read data for z acceleration 00182 } // End loop over bodies 00183 filePosition = file.tellg(); // Get position in file 00184 file.close(); // Close file 00185 } 00186 00187 //! Write target values to file 00188 void writeTarget(Bodies &bodies) { 00189 char fname[256]; // File name for saving direct calculation values 00190 sprintf(fname,"direct%4.4d",MPIRANK); // Set file name 00191 std::ofstream file(fname,std::ios::out | std::ios::app | std::ios::binary);// Open file 00192 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) { // Loop over bodies 00193 file << B->TRG[0] << std::endl; // Write data for potential 00194 file << B->TRG[1] << std::endl; // Write data for x acceleration 00195 file << B->TRG[2] << std::endl; // Write data for y acceleration 00196 file << B->TRG[3] << std::endl; // Write data for z acceleration 00197 } // End loop over bodies 00198 file.close(); // Close file 00199 } 00200 00201 //! Evaluate relative L2 norm error 00202 void evalError(Bodies &bodies, Bodies &bodies2, 00203 real &diff1, real &norm1, real &diff2, real &norm2) { 00204 B_iter B2 = bodies2.begin(); // Set iterator for bodies2 00205 for( B_iter B=bodies.begin(); B!=bodies.end(); ++B, ++B2 ) {// Loop over bodies & bodies2 00206 #ifdef DEBUG 00207 std::cout << B->ICELL << " " << B->TRG[0] << " " << B2->TRG[0] << std::endl;// Compare every element 00208 #endif 00209 diff1 += (B->TRG[0] - B2->TRG[0]) * (B->TRG[0] - B2->TRG[0]);// Difference of potential 00210 norm1 += B2->TRG[0] * B2->TRG[0]; // Value of potential 00211 diff2 += (B->TRG[1] - B2->TRG[1]) * (B->TRG[1] - B2->TRG[1]);// Difference of x acceleration 00212 diff2 += (B->TRG[2] - B2->TRG[2]) * (B->TRG[2] - B2->TRG[2]);// Difference of y acceleration 00213 diff2 += (B->TRG[3] - B2->TRG[3]) * (B->TRG[3] - B2->TRG[3]);// Difference of z acceleration 00214 norm2 += B2->TRG[1] * B2->TRG[1]; // Value of x acceleration 00215 norm2 += B2->TRG[2] * B2->TRG[2]; // Value of y acceleration 00216 norm2 += B2->TRG[3] * B2->TRG[3]; // Value of z acceleration 00217 } // End loop over bodies & bodies2 00218 } 00219 00220 //! Print relative L2 norm error 00221 void printError(real diff1, real norm1, real diff2, real norm2) { 00222 std::cout << "Error (pot) : " << std::sqrt(diff1/norm1) << std::endl; 00223 std::cout << "Error (acc) : " << std::sqrt(diff2/norm2) << std::endl; 00224 } 00225 }; 00226 00227 #endif
