Fix bug in the energy computation. Test works on MacOS X. But the results are...

Fix bug in the energy computation. Test works on MacOS X. But the results are strange. The force computations are more precise than the potential computation !!! We should investigate that.

UTests/utestSphericalDirectPeriodic.cpp

@@ -63,13 +63,14 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
             FReal physicalValue;
             FReal potential;
         };
-        FReal coeff = -1.0, value = 0.10, sum = 0.0;
+        FReal coeff = -1.0, value = 0.10, sum = 0.0, a= 0.0;
        TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
         for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
             FPoint position;
             loader.fillParticle(&position);
             value *= coeff ;
             sum += value ;
+            a = std::max(a,position.getX()*position.getX()+position.getY()*position.getY()+position.getZ()*position.getZ());
             // put in tree
             tree.insert(position, idxPart, value);
             // get copy
@@ -80,7 +81,7 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
             particles[idxPart].forces[1] = 0.0;
             particles[idxPart].forces[2] = 0.0;
         }
-
+        double CorErr = loader.getNumberOfParticles()*value/a;
         // Run FMM
         Print("Fmm...");
         FmmClass algo(&tree,PeriodicDeep);
@@ -153,7 +154,7 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
                     fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
                     fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
                     energy   += potentials[idxPart]*physicalValues[idxPart];
-                    energyD +=particles[indexPartOrig].potential*particles[indexPartOrig].potential;
+                    energyD +=particles[indexPartOrig].potential*particles[indexPartOrig].physicalValue;
 
                 }
             });
@@ -175,17 +176,37 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
 		printf("        RMSError %e\n",fz.getRMSError());
         FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm()  + fz.getRelativeL2Norm() *fz.getRelativeL2Norm()  );
 		printf(" Total L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
-		printf("  Energy Error =   %.12e\n",FMath::Abs(energy-energyD));
-
-        const FReal MaximumDiff = FReal(0.001);
-        uassert(potentialDiff.getRelativeL2Norm() < MaximumDiff);  // 1
-        uassert(potentialDiff.getRMSError() < MaximumDiff);  // 2
-        uassert(fx.getRelativeL2Norm()  < MaximumDiff);     // 3
-        uassert(fx.getRMSError() < MaximumDiff);                // 4
-        uassert(fy.getRelativeL2Norm()  < MaximumDiff);     // 5
-        uassert(fy.getRMSError() < MaximumDiff);                 // 6
-        uassert(fz.getRelativeL2Norm()  < MaximumDiff);    // 7
-        uassert(fz.getRMSError() < MaximumDiff);                // 8
+//
+		printf("  Energy Error  =   %.12e\n",FMath::Abs(energy-energyD));
+		printf("  Energy FMM    =   %.12e\n",FMath::Abs(energy));
+		printf("  Energy DIRECT =   %.12e\n",FMath::Abs(energyD));
+
+// ASSERT section
+		double epsilon = 1.0/FMath::pow2(DevP);
+		const FReal MaximumDiffPotential = FReal(CorErr*epsilon);
+		const FReal MaximumDiffForces     = FReal(10*CorErr*epsilon);
+		printf(" Criteria error - Epsilon  %e  \n",epsilon);
+//
+		Print("Test1 - Error Relative L2 norm Potential ");
+		uassert(potentialDiff.getRelativeL2Norm() < MaximumDiffPotential);    //1
+		Print("Test2 - Error RMS L2 norm Potential ");
+		uassert(potentialDiff.getRMSError() < MaximumDiffPotential);  //2
+		Print("Test3 - Error Relative L2 norm FX ");
+		uassert(fx.getRelativeL2Norm()  < MaximumDiffForces);                       //3
+		Print("Test4 - Error RMS L2 norm FX ");
+		uassert(fx.getRMSError() < MaximumDiffForces);                      //4
+		Print("Test5 - Error Relative L2 norm FY ");
+		uassert(fy.getRelativeL2Norm()  < MaximumDiffForces);                       //5
+		Print("Test6 - Error RMS L2 norm FY ");
+		uassert(fy.getRMSError() < MaximumDiffForces);                      //6
+		Print("Test7 - Error Relative L2 norm FZ ");
+		uassert(fz.getRelativeL2Norm()  < MaximumDiffForces);                      //8
+		Print("Test8 - Error RMS L2 norm FZ ");
+		uassert(fz.getRMSError() < MaximumDiffForces);                                           //8
+		Print("Test9 - Error Relative L2 norm F ");
+		uassert(L2error              < MaximumDiffForces);                                            //9   Total Force
+		Print("Test10 - Relative error Energy ");
+		uassert(FMath::Abs(energy-energyD) /FMath::Abs(energyD)< MaximumDiffPotential);                     //10  Total Energy
 
         delete[] particles;
     }