Commit 9db8d6e0 authored by berenger-bramas's avatar berenger-bramas

Remove a bug in M2L precompute in ALL spherical kernels.

git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@426 2616d619-271b-44dc-8df4-d4a8f33a7222
parent 834d15a0
......@@ -748,6 +748,7 @@ private:
((L_j_l.getImag() * I_l_j__l_k.getReal()) +
(L_j_l.getReal() * I_l_j__l_k.getImag())));
}
}//n
......
......@@ -54,7 +54,7 @@ protected:
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){
if(FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
const F3DPosition relativePos( FReal(idxX) * treeWidthAtLevel , FReal(idxY) * treeWidthAtLevel , FReal(idxZ) * treeWidthAtLevel );
const F3DPosition relativePos( FReal(-idxX) * treeWidthAtLevel , FReal(-idxY) * treeWidthAtLevel , FReal(-idxZ) * treeWidthAtLevel );
blasHarmonic.computeOuter(FSpherical(relativePos));
FComplexe* FRestrict fillTransfer = &preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)];
......@@ -121,7 +121,7 @@ public:
/** M2L with a cell and all the existing neighbors */
void M2L(CellClass* const FRestrict pole, const CellClass* distantNeighbors[343],
const int /*size*/, const int inLevel) {
const int /*size*/, const int inLevel) {
// For all neighbors compute M2L
for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
if( distantNeighbors[idxNeigh] ){
......
......@@ -69,7 +69,7 @@ protected:
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){
if(FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
const F3DPosition relativePos( FReal(idxX) * treeWidthAtLevel , FReal(idxY) * treeWidthAtLevel , FReal(idxZ) * treeWidthAtLevel );
const F3DPosition relativePos( FReal(-idxX) * treeWidthAtLevel , FReal(-idxY) * treeWidthAtLevel , FReal(-idxZ) * treeWidthAtLevel );
blasHarmonic.computeOuter(FSpherical(relativePos));
FComplexe* FRestrict fillTransfer = &preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)];
......
......@@ -48,7 +48,7 @@ protected:
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){
if(FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
const F3DPosition relativePos( FReal(idxX) * treeWidthAtLevel , FReal(idxY) * treeWidthAtLevel , FReal(idxZ) * treeWidthAtLevel );
const F3DPosition relativePos( FReal(-idxX) * treeWidthAtLevel , FReal(-idxY) * treeWidthAtLevel , FReal(-idxZ) * treeWidthAtLevel );
Parent::harmonic.computeOuter(FSpherical(relativePos));
FMemUtils::copyall<FComplexe>(&preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)], Parent::harmonic.result(), Parent::harmonic.getExpSize());
}
......
......@@ -366,7 +366,7 @@ protected:
new (&preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)]) RotationM2LTransfer(Parent::devP,devM2lP,Parent::harmonic.getExpSize());
if(FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
const F3DPosition relativePos( FReal(idxX) * treeWidthAtLevel , FReal(idxY) * treeWidthAtLevel , FReal(idxZ) * treeWidthAtLevel );
const F3DPosition relativePos( FReal(-idxX) * treeWidthAtLevel , FReal(-idxY) * treeWidthAtLevel , FReal(-idxZ) * treeWidthAtLevel );
preM2LTransitions[idxLevel + Parent::periodicLevels][indexM2LTransition(idxX,idxY,idxZ)].transfer_M2L_rotation_Fill(FSpherical(relativePos), rotation_Info);
}
}
......
......@@ -155,16 +155,6 @@ struct FMath{
return std::isfinite(value);
}
/** Compute a relative difference between two values */
template <class ValueClass>
static ValueClass RelativeDiff(const ValueClass& value1, const ValueClass& value2){
if(Abs(value1) > Abs(value2)){
return Abs((value2 - value1) / value1);
}
else{
return Abs((value2 - value1) / value2);
}
}
};
......
This diff is collapsed.
......@@ -97,7 +97,6 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
Print("Direct...");
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
//kernels.DIRECT_COMPUTATION_MUTUAL_SOFT(particles[idxTarget], particles[idxOther]);
kernels.directInteractionMutual(&particles[idxTarget], &particles[idxOther]);
}
}
......@@ -116,25 +115,13 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
while( leafIter.hasNotFinished() ){
const ParticleClass& other = particles[leafIter.data().getIndex()];
const FReal currentPotentialDiff = FMath::RelativeDiff(other.getPotential(),leafIter.data().getPotential());
if( potentialDiff < currentPotentialDiff ){
potentialDiff = currentPotentialDiff;
}
const FReal currentFx = FMath::RelativeDiff(other.getForces().getX() , leafIter.data().getForces().getX());
if( fx < currentFx ){
fx = currentFx;
}
const FReal currentFy = FMath::RelativeDiff(other.getForces().getY() , leafIter.data().getForces().getY());
if( fy < currentFy ){
fy = currentFy;
}
const FReal currentFz = FMath::RelativeDiff(other.getForces().getZ() , leafIter.data().getForces().getZ());
if( fz < currentFz ){
fz = currentFz;
}
potentialDiff += (other.getPotential(),leafIter.data().getPotential())/other.getPotential();
fx += FMath::Abs((other.getForces().getX()-leafIter.data().getForces().getX())/other.getForces().getX());
fy += FMath::Abs((other.getForces().getY()-leafIter.data().getForces().getY())/other.getForces().getY());
fz += FMath::Abs((other.getForces().getZ()-leafIter.data().getForces().getZ())/other.getForces().getZ());
leafIter.gotoNext();
}
......
......@@ -163,25 +163,13 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
other.getForces().getX(),other.getForces().getY(),other.getForces().getZ());// todo delete
const FReal currentPotentialDiff = FMath::RelativeDiff(other.getPotential(),leafIter.data().getPotential());
if( potentialDiff < currentPotentialDiff ){
potentialDiff = currentPotentialDiff;
}
potentialDiff += (other.getPotential(),leafIter.data().getPotential())/other.getPotential();
const FReal currentFx = FMath::RelativeDiff(other.getForces().getX() , leafIter.data().getForces().getX());
if( fx < currentFx ){
fx = currentFx;
}
fx += FMath::Abs((other.getForces().getX()-leafIter.data().getForces().getX())/other.getForces().getX());
const FReal currentFy = FMath::RelativeDiff(other.getForces().getY() , leafIter.data().getForces().getY());
if( fy < currentFy ){
fy = currentFy;
}
fy += FMath::Abs((other.getForces().getY()-leafIter.data().getForces().getY())/other.getForces().getY());
const FReal currentFz = FMath::RelativeDiff(other.getForces().getZ() , leafIter.data().getForces().getZ());
if( fz < currentFz ){
fz = currentFz;
}
fz += FMath::Abs((other.getForces().getZ()-leafIter.data().getForces().getZ())/other.getForces().getZ());
leafIter.gotoNext();
}
......@@ -328,25 +316,13 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
other.getForces().getX(),other.getForces().getY(),other.getForces().getZ());// todo delete
const FReal currentPotentialDiff = FMath::RelativeDiff(other.getPotential(),leafIter.data().getPotential());
if( potentialDiff < currentPotentialDiff ){
potentialDiff = currentPotentialDiff;
}
potentialDiff += (other.getPotential(),leafIter.data().getPotential())/other.getPotential();
const FReal currentFx = FMath::RelativeDiff(other.getForces().getX() , leafIter.data().getForces().getX());
if( fx < currentFx ){
fx = currentFx;
}
fx += FMath::Abs((other.getForces().getX()-leafIter.data().getForces().getX())/other.getForces().getX());
const FReal currentFy = FMath::RelativeDiff(other.getForces().getY() , leafIter.data().getForces().getY());
if( fy < currentFy ){
fy = currentFy;
}
fy += FMath::Abs((other.getForces().getY()-leafIter.data().getForces().getY())/other.getForces().getY());
const FReal currentFz = FMath::RelativeDiff(other.getForces().getZ() , leafIter.data().getForces().getZ());
if( fz < currentFz ){
fz = currentFz;
}
fz += FMath::Abs((other.getForces().getZ()-leafIter.data().getForces().getZ())/other.getForces().getZ());
leafIter.gotoNext();
}
......
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