Commit 56780480 authored by COULAUD Olivier's avatar COULAUD Olivier

IMprovenments for Utests

parent 93b0c92e
...@@ -23,19 +23,19 @@ ...@@ -23,19 +23,19 @@
#include <cstdio> #include <cstdio>
#include <cstdlib> #include <cstdlib>
#include "../../Src/Files/FFmaGenericLoader.hpp" #include "Utils/FParameters.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "../../Src/Kernels/Rotation/FRotationKernel.hpp" #include "Kernels/Rotation/FRotationKernel.hpp"
#include "../../Src/Kernels/Rotation/FRotationCell.hpp" #include "Kernels/Rotation/FRotationCell.hpp"
#include "../../Src/Components/FSimpleLeaf.hpp" #include "Components/FSimpleLeaf.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp" #include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Containers/FOctree.hpp" #include "Containers/FOctree.hpp"
#include "../../Src/Core/FFmmAlgorithmThread.hpp" #include "Core/FFmmAlgorithmThread.hpp"
/// \file RotationFMM.cpp /// \file RotationFMM.cpp
......
This diff is collapsed.
// =================================================================================== // ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner // Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr // olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM. // This software is a computer program whose purpose is to compute the FMM.
// //
...@@ -39,7 +39,7 @@ int main(void){\ ...@@ -39,7 +39,7 @@ int main(void){\
* Please refer to testUTest.cpp to see an example * Please refer to testUTest.cpp to see an example
* @warning Create a derived class that implement SetTests() and use TestClass() macro * @warning Create a derived class that implement SetTests() and use TestClass() macro
* *
* We recommand to have a look to a unit test to better understand how it works, * We recommend to have a look to a unit test to better understand how it works,
* as for example @example TestList * as for example @example TestList
*/ */
template <class TestClass> template <class TestClass>
...@@ -61,8 +61,8 @@ class FUTester{ ...@@ -61,8 +61,8 @@ class FUTester{
int currentTest; //< current processing test in the run int currentTest; //< current processing test in the run
int currentStep; //< current processing step in the run int currentStep; //< current processing step in the run
int failledSteps; //< number of failled step in the current test int failledSteps; //< number of failed step in the current test
int failledTests; //< number of failled tests int failledTests; //< number of failed tests
protected: protected:
/** Constructor */ /** Constructor */
...@@ -94,7 +94,7 @@ protected: ...@@ -94,7 +94,7 @@ protected:
*/ */
void AddTest(TestFunc inFunc){ void AddTest(TestFunc inFunc){
char buff[256]; char buff[256];
sprintf(buff,"Unamed Test number %d",totalTests+1); sprintf(buff,"Unnamed Test number %d",totalTests+1);
AddTest(inFunc,buff); AddTest(inFunc,buff);
} }
...@@ -119,11 +119,11 @@ protected: ...@@ -119,11 +119,11 @@ protected:
void Print(const Output& value){ void Print(const Output& value){
std::cout<< "--- Output from program : " << value << "\n"; std::cout<< "--- Output from program : " << value << "\n";
} }
/** /**
* To test * To test
* @param result the test result * @param result the test result
* if result is false test failled * if result is false test failed
*/ */
void uassert(const bool result){ void uassert(const bool result){
++currentStep; ++currentStep;
......
// =================================================================================== // ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner // Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr // olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM. // This software is a computer program whose purpose is to compute the FMM.
// //
...@@ -59,10 +59,10 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> { ...@@ -59,10 +59,10 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
// Warning in make test the exec dir it Build/UTests // Warning in make test the exec dir it Build/UTests
// Load particles // Load particles
const int NbLevels = 3; const int NbLevels = 4;
const int SizeSubLevels = 2; const int SizeSubLevels = 2;
const int PeriodicDeep = 1; const int PeriodicDeep = 2;
const int NbParticles = 6; const int NbParticles = 100;
FRandomLoader loader(NbParticles); FRandomLoader loader(NbParticles);
...@@ -78,19 +78,23 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> { ...@@ -78,19 +78,23 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
FReal physicalValue; FReal physicalValue;
FReal potential; FReal potential;
}; };
FReal coeff = -1.0, value = 0.10, sum = 0.0;
TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()]; TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){ for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
FPoint position; FPoint position;
loader.fillParticle(&position); loader.fillParticle(&position);
// put in tree // put in tree
tree.insert(position, idxPart, 0.10); value *= coeff ;
sum += value ;
// put in tree
tree.insert(position, idxPart, value);
// get copy // get copy
particles[idxPart].position = position; particles[idxPart].position = position;
particles[idxPart].physicalValue = 0.10; particles[idxPart].physicalValue = value;
particles[idxPart].potential = 0.0; particles[idxPart].potential = 0.0;
particles[idxPart].forces[0] = 0.0; particles[idxPart].forces[0] = 0.0;
particles[idxPart].forces[1] = 0.0; particles[idxPart].forces[1] = 0.0;
particles[idxPart].forces[2] = 0.0; particles[idxPart].forces[2] = 0.0;
} }
// Run FMM // Run FMM
...@@ -117,7 +121,8 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> { ...@@ -117,7 +121,8 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
particles[idxOther].position.getZ(),particles[idxOther].physicalValue, particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
&particles[idxOther].forces[0],&particles[idxOther].forces[1], &particles[idxOther].forces[0],&particles[idxOther].forces[1],
&particles[idxOther].forces[2],&particles[idxOther].potential); &particles[idxOther].forces[2],&particles[idxOther].potential);
}
}
for(int idxX = min.getX() ; idxX <= max.getX() ; ++idxX){ for(int idxX = min.getX() ; idxX <= max.getX() ; ++idxX){
for(int idxY = min.getY() ; idxY <= max.getY() ; ++idxY){ for(int idxY = min.getY() ; idxY <= max.getY() ; ++idxY){
for(int idxZ = min.getZ() ; idxZ <= max.getZ() ; ++idxZ){ for(int idxZ = min.getZ() ; idxZ <= max.getZ() ; ++idxZ){
...@@ -149,10 +154,13 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> { ...@@ -149,10 +154,13 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
Print("Compute Diff..."); Print("Compute Diff...");
FMath::FAccurater potentialDiff; FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz; FMath::FAccurater fx, fy, fz;
FReal energy= 0.0 , energyD = 0.0 ;
{ // Check that each particle has been summed with all other { // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){ tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const potentials = leaf->getTargets()->getPotentials(); const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const forcesX = leaf->getTargets()->getForcesX(); const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY(); const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ(); const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
...@@ -165,6 +173,8 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> { ...@@ -165,6 +173,8 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]); fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]); fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]); fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
energy += potentials[idxPart]*physicalValues[idxPart];
energyD +=particles[indexPartOrig].potential*particles[indexPartOrig].potential;
} }
}); });
} }
...@@ -172,31 +182,40 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> { ...@@ -172,31 +182,40 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
delete[] particles; delete[] particles;
// Print for information // Print for information
Print("Potential diff is = "); Print("Potential diff is = ");
Print(potentialDiff.getL2Norm()); printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm());
Print(potentialDiff.getInfNorm()); printf(" RMSError %e\n",potentialDiff.getRMSError());
Print("Fx diff is = "); Print("Fx diff is = ");
Print(fx.getL2Norm()); printf(" L2Norm %e\n",fx.getRelativeL2Norm());
Print(fx.getInfNorm()); printf(" RMSError %e\n",fx.getRMSError());
Print(fx.getRelativeL2Norm());
Print(fx.getRelativeInfNorm());
Print("Fy diff is = "); Print("Fy diff is = ");
Print(fy.getL2Norm()); printf(" L2Norm %e\n",fy.getRelativeL2Norm());
Print(fy.getInfNorm()); printf(" RMSError %e\n",fy.getRMSError());
Print("Fz diff is = "); Print("Fz diff is = ");
Print(fz.getL2Norm()); printf(" L2Norm %e\n",fz.getRelativeL2Norm());
Print(fz.getInfNorm()); printf(" RMSError %e\n",fz.getRMSError());
FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm() + fz.getRelativeL2Norm() *fz.getRelativeL2Norm() );
// Assert printf(" Total L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
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
//
const FReal MaximumDiffPotential = FReal(9e-4); const FReal MaximumDiffPotential = FReal(9e-4);
const FReal MaximumDiffForces = FReal(9e-3); const FReal MaximumDiffForces = FReal(9e-3);
uassert(potentialDiff.getL2Norm() < MaximumDiffPotential); uassert(potentialDiff.getL2Norm() < MaximumDiffPotential); //1
uassert(potentialDiff.getInfNorm() < MaximumDiffPotential); uassert(potentialDiff.getRMSError() < MaximumDiffPotential); //2
uassert(fx.getL2Norm() < MaximumDiffForces); uassert(fx.getL2Norm() < MaximumDiffForces); //3
uassert(fx.getInfNorm() < MaximumDiffForces); uassert(fx.getRMSError() < MaximumDiffForces); //4
uassert(fy.getL2Norm() < MaximumDiffForces); uassert(fy.getL2Norm() < MaximumDiffForces); //5
uassert(fy.getInfNorm() < MaximumDiffForces); uassert(fy.getRMSError() < MaximumDiffForces); //6
uassert(fz.getL2Norm() < MaximumDiffForces); uassert(fz.getL2Norm() < MaximumDiffForces); //8
uassert(fz.getInfNorm() < MaximumDiffForces); uassert(fz.getRMSError() < MaximumDiffForces); //8
} }
...@@ -220,8 +239,8 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> { ...@@ -220,8 +239,8 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
/** TestChebKernel */ /** TestChebKernel */
void TestChebKernel(){ void TestChebKernel(){
const unsigned int ORDER = 5; const unsigned int ORDER = 6;
const FReal epsilon = FReal(1e-5); const FReal epsilon = FReal(1e-6);
typedef FP2PParticleContainerIndexed<> ContainerClass; typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass; typedef FSimpleLeaf<ContainerClass> LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass; typedef FInterpMatrixKernelR MatrixKernelClass;
...@@ -230,6 +249,7 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> { ...@@ -230,6 +249,7 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass; typedef FChebKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FFmmAlgorithmPeriodic<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass; typedef FFmmAlgorithmPeriodic<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// run test // run test
std::cout <<" TEST 1 "<<std::endl;
RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(epsilon); RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(epsilon);
} }
...@@ -245,6 +265,7 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> { ...@@ -245,6 +265,7 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass; typedef FChebSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FFmmAlgorithmPeriodic<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass; typedef FFmmAlgorithmPeriodic<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// run test // run test
std::cout <<std::endl<<" TEST 2 "<<std::endl;
RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(epsilon); RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(epsilon);
} }
...@@ -256,8 +277,8 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> { ...@@ -256,8 +277,8 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
/** set test */ /** set test */
void SetTests(){ void SetTests(){
AddTest(&TestChebyshevDirect::TestChebKernel,"Test Chebyshev Kernel with one big SVD"); AddTest(&TestChebyshevDirect::TestChebKernel,"Test Chebyshev Kernel with one big SVD") ;
AddTest(&TestChebyshevDirect::TestChebSymKernel,"Test Chebyshev Kernel with 16 small SVDs and symmetries"); AddTest(&TestChebyshevDirect::TestChebSymKernel,"Test Chebyshev Kernel with 16 small SVDs and symmetries");
} }
}; };
......
...@@ -33,7 +33,7 @@ ...@@ -33,7 +33,7 @@
class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> { class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> {
/** Here we test only the P2P */ /** Here we test only the P2P */
void TestPeriodicFmm(){ void TestPeriodicFmm(){
static const int P = 9; static const int P = 18;
typedef FRotationCell<P> CellClass; typedef FRotationCell<P> CellClass;
typedef FP2PParticleContainerIndexed<> ContainerClass; typedef FP2PParticleContainerIndexed<> ContainerClass;
...@@ -45,9 +45,9 @@ class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> { ...@@ -45,9 +45,9 @@ class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> {
typedef FFmmAlgorithmPeriodic<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass; typedef FFmmAlgorithmPeriodic<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
// Parameters // Parameters
const int NbLevels = 4; const int NbLevels = 4;
const int SizeSubLevels = 2; const int SizeSubLevels = 2;
const int PeriodicDeep = -1; const int PeriodicDeep = 2;
const int NbParticles = 100; const int NbParticles = 100;
FRandomLoader loader(NbParticles); FRandomLoader loader(NbParticles);
...@@ -58,26 +58,32 @@ class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> { ...@@ -58,26 +58,32 @@ class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> {
FReal physicalValue; FReal physicalValue;
FReal potential; FReal potential;
}; };
FReal coeff = -1.0, value = 0.10, sum = 0.0;
TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()]; TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){ for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
FPoint position; FPoint position;
loader.fillParticle(&position); loader.fillParticle(&position);
value *= coeff ;
sum += value ;
// put in tree // put in tree
tree.insert(position, idxPart, 0.10); tree.insert(position, idxPart, value);
// get copy // get copy
particles[idxPart].position = position; particles[idxPart].position = position;
particles[idxPart].physicalValue = 0.10; particles[idxPart].physicalValue = value;
particles[idxPart].potential = 0.0; particles[idxPart].potential = 0.0;
particles[idxPart].forces[0] = 0.0; particles[idxPart].forces[0] = 0.0;
particles[idxPart].forces[1] = 0.0; particles[idxPart].forces[1] = 0.0;
particles[idxPart].forces[2] = 0.0; particles[idxPart].forces[2] = 0.0;
} }
if (FMath::Abs(sum)> 0.00001){
std::cerr << "Sum of charges is not equal zero!!! (sum=<<"<<sum<<" )"<<std::endl;
exit(-1);
}
// Run FMM // Run FMM
Print("Fmm..."); Print("Fmm...");
FmmClass algo(&tree,PeriodicDeep); FmmClass algo(&tree,PeriodicDeep);
KernelClass kernels( algo.extendedTreeHeight(), algo.extendedBoxWidth(), algo.extendedBoxCenter()); KernelClass *kernels = new KernelClass( algo.extendedTreeHeight(), algo.extendedBoxWidth(), algo.extendedBoxCenter());
algo.setKernel(&kernels); algo.setKernel(kernels);
algo.execute(); algo.execute();
// Run Direct // Run Direct
...@@ -127,11 +133,13 @@ class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> { ...@@ -127,11 +133,13 @@ class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> {
Print("Compute Diff..."); Print("Compute Diff...");
FMath::FAccurater potentialDiff; FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz; FMath::FAccurater fx, fy, fz;
{ // Check that each particle has been summed with all other FReal energy = 0.0 ;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){ tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const potentials = leaf->getTargets()->getPotentials(); const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX(); const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY(); const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ(); const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles(); const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
...@@ -143,38 +151,40 @@ class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> { ...@@ -143,38 +151,40 @@ class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> {
fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]); fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]); fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]); fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
} energy += potentials[idxPart]*physicalValues[idxPart];
}
}); });
} }
Print("Potential diff is = "); Print("Potential diff is = ");
printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm()); printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm());
printf(" RMSError %e\n",potentialDiff.getRMSError()); printf(" RMSError %e\n",potentialDiff.getRMSError());
Print("Fx diff is = "); Print("Fx diff is = ");
printf(" L2Norm %e\n",fx.getRelativeL2Norm()); printf(" L2Norm %e\n",fx.getRelativeL2Norm());
printf(" RMSError %e\n",fx.getRMSError()); printf(" RMSError %e\n",fx.getRMSError());
Print(fx.getRelativeL2Norm()); Print(fx.getRelativeL2Norm());
Print(fx.getRelativeInfNorm()); Print(fx.getRelativeInfNorm());
Print("Fy diff is = "); Print("Fy diff is = ");
printf(" L2Norm %e\n",fy.getRelativeL2Norm()); printf(" L2Norm %e\n",fy.getRelativeL2Norm());
printf(" RMSError %e\n",fy.getRMSError()); printf(" RMSError %e\n",fy.getRMSError());
Print("Fz diff is = "); Print("Fz diff is = ");
printf(" L2Norm %e\n",fz.getRelativeL2Norm()); printf(" L2Norm %e\n",fz.getRelativeL2Norm());
printf(" RMSError %e\n",fz.getRMSError()); printf(" RMSError %e\n",fz.getRMSError());
FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm() + fz.getRelativeL2Norm() *fz.getRelativeL2Norm() ); FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm() + fz.getRelativeL2Norm() *fz.getRelativeL2Norm() );
printf(" L2 Force Error= %e\n",FMath::Sqrt(L2error)) ; printf(" Total L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
const FReal MaximumDiff = FReal(0.0001); const FReal MaximumDiffPotential = FReal(9e-4);
const FReal MaximumDiffForces = FReal(9e-3);
uassert(potentialDiff.getRelativeL2Norm() < MaximumDiff); // 1 uassert(potentialDiff.getL2Norm() < MaximumDiffPotential); //1
uassert(potentialDiff.getRMSError() < MaximumDiff); // 2 uassert(potentialDiff.getRMSError() < MaximumDiffPotential); //2
uassert(fx.getRelativeL2Norm() < MaximumDiff); // 3 uassert(fx.getL2Norm() < MaximumDiffForces); //3
uassert(fx.getRMSError() < MaximumDiff); // 4 uassert(fx.getRMSError() < MaximumDiffForces); //4
uassert(fy.getRelativeL2Norm() < MaximumDiff); // 5 uassert(fy.getL2Norm() < MaximumDiffForces); //5
uassert(fy.getRMSError() < MaximumDiff); // 6 uassert(fy.getRMSError() < MaximumDiffForces); //6
uassert(fz.getRelativeL2Norm() < MaximumDiff); // 7 uassert(fz.getL2Norm() < MaximumDiffForces); //8
uassert(fz.getRMSError() < MaximumDiff); // 8 uassert(fz.getRMSError() < MaximumDiffForces); //8
delete[] particles; delete[] particles;
} }
......
...@@ -46,11 +46,11 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic> ...@@ -46,11 +46,11 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
typedef FFmmAlgorithmPeriodic<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass; typedef FFmmAlgorithmPeriodic<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
// Parameters // Parameters
const int NbLevels = 3; const int NbLevels = 4;
const int SizeSubLevels = 2; const int SizeSubLevels = 2;
const int PeriodicDeep = 3; const int PeriodicDeep = 1;
const int DevP = 9; const int DevP = 14;
const int NbParticles = 1; const int NbParticles = 100;
FSphericalCell::Init(DevP); FSphericalCell::Init(DevP);
...@@ -62,15 +62,18 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic> ...@@ -62,15 +62,18 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
FReal physicalValue; FReal physicalValue;
FReal potential; FReal potential;
}; };
TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()]; FReal coeff = -1.0, value = 0.10, sum = 0.0;
TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){ for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
FPoint position; FPoint position;
loader.fillParticle(&position); loader.fillParticle(&position);
value *= coeff ;
sum += value ;
// put in tree // put in tree
tree.insert(position, idxPart, 0.10); tree.insert(position, idxPart, value);
// get copy // get copy
particles[idxPart].position = position; particles[idxPart].position = position;
particles[idxPart].physicalValue = 0.10; particles[idxPart].physicalValue = value;
particles[idxPart].potential = 0.0; particles[idxPart].potential = 0.0;
particles[idxPart].forces[0] = 0.0; particles[idxPart].forces[0] = 0.0;
particles[idxPart].forces[1] = 0.0; particles[idxPart].forces[1] = 0.0;
...@@ -131,7 +134,7 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic> ...@@ -131,7 +134,7 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
Print("Compute Diff..."); Print("Compute Diff...");
FMath::FAccurater potentialDiff; FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz; FMath::FAccurater fx, fy, fz;
FReal energy = 0.0 ; FReal energy= 0.0 , energyD = 0.0 ;
{ // Check that each particle has been summed with all other { // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){ tree.forEachLeaf([&](LeafClass* leaf){
...@@ -149,29 +152,30 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic> ...@@ -149,29 +152,30 @@ class TestSphericalDirectPeriodic : public FUTester<TestSphericalDirectPeriodic>
fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]); fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]); fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]); fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
energy += potentials[idxPart]*physicalValues[idxPart]; energy += potentials[idxPart]*physicalValues[idxPart];
energyD +=particles[indexPartOrig].potential*particles[indexPartOrig].potential;
} }
}); });
} }
Print("Potential diff is = "); Print("Potential diff is = ");
printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm()); printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm());
printf(" RMSError %e\n",potentialDiff.getRMSError()); printf(" RMSError %e\n",potentialDiff.getRMSError());
Print("Fx diff is = "); Print("Fx diff is = ");
printf(" L2Norm %e\n",fx.getRelativeL2Norm()); printf(" L2Norm %e\n",fx.getRelativeL2Norm());
printf(" RMSError %e\n",fx.getRMSError()); printf(" RMSError %e\n",fx.getRMSError());
Print(fx.getRelativeL2Norm()); Print(fx.getRelativeL2Norm());
Print(fx.getRelativeInfNorm()); Print(fx.getRelativeInfNorm());
Print("Fy diff is = "); Print("Fy diff is = ");
printf(" L2Norm %e\n",fy.getRelativeL2Norm()); printf(" L2Norm %e\n",fy.getRelativeL2Norm());
printf(" RMSError %e\n",fy.getRMSError()); printf(" RMSError %e\n",fy.getRMSError());
Print("Fz diff is = "); Print("Fz diff is = ");
printf(" L2Norm %e\n",fz.getRelativeL2Norm()); printf(" L2Norm %e\n",fz.getRelativeL2Norm());
printf(" RMSError %e\n",fz.getRMSError()); printf(" RMSError %e\n",fz.getRMSError());
FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm() + fz.getRelativeL2Norm() *fz.getRelativeL2Norm() ); FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm() + fz.getRelativeL2Norm() *fz.getRelativeL2Norm() );
printf(" L2 Force Error= %e\n",FMath::Sqrt(L2error)) ; printf(" Total L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
printf(" Energy = %.12e\n",energy); printf(" Energy Error = %.12e\n",FMath::Abs(energy-energyD));
const FReal MaximumDiff = FReal(0.001); const FReal MaximumDiff = FReal(0.001);
uassert(potentialDiff.getRelativeL2Norm() < MaximumDiff); // 1 uassert(potentialDiff.getRelativeL2Norm() < MaximumDiff); // 1
......
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