Commit 5384d470 authored by berenger-bramas's avatar berenger-bramas

Try to find FMB kernel error (compare to FMB project and real application).

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@290 2616d619-271b-44dc-8df4-d4a8f33a7222
parent 37673773
......@@ -222,7 +222,7 @@ private:
// for each leafs
do{
FDEBUG(computationCounterL2P.tic());
kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
///kernels->L2P(octreeIterator.getCurrentCell(), octreeIterator.getCurrentListTargets());
FDEBUG(computationCounterL2P.tac());
// need the current particles and neighbors particles
const int counter = tree->getLeafsNeighborsWithIndex(neighbors, neighborsIndex, octreeIterator.getCurrentGlobalIndex(),heightMinusOne);
......
......@@ -40,7 +40,7 @@ class FFmbKernels : public FAbstractKernels<ParticleClass,CellClass,ContainerCla
protected:
// _GRAVITATIONAL_
static const int FMB_Info_eps_soft_square = 1;
static const int FMB_Info_eps_soft_square = 0;//1;
// Can be false or not in not blas kernels
static const int FMB_Info_up_to_P_in_M2L = true;
......@@ -1294,7 +1294,7 @@ public:
iterTarget.data().incForces( force_vector_tmp_x, force_vector_tmp_y, force_vector_tmp_z );
iterTarget.data().setPotential(expansion_Evaluate_local_with_Y_already_computed(local_exp));
iterTarget.data().incPotential(expansion_Evaluate_local_with_Y_already_computed(local_exp));
/*printf("[END] fx = %e \t fy = %e \t fz = %e \n\n",
iterTarget.data()->getForces().getX(),iterTarget.data()->getForces().getY(),iterTarget.data()->getForces().getZ());*/
......@@ -1372,10 +1372,6 @@ public:
iterSameBox.gotoNext();
}
//printf("x = %e \t y = %e \t z = %e \n",iterTarget.data()->getPosition().getX(),iterTarget.data()->getPosition().getY(),iterTarget.data()->getPosition().getZ());
//printf("\t P2P fx = %e \t fy = %e \t fz = %e \n",iterTarget.data()->getForces().getX(),iterTarget.data()->getForces().getY(),iterTarget.data()->getForces().getZ());
//printf("\t potential = %e \n",iterTarget.data()->getPotential());
iterTarget.data() = target;
iterTarget.gotoNext();
......@@ -1386,14 +1382,16 @@ public:
void DIRECT_COMPUTATION_MUTUAL_SOFT(ParticleClass& target, ParticleClass& source){
FReal dx = target.getPosition().getX() - source.getPosition().getX();
FReal dy = target.getPosition().getY() - source.getPosition().getY();
FReal dz = target.getPosition().getZ() - source.getPosition().getZ();
FReal dx = -(target.getPosition().getX() - source.getPosition().getX());
FReal dy = -(target.getPosition().getY() - source.getPosition().getY());
FReal dz = -(target.getPosition().getZ() - source.getPosition().getZ());
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz + FReal(FMB_Info_eps_soft_square));
FReal inv_distance = FMath::Sqrt(inv_square_distance);
inv_distance *= target.getPhysicalValue() * source.getPhysicalValue();
inv_square_distance *= inv_distance;
inv_square_distance *= target.getPhysicalValue() * source.getPhysicalValue();
dx *= inv_square_distance;
dy *= inv_square_distance;
......@@ -1404,15 +1402,14 @@ public:
dy,
dz
);
target.incPotential( inv_distance );
target.incPotential( inv_distance * source.getPhysicalValue() );
source.incForces(
(-dx),
(-dy),
(-dz)
);
source.incPotential( inv_distance );
source.incPotential( inv_distance * target.getPhysicalValue() );
}
......@@ -1454,10 +1451,6 @@ public:
iterSameBox.gotoNext();
}
//printf("x = %e \t y = %e \t z = %e \n",iterTarget.data()->getPosition().getX(),iterTarget.data()->getPosition().getY(),iterTarget.data()->getPosition().getZ());
//printf("\t P2P fx = %e \t fy = %e \t fz = %e \n",iterTarget.data()->getForces().getX(),iterTarget.data()->getForces().getY(),iterTarget.data()->getForces().getZ());
//printf("\t potential = %e \n",iterTarget.data()->getPotential());
iterTarget.data() = target;
iterTarget.gotoNext();
......@@ -1474,6 +1467,7 @@ public:
FReal inv_square_distance = FReal(1.0) / (dx*dx + dy*dy + dz*dz + FReal(FMB_Info_eps_soft_square));
FReal inv_distance = FMath::Sqrt(inv_square_distance);
inv_distance *= target.getPhysicalValue() * source.getPhysicalValue();
inv_square_distance *= inv_distance;
......
......@@ -39,7 +39,7 @@ int main(int argc, char ** argv){
typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FFmbKernels<ParticleClass, CellClass, ContainerClass > KernelClass;
typedef FFmmAlgorithmThread<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
///////////////////////What we do/////////////////////////////
std::cout << ">> This executable has to be used to test fmb algorithm.\n";
//////////////////////////////////////////////////////////////
......@@ -100,6 +100,9 @@ int main(int argc, char ** argv){
std::cout << "Done " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl;
{ // get sum forces&potential
FILE* fout = fopen("./res.temp.txt", "w");
FReal potential = 0;
F3DPosition forces;
typename OctreeClass::Iterator octreeIterator(&tree);
......@@ -110,11 +113,11 @@ int main(int argc, char ** argv){
potential += iter.data().getPotential() * iter.data().getPhysicalValue();
forces += iter.data().getForces();
//printf("x = %e y = %e z = %e \n",iter.data()->getPosition().getX(),iter.data()->getPosition().getY(),iter.data()->getPosition().getZ());
//printf("\t fx = %e fy = %e fz = %e \n",iter.data()->getForces().getX(),iter.data()->getForces().getY(),iter.data()->getForces().getZ());
//printf("\t\t Sum Forces ( %e , %e , %e)\n",
//forces.getX(),forces.getY(),forces.getZ());
fprintf(fout, "*** pos= (%f, %f, %f)\tv= %f\t\t\t\t\t\t\tforce= (%f, %f, %f)\t\t\t\tpotential energy=\n%f\n",
iter.data().getPosition().getX(),iter.data().getPosition().getY(),iter.data().getPosition().getZ(),
iter.data().getPhysicalValue(),
iter.data().getForces().getX(),iter.data().getForces().getY(),iter.data().getForces().getZ(),
iter.data().getPotential());
iter.gotoNext();
}
......@@ -122,6 +125,8 @@ int main(int argc, char ** argv){
std::cout << "Foces Sum x = " << forces.getX() << " y = " << forces.getY() << " z = " << forces.getZ() << std::endl;
std::cout << "Potential = " << potential << std::endl;
fclose(fout);
}
// -----------------------------------------------------
......
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