testUnifTensorialAlgorithm.cpp 13.6 KB
Newer Older
1
// ===================================================================================
2 3 4 5 6 7 8 9 10 11 12 13 14
// Copyright ScalFmm 2011 INRIA,
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
15
// ===================================================================================
16 17 18
// author P. Banchard
// Modifs
//  O. Coulaud
19
// ==== CMAKE =====
20
// @FUSE_FFT
21 22 23 24 25 26 27
// ================

#include <iostream>

#include <cstdio>
#include <cstdlib>

28
#include "Files/FFmaGenericLoader.hpp"
29 30


31 32 33
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Uniform/FUnifTensorialKernel.hpp"
34

35 36
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
37

38 39
#include "Utils/FParameters.hpp"
#include "Utils/FMemUtils.hpp"
40

41 42
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"
43

44 45
#include "Core/FFmmAlgorithm.hpp"
#include "Core/FFmmAlgorithmThread.hpp"
46 47 48 49 50 51 52 53 54

/**
 * This program runs the FMM Algorithm with the Uniform kernel and compares the results with a direct computation.
 */

// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
  const char* const filename       = FParameters::getStr(argc,argv,"-f", "../Data/test20k.fma");
55 56
  const unsigned int TreeHeight    = FParameters::getValue(argc, argv, "-depth", 3);
  const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, "-subdepth", 2);
57 58 59 60 61 62 63 64 65 66 67 68
  const unsigned int NbThreads     = FParameters::getValue(argc, argv, "-t", 1);

#ifdef _OPENMP
  omp_set_num_threads(NbThreads);
  std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
#else
  std::cout << "\n>> Sequential version.\n" << std::
#endif

    // init timer
    FTic time;

69 70

  // typedefs
71 72
  typedef FInterpMatrixKernel_R_IJ MatrixKernelClass;
//  typedef FInterpMatrixKernel_R_IJK MatrixKernelClass;
73

74
  // useful features of matrix kernel
75
  const KERNEL_FUNCTION_IDENTIFIER MK_ID = MatrixKernelClass::Identifier;
76 77
  const unsigned int NPV  = MatrixKernelClass::NPV;
  const unsigned int NPOT = MatrixKernelClass::NPOT;
78 79 80
  const unsigned int NRHS = MatrixKernelClass::NRHS;
  const unsigned int NLHS = MatrixKernelClass::NLHS;

81
  const double CoreWidth = 0.1;
82
  const MatrixKernelClass DirectMatrixKernel(CoreWidth);
83 84 85
  std::cout<< "Interaction kernel: ";
  if(MK_ID == ONE_OVER_R) std::cout<< "1/R";
  else if(MK_ID == R_IJ)
86 87
    std::cout<< "Ra_{,ij}" 
             << " with core width a=" << CoreWidth <<std::endl;
88
  else if(MK_ID == R_IJK)
89 90
    std::cout<< "Ra_{,ijk} (BEWARE! Forces NOT YET IMPLEMENTED)" 
             << " with core width a=" << CoreWidth <<std::endl;
91
  else std::runtime_error("Provide a valid matrix kernel!");
92

93 94 95
  // init particles position and physical value
  struct TestParticle{
    FPoint position;
96 97 98
    FReal forces[3][NPOT];
    FReal physicalValue[NPV];
    FReal potential[NPOT];
99 100 101
  };

  // open particle file
102 103
  FFmaGenericLoader loader(filename);

104 105 106 107 108 109 110
  if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");

  TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
  for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
    FPoint position;
    FReal physicalValue = 0.0;
    loader.fillParticle(&position,&physicalValue);
111

112 113
    // get copy
    particles[idxPart].position       = position;
114 115 116 117 118
    // Set physical values
    for(unsigned idxPV = 0; idxPV<NPV;++idxPV){
//    //   Either copy same physical value in each component
        particles[idxPart].physicalValue[idxPV]  = physicalValue;
    // ... or set random value
119
//      particles[idxPart].physicalValue[idxPV]  = physicalValue*FReal(drand48());
120 121 122 123 124 125
    }
    for(unsigned idxPot = 0; idxPot<NPOT;++idxPot){
      particles[idxPart].potential[idxPot]      = 0.0;
      particles[idxPart].forces[0][idxPot]      = 0.0;
      particles[idxPart].forces[1][idxPot]      = 0.0;
      particles[idxPart].forces[2][idxPot]      = 0.0;
126
    }
127 128 129 130 131 132 133 134 135 136
  }

  ////////////////////////////////////////////////////////////////////

  { // begin direct computation
    std::cout << "\nDirect Computation ... " << std::endl;
    time.tic();
    {
      for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
        for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
137 138 139 140 141 142 143 144
          if(MK_ID == R_IJ)
            FP2P::MutualParticlesRIJ(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
                                     particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
                                     particles[idxTarget].forces[0], particles[idxTarget].forces[1],
                                     particles[idxTarget].forces[2], particles[idxTarget].potential,
                                     particles[idxOther].position.getX(), particles[idxOther].position.getY(),
                                     particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
                                     particles[idxOther].forces[0], particles[idxOther].forces[1],
145 146 147 148 149 150 151 152 153 154 155 156
                                     particles[idxOther].forces[2], particles[idxOther].potential,
                                     &DirectMatrixKernel);
          else if(MK_ID == R_IJK)
            FP2P::MutualParticlesRIJK(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
                                      particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
                                      particles[idxTarget].forces[0], particles[idxTarget].forces[1],
                                      particles[idxTarget].forces[2], particles[idxTarget].potential,
                                      particles[idxOther].position.getX(), particles[idxOther].position.getY(),
                                      particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
                                      particles[idxOther].forces[0], particles[idxOther].forces[1],
                                      particles[idxOther].forces[2], particles[idxOther].potential,
                                      &DirectMatrixKernel);
157 158
          else 
            std::runtime_error("Provide a valid matrix kernel!");
159 160 161 162 163 164 165 166 167 168 169 170 171 172
        }
      }
    }
    time.tac();
    std::cout << "Done  " << "(@Direct Computation = "
              << time.elapsed() << "s)." << std::endl;

  } // end direct computation

  ////////////////////////////////////////////////////////////////////

  {	// begin Lagrange kernel

    // accuracy
173
    const unsigned int ORDER = 6 ;
174 175 176 177 178 179 180 181 182 183 184 185 186 187
    // set box width extension
    // ... either deduce from element size
    const FReal LeafCellWidth = FReal(loader.getBoxWidth()) / FReal(FMath::pow(2.,TreeHeight-1));
    const FReal ElementSize = LeafCellWidth / FReal(3.);  
    const FReal BoxWidthExtension = ElementSize; // depends on type of element
    // ... or set to arbitrary value (0. means no extension)
//    const FReal BoxWidthExtension = FReal(0.); 

    std::cout << "LeafCellWidth=" << LeafCellWidth 
              << ", BoxWidthExtension=" << BoxWidthExtension <<std::endl;

    // stop execution if interactions are homog and box extension is required 
    if(MatrixKernelClass::Type==HOMOGENEOUS && BoxWidthExtension>0.)
      throw std::runtime_error("Extension of box width is not yet supported for homogeneous kernels! Work-around: artificially set Type to NON_HOMOGENEOUS.");
188

189
    typedef FP2PParticleContainerIndexed<NRHS,NLHS> ContainerClass;
190 191 192 193 194 195 196 197 198 199 200 201 202 203

    typedef FSimpleLeaf< ContainerClass >  LeafClass;
    typedef FUnifCell<ORDER,NRHS,NLHS> CellClass;
    typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
    typedef FUnifTensorialKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
    typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
    //  typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;

    // init oct-tree
    OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());


    { // -----------------------------------------------------
      std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
204 205 206 207 208
                << " particles ..." << std::endl; 




209 210 211 212 213
      std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
      time.tic();

      for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
        // put in tree
214 215 216 217 218 219 220 221 222 223
        if(NPV==3) // R_IJ or IOR
          tree.insert(particles[idxPart].position, idxPart, 
                      particles[idxPart].physicalValue[0], particles[idxPart].physicalValue[1], particles[idxPart].physicalValue[2]);
        else if(NPV==9) // R_IJK
          tree.insert(particles[idxPart].position, idxPart, 
                      particles[idxPart].physicalValue[0], particles[idxPart].physicalValue[1], particles[idxPart].physicalValue[2],
                      particles[idxPart].physicalValue[3], particles[idxPart].physicalValue[4], particles[idxPart].physicalValue[5],
                      particles[idxPart].physicalValue[6], particles[idxPart].physicalValue[7], particles[idxPart].physicalValue[8]);
        else
          std::runtime_error("NPV not yet supported in test! Add new case.");
224 225 226 227 228 229 230 231 232 233
      }

      time.tac();
      std::cout << "Done  " << "(@Creating and Inserting Particles = "
                << time.elapsed() << "s)." << std::endl;
    } // -----------------------------------------------------

    { // -----------------------------------------------------
      std::cout << "\nLagrange/Uniform grid FMM (ORDER="<< ORDER << ") ... " << std::endl;
      time.tic();
234
      KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),BoxWidthExtension,CoreWidth);
235 236 237 238 239 240 241 242 243
      FmmClass algorithm(&tree, &kernels);
      algorithm.execute();
      time.tac();
      std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
    } // -----------------------------------------------------


    { // -----------------------------------------------------
      std::cout << "\nError computation ... " << std::endl;
244 245
      FMath::FAccurater potentialDiff[NPOT];
      FMath::FAccurater fx[NPOT], fy[NPOT], fz[NPOT];
246

247 248
      FReal checkPotential[20000][NPOT];
      FReal checkfx[20000][NPOT];
249 250 251 252

      { // Check that each particle has been summed with all other

        tree.forEachLeaf([&](LeafClass* leaf){
253
            for(unsigned idxPot = 0; idxPot<NPOT;++idxPot){
254

255 256 257 258
              const FReal*const potentials = leaf->getTargets()->getPotentials(idxPot);
              const FReal*const forcesX = leaf->getTargets()->getForcesX(idxPot);
              const FReal*const forcesY = leaf->getTargets()->getForcesY(idxPot);
              const FReal*const forcesZ = leaf->getTargets()->getForcesZ(idxPot);
259 260 261 262 263 264 265
              const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
              const FVector<int>& indexes = leaf->getTargets()->getIndexes();

              for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                const int indexPartOrig = indexes[idxPart];

                //PB: store potential in array[nbParticles]
266 267 268 269 270 271 272
                checkPotential[indexPartOrig][idxPot]=potentials[idxPart];              
                checkfx[indexPartOrig][idxPot]=forcesX[idxPart];              

                potentialDiff[idxPot].add(particles[indexPartOrig].potential[idxPot],potentials[idxPart]);
                fx[idxPot].add(particles[indexPartOrig].forces[0][idxPot],forcesX[idxPart]);
                fy[idxPot].add(particles[indexPartOrig].forces[1][idxPot],forcesY[idxPart]);
                fz[idxPot].add(particles[indexPartOrig].forces[2][idxPot],forcesZ[idxPart]);
273
              }
274
            }// NPOT
275 276 277 278 279
          });
      }


      // Print for information
280 281 282 283 284 285 286 287
      std::cout << "\nRelative Inf/L2 errors: " << std::endl;
      std::cout << "  Potential: " << std::endl;
      for(unsigned idxPot = 0; idxPot<NPOT;++idxPot) {
        std::cout << "    " << idxPot << ": "
                  << potentialDiff[idxPot].getRelativeInfNorm() << ", " 
                  << potentialDiff[idxPot].getRelativeL2Norm() 
                  << std::endl;
      }
288
      std::cout << std::endl;
289 290 291 292 293 294 295
      std::cout << "  Fx: " << std::endl; 
      for(unsigned idxPot = 0; idxPot<NPOT;++idxPot) {
        std::cout << "    " << idxPot << ": "
                  << fx[idxPot].getRelativeInfNorm() << ", " 
                  << fx[idxPot].getRelativeL2Norm()
                  << std::endl;
      }
296
      std::cout  << std::endl;
297 298 299 300 301 302 303
      std::cout << "  Fy: " << std::endl; 
      for(unsigned idxPot = 0; idxPot<NPOT;++idxPot) {
        std::cout << "    " << idxPot << ": "
                  << fy[idxPot].getRelativeInfNorm() << ", " 
                  << fy[idxPot].getRelativeL2Norm()
                  << std::endl;
      }
304
      std::cout  << std::endl;
305 306 307 308 309 310 311
      std::cout << "  Fz: " << std::endl; 
      for(unsigned idxPot = 0; idxPot<NPOT;++idxPot) {
        std::cout << "    " << idxPot << ": "
                  << fz[idxPot].getRelativeInfNorm() << ", " 
                  << fz[idxPot].getRelativeL2Norm()
                  << std::endl;
      }
312 313
      std::cout << std::endl;

314 315 316 317 318 319
    } // -----------------------------------------------------

  } // end Lagrange kernel

  return 0;
}