Mentions légales du service

Skip to content
Snippets Groups Projects
Commit c14155ac authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files

clean test cheb algorithm

parent 09d0ce91
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Tests/Kernels/testChebAlgorithm.cpp
+ 141
− 294
View file @ c14155ac
...@@ -40,7 +40,9 @@ ...@@ -40,7 +40,9 @@
#include "Containers/FVector.hpp" #include "Containers/FVector.hpp"
#include "Core/FFmmAlgorithm.hpp" #include "Core/FFmmAlgorithm.hpp"
#ifdef _OPENMP
#include "Core/FFmmAlgorithmThread.hpp" #include "Core/FFmmAlgorithmThread.hpp"
#endif
#include "../../Src/Utils/FParameterNames.hpp" #include "../../Src/Utils/FParameterNames.hpp"
...@@ -56,309 +58,154 @@ int main(int argc, char* argv[]) ...@@ -56,309 +58,154 @@ int main(int argc, char* argv[])
FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight, FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::NbThreads); FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::NbThreads);
const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma"); const char* const filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/test20k.fma");
const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5); const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2); const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1); const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
FTic time;
#ifdef _OPENMP #ifdef _OPENMP
omp_set_num_threads(NbThreads); omp_set_num_threads(NbThreads);
std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl; std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
#else #else
std::cout << "\n>> Sequential version.\n" << std:: std::cout << "\n>> Sequential version.\n" << std::
#endif #endif
// init timer // interaction kernel evaluator
FTic time; //typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
// interaction kernel evaluator
//typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
// init particles position and physical value
struct TestParticle{
FPoint position;
FReal forces[3];
FReal physicalValue;
FReal potential;
};
// open particle file
FFmaGenericLoader loader(filename);
if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");
TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
FPoint position;
FReal physicalValue = 0.0;
loader.fillParticle(&position,&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;
}
////////////////////////////////////////////////////////////////////
{ // begin direct computation
time.tic();
{
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential,&MatrixKernel);
}
}
}
time.tac();
printf("Elapsed Time for direct computation: %f\n",time.elapsed());
} // end direct computation
////////////////////////////////////////////////////////////////////
{ // begin Chebyshev kernel
// accuracy
const unsigned int ORDER = 7;
const FReal epsilon = FReal(1e-7);
// typedefs
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FChebCell<ORDER> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
//typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
{ // -----------------------------------------------------
std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
<< " particles ..." << std::endl;
std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
time.tic();
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// put in tree
tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
}
time.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = "
<< time.elapsed() << "s)." << std::endl;
} // -----------------------------------------------------
{ // -----------------------------------------------------
std::cout << "\nChebyshev FMM (ORDER="<< ORDER << ",EPS="<< epsilon <<") ... " << std::endl;
time.tic();
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
FmmClass algorithm(&tree, &kernels);
algorithm.execute();
time.tac();
std::cout << "Done " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
} // -----------------------------------------------------
{ // -----------------------------------------------------
std::cout << "\nError computation ... " << std::endl;
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 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]);
}
});
}
// Print for information
std::cout << "Potential " << potentialDiff << std::endl;
std::cout << "Fx " << fx << std::endl;
std::cout << "Fy " << fy << std::endl;
std::cout << "Fz " << fz << std::endl;
} // -----------------------------------------------------
} // end Chebyshev kernel
//// -----------------------------------------------------
//{ // cost of symmetric m2l opertors, weighted rank, etc.
// const unsigned int nnodes = ORDER*ORDER*ORDER;
// const SymmetryHandler<ORDER> *const SymHandler = kernels.getPtrToSymHandler();
// unsigned int expansionCounter[343];
// for (unsigned i=0; i<343; ++i) expansionCounter[i] = 0;
// for (unsigned i=0; i<343; ++i) if (SymHandler->pindices[i]) expansionCounter[SymHandler->pindices[i]]++;
//
// unsigned int overallCost = 0;
// unsigned int overallWeightedRank = 0;
// unsigned int nbExpansions = 0;
// for (unsigned i=0; i<343; ++i)
// if (expansionCounter[i]) {
// const unsigned int cost = (2*nnodes*SymHandler->LowRank[i]) * expansionCounter[i];
// const unsigned int weightedRank = SymHandler->LowRank[i] * expansionCounter[i];
// overallCost += cost;
// overallWeightedRank += weightedRank;
// nbExpansions += expansionCounter[i];
// std::cout << "expansionCounter[" << i << "] = " << expansionCounter[i]
// << "\tlow rank = " << SymHandler->LowRank[i]
// << "\t(2*nnodes*rank) * nb_exp = " << cost
// << std::endl;
// }
// std::cout << "=== Overall cost = " << overallCost << "\t Weighted rank = " << (double)overallWeightedRank / (double)nbExpansions << std::endl;
// if (nbExpansions!=316) std::cout << "Something went wrong, number of counted expansions = " << nbExpansions << std::endl;
//}
//// -----------------------------------------------------
// -----------------------------------------------------
// find first non empty leaf cell
/*if (FParameters::findParameter(argc,argv,"-dont_check_accuracy") == FParameters::NotFound) {
OctreeClass::Iterator iLeafs(&tree);
iLeafs.gotoBottomLeft();
const ContainerClass *const Targets = iLeafs.getCurrentListTargets();
const unsigned int NumTargets = Targets->getSize();
FReal* Potential = new FReal [NumTargets];
FBlas::setzero(NumTargets, Potential);
FReal* Force = new FReal [NumTargets * 3];
FBlas::setzero(NumTargets * 3, Force);
std::cout << "\nDirect computation of " << NumTargets << " target particles ..." << std::endl;
const MatrixKernelClass MatrixKernel; const MatrixKernelClass MatrixKernel;
do {
const ContainerClass *const Sources = iLeafs.getCurrentListSrc(); // init particles position and physical value
unsigned int counter = 0; struct TestParticle{
ContainerClass::ConstBasicIterator iTarget(*Targets); FPoint position;
while(iTarget.hasNotFinished()) { FReal forces[3];
const FReal wt = iTarget.data().getPhysicalValue(); FReal physicalValue;
ContainerClass::ConstBasicIterator iSource(*Sources); FReal potential;
while(iSource.hasNotFinished()) { };
if (&iTarget.data() != &iSource.data()) {
const FReal potential_value = MatrixKernel.evaluate(iTarget.data().getPosition(), // open particle file
iSource.data().getPosition()); FFmaGenericLoader loader(filename);
const FReal ws = iSource.data().getPhysicalValue(); if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");
// potential
Potential[counter] += potential_value * ws; TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
// force for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
if (MatrixKernelClass::Identifier == ONE_OVER_R) { // laplace force FPoint position;
FPoint force(iSource.data().getPosition() - iTarget.data().getPosition()); FReal physicalValue = 0.0;
force *= ((ws*wt) * (potential_value*potential_value*potential_value)); loader.fillParticle(&position,&physicalValue);
Force[counter*3 + 0] += force.getX(); // get copy
Force[counter*3 + 1] += force.getY(); particles[idxPart].position = position;
Force[counter*3 + 2] += force.getZ(); particles[idxPart].physicalValue = physicalValue;
} else if (MatrixKernelClass::Identifier == LEONARD_JONES_POTENTIAL) { // lenard-jones force particles[idxPart].potential = 0.0;
FPoint force(iSource.data().getPosition() - iTarget.data().getPosition()); particles[idxPart].forces[0] = 0.0;
const FReal one_over_r = FReal(1.) / FMath::Sqrt(force.getX()*force.getX() + particles[idxPart].forces[1] = 0.0;
force.getY()*force.getY() + particles[idxPart].forces[2] = 0.0;
force.getZ()*force.getZ()); // 1 + 15 + 5 = 21 flops
const FReal one_over_r3 = one_over_r * one_over_r * one_over_r;
const FReal one_over_r6 = one_over_r3 * one_over_r3;
const FReal one_over_r4 = one_over_r3 * one_over_r;
force *= ((ws*wt) * (FReal(12.)*one_over_r6*one_over_r4*one_over_r4 - FReal(6.)*one_over_r4*one_over_r4));
Force[counter*3 + 0] += force.getX();
Force[counter*3 + 1] += force.getY();
Force[counter*3 + 2] += force.getZ();
}
}
iSource.gotoNext();
}
counter++;
iTarget.gotoNext();
} }
} while(iLeafs.moveRight());
////////////////////////////////////////////////////////////////////
{ // begin direct computation
time.tic();
{
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential,&MatrixKernel);
}
}
}
time.tac();
printf("Elapsed Time for direct computation: %f\n",time.elapsed());
} // end direct computation
FReal* ApproxPotential = new FReal [NumTargets]; ////////////////////////////////////////////////////////////////////
FReal* ApproxForce = new FReal [NumTargets * 3];
{ // begin Chebyshev kernel
const unsigned int ORDER = 7;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FChebCell<ORDER> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
//typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
#ifdef _OPENMP
typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#else
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#endif
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
{ // -----------------------------------------------------
std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
<< " particles ..." << std::endl;
std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
time.tic();
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// put in tree
tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
}
time.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = "
<< time.elapsed() << "s)." << std::endl;
} // -----------------------------------------------------
{ // -----------------------------------------------------
std::cout << "\nChebyshev FMM (ORDER="<< ORDER << "... " << std::endl;
time.tic();
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
FmmClass algorithm(&tree, &kernels);
algorithm.execute();
time.tac();
std::cout << "Done " << "(@Algorithm = " << time.elapsed() << "s)." << std::endl;
} // -----------------------------------------------------
{ // -----------------------------------------------------
std::cout << "\nError computation ... " << std::endl;
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 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]);
}
});
}
// Print for information
std::cout << "Potential " << potentialDiff << std::endl;
std::cout << "Fx " << fx << std::endl;
std::cout << "Fy " << fy << std::endl;
std::cout << "Fz " << fz << std::endl;
} // -----------------------------------------------------
} // end Chebyshev kernel
unsigned int counter = 0;
ContainerClass::ConstBasicIterator iTarget(*Targets);
while(iTarget.hasNotFinished()) {
ApproxPotential[counter] = iTarget.data().getPotential();
ApproxForce[counter*3 + 0] = iTarget.data().getForces().getX();
ApproxForce[counter*3 + 1] = iTarget.data().getForces().getY();
ApproxForce[counter*3 + 2] = iTarget.data().getForces().getZ();
counter++;
iTarget.gotoNext();
}
std::cout << "\nPotential error:" << std::endl; return 0;
std::cout << "Relative L2 error = " << computeL2norm( NumTargets, Potential, ApproxPotential)
<< std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(NumTargets, Potential, ApproxPotential)
<< std::endl;
std::cout << "\nForce error:" << std::endl;
std::cout << "Relative L2 error = " << computeL2norm( NumTargets*3, Force, ApproxForce)
<< std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(NumTargets*3, Force, ApproxForce)
<< std::endl;
std::cout << std::endl;
// free memory
delete [] Potential;
delete [] ApproxPotential;
delete [] Force;
delete [] ApproxForce;
}*/
/*
// Check if particles are strictly within its containing cells
const FReal BoxWidthLeaf = loader.getBoxWidth() / FReal(FMath::pow(2, TreeHeight-1));
OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
const CellClass *const LeafCell = octreeIterator.getCurrentCell();
const FPoint LeafCellCenter(LeafCell->getCoordinate().getX() * BoxWidthLeaf + BoxWidthLeaf/2 + loader.getCenterOfBox().getX(),
LeafCell->getCoordinate().getY() * BoxWidthLeaf + BoxWidthLeaf/2 + loader.getCenterOfBox().getY(),
LeafCell->getCoordinate().getZ() * BoxWidthLeaf + BoxWidthLeaf/2 + loader.getCenterOfBox().getZ());
const ContainerClass *const Particles = octreeIterator.getCurrentListSrc();
ContainerClass::ConstBasicIterator particleIterator(*Particles);
while(particleIterator.hasNotFinished()) {
const FPoint distance(LeafCellCenter-particleIterator.data().getPosition());
std::cout << "center - particle = " << distance << " < " << BoxWidthLeaf/FReal(2.) << std::endl;
if (std::abs(distance.getX())>BoxWidthLeaf/FReal(2.) ||
std::abs(distance.getY())>BoxWidthLeaf/FReal(2.) ||
std::abs(distance.getZ())>BoxWidthLeaf/FReal(2.)) {
std::cout << "stop" << std::endl;
exit(-1);
}
particleIterator.gotoNext();
}
} while(octreeIterator.moveRight());
*/
return 0;
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment