New data files

Tests/noDist/statAdaptiveOctree.cpp

+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas
+// 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".
+// ===================================================================================
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+#include <fstream>
+#include <ctime>
+
+#include "Utils/FParameters.hpp"
+#include "Containers/FOctree.hpp"
+
+#include "adaptiveTree/FAdaptCell.hpp"
+#include "adaptiveTree/FAdaptTools.hpp"
+
+#include "Components/FSimpleIndexedLeaf.hpp"
+#include "Components/FBasicParticleContainer.hpp"
+
+
+#include "Utils/FMath.hpp"
+#include "Files/FFmaGenericLoader.hpp"
+
+
+/// \file  statisticsOnOctree.cpp
+//!
+//! \brief Driver to generate statistics on the octree for the particle distribution given by parameter -infile
+//!  \authors B. Bramas, O. Coulaud
+//!
+//!  This code gives you some statistics (Particles, Cells, P2P and M2L number of operators) on the octree
+//!    Those statistics are shown level by level.
+//!
+//! The statistics are
+//!
+//!  <b>for particles:</b> min/max particles per leaf, the average number, the variance and the average number of P2P neighbors.
+//!
+//! <b>For each level in the octree</b>
+//!
+//!  \arg       The number of cells, of adaptive cells (cell with more one child)
+//!  \arg          The average, min and max numbers of M2L operators and also its  variance.
+
+//!
+//!  <b> General arguments:</b>
+//!     \param   -help(-h)      to see the parameters available in this driver
+//!     \param   -depth    The depth of the octree
+//!     \param   -subdepth          Specifies the size of the sub octree
+//!
+//!     \param   -infile name   Name of the particles file. The file have to be in our FMA format
+//!     \param    -bin              if the input file in binary mode
+//!     \param   -outfile name Generic name  for output file  (without extension)
+//!
+//!  <b> Statistics options:</b>
+//!   \param -stat to build the statistivs on the octree
+//!   \param -histP build a file to generate histogram of particles per leaf. The data are store in file given by -outfile arguments and  .txt extension. (only if -stat is set)
+//   \param   -sM    s_min^M threshold for Multipole expansion (l+1)^2 for Spherical harmonics"
+//   \param   -sL    s_min^M threshold for local expansion  (l+1)^2 for Spherical harmonics"
+// Simply create particles and try the kernels
+//
+void usage() {
+	std::cout << "Driver to obtain statistics on the octree" << std::endl;
+	std::cout <<	 "Options  "<< std::endl
+			<<     "      -help       to see the parameters    " << std::endl
+			<<	     "      -depth        the depth of the octree   "<< std::endl
+			<<	     "      -subdepth   specifies the size of the sub octree   " << std::endl
+			<<     "      -infile name specifies the name of the particle distribution" << std::endl
+            <<    "        -bin              if the input file in binary mode"<< std::endl
+			<<     "      -outfile name  specifies the file for the diagnostics" << std::endl
+			<<     "      -histP   build the histogram of the particle number per leaf"<<std::endl
+			<<     "      -sM    s_min^M threshold for Multipole (l+1)^2 for Spherical harmonics"<<std::endl;
+}
+
+int main(int argc, char ** argv){
+	typedef FBasicParticleContainer<0>                                     ContainerClass;
+	typedef FSimpleIndexedLeaf<ContainerClass>                                  LeafClass;
+	typedef FAdaptCell<FBasicCell,LeafClass>                            CellClass;
+	typedef FOctree<CellClass, ContainerClass, LeafClass >       OctreeClass;
+	//
+	if(FParameters::existParameter(argc, argv, "-h")||FParameters::existParameter(argc, argv, "-help")|| (argc < 3 )){
+		usage() ;
+		exit(-1);
+	}
+	//
+	//   Octree parameters
+	//
+	const int NbLevels        = FParameters::getValue(argc,argv,"-depth", 5);
+	const int SizeSubLevels = FParameters::getValue(argc,argv,"subdepth", 3);
+	const int sminM            = FParameters::getValue(argc,argv,"-sM", 0);
+	const int sminL             = FParameters::getValue(argc,argv,"-sL", 0);
+	//
+	//  input and output  Files parameters
+	//
+	const char* const filename = FParameters::getStr(argc,argv,"-infile", "../Data/test20k.fma");
+	const std::string genericFileName(FParameters::getStr(argc,argv,"-outfile",   "output"));
+	//
+	std::cout << "Opening : " << filename << "\n";
+	bool binaryMode = false;
+	if(FParameters::existParameter(argc, argv, "-bin")){
+		binaryMode = true;
+	}
+	FFmaGenericLoader loader(filename,binaryMode);
+	if(!loader.isOpen()){
+		std::cout << "Loader Error, " << filename << " is missing\n";
+		return 1;
+	}
+	//
+	// -----------------------------------------------------
+	OctreeClass tree(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
+	//
+	// -----------------------------------------------------
+	//     Creating and Inserting particles in the tree
+	// -----------------------------------------------------
+	//
+	std::cout << "Tree box is cubic "<<std::endl
+			<< "         Centre:   "<< loader.getCenterOfBox() <<std::endl
+			<< "         Length:  "<< loader.getBoxWidth()       <<std::endl <<std::endl;
+	//
+	std::cout << "Creating and Inserting " << loader.getNumberOfParticles() << " particles ..." << std::endl;
+	std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
+	FPoint particlePosition, minPos, maxPos;
+	FReal physicalValue;
+	for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
+		loader.fillParticle(&particlePosition,&physicalValue);
+		//
+		minPos.setX(FMath::Min(minPos.getX(),particlePosition.getX())) ;
+		minPos.setY(FMath::Min(minPos.getY(),particlePosition.getY())) ;
+		minPos.setZ(FMath::Min(minPos.getZ(),particlePosition.getZ())) ;
+		maxPos.setX(FMath::Max(maxPos.getX(),particlePosition.getX())) ;
+		maxPos.setX(FMath::Max(maxPos.getY(),particlePosition.getY())) ;
+		maxPos.setX(FMath::Max(maxPos.getZ(),particlePosition.getZ())) ;
+		//
+		tree.insert(particlePosition );
+	}
+	std::cout << "Data are inside the box delimited by "<<std::endl
+			<< "         Min corner:  "<< minPos<<std::endl
+			<< "         Max corner:  "<< maxPos<<std::endl <<std::endl;
+	//
+	OctreeClass::Iterator octreeIterator(&tree);
+
+	//
+	// -----------------------------------------------------
+	//     Build information for adaptive tree
+	// -----------------------------------------------------
+	//
+	{
+		// Set global If for debug purpose
+		long int idCell  = setGlobalID(tree);
+
+		std::cout << " start build smin criteria " <<std::endl;
+		//
+		adaptiveTreeBuilSminC(tree,sminM,sminL) ;
+		//
+		//  Set Global id
+		//
+//		long int idCell  = setGlobalID(tree);
+		//
+		//  Build CA and FA  lists
+		std::cout << " start building CA and FA lists " <<std::endl;
+		//
+		adaptiveTreeBuildLists(tree) ;
+		//
+	}
+	//
+	// -----------------------------------------------------
+	//     Start statistics
+	// -----------------------------------------------------
+	//
+	int removeM = 0 ;
+	if(FParameters::existParameter(argc, argv, "-stat")){ // get stats
+		{    // get stats on the leaf level (Particles)
+			long int allLeaves =  (1 << (3* (NbLevels-1) )) ;
+			std::cout << std::endl<< "[STAT] Leaf level "  << " is  " << NbLevels << std::endl;
+			std::cout << "[STAT] potentials leafs number is " << allLeaves<< std::endl;
+
+			FReal averageParticles = 0.0, varianceParticles = 0.0 ;
+			int nbLeafs = 0,minParticles = 1000000.0, maxParticles = 0.0 ;
+			//
+			// Start to compute statistics on particles
+			//
+			int  nbPart ,nbTPart=0;
+			octreeIterator.gotoBottomLeft();
+			do{
+				nbPart                   = octreeIterator.getCurrentListTargets()->getNbParticles() ;
+				minParticles          =  FMath::Min(minParticles,nbPart) ;
+				maxParticles         =  FMath::Max(maxParticles,nbPart) ;
+				nbTPart              += nbPart;
+				varianceParticles += FReal(nbPart*nbPart) ;
+				++nbLeafs;
+				if(nbPart < sminM){
+					++removeM;
+				}
+			} while(octreeIterator.moveRight());
+			averageParticles   = nbTPart/FReal(nbLeafs);
+			varianceParticles  = varianceParticles/FReal(nbLeafs) - averageParticles*averageParticles;
+			//
+			std::cout.precision(4);
+			std::cout << "[STAT] Non empty leafs: " << nbLeafs << " % of non empty leaves: "<<100*static_cast<FReal>(nbLeafs)/static_cast<FReal>(allLeaves)<<" %" << std::endl;
+			std::cout << "[STAT]  Particles on leafs:"  << std::endl
+					<<   "[STAT]           Min:         "<< minParticles << std::endl
+					<<   "[STAT]           Max:        "<< maxParticles << std::endl
+					<<   "[STAT]           Average:  "<< averageParticles << std::endl
+					<<   "[STAT]           Variance: " << varianceParticles << std::endl;
+			std::cout << "[STAT]  number of P2M to remove: " << 		removeM <<	 std::endl;
+			//
+			//  Histogram of particles per leaf
+			//
+
+			if(FParameters::existParameter(argc, argv, "-histP")){
+				int size = maxParticles+1;
+				int * hist = new int [size] ;
+				memset(hist,0,(size)*sizeof(int));
+				octreeIterator.gotoBottomLeft();
+				do{
+					nbPart  = octreeIterator.getCurrentListTargets()->getNbParticles() ;
+					++hist[nbPart] ;
+				} while(octreeIterator.moveRight());
+				//
+				// write data
+				//
+				std::ofstream outfile( genericFileName + ".txt", std::ofstream::out);
+				if(!outfile) {
+					std::cout << "Cannot open file "<< std::endl;
+					exit(-1)	;
+				}	//
+				outfile << "# Particle histogram.   "<< size << " chunk" <<std::endl;
+				for(int i=0 ; i < size ; ++i){
+					outfile << i << "  " << hist[i] <<std::endl;
+				}
+				delete [] hist ;
+			}
+			FReal averageNeighbors = 0.0, varianceNeighbors =0.0 ;
+			int nbBox,minBox=30,maxBox=0;
+			octreeIterator.gotoBottomLeft();
+			ContainerClass*  neighborsP2P[27];
+			do{
+				//
+				//  P2P Neighbors
+				//
+				nbBox = tree.getLeafsNeighbors(neighborsP2P, octreeIterator.getCurrentGlobalCoordinate(),NbLevels-1) ;
+				// need the current particles and neighbors particles
+				minBox                      =  FMath::Min(minBox,nbBox) ;
+				maxBox                     =  FMath::Max(maxBox,nbBox) ;
+				averageNeighbors   += FReal(nbBox);
+				varianceNeighbors  += FReal(nbBox*nbBox) ;
+			} while(octreeIterator.moveRight());
+			//
+			averageNeighbors/=FReal(nbLeafs) ;
+			varianceNeighbors = varianceNeighbors/nbLeafs-averageNeighbors*averageNeighbors;
+			//
+			std::cout << "[STAT]  P2P Neighbors for each leaf " << std::endl
+					<< "[STAT]           Min:       " <<  minBox << std::endl
+					<< "[STAT]           Max:      " <<  maxBox << std::endl
+					<< "[STAT]           Average: " <<  averageNeighbors<< std::endl
+					<< "[STAT]           Variance: " <<  varianceNeighbors << std::endl<< std::endl;
+		}
+		//
+		//    ---------------  END LEAVES ---------------
+		//
+		{
+
+			long long int totalCells = 0;
+			long long int totalM2L = 0;
+			long long int totalM2ML2L = 0;
+
+			int nbCellsAtTop = 0;
+			int nbCellsAtBottom = 0;
+
+			octreeIterator.gotoBottomLeft();
+
+			for(int idxLevel = NbLevels - 1 ; idxLevel >= 1 ; --idxLevel){
+				removeM =0 ;
+				int nbCellsAtLevel = 0;
+				int nbChildAtLevel = 0, adaptiveCell=0 ,nbChildForMyCell;
+				int nbNeighborsAtLevel = 0;
+				//
+				int nbM2LNeighbors, minM2L=500,maxM2L=-1;
+				FReal averageM2LNeighbors=0.0, varianceM2LNeighbors=0.0	;
+				//
+				const CellClass* neighborsM2L[343];
+				do{
+					++nbCellsAtLevel;
+					// Check number of
+					if( idxLevel != NbLevels - 1 ){
+						nbChildForMyCell=0 ;
+						auto** child = octreeIterator.getCurrentChild();
+						auto  &  myCell = *(octreeIterator.getCurrentCell());
+						int nbPart = 0 ;
+						//						std::cout << "NB: ";
+						for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
+							if(child[idxChild]) {
+								++nbChildForMyCell;
+								nbPart += child[idxChild]->getnbPart();
+								//								std::cout << "  "<< child[idxChild]->getnbPart();
+							}
+						}
+						//						std::cout << std::endl;
+						octreeIterator.getCurrentCell()->addPart(nbPart);
+						if(octreeIterator.getCurrentCell()->getnbPart() < sminM){
+							++removeM;
+						}
+						nbChildAtLevel += nbChildForMyCell ;
+						if(nbChildForMyCell>1) {
+							++adaptiveCell ;
+						}
+						else
+						{myCell.setCelladaptive();}
+					}
+					const CellClass* neighbors[343];
+					nbNeighborsAtLevel += tree.getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(),idxLevel);
+					//
+					//  M2L Neighbors
+					//
+					nbM2LNeighbors = tree.getInteractionNeighbors(neighborsM2L, octreeIterator.getCurrentGlobalCoordinate(),idxLevel);
+					minM2L                          =  FMath::Min(minM2L,nbM2LNeighbors) ;
+					maxM2L                         =  FMath::Max(maxM2L,nbM2LNeighbors) ;
+					averageM2LNeighbors  += FReal(nbM2LNeighbors) ;
+					varianceM2LNeighbors += FReal(nbM2LNeighbors*nbM2LNeighbors) ;
+				} while(octreeIterator.moveRight());
+				//
+				averageM2LNeighbors/=FReal(nbCellsAtLevel) ;
+				varianceM2LNeighbors = varianceM2LNeighbors/nbCellsAtLevel-averageM2LNeighbors*averageM2LNeighbors;
+
+				std::cout << "[STAT] Level = " << idxLevel << std::endl
+						<< "[STAT]     >> Nb Cells =                                 \t " << nbCellsAtLevel << std::endl
+						<< "[STAT]     >> Nb Adaptive Cells =                   \t" << adaptiveCell    << "  Non Adaptive (1 son): " <<100*FReal(nbCellsAtLevel-adaptiveCell)/nbCellsAtLevel<<std::endl
+						<< "[STAT]     >> Number of M2M to remove: " << 		removeM <<	 std::endl
+						<< "[STAT]     >> Nb M2M/L2L interactions =        \t" << nbChildAtLevel << std::endl
+						<< "[STAT]     >> Average M2M/L2L interactions = \t" << FReal(nbChildAtLevel)/FReal(nbCellsAtLevel) << std::endl
+						<< "[STAT]     >> Nb M2L interactions =                 \t" << nbNeighborsAtLevel << std::endl;
+				std::cout << "[STAT]     >> M2L Neighbors for each leaf " << std::endl
+						<< "[STAT]             >>  Min:        " <<  minM2L << std::endl
+						<<  "[STAT]             >> Max:       " <<  maxM2L << std::endl
+						<<  "[STAT]             >> Average: " <<  averageM2LNeighbors<< std::endl
+						<<  "[STAT]             >> Variance: " <<  varianceM2LNeighbors << std::endl<< std::endl;
+
+				totalCells += (long long int)(nbCellsAtLevel);
+				totalM2L += (long long int)(nbNeighborsAtLevel);
+				totalM2ML2L += (long long int)(nbChildAtLevel);
+				nbCellsAtTop = nbCellsAtLevel;
+				if( idxLevel == NbLevels - 1 ) nbCellsAtBottom = nbCellsAtLevel;
+				std::cout << std::endl;
+				//
+				//  Go to next level
+				octreeIterator.moveUp();
+				octreeIterator.gotoLeft();
+				//
+			}
+			//
+			// Global statistics on the octree
+			//
+			std::cout << "[STAT] For all the tree\n";
+			std::cout << "[STAT] >> Total Nb Cells = " << totalCells-nbCellsAtTop << "\n";
+			std::cout << "[STAT] >> Total Nb M2M/L2L interactions = " << totalM2ML2L << "\n";
+			std::cout << "[STAT] >> Total Average M2M/L2L interactions = " << FReal(totalM2ML2L-nbCellsAtTop)/FReal(totalCells-nbCellsAtBottom) << "\n";
+			std::cout << "[STAT] >> Total Nb M2L interactions per cell = " << totalM2L << "\n";
+			std::cout << "[STAT] >> Total Average M2L interactions per cell = " << FReal(totalM2L)/FReal(totalCells) << "\n";
+
+			std::cout << "nbCellsAtTop " << nbCellsAtTop <<std::endl;
+			//		idCell = totalCells ;
+		}
+
+	}
+	//
+	//  Set Global id for tulip export
+	//
+	long int idCell  = setGlobalID(tree);
+	//
+	//
+	std::cout << " start export tulip " <<std::endl;
+
+	//
+	// Set Global indexes to save the octree in tulip format
+	//
+	// -----------------------------------------------------
+	std::ofstream tlp("aa.tlp", std::ofstream::out );
+
+	TulipExport( tlp, idCell, tree);
+
+	std::cout << " NVCells " << idCell <<  std::endl ;
+
+
+	octreeIterator.gotoTop() ;
+	for(int idxLevel = 1 ; idxLevel < NbLevels ;  ++idxLevel){
+		std::cout << "idxLevel: "<<idxLevel << "  Iterator Level    " << octreeIterator.level()<<  "  is leaves level: " << octreeIterator.isAtLeafLevel()	<<std::endl;
+		octreeIterator.moveDown() ;octreeIterator.gotoLeft();
+	}
+	std::cout << "Level max " <<  NbLevels <<std::endl;
+
+
+
+	std::ofstream file("aa.tree", std::ofstream::out );
+	//
+	octreeIterator.gotoTop() ;  // here we are at level 1 (first child)
+	//
+	////////////////////////////////////////////////////////////////////
+	//              Export adaptive tree in our format
+	////////////////////////////////////////////////////////////////////
+	//
+	// -----------------------------------------------------
+	//
+	std::cout << "Top of the octree " << octreeIterator.level() << std::endl ;
+	for(int idxLevel = 1 ; idxLevel < NbLevels ;  ++idxLevel){
+		file << std::endl << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<< std::endl;
+		file << "  Level " << idxLevel <<" OLevel  "<<  octreeIterator.level()<<  "  -- leave level " <<   std::boolalpha <<  octreeIterator.isAtLeafLevel() << std::endl;
+		do{
+			file << *(octreeIterator.getCurrentCell())<< std::endl ;
+		} while(octreeIterator.moveRight());
+		octreeIterator.moveDown() ;
+		octreeIterator.gotoLeft();
+	}
+	std::cout << "   END    " << std::endl;
+
+	// Check
+	octreeIterator.gotoBottomLeft();
+	do {
+		std::cout << " Level " <<octreeIterator.level() <<std::endl;
+	}while(octreeIterator.moveUp() );
+	//
+	return 0;
+}
+
+
+

UTests/utestChebyshevThread.cpp

+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
+// 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".
+// ===================================================================================
+
+// ==== CMAKE =====
+// @FUSE_BLAS
+// ================
+#include "ScalFmmConfig.h"
+#include "Utils/FGlobal.hpp"
+
+#include "Containers/FOctree.hpp"
+
+#include "Files/FFmaGenericLoader.hpp"
+
+#include "Core/FFmmAlgorithmThread.hpp"
+
+#include "FUTester.hpp"
+
+#include "Components/FSimpleLeaf.hpp"
+
+
+#include "Kernels/Chebyshev/FChebCell.hpp"
+#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
+#include "Kernels/Chebyshev/FChebKernel.hpp"
+#include "Kernels/Chebyshev/FChebSymKernel.hpp"
+
+#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
+/*
+  In this test we compare the Chebyschev fmm results and the direct results.
+ */
+
+
+/** the test class
+ *
+ */
+class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
+
+	///////////////////////////////////////////////////////////
+	// The tests!
+	///////////////////////////////////////////////////////////
+
+	template <class CellClass, class ContainerClass, class KernelClass, class MatrixKernelClass,
+	class LeafClass, class OctreeClass, class FmmClass>
+	void RunTest()	{
+		//
+#ifdef _OPENMP
+	std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
+#else
+	std::cout << "\n>> OpenMP test !!!\n" << std::endl;
+	exit(EXIT_FAILURE);
+#endif
+		//
+		// Load particles
+		//
+		if(sizeof(FReal) == sizeof(float) ) {
+			std::cerr << "No input data available for Float "<< std::endl;
+			exit(EXIT_FAILURE);
+		}
+		const std::string parFile( (sizeof(FReal) == sizeof(float))?
+				"Test/DirectFloatbfma":
+				"UTest/DirectDouble.bfma");
+		//
+		std::string filename(SCALFMMDataPath+parFile);
+		//
+		FFmaGenericLoader loader(filename);
+		Print("Number of particles:");
+		Print(loader.getNumberOfParticles());
+
+		const int NbLevels        = 4;
+		const int SizeSubLevels = 2;
+
+		// Create octree
+
+		FSize nbParticles = loader.getNumberOfParticles() ;
+		FmaR8W8Particle* const particles = new FmaR8W8Particle[nbParticles];
+
+		loader.fillParticle(particles,nbParticles);
+         //
+		//   Insert particle in the tree
+		//
+		OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+		for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
+			tree.insert(particles[idxPart].position , idxPart, particles[idxPart].physicalValue );
+		}
+		//
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+		// Run FMM computation
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+		Print("Fmm...");
+		KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+		FmmClass algo(&tree,&kernels);
+		algo.execute();
+		//0
+		FReal energy= 0.0 , energyD = 0.0 ;
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+		// Compute direct energy
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+
+		for(int idx = 0 ; idx <  loader.getNumberOfParticles()  ; ++idx){
+			energyD +=  particles[idx].potential*particles[idx].physicalValue ;
+		}
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+		// Compare
+		/////////////////////////////////////////////////////////////////////////////////////////////////
+		Print("Compute Diff...");
+		FMath::FAccurater potentialDiff;
+		FMath::FAccurater fx, fy, fz;
+		{ // Check that each particle has been summed with all other
+
+			tree.forEachLeaf([&](LeafClass* leaf){
+				const FReal*const potentials        = leaf->getTargets()->getPotentials();
+				const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
+				const FReal*const forcesX            = leaf->getTargets()->getForcesX();
+				const FReal*const forcesY            = leaf->getTargets()->getForcesY();
+				const FReal*const forcesZ            = leaf->getTargets()->getForcesZ();
+				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];
+					potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
+					fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
+					fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
+					fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
+					energy   += potentials[idxPart]*physicalValues[idxPart];
+				}
+			});
+		}
+
+		delete[] particles;
+
+		// Print for information
+
+		Print("Potential diff is = ");
+		printf("         Pot L2Norm     %e\n",potentialDiff.getL2Norm());
+		printf("         Pot RL2Norm   %e\n",potentialDiff.getRelativeL2Norm());
+		printf("         Pot RMSError   %e\n",potentialDiff.getRMSError());
+		Print("Fx diff is = ");
+		printf("         Fx L2Norm     %e\n",fx.getL2Norm());
+		printf("         Fx RL2Norm   %e\n",fx.getRelativeL2Norm());
+		printf("         Fx RMSError   %e\n",fx.getRMSError());
+		Print("Fy diff is = ");
+		printf("        Fy L2Norm     %e\n",fy.getL2Norm());
+		printf("        Fy RL2Norm   %e\n",fy.getRelativeL2Norm());
+		printf("        Fy RMSError   %e\n",fy.getRMSError());
+		Print("Fz diff is = ");
+		printf("        Fz L2Norm     %e\n",fz.getL2Norm());
+		printf("        Fz RL2Norm   %e\n",fz.getRelativeL2Norm());
+		printf("        Fz RMSError   %e\n",fz.getRMSError());
+		FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm()  + fz.getRelativeL2Norm() *fz.getRelativeL2Norm()  );
+		printf(" Total L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;