diff --git a/Src/Core/FFmmAlgorithm.hpp b/Src/Core/FFmmAlgorithm.hpp index f03447e421a214c2c5f98ecf1cc0366f043f0b25..8465c9cecab76ec0de3c3757bc33216e15c3ce0d 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 2b756fcd9395633238f67b0d7ce315d32267f8d6..c4782b6fef9651e9fd0b4b3bf0c2fcdd31b8279f 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 725f1c4c366fb0f7b20ac745484c6fa7e26c18f4..d7e2993fc910f0bfe72d8cebac31c855b2ffb26a 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); } // -----------------------------------------------------