Skip to content

GitLab

S
ScalFMM
Commit 1b21bf47 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files
Options

Continue the debug

parent b0feeacd
Hide whitespace changes
Inline Side-by-side
Showing with 141 additions and 103 deletions
Src/Adaptive/FAdaptChebSymKernel.hpp
View file @ 1b21bf47
......@@ -62,6 +62,7 @@ class FAdaptiveChebSymKernel : FChebSymKernel<CellClass, ContainerClass, MatrixK
nnodes = TensorTraits<ORDER>::nnodes};
const MatrixKernelClass *const MatrixKernel;
int sminM, sminL;
public:
......@@ -79,11 +80,11 @@ public:
// * runtime_error is thrown if the required file is not valid).
// */
FAdaptiveChebSymKernel(const int inTreeHeight, const FReal inBoxWidth,
const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel), MatrixKernel(inMatrixKernel)
const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel, const int &minM, const int &minL) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel), MatrixKernel(inMatrixKernel),sminM(minM),sminL(minM)
{}
// /** Copy constructor */
FAdaptiveChebSymKernel(const FAdaptiveChebSymKernel& other)
: KernelBaseClass(other), MatrixKernel(other.MatrixKernel)
: KernelBaseClass(other), MatrixKernel(other.MatrixKernel),sminM(other.sminM),sminL(other.sminL)
{ }
//
......@@ -319,14 +320,14 @@ public:
}
bool preferP2M(const ContainerClass* const particles) override {
return particles->getNbParticles() < 10;
return particles->getNbParticles() > this->sminM;
}
bool preferP2M(const int /*atLevel*/, const ContainerClass*const particles[], const int nbContainers) override {
int counterParticles = 0;
for(int idxContainer = 0 ; idxContainer < nbContainers ; ++idxContainer){
counterParticles += particles[idxContainer]->getNbParticles();
}
return counterParticles < 10;
return counterParticles >this->sminM;
}
};
......
Src/Adaptive/FAdaptiveKernelWrapper.hpp
View file @ 1b21bf47
......@@ -196,7 +196,7 @@ public:
FAdaptiveCell<CellClass, ContainerClass>* currentAdaptiveCell = nullptr;
int currentAdaptiveLevel = -1;
// If the current Adaptive cell is the current cell
std::cout << " M2L ->( " <<local->getGlobalId( ) << "): " ;
std::cout << " M2L cell (" <<local->getGlobalId( ) << "): " ;
if(local->isAdaptive()){
currentAdaptiveCell = local;
currentAdaptiveLevel = inLevel;
......@@ -220,20 +220,28 @@ public:
}
else if(currentAdaptiveCell->hasDevelopment()){
// If only current cell has development the neighbor has particles
std::cout << " P2L ->( " <<currentAdaptiveCell->getGlobalId( ) << "): " ;
for(int idxLeafSrc = 0 ; idxLeafSrc < distantNeighbors[idxNeigh]->getNbSubLeaves() ; ++idxLeafSrc){
kernel.P2L(currentAdaptiveCell->getRealCell(), currentAdaptiveLevel, distantNeighbors[idxNeigh]->getSubLeaf(idxLeafSrc));
}
}
else if(distantNeighbors[idxNeigh]->hasDevelopment()){
// If only current cell has particles the neighbor has development
std::cout << " M2P ->( " <<currentAdaptiveCell->getGlobalId( ) << "): " ;
for(int idxLeafTgt = 0 ; idxLeafTgt < currentAdaptiveCell->getNbSubLeaves() ; ++idxLeafTgt){
kernel.M2P(distantNeighbors[idxNeigh]->getRealCell(), currentAdaptiveLevel, currentAdaptiveCell->getSubLeaf(idxLeafTgt));
}
}
else{
// If both have particles
std::cout << " P2P: " ;
for(int idxLeafTgt = 0 ; idxLeafTgt < currentAdaptiveCell->getNbSubLeaves() ; ++idxLeafTgt){
for(int idxLeafSrc = 0 ; idxLeafSrc < distantNeighbors[idxNeigh]->getNbSubLeaves() ; ++idxLeafSrc){
std::cout << " ( " <<distantNeighbors[idxNeigh]->getSubLeaf(idxLeafSrc)<< " ) " ;
kernel.P2P(currentAdaptiveCell->getSubLeaf(idxLeafTgt), distantNeighbors[idxNeigh]->getSubLeaf(idxLeafSrc));
}
}
......@@ -301,10 +309,10 @@ public:
void L2L(const FAdaptiveCell<CellClass, ContainerClass>* const FRestrict local,
FAdaptiveCell<CellClass, ContainerClass>* FRestrict * const FRestrict child, const int inLevel) override {
// If there is something on this cell
std::cout << " L2L ( " <<local->getGlobalId( ) << "): " ;
std::cout << " L2L ( " <<local->getGlobalId( ) << "): ";
std::cout.flush();
if(local->isAdaptive() && local->hasDevelopment()){
std::cout << " ( adaptive) devel("<<std::boolalpha<< local->hasDevelopment()<<")" ;
std::cout << " ( adaptive) devel("<<std::boolalpha<< std::boolalpha<<local->hasDevelopment()<<")" ;
// We store the usual cell
CellClass* realChild[8] = {nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr};
......@@ -356,12 +364,12 @@ public:
/** We do a Local to Particles only if the local (leaf) cell has some development */
void L2P(const FAdaptiveCell<CellClass, ContainerClass>* const local, ContainerClass* const particles) override {
std::cout << " L2P ( " <<local->getGlobalId( ) << "): " << std::boolalpha<<local->hasDevelopment();
if(local->hasDevelopment()){
std::cout << " L2P ( " <<local->getGlobalId( ) << "): " ;
kernel.L2P(local->getRealCell(), particles);
std::cout << std::endl;
}
std::cout << std::endl;
}
/** This is a normal P2P */
......
Src/Kernels/Uniform/FAbstractUnifKernel.hpp
View file @ 1b21bf47
// ===================================================================================
// 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.
//
......@@ -156,8 +156,10 @@ public:
virtual void P2PRemote(const FTreeCoordinate& /*inPosition*/,
ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
ContainerClass* const inNeighbors[27], const int /*inSize*/) = 0;
ContainerClass* const FRestrict inTargets,
const ContainerClass* const FRestrict /*inSources*/,
ContainerClass* const inNeighbors[27],
const int /*inSize*/) = 0;
};
......
Tests/noDist/testAdaptiveChebSymFMM.cpp
View file @ 1b21bf47
......@@ -57,6 +57,7 @@
#include "Core/FFmmAlgorithm.hpp"
//#include "Core/FFmmAlgorithmThread.hpp"
//#include "Core/FFmmAlgorithmTask.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
......@@ -74,9 +75,33 @@ void usage() {
}
// 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);
for (int i = 0 ; i< argc ; ++i){
std::cout << argv[i] << " " ;
}
std::cout << std::endl<< std::endl;
// ---------------------------------------------
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);
//
// accuracy
const unsigned int P = 3 ;
const unsigned int P = 5 ;
// typedef FTestCell CellClass;
// typedef FAdaptiveTestKernel< CellClass, ContainerClass > KernelClass;
typedef FChebCell<P> CellClass;
......@@ -98,10 +123,8 @@ int main(int argc, char ** argv){
std::cout << ">> This executable has to be used to test the FMM algorithm.\n";
//////////////////////////////////////////////////////////////
//
const int NbLevels = FParameters::getValue(argc,argv,"-depth", 7);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-subdepth", 3);
const int sminM = FParameters::getValue(argc,argv,"-sM", P*P*P);
const int sminL = FParameters::getValue(argc,argv,"-sL", P*P*P);
const int sminM = FParameters::getValue(argc,argv,LocalOptionMinMultipoleThreshod.options, P*P*P);
const int sminL = FParameters::getValue(argc,argv,LocalOptionMinLocalThreshod.options, P*P*P);
//
......@@ -110,7 +133,6 @@ int main(int argc, char ** argv){
//////////////////////////////////////////////////////////////////////////////////
// Not Random Loader
//////////////////////////////////////////////////////////////////////////////////
const std::string fileName(FParameters::getStr(argc,argv,"-fin", "../Data/prolate50.out.fma"));
FFmaGenericLoader loader(fileName);
const long int NbPart = loader.getNumberOfParticles() ;
// Random Loader
......@@ -118,109 +140,114 @@ int main(int argc, char ** argv){
// FRandomLoader loader(NbPart, 1, FPoint(0.5,0.5,0.5), 1);
//////////////////////////////////////////////////////////////////////////////////
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Creating & Inserting " << NbPart << " particles ..." << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << 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();
FReal L= loader.getBoxWidth();
//FmaRParticle* particles= new FmaRParticle[NbPart];
FmaRWParticle<8,8>* const particles = new FmaRWParticle<8,8>[NbPart];
FPoint minPos(L,L,L), maxPos(-L,-L,-L);
counter.tic();
FReal L= loader.getBoxWidth();
//
FmaRWParticle<8,8>* const particles = new FmaRWParticle<8,8>[NbPart];
FPoint minPos(L,L,L), maxPos(-L,-L,-L);
//
loader.fillParticle(particles,NbPart);
for(int idxPart = 0 ; idxPart < NbPart; ++idxPart){
const FPoint PP(particles[idxPart].getPosition() ) ;
//
minPos.setX(FMath::Min(minPos.getX(),PP.getX())) ;
minPos.setY(FMath::Min(minPos.getY(),PP.getY())) ;
minPos.setZ(FMath::Min(minPos.getZ(),PP.getZ())) ;
maxPos.setX(FMath::Max(maxPos.getX(),PP.getX())) ;
maxPos.setY(FMath::Max(maxPos.getY(),PP.getY())) ;
maxPos.setZ(FMath::Max(maxPos.getZ(),PP.getZ())) ;
//
loader.fillParticle(particles,NbPart);
for(int idxPart = 0 ; idxPart < NbPart; ++idxPart){
const FPoint PP(particles[idxPart].getPosition() ) ;
//
minPos.setX(FMath::Min(minPos.getX(),PP.getX())) ;
minPos.setY(FMath::Min(minPos.getY(),PP.getY())) ;
minPos.setZ(FMath::Min(minPos.getZ(),PP.getZ())) ;
maxPos.setX(FMath::Max(maxPos.getX(),PP.getX())) ;
maxPos.setY(FMath::Max(maxPos.getY(),PP.getY())) ;
maxPos.setZ(FMath::Max(maxPos.getZ(),PP.getZ())) ;
//
tree.insert(PP, idxPart, particles[idxPart].getPhysicalValue());
counter.tac();
std::cout << "Data are inside the box delimited by "<<std::endl
<< " Min corner: "<< minPos<<std::endl
<< " Max corner: "<< maxPos<<std::endl <<std::endl;
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << " s)." << std::endl;
tree.insert(PP, idxPart, particles[idxPart].getPhysicalValue());
}
counter.tac();
std::cout << "Data are inside the box delimited by "<<std::endl
<< " Min corner: "<< minPos<<std::endl
<< " Max corner: "<< maxPos<<std::endl <<std::endl;
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << " s)." << std::endl;
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Working on particles ..." << std::endl;
// For debug purpose
// Set Global id
//
long int idCell = setGlobalID(tree);
counter.tic();
const MatrixKernelClass MatrixKernel;
KernelWrapperClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel); // FTestKernels FBasicKernels
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;
//
FReal energy= 0.0 , energyD = 0.0 ;
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compute direct energy
/////////////////////////////////////////////////////////////////////////////////////////////////
for(int idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){
energyD += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
}
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compare
/////////////////////////////////////////////////////////////////////////////////////////////////
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
{ // Check that each particle has been summed with all other
// std::cout << "indexPartOrig || DIRECT V fx || FMM V fx" << std::endl;
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];
// std::cout << indexPartOrig
// << " " << particles[indexPartOrig].getPotential() << " " << particles[indexPartOrig].getForces()[0]
// << " " << potentials[idxPart] << " " << forcesX[idxPart]
// << std::endl;
}
});
}
delete[] particles;
// 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;
//
FReal energy= 0.0 , energyD = 0.0 ;
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compute direct energy
/////////////////////////////////////////////////////////////////////////////////////////////////
for(int idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){
energyD += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
}
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compare
/////////////////////////////////////////////////////////////////////////////////////////////////
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
{ // Check that each particle has been summed with all other
// std::cout << "indexPartOrig || DIRECT V fx || FMM V fx" << std::endl;
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];
// std::cout << indexPartOrig
// << " " << particles[indexPartOrig].getPotential() << " " << particles[indexPartOrig].getForces()[0]
// << " " << potentials[idxPart] << " " << forcesX[idxPart]
// << std::endl;
}
});
}
delete[] particles;
// 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;
OctreeClass::Iterator octreeIterator(&tree);
......@@ -236,7 +263,7 @@ int main(int argc, char ** argv){
//
// Set Global id
//
long int idCell = setGlobalID(tree);
// long int idCell = setGlobalID(tree);
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
......@@ -251,7 +278,7 @@ int main(int argc, char ** argv){
std::cout << " Number of Cells: " << idCell <<std::endl;
//
std::cout << "Top of the octree " << octreeIterator.level() << std::endl ;
for(int idxLevel = 1; idxLevel < NbLevels ; ++idxLevel){
for(int idxLevel = 1; idxLevel < TreeHeight ; ++idxLevel){
file << std::endl << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<< std::endl;
file << " Level " << idxLevel <<" Level "<< octreeIterator.level()<< " -- leave level " << std::boolalpha << octreeIterator.isAtLeafLevel() << std::endl;
do{
......
Tests/noDist/testAdaptiveUnifFMM.cpp
View file @ 1b21bf47
......@@ -83,7 +83,7 @@ int main(int argc, char ** argv){
};
FHelpDescribeAndExit(argc, argv,
"Test Uniform kernel and compare it with the direct computation.",
"Test Adaptive kernel and compare it with the direct computation.",
FParameterDefinitions::OctreeHeight,FParameterDefinitions::NbThreads,
FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
LocalOptionMinMultipoleThreshod,LocalOptionMinLocalThreshod);
......
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