Commit d036973c authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files
parents 94e3dcc8 3ff16eec
This diff is collapsed.
This diff is collapsed.
......@@ -52,24 +52,30 @@ public:
char buffer[bufferSize];
file.getline(buffer, bufferSize);
int imcon ;
int tempi(0);
FReal tempf(0);
file >> tempi >> tempi >> this->nbParticles >> tempf;
FReal widthx, widthy, widthz;
file >> widthx >> tempf >> tempf;
file >> tempf >> widthy >> tempf;
file >> tempf >> tempf >> widthz;
this->centerOfBox.setPosition(0.0,0.0,0.0);
this->boxWidth = widthx;
file >> tempi >> imcon >> this->nbParticles >> tempf;
if(imcon >0 ) {
FReal widthx, widthy, widthz;
file >> widthx >> tempf >> tempf;
file >> tempf >> widthy >> tempf;
file >> tempf >> tempf >> widthz;
// this->centerOfBox.setPosition(0.0,0.0,0.0);
this->boxWidth = widthx;
}
else{ // Non periodic case.
this->boxWidth = tempf ;
}
this->centerOfBox.setPosition(0.0,0.0,0.0);
}
else {
this->boxWidth = 0;
this->nbParticles = 0;
}
std::cout << "boxWidth: "<< this->boxWidth <<std::endl ;
}
/**
* Default destructor, simply close the file
*/
......@@ -122,15 +128,17 @@ public:
FReal x, y, z, fx, fy, fz, vx, vy, vz;
int index;
char type[2];
std::string line;
file.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
file.read(type, 2);
file >> index;
std::getline(file, line); // needed to skip the end of the line in non periodic case
file >> x >> y >> z;
file >> vx >> vy >> vz;
file >> fx >> fy >> fz;
// std::cout << " x >> y >> z: " << x<< " " <<y<< " " <<z <<std::endl;
inParticle.setPosition(x,y,z);
inParticle.setForces(fx,fy,fz);
//inParticle.setForces(vx,vy,vz);
......
// ===================================================================================
// Logiciel initial: ScalFmm Version 0.5
// Co-auteurs : Olivier Coulaud, Bérenger Bramas.
// Propriétaires : INRIA.
// Copyright © 2011-2012, diffusé sous les termes et conditions d’une licence propriétaire.
// Initial software: ScalFmm Version 0.5
// Co-authors: Olivier Coulaud, Bérenger Bramas.
// Owners: INRIA.
// Copyright © 2011-2012, spread under the terms and conditions of a proprietary license.
// ===================================================================================
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>
#include "../../Src/Utils/FTic.hpp"
#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Containers/FOctree.hpp"
#include "../../Src/Containers/FVector.hpp"
#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Core/FFmmAlgorithmThread.hpp"
#include "../../Src/Core/FFmmAlgorithmTask.hpp"
#include "../../Src/Kernels/Spherical/FSphericalKernel.hpp"
#include "../../Src/Kernels/Spherical/FSphericalCell.hpp"
#include "../../Src/Kernels/Spherical/FSphericalParticle.hpp"
#include "../../Src/Components/FSimpleLeaf.hpp"
//#include "../../Src/Files/FFmaLoader.hpp"
#include "../../Src/Files/FEwalLoader.hpp"
/** This program show an example of use of
* the fmm basic algo
* it also check that each particles is little or longer
* related that each other
*/
class DLPolyParticle : public FSphericalParticle {
public:
// Type of particle
enum Type{
OW,
HW,
Undefined,
};
private:
Type type; //< current type
int index; //< current index in array
public:
// Basic constructor
DLPolyParticle() : type(Undefined), index(-1) {
}
Type getType() const{
return type;
}
void setType(const Type inType) {
type = inType;
}
int getIndex() const{
return index;
}
void setIndex( const int inIndex ){
index = inIndex;
}
};
//
// ------------------------------------------------------------------
//
// Simply create particles and try the kernels
int main(int argc, char ** argv){
typedef DLPolyParticle ParticleClass;
typedef FSphericalCell CellClass;
typedef FVector<ParticleClass> ContainerClass;
typedef FSimpleLeaf<ParticleClass, ContainerClass > LeafClass;
typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FSphericalKernel<ParticleClass, CellClass, ContainerClass > KernelClass;
typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
typedef FFmmAlgorithmThread<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClassThread;
typedef FFmmAlgorithmTask<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClassTask;
///////////////////////What we do/////////////////////////////
std::cout << ">> This executable has to be used to test Spherical algorithm.\n";
std::cout << ">> You can pass -sequential or -task (thread by default).\n";
//////////////////////////////////////////////////////////////
const int DevP = FParameters::getValue(argc,argv,"-p", 8);
const int NbLevels = FParameters::getValue(argc,argv,"-h", 5);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 3);
FTic counter;
//
const FReal coeff_MD= 138935.4835 / 418.4 ;
//
const char* const filename = FParameters::getStr(argc,argv,"-f", "../Data/test20k.fma");
std::cout << "Opening : " << filename << "\n";
FEwalLoader<ParticleClass> loader(filename);
if(!loader.isOpen()){
std::cout << "Loader Error, " << filename << " is missing\n";
return 1;
}
// -----------------------------------------------------
CellClass::Init(DevP);
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
// -----------------------------------------------------
std::cout << "Creating & Inserting " << loader.getNumberOfParticles() << " particles ..." << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
counter.tic();
loader.fillTree(tree);
counter.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << "s)." << std::endl;
// -----------------------------------------------------
std::cout << "Create kernel ..." << std::endl;
counter.tic();
KernelClass kernels(DevP, NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
counter.tac();
std::cout << "Done " << " in " << counter.elapsed() << "s)." << std::endl;
// -----------------------------------------------------
std::cout << "Working on particles ..." << std::endl;
if( FParameters::findParameter(argc,argv,"-sequential") != FParameters::NotFound){
FmmClass algo(&tree,&kernels);
counter.tic();
algo.execute();
}
else if( FParameters::findParameter(argc,argv,"-task") != FParameters::NotFound){
FmmClassTask algo(&tree,&kernels);
counter.tic();
algo.execute();
}
else {
FmmClassThread algo(&tree,&kernels);
counter.tic();
algo.execute();
}
counter.tac();
std::cout << "Done " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl;
{ // get sum forces&potential
FReal potential = 0;
FPoint forces;
OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
ContainerClass::ConstBasicIterator iter(*octreeIterator.getCurrentListTargets());
while( iter.hasNotFinished() ){
potential += iter.data().getPotential() * iter.data().getPhysicalValue();
forces += iter.data().getForces();
std::cout << " " << iter.data().getIndex()+1 << " \t "<< std::setprecision(5)<< iter.data().getPosition().getX() << " \t" << iter.data().getPosition().getY() << " \t" << iter.data().getPosition().getZ() << " Forces: \t"<< std::setprecision(8)<< iter.data().getForces().getX()*coeff_MD << " \t " << iter.data().getForces().getY()*coeff_MD << " \t " << iter.data().getForces().getZ()*coeff_MD<< std::endl;
iter.gotoNext();
}
} while(octreeIterator.moveRight());
std::cout << "Foces Sum x = " << forces.getX() << " y = " << forces.getY() << " z = " << forces.getZ() << std::endl;
std::cout << "Potential = " << potential*coeff_MD << std::endl;
}
std::cout << "Constante DL_POLY: " << coeff_MD << std::endl;
// -----------------------------------------------------
return 0;
}
......@@ -10,6 +10,7 @@
// ===================================================================================
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>
......@@ -90,6 +91,7 @@ int main(int argc, char ** argv){
const int PeriodicDeep = FParameters::getValue(argc,argv,"-per", 2);
const char* const filename = FParameters::getStr(argc,argv,"-f", "../Data/testEwal417.dt");
FTic counter;
const FReal coeff_MD= 138935.4835 / 418.4 ;
// -----------------------------------------------------
......@@ -145,6 +147,7 @@ int main(int argc, char ** argv){
// -----------------------------------------------------
{
FReal potential = 0;
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
......@@ -155,14 +158,17 @@ int main(int argc, char ** argv){
while( iter.hasNotFinished() ){
const ParticleClass& part = particles[iter.data().getIndex()];
std::cout << ">> index " << iter.data().getIndex() << " type " << iter.data().getType() << std::endl;
std::cout << "Good x " << part.getPosition().getX() << " y " << part.getPosition().getY() << " z " << part.getPosition().getZ() << std::endl;
std::cout << "FMM x " << iter.data().getPosition().getX() << " y " << iter.data().getPosition().getY() << " z " << iter.data().getPosition().getZ() << std::endl;
std::cout << "Good fx " <<part.getForces().getX() << " fy " << part.getForces().getY() << " fz " << part.getForces().getZ() << std::endl;
std::cout << "FMM fx " << iter.data().getForces().getX() << " fy " << iter.data().getForces().getY() << " fz " << iter.data().getForces().getZ() << std::endl;
std::cout << "GOOD physical value " << part.getPhysicalValue() << " potential " << part.getPotential() << std::endl;
std::cout << "FMM physical value " << iter.data().getPhysicalValue() << " potential " << iter.data().getPotential() << std::endl;
std::cout << "\n";
// std::cout << ">> index " << iter.data().getIndex() << " type " << iter.data().getType() << std::endl;
// std::cout << "Good x " << part.getPosition().getX() << " y " << part.getPosition().getY() << " z " << part.getPosition().getZ() << std::endl;
// std::cout << "FMM x " << iter.data().getPosition().getX() << " y " << iter.data().getPosition().getY() << " z " << iter.data().getPosition().getZ() << std::endl;
// std::cout << "Good fx " <<part.getForces().getX() << " fy " << part.getForces().getY() << " fz " << part.getForces().getZ() << std::endl;
// std::cout << "FMM fx " << iter.data().getForces().getX() << " fy " << iter.data().getForces().getY() << " fz " << iter.data().getForces().getZ() << std::endl;
// std::cout << "GOOD physical value " << part.getPhysicalValue() << " potential " << part.getPotential() << std::endl;
// std::cout << "FMM physical value " << iter.data().getPhysicalValue() << " potential " << iter.data().getPotential() << std::endl;
// std::cout << "\n";
potential += iter.data().getPotential() * iter.data().getPhysicalValue();
std::cout << " " << iter.data().getIndex()+1 << " \t "<< std::setprecision(5)<< iter.data().getPosition().getX() << " \t" << iter.data().getPosition().getY() << " \t" << iter.data().getPosition().getZ() << " Forces: \t"<< std::setprecision(8)<< iter.data().getForces().getX() << " \t " << iter.data().getForces().getY() << " \t " << iter.data().getForces().getZ()<< std::endl;
potentialDiff.add(part.getPotential(),iter.data().getPotential());
fx.add(part.getForces().getX(),iter.data().getForces().getX());
......@@ -172,7 +178,7 @@ int main(int argc, char ** argv){
iter.gotoNext();
}
} while(octreeIterator.moveRight());
printf("Potential diff is = \n");
printf("%f\n",potentialDiff.getL2Norm());
printf("%f\n",potentialDiff.getInfNorm());
......@@ -185,6 +191,8 @@ int main(int argc, char ** argv){
printf("Fz diff is = \n");
printf("%f\n",fz.getL2Norm());
printf("%f\n",fz.getInfNorm());
std::cout << std::endl<< std::endl<< "Potential = " << potential*coeff_MD << std::endl;
}
// -----------------------------------------------------
......
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