// ===================================================================================
// ===================================================================================
// 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".
// ===================================================================================
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include "ScalFmmConfig.h"
#include "../../Src/Utils/FTic.hpp"
#include "../../Src/Utils/FMath.hpp"
#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Utils/FIOVtk.hpp"
#include "../../Src/Containers/FOctree.hpp"
#include "../../Src/Containers/FVector.hpp"
#include "../../Src/Core/FFmmAlgorithmPeriodic.hpp"
#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Files/FDlpolyLoader.hpp"
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "../../Src/Kernels/Spherical/FSphericalKernel.hpp"
#include "../../Src/Kernels/Spherical/FSphericalCell.hpp"
#ifdef ScalFMM_USE_BLAS
// chebyshev kernel
#include "../../Src/Kernels/Chebyshev/FChebCell.hpp"
#include "../../Src/Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "../../Src/Kernels/Chebyshev/FChebKernel.hpp"
#include "../../Src/Kernels/Chebyshev/FChebSymKernel.hpp"
#endif
/** DLpoly particle is used in the gadget program
* here we try to make the same simulation
*/
struct MDParticle {
FPoint position;
FReal forces[3];
FReal physicalValue;
FReal potential;
int index;
};
// Simply create particles and try the kernels
int main(int argc, char ** argv){
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
#ifdef ScalFMM_USE_BLAS
// begin Chebyshev kernel
// accuracy
const unsigned int ORDER = 12;
// typedefs
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
#else
typedef FSphericalCell CellClass;
typedef FOctree< CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FSphericalKernel< CellClass, ContainerClass > KernelClass;
const int DevP = FParameters::getValue(argc,argv,"-P", 9);
#endif
typedef FFmmAlgorithmPeriodic<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
// typedef FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClassNoPer;
///////////////////////What we do/////////////////////////////
if( FParameters::existParameter(argc, argv, "-help")){
std::cout << ">> This executable has to be used to compute direct interaction either for periodic or non periodic system.\n";
std::cout << ">> options are -h H -sh SH [-per perdeep, -noper] -fin filenameIN (-bin) -fout filenameOUT \n";
std::cout << ">> Recommended files : ../Data/EwalTest_Periodic.run ../Data/EwalTest_NoPeriodic.run\n";
std::cout << " Options " << std::endl;
std::cout << " -per perDeep " << std::endl;
std::cout << " -noper no periodic boundary conditions " << std::endl;
std::cout << " -verbose : print index x y z fx fy fy Q and V" << std::endl;
std::cout << " -fout filenameOUT bunary output file " << std::endl;
exit(-1);
}
if(FParameters::existParameter(argc, argv, "-per") &&FParameters::existParameter(argc, argv, "-noper") ){
std::cerr <<" Error -per X and -noper are forbidden together " << std::endl;
exit(-1);
}
//////////////////////////////////////////////////////////////
const int NbLevels = FParameters::getValue(argc,argv,"-h", 4);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 2);
const int PeriodicDeep = FParameters::getValue(argc,argv,"-per", 3);
const std::string filenameIn(FParameters::getStr(argc,argv,"-fin", "../Data/forceNacl_128_dlpolyPer.bin"));
const char* const filenameOut = FParameters::getStr(argc,argv,"-fout", "periodicDirect.out");
// file for -saveError option
std::ofstream errorfile("outputEwaldError.txt",std::ios::out);
FTic counter;
// -----------------------------------------------------
// LOADER
// -----------------------------------------------------
std::cout << "Opening : " << filenameIn << "\n";
FDlpolyLoader *loader = nullptr ;
if(FParameters::existParameter(argc, argv, "-bin")){
loader = new FDlpolyBinLoader(filenameIn.c_str());
}
else {
loader = new FDlpolyAsciiLoader(filenameIn.c_str());
}
if(! loader->isOpen()){
std::cout << "Loader Error, " << filenameIn << " is missing\n";
return 1;
}
OctreeClass tree(NbLevels, SizeSubLevels, loader->getBoxWidth(), loader->getCenterOfBox());
// ---------------------------------------------------------------------------------
// Insert particles in the Octree
// --------------------------------------------------------------------------------- std::cout << "Creating & Inserting " << loader->getNumberOfParticles() << " particles ..." << std::endl;
std::cout << "\tWidth : " << loader->getBoxWidth() << " \t center x : " << loader->getCenterOfBox().getX()
<< " y : " << loader->getCenterOfBox().getY() << " z : " << loader->getCenterOfBox().getZ() << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
counter.tic();
FPoint electricMoment(0.0,0.0,0.0) ;
// const --> then shared
MDParticle * const particles = new MDParticle[loader->getNumberOfParticles()];
memset(particles, 0, sizeof(MDParticle) * loader->getNumberOfParticles()) ;
MDParticle* particlesDirect = nullptr;
particlesDirect = new MDParticle[loader->getNumberOfParticles()];
memset(particlesDirect, 0, sizeof(MDParticle) * loader->getNumberOfParticles()) ;
//
int nbParticles = static_cast<int>(loader->getNumberOfParticles());
double totalCharge = 0.0;
//
for(int idxPart = 0 ; idxPart < loader->getNumberOfParticles() ; ++idxPart){
//
loader->fillParticle(&particles[idxPart].position, particles[idxPart].forces,
&particles[idxPart].physicalValue,&particles[idxPart].index);
particlesDirect[idxPart].position = particles[idxPart].position;
particlesDirect[idxPart].index = particles[idxPart].index ;
particlesDirect[idxPart].physicalValue = particles[idxPart].physicalValue;
//
totalCharge += particles[idxPart].physicalValue ;
//
electricMoment.incX(particles[idxPart].physicalValue*particles[idxPart].position.getX() );
electricMoment.incY(particles[idxPart].physicalValue*particles[idxPart].position.getY() );
electricMoment.incZ(particles[idxPart].physicalValue*particles[idxPart].position.getZ() );
//
// reset forces and insert in the tree
//
tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
}
counter.tac();
std::cout << std::endl;
std::cout << "Total Charge = "<< totalCharge <<std::endl;
std::cout << "Electric Moment = "<< electricMoment <<std::endl;
std::cout << "Electric Moment norm = "<< electricMoment.norm2() <<std::endl;
std::cout << std::endl;
std::cout << std::endl;
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << " s)." << std::endl;
//
//
// ---------------------------------------------------------------------------------
FTreeCoordinate min, max;
if( FParameters::existParameter(argc, argv, "-per") ){
FmmClass algo(&tree,PeriodicDeep);
algo.repetitionsIntervals(&min, &max);
}
// ----------------------------------------------------------------------------------------------------------
// DIRECT COMPUTATION
// ----------------------------------------------------------------------------------------------------------
FReal denergy = 0.0;
//
// direct computation
//
{
printf("Compute direct:\n");
printf("Box [%d;%d][%d;%d][%d;%d]\n", min.getX(), max.getX(), min.getY(),
max.getY(), min.getZ(), max.getZ());
//
//
// particles : Input results
// particlesDirect : Direct results
//
counter.tic();
denergy = 0.0 ;
//#pragma omp parallel for shared(loader,nbParticles,min, max,particlesDirect) schedule(static,1)
#pragma omp parallel shared(loader,nbParticles,min, max,particlesDirect)
{
#pragma omp for schedule(static,1)
for(int idxTarget = 0 ; idxTarget < nbParticles ; ++idxTarget){
//
// compute with all other except itself
//
// Compute force and potential between particles[idxTarget] and particles inside the box
//
for(int idxOther = 0; idxOther < nbParticles ; ++idxOther){
if( idxOther != idxTarget ){
FP2P::NonMutualParticles(
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
particlesDirect[idxTarget].position.getX(), particlesDirect[idxTarget].position.getY(),
particlesDirect[idxTarget].position.getZ(),particlesDirect[idxTarget].physicalValue,
&particlesDirect[idxTarget].forces[0],&particlesDirect[idxTarget].forces[1],
&particlesDirect[idxTarget].forces[2],&particlesDirect[idxTarget].potential);
}
}
//
// compute with image boxes
//
for(int idxX = min.getX() ; idxX <= max.getX() ; ++idxX){
for(int idxY = min.getY() ; idxY <= max.getY() ; ++idxY){
for(int idxZ = min.getZ() ; idxZ <= max.getZ() ; ++idxZ){
{
if(idxX == 0 && idxY == 0 && idxZ == 0) continue;
const FPoint offset(loader->getBoxWidth() * FReal(idxX),
loader->getBoxWidth() * FReal(idxY),
loader->getBoxWidth() * FReal(idxZ));
//
for(int idxSource = 0 ; idxSource < nbParticles ; ++idxSource){
MDParticle source = particles[idxSource];
source.position += offset;
FP2P::NonMutualParticles(
source.position.getX(), source.position.getY(),source.position.getZ(),source.physicalValue,
particlesDirect[idxTarget].position.getX(), particlesDirect[idxTarget].position.getY(),particlesDirect[idxTarget].position.getZ(),particlesDirect[idxTarget].physicalValue,
&particlesDirect[idxTarget].forces[0],&particlesDirect[idxTarget].forces[1],&particlesDirect[idxTarget].forces[2],&particlesDirect[idxTarget].potential
);
}
}
}
}
} // End nested loop - compute with image boxes
} // end for
// Compute the energy
//#pragma omp for shared(nbParticles,particlesDirect) reduction(+:denergy)
#pragma omp for reduction(+:denergy)
for(int idx = 0 ; idx < nbParticles ; ++idx){
denergy += particlesDirect[idx].potential*particlesDirect[idx].physicalValue ;
}
} // end pragma parallel
denergy *= 0.5 ;
counter.tac();
//
printf("Energy DIRECT= %.14e\n",denergy);
//
std::cout << "Done " << "(@DIRECT Algorithm = " << counter.elapsed() << " s)." << std::endl;
std::cout << "\n"<< "END DIRECT "
<< "-------------------------------------------------------------------------"
<< std::endl << std::endl ;
} // END DIRECT
//
// ----------------------------------------------------------------
// Save direct computation in binary format
// write binary output file
//
std::cout << "Generate " << filenameOut <<" from input file" << std::endl;
std::fstream fileout(filenameOut,std::ifstream::in|std::ifstream::out| std::ios::binary| std::ios::trunc);
if(!fileout) {
std::cout << "Cannot open file."<< std::endl;
return 1;
}
//
std::cout << " nbParticles: " << nbParticles <<" " << sizeof(nbParticles) <<std::endl;
std::cout << " denergy: " << denergy <<" " << sizeof(denergy) <<std::endl;
std::cout << " Box size: " << loader->getBoxWidth() << " " << sizeof(loader->getBoxWidth())<<std::endl;
double boxsize[3] ; boxsize[0] = boxsize[1]= boxsize[2]=loader->getBoxWidth();
int PER[4] ;
if( FParameters::existParameter(argc, argv, "-noper") ){
PER[0] = 0 ;
PER[1] = PER[2] = PER[3] = 0 ;
}
else {
PER[0] = 1 ;
PER[1] = PER[2] = PER[3] = PeriodicDeep ;
}
fileout.write((char*)&nbParticles,sizeof(nbParticles));
fileout.write((char*)&boxsize,sizeof(double)*3);
fileout.write((char*)&PER,sizeof(int)*4);
fileout.write((char*)&denergy,sizeof(denergy));
//
fileout.write ((char*)&particlesDirect[0], sizeof(MDParticle)*nbParticles);
fileout.flush();
//
// end generate
// -----------------------------------------------------
//
if(FParameters::existParameter(argc, argv, "-verbose")){
denergy = 0 ;
for(int idx = 0 ; idx < nbParticles ; ++idx){
std::cout << ">> index " << particlesDirect[idx].index << std::endl;
std::cout << " x " << particlesDirect[idx].position.getX() << " y " << particlesDirect[idx].position.getY() << " z " << particlesDirect[idx].position.getZ() << std::endl;
std::cout << " fx " << particlesDirect[idx].forces[0] << " fy " << particlesDirect[idx].forces[1] << " fz " << particlesDirect[idx].forces[2] << std::endl;
std::cout << " Q " << particlesDirect[idx].physicalValue << " V " << particlesDirect[idx].potential << std::endl;
std::cout << "\n";
denergy += particlesDirect[idx].potential*particlesDirect[idx].physicalValue ;
}
}
std::cout << " ENERGY " << denergy*0.5 << std::endl;
delete[] particles;
delete[] particlesDirect;
return 0;
}