From 5367316b28b03d279843e4ece8056ecf5ac2a45b Mon Sep 17 00:00:00 2001 From: Berenger Bramas Date: Fri, 6 Nov 2015 21:58:36 +0100 Subject: [PATCH] add a tsm utest on the chebyshev kernel --- UTests/utestChebyshevDirectTsm.cpp | 275 +++++++++++++++++++++++++++++ 1 file changed, 275 insertions(+) create mode 100644 UTests/utestChebyshevDirectTsm.cpp diff --git a/UTests/utestChebyshevDirectTsm.cpp b/UTests/utestChebyshevDirectTsm.cpp new file mode 100644 index 00000000..eccd79f3 --- /dev/null +++ b/UTests/utestChebyshevDirectTsm.cpp @@ -0,0 +1,275 @@ + +// =================================================================================== +// Copyright ScalFmm 2011 INRIA, +// 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 "Utils/FGlobal.hpp" + +#include "Containers/FOctree.hpp" +#include "Containers/FVector.hpp" + +#include "Kernels/Chebyshev/FChebCell.hpp" +#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp" +#include "Kernels/Interpolation/FInterpMatrixKernel.hpp" + +#include "Components/FTypedLeaf.hpp" +#include "Extensions/FExtendCellType.hpp" +#include "Kernels/Chebyshev/FChebSymKernel.hpp" + +#include "Files/FRandomLoader.hpp" +#include "Files/FFmaGenericLoader.hpp" + +#include "Core/FFmmAlgorithmThreadTsm.hpp" +#include "Core/FFmmAlgorithmTsm.hpp" + +#include "FUTester.hpp" + + +/** the test class the rotation and target source model. + * + */ +class TestChebyshevDirectTsm : public FUTester { + /** The test method to factorize all the test based on different kernels */ + template + void RunTest(){ + const MatrixKernelClass MatrixKernel; + // Warning in make test the exec dir it Build/UTests + // Load particles + const int nbSources = 5000; + const int nbTargets = 5000; + + FRandomLoader loader(nbSources + nbTargets); + + Print("Number of particles:"); + Print(loader.getNumberOfParticles()); + + const int NbLevels = 4; + const int SizeSubLevels = 3; + + // Create octree + OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox()); + + const FReal physicalValue = 0.10; + // + FmaRWParticle* const particlesTargets = new FmaRWParticle[nbTargets]; + for(FSize idxPart = 0 ; idxPart < nbTargets ; ++idxPart){ + FPoint position; + loader.fillParticle(&position); + // put in tree + tree.insert(position, FParticleTypeTarget, idxPart, physicalValue); + // get copy + particlesTargets[idxPart].setPosition(position); + *(particlesTargets[idxPart].setPhysicalValue()) = physicalValue; + *(particlesTargets[idxPart].setPotential()) = 0.0; + particlesTargets[idxPart].setForces()[0] = 0.0; + particlesTargets[idxPart].setForces()[1] = 0.0; + particlesTargets[idxPart].setForces()[2] = 0.0; + } + + FmaRWParticle* const particlesSources = new FmaRWParticle[nbSources]; + for(FSize idxPart = 0 ; idxPart < nbSources ; ++idxPart){ + FPoint position; + loader.fillParticle(&position); + // put in tree + tree.insert(position, FParticleTypeSource, idxPart, physicalValue); + // get copy + particlesSources[idxPart].setPosition(position); + *(particlesSources[idxPart].setPhysicalValue()) = physicalValue; + } + + + // Run FMM + Print("Fmm..."); + //KernelClass kernels(NbLevels,loader.getBoxWidth()); + KernelClass kernels(NbLevels,loader.getBoxWidth(), loader.getCenterOfBox(), &MatrixKernel); + FmmClass algo(&tree,&kernels); + algo.execute(); + // + // + + // Run direct computation + Print("Direct..."); + for(int idxTarget = 0 ; idxTarget < nbTargets ; ++idxTarget){ + for(int idxOther = 0 ; idxOther < nbSources ; ++idxOther){ + FP2P::NonMutualParticles( + particlesSources[idxOther].getPosition().getX(), particlesSources[idxOther].getPosition().getY(), + particlesSources[idxOther].getPosition().getZ(),particlesSources[idxOther].getPhysicalValue(), + particlesTargets[idxTarget].getPosition().getX(), particlesTargets[idxTarget].getPosition().getY(), + particlesTargets[idxTarget].getPosition().getZ(),particlesTargets[idxTarget].getPhysicalValue(), + &particlesTargets[idxTarget].setForces()[0],&particlesTargets[idxTarget].setForces()[1], + &particlesTargets[idxTarget].setForces()[2],particlesTargets[idxTarget].setPotential(), + &MatrixKernel); + } + } + + // + // Assert + ///////////////////////////////////////////////////////////////////////////////////////////////// + // Compute direct energy + ///////////////////////////////////////////////////////////////////////////////////////////////// + FReal energy= 0.0 , energyD = 0.0 ; + for(int idx = 0 ; idx < nbTargets ; ++idx){ + energyD += particlesTargets[idx].getPotential()*particlesTargets[idx].getPhysicalValue() ; + } + ///////////////////////////////////////////////////////////////////////////////////////////////// + // 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){ + if( leaf->getTargets()->getNbParticles() ){ + 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 FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles(); + const FVector& indexes = leaf->getTargets()->getIndexes(); + + for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){ + const FSize indexPartOrig = indexes[idxPart]; + potentialDiff.add(particlesTargets[indexPartOrig].getPotential(),potentials[idxPart]); + fx.add(particlesTargets[indexPartOrig].getForces()[0],forcesX[idxPart]); + fy.add(particlesTargets[indexPartOrig].getForces()[1],forcesY[idxPart]); + fz.add(particlesTargets[indexPartOrig].getForces()[2],forcesZ[idxPart]); + energy += potentials[idxPart]*physicalValues[idxPart]; + } + } + }); + } + + delete[] particlesTargets; + delete[] particlesSources; + // + // 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"; + } + } + + /////////////////////////////////////////////////////////// + // The tests! + /////////////////////////////////////////////////////////// + + static const int P = 5; + + /** Chebyshev */ + void TestChebyshev(){ + typedef double FReal; + typedef FInterpMatrixKernelR MatrixKernelClass; + typedef FTypedChebCell CellClass; + typedef FP2PParticleContainerIndexed ContainerClass; + + typedef FChebSymKernel KernelClass; + + typedef FTypedLeaf LeafClass; + typedef FOctree OctreeClass; + + typedef FFmmAlgorithmTsm FmmClass; + + RunTest(); + } + + void TestChebyshevThread(){ + typedef double FReal; + typedef FInterpMatrixKernelR MatrixKernelClass; + typedef FTypedChebCell CellClass; + typedef FP2PParticleContainerIndexed ContainerClass; + + typedef FChebSymKernel KernelClass; + + typedef FTypedLeaf LeafClass; + typedef FOctree OctreeClass; + + typedef FFmmAlgorithmThreadTsm FmmClass; + + RunTest(); + } + + /////////////////////////////////////////////////////////// + // Set the tests! + /////////////////////////////////////////////////////////// + + /** set test */ + void SetTests(){ + AddTest(&TestChebyshevDirectTsm::TestChebyshev,"Test Chebyshev Kernel TSM"); + AddTest(&TestChebyshevDirectTsm::TestChebyshevThread,"Test Chebyshev Kernel TSM thread"); + } +}; + + +// You must do this +TestClass(TestChebyshevDirectTsm) + + + + -- GitLab