diff --git a/Tests/Kernels/testAdaptiveChebSymFMM.cpp b/Tests/Kernels/testAdaptiveChebSymFMM.cpp new file mode 100644 index 0000000000000000000000000000000000000000..c14b26624387faeea662eec68523a78a36a44d39 --- /dev/null +++ b/Tests/Kernels/testAdaptiveChebSymFMM.cpp @@ -0,0 +1,189 @@ +// =================================================================================== +// 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 <iostream> +#include <cstdio> + +#include "Utils/FParameters.hpp" +#include "Utils/FTic.hpp" + +#include "Containers/FOctree.hpp" +#include "Components/FSimpleLeaf.hpp" + +#include "Utils/FPoint.hpp" + +#include "Files/FFmaGenericLoader.hpp" +#include "Files/FRandomLoader.hpp" + +#include "Components/FBasicKernels.hpp" +#include "Components/FSimpleIndexedLeaf.hpp" +#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp" + +#include "Kernels/Interpolation/FInterpMatrixKernel.hpp" +#include "Kernels/Chebyshev/FChebCell.hpp" + +#include "Adaptive/FAdaptiveCell.hpp" +#include "Adaptive/FAdaptiveKernelWrapper.hpp" +#include "Adaptive/FAbstractAdaptiveKernel.hpp" +#include "Adaptive/FAdaptChebSymKernel.hpp" + +#include "Kernels/Interpolation/FInterpMatrixKernel.hpp" +#include "Kernels/Chebyshev/FChebCell.hpp" +#include "Adaptive/FAdaptTools.hpp" + +#include "Adaptive/FAdaptivePrintKernel.hpp" + +#include "Core/FFmmAlgorithm.hpp" +#include "Utils/FParameterNames.hpp" + +/** This program show an example of use of the fmm basic algo + * it also check that each particles is impacted each other particles + */ + + +// Simply create particles and try the kernels +int main(int argc, char ** argv){ + // + const FParameterNames LocalOptionMinMultipoleThreshod {{"-sM"}," s_min^M threshold for Multipole (l+1)^2 for Spherical harmonic."}; + const FParameterNames LocalOptionMinLocalThreshod {{"-SL"}," s_min^L threshold for Local (l+1)^2 for Spherical harmonics."}; + + FHelpDescribeAndExit(argc, argv, + "Test Adaptive kernel and compare it with the direct computation.", + FParameterDefinitions::OctreeHeight,FParameterDefinitions::NbThreads, + FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile, + LocalOptionMinMultipoleThreshod,LocalOptionMinLocalThreshod); + + const unsigned int P = 5 ; + typedef FChebCell<P> CellClass; + typedef FP2PParticleContainerIndexed<> ContainerClass; + typedef FSimpleLeaf<ContainerClass> LeafClass; + typedef FInterpMatrixKernelR MatrixKernelClass; + typedef FAdaptiveChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,P> KernelClass; + typedef FAdaptiveCell< CellClass, ContainerClass > CellWrapperClass; + typedef FAdaptiveKernelWrapper< KernelClass, CellClass, ContainerClass > KernelWrapperClass; + typedef FOctree< CellWrapperClass, ContainerClass , LeafClass > OctreeClass; + typedef FFmmAlgorithm<OctreeClass, CellWrapperClass, ContainerClass, KernelWrapperClass, LeafClass > FmmClass; + + FTic counter; + + ////////////////////////////////////////////////////////////////////////////////// + const int sminM = FParameters::getValue(argc,argv,LocalOptionMinMultipoleThreshod.options, P*P*P); + const int sminL = FParameters::getValue(argc,argv,LocalOptionMinLocalThreshod.options, P*P*P); + const std::string fileName(FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/noDistprolate50.out.fma")); + const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 3); + const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2); + + FFmaGenericLoader loader(fileName); + const long int NbPart = loader.getNumberOfParticles() ; + ////////////////////////////////////////////////////////////////////////////////// + + OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox()); + + ////////////////////////////////////////////////////////////////////////////////// + + std::cout << "Creating & Inserting " << NbPart << " particles ..." << std::endl; + std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl; + std::cout << " criteria SM: "<< sminM <<std::endl + << " criteria SL: "<< sminL <<std::endl <<std::endl; + + counter.tic(); + + FmaRWParticle<8,8>* const particles = new FmaRWParticle<8,8>[NbPart]; + loader.fillParticle(particles,NbPart); + + for(int idxPart = 0 ; idxPart < NbPart; ++idxPart){ + const FPoint PP(particles[idxPart].getPosition() ) ; + tree.insert(PP, idxPart, particles[idxPart].getPhysicalValue()); + } + + counter.tac(); + std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << " s)." << std::endl; + + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + + std::cout << "Working on particles ..." << std::endl; + + counter.tic(); + + const MatrixKernelClass MatrixKernel; + KernelWrapperClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel,sminM,sminL); // FTestKernels FBasicKernels + FmmClass algo(&tree,&kernels); //FFmmAlgorithm FFmmAlgorithmThread + algo.execute(); + + counter.tac(); + std::cout << "Done " << "(@Algorithm = " << counter.elapsed() << " s)." << std::endl; + + ///////////////////////////////////////////////////////////////////////////////////////////////// + // Compute direct energy + ///////////////////////////////////////////////////////////////////////////////////////////////// + FReal energyD = 0.0 ; + for(int idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){ + energyD += particles[idx].getPotential()*particles[idx].getPhysicalValue() ; + } + + ///////////////////////////////////////////////////////////////////////////////////////////////// + // Compare + ///////////////////////////////////////////////////////////////////////////////////////////////// + { + FMath::FAccurater potentialDiff; + FMath::FAccurater fx, fy, fz; + FReal energy= 0.0; + { // Check that each particle has been summed with all other + + // std::cout << "indexPartOrig || DIRECT V fx || FMM V fx" << std::endl; + + tree.forEachCellLeaf([&](CellWrapperClass* cell, 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].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]); + energy += potentials[idxPart]*physicalValues[idxPart]; + + } + }); + } + + // Print for information + std::cout << "Energy [relative L2 error] " << FMath::Abs(energy-energyD) /energyD << std::endl; + std::cout << "Potential " << potentialDiff << std::endl; + std::cout << "Fx " << fx << std::endl; + std::cout << "Fy " << fy << std::endl; + std::cout << "Fz " << fz << std::endl; + } + + delete[] particles; + + return 0; +} + + + + diff --git a/Tests/noDist/testAdaptiveChebSymFMM.cpp b/Tests/noDist/testAdaptiveChebSymFMMDebug.cpp similarity index 100% rename from Tests/noDist/testAdaptiveChebSymFMM.cpp rename to Tests/noDist/testAdaptiveChebSymFMMDebug.cpp