Skip to content

S
ScalFMM
Commit 96d44b04 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files
Options

Add cpmments, code dimplifications, ... But still a bug with spherical kernels

parent ee6289b7
Hide whitespace changes
Inline Side-by-side
Showing with 87 additions and 65 deletions
Src/Kernels/Rotation/FRotationKernel.hpp
View file @ 96d44b04
......@@ -319,7 +319,7 @@ class FRotationKernel : public FAbstractKernels<CellClass,ContainerClass> {
{
int index_lm = index_P0;
for(int m = 0 ; m <= P ; ++m, ++index_lm ){
const FReal mphi = (sph.getAzimuth() + FMath::FPiDiv2) * FReal(m);
const FReal mphi = (sph.getPhiInO2PI() + FMath::FPiDiv2) * FReal(m);
// O_{l,m}( \alpha, \beta + \phi ) = e^{-i \phi m} O_{l,m}( \alpha, \beta )
rotationExpMinusImPhi[idxChild][index_lm].setRealImag(FMath::Cos(-mphi), FMath::Sin(-mphi));
// M_{l,m}( \alpha, \beta + \phi ) = e^{i \phi m} M_{l,m}( \alpha, \beta )
......@@ -418,7 +418,7 @@ class FRotationKernel : public FAbstractKernels<CellClass,ContainerClass> {
{
int index_lm = index_P0;
for(int m = 0 ; m <= P ; ++m, ++index_lm ){
const FReal mphi = (sph.getAzimuth() + FMath::FPiDiv2) * FReal(m);
const FReal mphi = (sph.getPhiInO2PI() + FMath::FPiDiv2) * FReal(m);
// O_{l,m}( \alpha, \beta + \phi ) = e^{-i \phi m} O_{l,m}( \alpha, \beta )
rotationM2LExpMinusImPhi[position][index_lm].setRealImag(FMath::Cos(-mphi), FMath::Sin(-mphi));
// M_{l,m}( \alpha, \beta + \phi ) = e^{i \phi m} M_{l,m}( \alpha, \beta )
......
Src/Kernels/Rotation/FRotationOriginalKernel.hpp
View file @ 96d44b04
......@@ -508,7 +508,7 @@ public:
// rotate it forward
const FSpherical sph = getSphericalChild(inLevel, idxChild);
rotateMultipole(source_w, sph.getAzimuth(), sph.getInclination());
rotateMultipole(source_w, sph.getPhiInO2PI(), sph.getInclination());
const FReal b = -sph.getR();
// Translate it
......@@ -528,7 +528,7 @@ public:
}
// Rotate it back
deRotateMultipole(target_w, sph.getAzimuth(), sph.getInclination());
deRotateMultipole(target_w, sph.getPhiInO2PI(), sph.getInclination());
// Sum the result
FMemUtils::addall( inPole->getMultipole(), target_w, SizeArray);
......@@ -559,7 +559,7 @@ public:
// Rotate
const FSpherical sph = getSphericalInteraction(inLevel, idxNeigh);
rotateMultipole(source_w, sph.getAzimuth(), sph.getInclination());
rotateMultipole(source_w, sph.getPhiInO2PI(), sph.getInclination());
const FReal b = sph.getR();
// Transfer to u
......@@ -580,7 +580,7 @@ public:
}
// Rotate it back
deRotateTaylor(target_u, sph.getAzimuth(), sph.getInclination());
deRotateTaylor(target_u, sph.getPhiInO2PI(), sph.getInclination());
// Sum
FMemUtils::addall(inLocal->getLocal(), target_u, SizeArray);
......@@ -611,7 +611,7 @@ public:
// Rotate
const FSpherical sph = getSphericalChild(inLevel, idxChild);
rotateTaylor(source_u, sph.getAzimuth(), sph.getInclination());
rotateTaylor(source_u, sph.getPhiInO2PI(), sph.getInclination());
const FReal b = sph.getR();
// Translate
......@@ -631,7 +631,7 @@ public:
}
// Rotate
deRotateTaylor(target_u, sph.getAzimuth(), sph.getInclination());
deRotateTaylor(target_u, sph.getPhiInO2PI(), sph.getInclination());
// Sum in child
FMemUtils::addall(inChildren[idxChild]->getLocal(), target_u, SizeArray);
......
Src/Kernels/Spherical/FAbstractSphericalKernel.hpp
View file @ 96d44b04
......@@ -536,38 +536,41 @@ private:
// The -(- 1.0) computing is not the most efficient programming ...
const FReal signe = 1.0;
//if( FMath::Epsilon < spherical.getR()){
// POURQUOI diviser par spherical.getR()
force_vector_in_local_base_x = ( force_vector_in_local_base_x * signe / spherical.getR()); // ????
// force_vector_in_local_base_x = ( force_vector_in_local_base_x ); // Better
//
// The derivative wrt r is equivalent to use classical outer expansion multiplied by l and the result is divided by r
force_vector_in_local_base_x = ( force_vector_in_local_base_x * signe / spherical.getR());
// classical definition of the derivatives
force_vector_in_local_base_y = ( force_vector_in_local_base_y * signe / spherical.getR());
force_vector_in_local_base_z = ( force_vector_in_local_base_z * signe / (spherical.getR() * spherical.getSinTheta()));
//}
/////////////////////////////////////////////////////////////////////
// ???????????????????
//spherical_position_Set_ph
FReal ph = FMath::Fmod(spherical.getPhi(), FReal(2)*FMath::FPi);
if (ph > M_PI) ph -= FReal(2) * FMath::FPi;
if (ph < -M_PI + FMath::Epsilon) ph += FReal(2) * FMath::FPi;
//spherical_position_Set_th
FReal th = FMath::Fmod(spherical.getTheta(), FReal(2) * FMath::FPi);
if (th < 0.0) th += 2*FMath::FPi;
if (th > FMath::FPi){
th = 2*FMath::FPi - th;
//spherical_position_Set_ph(p, spherical_position_Get_ph(p) + M_PI);
ph = FMath::Fmod(ph + FMath::FPi, 2*FMath::FPi);
if (ph > M_PI) ph -= 2*FMath::FPi;
if (ph < -M_PI + FMath::Epsilon) ph += 2 * FMath::FPi;
}
//spherical_position_Set_r
//FReal rh = spherical.r;
//
//spherical_position_Set_ph
FReal ph = spherical.getPhiInO2PI();
// FReal ph = FMath::Fmod(spherical.getPhi(), FReal(2)*FMath::FPi);
// if (ph > M_PI) ph -= FReal(2) * FMath::FPi;
// if (ph < -M_PI + FMath::Epsilon) ph += FReal(2) * FMath::FPi;
//
// //spherical_position_Set_th
// FReal th = FMath::Fmod(spherical.getTheta(), FReal(2) * FMath::FPi);
// //FReal th = spherical.getTheta();
// if (th < 0.0) th += 2*FMath::FPi;
// if (th > FMath::FPi){
// th = 2*FMath::FPi - th;
// //spherical_position_Set_ph(p, spherical_position_Get_ph(p) + M_PI);
// ph = FMath::Fmod(ph + FMath::FPi, 2*FMath::FPi);
// if (ph > M_PI) ph -= 2*FMath::FPi;
// if (ph < -M_PI + FMath::Epsilon) ph += 2 * FMath::FPi;
// }
// //spherical_position_Set_r
// //FReal rh = spherical.r;
FAssertLF(spherical.getR() >= 0 , "R should be < 0!");
const FReal cos_theta = FMath::Cos(th);
const FReal cos_theta = spherical.getCosTheta(); //FMath::Cos(th);
// const FReal cos_theta = FMath::Cos(th);
const FReal cos_phi = FMath::Cos(ph);
const FReal sin_theta = FMath::Sin(th);
const FReal sin_phi = FMath::Sin(ph);
const FReal sin_theta = spherical.getSinTheta() ;// FMath::Sin(th);
// const FReal sin_theta = FMath::Sin(th);
const FReal sin_phi = FMath::Sin(ph);
//
//
// Formulae below are OK
......
Src/Utils/FMath.cpp
View file @ 96d44b04
......@@ -17,6 +17,7 @@
// Constant values
const FReal FMath::FPi = FReal(M_PI);
const FReal FMath::FTwoPi = FReal(2.0*M_PI);
const FReal FMath::FPiDiv2 = FReal(M_PI_2);
const FReal FMath::Epsilon = FReal(0.00000000000000000001);
......
Src/Utils/FMath.hpp
View file @ 96d44b04
......@@ -19,7 +19,6 @@
#include <cmath>
#include <limits>
#include <limits>
#include "FGlobal.hpp"
......@@ -32,6 +31,7 @@
*/
struct FMath{
static const FReal FPi; //< Pi constant
static const FReal FTwoPi; //< 2 Pi constant
static const FReal FPiDiv2; //< Pi/2 constant
static const FReal Epsilon; //< Epsilon
......@@ -134,7 +134,8 @@ struct FMath{
return log2(inValue);
}
/** To get atan2 of a 2 FReal */
/** To get atan2 of a 2 FReal, The return value is given in radians and is in the
range -pi to pi, inclusive. */
static float Atan2(const float inValue1,const float inValue2){
return atan2f(inValue1,inValue2);
}
......@@ -158,10 +159,11 @@ struct FMath{
return cos(inValue);
}
/** To get acos of a FReal */
/** To get arccos of a float. The result is in the range [0, pi]*/
static float ACos(const float inValue){
return acosf(inValue);
}
/** To get arccos of a double. The result is in the range [0, pi]*/
static double ACos(const double inValue){
return acos(inValue);
}
......@@ -170,6 +172,7 @@ struct FMath{
static float Fmod(const float inValue1,const float inValue2){
return fmodf(inValue1,inValue2);
}
/** return the floating-point remainder of inValue1 / inValue2 */
static double Fmod(const double inValue1,const double inValue2){
return fmod(inValue1,inValue2);
}
......
UTests/utestRotationDirectDebug.cpp
View file @ 96d44b04
......@@ -102,13 +102,13 @@ class TestRotationDirect : public FUTester<TestRotationDirect> {
// Run FMM
Print("Fmm...");
// Print("Fmm...");
KernelClass kernels(NbLevels,boxWidth, boxCenter);
FmmClass algo(&tree,&kernels);
algo.execute();
// Run direct computation
Print("Direct...");
// Print("Direct...");
for(int idxTarget = 0 ; idxTarget < nbParticles ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < nbParticles ; ++idxOther){
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
......@@ -123,7 +123,7 @@ class TestRotationDirect : public FUTester<TestRotationDirect> {
}
// Compare
Print("Compute Diff...");
// Print("Compute Diff...");
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
{ // Check that each particle has been summed with all other
......@@ -143,9 +143,13 @@ class TestRotationDirect : public FUTester<TestRotationDirect> {
fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
std::cout << indexPartOrig << " x " << particles[indexPartOrig].forces[0]<<" "<<forcesX[idxPart] << std::endl;
std::cout << indexPartOrig << " y " << particles[indexPartOrig].forces[1]<<" "<<forcesY[idxPart] << std::endl;
std::cout << indexPartOrig << " z " << particles[indexPartOrig].forces[2]<<" "<<forcesZ[idxPart] << std::endl;
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;
printf("Printf : forces x %e y %e z %e\n", forcesX[idxPart],forcesY[idxPart],forcesZ[idxPart]);//TODO delete
// std::cout << indexPartOrig << " x " << particles[indexPartOrig].forces[0]<<" "<<forcesX[idxPart] << std::endl;
// std::cout << indexPartOrig << " y " << particles[indexPartOrig].forces[1]<<" "<<forcesY[idxPart] << std::endl;
// std::cout << indexPartOrig << " z " << particles[indexPartOrig].forces[2]<<" "<<forcesZ[idxPart] << std::endl;
}
});
......@@ -217,7 +221,7 @@ class TestRotationDirect : public FUTester<TestRotationDirect> {
// The tests!
///////////////////////////////////////////////////////////
static const int P = 5;
static const int P = 2;
/** Rotation */
void TestRotation(){
......
UTests/utestSphericalDirectDebug.cpp
View file @ 96d44b04
......@@ -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 = 5;
const int DevP = 2;
// Warning in make test the exec dir it Build/UTests
// Load particles
const int nbParticles = 2;
......@@ -88,21 +88,22 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
Print(nbParticles);
for(int idxLeafX = 0 ; idxLeafX < dimGrid ; ++idxLeafX){
/*for(int idxLeafY = 0 ; idxLeafY < dimGrid ; ++idxLeafY)*/{
/*for(int idxLeafZ = 0 ; idxLeafZ < dimGrid ; ++idxLeafZ)*/{
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 + quarterDimLeaf,
// FReal(idxLeafY)*dimLeaf + quarterDimLeaf,
// FReal(idxLeafZ)*dimLeaf + quarterDimLeaf);
/*particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + quarterDimLeaf,
particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
FReal(idxLeafY)*dimLeaf + quarterDimLeaf,
FReal(idxLeafZ)*dimLeaf + quarterDimLeaf);
/* particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + quarterDimLeaf,
2*quarterDimLeaf,
2*quarterDimLeaf);*/
/*particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
quarterDimLeaf,
quarterDimLeaf);*/
2*quarterDimLeaf);
particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
3*quarterDimLeaf,
quarterDimLeaf,
quarterDimLeaf);
particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
quarterDimLeaf,
quarterDimLeaf);*/
// Create octree
OctreeClass tree(NbLevels, SizeSubLevels, boxWidth, boxCenter);
......@@ -207,16 +208,26 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
Print(fz.getInfNorm());
// Assert
const FReal MaximumDiff = FReal(0.0001);
if(fx.getL2Norm() > MaximumDiff || fx.getInfNorm() > MaximumDiff){
std::cout << "Error in X " << fx.getL2Norm() << " " << fx.getInfNorm() << std::endl;
}
if(fy.getL2Norm() > MaximumDiff || fy.getInfNorm() > MaximumDiff){
std::cout << "Error in Y " << fy.getL2Norm() << " " << fy.getInfNorm() << std::endl;
}
if(fz.getL2Norm() > MaximumDiff || fz.getInfNorm() > MaximumDiff){
std::cout << "Error in Z " << fz.getL2Norm() << " " << fz.getInfNorm() << std::endl;
}
// Assert
const FReal MaximumDiff = FReal(0.0001);
uassert(potentialDiff.getL2Norm() < MaximumDiff);
uassert(potentialDiff.getInfNorm() < MaximumDiff);
uassert(fx.getL2Norm() < MaximumDiff);
uassert(fx.getInfNorm() < MaximumDiff);
uassert(fy.getL2Norm() < MaximumDiff);
uassert(fy.getInfNorm() < MaximumDiff);
uassert(fz.getL2Norm() < MaximumDiff);
uassert(fz.getInfNorm() < MaximumDiff);
// const FReal MaximumDiff = FReal(0.0001);
// if(fx.getL2Norm() > MaximumDiff || fx.getInfNorm() > MaximumDiff){
// std::cout << "Error in X " << fx.getL2Norm() << " " << fx.getInfNorm() << std::endl;
// }
// if(fy.getL2Norm() > MaximumDiff || fy.getInfNorm() > MaximumDiff){
// std::cout << "Error in Y " << fy.getL2Norm() << " " << fy.getInfNorm() << std::endl;
// }
// if(fz.getL2Norm() > MaximumDiff || fz.getInfNorm() > MaximumDiff){
// std::cout << "Error in Z " << fz.getL2Norm() << " " << fz.getInfNorm() << 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