Commit 88282e57 authored by Florent Pruvost's avatar Florent Pruvost
Browse files
parents 0f9d37d0 2be6b19c
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas
// 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_FFT
// ==== Git =====
// @SCALFMM_PRIVATE
// ================
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <string>
#include "ScalFmmConfig.h"
#include "Files/FFmaGenericLoader.hpp"
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "Utils/FParameters.hpp"
#include "Containers/FOctree.hpp"
#ifdef _OPENMP
#include "Core/FFmmAlgorithmThread.hpp"
#else
#include "Core/FFmmAlgorithm.hpp"
#endif
#include "Utils/FParameterNames.hpp"
/**
* This program runs the FMM Algorithm with the uniform interpolation kernel and compares the results with a direct computation.
*/
/// \file LagrangeInterpolationFMM.cpp
//!
//! \brief This program runs the FMM Algorithm with the interpolation kernel based on uniform (grid points) interpolation (1/r kernel)
//! \authors B. Bramas, O. Coulaud
//!
//! This code is a short example to use the Interpolation approach for the 1/r kernel
// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
FHelpDescribeAndExit(argc, argv,
"Driver for Lagrange interpolation kernel (1/r kernel).",
FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,FParameterDefinitions::OutputFile,
FParameterDefinitions::NbThreads);
const std::string defaultFile(/*SCALFMMDataPath+*/"Data/unitCubeXYZQ100.bfma" );
const std::string filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFile.c_str());
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, 1);
#ifdef _OPENMP
omp_set_num_threads(NbThreads);
std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
#else
std::cout << "\n>> Sequential version.\n" << std::endl;
#endif
//
std::cout << "Parameters "<< std::endl
<< " Octree Depth "<< TreeHeight <<std::endl
<< " SubOctree depth " << SubTreeHeight <<std::endl
<< " Input file name: " <<filename <<std::endl
<< " Thread number: " << NbThreads <<std::endl
<<std::endl;
//
// init timer
FTic time;
// open particle file
////////////////////////////////////////////////////////////////////
//
FFmaGenericLoader loader(filename);
//
////////////////////////////////////////////////////////////////////
// begin Lagrange kernel
// accuracy
const unsigned int ORDER = 7;
// typedefs
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FUnifCell<ORDER> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
//
typedef FInterpMatrixKernelR MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
typedef FUnifKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
//
#ifdef _OPENMP
typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#else
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#endif
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
{ // -----------------------------------------------------
std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
<< " particles ..." << std::endl;
std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
time.tic();
//
FPoint position;
FReal physicalValue = 0.0;
//
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
//
// Read particle per particle from file
loader.fillParticle(&position,&physicalValue);
//
// put particle in octree
tree.insert(position, idxPart, physicalValue);
}
time.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = "
<< time.elapsed() << " s) ." << std::endl;
} // -----------------------------------------------------
{ // -----------------------------------------------------
std::cout << "\nLagrange FMM (ORDER="<< ORDER << ") ... " << std::endl;
time.tic();
//
KernelClass *kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
//
FmmClass algo(&tree, kernels);
//
algo.execute(); // Here the call of the FMM algorithm
//
time.tac();
std::cout << "Timers Far Field \n"
<< "P2M " << algo.getTime(FAlgorithmTimers::P2MTimer) << " seconds\n"
<< "M2M " << algo.getTime(FAlgorithmTimers::M2MTimer) << " seconds\n"
<< "M2L " << algo.getTime(FAlgorithmTimers::M2LTimer) << " seconds\n"
<< "L2L " << algo.getTime(FAlgorithmTimers::L2LTimer) << " seconds\n"
<< "P2P and L2P " << algo.getTime(FAlgorithmTimers::NearTimer) << " seconds\n"
<< std::endl;
std::cout << "Done " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
}
// -----------------------------------------------------
//
// Some output
//
//
{ // -----------------------------------------------------
long int 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([&](LeafClass* 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 int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FVector<int>& indexes = leaf->getTargets()->getIndexes();
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const int 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 writer(name) ;
//
int NbPoints = loader.getNumberOfParticles();
FReal * particles ;
particles = new FReal[8*NbPoints] ;
memset(particles,0,8*NbPoints*sizeof(FReal));
int j = 0 ;
tree.forEachLeaf([&](LeafClass* 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<int>& 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 int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
j = 8*indexes[idxPart];
particles[j] = posX[idxPart] ;
particles[j+1] = posY[idxPart] ;
particles[j+2] = posZ[idxPart] ;
particles[j+3] = physicalValues[idxPart] ;
particles[j+4] = potentials[idxPart] ;
particles[j+5] = forcesX[idxPart] ;
particles[j+6] = forcesY[idxPart] ;
particles[j+7] = forcesZ[idxPart] ;
}
});
writer.writeHeader( loader.getCenterOfBox(), loader.getBoxWidth() , NbPoints, sizeof(FReal), 8) ;
writer.writeArrayOfReal(particles, 8 , NbPoints);
}
return 0;
}
......@@ -18,6 +18,8 @@
#include "Utils/FPoint.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FGenerateDistribution.hpp"
#include "Files/FExportWriter.hpp"
#include "Utils/FParameterNames.hpp"
//
......@@ -35,7 +37,6 @@
//! \param -fout name: generic name for files (without extension) and save data
//! with following format in name.fma or name.bfma if -bin is set"
//! \param -bin save output in binary mode (name file name.bfma
//! \param -visufmt format for the visu file (vtk, vtp, cvs or cosmo). vtp is the default
//!
//!
//! \b examples
......@@ -44,6 +45,7 @@
//!
//! changeFormat -fin unitCubeXYZQ100.fma -fout unitCubeXYZQ100 -bin
// \param -visufmt format for the visu file (vtk, vtp, cvs or cosmo). vtp is the default
int main(int argc, char ** argv){
......@@ -51,8 +53,7 @@ int main(int argc, char ** argv){
"Driver to change the format of the input file. "
"fdlpoly is not supported for now.",
FParameterDefinitions::InputFile, FParameterDefinitions::OutputFile,
FParameterDefinitions::OutputVisuFile,FParameterDefinitions::FormatVisuFile,
FParameterDefinitions::OutputBinFormat);
FParameterDefinitions::OutputVisuFile);
FSize NbPoints;
FReal * particles = nullptr ;
......@@ -86,19 +87,8 @@ int main(int argc, char ** argv){
// Generate file for ScalFMM FMAGenericLoader
//
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output"));
std::string ext(".");
if(name.find(ext) !=std::string::npos) {
std::cout << "No file with extension permitted for output name : " << name << std::endl;
exit(-1);
}
if( FParameters::existParameter(argc, argv, FParameterDefinitions::OutputBinFormat.options)){
name += ".bfma";
}
else {
name += ".fma";
}
FFmaGenericWriter writer(name) ;
std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.fma"));
FFmaGenericWriter writer(name) ;
writer.writeHeader( loader->getCenterOfBox(), loader->getBoxWidth() , NbPoints, sizeof(FReal), nbData) ;
writer.writeArrayOfReal(particles, nbData, NbPoints);
}
......@@ -106,53 +96,8 @@ int main(int argc, char ** argv){
// Generate file for visualization purpose
//
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputVisuFile.options)){
std::string outfilename(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output"));
std::string visufile(""), fmt(FParameters::getStr(argc,argv,FParameterDefinitions::OutputVisuFile.options, "vtp"));
if( fmt == "vtp" ){
visufile = outfilename + ".vtp" ;
}
else if( fmt == "vtk" ){
visufile = outfilename + ".vtk" ;
}
else if( fmt == "cosmo" ){
if(nbData !=4) {
std::cerr << "Cosmos export accept only 4 data per particles. here: "<<nbData<<std::endl;
std::exit(EXIT_FAILURE);
}
visufile = outfilename + ".cosmo" ;
}
else {
visufile = outfilename + ".csv" ;
}
std::ofstream file( visufile, std::ofstream::out);
if(!file) {
std::cout << "Cannot open file."<< std::endl;
exit(-1) ;
} //
//
// Export data in cvs format
//
if( fmt == "vtp" ){
std::cout << "Writes in XML VTP format (visualization) in file "<< visufile <<std::endl ;
if(nbData==4){
exportVTKxml( file, particles, NbPoints) ;
}
else {
exportVTKxml( file, particles, NbPoints,nbData) ;
}
}
else if( fmt == "vtk" ){
std::cout << "Writes in VTK format (visualization) in file "<< visufile <<std::endl ;
exportVTK( file, particles, NbPoints,nbData) ;
}
else if( fmt == "cosmo" ){
std::cout << "Writes in COSMO format (visualization) in file "<< visufile <<std::endl ;
exportCOSMOS( file, particles, NbPoints) ;
}
else {
std::cout << "Writes in CVS format (visualization) in file "<<visufile<<std::endl ;
exportCVS( file, particles, NbPoints,nbData) ;
}
std::string outfilename(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.vtp"));
driverExportData(outfilename, particles , NbPoints,loader->getNbRecordPerline() );
}
//
delete particles ;
......
......@@ -16,6 +16,7 @@
#include "Utils/FPoint.hpp"
#include "Utils/FGenerateDistribution.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "Files/FExportWriter.hpp"
#include "Utils/FParameterNames.hpp"
......@@ -74,6 +75,10 @@
int main(int argc, char ** argv){
const FParameterNames LocalOptionEllipsoid = {
{"-ellipsoid"} ,
" non uniform distribution on an ellipsoid of aspect ratio given by \n -ar a:b:c with a, b and c > 0\n"
};
FHelpDescribeAndExit(argc, argv,
">> Driver to generate N points (non)uniformly distributed on a given geometry.\n"
"Options \n"
......@@ -97,12 +102,13 @@ int main(int argc, char ** argv){
" -radius R - default value 10.0\n"
" Physical values\n"
" -charge generate physical values between -1 and 1 otherwise generate between 0 and 1 \n"
" -zeromean the average of the physical values is zero \n"
" Output \n"
" -filename name: generic name for files (without extension) and save data\n"
" with following format in name.fma or name.bfma in -bin is set\n"
" -visufmt vtk, vtp, cosmo or cvs format.",
FParameterDefinitions::InputFile, FParameterDefinitions::OutputBinFormat, FParameterDefinitions::NbParticles);
" -zeromean the average of the physical values is zero \n",
// " Output \n"
// " -filename name: generic name for files (without extension) and save data\n"
// " with following format in name.fma or name.bfma in -bin is set\n"
// " -visufmt vtk, vtp, cosmo or cvs format.",
FParameterDefinitions::OutputFile,
FParameterDefinitions::NbParticles,FParameterDefinitions::OutputVisuFile,LocalOptionEllipsoid);
......@@ -188,14 +194,16 @@ int main(int argc, char ** argv){
std::cout << "A: "<<A<<" B "<< B << " C: " << C<<std::endl;
unifRandonPointsOnProlate(NbPoints,A,C,particles);
BoxWith = 2.0*C;
}
} //const int NbPoints = FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options, 20000);
else if(FParameters::existParameter(argc, argv, "-ellipsoid")){
// else if(FParameters::existParameter(argc, argv, "-ellipsoid")){
std::string dd(":"),aspectRatio = FParameters::getStr(argc,argv,"-ar", "1:1:2");
// std::string dd(":"),aspectRatio = FParameters::getStr(argc,argv,"-ar", "1:1:2");
FReal A,B,C ;
size_t pos = aspectRatio.find(":"); aspectRatio.replace(pos,1," ");
pos = aspectRatio.find(":"); aspectRatio.replace(pos,1," ");
std::stringstream ss(aspectRatio); ss >>A >> B >> C ;
// std::cout << "A: "<<A<<" B "<< B << " C: " << C<<std::endl;
std::cout << "A: "<<A<<" B "<< B << " C: " << C<<std::endl;
nonunifRandonPointsOnElipsoid(NbPoints,A,B,C,particles);
BoxWith = 2.0*FMath::Max( A,FMath::Max( B,C)) ;
}
......@@ -220,14 +228,7 @@ int main(int argc, char ** argv){
BoxWith += 2*extraRadius ;
}
std::string name(genericFileName);
if( FParameters::existParameter(argc, argv, FParameterDefinitions::OutputBinFormat.options)){
name += ".bfma";
}
else {
name += ".fma";
}
std::cout << "Write "<< NbPoints <<" Particles in file " << name <<std::endl;
std::cout << "Write "<< NbPoints <<" Particles in file " << name << std::endl;
FFmaGenericWriter writer(name) ;
writer.writeHeader(Centre,BoxWith, NbPoints, *ppart) ;
writer.writeArrayOfParticles(ppart, NbPoints);
......@@ -237,43 +238,8 @@ int main(int argc, char ** argv){
// Generate file for visualization
//
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputVisuFile.options)){
std::string visufile(""), fmt(FParameters::getStr(argc,argv,FParameterDefinitions::OutputVisuFile.options, "vtp"));
if( fmt == "vtp" ){
visufile = genericFileName + ".vtp" ;
}
else if( fmt == "vtk" ){
visufile = genericFileName + ".vtk" ;
}
else if( fmt == "cosmo" ){
visufile = genericFileName + ".cosmo" ;
}
else {
visufile = genericFileName + ".csv" ;
}
std::ofstream file( visufile.c_str(), std::ofstream::out);
if(!file) {
std::cout << "Cannot open file."<< std::endl;
exit(-1) ;
} //
//
// Export data in VTK format
//
if( fmt == "vtp" ){
std::cout << "Writes in XML VTP format (visualization) in file "<< visufile <<std::endl ;
exportVTKxml( file, particles, NbPoints) ;
}
else if( fmt == "vtk" ){
std::cout << "Writes in VTK format (visualization) in file "<< visufile <<std::endl ;
exportVTK( file, particles, NbPoints) ;
}
else if( fmt == "cosmo" ){
std::cout << "Writes in COSMO format (visualization) in file "<< visufile <<std::endl ;
exportCOSMOS( file, particles, NbPoints) ;
}
else {
std::cout << "Writes in CVS format (visualization) in file "<<visufile<<std::endl ;
exportCVS( file, particles, NbPoints) ;
}
std::string visufile(FParameters::getStr(argc,argv,FParameterDefinitions::OutputVisuFile.options, "output.vtp"));
driverExportData(visufile, particles , NbPoints);
}
//
delete [] particles ;
......
......@@ -6,10 +6,18 @@ Copyright (c) 2011-2014 Inria, All rights reserved.
This file contains the main features as well as overviews of specific
bug fixes (and other actions) for each version of ScalFMM since
version 1.1
1.3-
-----
- Add interpolation FMM based on uniform grid points.
- Add blocked version of the algorithm to increase the granularity (task-based approach)
-
1.2.1
-----
- Bug fix : Support for huge MPI message in tree construction and Parallel QuickSort (count can be greater than Int32.MaxValue)
- Bug fix : Data sharing attribute clauses for omp in Core/FAlgorithmThreadProc.hpp
1.2-
1.2
-----
- New FMA format to read/write particles
- Add examples repository
......@@ -23,7 +31,3 @@ version 1.1
- Add SSE and AVX code for 1/r kernel
- CMake improvements
1.2.1
-----
- Bug fix : Support for huge MPI message in tree construction and Parallel QuickSort (count can be greater than Int32.MaxValue)
- Bug fix : Data sharing attribute clauses for omp in Core/FAlgorithmThreadProc.hpp
......@@ -117,9 +117,9 @@ public:
// To get access to descriptor
friend struct FTestCellDescriptor;
friend struct FTestCell_Alignement;
};
#endif //FTESTCELL_HPP
// ===================================================================================
// Copyright ScalFmm 2014 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".
// ===================================================================================
#ifndef FEXPORTDATA_HPP
#define FEXPORTDATA_HPP
// @author O. Coulaud
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <string>
//! \fn void exportCVS(std::ofstream& file, const FReal * particles , const int N, const int nbDataPerParticle=4)
//! \brief Export particles in CVS Format
//!
//! Export particles in CVS Format as follow
//! x , y , z , physicalValue, P, FX, FY, FY
//! It is useful to plot the distribution with paraView
//!
//! @param file stream to save the data
//! @param N number of particles
//! @param particles array of particles of type FReal (float or double) Its size is N*nbDataPerParticle
//! @param nbDataPerParticle number of values per particles (default value 4)
//!
void exportCVS(std::ofstream& file, const FReal * particles , const int N, const int nbDataPerParticle=4){
int j = 0;
if (nbDataPerParticle==4){
file << " x , y , z, q " <<std::endl;
}
else {
file << " x , y , z, q P FX FY FZ" <<std::endl;
}
for(int i = 0 ; i< N; ++i, j+=nbDataPerParticle){
file << particles[j] ;
for (int k = 1 ; k< nbDataPerParticle ; ++k) {
file << " , " << particles[j+k] ;
}
file << std::endl;
}
}
//
//! \fn void exportCOSMOS(std::ofstream& file, const FReal * particles, const int N )
//! \brief Export particles in COSMOS Format
//!
//! Export particles in CVS Format as follow
//! x , y , z , 0.0, 0.0, 0.0, physicalValue
//!
//! @param file stream to save the data
//! @param particles array of particles of type FReal (float or double) Its size is 4*N (X,Y,Z,M)
//! @param N number of particles
//!
void exportCOSMOS(std::ofstream& file, const FReal * particles , const int N){