Merge branch 'develop' of gitlab.inria.fr:solverstack/ScalFMM into develop

Examples/MPIInterpolationFMM.hpp

@@ -70,16 +70,19 @@ using FmmClassProcPER  = FFmmAlgorithmThreadProcPeriodic<FReal,OctreeClass,CellC
 int main(int argc, char* argv[])
 {
   ///////// PARAMETERS HANDLING //////////////////////////////////////
-  //
-    FHelpDescribeAndExit(argc, argv,
-                       "Driver for Interpolation Laplace kernel using MPI. "
+  const FParameterNames  localIncreaseBox = { {"ratio","-L"}, "Increase teh Box size by a factor L:= ratio"};
+  FHelpDescribeAndExit(argc, argv,
+                       "Driver for Chebyshev Interpolation kernel using MPI  (1/r kernel).\n "
                        "Usully run using : mpirun -np nb_proc_needed ./ChebyshevInterpolationAlgorithm [params].",
                        FParameterDefinitions::OctreeHeight,
                        FParameterDefinitions::OctreeSubHeight,
                        FParameterDefinitions::InputFile,
                        FParameterDefinitions::OutputFile,
                        FParameterDefinitions::NbThreads,
-                       FParameterDefinitions::PeriodicityNbLevels);
+                       FParameterDefinitions::PeriodicityNbLevels,
+		       localIncreaseBox
+		       )
+;
 
   const std::string defaultFile("../Data/test20k.fma");
   const std::string  filename      = FParameters::getStr(argc,argv,    FParameterDefinitions::InputFile.options, defaultFile.c_str());
@@ -121,16 +124,25 @@ int main(int argc, char* argv[])
   if(!loader.isOpen()) {
       throw std::runtime_error("Particle file couldn't be opened!") ;
     }
-
+  auto boxWidth = loader.getBoxWidth() ;
+  //
+  if(FParameters::existParameter(argc, argv, localIncreaseBox.options)){
+    FReal ratio=  FParameters::getValue(argc, argv, localIncreaseBox.options, 1.0);
+    boxWidth *= ratio;
+  }
+ 
   // Initialize empty oct-tree
-  OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
+  OctreeClass tree(TreeHeight, SubTreeHeight, boxWidth, loader.getCenterOfBox());
 
   FSize localParticlesNumber = 0 ;
 
   // -----------------------------------------------------
     if(app.global().processId() == 0){
         std::cout << "Loading & Inserting " << loader.getNumberOfParticles()
-                  << " particles ..." << std::endl;
+                  << " particles ..." << std::endl
+		  <<" Box: "<< std::endl
+		  << "    width  " << boxWidth << std::endl
+		  << "    Centre " << loader.getCenterOfBox()<< std::endl;
         std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
       }
     time.tic();
@@ -221,7 +233,7 @@ int main(int argc, char* argv[])
 
 
     // non periodic FMM algorithm
-    std::unique_ptr<KernelClass> kernelsNoPer(new KernelClass(TreeHeight, loader.getBoxWidth(),
+    std::unique_ptr<KernelClass> kernelsNoPer(new KernelClass(TreeHeight, boxWidth,
                                                               loader.getCenterOfBox(),&MatrixKernel));
     FmmClassProc    algoNoPer(app.global(),&tree, kernelsNoPer.get());
     //
@@ -240,13 +252,12 @@ int main(int argc, char* argv[])
         //timer      = &algoPer ;
       }
     //
-    // FMM exectution  FFmmFarField
-    algorithm->execute();
+    // FMM exectution  FFmmFarField FFmmNearField
+    algorithm->execute(FFmmNearField);
     //
 
     time.tac();
-    double timeUsed = time.elapsed();
-    double minTime,maxTime;
+    timeUsed = time.elapsed();
     std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << " s)." << std::endl;
     MPI_Reduce(&timeUsed,&minTime,1,MPI_DOUBLE,MPI_MIN,0,app.global().getComm());
     MPI_Reduce(&timeUsed,&maxTime,1,MPI_DOUBLE,MPI_MAX,0,app.global().getComm());
@@ -270,7 +281,9 @@ FReal energy =0.0 ;
 //
 std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
 std::cout << std::scientific;
-std::cout.precision(10) ;
+std::cout.precision(15) ;
+ 
+ FReal TotalPhysicalValue=0.0 ;
 
 tree.forEachLeaf([&](LeafClass* leaf){
   const FReal*const posX = leaf->getTargets()->getPositions()[0];
@@ -294,11 +307,12 @@ tree.forEachLeaf([&](LeafClass* leaf){
                           << "   Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
         }
       energy += potentials[idxPart]*physicalValues[idxPart] ;
+      TotalPhysicalValue += physicalValues[idxPart]  ;
     }
 });
 FReal gloEnergy = app.global().reduceSum(energy);
 if(0 == app.global().processId()){
-  std::cout <<std::endl << "Proc "<< app.global().processId() << " Energy: "<< gloEnergy <<std::endl;
+     std::cout <<std::endl<<"aboveRoot: " << aboveTree << "  Energy: "<< energy<<"  TotalPhysicalValue: " << TotalPhysicalValue<< std::endl; std::cout <<std::endl << "Proc "<< app.global().processId() << " Energy: "<< gloEnergy <<std::endl;
   std::cout <<std::endl <<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
 }
 }