Attention une mise à jour du service Gitlab va être effectuée le mardi 18 janvier (et non lundi 17 comme annoncé précédemment) entre 18h00 et 18h30. Cette mise à jour va générer une interruption du service dont nous ne maîtrisons pas complètement la durée mais qui ne devrait pas excéder quelques minutes.

Commit 8d61751c authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files

Test for the periodic case

parent a50a6eaf
......@@ -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