diff --git a/Data/UTest/unitCubeRef20kDouble.bfma b/Data/UTest/unitCubeRef20kDouble.bfma new file mode 100644 index 0000000000000000000000000000000000000000..48ed9ddd4f0e01e6654fdb03b93ebe7d04f1fba5 Binary files /dev/null and b/Data/UTest/unitCubeRef20kDouble.bfma differ diff --git a/Tests/noDist/statAdaptiveOctree.cpp b/Tests/noDist/statAdaptiveOctree.cpp new file mode 100755 index 0000000000000000000000000000000000000000..167dd377ed14cc56ddf036931d267f93055cb805 --- /dev/null +++ b/Tests/noDist/statAdaptiveOctree.cpp @@ -0,0 +1,433 @@ +// =================================================================================== +// 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; +} + + + diff --git a/UTests/utestChebyshevThread.cpp b/UTests/utestChebyshevThread.cpp new file mode 100755 index 0000000000000000000000000000000000000000..6be6e4a889eb3b25b69bcf62e8c9e54efcf6c35c --- /dev/null +++ b/UTests/utestChebyshevThread.cpp @@ -0,0 +1,261 @@ +// =================================================================================== +// 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)) ; + printf(" Energy Error = %.12e\n",FMath::Abs(energy-energyD)); + printf(" Energy FMM = %.12e\n",FMath::Abs(energy)); + printf(" Energy DIRECT = %.12e\n",FMath::Abs(energyD)); + + // Assert + const FReal MaximumDiffPotential = FReal(9e-3); + const FReal MaximumDiffForces = FReal(9e-2); + + Print("Test1 - Error Relative L2 norm Potential "); + uassert(potentialDiff.getRelativeL2Norm() < MaximumDiffPotential); //1 + Print("Test2 - Error RMS L2 norm Potential "); + uassert(potentialDiff.getRMSError() < MaximumDiffPotential); //2 + Print("Test3 - Error Relative L2 norm FX "); + uassert(fx.getRelativeL2Norm() < MaximumDiffForces); //3 + Print("Test4 - Error RMS L2 norm FX "); + uassert(fx.getRMSError() < MaximumDiffForces); //4 + Print("Test5 - Error Relative L2 norm FY "); + uassert(fy.getRelativeL2Norm() < MaximumDiffForces); //5 + Print("Test6 - Error RMS L2 norm FY "); + uassert(fy.getRMSError() < MaximumDiffForces); //6 + Print("Test7 - Error Relative L2 norm FZ "); + uassert(fz.getRelativeL2Norm() < MaximumDiffForces); //8 + Print("Test8 - Error RMS L2 norm FZ "); + uassert(fz.getRMSError() < MaximumDiffForces); //8 + Print("Test9 - Error Relative L2 norm F "); + uassert(L2error < MaximumDiffForces); //9 Total Force + Print("Test10 - Relative error Energy "); + uassert(FMath::Abs(energy-energyD) /energyD< MaximumDiffPotential); //10 Total Energy + + + } + + /** If memstas is running print the memory used */ + void PostTest() { + if( FMemStats::controler.isUsed() ){ + std::cout << "Memory used at the end " << FMemStats::controler.getCurrentAllocated() + << " Bytes (" << FMemStats::controler.getCurrentAllocatedMB() << "MB)\n"; + std::cout << "Max memory used " << FMemStats::controler.getMaxAllocated() + << " Bytes (" << FMemStats::controler.getMaxAllocatedMB() << "MB)\n"; + std::cout << "Total memory used " << FMemStats::controler.getTotalAllocated() + << " Bytes (" << FMemStats::controler.getTotalAllocatedMB() << "MB)\n"; + } + } + + + /////////////////////////////////////////////////////////// + // Set the tests! + /////////////////////////////////////////////////////////// + + + /** TestChebKernel */ + void TestChebKernel(){ + const unsigned int ORDER = 6; + typedef FP2PParticleContainerIndexed<> ContainerClass; + typedef FSimpleLeaf<ContainerClass> LeafClass; + typedef FInterpMatrixKernelR MatrixKernelClass; + typedef FChebCell<ORDER> CellClass; + typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass; + typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass; + typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass; + // run test + RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(); + } + + /** TestChebSymKernel */ + void TestChebSymKernel(){ + const unsigned int ORDER = 6; + typedef FP2PParticleContainerIndexed<> ContainerClass; + typedef FSimpleLeaf<ContainerClass> LeafClass; + typedef FInterpMatrixKernelR MatrixKernelClass; + typedef FChebCell<ORDER> CellClass; + typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass; + typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass; + typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass; + // run test + RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(); + } + + + + /////////////////////////////////////////////////////////// + // Set the tests! + /////////////////////////////////////////////////////////// + + /** set test */ + void SetTests(){ + AddTest(&TestChebyshevDirect::TestChebKernel,"Test Chebyshev Kernel with one big SVD"); + AddTest(&TestChebyshevDirect::TestChebSymKernel,"Test Chebyshev Kernel with 16 small SVDs and symmetries"); + } +}; + + +// You must do this +TestClass(TestChebyshevDirect) + + + + diff --git a/UTests/utestLagrangeThread.cpp b/UTests/utestLagrangeThread.cpp new file mode 100755 index 0000000000000000000000000000000000000000..03736c7dab226a1f526d6998501515af266e307f --- /dev/null +++ b/UTests/utestLagrangeThread.cpp @@ -0,0 +1,246 @@ +// =================================================================================== +// 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_FFT +// ================ + +#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 "../../Src/Kernels/Uniform/FUnifCell.hpp" +#include "../../Src/Kernels/Interpolation/FInterpMatrixKernel.hpp" +#include "../../Src/Kernels/Uniform/FUnifKernel.hpp" + +#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp" +/* + In this test we compare the TestLagrange fmm results with the direct results. + */ + + +/** the test class + * + */ +class TestLagrange : public FUTester<TestLagrange> { + + /////////////////////////////////////////////////////////// + // 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)) ; + printf(" Energy Error = %.12e\n",FMath::Abs(energy-energyD)); + printf(" Energy FMM = %.12e\n",FMath::Abs(energy)); + printf(" Energy DIRECT = %.12e\n",FMath::Abs(energyD)); + + // Assert + const FReal MaximumDiffPotential = FReal(9e-3); + const FReal MaximumDiffForces = FReal(9e-2); + + Print("Test1 - Error Relative L2 norm Potential "); + uassert(potentialDiff.getRelativeL2Norm() < MaximumDiffPotential); //1 + Print("Test2 - Error RMS L2 norm Potential "); + uassert(potentialDiff.getRMSError() < MaximumDiffPotential); //2 + Print("Test3 - Error Relative L2 norm FX "); + uassert(fx.getRelativeL2Norm() < MaximumDiffForces); //3 + Print("Test4 - Error RMS L2 norm FX "); + uassert(fx.getRMSError() < MaximumDiffForces); //4 + Print("Test5 - Error Relative L2 norm FY "); + uassert(fy.getRelativeL2Norm() < MaximumDiffForces); //5 + Print("Test6 - Error RMS L2 norm FY "); + uassert(fy.getRMSError() < MaximumDiffForces); //6 + Print("Test7 - Error Relative L2 norm FZ "); + uassert(fz.getRelativeL2Norm() < MaximumDiffForces); //8 + Print("Test8 - Error RMS L2 norm FZ "); + uassert(fz.getRMSError() < MaximumDiffForces); //8 + Print("Test9 - Error Relative L2 norm F "); + uassert(L2error < MaximumDiffForces); //9 Total Force + Print("Test10 - Relative error Energy "); + uassert(FMath::Abs(energy-energyD) /energyD< MaximumDiffPotential); //10 Total Energy + + + } + + /** If memstas is running print the memory used */ + void PostTest() { + if( FMemStats::controler.isUsed() ){ + std::cout << "Memory used at the end " << FMemStats::controler.getCurrentAllocated() + << " Bytes (" << FMemStats::controler.getCurrentAllocatedMB() << "MB)\n"; + std::cout << "Max memory used " << FMemStats::controler.getMaxAllocated() + << " Bytes (" << FMemStats::controler.getMaxAllocatedMB() << "MB)\n"; + std::cout << "Total memory used " << FMemStats::controler.getTotalAllocated() + << " Bytes (" << FMemStats::controler.getTotalAllocatedMB() << "MB)\n"; + } + } + + + /////////////////////////////////////////////////////////// + // Set the tests! + /////////////////////////////////////////////////////////// + + + /** TestUnifKernel */ + void TestUnifKernel(){ + const unsigned int ORDER = 6; + // typedefs + typedef FP2PParticleContainerIndexed<> ContainerClass; + typedef FSimpleLeaf< ContainerClass > LeafClass; + typedef FInterpMatrixKernelR MatrixKernelClass; + typedef FUnifCell<ORDER> CellClass; + typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass; + typedef FUnifKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass; + typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass; + + // run test + RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(); + } + + /////////////////////////////////////////////////////////// + // Set the tests! + /////////////////////////////////////////////////////////// + + /** set test */ + void SetTests(){ + AddTest(&TestLagrange::TestUnifKernel,"Test Lagrange Kernel "); + } +}; + + +// You must do this +TestClass(TestLagrange) + + + +