Commit ed193e2d authored by COULAUD Olivier's avatar COULAUD Olivier

Now test uestSphericalDirect works

Rename utestChebyshevMultiRhs in utestInterpolationMultiRhs and test on Lagrange works but not Chebyschev one.
parent 7bc2123a
......@@ -30,6 +30,7 @@ class FTreeCoordinate;
// ==== CMAKE =====
// @FUSE_BLAS
// @FUSE_FFT
// ================
// for verbosity only!!!
......@@ -221,7 +222,7 @@ public:
// Target cell: local
const FReal localCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, localLevel)));
const FPoint localCellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel));
std::cout << " call P2L localLevel "<< localLevel << " localCellCenter "<< localCellCenter <<std::endl;
// std::cout << " call P2L localLevel "<< localLevel << " localCellCenter "<< localCellCenter <<std::endl;
// interpolation points of target (X) cell
FPoint X[nnodes];
FUnifTensor<order>::setRoots(localCellCenter, localCellWidth, X);
......@@ -458,7 +459,7 @@ public:
for(int idxContainer = 0 ; idxContainer < nbContainers ; ++idxContainer){
counterParticles += particles[idxContainer]->getNbParticles();
}
std::cout << " Part("<<counterParticles<< ") ";
// std::cout << " Part("<<counterParticles<< ") ";
return counterParticles >this->sminM;
}
};
......
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// 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.
//
......@@ -24,16 +24,16 @@
* To chose which operation has to be performed.
*/
enum FFmmOperations {
FFmmP2P = (1 << 0),
FFmmP2M = (1 << 1),
FFmmP2P = (1 << 0),
FFmmP2M = (1 << 1),
FFmmM2M = (1 << 2),
FFmmM2L = (1 << 3),
FFmmL2L = (1 << 4),
FFmmL2P = (1 << 5),
FFmmM2L = (1 << 3),
FFmmL2L = (1 << 4),
FFmmL2P = (1 << 5),
//
FFmmNearField = FFmmP2P,
FFmmFarField = (FFmmP2M|FFmmM2M|FFmmM2L|FFmmL2L|FFmmL2P),
//
FFmmNearAndFarFields = (FFmmNearField|FFmmFarField)
};
......@@ -52,12 +52,15 @@ protected:
int nbLevelsInTree;
void setNbLevelsInTree(const int inNbLevelsInTree){
nbLevelsInTree = inNbLevelsInTree;
nbLevelsInTree = inNbLevelsInTree;
lowerWorkingLevel = nbLevelsInTree;
}
void validateLevels() const {
std::cout << "upperWorkingLevel: "<< FAbstractAlgorithm::upperWorkingLevel << std::endl
<< "lowerWorkingLevel "<< FAbstractAlgorithm::lowerWorkingLevel << std::endl ;
FAssertLF(FAbstractAlgorithm::upperWorkingLevel <= FAbstractAlgorithm::lowerWorkingLevel);
std::cout << "End assert 1" << std::endl ;
FAssertLF(2 <= FAbstractAlgorithm::upperWorkingLevel);
}
......
......@@ -17,6 +17,7 @@
// ==== CMAKE =====
// @FUSE_BLAS
// ==============
#include <array>
#include "ScalFmmConfig.h"
#include "Utils/FGlobal.hpp"
......@@ -35,29 +36,33 @@
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
/*
In this test we compare the spherical FMM results and the direct results.
*/
*/
/** the test class
*
*/
class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
class TestInterpolationKernel : public FUTester<TestInterpolationKernel> {
///////////////////////////////////////////////////////////
// The tests!
///////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////
// The tests!
///////////////////////////////////////////////////////////
template <class CellClass, class ContainerClass, class KernelClass, class MatrixKernelClass,
class LeafClass, class OctreeClass, class FmmClass, const int NVals>
void RunTest() {
// Warning in make test the exec dir it Build/UTests
// Load particles
template <class CellClass, class ContainerClass, class KernelClass, class MatrixKernelClass,
class LeafClass, class OctreeClass, class FmmClass, const int NVals>
void RunTest() {
// Warning in make test the exec dir it Build/UTests
// Load particles
//
// Load particles
//
......@@ -72,101 +77,113 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
std::string filename(SCALFMMDataPath+parFile);
//
FFmaGenericLoader loader(filename);
if(!loader.isOpen()){
Print("Cannot open particles file.");
uassert(false);
return;
}
Print("Number of particles:");
Print(loader.getNumberOfParticles());
const int NbLevels = 4;
const int SizeSubLevels = 2;
// Create Matrix Kernel
const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
//
if(!loader.isOpen()){
Print("Cannot open particles file.");
uassert(false);
return;
}
Print("Number of particles:");
Print(loader.getNumberOfParticles());
const int NbLevels = 4;
const int SizeSubLevels = 2;
// std::cout << "\nInterpolation FMM (ORDER="<< ORDER << ") ... " << std::endl;
// Create Matrix Kernel
const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
//
FSize nbParticles = loader.getNumberOfParticles() ;
FmaRWParticle<8,8>* const particles = new FmaRWParticle<8,8>[nbParticles];
loader.fillParticle(particles,nbParticles);
//
//
// Create octree
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
// Insert particle in the tree
//
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
tree.insert(particles[idxPart].getPosition() , idxPart, particles[idxPart].getPhysicalValue() ,0.0,0.0,0.0);
// For each particles we associate Nvals charge ( q,0,0,0)
//
/* for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
double q = particles[idxPart].getPhysicalValue();
tree.insert(particles[idxPart].getPosition() , idxPart, q,q);//,0.0,0.0,0.0);
}*/
for(int idxPart = 0 ; idxPart < nbParticles ; ++idxPart){
// // Convert FReal[NVALS] to std::array<FReal,NVALS>
std::array<FReal, (1+4*1)*NVals> physicalState;
for(int idxVals = 0 ; idxVals < NVals ; ++idxVals){
double q = particles[idxPart].getPhysicalValue();
physicalState[0*NVals+idxVals]= particles[idxPart].getPhysicalValue();
physicalState[1*NVals+idxVals]=0.0;
physicalState[2*NVals+idxVals]=0.0;
physicalState[3*NVals+idxVals]=0.0;
physicalState[4*NVals+idxVals]=0.0;
}
// put in tree
tree.insert(particles[idxPart].getPosition(), idxPart, physicalState);
}
// //
// // Create octree
// for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// FPoint position;
// FReal physicalValue;
// loader.fillParticle(&position,&physicalValue);
// // put in tree
// tree.insert(position, idxPart, physicalValue);
// // get copy
// particles[idxPart].position = position;
// particles[idxPart].physicalValue = physicalValue;
// particles[idxPart].potential = 0.0;
// particles[idxPart].forces[0] = 0.0;
// particles[idxPart].forces[1] = 0.0;
// particles[idxPart].forces[2] = 0.0;
// }
// Run FMM
Print("Fmm...");
KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
FmmClass algo(&tree,&kernels);
algo.execute();
//
FReal energy= 0.0 , energyD = 0.0 ;
// Run direct computation
Print("Direct...");
for(int idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){
*(particles[idx].setPotential()) *= NVals ;
particles[idx].getForces()[0] *= NVals ;
particles[idx].getForces()[1] *= NVals ;
particles[idx].getForces()[2] *= NVals ;
energyD += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
}
//
// Compare
Print("Compute Diff...");
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
{ // Check that each particle has been summed with all other
//
// Run FMM
Print("Fmm...");
KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
FmmClass algo(&tree,&kernels);
algo.execute();
//
FReal energy= 0.0 , energyD = 0.0 ;
for(int idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){
energyD += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
}
//
// Compare
Print("Compute Diff...");
FMath::FAccurater potentialDiff[NVals];
FMath::FAccurater fx, fy, fz;
{
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].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];
//
for(int idxVals = 0 ; idxVals < NVals ; ++idxVals){
const FReal* const physicalValues = leaf->getTargets()->getPhysicalValues(idxVals);
const FReal*const potentials = leaf->getTargets()->getPotentials(idxVals);
const FReal*const forcesX = leaf->getTargets()->getForcesX(idxVals);
const FReal*const forcesY = leaf->getTargets()->getForcesY(idxVals);
const FReal*const forcesZ = leaf->getTargets()->getForcesZ(idxVals);
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];
// std::cout << " index "<< indexPartOrig << " " << particles[indexPartOrig].getPotential() << " " << potentials[idxPart] << std::endl;
potentialDiff[idxVals].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];
}
}
});
}
}
delete[] particles;
// Print for information
energy /=NVals;
// Print for information
double errorPotRL2=0.0, errorPotRMS=0.0;
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());
for(int idxVals = 0 ; idxVals < NVals ; ++idxVals){
printf(" Charge: %d\n", idxVals);
printf(" Pot L2Norm %e\n",potentialDiff[idxVals].getL2Norm());
printf(" Pot RL2Norm %e\n",potentialDiff[idxVals].getRelativeL2Norm());
printf(" Pot RMSError %e\n",potentialDiff[idxVals].getRMSError());
errorPotRL2 = std::max(errorPotRL2, potentialDiff[idxVals].getRelativeL2Norm());
errorPotRMS = std::max(errorPotRMS, potentialDiff[idxVals].getRMSError());
}
Print("Fx diff is = ");
printf(" Fx L2Norm %e\n",fx.getL2Norm());
printf(" Fx RL2Norm %e\n",fx.getRelativeL2Norm());
......@@ -185,14 +202,14 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
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);
// 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
uassert(errorPotRL2 < MaximumDiffPotential); //1
Print("Test2 - Error RMS L2 norm Potential ");
uassert(potentialDiff.getRMSError() < MaximumDiffPotential); //2
uassert(errorPotRMS< MaximumDiffPotential); //2
Print("Test3 - Error Relative L2 norm FX ");
uassert(fx.getRelativeL2Norm() < MaximumDiffForces); //3
Print("Test4 - Error RMS L2 norm FX ");
......@@ -211,92 +228,96 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
uassert(FMath::Abs(energy-energyD) /energyD< MaximumDiffPotential); //10 Total Energy
// Compute multipole local rhs diff
FMath::FAccurater localDiff;
FMath::FAccurater multiPoleDiff;
tree.forEachCell([&](CellClass* cell){
for( int idxRhs = 1 ; idxRhs < NVals ; ++idxRhs){
localDiff.add(cell->getLocal(0), cell->getLocal(idxRhs), cell->getVectorSize());
multiPoleDiff.add(cell->getMultipole(0), cell->getMultipole(idxRhs), cell->getVectorSize());
}
});
Print("Local diff is = ");
Print(localDiff.getL2Norm());
Print(localDiff.getInfNorm());
Print("Multipole diff is = ");
Print(multiPoleDiff.getL2Norm());
Print(multiPoleDiff.getInfNorm());
uassert(localDiff.getL2Norm() < 1e-10);
uassert(localDiff.getInfNorm() < 1e-10);
uassert(multiPoleDiff.getL2Norm() < 1e-10);
uassert(multiPoleDiff.getInfNorm() < 1e-10);
}
/** 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 int NVals = 4;
const unsigned int ORDER = 6 ;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER, 1, 1, NVals> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER, NVals> KernelClass;
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// run test
RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass, NVals>();
}
/** TestChebSymKernel */
void TestChebSymKernel(){
const int NVals = 4;
const unsigned int ORDER = 6;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER, 1, 1, NVals> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER, NVals> KernelClass;
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// run test
RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass, NVals>();
}
///////////////////////////////////////////////////////////
// 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");
}
// Compute multipole local rhs diff
FMath::FAccurater localDiff;
FMath::FAccurater multiPoleDiff;
tree.forEachCell([&](CellClass* cell){
for( int idxRhs = 1 ; idxRhs < NVals ; ++idxRhs){
localDiff.add(cell->getLocal(0), cell->getLocal(idxRhs), cell->getVectorSize());
multiPoleDiff.add(cell->getMultipole(0), cell->getMultipole(idxRhs), cell->getVectorSize());
}
});
Print("Local diff is = ");
Print(localDiff.getL2Norm());
Print(localDiff.getInfNorm());
Print("Multipole diff is = ");
Print(multiPoleDiff.getL2Norm());
Print(multiPoleDiff.getInfNorm());
uassert(localDiff.getL2Norm() < 1e-10);
uassert(localDiff.getInfNorm() < 1e-10);
uassert(multiPoleDiff.getL2Norm() < 1e-10);
uassert(multiPoleDiff.getInfNorm() < 1e-10);
}
/** 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 int NVals = 1;
const unsigned int ORDER = 6 ;
// run test
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FP2PParticleContainerIndexed<1,1,NVals> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FUnifCell<ORDER,1,1,NVals> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FUnifKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER,NVals> KernelClass;
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass, NVals>();
}
/** TestChebSymKernel */
void TestChebSymKernel(){
const int NVals = 2;
const unsigned int ORDER = 6;
typedef FP2PParticleContainerIndexed<1,1,NVals> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER, 1, 1, NVals> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER, NVals> KernelClass;
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// run test
RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass, NVals>();
}
///////////////////////////////////////////////////////////
// Set the tests!
///////////////////////////////////////////////////////////
/** set test */
void SetTests(){
AddTest(&TestInterpolationKernel::TestUnifKernel,"Test Lagrange/Uniform grid FMM");
AddTest(&TestInterpolationKernel::TestChebSymKernel,"Test Symmetric Chebyshev Kernel with 16 small SVDs and symmetries");
}
};
// You must do this
TestClass(TestChebyshevDirect)
TestClass(TestInterpolationKernel)
......
......@@ -330,12 +330,12 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
OctreeClass, FmmClass, 24>(true);
RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
OctreeClass, FmmClass, 26>(true);
RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
OctreeClass, FmmClass, 28>(true);
RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
OctreeClass, FmmClass, 30>(true);
RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
OctreeClass, FmmClass, 32>(true);
// RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
// OctreeClass, FmmClass, 28>(true);
// RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
// OctreeClass, FmmClass, 30>(true);
// RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
// OctreeClass, FmmClass, 32>(true);
}
#endif
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment