Commit bd07c935 authored by berenger-bramas's avatar berenger-bramas
Browse files

FMM Algorithim for sources != targets model

In this version only what is needed is computed.
So the cell has a trace to know if there are
sources or targets inside them.

Need to clean the code
And rename file/class linked to this system.

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@37 2616d619-271b-44dc-8df4-d4a8f33a7222
parent e9168d2a
......@@ -43,6 +43,28 @@ public:
* @param inPosition the position of the current cell
*/
virtual void setPosition(const F3DPosition& inPosition) = 0;
/** Because the system can run in ToS mode
* a cell has to express if it has sources
* @return true if there are sources particules inside
*/
virtual bool hasSourcesChild() const = 0;
/** Because the system can run in ToS mode
* a cell has to express if it has targets
* @return true if there are targets particules inside
*/
virtual bool hasTargetsChild() const = 0;
/**
* This function make the cell containing sources
*/
virtual void setSourcesChildTrue() = 0;
/**
* This function make the cell containing targets
*/
virtual void setTargetsChildTrue() = 0;
};
......
#ifndef FFMMALGORITHMTOS_HPP
#define FFMMALGORITHMTOS_HPP
// /!\ Please, you must read the license at the bottom of this page
#include "../Utils/FAssertable.hpp"
#include "../Utils/FDebug.hpp"
#include "../Utils/FTrace.hpp"
#include "../Utils/FTic.hpp"
#include "../Containers/FOctree.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FFMMAlgorithmToS
* @brief
* Please read the license
*
* This class is a basic FMM algorithm
* It just iterates on a tree and call the kernels with good arguments.
*
* Of course this class does not deallocate pointer given in arguements.
*/
template<template< class ParticuleClass, class CellClass, int OctreeHeight> class KernelClass,
class ParticuleClass, class CellClass,
template<class ParticuleClass> class LeafClass,
int OctreeHeight, int SubtreeHeight>
class FFMMAlgorithmToS : protected FAssertable{
// To reduce the size of variable type based on foctree in this file
typedef FOctree<ParticuleClass, CellClass, LeafClass, OctreeHeight, SubtreeHeight> Octree;
typedef typename FOctree<ParticuleClass, CellClass,LeafClass, OctreeHeight, SubtreeHeight>::Iterator FOctreeIterator;
Octree* const tree; //< The octree to work on
KernelClass<ParticuleClass, CellClass, OctreeHeight>* const kernels; //< The kernels
FDEBUG(FTic counterTime); //< In case of debug: to count the elapsed time
FDEBUG(FTic computationCounter); //< In case of debug: to count computation time
public:
/** The constructor need the octree and the kernels used for computation
* @param inTree the octree to work on
* @param inKernels the kernels to call
* An assert is launched if one of the arguments is null
*/
FFMMAlgorithmToS(Octree* const inTree, KernelClass<ParticuleClass,CellClass,OctreeHeight>* const inKernels)
: tree(inTree) , kernels(inKernels) {
assert(tree, "tree cannot be null", __LINE__, __FILE__);
assert(kernels, "kernels cannot be null", __LINE__, __FILE__);
FDEBUG(FDebug::Controller << "FFMMAlgorithmToS\n");
}
/** Default destructor */
virtual ~FFMMAlgorithmToS(){
}
/**
* To execute the fmm algorithm
* Call this function to run the complete algorithm
*/
void execute(){
FTRACE( FTrace::Controller.enterFunction(FTrace::FMM, __FUNCTION__ , __FILE__ , __LINE__) );
kernels->init();
bottomPass();
upwardPass();
downardPass();
directPass();
FTRACE( FTrace::Controller.leaveFunction(FTrace::FMM) );
}
/** P2M */
void bottomPass(){
FTRACE( FTrace::Controller.enterFunction(FTrace::FMM, __FUNCTION__ , __FILE__ , __LINE__) );
FDEBUG( FDebug::Controller.write("\tStart Bottom Pass\n").write(FDebug::Flush) );
FDEBUG( counterTime.tic() );
FDEBUG( double totalComputation = 0 );
FOctreeIterator octreeIterator(tree);
// Iterate on leafs
octreeIterator.gotoBottomLeft();
do{
// We need the current cell that represent the leaf
// and the list of particules
FDEBUG(computationCounter.tic());
FList<ParticuleClass*>* const sources = octreeIterator.getCurrentListSources();
if(sources->getSize()){
octreeIterator.getCurrentCell()->setSourcesChildTrue();
kernels->P2M( octreeIterator.getCurrentCell() , sources);
}
if(octreeIterator.getCurrentListTargets()->getSize()){
octreeIterator.getCurrentCell()->setTargetsChildTrue();
}
FDEBUG(computationCounter.tac());
FDEBUG(totalComputation += computationCounter.elapsed());
} while(octreeIterator.moveRight());
FDEBUG( counterTime.tac() );
FDEBUG( FDebug::Controller << "\tFinished (" << counterTime.elapsed() << "s)\n" );
FDEBUG( FDebug::Controller << "\t\t Computation : " << totalComputation << " s\n" );
FTRACE( FTrace::Controller.leaveFunction(FTrace::FMM) );
}
/** M2M */
void upwardPass(){
FTRACE( FTrace::Controller.enterFunction(FTrace::FMM, __FUNCTION__ , __FILE__ , __LINE__) );
FDEBUG( FDebug::Controller.write("\tStart Upward Pass\n").write(FDebug::Flush); );
FDEBUG( counterTime.tic() );
FDEBUG( double totalComputation = 0 );
// Start from leal level - 1
FOctreeIterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
octreeIterator.moveUp();
FOctreeIterator avoidGotoLeftIterator(octreeIterator);
// for each levels
for(int idxLevel = OctreeHeight - 2 ; idxLevel > 1 ; --idxLevel ){
// for each cells
do{
// We need the current cell and the child
// child is an array (of 8 child) that may be null
FDEBUG(computationCounter.tic());
CellClass* potentialChild[8];
CellClass** const realChild = octreeIterator.getCurrentChild();
CellClass* const currentCell = octreeIterator.getCurrentCell();
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
potentialChild[idxChild] = 0;
if(realChild[idxChild]){
if(realChild[idxChild]->hasSourcesChild()){
currentCell->setSourcesChildTrue();
potentialChild[idxChild] = realChild[idxChild];
}
if(realChild[idxChild]->hasTargetsChild()){
currentCell->setTargetsChildTrue();
}
}
}
kernels->M2M( currentCell , potentialChild, idxLevel);
FDEBUG(computationCounter.tac());
FDEBUG(totalComputation += computationCounter.elapsed());
} while(octreeIterator.moveRight());
avoidGotoLeftIterator.moveUp();
octreeIterator = avoidGotoLeftIterator;// equal octreeIterator.moveUp(); octreeIterator.gotoLeft();
}
FDEBUG( counterTime.tac() );
FDEBUG( FDebug::Controller << "\tFinished (" << counterTime.elapsed() << "s)\n" );
FDEBUG( FDebug::Controller << "\t\t Computation : " << totalComputation << " s\n" );
FTRACE( FTrace::Controller.leaveFunction(FTrace::FMM) );
}
/** M2L L2L */
void downardPass(){
FTRACE( FTrace::Controller.enterFunction(FTrace::FMM, __FUNCTION__ , __FILE__ , __LINE__) );
FDEBUG( FDebug::Controller.write("\tStart Downward Pass (M2L)\n").write(FDebug::Flush); );
FDEBUG( counterTime.tic() );
FDEBUG( double totalComputation = 0 );
{ // first M2L
FOctreeIterator octreeIterator(tree);
octreeIterator.moveDown();
FOctreeIterator avoidGotoLeftIterator(octreeIterator);
CellClass* neighbors[208];
// for each levels
for(int idxLevel = 2 ; idxLevel < OctreeHeight ; ++idxLevel ){
// for each cells
do{
FDEBUG(computationCounter.tic());
CellClass* const currentCell = octreeIterator.getCurrentCell();
if(currentCell->hasTargetsChild()){
const int counter = tree->getDistantNeighbors(neighbors, octreeIterator.getCurrentGlobalIndex(),idxLevel);
int offsetTargetNeighbors = 0;
for(int idxRealNeighbors = 0 ; idxRealNeighbors < counter ; ++idxRealNeighbors, ++offsetTargetNeighbors){
if(neighbors[idxRealNeighbors]->hasSourcesChild()){
if(idxRealNeighbors != offsetTargetNeighbors){
neighbors[offsetTargetNeighbors] = neighbors[idxRealNeighbors];
}
}
else{
--offsetTargetNeighbors;
}
}
if(offsetTargetNeighbors){
kernels->M2L( currentCell , neighbors, offsetTargetNeighbors, idxLevel);
}
}
FDEBUG(computationCounter.tac());
FDEBUG(totalComputation += computationCounter.elapsed());
} while(octreeIterator.moveRight());
avoidGotoLeftIterator.moveDown();
octreeIterator = avoidGotoLeftIterator;
}
}
FDEBUG( counterTime.tac() );
FDEBUG( FDebug::Controller << "\tFinished (" << counterTime.elapsed() << "s)\n" );
FDEBUG( FDebug::Controller << "\t\t Computation : " << totalComputation << " s\n" );
FDEBUG( FDebug::Controller.write("\tStart Downward Pass (L2L)\n").write(FDebug::Flush); );
FDEBUG( counterTime.tic() );
FDEBUG( totalComputation = 0 );
{ // second L2L
FOctreeIterator octreeIterator(tree);
octreeIterator.moveDown();
FOctreeIterator avoidGotoLeftIterator(octreeIterator);
const int heightMinusOne = OctreeHeight - 1;
// for each levels exepted leaf level
for(int idxLevel = 2 ; idxLevel < heightMinusOne ; ++idxLevel ){
// for each cells
do{
FDEBUG(computationCounter.tic());
CellClass* potentialChild[8];
CellClass** const realChild = octreeIterator.getCurrentChild();
CellClass* const currentCell = octreeIterator.getCurrentCell();
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
if(realChild[idxChild] && realChild[idxChild]->hasTargetsChild()){
potentialChild[idxChild] = realChild[idxChild];
}
else{
potentialChild[idxChild] = 0;
}
}
kernels->L2L( currentCell , potentialChild, idxLevel);
FDEBUG(computationCounter.tac());
FDEBUG(totalComputation += computationCounter.elapsed());
} while(octreeIterator.moveRight());
avoidGotoLeftIterator.moveDown();
octreeIterator = avoidGotoLeftIterator;
}
}
FDEBUG( counterTime.tac() );
FDEBUG( FDebug::Controller << "\tFinished (" << counterTime.elapsed() << "s)\n" );
FDEBUG( FDebug::Controller << "\t\t Computation : " << totalComputation << " s\n" );
FTRACE( FTrace::Controller.leaveFunction(FTrace::FMM) );
}
/** P2P */
void directPass(){
FTRACE( FTrace::Controller.enterFunction(FTrace::FMM, __FUNCTION__ , __FILE__ , __LINE__) );
FDEBUG( FDebug::Controller.write("\tStart Direct Pass\n").write(FDebug::Flush); );
FDEBUG( counterTime.tic() );
FDEBUG( double totalComputation = 0 );
const int heightMinusOne = OctreeHeight - 1;
FOctreeIterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
// There is a maximum of 26 neighbors
FList<ParticuleClass*>* neighbors[26];
// for each leafs
do{
FDEBUG(computationCounter.tic());
kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
// need the current particules and neighbors particules
const int counter = tree->getLeafsNeighbors(neighbors, octreeIterator.getCurrentGlobalIndex(),heightMinusOne);
kernels->P2P( octreeIterator.getCurrentListTargets(), octreeIterator.getCurrentListSources() , neighbors, counter);
FDEBUG(computationCounter.tac());
FDEBUG(totalComputation += computationCounter.elapsed());
} while(octreeIterator.moveRight());
FDEBUG( counterTime.tac() );
FDEBUG( FDebug::Controller << "\tFinished (" << counterTime.elapsed() << "s)\n" );
FDEBUG( FDebug::Controller << "\t\t Computation : " << totalComputation << " s\n" );
FTRACE( FTrace::Controller.leaveFunction(FTrace::FMM) );
}
};
#endif //FFMMALGORITHMTOS_HPP
// [--LICENSE--]
......@@ -108,7 +108,7 @@ void ValidateFMMAlgo(FOctree<ParticuleClass, CellClass, LeafClass, TreeHeight ,
std::cout << "Check Result\n";
int NbPart = 0;
{ // Check that each particule has been summed with all other
typename FOctree<FTestParticule, FTestCell, LeafClass, TreeHeight, SizeSubLevels>::Iterator octreeIterator(tree);
typename FOctree<ParticuleClass, CellClass, LeafClass, TreeHeight, SizeSubLevels>::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
do{
if(octreeIterator.getCurrentCell()->getDataUp() != octreeIterator.getCurrentListSources()->getSize() ){
......@@ -118,7 +118,7 @@ void ValidateFMMAlgo(FOctree<ParticuleClass, CellClass, LeafClass, TreeHeight ,
} while(octreeIterator.moveRight());
}
{ // Ceck if there is number of NbPart summed at level 1
typename FOctree<FTestParticule, FTestCell, LeafClass, TreeHeight, SizeSubLevels>::Iterator octreeIterator(tree);
typename FOctree<ParticuleClass, CellClass, LeafClass, TreeHeight, SizeSubLevels>::Iterator octreeIterator(tree);
octreeIterator.moveDown();
long res = 0;
do{
......@@ -129,7 +129,7 @@ void ValidateFMMAlgo(FOctree<ParticuleClass, CellClass, LeafClass, TreeHeight ,
}
}
{ // Ceck if there is number of NbPart summed at level 1
typename FOctree<FTestParticule, FTestCell, LeafClass, TreeHeight, SizeSubLevels>::Iterator octreeIterator(tree);
typename FOctree<ParticuleClass, CellClass, LeafClass, TreeHeight, SizeSubLevels>::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
for(int idxLevel = TreeHeight - 1 ; idxLevel > 1 ; --idxLevel ){
long res = 0;
......@@ -144,14 +144,18 @@ void ValidateFMMAlgo(FOctree<ParticuleClass, CellClass, LeafClass, TreeHeight ,
}
}
{ // Check that each particule has been summed with all other
typename FOctree<FTestParticule, FTestCell, LeafClass, TreeHeight, SizeSubLevels>::Iterator octreeIterator(tree);
typename FOctree<ParticuleClass, CellClass, LeafClass, TreeHeight, SizeSubLevels>::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
do{
typename FList<FTestParticule*>::BasicIterator iter(*octreeIterator.getCurrentListTargets());
typename FList<ParticuleClass*>::BasicIterator iter(*octreeIterator.getCurrentListTargets());
const bool isUsingToS = (octreeIterator.getCurrentListTargets() != octreeIterator.getCurrentListSources());
while( iter.isValide() ){
// If a particules has been impacted by less than NbPart - 1 (the current particule)
// there is a problem
if(iter.value()->getDataDown() != NbPart - 1){
if( (!isUsingToS && iter.value()->getDataDown() != NbPart - 1) ||
(isUsingToS && iter.value()->getDataDown() != NbPart) ){
std::cout << "Problem L2P + P2P : " << iter.value()->getDataDown() << "\n";
}
iter.progress();
......
#ifndef FEXTENDCELLTYPE_HPP
#define FEXTENDCELLTYPE_HPP
// /!\ Please, you must read the license at the bottom of this page
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FExtendCellType
* Please read the license
* This class is an extenssion.
* It proposes a target/source extenssion for cell.
* Because cells may have child that contains only
* sources or targets (in ToS system) then it is important
* to not compute for nothing.
*/
class FExtendCellType {
protected:
/** Particule potential type */
static const int Neither = 0;
static const int ContainsSources = 1;
static const int ContainsTargets = 2;
/** Current type */
int type;
public:
/** Default constructor */
FExtendCellType() : type(Neither) {
}
/** Copy constructor */
FExtendCellType(const FExtendCellType& other) : type(other.type) {
}
/** Destructor */
virtual ~FExtendCellType(){
}
/** Copy operator */
FExtendCellType& operator=(const FExtendCellType& other) {
this->type = other.type;
return *this;
}
/** To know if a cell has sources */
bool hasSourcesChild() const {
return this->type & ContainsSources;
}
/** To know if a cell has targets */
bool hasTargetsChild() const {
return this->type & ContainsTargets;
}
/** To set cell as sources container */
void setSourcesChildTrue() {
this->type |= ContainsSources;
}
/** To set cell as targets container */
void setTargetsChildTrue() {
this->type |= ContainsTargets;
}
};
#endif //FEXTENDCELLTYPE_HPP
// [--LICENSE--]
......@@ -115,6 +115,7 @@ public:
inParticule->setPosition(x,y,z);
inParticule->setPhysicalValue(data);
if(isTarget) inParticule->setAsTarget();
else inParticule->setAsSource();
}
};
......
......@@ -84,7 +84,7 @@ int main(int argc, char ** argv){
// FTestKernels FBasicKernels
FTestKernels<FTestParticule, FTestCell, NbLevels> kernels;
//FFMMAlgorithm FFMMAlgorithmThreaded FFMMAlgorithmArray FFMMAlgorithmTask
FFMMAlgorithmArray<FTestKernels, FTestParticule, FTestCell, FSimpleLeaf, NbLevels, SizeSubLevels> algo(&tree,&kernels);
FFMMAlgorithm<FTestKernels, FTestParticule, FTestCell, FSimpleLeaf, NbLevels, SizeSubLevels> algo(&tree,&kernels);
algo.execute();
counter.tac();
......
// /!\ Please, you must read the license at the bottom of this page
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include "../Sources/Utils/FTic.hpp"
#include "../Sources/Containers/FOctree.hpp"
#include "../Sources/Containers/FList.hpp"
#include "../Sources/Core/FTypedLeaf.hpp"
#include "../Sources/Utils/F3DPosition.hpp"
#include "../Sources/Core/FTestParticule.hpp"
#include "../Sources/Core/FTestCell.hpp"
#include "../Sources/Core/FTestKernels.hpp"
#include "../Sources/Extenssions/FExtendParticuleType.hpp"
#include "../Sources/Extenssions/FExtendCellType.hpp"
#include "../Sources/Core/FFMMAlgorithmToS.hpp"
#include "../Sources/Core/FFMMAlgorithm.hpp"
#include "../Sources/Core/FFMMAlgorithmArray.hpp"
#include "../Sources/Core/FFMMAlgorithmThreaded.hpp"
#include "../Sources/Core/FFMMAlgorithmTask.hpp"
#include "../Sources/Core/FBasicKernels.hpp"
/** This program show an example of use of
* the fmm basic algo
* it also check that each particules is impacted each other particules
*/
class FTestParticuleToS : public FTestParticule, public FExtendParticuleType {
};
class FTestCellToS: public FTestCell , public FExtendCellType{
};
// Simply create particules and try the kernels
int main(int argc, char ** argv){
///////////////////////What we do/////////////////////////////
std::cout << ">> This executable has to be used to test the FMM algorithm.\n";
//////////////////////////////////////////////////////////////
const int NbLevels = 10;//10;
const int SizeSubLevels = 3;//3
const long NbPart = 2000000;//2000000
FTestParticuleToS* particules = new FTestParticuleToS[NbPart];
FTic counter;
srand ( 1 ); // volontary set seed to constant
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Creating " << NbPart << " particules ..." << std::endl;
counter.tic();
for(long idxPart = 0 ; idxPart < NbPart ; ++idxPart){
particules[idxPart].setPosition(FReal(rand())/RAND_MAX,FReal(rand())/RAND_MAX,FReal(rand())/RAND_MAX);
if(rand() > RAND_MAX/2) particules[idxPart].setAsTarget();
else particules[idxPart].setAsSource();
}
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
FOctree<FTestParticuleToS, FTestCellToS, FTypedLeaf, NbLevels, SizeSubLevels> tree(1.0,F3DPosition(0.5,0.5,0.5));
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Inserting particules ..." << std::endl;
counter.tic();
for(long idxPart = 0 ; idxPart < NbPart ; ++idxPart){
tree.insert(&particules[idxPart]);
}
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Working on particules ..." << std::endl;
counter.tic();
// FTestKernels FBasicKernels
FTestKernels<FTestParticuleToS, FTestCellToS, NbLevels> kernels;
//FFMMAlgorithm FFMMAlgorithmThreaded FFMMAlgorithmArray FFMMAlgorithmTask
FFMMAlgorithm<FTestKernels, FTestParticuleToS, FTestCellToS, FTypedLeaf, NbLevels, SizeSubLevels> algo(&tree,&kernels);
algo.execute();
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
ValidateFMMAlgo(&tree);
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Deleting particules ..." << std::endl;
counter.tic();
for(long idxPart = 0 ; idxPart < NbPart ; ++idxPart){
particules[idxPart].~FTestParticuleToS();
}
delete [] particules;
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
return 0;
}
// [--LICENSE--]
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