Commit 22bd92f6 authored by COULAUD Olivier's avatar COULAUD Olivier

Add debug information for spherical kernels

parent bbeed7db
......@@ -35,10 +35,11 @@
class FEwalLoader : public FAbstractLoader {
protected:
std::ifstream file; //< The file to read
FPoint centerOfBox; //< The center of box read from file
FReal boxWidth; //< the box width read from file
int nbParticles; //< the number of particles read from file
int levcfg ; //< DL_POLY CONFIG file key. 0,1 or 2
FPoint centerOfBox; //< The center of box read from file
FReal boxWidth; //< the box width read from file
int nbParticles; //< the number of particles read from file
int levcfg ; //< DL_POLY CONFIG file key. 0,1 or 2
FReal energy ;
public:
enum Type{
OW,
......@@ -66,7 +67,7 @@ public:
int imcon ;
//int tempi(0);
FReal tempf(0);
file >> levcfg >> imcon >> this->nbParticles;
file >> levcfg >> imcon >> this->nbParticles >> this->energy;
// Periodic case
if( imcon > 0 ) {
FReal widthx, widthy, widthz;
......@@ -83,7 +84,7 @@ public:
this->centerOfBox.setPosition(0.0,0.0,0.0);
}
else {
this->boxWidth = 0;
this->boxWidth = 0;
this->nbParticles = 0;
}
}
......@@ -127,7 +128,7 @@ public:
}
FReal getEnergy() const{
return static_cast<FReal>(0.0);
return this->energy;
}
/**
......@@ -232,7 +233,7 @@ public:
centerOfBox.setPosition(0.0,0.0,0.0);
}
else {
this->boxWidth = 0;
this->boxWidth = 0;
this->nbParticles = 0;
}
}
......
......@@ -16,6 +16,7 @@
#ifndef FABSTRACTSPHERICALKERNEL_HPP
#define FABSTRACTSPHERICALKERNEL_HPP
#include <iostream>
#include "../../Components/FAbstractKernels.hpp"
#include "../../Utils/FGlobal.hpp"
......@@ -491,6 +492,11 @@ private:
const FSpherical spherical(particlePosition - local_position);
harmonic.computeInnerTheta( spherical );
// std::cout << " L2P:"<<std::endl
// << " Centre: " << local_position <<std::endl
// << " PArt: " << particlePosition <<std::endl
// << " Diff " << particlePosition - local_position<<std::endl
// << " Spherical Diff " << spherical <<std::endl ;
int index_j_k = 1;
......
......@@ -119,7 +119,7 @@ public:
*/
template <class StreamClass>
friend StreamClass& operator<<(StreamClass& output, const FSpherical& inPosition){
output << "(" << inPosition.getR() << ", " << inPosition.getPhi() << ", " << inPosition.getTheta() <<")";
output << "(" << inPosition.getR() << ", " << inPosition.getTheta() << ", " << inPosition.getPhi() <<")";
return output; // for multiple << operators.
}
......
......@@ -71,19 +71,19 @@ class TestRotationDirect : public FUTester<TestRotationDirect> {
TestParticle* const particles = new TestParticle[nbParticles];
particles[0].position = FPoint(2*quarterDimLeaf, quarterDimLeaf, quarterDimLeaf);
particles[0].physicalValue = 0.50;
particles[0].physicalValue = 1.0;
particles[1].position = FPoint(2*quarterDimLeaf, quarterDimLeaf, quarterDimLeaf);
particles[1].physicalValue = -0.10;
particles[1].physicalValue = 1.0;
Print("Number of particles:");
Print(nbParticles);
for(int idxLeafX = 0 ; idxLeafX < 2/*dimGrid*/ ; ++idxLeafX){
for(int idxLeafX = 0 ; idxLeafX < dimGrid ; ++idxLeafX){
for(int idxLeafY = 0 ; idxLeafY < 1/*dimGrid*/ ; ++idxLeafY){
for(int idxLeafZ = 0 ; idxLeafZ < 1/*dimGrid*/ ; ++idxLeafZ){
particles[0].position = FPoint((idxLeafX+1)*quarterDimLeaf, quarterDimLeaf, quarterDimLeaf);
particles[1].position = FPoint(FReal(2)*dimLeaf + 2*quarterDimLeaf,
// particles[0].position = FPoint((idxLeafX+1)*quarterDimLeaf, quarterDimLeaf, quarterDimLeaf);
particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + quarterDimLeaf,
FReal(idxLeafY)*dimLeaf + quarterDimLeaf,
FReal(idxLeafZ)*dimLeaf + quarterDimLeaf);
......@@ -97,8 +97,9 @@ class TestRotationDirect : public FUTester<TestRotationDirect> {
particles[idxPart].forces[0] = 0.0;
particles[idxPart].forces[1] = 0.0;
particles[idxPart].forces[2] = 0.0;
std::cout << idxPart << " " << particles[idxPart].position << std::endl;
// FSpherical PP(particles[idxPart].position) ;
std::cout << idxPart << " cart " << particles[idxPart].position << std::endl;
// std::cout << idxPart << " sph " << FSpherical(particles[idxPart].position) << std::endl;
}
......@@ -222,7 +223,7 @@ class TestRotationDirect : public FUTester<TestRotationDirect> {
// The tests!
///////////////////////////////////////////////////////////
static const int P = 1;
static const int P = 2;
/** Rotation */
void TestRotation(){
......
......@@ -48,7 +48,7 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
template < class CellClass, class ContainerClass, class KernelClass, class LeafClass,
class OctreeClass, class FmmClass>
void RunTest(const bool isBlasKernel){
const int DevP = 2;
const int DevP = 26;
// Warning in make test the exec dir it Build/UTests
// Load particles
const int nbParticles = 2;
......@@ -80,19 +80,20 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
TestParticle* const particles = new TestParticle[nbParticles];
particles[0].position = FPoint(quarterDimLeaf, quarterDimLeaf, quarterDimLeaf);
particles[0].physicalValue = 0.50;
particles[1].position = FPoint(2*quarterDimLeaf, quarterDimLeaf, quarterDimLeaf);
particles[1].physicalValue = -0.10;
particles[0].physicalValue = 1;
// particles[1].position = FPoint(2*quarterDimLeaf, quarterDimLeaf, quarterDimLeaf);
particles[1].physicalValue = -1;
Print("Number of particles:");
Print(nbParticles);
for(int idxLeafX = 0 ; idxLeafX < dimGrid ; ++idxLeafX){
for(int idxLeafY = 0 ; idxLeafY < dimGrid ; ++idxLeafY){
for(int idxLeafZ = 0 ; idxLeafZ < dimGrid ; ++idxLeafZ){
int idxLeafY = 0,idxLeafZ = 0 ;
std::cout << "\n ------ Loop starts ---"<< std::endl ;
for(int idxLeafX = 2 ; idxLeafX < dimGrid ; ++idxLeafX){
/*for(int idxLeafY = 0 ; idxLeafY < dimGrid ; ++idxLeafY)*/{
/* for(int idxLeafZ = 0 ; idxLeafZ < dimGrid ; ++idxLeafZ)*/{
//std::cout << "Shift : " << idxLeafX << " " << idxLeafY << " " << idxLeafZ << std::endl;
particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 3*quarterDimLeaf,
FReal(idxLeafY)*dimLeaf + quarterDimLeaf,
FReal(idxLeafZ)*dimLeaf + quarterDimLeaf);
/* particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + quarterDimLeaf,
......@@ -162,30 +163,32 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
std::cout << indexPartOrig << " x Direct " << particles[indexPartOrig].forces[0]<<" FMM "<<forcesX[idxPart] << std::endl;
std::cout << indexPartOrig << " y Direct " << particles[indexPartOrig].forces[1]<<" FMM "<<forcesY[idxPart] << std::endl;
std::cout << indexPartOrig << " z Direct " << particles[indexPartOrig].forces[2]<<" FMM "<<forcesZ[idxPart] << std::endl;
std::cout <<"indx: " << indexPartOrig << " Potential Direct " << particles[indexPartOrig].potential<<" FMM "<<potentials[idxPart] << std::endl;
std::cout <<"indx: " << indexPartOrig << " Fx Direct " << particles[indexPartOrig].forces[0]<<" FMM "<<forcesX[idxPart] << std::endl;
std::cout <<"indx: " << indexPartOrig << " Fy Direct " << particles[indexPartOrig].forces[1]<<" FMM "<<forcesY[idxPart] << std::endl;
std::cout <<"indx: " << indexPartOrig << " Fz Direct " << particles[indexPartOrig].forces[2]<<" FMM "<<forcesZ[idxPart] << std::endl;
printf("Printf : forces x %e y %e z %e\n", forcesX[idxPart],forcesY[idxPart],forcesZ[idxPart]);//TODO delete
}
});
tree.forEachCell([&](CellClass* cell){
std::cout << "Multipole:\n";
std::cout << "cell " << cell->getMortonIndex() << "\n Multipole: \n";
int index_j_k = 0;
for (int j = 0 ; j <= DevP ; ++j ){
std::cout <<"[" << j << "]\n";
std::cout <<"[" << j << "] ----- level\n";
for (int k=0; k<=j ;++k, ++index_j_k){
std::cout << "[" << k << "] " << cell->getMultipole()[index_j_k].getReal() << " i" << cell->getMultipole()[index_j_k].getImag() << " ";
std::cout << "[" << k << "] ( " << cell->getMultipole()[index_j_k].getReal() << " , " << cell->getMultipole()[index_j_k].getImag() << " ) ";
}
std::cout << "\n";
}
std::cout << "\n";
std::cout << "Local:\n";
std::cout << " Local:\n";
index_j_k = 0;
for (int j = 0 ; j <= DevP ; ++j ){
std::cout <<"[" << j << "]\n";
std::cout <<"[" << j << "] ----- level \n";
for (int k=0; k<=j ;++k, ++index_j_k){
std::cout << "[" << k << "] " << cell->getLocal()[index_j_k].getReal() << " i" << cell->getLocal()[index_j_k].getImag() << " ";
std::cout << "[" << k << "] ( " << cell->getLocal()[index_j_k].getReal() << " , " << cell->getLocal()[index_j_k].getImag() << " ) ";
}
std::cout << "\n";
}
......
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