sharedMemoryInterpolationFMM.hpp 10 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
// -*-c++-*-
// See LICENCE file at project root

// ==== CMAKE =====
// @FUSE_FFT
// @FUSE_BLAS
//  ==== Git =====

// ================

/** \brief Uniform FMM example
 *
 * \file
 * \authors B. Bramas, O. Coulaud
 *
 * This program runs the FMM Algorithm with the interpolation kernel based on
 * either uniform (grid points) or Chebychev interpolation (1/r kernel). 
 * It then compares the results with a direct computation.
 */

#include <iostream>
#include <iomanip>
#include <memory>

#include <cstdio>
#include <cstdlib>
#include <string>


#include "ScalFmmConfig.h"
#include "Utils/FGlobal.hpp"

#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"

#include "Files/FFmaGenericLoader.hpp"
// Leaves
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
// Octree
#include "Containers/FOctree.hpp"
//
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
//
//
// Order of the Interpolation approximation
47
static constexpr unsigned ORDER = 6 ;
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
using FReal                 = double;
//   1/r kernel
//
using MatrixKernelClass     = FInterpMatrixKernelR<FReal> ;

//
/// definition of the common tree structure
using CellClass      = FInterpolationCell<FReal, ORDER>;
//using CellUpClass    = typename CellClass::multipole_t;
//using CellDownClass  = typename CellClass::local_expansion_t;
//using CellSymbClass  = FSymbolicData;  
using ContainerClass = FP2PParticleContainerIndexed<FReal>;
using LeafClass      = FSimpleLeaf<FReal,  ContainerClass >   ;
using OctreeClass    = FOctree<FReal, CellClass,ContainerClass,LeafClass>;
using KernelClass    = FInterpolationKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> ;


#ifdef _OPENMP
66 67
//#include "Core/FFmmAlgorithmThread.hpp"
#include "Core/FFmmAlgorithmPeriodic.hpp"
68 69 70 71 72 73 74 75
#include "Core/FFmmAlgorithmSectionTask.hpp"
#else
#include "Core/FFmmAlgorithm.hpp"
#endif


#ifdef _OPENMP
//using FmmClass = FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
76 77
using FmmClass    = FFmmAlgorithmSectionTask<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> ;
using FmmClassPer = FFmmAlgorithmPeriodic<FReal,OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass>;
78
#else
79
using FmmClass  = FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass>;
80 81 82 83 84 85 86 87 88 89 90 91 92 93
#endif



// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
  FHelpDescribeAndExit(
		       argc, argv,
		       "Driver for Lagrange interpolation kernel  (1/r kernel).",
		       FParameterDefinitions::OctreeHeight,
		       FParameterDefinitions::OctreeSubHeight,
		       FParameterDefinitions::InputFile,
		       FParameterDefinitions::OutputFile,
94 95
		       FParameterDefinitions::NbThreads,
		       FParameterDefinitions::PeriodicityNbLevels
96 97 98 99
		       );


  const std::string defaultFile(SCALFMMDataPath+"unitCubeXYZQ100.bfma" );
100 101 102
  const std::string filename       = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFile.c_str());
  const int TreeHeight    = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
  const int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
103 104 105 106
  bool periodicCondition = false ;
  if(FParameters::existParameter(argc, argv, FParameterDefinitions::PeriodicityNbLevels.options)){
      periodicCondition = true;
    }
107
  const unsigned int aboveTree = FParameters::getValue(argc, argv, FParameterDefinitions::PeriodicityNbLevels.options, -1);
108 109

#ifdef _OPENMP
110 111 112 113
  const unsigned int NbThreads = periodicCondition ? 1 : FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, omp_get_max_threads());

   omp_set_num_threads(NbThreads);

114 115
  std::cout << "\n>> Using " << NbThreads << " threads.\n" << std::endl;
#else
116
  const int NbThreads =  1;
117 118 119 120 121 122 123
  std::cout << "\n>> Sequential version.\n" << std::endl;
#endif
  //
  {
    std::string indent("    ");
    auto w = std::setw(18);
    std::cout << "Parameters" << std::endl << std::left
124 125 126 127 128 129 130 131 132
              << indent << w << "Octree Depth: "     << TreeHeight    << std::endl
              << indent << w << "SubOctree depth: "  << SubTreeHeight << std::endl;
    if(periodicCondition){
        std::cout << indent << w << "AboveTree    "<< aboveTree <<std::endl;

      }
    std::cout << indent << w << "Input file  name: " << filename      << std::endl
              << indent << w << "Thread number: "    << NbThreads     << std::endl
              << std::endl;
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
  }
  //
  // init timer
  FTic time;

  // open particle file
  ////////////////////////////////////////////////////////////////////
  //
  FFmaGenericLoader<FReal> loader(filename);
  //
  // init oct-tree
  OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
  //
  // Read particles and insert them in octree
  { // -----------------------------------------------------
    std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
	      << " particles ..." << std::endl;
    std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
    time.tic();
    //
    FPoint<FReal> position;
    FReal physicalValue = 0.0;
    //
    for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
      //
      // Read particle per particle from file
      loader.fillParticle(&position,&physicalValue);
      //
      // put particle in octree
      tree.insert(position, idxPart, physicalValue);
    }

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

175
 // { // -----------------------------------------------------
176 177 178 179 180
    std::cout << "\n" << interpolationType << "  FMM (ORDER= "<< ORDER << ") ... " << std::endl;

    const MatrixKernelClass    MatrixKernel;
    time.tic();
    //
181 182 183 184 185 186 187 188
    std::unique_ptr<KernelClass> kernelsNoPer(new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel));
    //
    FmmClass algoNoPer(&tree, kernelsNoPer.get());
    // periodic FMM algorithm
    FmmClassPer algoPer(&tree, aboveTree);
    KernelClass kernelsPer(algoPer.extendedTreeHeight(), algoPer.extendedBoxWidth(),
                           algoPer.extendedBoxCenter(),&MatrixKernel);
    algoPer.setKernel(&kernelsPer);
189
    //
190 191 192 193 194 195 196 197 198 199 200
    FAbstractAlgorithm * algorithm  = nullptr;
    FAlgorithmTimers   * timer      = nullptr;
    if(! periodicCondition) { // Non periodic case
        algorithm  = &algoNoPer ;
        timer      = &algoNoPer ;
      }
    else {                    // Periodic case
        algorithm  = &algoPer ;
        timer      = &algoPer ;

      }
COULAUD Olivier's avatar
COULAUD Olivier committed
201

202
    //
203
    algorithm->execute();   // Here the call of the FMM algorithm
204 205 206
    //
    time.tac();
    std::cout << "Timers Far Field \n"
207 208 209 210 211
              << "P2M " << timer->getTime(FAlgorithmTimers::P2MTimer) << " seconds\n"
              << "M2M " << timer->getTime(FAlgorithmTimers::M2MTimer) << " seconds\n"
              << "M2L " << timer->getTime(FAlgorithmTimers::M2LTimer) << " seconds\n"
              << "L2L " << timer->getTime(FAlgorithmTimers::L2LTimer) << " seconds\n"
              << "P2P and L2P " << timer->getTime(FAlgorithmTimers::NearTimer) << " seconds\n"
212 213 214 215
	      << std::endl;


    std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
216
  //}
217 218 219 220 221 222 223 224 225 226 227 228 229
  // -----------------------------------------------------
  //
  // Some output
  //
  //
  { // -----------------------------------------------------
    FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= loader.getNumberOfParticles() -1; ;
    FReal energy =0.0 ;
    //
    //   Loop over all leaves
    //
    std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
    std::cout << std::scientific;
230 231 232
    std::cout.precision(15) ;
//     std::size_t count =0 , numLeaf =0 ;
    FReal TotalPhysicalValue=0.0 ;
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
    tree.forEachLeaf([&](LeafClass* leaf){
	const FReal*const posX = leaf->getTargets()->getPositions()[0];
	const FReal*const posY = leaf->getTargets()->getPositions()[1];
	const FReal*const posZ = leaf->getTargets()->getPositions()[2];

	const FReal*const potentials = leaf->getTargets()->getPotentials();
	const FReal*const forcesX = leaf->getTargets()->getForcesX();
	const FReal*const forcesY = leaf->getTargets()->getForcesY();
	const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
	const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
	const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();

	const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();

	for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
	  const FSize indexPartOrig = indexes[idxPart];
	  if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3)  ) {
	    std::cout << "Index "<< indexPartOrig <<"  potential  " << potentials[idxPart]
		      << " Pos "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
		      << "   Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
	  }
	  energy += potentials[idxPart]*physicalValues[idxPart] ;
255 256 257
	  TotalPhysicalValue += physicalValues[idxPart]  ;
	//  ++count ;

258
	}
259 260
//	std::cout << "numLeaf "<< numLeaf << "  " << energy <<"  " << count<<std::endl<<std::endl ;
//	++numLeaf ;
261
      });
262
    std::cout <<std::endl<<"aboveRoot: " << aboveTree << "  Energy: "<< energy<<"  TotalPhysicalValue: " << TotalPhysicalValue<< std::endl;
263 264 265
    std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;

  }
266
  //
267
  // -----------------------------------------------------
268
  //
269 270
  if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
    std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options,   "output.fma"));
271
    
272
    FFmaGenericWriter<FReal> writer(name) ;
273
    writer.writeDataFromOctree(&tree,loader.getNumberOfParticles());
274 275 276 277 278 279

  }


  return 0;
}