Skip to content

S
ScalFMM
Commit 4c309ff5 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files
Options

Modifs Taylor

parent f5c31024
Hide whitespace changes
Inline Side-by-side
Showing with 68 additions and 89 deletions
Src/Kernels/Taylor/FTaylorKernel_Opt.hpp
View file @ 4c309ff5
......@@ -510,48 +510,13 @@ public:
FReal xc = cellCenter.getX(), yc = cellCenter.getY(), zc = cellCenter.getZ() ;
FReal dx[3] ;
for(int idPart=0 ; idPart<nbPart ; ++idPart){
// // printf("P2M :: part : X: %f, Y: %f, Z:%f Q %f\n", posX[idPart] , posY[idPart] , posZ[idPart] ,phyValue[idPart]);
// //#ifndef OC
// // compute the distance to the centre
// FReal dx = xc - posX[idPart] ;
// FReal dy = yc - posY[idPart] ;
// FReal dz = zc - posZ[idPart] ;
// // Precompute the arrays of dx^i
// arrayDX[0] = 1.0 ;
// arrayDY[0] = 1.0 ;
// arrayDZ[0] = 1.0 ;
// for (int i = 1 ; i <= P ; ++i) {
// arrayDX[i] = dx * arrayDX[i-1] ;
// arrayDY[i] = dy * arrayDY[i-1] ;
// arrayDZ[i] = dz * arrayDZ[i-1] ;
// }
// //
// //Iterating over MutlipoleVector
// //
// // a = 0, b = 0, c = 0;
// for(int i=0,a = 0, b = 0, c = 0 ; i<SizeVector ; ++i)
// {
// //
// // update needed values
// //
// facto = fact3int(a,b,c);
// coeff = factorials[a+b+c]/(facto*facto);
// //
// //Computation
// //
// multipole[i] += coeff*phyValue[idPart]*arrayDX[a]*arrayDY[b]*arrayDZ[c];
// // multipole[i] += phyValue[idPart]*arrayDX[a]*arrayDY[b]*arrayDZ[c];
// // std::cout << " (a,b,c), ("<<a<<" , " << b << " , "<< c << " ) - j " << i <<std::endl;
// incPowers(&a,&b,&c); //inc powers of expansion
// } // end loop on multipoles
// //
// printf("P2M :: part : X: %f, Y: %f, Z:%f Q %f\n", posX[idPart] , posY[idPart] , posZ[idPart] ,phyValue[idPart]);
//
dx[0] = xc - posX[idPart] ;
dx[1] = yc - posY[idPart] ;
dx[2] = zc - posZ[idPart] ;
multipole2[0] = phyValue[idPart] ;
// printf(" : DX: %f, DY: %f, DZ:%f \n",dx[0] ,dx[01] ,dx[2] );
multipole2[0] = phyValue[idPart] ;
//
int leading[3] = {0,0,0 } ;
for (int k=1, t=1, tail=1; k <= P; ++k, tail=t)
......@@ -571,12 +536,15 @@ public:
for( i=0 ; i < SizeVector ; ++i)
{
multipole[i] += multipole2[i] ;
multipole2[i] = 0.0;
}
} // end loop on particles
// Multiply by the coefficient
for( i=0 ; i < SizeVector ; ++i)
{
multipole[i] *=_coeffPoly[i] ;
// printf(" M[%d] %f\n", i,multipole[i] );
}
// //
// // Nouveau parcours graded lexicographic order.
......@@ -760,12 +728,12 @@ public:
this->computeFullDerivative(dx,dy,dz,this->_PsiVector);
const FReal * multipole = distantNeighbors[idxNeigh]->getMultipole();
FReal res = FReal(0);
for(int k=0; k<SizeVector ; ++k)
{
res += _PsiVector[k]*multipole[k];
}
printf("res : %.10f size : %d\n",res,sizeof(res));
// FReal res = FReal(0);
// for(int k=0; k<SizeVector ; ++k)
// {
// res += _PsiVector[k]*multipole[k];
// }
// printf("res : %.10f size : %d\n",res,sizeof(res));
//Iteration over Multipole / Local
int al=0,bl=0,cl=0; // For local array
int am,bm,cm; // For distant array
......@@ -787,7 +755,7 @@ public:
for(int j = 0 ; j < SizeVector ; ++j){ //corresponding powers am,bm,cm
int idxPsi = powerToIdx(al+am,bl+bm,cl+cm);
tmp += this->_PsiVector[idxPsi]*multipole[j];
if((al+bl+cl)==0) {printf("coeffL:%f (%d,%d,%d)--> mult x psi ::> %f * %f tmp: %f\n",coeffL,am,bm,cm,multipole[j],_PsiVector[idxPsi],tmp);}
// if((al+bl+cl)==0) {printf("coeffL:%f (%d,%d,%d)--> mult x psi ::> %f * %f tmp: %f\n",coeffL,am,bm,cm,multipole[j],_PsiVector[idxPsi],tmp);}
//updating a,b,c
incPowers(&am,&bm,&cm);
}
......@@ -808,10 +776,10 @@ public:
// incPowers(&x,&y,&z);
// }
int x=0,y=0,z=0;
printf("Cell in (%f,%f,%f), leading to (dx,dy,dz) = (%f,%f,%f)\n",locCenter.getX(),locCenter.getY(),locCenter.getZ(),dx,dy,dz);
// printf("Cell in (%f,%f,%f), leading to (dx,dy,dz) = (%f,%f,%f)\n",locCenter.getX(),locCenter.getY(),locCenter.getZ(),dx,dy,dz);
for(int m=0 ; m < SizeVector ; ++m)
{
printf("Cell in (%f,%f,%f) :: (%d,%d,%d) coeffL:%f ::> %.8f \n",dx,dy,dz,x,y,z,factorials[x+y+z]/(fact3int(x,y,z)*fact3int(x,y,z)),iterLocal[m]);
// printf("Cell in (%f,%f,%f) :: (%d,%d,%d) coeffL:%f ::> %.8f \n",dx,dy,dz,x,y,z,factorials[x+y+z]/(fact3int(x,y,z)*fact3int(x,y,z)),iterLocal[m]);
incPowers(&x,&y,&z);
}
}
......@@ -851,7 +819,7 @@ public:
// if child exists
if(childCell[idxChild]){
FReal* FRestrict childExpansion = childCell[idxChild]->getLocal() ;
const FPoint& childCenter =getCellCenter(childCell[idxChild]->getCoordinate(),inLevel+1);
const FPoint& childCenter = getCellCenter(childCell[idxChild]->getCoordinate(),inLevel+1);
//Set the distance between centers of cells
// Child - father
......@@ -911,7 +879,7 @@ public:
{
// printf("L2P %%%%%%%%%%%%%%%%%%%%%%%%%%\n");
FPoint locCenter = getLeafCenter(local->getCoordinate());
printf("L2P :: locCenter (%f,%f,%f)\n",locCenter.getX(),locCenter.getY(),locCenter.getZ());
// printf("L2P :: locCenter (%f,%f,%f)\n",locCenter.getX(),locCenter.getY(),locCenter.getZ());
//Iterator over particles
int nbPart = particles->getNbParticles();
......@@ -942,7 +910,7 @@ public:
FReal dx = posX[i] - locCenter.getX();
FReal dy = posY[i] - locCenter.getY();
FReal dz = posZ[i] - locCenter.getZ();
printf("Distance : (%f,%f,%f)\n",dx,dy,dz);
// printf("Distance : (%f,%f,%f)\n",dx,dy,dz);
// Precompute an arrays of Array[i] = dx^(i-1)
arrayDX[0] = 0.0 ;
arrayDY[0] = 0.0 ;
......@@ -972,12 +940,13 @@ public:
//
incPowers(&a,&b,&c);
}
targetsPotentials[i] = partPhyValue*locPot ;
targetsPotentials[i] =/* partPhyValue*/locPot ;
forceX[i] = partPhyValue*locForceX ;
forceY[i] = partPhyValue*locForceY ;
forceZ[i] = partPhyValue*locForceZ ;
// printf(" L2P part %d %f %f %f %f\n",i,targetsPotentials[i],forceX[i],forceY[i], forceZ[i]);
}
// printf(" L2P end\n");
}
/**
......
Tests/Kernels/testNewCompareKernels.cpp
View file @ 4c309ff5
......@@ -55,7 +55,7 @@
//
// taylor kernel
#include "../../Src/Kernels/Taylor/FTaylorCell.hpp"
#include "../../Src/Kernels/Taylor/FTaylorKernel.hpp"
#include "../../Src/Kernels/Taylor/FTaylorKernel_Opt.hpp"
//
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"
......@@ -74,7 +74,7 @@ int main(int argc, char* argv[])
const unsigned int TreeHeight = FParameters::getValue(argc, argv, "-h", 5);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, "-sh", 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, "-t", omp_get_max_threads());
const int DevP = FParameters::getValue(argc, argv, "-p", 6);
const int DevP = FParameters::getValue(argc, argv, "-p", 7);
#ifdef _OPENMP
omp_set_num_threads(NbThreads);
std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
......@@ -197,14 +197,14 @@ int main(int argc, char* argv[])
// Print for information
std::cout << "Potential " << potentialDiff << std::endl;
std::cout << "Fx " << potentialDiff << std::endl;
std::cout << "Fy " << potentialDiff << std::endl;
std::cout << "Fz " << potentialDiff << std::endl;
std::cout << "Fx " << fx << std::endl;
std::cout << "Fy " << fy << std::endl;
std::cout << "Fz " << fz << std::endl;
} // end Chebyshev kernel
//
////////////////////////////////////////////////////////////////////
/:
//
{ // begin FFmaBlas kernel
// accuracy
......@@ -276,15 +276,15 @@ int main(int argc, char* argv[])
// Print for information
std::cout << "Potential " << potentialDiff << std::endl;
std::cout << "Fx " << potentialDiff << std::endl;
std::cout << "Fy " << potentialDiff << std::endl;
std::cout << "Fz " << potentialDiff << std::endl;
std::cout << "Fx " << fx << std::endl;
std::cout << "Fy " << fy << std::endl;
std::cout << "Fz " << fz << std::endl;
} // end FFmaBlas kernel
#endif
////////////////////////////////////////////////////////////////////
{ // begin SphericalRotation kernel
// accuracy
// accuracy
// const int DevP = FParameters::getValue(argc, argv, "-p", 6);
// typedefs
......@@ -352,9 +352,9 @@ int main(int argc, char* argv[])
// Print for information
std::cout << "Potential " << potentialDiff << std::endl;
std::cout << "Fx " << potentialDiff << std::endl;
std::cout << "Fy " << potentialDiff << std::endl;
std::cout << "Fz " << potentialDiff << std::endl;
std::cout << "Fx " << fx << std::endl;
std::cout << "Fy " << fy << std::endl;
std::cout << "Fz " << fz << std::endl;
} // end FFSphericalRotation kernel
//
////////////////////////////////////////////////////////////////////
......@@ -414,22 +414,21 @@ int main(int argc, char* argv[])
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].potential,potentials[idxPart]);
fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
}
}
});
}
// Print for information
std::cout << "Potential " << potentialDiff << std::endl;
std::cout << "Fx " << potentialDiff << std::endl;
std::cout << "Fy " << potentialDiff << std::endl;
std::cout << "Fz " << potentialDiff << std::endl;
std::cout << "Fx " << fx << std::endl;
std::cout << "Fy " << fy << std::endl;
std::cout << "Fz " << fz << std::endl;
} // end FFTaylor kernel
delete[] particles;
......
Tests/Kernels/testTaylorSimple.cpp
View file @ 4c309ff5
......@@ -20,12 +20,14 @@
#include "../../Src/Core/FFmmAlgorithmTask.hpp"
#include "../../Src/Kernels/Taylor/FTaylorCell.hpp"
#include "../../Src/Kernels/Taylor/FTaylorKernel.hpp"
#include "../../Src/Kernels/Taylor/FTaylorKernel_Opt.hpp"
#include "../../Src/Files/FFmaLoader.hpp"
int main(int argc,char* argv[]){
static const int P = 8;
const int P = 10;
const int NbLevels = FParameters::getValue(argc, argv, "-h", 3);
const int SizeSubLevels = FParameters::getValue(argc, argv, "-sh", 2);
static const int order = 1;
FPoint rootCenter(FReal(0.0),FReal(0.0),FReal(0.0));
FReal boxWidth = FReal(8);
......@@ -42,19 +44,20 @@ int main(int argc,char* argv[]){
typedef FFmmAlgorithmTask<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClassTask;
const int NbLevels = 4;
const int SizeSubLevels = 1;
FTic counter;
OctreeClass tree(NbLevels, SizeSubLevels, boxWidth, rootCenter);
FPoint part1Pos = FPoint(FReal(3.25),FReal(0.75),FReal(0.75));
FReal physVal1 = 2;
FReal physVal1 = -1;
FPoint part3Pos = FPoint(FReal(2.25),FReal(1.75),FReal(0.75));
FReal physVal3 = -1;
FPoint part2Pos = FPoint(FReal(-3.25),FReal(0.75),FReal(0.75));
FPoint part2Pos = FPoint(FReal(-3.25),FReal(1.75),FReal(0.5));
FReal physVal2 = 2;
tree.insert(part1Pos,physVal1);
tree.insert(part3Pos,physVal3);
tree.insert(part2Pos,physVal2);
......@@ -67,11 +70,12 @@ int main(int argc,char* argv[]){
{ // get sum forces&potential
FReal potential = 0;
FReal Energy = 0.0;
tree.forEachLeaf([&](LeafClass* leaf){
FReal fx = 0.0, fy = 0.0, fz = 0.0;
const FReal * FRestrict charges = leaf->getTargets()->getPhysicalValues();
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
......@@ -79,29 +83,36 @@ int main(int argc,char* argv[]){
const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
potential += potentials[idxPart];
Energy += charges[idxPart]*potentials[idxPart];
fx += forcesX[idxPart];
fy += forcesY[idxPart];
fz += forcesZ[idxPart];
printf("part : %f, fx = %f, fy = %f, fz = %f\n",leaf->getTargets()->getPositions()[0][0],fx,fy,fz);
printf(" part : %f, p= %f fx = %f, fy = %f, fz = %f\n",leaf->getTargets()->getPositions()[0][0],Energy,fx,fy,fz);
}
});
//
FReal potTheo ;
//
FReal dx = part1Pos.getX() - part2Pos.getX();
FReal dy = part1Pos.getY() - part2Pos.getY();
FReal dz = part1Pos.getZ() - part2Pos.getZ();
FReal dx1 = part3Pos.getX() - part2Pos.getX();
FReal dy1 = part3Pos.getY() - part2Pos.getY();
FReal dz1 = part3Pos.getZ() - part2Pos.getZ();
FReal dx2 = part1Pos.getX() - part3Pos.getX();
FReal dy2 = part1Pos.getY() - part3Pos.getY();
FReal dz2 = part1Pos.getZ() - part3Pos.getZ();
std::cout << dx <<" " <<dy <<" "<< dz <<" " <<std::endl;
std::cout << dx1 <<" " <<dy1 <<" "<< dz1 <<" " <<std::endl;
std::cout << dx2 <<" " <<dy2 <<" "<< dz2 <<" " <<std::endl;
//
FReal potTheo = physVal2*physVal1*FMath::Sqrt(FReal(1) / (dx*dx + dy*dy + dz*dz));
potTheo = physVal2*physVal1 / FMath::Sqrt(dx*dx + dy*dy + dz*dz)+ physVal2*physVal3 / FMath::Sqrt(dx1*dx1 + dy1*dy1 + dz1*dz1) + physVal1*physVal3/ FMath::Sqrt(dx2*dx2 + dy2*dy2 + dz2*dz2);
FReal coeffa = physVal2*physVal1*(FMath::Sqrt(FReal(1.0) / (dx*dx + dy*dy + dz*dz)) ) / (dx*dx + dy*dy + dz*dz);
potential *=FReal(0.5) ;
printf("Exact potential : %f Computed potential : %f Error: %e \n",potTheo, potential,std::abs(potTheo- potential));
printf("Exact Force : %f %f : %f \n",dx*coeffa,dy*coeffa,dz*coeffa);
Energy *=FReal(0.5) ;
printf("Exact potential : %f Computed potential : %f Error: %e \n",potTheo, Energy,std::abs(potTheo- Energy));
printf("Exact Force : %f %f : %f \n",dx*coeffa,dy*coeffa,dz*coeffa);
std::cout << Energy/potTheo << std::endl;
}
return 0;
}
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