Commit 987cee7d authored by berenger-bramas's avatar berenger-bramas
Browse files

Easy "P" cutomization

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@24 2616d619-271b-44dc-8df4-d4a8f33a7222
parent 3bac8835
......@@ -2,6 +2,8 @@
#define FLIST_HPP
// /!\ Please, you must read the license at the bottom of this page
#include "../Utils/FGlobal.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FList
......
......@@ -16,6 +16,9 @@
#include <iostream>
// P is a input parameter
static const int FMB_Info_P = 3;
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FAbstractFmbKernels
......@@ -36,8 +39,6 @@
template< class ParticuleClass, class CellClass, int TreeHeight>
class FAbstractFmbKernels : public FAbstractKernels<ParticuleClass,CellClass, TreeHeight> {
protected:
// P is a input parameter
static const int FMB_Info_P = 2;
// _GRAVITATIONAL_
static const int FMB_Info_eps_soft_square = 1;
......@@ -284,25 +285,29 @@ protected:
ptrCosSin->setReal( FMath::Sin(angle + FMath::FPiDiv2) );
ptrCosSin->setImag( FMath::Sin(angle) );
//std::cout<< idxl << "=" << ptrCosSin->getReal() << "/" << ptrCosSin->getImag() << " (" << idxl << "/" << inSphere.phi << "/" << PiArray[lMod4] << ")\n";
//printf("%d=%f/%f (%d/%f/%f)\n",idxl,ptrCosSin->getReal(),ptrCosSin->getImag(),idxl,inSphere.phi,this->PiArrayInner[idxlMod4]);
}
legendreFunction(inSphere.cosTheta, inSphere.sinTheta, this->legendre);
//printf("FMB_Info_M2L_exp_size=%d\n",FMB_Info_M2L_exp_size);
//for(int temp = 0 ; temp < FMB_Info_M2L_exp_size ; ++temp){
// printf("%f\n",this->legendre[temp]);
//}
/*printf("FMB_Info_M2L_exp_size=%d\n",FMB_Info_M2L_exp_size);
for(int temp = 0 ; temp < FMB_Info_M2L_exp_size ; ++temp){
printf("%f\n",this->legendre[temp]);
}*/
FComplexe* currentResult = outResults;
int idxLegendre = 0;//ptr_associated_Legendre_function_Array
int idxSphereHarmoCoef = 0;
FReal idxRl = 1.0 ;
//printf("lmax = %d\n",LMax);
for(int idxl = 0; idxl <= this->LMax ; ++idxl, idxRl *= inSphere.r){
for(int idxm = 0 ; idxm <= idxl ; ++idxm, ++currentResult, ++idxSphereHarmoCoef, ++idxLegendre){
const FReal magnitude = this->sphereHarmoInnerCoef[idxSphereHarmoCoef] * idxRl * legendre[idxLegendre];
currentResult->setReal( magnitude * this->cosSin[idxm].getReal() );
currentResult->setImag( magnitude * this->cosSin[idxm].getImag() );
//printf("l = %d m = %d\n",idxl,idxm);
//printf("magnitude=%f idxRl=%f sphereHarmoInnerCoef=%f real=%f imag=%f\n",magnitude,idxRl,this->sphereHarmoInnerCoef[idxSphereHarmoCoef],currentResult->getReal(),currentResult->getImag());
}
}
......@@ -504,7 +509,7 @@ protected:
void precomputeM2M(){
FReal treeWidthAtLevel = this->treeWidthAtRoot/2;
for(int idxLevel = 0 ; idxLevel < TreeHeight ; ++idxLevel ){
for(int idxLevel = 0 ; idxLevel < TreeHeight - 1 ; ++idxLevel ){
const F3DPosition father(treeWidthAtLevel,treeWidthAtLevel,treeWidthAtLevel);
treeWidthAtLevel /= 2;
......@@ -533,11 +538,11 @@ protected:
positionToSphere(L2LVector,&sphericalL2L);
harmonicInner(sphericalL2L,this->transitionL2L[idxLevel][idxChild]);
//std::cout << "[M2M_vector]" << idxLevel << "/" << idxChild << " = " << M2MVector.getX() << "/" << M2MVector.getY() << "/" << M2MVector.getZ() << "\n";
//std::cout << "[M2M_vectorSpherical]" << idxLevel << "/" << idxChild << " = " << sphericalM2M.r << "/" << sphericalM2M.cosTheta << "/" << sphericalM2M.sinTheta << "/" << sphericalM2M.phi << "\n";
//printf("[M2M_vector]%d/%d = %f/%f/%f\n", idxLevel , idxChild , M2MVector.getX() , M2MVector.getY() , M2MVector.getZ() );
//printf("[M2M_vectorSpherical]%d/%d = %f/%f/%f/%f\n", idxLevel , idxChild , sphericalM2M.r , sphericalM2M.cosTheta , sphericalM2M.sinTheta , sphericalM2M.phi );
//for(int idxExpSize = 0 ; idxExpSize < FMB_Info_exp_size ; ++idxExpSize){
//std::cout << "transitionL2L[" << idxLevel << "][" << idxChild << "][" << idxExpSize << "]=" << this->transitionL2L[idxLevel][idxChild][idxExpSize].getReal()<<"/"<<this->transitionL2L[idxLevel][idxChild][idxExpSize].getImag()<< "\n";
//std::cout << "transitionM2M[" << idxLevel << "][" << idxChild << "][" << idxExpSize << "]=" << this->transitionM2M[idxLevel][idxChild][idxExpSize].getReal()<<"/"<<this->transitionM2M[idxLevel][idxChild][idxExpSize].getImag()<< "\n";
//printf("transitionM2M[%d][%d][%d]=%f/%f\n", idxLevel , idxChild , idxExpSize , this->transitionM2M[idxLevel][idxChild][idxExpSize].getReal(),this->transitionM2M[idxLevel][idxChild][idxExpSize].getImag());
//}
}
}
......
......@@ -6,6 +6,8 @@
#include "../Utils/FComplexe.hpp"
#include "FAbstractFmbKernels.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class FExtendFmbCell
......@@ -14,10 +16,9 @@
* This class is an extenssion.
* It is needed by the Fmb Kernels.
*/
template <int P>
class FExtendFmbCell {
protected:
static const int FMB_Info_P = P; //< P >> FMB_Info.P
// FMB_Info_P is declared in FAbstractFmbKernels
static const int MultipoleSize = int(((FMB_Info_P)+1) * ((FMB_Info_P)+2) * 0.5); //< The size of the multipole
FComplexe multipole_exp[MultipoleSize]; //< For multipole extenssion
......
......@@ -53,7 +53,7 @@ class FFmbKernelsForces : public FAbstractFmbKernels<ParticuleClass,CellClass, T
const FComplexe* p_Y_theta_derivated_term = FAbstractFmbKernels<ParticuleClass,CellClass,TreeHeight>::current_thread_Y_theta_derivated+1;
const FComplexe* p_local_exp_term = local->getLocal()+1;
for (int j = 1 ; j <= FAbstractFmbKernels<ParticuleClass,CellClass,TreeHeight>::FMB_Info_P ; ++j ){
for (int j = 1 ; j <= FMB_Info_P ; ++j ){
FComplexe exp_term_aux;
// k=0:
......@@ -197,7 +197,7 @@ class FFmbKernelsForces : public FAbstractFmbKernels<ParticuleClass,CellClass, T
FReal result = 0.0;
FComplexe* p_Y_term = FAbstractFmbKernels<ParticuleClass,CellClass,TreeHeight>::current_thread_Y;
for(int j = 0 ; j<= FAbstractFmbKernels<ParticuleClass,CellClass,TreeHeight>::FMB_Info_P ; ++j){
for(int j = 0 ; j<= FMB_Info_P ; ++j){
// k=0
(*p_Y_term) *= (*local_exp);
result += p_Y_term->getReal();
......
......@@ -60,7 +60,7 @@ public:
FReal result = 0.0;
FComplexe* p_Y_term = FAbstractFmbKernels<ParticuleClass,CellClass,TreeHeight>::current_thread_Y;
for(int j = 0 ; j<= FAbstractFmbKernels<ParticuleClass,CellClass,TreeHeight>::FMB_Info_P ; ++j){
for(int j = 0 ; j<= FMB_Info_P ; ++j){
// k=0
(*p_Y_term) *= (*local_exp);
result += p_Y_term->getReal();
......
......@@ -63,7 +63,7 @@ public:
const FComplexe* p_Y_theta_derivated_term = FAbstractFmbKernels<ParticuleClass,CellClass,TreeHeight>::current_thread_Y_theta_derivated+1;
const FComplexe* p_local_exp_term = local->getLocal()+1;
for (int j = 1 ; j <= FAbstractFmbKernels<ParticuleClass,CellClass,TreeHeight>::FMB_Info_P ; ++j ){
for (int j = 1 ; j <= FMB_Info_P ; ++j ){
FComplexe exp_term_aux;
// k=0:
......@@ -258,7 +258,7 @@ public:
FReal result = 0.0;
FComplexe* p_Y_term = FAbstractFmbKernels<ParticuleClass,CellClass,TreeHeight>::current_thread_Y;
for(int j = 0 ; j<= FAbstractFmbKernels<ParticuleClass,CellClass,TreeHeight>::FMB_Info_P ; ++j){
for(int j = 0 ; j<= FMB_Info_P ; ++j){
// k=0
(*p_Y_term) *= (*local_exp);
result += p_Y_term->getReal();
......
......@@ -4,7 +4,7 @@
// To get memcpy
#include <cstring>
#include "FGlobal.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
......
......@@ -3,6 +3,7 @@
// /!\ Please, you must read the license at the bottom of this page
#include <math.h>
#include "FGlobal.hpp"
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
......
......@@ -45,7 +45,7 @@ public:
/** Custom cell
*
*/
class FmbCell : public FBasicCell, public FExtendFmbCell<2> {
class FmbCell : public FBasicCell, public FExtendFmbCell {
public:
};
......@@ -55,7 +55,7 @@ int main(int , char ** ){
const int NbLevels = 9;//10;
const int SizeSubLevels = 3;//3
FTic counter;
const char* const filename = "testFMAlgorithm.fma"; //"testLoaderFMA.fma" "testFMAlgorithm.fma"
const char* const filename = "testLoaderFMA.fma"; //"testLoaderFMA.fma" "testFMAlgorithm.fma"
FFMALoader<FmbParticule> loader(filename);
if(!loader.isValide()){
......
......@@ -26,7 +26,7 @@
*/
int main(int , char ** ){
const int NbLevels = 10;
const int NbLevels = 9;
const int NbSubLevels = 3;
const long NbPart = 2E6;
FList<FBasicParticule*> particules;
......@@ -59,25 +59,45 @@ int main(int , char ** ){
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
// -----------------------------------------------------
std::cout << "Itering on particules ..." << std::endl;
counter.tic();
FOctree<FBasicParticule, FBasicCell, NbLevels, NbSubLevels>::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
for(int idx = 0 ; idx < NbLevels - 1; ++idx ){
int counter = 0;
do{
++counter;
//counter += octreeIterator.getCurrentList()->getSize();
} while(octreeIterator.moveRight());
octreeIterator.moveUp();
octreeIterator.gotoLeft();
std::cout << "Next level (" << counter << ")...\n";
{
std::cout << "Itering on particules ..." << std::endl;
counter.tic();
FOctree<FBasicParticule, FBasicCell, NbLevels, NbSubLevels>::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
for(int idx = 0 ; idx < NbLevels - 1; ++idx ){
int counter = 0;
do{
++counter;
//counter += octreeIterator.getCurrentList()->getSize();
} while(octreeIterator.moveRight());
octreeIterator.moveUp();
octreeIterator.gotoLeft();
std::cout << "Cells at this level " << counter << " ...\n";
}
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
}
// -----------------------------------------------------
{
std::cout << "Itering on particules fast ..." << std::endl;
counter.tic();
FOctree<FBasicParticule, FBasicCell, NbLevels, NbSubLevels>::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
FOctree<FBasicParticule, FBasicCell, NbLevels, NbSubLevels>::Iterator avoidGoLeft(octreeIterator);
for(int idx = 0 ; idx < NbLevels - 1; ++idx ){
int counter = 0;
do{
} while(octreeIterator.moveRight());
avoidGoLeft.moveUp();
octreeIterator = avoidGoLeft;
}
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
}
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
// -----------------------------------------------------
std::cout << "Deleting particules ..." << std::endl;
counter.tic();
......
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