From 5384d470a457283e170134de1a6c3e17620bb99a Mon Sep 17 00:00:00 2001
From: berenger-bramas <berenger-bramas@2616d619-271b-44dc-8df4-d4a8f33a7222>
Date: Mon, 9 Jan 2012 08:44:40 +0000
Subject: [PATCH] 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
---
Src/Core/FFmmAlgorithm.hpp | 2 +-
Src/Fmb/FFmbKernels.hpp | 28 +++++++++++-----------------
Tests/testFmbAlgorithm.cpp | 17 +++++++++++------
3 files changed, 23 insertions(+), 24 deletions(-)
diff --git a/Src/Core/FFmmAlgorithm.hpp b/Src/Core/FFmmAlgorithm.hpp
index f03447e42..8465c9cec 100644
--- a/Src/Core/FFmmAlgorithm.hpp
+++ b/Src/Core/FFmmAlgorithm.hpp
@@ -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);
diff --git a/Src/Fmb/FFmbKernels.hpp b/Src/Fmb/FFmbKernels.hpp
index 2b756fcd9..c4782b6fe 100644
--- a/Src/Fmb/FFmbKernels.hpp
+++ b/Src/Fmb/FFmbKernels.hpp
@@ -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;
diff --git a/Tests/testFmbAlgorithm.cpp b/Tests/testFmbAlgorithm.cpp
index 725f1c4c3..d7e2993fc 100644
--- a/Tests/testFmbAlgorithm.cpp
+++ b/Tests/testFmbAlgorithm.cpp
@@ -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);
}
// -----------------------------------------------------
--
GitLab