NEw files to check Algorithm with Interpolation method

Utils/src/ChebyshevInterpolationCmpAlgo.cpp

+// ==== CMAKE =====
+// @FUSE_BLAS
+// ================
+// Keep in private GIT
+//
+//
+/** \brief Chebyshev FMM example
+ *
+ * \file
+ * \authors B. Bramas, O. Coulaud
+ *
+ * This program runs the FMM Algorithm with the interpolation kernel based on
+ * Chebyshev (grid points) interpolation (1/r kernel). It then compares the
+ * results with a direct computation.
+ */
+#include <string> 
+#include "Kernels/Chebyshev/FChebCell.hpp"
+#include "Kernels/Chebyshev/FChebSymKernel.hpp"
+//
+template<typename FReal, int ORDER> 
+using FInterpolationCell =  FChebCell<FReal, ORDER>;
+
+template<typename FReal, typename GroupCellClass,
+	 typename GroupContainerClass,
+	 typename MatrixKernelClass, int ORDER>  
+using FInterpolationKernel = FChebSymKernel<FReal,
+					    GroupCellClass,
+					    GroupContainerClass,
+					    MatrixKernelClass,
+					    ORDER> ;
+const std::string interpolationType("Chebyshev interpolation");
+
+#include "sharedMemoryInterpolatiomCmpAlgos.hpp" 

Utils/src/UniformInterpolationCmpAlgo.cpp

+// ==== CMAKE =====
+// @FUSE_BLAS
+// @FUSE_FFT
+// ================
+// Keep in private GIT
+//
+//
+/** \brief Uniform FMM example
+ *
+ * \file
+ * \authors B. Bramas, O. Coulaud
+ *
+ * This program runs the FMM Algorithm with the interpolation kernel based on
+ * uniform (grid points) interpolation (1/r kernel). It then compares the
+ * results with a direct computation.
+ */
+#include <string> 
+#include "Kernels/Uniform/FUnifCell.hpp"
+#include "Kernels/Uniform/FUnifKernel.hpp"
+//
+template<typename FReal, int ORDER> 
+using FInterpolationCell =  FUnifCell<FReal, ORDER>;
+
+template<typename FReal, typename GroupCellClass,
+	 typename GroupContainerClass,
+	 typename MatrixKernelClass, int ORDER>  
+using FInterpolationKernel = FUnifKernel<FReal,
+					    GroupCellClass,
+					    GroupContainerClass,
+					    MatrixKernelClass,
+					    ORDER> ;
+const std::string interpolationType("Uniform interpolation");
+
+#include "sharedMemoryInterpolatiomCmpAlgos.hpp" 

Utils/src/sharedMemoryInterpolatiomCmpAlgos.hpp

+// -*-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"
+//
+#ifdef _OPENMP
+#include "Core/FFmmAlgorithmTask.hpp"
+#include "Core/FFmmAlgorithmNewTask.hpp"
+#include "Core/FFmmAlgorithmSectionTask.hpp"
+#else
+#include "Core/FFmmAlgorithm.hpp"
+#endif
+//
+// Order of the Interpolation approximation
+static constexpr unsigned ORDER = 6 ;
+using FReal                 = double;
+//   1/r kernel
+//
+using MatrixKernelType     = FInterpMatrixKernelR<FReal> ;
+
+//
+/// definition of the common tree structure
+using CellType      = FInterpolationCell<FReal, ORDER>;
+//using CellUpType    = typename CellType::multipole_t;
+//using CellDownType  = typename CellType::local_expansion_t;
+//using CellSymbType  = FSymbolicData;  
+using ContainerType = FP2PParticleContainerIndexed<FReal>;
+using LeafType      = FSimpleLeaf<FReal,  ContainerType >   ;
+using OctreeType     = FOctree<FReal, CellType,ContainerType,LeafType>;
+using KernelType    = FInterpolationKernel<FReal,CellType,ContainerType,MatrixKernelType,ORDER> ;
+
+
+#ifdef _OPENMP
+using TaskFmmAlgo        = FFmmAlgorithmTask<OctreeType,CellType,ContainerType,KernelType,LeafType>;
+using TaskNewFmmAlgo     = FFmmAlgorithmNewTask<OctreeType,CellType,ContainerType,KernelType,LeafType>;
+using SectionTaskFmmAlgo = FFmmAlgorithmSectionTask<OctreeType,CellType,ContainerType,KernelType,LeafType> ;
+#else
+using FmmType  = FFmmAlgorithm<OctreeType,CellType,ContainerType,KernelType,LeafType>;
+#endif
+
+
+
+// Simply create particles and try the kernels
+int main(int argc, char* argv[]) {
+  const FParameterNames LocalOptionAlgo= { {"-algo"} , " Algorithm to run (task, newtask, sectiontask)\n"};
+  const FParameterNames LocalOptionCmp = {
+    {"-cmp"} , "Use to check the result with the exact solution given in the input file\n" };
+  FHelpDescribeAndExit(
+		       argc, argv,
+		       "Driver for Lagrange Or Chebychev interpolation kernel  (1/r kernel).",
+		       FParameterDefinitions::OctreeHeight,
+		       FParameterDefinitions::OctreeSubHeight,
+		       FParameterDefinitions::InputFile,
+		       FParameterDefinitions::OutputFile,
+		       FParameterDefinitions::NbThreads,
+		       LocalOptionAlgo,LocalOptionCmp
+		       );
+  
+
+  const std::string defaultFile(SCALFMMDataPath+"unitCubeXYZQ100.bfma" );
+  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);
+  
+#ifdef _OPENMP
+  const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, omp_get_max_threads());
+  omp_set_num_threads(NbThreads);
+  std::cout << "\n>> Using " << NbThreads << " threads.\n" << std::endl;
+#else
+  const int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
+  std::cout << "\n>> Sequential version.\n" << std::endl;
+#endif
+  //
+  {
+    std::string indent("    ");
+    auto w = std::setw(18);
+    std::cout << "Parameters" << std::endl << std::left
+	      << indent << w << "Octree Depth: "     << TreeHeight    << std::endl
+	      << indent << w << "SubOctree depth: "  << SubTreeHeight << std::endl
+	      << indent << w << "Input file  name: " << filename      << std::endl
+	      << indent << w << "Thread number: "    << NbThreads     << std::endl
+	      << std::endl;
+  }
+  //
+  // init timer
+  FTic time;
+  
+  // open particle file
+  ////////////////////////////////////////////////////////////////////
+  //
+  FFmaGenericLoader<FReal> loader(filename);
+  if (loader.getNbRecordPerline() !=8 ){
+    std::cerr << "File should contain 8 data to read (x,y,z,q,p,fx,fy,fz)\n";
+    std::exit(EXIT_FAILURE);
+  }
+  FSize nbParticles = loader.getNumberOfParticles() ;
+  FmaRWParticle<FReal,8,8> * const particles = new FmaRWParticle<FReal, 8,8>[nbParticles];
+  //
+  //
+  // open particle file
+  ////////////////////////////////////////////////////////////////////
+  //
+  loader.fillParticle(particles,nbParticles);
+  time.tac();
+  //
+  // init oct-tree
+  OctreeType tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
+  FReal  energyD = 0.0 ;
+  //
+  // 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();
+    //
+    //
+    //
+    /////////////////////////////////////////////////////////////////////////////////////////////////
+    // Compute direct energy
+    /////////////////////////////////////////////////////////////////////////////////////////////////
+    
+    for(int idx = 0 ; idx <  nbParticles  ; ++idx){
+      tree.insert(particles[idx].getPosition() , idx, particles[idx].getPhysicalValue() );
+      energyD +=  particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
+    } 
+    time.tac();
+    std::cout << "Done  " << "(@Creating and Inserting Particles = "
+	      << time.elapsed() << " s) ." << std::endl;
+  }
+  ////////////////////////////////////////////////////////////////////
+  //
+  //    Execute FMM Algorithm
+  //
+  ////////////////////////////////////////////////////////////////////
+
+  { // -----------------------------------------------------
+    std::cout << "\n" << interpolationType << "  FMM (ORDER= "<< ORDER << ") ... " << std::endl;
+    
+    const MatrixKernelType    MatrixKernel;
+    time.tic();
+    //
+    std::unique_ptr<KernelType> kernels(new KernelType(TreeHeight, loader.getBoxWidth(), 
+							 loader.getCenterOfBox(),&MatrixKernel));
+    std::string  algoStr  = FParameters::getStr(argc,argv,"-algo",  "task");
+    //
+    TaskFmmAlgo         algo1(&tree, kernels.get() );
+    TaskNewFmmAlgo      algo2(&tree, kernels.get() );
+    SectionTaskFmmAlgo  algo3(&tree, kernels.get() );
+    
+    
+    FAbstractAlgorithm * algo  = nullptr;
+    FAlgorithmTimers   * timer = nullptr;
+    if( "task" == algoStr) {
+      algo  = &algo1 ;
+      timer = &algo1 ;
+    } else if( "newtask" == algoStr) {
+      algo  = &algo2 ;
+      timer = &algo2 ;
+    } else if ( "sectiontask" == algoStr ) {
+      algo  = &algo3 ;
+      timer = &algo3 ;
+    } else {
+      std::cout << "Unknown algorithm: " << algoStr << std::endl;
+      std::exit(EXIT_FAILURE);
+    }
+    std::cout << "Algorithm to check: "<< algoStr <<std::endl;
+    time.tic();
+    
+    //
+    algo->execute();   // Here the call of the FMM algorithm
+    //
+    time.tac();
+    std::cout << "Timers Far Field \n"
+	      << "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"
+	      << "L2P " << timer->getTime(FAlgorithmTimers::L2PTimer) << " seconds\n"
+	      << "P2P and L2P " << timer->getTime(FAlgorithmTimers::NearTimer) << " seconds\n"
+	      << std::endl;
+    
+    
+    std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
+  }
+  // -----------------------------------------------------
+  //
+  // 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;
+    std::cout.precision(10) ;
+
+    tree.forEachLeaf([&](LeafType* 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] ;
+	}
+      });
+    std::cout <<std::endl<<"Energy: "<< energy<<std::endl;
+    std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
+
+  }
+  // -----------------------------------------------------
+  if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
+    std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options,   "output.fma"));
+    FFmaGenericWriter<FReal> writer(name) ;
+    //
+    FSize NbPoints = loader.getNumberOfParticles();
+    FReal * computedParticles = new FReal[8*NbPoints]() ;
+    memset(computedParticles,0,8*NbPoints*sizeof(FReal));
+    FSize j = 0 ;
+    tree.forEachLeaf([&](LeafType* leaf){
+	//
+	// Input
+	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 physicalValues = leaf->getTargets()->getPhysicalValues();
+	const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
+	//
+	// Computed data
+	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();
+	for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
+	  j = 8*indexes[idxPart];
+	  computedParticles[j]      = posX[idxPart] ;
+	  computedParticles[j+1]  = posY[idxPart] ;
+	  computedParticles[j+2]  = posZ[idxPart] ;
+	  computedParticles[j+3]  = physicalValues[idxPart] ;
+	  computedParticles[j+4]  = potentials[idxPart] ;
+	  computedParticles[j+5]  =  forcesX[idxPart] ;
+	  computedParticles[j+6]  =  forcesY[idxPart] ;
+	  computedParticles[j+7]  =  forcesZ[idxPart] ;
+	}
+      });
+
+    writer.writeHeader( loader.getCenterOfBox(), loader.getBoxWidth() ,  NbPoints, sizeof(FReal), 8) ;
+    writer.writeArrayOfReal(computedParticles,  8 , NbPoints);
+
+    delete[] computedParticles;
+
+    //
+    std::string name1( "output.fma");
+    //
+    FFmaGenericWriter<FReal> writer1(name1) ;
+    writer1.writeDistributionOfParticlesFromOctree(&tree,NbPoints) ;
+  }
+  if(FParameters::existParameter(argc, argv, "-cmp") ){
+    std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
+    std::cout << std::scientific;
+    std::cout.precision(10) ;
+
+    printf("Compute Diff...");
+    FMath::FAccurater<FReal> potentialDiff;
+    FMath::FAccurater<FReal> fx, fy, fz, f;
+    FReal energy = 0.0;
+    { // Check that each particle has been summed with all other
+      
+      tree.forEachLeaf([&](LeafType* leaf){
+	  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];
+	    potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
+	    fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
+	    fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
+	    fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
+	    f.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
+	    f.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
+	    f.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
+	    energy   += potentials[idxPart]*physicalValues[idxPart];
+	  }
+	});
+
+      std::cout << energy << " " << energyD << std::endl;
+      delete[] particles;
+
+      f.setNbElements(nbParticles);
+      std::cout << "FChebSymKernel Energy "  << FMath::Abs(energy-energyD) <<  "  Relative     "<< FMath::Abs(energy-energyD) / FMath::Abs(energyD) <<std::endl;
+      std::cout << " Potential Error " << potentialDiff << std::endl;
+      std::cout << " Fx Error " << fx << std::endl;
+      std::cout << " Fy Error " << fy << std::endl;
+      std::cout << " Fz Error " << fz << std::endl;
+      std::cout << " F  Error " << f << std::endl;
+      std::cout <<std::endl<<"Energy: "<< energy<<std::endl;
+      std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
+    }
+    
+    // -----------------------------------------------------
+    
+  }
+
+
+  return 0;
+}