Skip to content

GitLab

Commit cd45b491 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre
Browse files
Options

Add test for adaptive Uniform kernel and update FUnifCell.

parent f12dc2f9
Hide whitespace changes
Inline Side-by-side
Showing with 534 additions and 0 deletions
Src/AdaptiveTree/FAdaptUnifKernel.hpp 0 → 100755
View file @ cd45b491
#ifndef FADAPTUNIFKERNEL_HPP
#define FADAPTUNIFKERNEL_HPP
// ===================================================================================
// Copyright ScalFmm 2011 INRIA,
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#include "Utils/FGlobal.hpp"
#include "Utils/FTrace.hpp"
#include "Utils/FPoint.hpp"
#include "Adaptative/FAdaptiveCell.hpp"
#include "Adaptative/FAdaptiveKernelWrapper.hpp"
#include "Adaptative/FAbstractAdaptiveKernel.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Kernels/Uniform/FUnifM2LHandler.hpp"
class FTreeCoordinate;
// ==== CMAKE =====
// @FUSE_BLAS
// ================
// for verbosity only!!!
//#define COUNT_BLOCKED_INTERACTIONS
// if timings should be logged
//#define LOG_TIMINGS
/**
* @author O. Coulaud
* @class FAdaptUnifKernel
* @brief
* Please read the license
*
* This kernels implement the Lagrange interpolation based FMM operators.
* It implements all interfaces (P2P, P2M, M2M, M2L, L2L, L2P) which are
* required by the FFmmAlgorithm and FFmmAlgorithmThread.
*
* @tparam CellClass Type of cell
* @tparam ContainerClass Type of container to store particles
* @tparam MatrixKernelClass Type of matrix kernel function
* @tparam ORDER Chebyshev interpolation order
*/
template< class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
class FAdaptiveUnifKernel : FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
, public FAbstractAdaptiveKernel<CellClass, ContainerClass> {
//
typedef FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
/// Needed for M2L operator
typedef FUnifM2LHandler<ORDER,MatrixKernelClass::Type> M2LHandlerClass;
const M2LHandlerClass M2LHandler;
public:
using KernelBaseClass::P2M;
using KernelBaseClass::M2M;
using KernelBaseClass::M2L;
using KernelBaseClass::finishedLevelM2L;
using KernelBaseClass::L2L;
using KernelBaseClass::L2P;
using KernelBaseClass::P2P;
using KernelBaseClass::P2PRemote;
// /**
// * The constructor initializes all constant attributes and it reads the
// * precomputed and compressed M2L operators from a binary file (an
// * runtime_error is thrown if the required file is not valid).
// */
FAdaptiveUnifKernel(const int inTreeHeight, const FReal inBoxWidth,
const FPoint& inBoxCenter) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter), M2LHandler(KernelBaseClass::MatrixKernel.getPtr(), inTreeHeight, inBoxWidth)
{}
// /** Copy constructor */
FAdaptiveUnifKernel(const FAdaptiveUnifKernel& other)
: KernelBaseClass(other), M2LHandler(other.M2LHandler)
{ }
//
// /** Destructor */
~FAdaptiveUnifKernel()
{
//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);
//
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy
KernelBaseClass::Interpolator->applyP2M(LeafCellCenter, BoxWidth,
pole->getMultipole(idxRhs), particles);
// 2) apply Discrete Fourier Transform
M2LHandler.applyZeroPaddingAndDFT(pole->getMultipole(idxRhs),
pole->getTransformedMultipole(idxRhs));
}
}
void M2M(CellClass* const pole, const int /*poleLevel*/, const CellClass* const subCell, const int /*subCellLevel*/) override {
// pole->setDataUp(pole->getDataUp() + subCell->getDataUp());
}
void P2L(CellClass* const local, const int /*localLevel*/, const ContainerClass* const particles) override {
// local->setDataDown(local->getDataDown() + particles->getNbParticles());
}
void M2L(CellClass* const local, const int /*localLevel*/, const CellClass* const aNeighbor, const int /*neighborLevel*/) override {
// local->setDataDown(local->getDataDown() + aNeighbor->getDataUp());
}
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 L2L(const CellClass* const local, const int /*localLevel*/, CellClass* const subCell, const int /*subCellLevel*/) override {
// subCell->setDataDown(local->getDataDown() + subCell->getDataDown());
}
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 P2P(ContainerClass* target, const ContainerClass* sources) override {
// long long int*const particlesAttributes = target->getDataDown();
// for(int idxPart = 0 ; idxPart < target->getNbParticles() ; ++idxPart){
// particlesAttributes[idxPart] += sources->getNbParticles();
// }
}
bool preferP2M(const ContainerClass* const particles) override {
return particles->getNbParticles() < 10;
}
bool preferP2M(const int /*atLevel*/, const ContainerClass*const particles[], const int nbContainers) override {
int counterParticles = 0;
for(int idxContainer = 0 ; idxContainer < nbContainers ; ++idxContainer){
counterParticles += particles[idxContainer]->getNbParticles();
}
return counterParticles < 10;
}
};
//
//template < class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
//class FAdaptUnifKernel
// : public FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
//{
// typedef FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
//
//#ifdef LOG_TIMINGS
// FTic time;
// FReal t_m2l_1, t_m2l_2, t_m2l_3;
//#endif
//
//public:
// /**
// * The constructor initializes all constant attributes and it reads the
// * precomputed and compressed M2L operators from a binary file (an
// * runtime_error is thrown if the required file is not valid).
// */
// FAdaptUnifKernel(const int inTreeHeight,
// const FReal inBoxWidth,
// const FPoint& inBoxCenter)
//: KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter)
//{
//
//#ifdef LOG_TIMINGS
// t_m2l_1 = FReal(0.);
// t_m2l_2 = FReal(0.);
// t_m2l_3 = FReal(0.);
//#endif
//}
//
//
// /** Copy constructor */
// FAdaptUnifKernel(const FAdaptUnifKernel& other)
// : KernelBaseClass(other)
// { }
//
//
//
// /** Destructor */
// ~FAdaptUnifKernel()
// {
// this->~KernelBaseClass() ;
//#ifdef LOG_TIMINGS
// std::cout << "- Permutation took " << t_m2l_1 << "s"
// << "\n- GEMMT and GEMM took " << t_m2l_2 << "s"
// << "\n- Unpermutation took " << t_m2l_3 << "s"
// << std::endl;
//#endif
// }
//
//
// void P2MAdapt(CellClass* const ParentCell, const int &level)
// {
// const FPoint LeafCellCenter(KernelBaseClass::getLeafCellCenter(ParentCell->getCoordinate()));
// const FReal BoxWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-level);
// //
// for(int i = 0 ; i <ParentCell->getLeavesSize(); ++i ){
// //
// for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// KernelBaseClass::Interpolator->applyP2M(LeafCellCenter, BoxWidth,
// ParentCell->getMultipole(idxRhs), ParentCell->getLeaf(i)->getSrc());
// }
// }
// }
// void M2MAdapt(CellClass* const FRestrict ParentCell, const int &TreeLevel, const int &numberOfM2M,
// const int * FRestrict ChildLevel , const CellClass*const FRestrict *const FRestrict ChildCells)
// {
// for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// // // apply Sy
// for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
// if (ChildCells[ChildIndex]){
// // KernelBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs), ParentCell->getMultipole(idxRhs));
// }
// }
// }
// }
//
//
//
// void M2L(CellClass* const FRestrict TargetCell,
// const CellClass* SourceCells[343],
// const int /*NumSourceCells*/,
// const int TreeLevel)
// {
//
// }
//
//
// void L2L(const CellClass* const FRestrict ParentCell,
// CellClass* FRestrict *const FRestrict ChildCells,
// const int /*TreeLevel*/)
// {
// // for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// // // apply Sx
// // for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
// // if (ChildCells[ChildIndex]){
// // AbstractBaseClass::Interpolator->applyL2L(ChildIndex, ParentCell->getLocal(idxRhs), ChildCells[ChildIndex]->getLocal(idxRhs));
// // }
// // }
// // }
// }
//
// void L2P(const CellClass* const LeafCell,
// ContainerClass* const TargetParticles)
// {
// KernelBaseClass::L2P(LeafCell,TargetParticles) ;
// }
//
// // void P2P(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions
// // ContainerClass* const FRestrict TargetParticles,
// // const ContainerClass* const FRestrict /*SourceParticles*/,
// // ContainerClass* const NeighborSourceParticles[27],
// // const int /* size */)
// // {
// // DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
// // }
// //
// //
// // void P2PRemote(const FTreeCoordinate& /*inPosition*/,
// // ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
// // ContainerClass* const inNeighbors[27], const int /*inSize*/){
// // DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
// // }
//
//};
//
//
#endif //FADAPTUNIFKERNELS_HPP
// [--END--]
Src/Kernels/Uniform/FUnifCell.hpp
View file @ cd45b491
......@@ -170,6 +170,21 @@ public:
return (NRHS+NLHS)*NVALS*VectorSize * (int) sizeof(FReal) + (NRHS+NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplexe);
}
template <class StreamClass>
friend StreamClass& operator<<(StreamClass& output, const FUnifCell<ORDER, NRHS, NLHS, NVALS>& cell){
output <<" Multipole exp NRHS " << NRHS <<" NVALS " <<NVALS << " VectorSize " << cell.getVectorSize() << std::endl;
for (int rhs= 0 ; rhs < NRHS ; ++rhs) {
const FReal* pole = cell.getMultipole(rhs);
for (int val= 0 ; val < NVALS ; ++val) {
output<< " val : " << val << " exp: " ;
for (int i= 0 ; i < cell.getVectorSize() ; ++i) {
output<< pole[i] << " ";
}
output << std::endl;
}
}
return output;
}
};
......
Tests/noDist/testAdaptiveUnifFMM.cpp 0 → 100644
View file @ cd45b491
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
// ==== CMAKE =====
// @FUSE_BLAS
// ================
#include <iostream>
#include <cstdio>
#include "Utils/FParameters.hpp"
#include "Utils/FTic.hpp"
#include "Containers/FOctree.hpp"
//#include "Containers/FVector.hpp"
//#include "Components/FSimpleLeaf.hpp"
#include "Utils/FPoint.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "Files/FRandomLoader.hpp"
#include "Components/FBasicKernels.hpp"
#include "Components/FSimpleIndexedLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "Adaptative/FAdaptiveCell.hpp"
#include "Adaptative/FAdaptiveKernelWrapper.hpp"
#include "Adaptative/FAbstractAdaptiveKernel.hpp"
//
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Uniform/FUnifCell.hpp"
#include "AdaptiveTree/FAdaptUnifKernel.hpp"
#include "AdaptiveTree/FAdaptTools.hpp"
//
//
#include "Core/FFmmAlgorithm.hpp"
//#include "Core/FFmmAlgorithmThread.hpp"
//#include "Core/FFmmAlgorithmTask.hpp"
/** This program show an example of use of the fmm basic algo
* it also check that each particles is impacted each other particles
*/
void usage() {
std::cout << "Driver to obtain statistics on the octree" << std::endl;
std::cout << "Options "<< std::endl
<< " -help to see the parameters " << std::endl
<< " -depth the depth of the octree "<< std::endl
<< " -subdepth specifies the size of the sub octree " << std::endl
<< " -fin name specifies the name of the particle distribution" << std::endl
<< " -sM s_min^M threshold for Multipole (l+1)^2 for Spherical harmonics"<<std::endl
<< " -sL s_min^L threshold for Local (l+1)^2 for Spherical harmonics"<<std::endl;
}
// Simply create particles and try the kernels
int main(int argc, char ** argv){
//
// accuracy
const unsigned int P = 3 ;
// typedef FTestCell CellClass;
// typedef FAdaptiveTestKernel< CellClass, ContainerClass > KernelClass;
typedef FUnifCell<P> CellClass;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleIndexedLeaf<ContainerClass> LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
//
typedef FAdaptiveUnifKernel<CellClass,ContainerClass,MatrixKernelClass,P> KernelClass;
//
//
typedef FAdaptiveCell< CellClass, ContainerClass > CellWrapperClass;
typedef FAdaptiveKernelWrapper< KernelClass, CellClass, ContainerClass > KernelWrapperClass;
typedef FOctree< CellWrapperClass, ContainerClass , LeafClass > OctreeClass;
// FFmmAlgorithmTask FFmmAlgorithmThread
typedef FFmmAlgorithm<OctreeClass, CellWrapperClass, ContainerClass, KernelWrapperClass, LeafClass > FmmClass;
///////////////////////What we do/////////////////////////////
std::cout << ">> This executable has to be used to test the FMM algorithm.\n";
//////////////////////////////////////////////////////////////
//
const int NbLevels = FParameters::getValue(argc,argv,"-depth", 7);
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);
//
FTic counter;
//////////////////////////////////////////////////////////////////////////////////
// Not Random Loader
//////////////////////////////////////////////////////////////////////////////////
const std::string fileName(FParameters::getStr(argc,argv,"-fin", "../Data/prolate50.fma"));
FFmaGenericLoader loader(fileName);
const long int NbPart = loader.getNumberOfParticles() ;
// Random Loader
//const int NbPart = FParameters::getValue(argc,argv,"-nb", 2000000);
// FRandomLoader loader(NbPart, 1, FPoint(0.5,0.5,0.5), 1);
//////////////////////////////////////////////////////////////////////////////////
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Creating & Inserting " << NbPart << " particles ..." << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
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];
FPoint minPos(L,L,L), maxPos(-L,-L,-L);
//
loader.fillParticle(particles,NbPart);
for(int idxPart = 0 ; idxPart < NbPart; ++idxPart){
const FPoint PP(particles[idxPart].getPosition() ) ;
//
minPos.setX(FMath::Min(minPos.getX(),PP.getX())) ;
minPos.setY(FMath::Min(minPos.getY(),PP.getY())) ;
minPos.setZ(FMath::Min(minPos.getZ(),PP.getZ())) ;
maxPos.setX(FMath::Max(maxPos.getX(),PP.getX())) ;
maxPos.setY(FMath::Max(maxPos.getY(),PP.getY())) ;
maxPos.setZ(FMath::Max(maxPos.getZ(),PP.getZ())) ;
//
tree.insert(PP, idxPart, particles[idxPart].getPhysicalValue());
}
counter.tac();
std::cout << "Data are inside the box delimited by "<<std::endl
<< " Min corner: "<< minPos<<std::endl
<< " Max corner: "<< maxPos<<std::endl <<std::endl;
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << " s)." << std::endl;
}
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Working on particles ..." << std::endl;
counter.tic();
KernelWrapperClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());; // FTestKernels FBasicKernels
FmmClass algo(&tree,&kernels); //FFmmAlgorithm FFmmAlgorithmThread
algo.execute();
counter.tac();
std::cout << "Done " << "(@Algorithm = " << counter.elapsed() << " s)." << std::endl;
OctreeClass::Iterator octreeIterator(&tree);
std::ofstream file("aa.tree", std::ofstream::out );
//
////////////////////////////////////////////////////////////////////
// Export adaptive tree in our format
////////////////////////////////////////////////////////////////////
//
// -----------------------------------------------------
//
//
// Set Global id
//
long int idCell = setGlobalID(tree);
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
tree.forEachCellLeaf([&](CellWrapperClass* cell, LeafClass* leaf){
file << "Cell Id " << cell->getGlobalId( ) << " Nb particles "<< leaf->getSrc()->getNbParticles()<<std::endl;
});
octreeIterator.gotoTop() ; // here we are at level 1 (first child)
// octreeIterator.moveDown() ;
octreeIterator.gotoLeft();
// octreeIterator.moveDown() ; // We are at the levell 2
std::cout << " Number of Cells: " << idCell <<std::endl;
//
std::cout << "Top of the octree " << octreeIterator.level() << std::endl ;
for(int idxLevel = 1; idxLevel < NbLevels ; ++idxLevel){
file << std::endl << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<< std::endl;
file << " Level " << idxLevel <<" Level "<< octreeIterator.level()<< " -- leave level " << std::boolalpha << octreeIterator.isAtLeafLevel() << std::endl;
do{
if(octreeIterator.getCurrentCell()->hasDevelopment()){
file <<"Cell id "<< octreeIterator.getCurrentCell()->getGlobalId( ) << " "<<*(octreeIterator.getCurrentCell())<< std::endl ;
}
} while(octreeIterator.moveRight());
octreeIterator.moveDown() ;
octreeIterator.gotoLeft();
}
std::cout << " END " << std::endl;
// Check
octreeIterator.gotoBottomLeft();
do {
std::cout << " Level " <<octreeIterator.level() <<std::endl;
}while(octreeIterator.moveUp() );
std::cout << " RETURN 0 " << 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