Commit 6f55fa9a authored by BRAMAS Berenger's avatar BRAMAS Berenger

Create the periodic system using normal FMM functions for all kernels, add...

Create the periodic system using normal FMM functions for all kernels, add unit test for all 3 kernels, be aware that current particles are randomly generated and not real system for periodic simulations
parent 74ca87da
This diff is collapsed.
This diff is collapsed.
......@@ -61,7 +61,8 @@ class FFmmAlgorithmThreadProcPeriodic : protected FAssertable {
const FMpi::FComm& comm; //< MPI comm
CellClass root; //< root of tree needed by the periodicity
int periodicLimit; //< the upper periodic limit
int nbLevelAboveRoot; //< the upper periodic limit
int offsetHeight;
typename OctreeClass::Iterator* iterArray;
int numberOfLeafs; //< To store the size at the previous level
......@@ -89,6 +90,33 @@ class FFmmAlgorithmThreadProcPeriodic : protected FAssertable {
public:
/** To know how many times the box is repeated in each direction
* -x +x -y +y -z +z
*/
int repeatedBox() const {
return nbLevelAboveRoot == -1 ? 3 : 7 * (1<<(nbLevelAboveRoot));
}
/** This function has to be used to init the kernel with correct args
* it return the box seen from a kernel point of view from the periodicity the user ask for
* this is computed using the originalBoxWidth given in parameter
* @param originalBoxWidth the real system size
* @return the size the kernel should use
*/
int boxwidthForKernel(const FReal originalBoxWidth) const {
return originalBoxWidth * ((1<<(offsetHeight+1)) - 1);
}
/** This function has to be used to init the kernel with correct args
* it return the tree heigh seen from a kernel point of view from the periodicity the user ask for
* this is computed using the originalTreeHeight given in parameter
* @param originalTreeHeight the real tree heigh
* @return the heigh the kernel should use
*/
int treeHeightForKernel(const FReal originalTreeHeight) const {
return OctreeHeight + offsetHeight;
}
Interval& getWorkingInterval(const int level){
return getWorkingInterval(level, idProcess);
}
......@@ -102,13 +130,14 @@ public:
* @param inKernels the kernels to call
* An assert is launched if one of the arguments is null
*/
FFmmAlgorithmThreadProcPeriodic(const FMpi::FComm& inComm, OctreeClass* const inTree, KernelClass* const inKernels, const int inPeriodicLimit = 5)
: tree(inTree) , kernels(0), comm(inComm), periodicLimit(inPeriodicLimit), numberOfLeafs(0),
FFmmAlgorithmThreadProcPeriodic(const FMpi::FComm& inComm, OctreeClass* const inTree, KernelClass* const inKernels, const int inUpperLevel = 5)
: tree(inTree) , kernels(0), comm(inComm), nbLevelAboveRoot(inUpperLevel), offsetHeight(inUpperLevel+2), numberOfLeafs(0),
MaxThreads(omp_get_max_threads()), nbProcess(inComm.processCount()), idProcess(inComm.processId()),
OctreeHeight(tree->getHeight()),intervals(new Interval[inComm.processCount()]),
workingIntervalsPerLevel(new Interval[inComm.processCount() * tree->getHeight()]) {
fassert(tree, "tree cannot be null", __LINE__, __FILE__);
fassert(-1 <= inUpperLevel, "inUpperLevel cannot be < -1", __LINE__, __FILE__);
this->kernels = new KernelClass*[MaxThreads];
for(int idxThread = 0 ; idxThread < MaxThreads ; ++idxThread){
......@@ -277,7 +306,7 @@ private:
int firstProcThatSend = idProcess + 1;
// for each levels
for(int idxLevel = OctreeHeight - 2 ; idxLevel >= 1 ; --idxLevel ){
for(int idxLevel = OctreeHeight - 2 ; idxLevel > 0 ; --idxLevel ){
// No more work for me
if(idProcess != 0
&& getWorkingInterval((idxLevel+1), idProcess).max <= getWorkingInterval((idxLevel+1), idProcess - 1).max){
......@@ -1430,62 +1459,78 @@ private:
/** Periodicity */
void processPeriodicLevels(){
if( !periodicLimit ){
return;
}
// If nb level == -1 nothing to do
if( nbLevelAboveRoot == -1 ){
return;
}
CellClass upperCells[periodicLimit];
upperCells[0] = root;
// in other situation, we have to compute M2L a level 1
// but also
// the cell at level >= 0
CellClass upperCells[nbLevelAboveRoot+1];
upperCells[nbLevelAboveRoot] = root;
// Then M2M from level 0 to level -LIMITE
{
CellClass* virtualChild[8];
for(int idxLevel = 1 ; idxLevel < periodicLimit ; ++idxLevel){
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
virtualChild[idxChild] = &upperCells[idxLevel-1];
// Then M2M from level 0 to level -LIMITE
{
CellClass* virtualChild[8];
for(int idxLevel = nbLevelAboveRoot - 1 ; idxLevel >= 0 ; --idxLevel){
FMemUtils::setall(virtualChild, &upperCells[idxLevel+1], 8);
kernels[0]->M2M( &upperCells[idxLevel], virtualChild, idxLevel + 2);
}
kernels[0]->M2M( &upperCells[idxLevel], virtualChild, -idxLevel);
}
}
// Then M2L at all level
{
// We say that we are in the child index 0
// So we can compute one time the relative indexes
const CellClass* neighbors[343];
memset(neighbors, 0, sizeof(CellClass*) * 343);
for(int idxX = -2 ; idxX <= 3 ; ++idxX){
for(int idxY = -2 ; idxY <= 3 ; ++idxY){
for(int idxZ = -2 ; idxZ <= 3 ; ++idxZ){
if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
neighbors[(((idxX+3)*7) + (idxY+3))*7 + (idxZ + 3)] = reinterpret_cast<const CellClass*>(~0);
// Then M2L at all level
{
// We say that we are in the child index 0
// So we can compute one time the relative indexes
const CellClass* neighbors[343];
memset(neighbors, 0, sizeof(CellClass*) * 343);
for(int idxX = -2 ; idxX <= 3 ; ++idxX){
for(int idxY = -2 ; idxY <= 3 ; ++idxY){
for(int idxZ = -2 ; idxZ <= 3 ; ++idxZ){
if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
neighbors[(((idxX+3)*7) + (idxY+3))*7 + (idxZ + 3)] = reinterpret_cast<const CellClass*>(~0);
}
}
}
}
}
const int counter = 189;
for(int idxLevel = 0 ; idxLevel < periodicLimit ; ++idxLevel ){
for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
if(neighbors[idxNeigh]){
neighbors[idxNeigh] = &upperCells[idxLevel];
for(int idxLevel = nbLevelAboveRoot ; idxLevel > 0 ; --idxLevel ){
for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
if(neighbors[idxNeigh]){
neighbors[idxNeigh] = &upperCells[idxLevel];
}
}
kernels[0]->M2L( &upperCells[idxLevel] , neighbors, 189, idxLevel + 2);
}
// We had the last level border
memset(neighbors, 0, sizeof(CellClass*) * 343);
for(int idxX = -3 ; idxX <= 3 ; ++idxX){
for(int idxY = -3 ; idxY <= 3 ; ++idxY){
for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ){
if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
neighbors[(((idxX+3)*7) + (idxY+3))*7 + (idxZ + 3)] = &upperCells[0];
}
}
}
}
kernels[0]->M2L( &upperCells[idxLevel] , neighbors, counter, -idxLevel);
kernels->M2L( &upperCells[0] , neighbors, 343-27, 2);
}
}
// Finally L2L until level 0
{
CellClass* virtualChild[8];
memset(virtualChild, 0, sizeof(CellClass*) * 8);
for(int idxLevel = periodicLimit - 1 ; idxLevel > 0 ; --idxLevel){
virtualChild[0] = &upperCells[idxLevel-1];
kernels[0]->L2L( &upperCells[idxLevel], virtualChild, -idxLevel);
// Finally L2L until level 0
{
CellClass* virtualChild[8];
memset(virtualChild, 0, sizeof(CellClass*) * 8);
for(int idxLevel = 0 ; idxLevel < nbLevelAboveRoot ; ++idxLevel){
virtualChild[0] = &upperCells[idxLevel+1];
kernels[0]->L2L( &upperCells[idxLevel], virtualChild, idxLevel + 2);
}
}
}
root = upperCells[0];
}
root = upperCells[nbLevelAboveRoot];
}
};
......
#ifndef FGLOBALPERIODIC_HPP
#define FGLOBALPERIODIC_HPP
///////////////////////////////////////////////////////
// Periodic condition definition
///////////////////////////////////////////////////////
enum PeriodicCondition {
DirPlusX = 1 << 0,
DirMinusX = 1 << 1,
DirPlusY = 1 << 2,
DirMinusY = 1 << 3,
DirPlusZ = 1 << 4,
DirMinusZ = 1 << 5,
DirX = (DirPlusX | DirMinusX),
DirY = (DirPlusY | DirMinusY),
DirZ = (DirPlusZ | DirMinusZ),
AllDirs = (DirX | DirY | DirZ)
};
bool testPeriodicCondition(const int conditions, const PeriodicCondition testConditions) {
return (conditions & testConditions) == testConditions;
}
#endif // FGLOBALPERIODIC_HPP
......@@ -86,8 +86,8 @@ int main(int argc, char ** argv){
const int NbLevels = FParameters::getValue(argc,argv,"-h", 4);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 2);
const int DevP = FParameters::getValue(argc,argv,"-P", 5);
const int BoundaryDeep = FParameters::getValue(argc,argv,"-bd", 5);
const int DevP = FParameters::getValue(argc,argv,"-P", 9);
const int PeriodicDeep = FParameters::getValue(argc,argv,"-per", 2);
const char* const filename = FParameters::getStr(argc,argv,"-f", "../Data/testEwal417.dt");
FTic counter;
......@@ -119,6 +119,7 @@ int main(int argc, char ** argv){
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
loader.fillParticle(particles[idxPart]);
// reset forces and insert in the tree
particles[idxPart].setIndex(idxPart);
ParticleClass part = particles[idxPart];
part.setForces(0,0,0);
part.setPotential(0);
......@@ -130,46 +131,61 @@ int main(int argc, char ** argv){
// -----------------------------------------------------
std::cout << "Create kernel ..." << std::endl;
std::cout << "Create kernel & run simu ..." << std::endl;
counter.tic();
KernelClass kernels(DevP, NbLevels,loader.getBoxWidth(),loader.getCenterOfBox(), BoundaryDeep);
counter.tac();
std::cout << "Done " << " in " << counter.elapsed() << "s)." << std::endl;
// -----------------------------------------------------
std::cout << "Working on particles ..." << std::endl;
FmmClass algo(&tree,&kernels,BoundaryDeep);
counter.tic();
FmmClass algo(&tree,PeriodicDeep);
KernelClass kernels( DevP, algo.treeHeightForKernel(),
algo.boxwidthForKernel(loader.getBoxWidth()), loader.getCenterOfBox());
algo.setKernel(&kernels);
algo.execute();
counter.tac();
std::cout << "Done " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl;
// -----------------------------------------------------
{
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
ContainerClass::ConstBasicIterator iter(*octreeIterator.getCurrentListTargets());
while( iter.hasNotFinished() ){
std::cout << ">> index " << iter.data().getIndex() << " type " << iter.data().getType() << std::endl;
std::cout << "x " << iter.data().getPosition().getX() << " y " << iter.data().getPosition().getY() << " z " << iter.data().getPosition().getZ() << std::endl;
std::cout << "fx " << iter.data().getForces().getX() << " fy " << iter.data().getForces().getY() << " fz " << iter.data().getForces().getZ() << std::endl;
std::cout << "physical value " << iter.data().getPhysicalValue() << " potential " << iter.data().getPotential() << std::endl;
const ParticleClass& part = particles[iter.data().getIndex()];
std::cout << "x " << part.getPosition().getX() << " y " << part.getPosition().getY() << " z " << part.getPosition().getZ() << std::endl;
std::cout << "fx " <<part.getForces().getX() << " fy " << part.getForces().getY() << " fz " << part.getForces().getZ() << std::endl;
std::cout << "physical value " << part.getPhysicalValue() << " potential " << part.getPotential() << std::endl;
std::cout << ">> index " << iter.data().getIndex() << " type " << iter.data().getType() << std::endl;
std::cout << "Good x " << part.getPosition().getX() << " y " << part.getPosition().getY() << " z " << part.getPosition().getZ() << std::endl;
std::cout << "FMM x " << iter.data().getPosition().getX() << " y " << iter.data().getPosition().getY() << " z " << iter.data().getPosition().getZ() << std::endl;
std::cout << "Good fx " <<part.getForces().getX() << " fy " << part.getForces().getY() << " fz " << part.getForces().getZ() << std::endl;
std::cout << "FMM fx " << iter.data().getForces().getX() << " fy " << iter.data().getForces().getY() << " fz " << iter.data().getForces().getZ() << std::endl;
std::cout << "GOOD physical value " << part.getPhysicalValue() << " potential " << part.getPotential() << std::endl;
std::cout << "FMM physical value " << iter.data().getPhysicalValue() << " potential " << iter.data().getPotential() << std::endl;
std::cout << "\n";
potentialDiff.add(part.getPotential(),iter.data().getPotential());
fx.add(part.getForces().getX(),iter.data().getForces().getX());
fy.add(part.getForces().getY(),iter.data().getForces().getY());
fz.add(part.getForces().getZ(),iter.data().getForces().getZ());
iter.gotoNext();
}
} while(octreeIterator.moveRight());
printf("Potential diff is = \n");
printf("%f\n",potentialDiff.getL2Norm());
printf("%f\n",potentialDiff.getInfNorm());
printf("Fx diff is = \n");
printf("%f\n",fx.getL2Norm());
printf("%f\n",fx.getInfNorm());
printf("Fy diff is = \n");
printf("%f\n",fy.getL2Norm());
printf("%f\n",fy.getInfNorm());
printf("Fz diff is = \n");
printf("%f\n",fz.getL2Norm());
printf("%f\n",fz.getInfNorm());
}
// -----------------------------------------------------
......
......@@ -18,6 +18,8 @@
#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Utils/FTic.hpp"
#include "../../Src/Files/FRandomLoader.hpp"
#include "../../Src/Containers/FOctree.hpp"
#include "../../Src/Containers/FVector.hpp"
......@@ -54,50 +56,29 @@ int main(int argc, char ** argv){
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);
const long NbParticles = FParameters::getValue(argc,argv,"-nb", 1000);
const int PeriodicDeep = FParameters::getValue(argc,argv,"-per", 2);
// choose in +x dir or -/+x dir or all dirs
int PeriodicDirs = (FParameters::existParameter(argc,argv,"-x")?DirMinusX:0) |
(FParameters::existParameter(argc,argv,"+x")?DirPlusX:0) |
(FParameters::existParameter(argc,argv,"-y")?DirMinusY:0) |
(FParameters::existParameter(argc,argv,"+y")?DirPlusY:0) |
(FParameters::existParameter(argc,argv,"-z")?DirMinusZ:0) |
(FParameters::existParameter(argc,argv,"+z")?DirPlusZ:0);
if( PeriodicDirs == 0 ) PeriodicDirs = AllDirs;
FTic counter;
srand ( 1 ); // volontary set seed to constant
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
const FReal BoxWidth = 1.0;
const FReal BoxCenter = 0.5;
OctreeClass tree(NbLevels, SizeSubLevels, BoxWidth, FPoint(BoxCenter,BoxCenter,BoxCenter));
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Creating & Inserting " << NbPartPerBoxes << " particles per boxes ..." << std::endl;
std::cout << "Creating & Inserting " << NbParticles << " particles ..." << 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 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);
}
}
}
}
}
FRandomLoader<ParticleClass> loader(NbParticles);
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
loader.fillTree(tree);
counter.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << "s)." << std::endl;
......@@ -108,11 +89,9 @@ int main(int argc, char ** argv){
std::cout << "Working on particles ..." << std::endl;
counter.tic();
const int PeriodicDeep = 4;
// FTestKernels FBasicKernels
KernelClass kernels;
//FFmmAlgorithm FFmmAlgorithmThread
FmmClass algo( &tree, &kernels, PeriodicDeep);
FmmClass algo( &tree, PeriodicDeep, PeriodicDirs);
algo.setKernel(&kernels);
algo.execute();
counter.tac();
......@@ -121,8 +100,6 @@ int main(int argc, char ** argv){
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
//ValidateFMMAlgo<OctreeClass, ParticleClass, CellClass, ContainerClass, LeafClass>(&tree);
{ // Check that each particle has been summed with all other
long long int counterNbPart = 0;
......@@ -130,92 +107,34 @@ int main(int argc, char ** argv){
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() <<
if(octreeIterator.getCurrentCell()->getDataUp() != octreeIterator.getCurrentListSrc()->getSize() ){
std::cout << "Problem P2M 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;
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( counterNbPart != NbParticles){
std::cout << "Problem global nb part, counter = " << counterNbPart << " created = " << NbParticles << std::endl;
}
}
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;
}
{
const FTreeCoordinate repetitions = algo.repetitions();
const int totalRepeatedBox = repetitions.getX() * repetitions.getY() * repetitions.getZ();
std::cout << "The box is repeated " << repetitions.getX() <<" "<< repetitions.getY()<<" "<<
repetitions.getZ() << " there are " << totalRepeatedBox << " boxes in total\n";
const long long NbParticlesEntireSystem = NbParticles * totalRepeatedBox;
std::cout << "The total number of particles is " << NbParticlesEntireSystem << "\n";
for( int idxLevel = PeriodicDeep - 1 ; idxLevel >= 0 ; --idxLevel ){
counterL2L += counterM2L[idxLevel];
}
}
{
long long int particlesPerBox = NbPartPerBoxes * FMath::pow(8,NbLevels-2);
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
OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
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 <<
if( NbParticlesEntireSystem - 1 != iter.data().getDataDown()){
std::cout << "P2P probleme, should be " << NbParticlesEntireSystem - 1 <<
" iter.data().getDataDown() "<< iter.data().getDataDown() << std::endl;
}
......
......@@ -25,6 +25,7 @@
#include "../../Src/Containers/FVector.hpp"
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Files/FRandomLoader.hpp"
#include "../../Src/Utils/FPoint.hpp"
......@@ -62,7 +63,8 @@ int main(int argc, char ** argv){
const int NbLevels = FParameters::getValue(argc,argv,"-h", 7);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 3);
const int NbPartPerBoxesPerProc = FParameters::getValue(argc,argv,"-nb", 1);
const long NbParticles = FParameters::getValue(argc,argv,"-nb", 1000);
const int PeriodicDeep = FParameters::getValue(argc,argv,"-per", 2);
FMpi app(argc, argv);
......@@ -70,16 +72,6 @@ int main(int argc, char ** argv){
FTic counter;
srand ( 1 ); // volontary set seed to constant
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
const FReal BoxWidth = 1.0;
const FReal BoxCenter = 0.5;
OctreeClass tree(NbLevels, SizeSubLevels, BoxWidth, FPoint(BoxCenter,BoxCenter,BoxCenter));
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
......@@ -87,27 +79,14 @@ int main(int argc, char ** argv){
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
counter.tic();
int NbPart = 0;
FRandomLoader<ParticleClass> loader(NbParticles);
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
{
const int NbSmallBoxesPerSide = (1 << (NbLevels-1));
const FReal SmallBoxWidth = BoxWidth / FReal(NbSmallBoxesPerSide);
const FReal SmallBoxWidthDiv2 = SmallBoxWidth / 2;
NbPart = NbSmallBoxesPerSide * NbSmallBoxesPerSide * NbSmallBoxesPerSide * NbPartPerBoxesPerProc;
FTestParticle*const particles = new FTestParticle[NbPart];
int indexPart = 0;
for(int idxX = 0 ; idxX < NbSmallBoxesPerSide ; ++idxX){
for(int idxY = 0 ; idxY < NbSmallBoxesPerSide ; ++idxY){
for(int idxZ = 0 ; idxZ < NbSmallBoxesPerSide ; ++idxZ){
for(int idxPart = 0 ; idxPart < NbPartPerBoxesPerProc ; ++idxPart ){
particles[indexPart++].setPosition(FReal(idxX)*SmallBoxWidth + SmallBoxWidthDiv2,
FReal(idxY)*SmallBoxWidth + SmallBoxWidthDiv2,
FReal(idxZ)*SmallBoxWidth + SmallBoxWidthDiv2);
}
}
}
FTestParticle*const particles = new FTestParticle[NbParticles];
for(int idx = 0 ; idx < NbParticles ; ++idx){
loader.fillParticle(particles[idx]);
}
FMpiTreeBuilder<ParticleClass>::ArrayToTree(app.global(), particles, NbPart, FPoint(BoxCenter,BoxCenter,BoxCenter), BoxWidth, tree);
......@@ -123,10 +102,7 @@ int main(int argc, char ** argv){
std::cout << "Working on particles ..." << std::endl;
counter.tic();
const int PeriodicDeep = 3;
// FTestKernels FBasicKernels
KernelClass kernels;
//FFmmAlgorithm FFmmAlgorithmThread
FmmClass algo( app.global(), &tree, &kernels, PeriodicDeep);
algo.execute();
......@@ -136,123 +112,21 @@ int main(int argc, char ** argv){
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
{ // Check that each particle has been summed with all other
OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
if( octreeIterator.getCurrentListTargets()->getSize() != NbPartPerBoxes){
std::cout << "Incorrect number of particles per leaf, should be " << NbPartPerBoxes <<
" iter.data().getDataDown() "<< octreeIterator.getCurrentListTargets()->getSize() << std::endl;
}
} while(octreeIterator.moveRight());
}
{ // Check that each particle has been summed with all other
long long int counterNbPart = 0;
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());
const long long int allPart = app.global().reduceSum( counterNbPart );
if( app.global().processId() == 0 && allPart != NbPart * app.global().processCount()){
std::cout << "Problem global nb part, counter = " << allPart << " created = " << NbPart * app.global().processCount() << std::endl;
}
}
{ // Ceck if there is number of NbPart summed at level 1
long long particlesPerBox = NbPartPerBoxes;
OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
for(int idxLevel = NbLevels - 1 ; idxLevel >= 1 ; --idxLevel ){
if(algo.hasWorkAtLevel(idxLevel)){
while(octreeIterator.getCurrentGlobalIndex() != algo.getWorkingInterval(idxLevel).min){
octreeIterator.moveRight();
}
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;
}
}
{
long long int counterL2L = 0;
if( PeriodicDeep ){
long long int particlesPerBox = NbPartPerBoxes * FMath::pow(8,NbLevels-1);
long long int counterUp[PeriodicDeep];
for( int idxLevel = 0 ; idxLevel < PeriodicDeep ; ++idxLevel ){