Commit 7e96583f authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Implemented adaptive routines for AdaptChebSym FMM kernel, Switched base...

Implemented adaptive routines for AdaptChebSym FMM kernel, Switched base kernel from dense to classic (fft accelerated) Unif FMM kernel in AdaptUnif FMM kernel.
parent dec16fcc
......@@ -58,6 +58,11 @@ class FAdaptiveChebSymKernel : FChebSymKernel<CellClass, ContainerClass, MatrixK
, public FAbstractAdaptiveKernel<CellClass, ContainerClass> {
//
typedef FChebSymKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes};
const MatrixKernelClass *const MatrixKernel;
public:
using KernelBaseClass::P2M;
......@@ -74,11 +79,11 @@ public:
// * runtime_error is thrown if the required file is not valid).
// */
FAdaptiveChebSymKernel(const int inTreeHeight, const FReal inBoxWidth,
const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel)
const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel), MatrixKernel(inMatrixKernel)
{}
// /** Copy constructor */
FAdaptiveChebSymKernel(const FAdaptiveChebSymKernel& other)
: KernelBaseClass(other)
: KernelBaseClass(other), MatrixKernel(other.MatrixKernel)
{ }
//
......@@ -88,47 +93,222 @@ public:
//this->~KernelBaseClass() ;
}
void P2M(CellClass* const pole, const int cellLevel, const ContainerClass* const particles) override {
//pole->setDataUp(pole->getDataUp() + particles->getNbParticles());
//
const FPoint LeafCellCenter(KernelBaseClass::getLeafCellCenter(pole->getCoordinate()));
const FReal BoxWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-cellLevel);
//
const FPoint CellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),cellLevel));
const FReal BoxWidth = KernelBaseClass::BoxWidth / FMath::pow(2.0,cellLevel);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
KernelBaseClass::Interpolator->applyP2M(LeafCellCenter, BoxWidth,
pole->getMultipole(idxRhs), particles);
// 1) apply Sy
KernelBaseClass::Interpolator->applyP2M(CellCenter, BoxWidth,
pole->getMultipole(idxRhs), particles);
}
}
void M2M(CellClass* const pole, const int /*poleLevel*/, const CellClass* const subCell, const int /*subCellLevel*/) override {
// pole->setDataUp(pole->getDataUp() + subCell->getDataUp());
void M2M(CellClass* const pole, const int poleLevel, const CellClass* const subCell, const int subCellLevel) override {
const FPoint subCellCenter(KernelBaseClass::getCellCenter(subCell->getCoordinate(),subCellLevel));
const FReal subCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,subCellLevel)));
const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel));
const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel)));
////////////////////////////////////////////////////////////////////////////
/// p^6 version
// allocate memory
FReal* subChildParentInterpolator = new FReal [nnodes * nnodes];
// set child info
FPoint ChildRoots[nnodes], localChildRoots[nnodes];
FChebTensor<ORDER>::setRoots(subCellCenter, subCellWidth, ChildRoots);
// map global position of roots to local position in parent cell
const map_glob_loc map(poleCellCenter, poleCellWidth);
for (unsigned int n=0; n<nnodes; ++n)
map(ChildRoots[n], localChildRoots[n]);
// assemble child - parent - interpolator
KernelBaseClass::Interpolator->assembleInterpolator(nnodes, localChildRoots, subChildParentInterpolator);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy (using tensor product M2M with an interpolator computed on the fly)
// Do NOT reset multipole expansion !
//FBlas::scal(nnodes, FReal(0.), pole->getMultipole(idxRhs));
/// p^6 version
FBlas::gemtva(nnodes, nnodes, FReal(1.),
subChildParentInterpolator,
const_cast<FReal*>(subCell->getMultipole(idxRhs)), pole->getMultipole(idxRhs));
}// NVALS
}
void P2L(CellClass* const local, const int /*localLevel*/, const ContainerClass* const particles) override {
// local->setDataDown(local->getDataDown() + particles->getNbParticles());
void P2L(CellClass* const local, const int localLevel, const ContainerClass* const particles) override {
// Target cell: local
const FReal localCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, localLevel)));
const FPoint localCellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel));
// interpolation points of target (X) cell
FPoint X[nnodes];
FChebTensor<order>::setRoots(localCellCenter, localCellWidth, X);
// read positions
const FReal*const positionsX = particles->getPositions()[0];
const FReal*const positionsY = particles->getPositions()[1];
const FReal*const positionsZ = particles->getPositions()[2];
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// read physicalValue
const FReal*const physicalValues = particles->getPhysicalValues();
// apply P2L
for (unsigned int idxPart=0; idxPart<particles->getNbParticles(); ++idxPart){
const FPoint y = FPoint(positionsX[idxPart],
positionsY[idxPart],
positionsZ[idxPart]);
for (unsigned int m=0; m<nnodes; ++m)
local->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], y) * physicalValues[idxPart];
}
}// NVALS
}
void M2L(CellClass* const local, const int /*localLevel*/, const CellClass* const aNeighbor, const int /*neighborLevel*/) override {
// local->setDataDown(local->getDataDown() + aNeighbor->getDataUp());
void M2L(CellClass* const local, const int localLevel, const CellClass* const pole, const int poleLevel) override {
// Source cell: pole
const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel)));
const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel));
// Target cell: local
const FReal localCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, localLevel)));
const FPoint localCellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel));
// interpolation points of source (Y) and target (X) cell
FPoint X[nnodes], Y[nnodes];
FChebTensor<order>::setRoots(poleCellCenter, poleCellWidth, Y);
FChebTensor<order>::setRoots(localCellCenter, localCellWidth, X);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// Dense M2L
const FReal *const MultipoleExpansion = pole->getMultipole(idxRhs);
for (unsigned int m=0; m<nnodes; ++m)
for (unsigned int n=0; n<nnodes; ++n){
local->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], Y[n]) * MultipoleExpansion[n];
}
}
}
void M2P(const CellClass* const pole, const int /*poleLevel*/, ContainerClass* const particles) override {
// long long int*const particlesAttributes = particles->getDataDown();
// for(int idxPart = 0 ; idxPart < particles->getNbParticles() ; ++idxPart){
// particlesAttributes[idxPart] += pole->getDataUp();
// }
void M2P(const CellClass* const pole, const int poleLevel, ContainerClass* const particles) override {
// Source cell: pole
const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel)));
const FPoint poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel));
// interpolation points of source (Y) cell
FPoint Y[nnodes];
FChebTensor<order>::setRoots(poleCellCenter, poleCellWidth, Y);
// read positions
const FReal*const positionsX = particles->getPositions()[0];
const FReal*const positionsY = particles->getPositions()[1];
const FReal*const positionsZ = particles->getPositions()[2];
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
const FReal *const MultipoleExpansion = pole->getMultipole(idxRhs);
// apply M2P
for (unsigned int idxPart=0; idxPart<particles->getNbParticles(); ++idxPart){
const FPoint x = FPoint(positionsX[idxPart],positionsY[idxPart],positionsZ[idxPart]);
FReal targetValue=0.;
{
for (unsigned int n=0; n<nnodes; ++n)
targetValue +=
MatrixKernel->evaluate(x, Y[n]) * MultipoleExpansion[/*idxLhs*nnodes+*/n];
}
// get potential
FReal*const potentials = particles->getPotentials(/*idxPot*/);
// add contribution to potential
potentials[idxPart] += (targetValue);
}// Particles
}// NVALS
}
void L2L(const CellClass* const local, const int /*localLevel*/, CellClass* const subCell, const int /*subCellLevel*/) override {
// subCell->setDataDown(local->getDataDown() + subCell->getDataDown());
void L2L(const CellClass* const local, const int localLevel, CellClass* const subCell, const int subCellLevel) override {
const FPoint subCellCenter(KernelBaseClass::getCellCenter(subCell->getCoordinate(),subCellLevel));
const FReal subCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,subCellLevel)));
const FPoint localCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel));
const FReal localWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,localLevel)));
////////////////////////////////////////////////////////////////////////////
/// p^6 version
// allocate memory
FReal* subChildParentInterpolator = new FReal [nnodes * nnodes];
// set child info
FPoint ChildRoots[nnodes], localChildRoots[nnodes];
FChebTensor<ORDER>::setRoots(subCellCenter, subCellWidth, ChildRoots);
// map global position of roots to local position in parent cell
const map_glob_loc map(localCenter, localWidth);
for (unsigned int n=0; n<nnodes; ++n)
map(ChildRoots[n], localChildRoots[n]);
// assemble child - parent - interpolator
KernelBaseClass::Interpolator->assembleInterpolator(nnodes, localChildRoots, subChildParentInterpolator);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 2) apply Sx
/// p^6 version
FBlas::gemva(nnodes, nnodes, FReal(1.),
subChildParentInterpolator,
const_cast<FReal*>(local->getLocal(idxRhs)), subCell->getLocal(idxRhs));
}// NVALS
}
void L2P(const CellClass* const local, const int /*cellLevel*/, ContainerClass* const particles) override {
// long long int*const particlesAttributes = particles->getDataDown();
// for(int idxPart = 0 ; idxPart < particles->getNbParticles() ; ++idxPart){
// particlesAttributes[idxPart] += local->getDataDown();
// }
void L2P(const CellClass* const local, const int cellLevel, ContainerClass* const particles) override {
const FPoint CellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),cellLevel));
const FReal BoxWidth = KernelBaseClass::BoxWidth / FMath::pow(2.0,cellLevel);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 2.a) apply Sx
KernelBaseClass::Interpolator->applyL2P(CellCenter, BoxWidth,
local->getLocal(idxRhs), particles);
// 2.b) apply Px (grad Sx)
KernelBaseClass::Interpolator->applyL2PGradient(CellCenter, BoxWidth,
local->getLocal(idxRhs), particles);
}
}
void P2P(ContainerClass* target, const ContainerClass* sources) override {
......
......@@ -22,8 +22,7 @@
#include "Adaptative/FAdaptiveCell.hpp"
#include "Adaptative/FAdaptiveKernelWrapper.hpp"
#include "Adaptative/FAbstractAdaptiveKernel.hpp"
//#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Kernels/Uniform/FUnifDenseKernel.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Kernels/Uniform/FUnifM2LHandler.hpp"
class FTreeCoordinate;
......@@ -55,15 +54,21 @@ class FTreeCoordinate;
*/
template< class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
class FAdaptiveUnifKernel : FUnifDenseKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
class FAdaptiveUnifKernel : FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
, public FAbstractAdaptiveKernel<CellClass, ContainerClass> {
//
typedef FUnifDenseKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
typedef FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
enum {order = ORDER,
nnodes = KernelBaseClass::nnodes};
nnodes = TensorTraits<ORDER>::nnodes};
/// Needed for M2L operator
typedef FUnifM2LHandler<ORDER,MatrixKernelClass::Type> M2LHandlerClass;
const M2LHandlerClass M2LHandler;
// // If we choose to pre-assemble adaptive M2L operators
// // then we need to provide an adaptive M2L handler
// // and the transfer in Fourier space is straightforward.
// typedef FUnifM2LHandler<ORDER,MatrixKernelClass::Type> M2LHandlerClass;
// const M2LHandlerClass M2LHandler;
const MatrixKernelClass *const MatrixKernel;
public:
......@@ -82,11 +87,11 @@ public:
// * runtime_error is thrown if the required file is not valid).
// */
FAdaptiveUnifKernel(const int inTreeHeight, const FReal inBoxWidth,
const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel), M2LHandler(inMatrixKernel, inTreeHeight, inBoxWidth), MatrixKernel(inMatrixKernel)
const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel)/*, M2LHandler(inMatrixKernel, inTreeHeight, inBoxWidth)*/, MatrixKernel(inMatrixKernel)
{}
// /** Copy constructor */
FAdaptiveUnifKernel(const FAdaptiveUnifKernel& other)
: KernelBaseClass(other), M2LHandler(other.M2LHandler), MatrixKernel(other.MatrixKernel)
: KernelBaseClass(other)/*, M2LHandler(other.M2LHandler)*/, MatrixKernel(other.MatrixKernel)
{ }
//
......@@ -435,7 +440,6 @@ public:
// for(int idxPart = 0 ; idxPart < target->getNbParticles() ; ++idxPart){
// particlesAttributes[idxPart] += sources->getNbParticles();
// }
std::cout << "AdapP2P" << std::endl;
}
bool preferP2M(const ContainerClass* const particles) override {
......
......@@ -72,6 +72,27 @@ protected:
BoxCorner.getZ() + (FReal(Coordinate.getZ()) + FReal(.5)) * BoxWidthLeaf);
}
/**
* @brief Return the position of the center of a cell from its tree
* coordinate
* @param FTreeCoordinate
* @param inLevel the current level of Cell
*/
FPoint getCellCenter(const FTreeCoordinate coordinate, int inLevel)
{
//Set the boxes width needed
FReal widthAtCurrentLevel = BoxWidthLeaf*FReal(1 << (TreeHeight-(inLevel+1)));
FReal widthAtCurrentLevelDiv2 = widthAtCurrentLevel/FReal(2.);
//Set the center real coordinates from box corner and widths.
FReal X = BoxCorner.getX() + FReal(coordinate.getX())*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
FReal Y = BoxCorner.getY() + FReal(coordinate.getY())*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
FReal Z = BoxCorner.getZ() + FReal(coordinate.getZ())*widthAtCurrentLevel + widthAtCurrentLevelDiv2;
return FPoint(X,Y,Z);
}
public:
/**
* The constructor initializes all constant attributes and it reads the
......
......@@ -98,7 +98,7 @@ public:
{
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy
FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxRhs));
//FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxRhs));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
......
......@@ -98,7 +98,7 @@ int main(int argc, char ** argv){
//////////////////////////////////////////////////////////////
//
const int NbLevels = FParameters::getValue(argc,argv,"-depth", 7);
const int SizeSubLevels = FParameters::getValue(argc,argv,"subdepth", 3);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-subdepth", 3);
const int sminM = FParameters::getValue(argc,argv,"-sM", P*P*P);
const int sminL = FParameters::getValue(argc,argv,"-sL", P*P*P);
//
......@@ -109,7 +109,7 @@ int main(int argc, char ** argv){
//////////////////////////////////////////////////////////////////////////////////
// Not Random Loader
//////////////////////////////////////////////////////////////////////////////////
const std::string fileName(FParameters::getStr(argc,argv,"-fin", "../Data/prolate50.fma"));
const std::string fileName(FParameters::getStr(argc,argv,"-fin", "../Data/prolate50.out.fma"));
FFmaGenericLoader loader(fileName);
const long int NbPart = loader.getNumberOfParticles() ;
// Random Loader
......@@ -127,10 +127,12 @@ int main(int argc, char ** argv){
std::cout << " criteria SM: "<< sminM <<std::endl
<< " criteria SL: "<< sminL <<std::endl <<std::endl;
//
{
counter.tic();
FReal L= loader.getBoxWidth();
FmaRParticle* particles= new FmaRParticle[NbPart];
//FmaRParticle* particles= new FmaRParticle[NbPart];
FmaRWParticle<8,8>* const particles = new FmaRWParticle<8,8>[NbPart];
FPoint minPos(L,L,L), maxPos(-L,-L,-L);
//
loader.fillParticle(particles,NbPart);
......@@ -146,7 +148,7 @@ int main(int argc, char ** argv){
maxPos.setZ(FMath::Max(maxPos.getZ(),PP.getZ())) ;
//
tree.insert(PP, idxPart, particles[idxPart].getPhysicalValue());
}
counter.tac();
......@@ -167,6 +169,58 @@ int main(int argc, char ** argv){
counter.tac();
std::cout << "Done " << "(@Algorithm = " << counter.elapsed() << " s)." << std::endl;
//
FReal energy= 0.0 , energyD = 0.0 ;
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compute direct energy
/////////////////////////////////////////////////////////////////////////////////////////////////
for(int idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){
energyD += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
}
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compare
/////////////////////////////////////////////////////////////////////////////////////////////////
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
{ // Check that each particle has been summed with all other
// std::cout << "indexPartOrig || DIRECT V fx || FMM V fx" << std::endl;
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
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].getPotential(),potentials[idxPart]);
fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
energy += potentials[idxPart]*physicalValues[idxPart];
// std::cout << indexPartOrig
// << " " << particles[indexPartOrig].getPotential() << " " << particles[indexPartOrig].getForces()[0]
// << " " << potentials[idxPart] << " " << forcesX[idxPart]
// << std::endl;
}
});
}
delete[] particles;
// Print for information
std::cout << "Potential " << potentialDiff << std::endl;
std::cout << "Fx " << fx << std::endl;
std::cout << "Fy " << fy << std::endl;
std::cout << "Fz " << fz << std::endl;
OctreeClass::Iterator octreeIterator(&tree);
std::ofstream file("aa.tree", std::ofstream::out );
......
......@@ -94,7 +94,7 @@ int main(int argc, char ** argv){
//////////////////////////////////////////////////////////////
//
const int NbLevels = FParameters::getValue(argc,argv,"-depth", 3);
const int SizeSubLevels = FParameters::getValue(argc,argv,"subdepth", 2);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-subdepth", 2);
const int sminM = FParameters::getValue(argc,argv,"-sM", P*P*P);
const int sminL = FParameters::getValue(argc,argv,"-sL", P*P*P);
//
......
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