-
Quentin Khan authored
This allows altering the LICENCE text without changing all files. It should avoid merge conflicts in the future.
Quentin Khan authoredThis allows altering the LICENCE text without changing all files. It should avoid merge conflicts in the future.
compareAllPoissonKernels.cpp 39.95 KiB
// See LICENCE file at project root
//
// ==== CMAKE =====
//
// ================
// Keep in private GIT
//
#include <iostream>
#include <stdexcept>
#include <string>
#include <cstdlib>
#include <cstdio>
//
#include "ScalFmmConfig.h"
#include "Utils/FTic.hpp"
#include "Utils/FParameters.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"
#include "Core/FFmmAlgorithmThread.hpp"
#ifdef SCALFMM_USE_BLAS
// chebyshev kernel
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
#endif
//
// spherical kernel
#include "Kernels/Spherical/FSphericalCell.hpp"
#ifdef SCALFMM_USE_BLAS
#include "Kernels/Spherical/FSphericalBlasKernel.hpp"
#include "Kernels/Spherical/FSphericalBlockBlasKernel.hpp"
#endif
//
// taylor kernel
#include "Kernels/Taylor/FTaylorCell.hpp"
#include "Kernels/Taylor/FTaylorKernel.hpp"
//
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
//Classical Spherical kernel
#include "Kernels/Spherical/FSphericalCell.hpp"
#include "Kernels/Spherical/FSphericalKernel.hpp"
//Rotation kernel
#include "Kernels/Rotation/FRotationKernel.hpp"
#include "Kernels/Rotation/FRotationCell.hpp"
#ifdef SCALFMM_USE_FFT
// Uniform grid kernel
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
#endif
#include "Utils/FParameterNames.hpp"
/**
* This program compares two different kernels, eg., the Chebyshev kernel with
* the SphericalBlas kernel.
*/
//
/// \file compareAllPoissonKernels.cpp
//!
//! \brief compareAllPoissonKernels: Driver to compare all different implementations for 1/r kernel.
//!
//! \details This driver compare all different implementations provided by our library for the classical Poisson kernel (1/r)<br>
//! We check the following kernels <br>
//! - Spherical expansion; Spherical with rotation optimization
//! - Taylor expansion; Spherical with rotation optimization
//! - if BLAS is activated
//! - Blas And Block Blas optiization of the M2L operator
//! - Chebychev and symetric Chebychev interpolation
//! - if FFT is activated: interpolation on uniform grid
//!<br>
// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
FHelpDescribeAndExit(argc, argv,
"Driver for testing different approximations for the 1/r kernel.",
FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
FParameterDefinitions::NbThreads, FParameterDefinitions::SHDevelopment);
typedef double FReal;
// get info from commande line
const std::string filename(FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, "../Data/UTest/unitCubeRef20kDouble.bfma"));
const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, omp_get_max_threads());
const int DevP = FParameters::getValue(argc, argv, FParameterDefinitions::SHDevelopment.options, 11);
//
#ifdef _OPENMP
omp_set_num_threads(NbThreads);
#else
std::cout << "\n>> Sequential version.\n" << std::endl;
#endif
std::cout << "Parameters "<< std::endl
<< " Octree Depth \t"<< TreeHeight <<std::endl
<< " SubOctree depth \t"<< SubTreeHeight <<std::endl
<< " Input file name: \t" <<filename <<std::endl
<< " Thread number: \t" << NbThreads <<std::endl
<<std::endl;
// init timer
FTic time;
FFmaGenericLoader<FReal> loader(filename);
// if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");
//
FSize nbParticles = loader.getNumberOfParticles() ;
FmaRWParticle<FReal, 8,8>* const particles = new FmaRWParticle<FReal, 8,8>[nbParticles];
//
loader.fillParticle(particles,nbParticles);
//
////////////////////////////////////////////////////////////////////
// Compute direct energy
FReal energyD =0.0, totPhysicalValue =0.0;
#pragma omp parallel for reduction(+:energyD,totPhysicalValue)
for(FSize idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){
energyD += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
totPhysicalValue += particles[idx].getPhysicalValue() ;
}
std::cout << " Total Physical values: "<< totPhysicalValue <<std::endl;
std::cout << " Energy of the system: "<< energyD <<std::endl;
////////////////////////////////////////////////////////////////////
#ifdef SCALFMM_USE_BLAS
{ // begin Chebyshev kernel
// accuracy
const unsigned int ORDER = 7;
std::cout << "\nFChebKernel FMM ... ORDER: " << ORDER <<std::endl;
// typedefs
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
typedef FChebCell<FReal, ORDER> CellClass;
typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebKernel<FReal, CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
// Create Matrix Kernel
const MatrixKernelClass MatrixKernel;
{ // -----------------------------------------------------
time.tic();
for(FSize idxPart = 0 ; idxPart < nbParticles; ++idxPart){
// put in tree
tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
}
time.tac();
std::cout << "(FChebKernel @ Inserting Particles = "<< time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
{ // -----------------------------------------------------
time.tic();
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
FmmClass algorithm(&tree, &kernels);
algorithm.execute();
time.tac();
std::cout <<"(FChebKernel @Algorithm = " << time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
FReal energy = 0.0;
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize 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] ;
}
});
}
// Print for information
std::cout << "FChebKernel Energy " << FMath::Abs(energy-energyD) /energyD << std::endl;
std::cout << "FChebKernel Potential " << potentialDiff << std::endl;
std::cout << "FChebKernel Fx " << fx << std::endl;
std::cout << "FChebKernel Fy " << fy << std::endl;
std::cout << "FChebKernel Fz " << fz << std::endl;
} // end Chebyshev kernel
{ // begin ChebSymKernel kernel
// accuracy
const unsigned int ORDER = 7;
std::cout << "\nFChebSymKernel FMM ... ORDER: " << ORDER <<std::endl;
// typedefs
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
typedef FChebCell<FReal, ORDER> CellClass;
typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebSymKernel<FReal, CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
// Create Matrix Kernel
const MatrixKernelClass MatrixKernel;
{ // -----------------------------------------------------
time.tic();
for(FSize idxPart = 0 ; idxPart < nbParticles; ++idxPart){
// put in tree
tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
}
time.tac();
std::cout << "(FChebSymKernel @Inserting Particles = "<< time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
{ // -----------------------------------------------------
time.tic();
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
FmmClass algorithm(&tree, &kernels);
algorithm.execute();
time.tac();
std::cout << "(FChebSymKernel @Algorithm = " << time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
FReal energy = 0.0;
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize 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] ;
}
});
}
// Print for information
std::cout << "FChebSymKernel Energy " << FMath::Abs(energy-energyD) /energyD << std::endl;
std::cout << "FChebSymKernel Potential " << potentialDiff << std::endl;
std::cout << "FChebSymKernel Fx " << fx << std::endl;
std::cout << "FChebSymKernel Fy " << fy << std::endl;
std::cout << "FChebSymKernel Fz " << fz << std::endl;
} // end Chebyshev kernel
//
////////////////////////////////////////////////////////////////////
//
{ // begin FFmaBlas kernel FSphericalBlockBlasKernel
std::cout << "\nFFmaBlas FMM ... P: " <<DevP << std::endl;
// typedefs
typedef FSphericalCell<FReal> CellClass;
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
typedef FOctree<FReal, CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FSphericalBlasKernel<FReal, CellClass, ContainerClass > KernelClass;
typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
// init cell class and oct-tree
CellClass::Init(DevP, true); // only for blas
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
{ // -----------------------------------------------------
time.tic();
for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// put in tree
tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
}
time.tac();
std::cout << "( FFmaBlas@Inserting Particles = " << time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
// -----------------------------------------------------
time.tic();
KernelClass kernels(DevP, TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
FmmClass algorithm(&tree, &kernels);
algorithm.execute();
time.tac();
std::cout << "(FFmaBlas @Algorithm = " << time.elapsed() << " s)." << std::endl;
// -----------------------------------------------------
FReal energy = 0.0;
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize 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] ;
}
});
}
// Print for information
std::cout << "FFmaBlas Energy " << FMath::Abs(energy-energyD) /energyD << std::endl;
std::cout << "FFmaBlas Potential " << potentialDiff << std::endl;
std::cout << "FFmaBlas Fx " << fx << std::endl;
std::cout << "FFmaBlas Fy " << fy << std::endl;
std::cout << "FFmaBlas Fz " << fz << std::endl;
} // end FFmaBlas kernel
//
////////////////////////////////////////////////////////////////////
//
{ // begin FFmaBlockBlas kernel
std::cout << "\nFFmaBlockBlas FMM ... P: " <<DevP << std::endl;
// typedefs
typedef FSphericalCell<FReal> CellClass;
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
typedef FOctree<FReal, CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FSphericalBlockBlasKernel<FReal, CellClass, ContainerClass > KernelClass;
typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
// init cell class and oct-tree
CellClass::Init(DevP, true); // only for blas
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
{ // -----------------------------------------------------
time.tic();
for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// put in tree
tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
}
time.tac();
std::cout << "( FFmaBlockBlas@Inserting Particles = " << time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
// -----------------------------------------------------
time.tic();
KernelClass kernels(DevP, TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
FmmClass algorithm(&tree, &kernels);
algorithm.execute();
time.tac();
std::cout << "(FFmaBlockBlas @Algorithm = " << time.elapsed() << " s)." << std::endl;
// -----------------------------------------------------
FReal energy = 0.0;
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize 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] ;
}
});
}
// Print for information
std::cout << "FFmaBlockBlas Energy " << FMath::Abs(energy-energyD) /energyD << std::endl;
std::cout << "FFmaBlockBlas Potential " << potentialDiff << std::endl;
std::cout << "FFmaBlockBlas Fx " << fx << std::endl;
std::cout << "FFmaBlockBlas Fy " << fy << std::endl;
std::cout << "FFmaBlockBlas Fz " << fz << std::endl;
} // end FFmaBlas kernel
#endif
#ifdef SCALFMM_USE_FFT
//
////////////////////////////////////////////////////////////////////
//
{ // begin Lagrange/Uniform Grid kernel
// TODO
// accuracy
const unsigned int ORDER = 7;
std::cout << "\nLagrange FMM ... ORDER " << ORDER <<std::endl;
// typedefs
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
typedef FUnifCell<FReal, ORDER> CellClass;
typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FUnifKernel<FReal, CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
// Create Matrix Kernel
const MatrixKernelClass MatrixKernel;
{ // -----------------------------------------------------
time.tic();
for(FSize idxPart = 0 ; idxPart <nbParticles ; ++idxPart){
// put in tree
tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
}
time.tac();
std::cout << "(Lagrange @Inserting Particles = " << time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
{ // -----------------------------------------------------
time.tic();
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
FmmClass algorithm(&tree, &kernels);
algorithm.execute();
time.tac();
std::cout << "(Lagrange @Algorithm = " << time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
FReal energy = 0.0;
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize 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] ;
}
});
}
// Print for information
std::cout << "Lagrange Energy " << FMath::Abs(energy-energyD) /energyD << std::endl;
std::cout << "Lagrange Potential " << potentialDiff << std::endl;
std::cout << "Lagrange Fx " << fx << std::endl;
std::cout << "Lagrange Fy " << fy << std::endl;
std::cout << "Lagrange Fz " << fz << std::endl;
} // end Lagrange/Uniform Grid kernel
#endif
//
// Spherical approximation
//
{
//const static int P = 10;
typedef FSphericalCell<FReal> CellClass;
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
typedef FOctree<FReal, CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FSphericalKernel<FReal, CellClass, ContainerClass > KernelClass;
typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
#ifndef SCALFMM_USE_BLAS
CellClass::Init(DevP, true);
#endif
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
std::cout << "\nFFmaSpherical FMM ... P: " << DevP<< std::endl;
{ // -----------------------------------------------------
time.tic();
for(FSize idxPart = 0 ; idxPart < nbParticles; ++idxPart){
// put in tree
tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
}
time.tac();
std::cout << "(FFmaSpherical @Inserting Particles = "<< time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
// -----------------------------------------------------
time.tic();
CellClass::Init(DevP);
KernelClass *kernels = new KernelClass(DevP, TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
FmmClass algorithm(&tree, kernels);
algorithm.execute();
time.tac();
std::cout << "(FFmaSpherical @Algorithm = " << time.elapsed() << " s)." << std::endl;
// -----------------------------------------------------
FReal energy = 0.0;
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize 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] ;
}
});
}
// Print for information
std::cout << "FFmaSpherical Energy " << FMath::Abs(energy-energyD) /energyD << std::endl;
std::cout << "FFmaSpherical Potential " << potentialDiff << std::endl;
std::cout << "FFmaSpherical Fx " << fx << std::endl;
std::cout << "FFmaSpherical Fy " << fy << std::endl;
std::cout << "FFmaSpherical Fz " << fz << std::endl;
}
//
// Spherical Rotation
//
{
const static int P = 11;
typedef FRotationCell<FReal, P> CellClass;
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
typedef FOctree<FReal, CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FRotationKernel<FReal, CellClass, ContainerClass , P> KernelClass;
typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
std::cout << "\nFFmaRotation FMM ... P: " << P<< std::endl;
{ // -----------------------------------------------------
time.tic();
for(FSize idxPart = 0 ; idxPart < nbParticles; ++idxPart){
// put in tree
tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
}
time.tac();
std::cout << "(FFmaRotation @Inserting Particles = "<< time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
// -----------------------------------------------------
time.tic();
KernelClass *kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
FmmClass algorithm(&tree, kernels);
algorithm.execute();
time.tac();
std::cout << "(FFmaRotation @Algorithm = " << time.elapsed() << " s)." << std::endl;
// -----------------------------------------------------
FReal energy = 0.0;
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize 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] ;
}
});
}
// Print for information
std::cout << "FFmaRotation Energy " << FMath::Abs(energy-energyD) /energyD << std::endl;
std::cout << "FFmaRotation Potential " << potentialDiff << std::endl;
std::cout << "FFmaRotation Fx " << fx << std::endl;
std::cout << "FFmaRotation Fy " << fy << std::endl;
std::cout << "FFmaRotation Fz " << fz << std::endl;
}
////////////////////////////////////////////////////////////////////
{ // begin Taylor kernel
// accuracy
const unsigned int ORDER = 10;
// typedefs
typedef FTaylorCell<FReal, ORDER,1> CellClass;
std::cout << "\nFFmaTaylor FMM ... ORDER: " << ORDER << std::endl;
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
typedef FOctree<FReal, CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FTaylorKernel<FReal, CellClass,ContainerClass,ORDER,1> KernelClass;
typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
// init cell class and oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
{ // -----------------------------------------------------
time.tic();
for(FSize idxPart = 0 ; idxPart <nbParticles ; ++idxPart){
// put in tree
tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue());
}
time.tac();
std::cout <<"(FFmaTaylor @Inserting Particles = " << time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
// -----------------------------------------------------
time.tic();
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
FmmClass algorithm(&tree, &kernels);
algorithm.execute();
time.tac();
std::cout << "(FFmaTaylor @Algorithm = " << time.elapsed() << " s)." << std::endl;
// -----------------------------------------------------
FReal energy = 0.0;
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize 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] ;
}
});
}
// Print for information
std::cout << "FFmaTaylor Energy " << FMath::Abs(energy-energyD) /energyD << std::endl;
std::cout << "FFmaTaylor Potential " << potentialDiff << std::endl;
std::cout << "FFmaTaylor Fx " << fx << std::endl;
std::cout << "FFmaTaylor Fy " << fy << std::endl;
std::cout << "FFmaTaylor Fz " << fz << std::endl;
} // end FFTaylor kernel
delete[] particles;
return 0;
}