Skip to content

GitLab

Attention une mise à jour du serveur va être effectuée le lundi 17 mai entre 13h et 13h30. Cette mise à jour va générer une interruption du service de quelques minutes.

Commit a3ca736e authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files
Options

mv these files in SalFMM-Tools project

parent f5a3da89
Hide whitespace changes
Inline Side-by-side
Showing with 0 additions and 454 deletions
Utils/src/ChebyshevInterpolationCmpAlgo.cpp deleted 100644 → 0
View file @ f5a3da89
// ==== CMAKE =====
// @FUSE_BLAS
// ================
// Keep in private GIT
//
//
/** \brief Chebyshev FMM example
*
* \file
* \authors B. Bramas, O. Coulaud
*
* This program runs the FMM Algorithm with the interpolation kernel based on
* Chebyshev (grid points) interpolation (1/r kernel). It then compares the
* results with a direct computation.
*/
#include <string>
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
//
template<typename FReal, int ORDER>
using FInterpolationCell = FChebCell<FReal, ORDER>;
template<typename FReal, typename GroupCellClass,
typename GroupContainerClass,
typename MatrixKernelClass, int ORDER>
using FInterpolationKernel = FChebSymKernel<FReal,
GroupCellClass,
GroupContainerClass,
MatrixKernelClass,
ORDER> ;
const std::string interpolationType("Chebyshev interpolation");
#include "sharedMemoryInterpolatiomCmpAlgos.hpp"
Utils/src/UniformInterpolationCmpAlgo.cpp deleted 100644 → 0
View file @ f5a3da89
// ==== CMAKE =====
// @FUSE_BLAS
// @FUSE_FFT
// ================
// Keep in private GIT
//
//
/** \brief Uniform FMM example
*
* \file
* \authors B. Bramas, O. Coulaud
*
* This program runs the FMM Algorithm with the interpolation kernel based on
* uniform (grid points) interpolation (1/r kernel). It then compares the
* results with a direct computation.
*/
#include <string>
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
//
template<typename FReal, int ORDER>
using FInterpolationCell = FUnifCell<FReal, ORDER>;
template<typename FReal, typename GroupCellClass,
typename GroupContainerClass,
typename MatrixKernelClass, int ORDER>
using FInterpolationKernel = FUnifKernel<FReal,
GroupCellClass,
GroupContainerClass,
MatrixKernelClass,
ORDER> ;
const std::string interpolationType("Uniform interpolation");
#include "sharedMemoryInterpolatiomCmpAlgos.hpp"
Utils/src/sharedMemoryInterpolatiomCmpAlgos.hpp deleted 100644 → 0
View file @ f5a3da89
// -*-c++-*-
// See LICENCE file at project root
// ==== CMAKE =====
// @FUSE_FFT
// @FUSE_BLAS
// ==== Git =====
// ================
/** \brief Uniform FMM example
*
* \file
* \authors B. Bramas, O. Coulaud
*
* This program runs the FMM Algorithm with the interpolation kernel based on
* either uniform (grid points) or Chebychev interpolation (1/r kernel).
* It then compares the results with a direct computation.
*/
#include <iostream>
#include <iomanip>
#include <memory>
#include <cstdio>
#include <cstdlib>
#include <string>
#include "ScalFmmConfig.h"
#include "Utils/FGlobal.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
#include "Files/FFmaGenericLoader.hpp"
// Leaves
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
// Octree
#include "Containers/FOctree.hpp"
//
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
//
#ifdef _OPENMP
#include "Core/FFmmAlgorithmTask.hpp"
#include "Core/FFmmAlgorithmNewTask.hpp"
#include "Core/FFmmAlgorithmSectionTask.hpp"
#include "Core/FFmmAlgorithmOmp4.hpp"
#else
#include "Core/FFmmAlgorithm.hpp"
#endif
//
// Order of the Interpolation approximation
static constexpr unsigned ORDER = 6 ;
using FReal = double;
// 1/r kernel
//
using MatrixKernelType = FInterpMatrixKernelR<FReal> ;
//
/// definition of the common tree structure
using CellType = FInterpolationCell<FReal, ORDER>;
//using CellUpType = typename CellType::multipole_t;
//using CellDownType = typename CellType::local_expansion_t;
//using CellSymbType = FSymbolicData;
using ContainerType = FP2PParticleContainerIndexed<FReal>;
using LeafType = FSimpleLeaf<FReal, ContainerType > ;
using OctreeType = FOctree<FReal, CellType,ContainerType,LeafType>;
using KernelType = FInterpolationKernel<FReal,CellType,ContainerType,MatrixKernelType,ORDER> ;
#ifdef _OPENMP
using TaskFmmAlgo = FFmmAlgorithmTask<OctreeType,CellType,ContainerType,KernelType,LeafType>;
using TaskNewFmmAlgo = FFmmAlgorithmNewTask<OctreeType,CellType,ContainerType,KernelType,LeafType>;
using SectionTaskFmmAlgo = FFmmAlgorithmSectionTask<OctreeType,CellType,ContainerType,KernelType,LeafType> ;
using Omp4FmmAlgo = FFmmAlgorithmOmp4<OctreeType,CellType,ContainerType,KernelType,LeafType>;
#else
using FmmType = FFmmAlgorithm<OctreeType,CellType,ContainerType,KernelType,LeafType>;
#endif
// Simply create particles and try the kernels
int main(int argc, char* argv[]) {
const FParameterNames LocalOptionAlgo= { {"-algo"} , " Algorithm to run (task, newtask, sectiontask, omp4)\n"};
const FParameterNames LocalOptionCmp = {
{"-cmp"} , "Use to check the result with the exact solution given in the input file\n" };
FHelpDescribeAndExit(
argc, argv,
"Driver for Lagrange Or Chebychev interpolation kernel (1/r kernel).",
FParameterDefinitions::OctreeHeight,
FParameterDefinitions::OctreeSubHeight,
FParameterDefinitions::InputFile,
FParameterDefinitions::OutputFile,
FParameterDefinitions::NbThreads,
LocalOptionAlgo,LocalOptionCmp
);
const std::string defaultFile(SCALFMMDataPath+"unitCubeXYZQ100.bfma" );
const std::string filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFile.c_str());
const int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
const int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
#ifdef _OPENMP
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, omp_get_max_threads());
omp_set_num_threads(NbThreads);
std::cout << "\n>> Using " << NbThreads << " threads.\n" << std::endl;
#else
const int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
std::cout << "\n>> Sequential version.\n" << std::endl;
#endif
//
{
std::string indent(" ");
auto w = std::setw(18);
std::cout << "Parameters" << std::endl << std::left
<< indent << w << "Octree Depth: " << TreeHeight << std::endl
<< indent << w << "SubOctree depth: " << SubTreeHeight << std::endl
<< indent << w << "Input file name: " << filename << std::endl
<< indent << w << "Thread number: " << NbThreads << std::endl
<< std::endl;
}
//
// init timer
FTic time;
// open particle file
////////////////////////////////////////////////////////////////////
//
FFmaGenericLoader<FReal> loader(filename);
if (loader.getNbRecordPerline() !=8 && FParameters::existParameter(argc, argv, "-cmp") ){
std::cerr << "File should contain 8 data to read (x,y,z,q,p,fx,fy,fz)\n";
std::exit(EXIT_FAILURE);
}
FSize nbParticles = loader.getNumberOfParticles() ;
FmaRWParticle<FReal,8,8> * const particles = new FmaRWParticle<FReal, 8,8>[nbParticles];
//
//
// open particle file
////////////////////////////////////////////////////////////////////
//
if(loader.getNbRecordPerline() == 8 ){
loader.fillParticle(particles,nbParticles);
}
else {
FmaRWParticle<FReal,4,4> * const part = new FmaRWParticle<FReal, 4,4>[nbParticles];
loader.fillParticle(part,nbParticles);
for(int idx = 0 ; idx < nbParticles ; ++idx){
auto p = part[idx].getPosition() ;
particles[idx].setPosition(p) ;
particles[idx].setPhysicalValue(part[idx].getPhysicalValue() );
}
}
time.tac();
//
// init oct-tree
OctreeType tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
FReal energyD = 0.0 ;
//
// Read particles and insert them in octree
{ // -----------------------------------------------------
std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
<< " particles ..." << std::endl;
std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
time.tic();
//
//
//
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compute direct energy
/////////////////////////////////////////////////////////////////////////////////////////////////
for(int idx = 0 ; idx < nbParticles ; ++idx){
tree.insert(particles[idx].getPosition() , idx, particles[idx].getPhysicalValue() );
energyD += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
}
time.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = "
<< time.elapsed() << " s) ." << std::endl;
}
////////////////////////////////////////////////////////////////////
//
// Execute FMM Algorithm
//
////////////////////////////////////////////////////////////////////
{ // -----------------------------------------------------
std::cout << "\n" << interpolationType << " FMM (ORDER= "<< ORDER << ") ... " << std::endl;
const MatrixKernelType MatrixKernel;
time.tic();
//
std::unique_ptr<KernelType> kernels(new KernelType(TreeHeight, loader.getBoxWidth(),
loader.getCenterOfBox(),&MatrixKernel));
std::string algoStr = FParameters::getStr(argc,argv,"-algo", "task");
//
TaskFmmAlgo algo1(&tree, kernels.get() );
TaskNewFmmAlgo algo2(&tree, kernels.get() );
SectionTaskFmmAlgo algo3(&tree, kernels.get() );
Omp4FmmAlgo algo4(&tree, kernels.get() );
FAbstractAlgorithm * algo = nullptr;
FAlgorithmTimers * timer = nullptr;
if( "task" == algoStr) {
algo = &algo1 ;
timer = &algo1 ;
} else if( "newtask" == algoStr) {
algo = &algo2 ;
timer = &algo2 ;
} else if ( "sectiontask" == algoStr ) {
algo = &algo3 ;
timer = &algo3 ;
} else if ( "omp4" == algoStr ) {
algo = &algo4 ;
timer = &algo4 ;
} else {
std::cout << "Unknown algorithm: " << algoStr << std::endl;
std::exit(EXIT_FAILURE);
}
std::cout << "Algorithm to check: "<< algoStr <<std::endl;
time.tic();
//
algo->execute(); // Here the call of the FMM algorithm
//
time.tac();
std::cout << "Timers Far Field \n"
<< "P2M " << timer->getTime(FAlgorithmTimers::P2MTimer) << " seconds\n"
<< "M2M " << timer->getTime(FAlgorithmTimers::M2MTimer) << " seconds\n"
<< "M2L " << timer->getTime(FAlgorithmTimers::M2LTimer) << " seconds\n"
<< "L2L " << timer->getTime(FAlgorithmTimers::L2LTimer) << " seconds\n"
<< "L2P " << timer->getTime(FAlgorithmTimers::L2PTimer) << " seconds\n"
<< "P2P and L2P " << timer->getTime(FAlgorithmTimers::NearTimer) << " seconds\n"
<< std::endl;
std::cout << "Done " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
}
// -----------------------------------------------------
//
// Some output
//
//
{ // -----------------------------------------------------
FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= loader.getNumberOfParticles() -1; ;
FReal energy =0.0 ;
//
// Loop over all leaves
//
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
std::cout << std::scientific;
std::cout.precision(10) ;
tree.forEachLeaf([&](LeafType* leaf){
const FReal*const posX = leaf->getTargets()->getPositions()[0];
const FReal*const posY = leaf->getTargets()->getPositions()[1];
const FReal*const posZ = leaf->getTargets()->getPositions()[2];
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 FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize indexPartOrig = indexes[idxPart];
if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3) ) {
std::cout << "Index "<< indexPartOrig <<" potential " << potentials[idxPart]
<< " Pos "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
<< " Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
}
energy += potentials[idxPart]*physicalValues[idxPart] ;
}
});
std::cout <<std::endl<<"Energy: "<< energy<<std::endl;
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
}
// -----------------------------------------------------
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.fma"));
FFmaGenericWriter<FReal> writer(name) ;
//
FSize NbPoints = loader.getNumberOfParticles();
FReal * computedParticles = new FReal[8*NbPoints]() ;
memset(computedParticles,0,8*NbPoints*sizeof(FReal));
FSize j = 0 ;
tree.forEachLeaf([&](LeafType* leaf){
//
// Input
const FReal*const posX = leaf->getTargets()->getPositions()[0];
const FReal*const posY = leaf->getTargets()->getPositions()[1];
const FReal*const posZ = leaf->getTargets()->getPositions()[2];
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
//
// Computed data
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();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
j = 8*indexes[idxPart];
computedParticles[j] = posX[idxPart] ;
computedParticles[j+1] = posY[idxPart] ;
computedParticles[j+2] = posZ[idxPart] ;
computedParticles[j+3] = physicalValues[idxPart] ;
computedParticles[j+4] = potentials[idxPart] ;
computedParticles[j+5] = forcesX[idxPart] ;
computedParticles[j+6] = forcesY[idxPart] ;
computedParticles[j+7] = forcesZ[idxPart] ;
}
});
writer.writeHeader( loader.getCenterOfBox(), loader.getBoxWidth() , NbPoints, sizeof(FReal), 8) ;
writer.writeArrayOfReal(computedParticles, 8 , NbPoints);
delete[] computedParticles;
//
std::string name1( "output.fma");
//
FFmaGenericWriter<FReal> writer1(name1) ;
writer1.writeDistributionOfParticlesFromOctree(&tree,NbPoints) ;
}
if(FParameters::existParameter(argc, argv, "-cmp") ){
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
std::cout << std::scientific;
std::cout.precision(10) ;
printf("Compute Diff...");
FMath::FAccurater<FReal> potentialDiff;
FMath::FAccurater<FReal> fx, fy, fz, f;
FReal energy = 0.0;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafType* leaf){
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 FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
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]);
f.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
f.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
f.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
energy += potentials[idxPart]*physicalValues[idxPart];
}
});
std::cout << energy << " " << energyD << std::endl;
delete[] particles;
f.setNbElements(nbParticles);
std::cout << "FChebSymKernel Energy " << FMath::Abs(energy-energyD) << " Relative "<< FMath::Abs(energy-energyD) / FMath::Abs(energyD) <<std::endl;
std::cout << " Potential Error " << potentialDiff << std::endl;
std::cout << " Fx Error " << fx << std::endl;
std::cout << " Fy Error " << fy << std::endl;
std::cout << " Fz Error " << fz << std::endl;
std::cout << " F Error " << f << std::endl;
std::cout <<std::endl<<"Energy: "<< energy<<std::endl;
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<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