Commit 86e28572 authored by berenger-bramas's avatar berenger-bramas

Finalize the periodic system and add some tests to validate.

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@257 2616d619-271b-44dc-8df4-d4a8f33a7222
parent 32021c8d
......@@ -100,7 +100,7 @@ public:
for(int indexLevel = 0; indexLevel < this->height; ++indexLevel ){
this->boxWidthAtLevel[indexLevel] = tempWidth;
tempWidth /= FReal(2.0);
}
}
}
/** Desctructor */
......
......@@ -30,6 +30,7 @@ class FFmmAlgorithmPeriodic : protected FAssertable{
KernelClass* const kernels; //< The kernels
const int OctreeHeight;
const int periodicLimit;
public:
/** The constructor need the octree and the kernels used for computation
......@@ -37,8 +38,8 @@ public:
* @param inKernels the kernels to call
* An assert is launched if one of the arguments is null
*/
FFmmAlgorithmPeriodic(OctreeClass* const inTree, KernelClass* const inKernels)
: tree(inTree) , kernels(inKernels), OctreeHeight(tree->getHeight()) {
FFmmAlgorithmPeriodic(OctreeClass* const inTree, KernelClass* const inKernels, const int inPeriodicLimit = 0)
: tree(inTree) , kernels(inKernels), OctreeHeight(tree->getHeight()), periodicLimit(inPeriodicLimit) {
fassert(tree, "tree cannot be null", __LINE__, __FILE__);
fassert(kernels, "kernels cannot be null", __LINE__, __FILE__);
......@@ -242,8 +243,11 @@ private:
/** Periodicity */
void processPeriodicLevels(){
const int PeriodicLimit = 10;
CellClass upperCells[PeriodicLimit];
if( !periodicLimit ){
return;
}
CellClass upperCells[periodicLimit];
// First M2M from level 1 to level 0
{
......@@ -254,7 +258,7 @@ private:
// Then M2M from level 0 to level -LIMITE
{
CellClass* virtualChild[8];
for(int idxLevel = 1 ; idxLevel < PeriodicLimit ; ++idxLevel){
for(int idxLevel = 1 ; idxLevel < periodicLimit ; ++idxLevel){
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
virtualChild[idxChild] = &upperCells[idxLevel-1];
}
......@@ -282,7 +286,7 @@ private:
const CellClass* neighbors[189];
const int counter = 189;
for(int idxLevel = 0 ; idxLevel < PeriodicLimit ; ++idxLevel ){
for(int idxLevel = 0 ; idxLevel < periodicLimit ; ++idxLevel ){
for(int idxNeigh = 0 ; idxNeigh < 189 ; ++idxNeigh){
neighbors[idxNeigh] = &upperCells[idxLevel];
}
......@@ -295,7 +299,7 @@ private:
{
CellClass* virtualChild[8];
memset(virtualChild, 0, sizeof(CellClass*) * 8);
for(int idxLevel = PeriodicLimit - 1 ; idxLevel > 0 ; --idxLevel){
for(int idxLevel = periodicLimit - 1 ; idxLevel > 0 ; --idxLevel){
virtualChild[0] = &upperCells[idxLevel-1];
kernels->L2L( &upperCells[idxLevel], virtualChild, -idxLevel);
}
......
......@@ -85,10 +85,9 @@ 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,"-h", 9);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 3);
const long NbPart = FParameters::getValue(argc,argv,"-nb", 2000000);
const FReal FRandMax = FReal(RAND_MAX);
const int NbLevels = FParameters::getValue(argc,argv,"-h", 7);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 3);
const long NbPartPerBoxes = FParameters::getValue(argc,argv,"-nb", 3);
FTic counter;
......@@ -97,20 +96,39 @@ int main(int argc, char ** argv){
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
OctreeClass tree(NbLevels, SizeSubLevels, 1.0, F3DPosition(0.5,0.5,0.5));
const FReal BoxWidth = 1.0;
const FReal BoxCenter = 0.5;
OctreeClass tree(NbLevels, SizeSubLevels, BoxWidth, F3DPosition(BoxCenter,BoxCenter,BoxCenter));
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Creating & Inserting " << NbPart << " particles ..." << std::endl;
std::cout << "Creating & Inserting " << NbPartPerBoxes << " particles per boxes ..." << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
counter.tic();
long NbPart = 0;
{
const long NbSmallBoxesPerSide = (1 << (NbLevels-1));
const FReal SmallBoxWidth = BoxWidth / FReal(NbSmallBoxesPerSide);
const FReal SmallBoxWidthDiv2 = SmallBoxWidth / 2;
NbPart = NbSmallBoxesPerSide * NbSmallBoxesPerSide * NbSmallBoxesPerSide * NbPartPerBoxes;
FTestParticle particleToFill;
for(int idxPart = 0 ; idxPart < NbPart ; ++idxPart){
particleToFill.setPosition(FReal(rand())/FRandMax,FReal(rand())/FRandMax,FReal(rand())/FRandMax);
tree.insert(particleToFill);
for(int idxX = 0 ; idxX < NbSmallBoxesPerSide ; ++idxX){
for(int idxY = 0 ; idxY < NbSmallBoxesPerSide ; ++idxY){
for(int idxZ = 0 ; idxZ < NbSmallBoxesPerSide ; ++idxZ){
particleToFill.setPosition(FReal(idxX)*SmallBoxWidth + SmallBoxWidthDiv2,
FReal(idxY)*SmallBoxWidth + SmallBoxWidthDiv2,
FReal(idxZ)*SmallBoxWidth + SmallBoxWidthDiv2);
for(int idxPart = 0 ; idxPart < NbPartPerBoxes ; ++idxPart){
tree.insert(particleToFill);
}
}
}
}
}
......@@ -123,10 +141,11 @@ int main(int argc, char ** argv){
std::cout << "Working on particles ..." << std::endl;
counter.tic();
const int PeriodicDeep = 1;
// FTestKernels FBasicKernels
KernelClass kernels;
//FFmmAlgorithm FFmmAlgorithmThread
FmmClass algo(&tree,&kernels);
FmmClass algo( &tree, &kernels, PeriodicDeep);
algo.execute();
counter.tac();
......@@ -137,6 +156,107 @@ int main(int argc, char ** argv){
//ValidateFMMAlgo<OctreeClass, ParticleClass, CellClass, ContainerClass, LeafClass>(&tree);
{ // Check that each particle has been summed with all other
long counterNbPart = 0;
typename OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
// on each leaf we should have the same number of particles
if(octreeIterator.getCurrentCell()->getDataUp() != octreeIterator.getCurrentListSrc()->getSize()
|| NbPartPerBoxes != octreeIterator.getCurrentCell()->getDataUp() ){
std::cout << "Problem P2M NbPartPerBoxes = " << NbPartPerBoxes <<
" Data up = " << octreeIterator.getCurrentCell()->getDataUp() <<
" Size = " << octreeIterator.getCurrentListSrc()->getSize() << "\n";
}
// we also count the number of particles.
counterNbPart += octreeIterator.getCurrentListSrc()->getSize();
} while(octreeIterator.moveRight());
if( counterNbPart != NbPart){
std::cout << "Problem global nb part, counter = " << counterNbPart << " created = " << NbPart << std::endl;
}
}
{ // Ceck if there is number of NbPart summed at level 1
long long particlesPerBox = NbPartPerBoxes;
typename OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
for(int idxLevel = NbLevels - 1 ; idxLevel >= 1 ; --idxLevel ){
do{
if(octreeIterator.getCurrentCell()->getDataUp() != particlesPerBox){
std::cout << "Problem M2M particlesPerBox = " << particlesPerBox <<
" Data up = " << octreeIterator.getCurrentCell()->getDataUp() <<
" level = " << idxLevel << "\n";
}
} while(octreeIterator.moveRight());
octreeIterator.moveUp();
octreeIterator.gotoLeft();
particlesPerBox *= 8;
}
}
if( PeriodicDeep ){
long long int counterL2L = 0;
{
long long int particlesPerBox = NbPartPerBoxes * FMath::pow(8,NbLevels-1);
long long int counterUp[PeriodicDeep];
for( int idxLevel = 0 ; idxLevel < PeriodicDeep ; ++idxLevel ){
counterUp[idxLevel] = particlesPerBox;
particlesPerBox *= 8;
}
long long int counterM2L[PeriodicDeep];
for( int idxLevel = 0 ; idxLevel < PeriodicDeep ; ++idxLevel ){
counterM2L[idxLevel] = counterUp[idxLevel] * 189;
}
for( int idxLevel = PeriodicDeep - 1 ; idxLevel >= 0 ; --idxLevel ){
counterL2L += counterM2L[idxLevel];
}
}
{
long long int particlesPerBox = NbPartPerBoxes * FMath::pow(8,NbLevels-2);
typename OctreeClass::Iterator octreeIterator(&tree);
for(int idxLevel = 1 ; idxLevel < NbLevels ; ++idxLevel ){
counterL2L = particlesPerBox * 189 + counterL2L;
do{
if(octreeIterator.getCurrentCell()->getDataDown() != counterL2L){
std::cout << "Problem L2L counterL2L = " << counterL2L <<
" Data Down = " << octreeIterator.getCurrentCell()->getDataDown() <<
" level = " << idxLevel << "\n";
}
} while(octreeIterator.moveRight());
octreeIterator.gotoLeft();
octreeIterator.moveDown();
particlesPerBox /= 8;
}
}
}
{ // Check that each particle has been summed with all other
typename OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
typename ContainerClass::BasicIterator iter(*octreeIterator.getCurrentListTargets());
const long long int sumParticles = octreeIterator.getCurrentCell()->getDataDown() + (26 * NbPartPerBoxes) + (NbPartPerBoxes - 1);
while( iter.hasNotFinished() ){
if( sumParticles != iter.data().getDataDown()){
std::cout << "P2P probleme, should be " << sumParticles <<
" iter.data().getDataDown() "<< iter.data().getDataDown() << std::endl;
}
iter.gotoNext();
}
} while(octreeIterator.moveRight());
}
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
......
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