From 567804800f40f4059f310a49a1692b0e614b7cba Mon Sep 17 00:00:00 2001 From: Olivier Coulaud Date: Wed, 7 May 2014 11:51:10 +0200 Subject: [PATCH] IMprovenments for Utests --- Examples/RotationFMM.cpp | 16 +- Src/Utils/FPoint.hpp | 601 ++++++++++++------------ UTests/FUTester.hpp | 16 +- UTests/utestChebyshevDirectPeriodic.cpp | 91 ++-- UTests/utestRotationDirectPeriodic.cpp | 70 +-- UTests/utestSphericalDirectPeriodic.cpp | 44 +- 6 files changed, 444 insertions(+), 394 deletions(-) diff --git a/Examples/RotationFMM.cpp b/Examples/RotationFMM.cpp index 86be53a8..25b5d827 100755 --- a/Examples/RotationFMM.cpp +++ b/Examples/RotationFMM.cpp @@ -23,19 +23,19 @@ #include #include -#include "../../Src/Files/FFmaGenericLoader.hpp" +#include "Utils/FParameters.hpp" +#include "Files/FFmaGenericLoader.hpp" -#include "../../Src/Kernels/Rotation/FRotationKernel.hpp" -#include "../../Src/Kernels/Rotation/FRotationCell.hpp" +#include "Kernels/Rotation/FRotationKernel.hpp" +#include "Kernels/Rotation/FRotationCell.hpp" -#include "../../Src/Components/FSimpleLeaf.hpp" -#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp" +#include "Components/FSimpleLeaf.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 diff --git a/Src/Utils/FPoint.hpp b/Src/Utils/FPoint.hpp index fe84ba8b..e8f4b6ae 100755 --- a/Src/Utils/FPoint.hpp +++ b/Src/Utils/FPoint.hpp @@ -22,307 +22,322 @@ #include #include +#include "ScalFmmConfig.h" #include "FMath.hpp" #include "FGlobal.hpp" /** -* @author Berenger Bramas (berenger.bramas@inria.fr) -* @class FPoint -* Please read the license -* -* This class is a 3D vector. It can be used as a position -* or as a 3d forces vector etc. -*/ + * @author Berenger Bramas (berenger.bramas@inria.fr) + * @class FPoint + * Please read the license + * + * This class is a 3D vector. It can be used as a position + * or as a 3d forces vector etc. + */ class FPoint { private: - FReal data[3]; //< all positions x y z + FReal data[3]; //< all positions x y z public: - /** Default constructor (sets position to 0/0/0) */ - FPoint(){ - data[0] = data[1] = data[2] = FReal(0); - } - - /** Constructor from an array */ - explicit FPoint(const FReal inPosition[3]){ - data[0] = inPosition[0]; - data[1] = inPosition[1]; - data[2] = inPosition[2]; - } - - /** Constructor from values */ - explicit FPoint(const FReal inX,const FReal inY,const FReal inZ){ - data[0] = inX; - data[1] = inY; - data[2] = inZ; - } - - /** - * Copdata[1] constructor - * @param other the source class to copy - */ - FPoint(const FPoint& other) { - data[0] = other.data[0]; - data[1] = other.data[1]; - data[2] = other.data[2]; - } - - /** - * Assignment operator - * @param other the source class to copy - */ - FPoint(const FPoint& other, const FReal addset) { - data[0] = other.data[0] + addset; - data[1] = other.data[1] + addset; - data[2] = other.data[2] + addset; - } - - /** - * Copdata[1] constructor - * @param other the source class to copy - * @return this a reference to the current class - */ - FPoint& operator=(const FPoint& other){ - this->data[0] = other.data[0]; - this->data[1] = other.data[1]; - this->data[2] = other.data[2]; - return *this; - } - - /** - * Position setter - * @param other the source class to copy - * @return this a reference to the current class - */ - void setPosition(const FReal inX,const FReal inY,const FReal inZ){ - this->data[0] = inX; - this->data[1] = inY; - this->data[2] = inZ; - } - - /** - * Get x - * @return this->data[0] - */ - FReal getX() const{ - return this->data[0]; - } - - /** - * Get y - * @return this->data[1] - */ - FReal getY() const{ - return this->data[1]; - } - - /** - * Get z - * @return this->data[2] - */ - FReal getZ() const{ - return this->data[2]; - } - - /** - * Set x - * @param the new x - */ - void setX(const FReal inX){ - this->data[0] = inX; - } - - /** - * Set y - * @param the new y - */ - void setY(const FReal inY){ - this->data[1] = inY; - } - - /** - * Set z - * @param the new z - */ - void setZ(const FReal inZ){ - this->data[2] = inZ; - } - - /** - * Add to the x-dimension the inX value - * @param inXthe increment in x - */ - void incX(const FReal inX){ - this->data[0] += inX; - } - - /** - * Add to the y-dimension the inY value - * @param in<<<<<data[1] += inY; - } - - /** - * Add to z-dimension the inZ value - * @param inZ the increment in z - */ - void incZ(const FReal inZ){ - this->data[2] += inZ; - } - /** - * Get a pointer on the coordinate of FPoint - * @return the data value array - */ - FReal * getDataValue(){ - return this->data ; - } - /** - *Compute the distance to the origin - * @return the norm of the Fpoint - */ - FReal norm() const { - return FMath::Sqrt(this->data[0]*this->data[0]+this->data[1]*this->data[1] - +this->data[2]*this->data[2]) ; - } - - /** - *Compute the distance to the origin - * @return the square norm of the Fpoint - */ - FReal norm2() const { - return (this->data[0]*this->data[0]+this->data[1]*this->data[1] - +this->data[2]*this->data[2]) ; - } - - /** - * Subtract to all dim the inValue - * @param inValue the value to substract - * @return the current object after being subtracted - */ - FPoint& operator-=(const FReal inValue){ - this->data[0] -= inValue; - this->data[1] -= inValue; - this->data[2] -= inValue; - return *this; - } - - /** - * Affect to all dim the inValue - * @param inValue the value to affect - * @return the current object after being affected - */ - FPoint& operator+=(const FReal inValue){ - this->data[0] += inValue; - this->data[1] += inValue; - this->data[2] += inValue; - return *this; - } - - /** - * Subtract to all dim the other position - * @param other the value to substract - * @return the current object after being subtracted - */ - FPoint& operator-=(const FPoint& other){ - this->data[0] -= other.data[0]; - this->data[1] -= other.data[1]; - this->data[2] -= other.data[2]; - return *this; - } - - /** - * Affect to all dim the other position - * @param other the value to afect - * @return the current object after being affected - */ - FPoint& operator+=(const FPoint& other){ - this->data[0] += other.data[0]; - this->data[1] += other.data[1]; - this->data[2] += other.data[2]; - return *this; - } - - /** - * Affect to all dim the other position - * @param other the value to afect - * @return the current object after being affected - */ - FPoint& operator*=(const FReal value){ - this->data[0] *= value; - this->data[1] *= value; - this->data[2] *= value; - return *this; - } - - /** - * Operator F3Position minus FReal - * This substract inValue to all dimensions of the inPosition - * @param inPosition the position to compute - * @param inValue the value to decrease/substract position - * @return the resulting position - */ - friend inline FPoint operator-(const FPoint& inPosition, const FReal inValue){ - return FPoint(inPosition, -inValue); - } - - /** - * Operator F3Position plus FReal - * This affect from inValue all dimensions of the inPosition - * @param inPosition the position to compute - * @param inValue the value to increase/affect position - * @return the resulting position - */ - friend inline FPoint operator+(const FPoint& inPosition, const FReal inValue){ - return FPoint(inPosition, inValue); - } - - /** - * Operator F3Position minus F3Position - * This substract one from anther - * @param inPosition the position to reduce - * @param inOther the position to decrease/substract inPosition - * @return the resulting position - */ - friend inline FPoint operator-(const FPoint& inPosition, const FPoint& inOther){ - return FPoint(inPosition.data[0] - inOther.data[0], inPosition.data[1] - inOther.data[1], inPosition.data[2] - inOther.data[2]); - } - - /** - * Operator F3Position plus F3Position - * This substract one from anther - * @param inPosition the position to reduce - * @param inOther the position to increase inPosition - * @return the resulting position - */ - friend inline FPoint operator+(const FPoint& inPosition, const FPoint& inOther){ - return FPoint(inPosition.data[0] + inOther.data[0], inPosition.data[1] + inOther.data[1], inPosition.data[2] + inOther.data[2]); - } - - /** - * Operator stream FPoint to std::ostream - * This can be used to simpldata[1] write out a position - * @param[in,out] output where to write the position - * @param[in] inPosition the position to write out - * @return the output for multiple << operators - */ - template - friend StreamClass& operator<<(StreamClass& output, const FPoint& inPosition){ - output << "(" << inPosition.getX() << ", " << inPosition.getY() << ", " << inPosition.getZ() <<")"; - return output; // for multiple << operators. - } - - /** Save current object */ - template - void save(BufferWriterClass& buffer) const { - buffer << data[0] << data[1] << data[2]; - } - /** Retrieve current object */ - template - void restore(BufferReaderClass& buffer) { - buffer >> data[0] >> data[1] >> data[2]; - } + /** Default constructor (sets position to 0/0/0) */ + FPoint(){ + data[0] = data[1] = data[2] = FReal(0); + } + + /** Constructor from an array */ + explicit FPoint(const FReal inPosition[3]){ + data[0] = inPosition[0]; + data[1] = inPosition[1]; + data[2] = inPosition[2]; + } + + /** Constructor from values */ + explicit FPoint(const FReal inX,const FReal inY,const FReal inZ){ + data[0] = inX; + data[1] = inY; + data[2] = inZ; + } + + /** + * Copdata[1] constructor + * @param other the source class to copy + */ + FPoint(const FPoint& other) { + data[0] = other.data[0]; + data[1] = other.data[1]; + data[2] = other.data[2]; + } + + /** + * Assignment operator + * @param other the source class to copy + */ + FPoint(const FPoint& other, const FReal addset) { + data[0] = other.data[0] + addset; + data[1] = other.data[1] + addset; + data[2] = other.data[2] + addset; + } + + /** + * Copdata[1] constructor + * @param other the source class to copy + * @return this a reference to the current class + */ + FPoint& operator=(const FPoint& other){ + this->data[0] = other.data[0]; + this->data[1] = other.data[1]; + this->data[2] = other.data[2]; + return *this; + } + + /** + * Position setter + * @param other the source class to copy + * @return this a reference to the current class + */ + void setPosition(const FReal inX,const FReal inY,const FReal inZ){ + this->data[0] = inX; + this->data[1] = inY; + this->data[2] = inZ; + } + + /** + * Get x + * @return this->data[0] + */ + FReal getX() const{ + return this->data[0]; + } + + /** + * Get y + * @return this->data[1] + */ + FReal getY() const{ + return this->data[1]; + } + + /** + * Get z + * @return this->data[2] + */ + FReal getZ() const{ + return this->data[2]; + } + + /** + * Set x + * @param the new x + */ + void setX(const FReal inX){ + this->data[0] = inX; + } + + /** + * Set y + * @param the new y + */ + void setY(const FReal inY){ + this->data[1] = inY; + } + + /** + * Set z + * @param the new z + */ + void setZ(const FReal inZ){ + this->data[2] = inZ; + } + + /** + * Add to the x-dimension the inX value + * @param inXthe increment in x + */ + void incX(const FReal inX){ + this->data[0] += inX; + } + + /** + * Add to the y-dimension the inY value + * @param in<<<<<data[1] += inY; + } + + /** + * Add to z-dimension the inZ value + * @param inZ the increment in z + */ + void incZ(const FReal inZ){ + this->data[2] += inZ; + } + /** + * Get a pointer on the coordinate of FPoint + * @return the data value array + */ + FReal * getDataValue(){ + return this->data ; + } + /** + *Compute the distance to the origin + * @return the norm of the Fpoint + */ + FReal norm() const { + return FMath::Sqrt(this->data[0]*this->data[0]+this->data[1]*this->data[1] + +this->data[2]*this->data[2]) ; + } + + /** + *Compute the distance to the origin + * @return the square norm of the Fpoint + */ + FReal norm2() const { + return (this->data[0]*this->data[0]+this->data[1]*this->data[1] + +this->data[2]*this->data[2]) ; + } + + /** + * Subtract to all dim the inValue + * @param inValue the value to substract + * @return the current object after being subtracted + */ + FPoint& operator-=(const FReal inValue){ + this->data[0] -= inValue; + this->data[1] -= inValue; + this->data[2] -= inValue; + return *this; + } + + /** + * Affect to all dim the inValue + * @param inValue the value to affect + * @return the current object after being affected + */ + FPoint& operator+=(const FReal inValue){ + this->data[0] += inValue; + this->data[1] += inValue; + this->data[2] += inValue; + return *this; + } + + /** + * Subtract to all dim the other position + * @param other the value to substract + * @return the current object after being subtracted + */ + FPoint& operator-=(const FPoint& other){ + this->data[0] -= other.data[0]; + this->data[1] -= other.data[1]; + this->data[2] -= other.data[2]; + return *this; + } + + /** + * Affect to all dim the other position + * @param other the value to afect + * @return the current object after being affected + */ + FPoint& operator+=(const FPoint& other){ + this->data[0] += other.data[0]; + this->data[1] += other.data[1]; + this->data[2] += other.data[2]; + return *this; + } + + /** + * Affect to all dim the other position + * @param other the value to afect + * @return the current object after being affected + */ + FPoint& operator*=(const FReal value){ + this->data[0] *= value; + this->data[1] *= value; + this->data[2] *= value; + return *this; + } + + /** + * Operator F3Position minus FReal + * This substract inValue to all dimensions of the inPosition + * @param inPosition the position to compute + * @param inValue the value to decrease/substract position + * @return the resulting position + */ + friend inline FPoint operator-(const FPoint& inPosition, const FReal inValue){ + return FPoint(inPosition, -inValue); + } + + /** + * Operator F3Position plus FReal + * This affect from inValue all dimensions of the inPosition + * @param inPosition the position to compute + * @param inValue the value to increase/affect position + * @return the resulting position + */ + friend inline FPoint operator+(const FPoint& inPosition, const FReal inValue){ + return FPoint(inPosition, inValue); + } + + /** + * Operator F3Position minus F3Position + * This substract one from anther + * @param inPosition the position to reduce + * @param inOther the position to decrease/substract inPosition + * @return the resulting position + */ + friend inline FPoint operator-(const FPoint& inPosition, const FPoint& inOther){ + return FPoint(inPosition.data[0] - inOther.data[0], inPosition.data[1] - inOther.data[1], inPosition.data[2] - inOther.data[2]); + } + + /** + * Operator F3Position plus F3Position + * This substract one from anther + * @param inPosition the position to reduce + * @param inOther the position to increase inPosition + * @return the resulting position + */ + friend inline FPoint operator+(const FPoint& inPosition, const FPoint& inOther){ + return FPoint(inPosition.data[0] + inOther.data[0], inPosition.data[1] + inOther.data[1], inPosition.data[2] + inOther.data[2]); + } + + /** + * Operator stream FPoint to std::ostream + * This can be used to simpldata[1] write out a position + * @param[in,out] output where to write the position + * @param[in] inPosition the position to write out + * @return the output for multiple << operators + */ + template + friend StreamClass& operator<<(StreamClass& output, const FPoint& inPosition){ + output << "(" << inPosition.getX() << ", " << inPosition.getY() << ", " << inPosition.getZ() <<")"; + return output; // for multiple << operators. + } + /** + * Operator stream FPoint to std::istream + * This can be used to simpldata[1] write out a position + * @param[in,out] input where to write the position + * @param[out] outPosition the position to write out + * @return the input for multiple << operators + */ + template + friend StreamClass& operator>>(StreamClass& input, FPoint& outPosition){ + FReal x,y,z; + input >> x>> y>> z ; + outPosition.setPosition(x,y,z); + return input; // for multiple << operators. + } + + /** Save current object */ + template + void save(BufferWriterClass& buffer) const { + buffer << data[0] << data[1] << data[2]; + } + /** Retrieve current object */ + template + void restore(BufferReaderClass& buffer) { + buffer >> data[0] >> data[1] >> data[2]; + } }; diff --git a/UTests/FUTester.hpp b/UTests/FUTester.hpp index c8c97299..11e023a2 100755 --- a/UTests/FUTester.hpp +++ b/UTests/FUTester.hpp @@ -1,5 +1,5 @@ // =================================================================================== -// 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 // This software is a computer program whose purpose is to compute the FMM. // @@ -39,7 +39,7 @@ int main(void){\ * Please refer to testUTest.cpp to see an example * @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 */ template @@ -61,8 +61,8 @@ class FUTester{ int currentTest; //< current processing test in the run int currentStep; //< current processing step in the run - int failledSteps; //< number of failled step in the current test - int failledTests; //< number of failled tests + int failledSteps; //< number of failed step in the current test + int failledTests; //< number of failed tests protected: /** Constructor */ @@ -94,7 +94,7 @@ protected: */ void AddTest(TestFunc inFunc){ char buff[256]; - sprintf(buff,"Unamed Test number %d",totalTests+1); + sprintf(buff,"Unnamed Test number %d",totalTests+1); AddTest(inFunc,buff); } @@ -119,11 +119,11 @@ protected: void Print(const Output& value){ std::cout<< "--- Output from program : " << value << "\n"; } - - /** + + /** * To test * @param result the test result - * if result is false test failled + * if result is false test failed */ void uassert(const bool result){ ++currentStep; diff --git a/UTests/utestChebyshevDirectPeriodic.cpp b/UTests/utestChebyshevDirectPeriodic.cpp index ddddfb8a..af79ba33 100755 --- a/UTests/utestChebyshevDirectPeriodic.cpp +++ b/UTests/utestChebyshevDirectPeriodic.cpp @@ -1,5 +1,5 @@ // =================================================================================== -// 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 // This software is a computer program whose purpose is to compute the FMM. // @@ -59,10 +59,10 @@ class TestChebyshevDirect : public FUTester { // Warning in make test the exec dir it Build/UTests // Load particles - const int NbLevels = 3; + const int NbLevels = 4; const int SizeSubLevels = 2; - const int PeriodicDeep = 1; - const int NbParticles = 6; + const int PeriodicDeep = 2; + const int NbParticles = 100; FRandomLoader loader(NbParticles); @@ -78,19 +78,23 @@ class TestChebyshevDirect : public FUTester { FReal physicalValue; FReal potential; }; + 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){ FPoint position; loader.fillParticle(&position); // put in tree - tree.insert(position, idxPart, 0.10); + value *= coeff ; + sum += value ; + // put in tree + tree.insert(position, idxPart, value); // get copy - particles[idxPart].position = position; - particles[idxPart].physicalValue = 0.10; - particles[idxPart].potential = 0.0; - particles[idxPart].forces[0] = 0.0; - particles[idxPart].forces[1] = 0.0; - particles[idxPart].forces[2] = 0.0; + particles[idxPart].position = position; + particles[idxPart].physicalValue = value; + particles[idxPart].potential = 0.0; + particles[idxPart].forces[0] = 0.0; + particles[idxPart].forces[1] = 0.0; + particles[idxPart].forces[2] = 0.0; } // Run FMM @@ -117,7 +121,8 @@ class TestChebyshevDirect : public FUTester { particles[idxOther].position.getZ(),particles[idxOther].physicalValue, &particles[idxOther].forces[0],&particles[idxOther].forces[1], &particles[idxOther].forces[2],&particles[idxOther].potential); - } + + } for(int idxX = min.getX() ; idxX <= max.getX() ; ++idxX){ for(int idxY = min.getY() ; idxY <= max.getY() ; ++idxY){ for(int idxZ = min.getZ() ; idxZ <= max.getZ() ; ++idxZ){ @@ -149,10 +154,13 @@ class TestChebyshevDirect : public FUTester { Print("Compute Diff..."); FMath::FAccurater potentialDiff; FMath::FAccurater fx, fy, fz; + FReal energy= 0.0 , energyD = 0.0 ; + { // Check that each particle has been summed with all other 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 forcesY = leaf->getTargets()->getForcesY(); const FReal*const forcesZ = leaf->getTargets()->getForcesZ(); @@ -165,6 +173,8 @@ class TestChebyshevDirect : public FUTester { fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]); 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; } }); } @@ -172,31 +182,40 @@ class TestChebyshevDirect : public FUTester { delete[] particles; // Print for information + Print("Potential diff is = "); - Print(potentialDiff.getL2Norm()); - Print(potentialDiff.getInfNorm()); + printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm()); + printf(" RMSError %e\n",potentialDiff.getRMSError()); Print("Fx diff is = "); - Print(fx.getL2Norm()); - Print(fx.getInfNorm()); + printf(" L2Norm %e\n",fx.getRelativeL2Norm()); + printf(" RMSError %e\n",fx.getRMSError()); + Print(fx.getRelativeL2Norm()); + Print(fx.getRelativeInfNorm()); Print("Fy diff is = "); - Print(fy.getL2Norm()); - Print(fy.getInfNorm()); + printf(" L2Norm %e\n",fy.getRelativeL2Norm()); + printf(" RMSError %e\n",fy.getRMSError()); Print("Fz diff is = "); - Print(fz.getL2Norm()); - Print(fz.getInfNorm()); - - // Assert + printf(" L2Norm %e\n",fz.getRelativeL2Norm()); + 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)); + 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 MaximumDiffForces = FReal(9e-3); - uassert(potentialDiff.getL2Norm() < MaximumDiffPotential); - uassert(potentialDiff.getInfNorm() < MaximumDiffPotential); - uassert(fx.getL2Norm() < MaximumDiffForces); - uassert(fx.getInfNorm() < MaximumDiffForces); - uassert(fy.getL2Norm() < MaximumDiffForces); - uassert(fy.getInfNorm() < MaximumDiffForces); - uassert(fz.getL2Norm() < MaximumDiffForces); - uassert(fz.getInfNorm() < MaximumDiffForces); + uassert(potentialDiff.getL2Norm() < MaximumDiffPotential); //1 + uassert(potentialDiff.getRMSError() < MaximumDiffPotential); //2 + uassert(fx.getL2Norm() < MaximumDiffForces); //3 + uassert(fx.getRMSError() < MaximumDiffForces); //4 + uassert(fy.getL2Norm() < MaximumDiffForces); //5 + uassert(fy.getRMSError() < MaximumDiffForces); //6 + uassert(fz.getL2Norm() < MaximumDiffForces); //8 + uassert(fz.getRMSError() < MaximumDiffForces); //8 } @@ -220,8 +239,8 @@ class TestChebyshevDirect : public FUTester { /** TestChebKernel */ void TestChebKernel(){ - const unsigned int ORDER = 5; - const FReal epsilon = FReal(1e-5); + const unsigned int ORDER = 6; + const FReal epsilon = FReal(1e-6); typedef FP2PParticleContainerIndexed<> ContainerClass; typedef FSimpleLeaf LeafClass; typedef FInterpMatrixKernelR MatrixKernelClass; @@ -230,6 +249,7 @@ class TestChebyshevDirect : public FUTester { typedef FChebKernel KernelClass; typedef FFmmAlgorithmPeriodic FmmClass; // run test + std::cout <<" TEST 1 "<(epsilon); } @@ -245,6 +265,7 @@ class TestChebyshevDirect : public FUTester { typedef FChebSymKernel KernelClass; typedef FFmmAlgorithmPeriodic FmmClass; // run test + std::cout <(epsilon); } @@ -256,8 +277,8 @@ class TestChebyshevDirect : public FUTester { /** set test */ void SetTests(){ - AddTest(&TestChebyshevDirect::TestChebKernel,"Test Chebyshev Kernel with one big SVD"); - AddTest(&TestChebyshevDirect::TestChebSymKernel,"Test Chebyshev Kernel with 16 small SVDs and symmetries"); + AddTest(&TestChebyshevDirect::TestChebKernel,"Test Chebyshev Kernel with one big SVD") ; + AddTest(&TestChebyshevDirect::TestChebSymKernel,"Test Chebyshev Kernel with 16 small SVDs and symmetries"); } }; diff --git a/UTests/utestRotationDirectPeriodic.cpp b/UTests/utestRotationDirectPeriodic.cpp index 1a65a53a..a3767ab5 100755 --- a/UTests/utestRotationDirectPeriodic.cpp +++ b/UTests/utestRotationDirectPeriodic.cpp @@ -33,7 +33,7 @@ class TestRotationDirectPeriodic : public FUTester { /** Here we test only the P2P */ void TestPeriodicFmm(){ - static const int P = 9; + static const int P = 18; typedef FRotationCell

CellClass; typedef FP2PParticleContainerIndexed<> ContainerClass; @@ -45,9 +45,9 @@ class TestRotationDirectPeriodic : public FUTester { typedef FFmmAlgorithmPeriodic FmmClass; // Parameters - const int NbLevels = 4; + const int NbLevels = 4; const int SizeSubLevels = 2; - const int PeriodicDeep = -1; + const int PeriodicDeep = 2; const int NbParticles = 100; FRandomLoader loader(NbParticles); @@ -58,26 +58,32 @@ class TestRotationDirectPeriodic : public FUTester { FReal physicalValue; FReal potential; }; + 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){ FPoint position; loader.fillParticle(&position); + value *= coeff ; + sum += value ; // put in tree - tree.insert(position, idxPart, 0.10); + tree.insert(position, idxPart, value); // get copy - particles[idxPart].position = position; - particles[idxPart].physicalValue = 0.10; + particles[idxPart].position = position; + particles[idxPart].physicalValue = value; particles[idxPart].potential = 0.0; particles[idxPart].forces[0] = 0.0; particles[idxPart].forces[1] = 0.0; particles[idxPart].forces[2] = 0.0; } - + if (FMath::Abs(sum)> 0.00001){ + std::cerr << "Sum of charges is not equal zero!!! (sum=<<"< { Print("Compute Diff..."); FMath::FAccurater potentialDiff; 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){ 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 forcesZ = leaf->getTargets()->getForcesZ(); const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles(); @@ -143,38 +151,40 @@ class TestRotationDirectPeriodic : public FUTester { fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]); fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]); fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]); - } + energy += potentials[idxPart]*physicalValues[idxPart]; + } }); } Print("Potential diff is = "); - printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm()); - printf(" RMSError %e\n",potentialDiff.getRMSError()); + printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm()); + printf(" RMSError %e\n",potentialDiff.getRMSError()); Print("Fx diff is = "); - printf(" L2Norm %e\n",fx.getRelativeL2Norm()); - printf(" RMSError %e\n",fx.getRMSError()); + printf(" L2Norm %e\n",fx.getRelativeL2Norm()); + printf(" RMSError %e\n",fx.getRMSError()); Print(fx.getRelativeL2Norm()); Print(fx.getRelativeInfNorm()); Print("Fy diff is = "); - printf(" L2Norm %e\n",fy.getRelativeL2Norm()); - printf(" RMSError %e\n",fy.getRMSError()); + printf(" L2Norm %e\n",fy.getRelativeL2Norm()); + printf(" RMSError %e\n",fy.getRMSError()); Print("Fz diff is = "); - printf(" L2Norm %e\n",fz.getRelativeL2Norm()); - printf(" RMSError %e\n",fz.getRMSError()); + printf(" L2Norm %e\n",fz.getRelativeL2Norm()); + printf(" RMSError %e\n",fz.getRMSError()); 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.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 + uassert(potentialDiff.getL2Norm() < MaximumDiffPotential); //1 + uassert(potentialDiff.getRMSError() < MaximumDiffPotential); //2 + uassert(fx.getL2Norm() < MaximumDiffForces); //3 + uassert(fx.getRMSError() < MaximumDiffForces); //4 + uassert(fy.getL2Norm() < MaximumDiffForces); //5 + uassert(fy.getRMSError() < MaximumDiffForces); //6 + uassert(fz.getL2Norm() < MaximumDiffForces); //8 + uassert(fz.getRMSError() < MaximumDiffForces); //8 delete[] particles; } diff --git a/UTests/utestSphericalDirectPeriodic.cpp b/UTests/utestSphericalDirectPeriodic.cpp index d0557d88..1c82d776 100755 --- a/UTests/utestSphericalDirectPeriodic.cpp +++ b/UTests/utestSphericalDirectPeriodic.cpp @@ -46,11 +46,11 @@ class TestSphericalDirectPeriodic : public FUTester typedef FFmmAlgorithmPeriodic FmmClass; // Parameters - const int NbLevels = 3; + const int NbLevels = 4; const int SizeSubLevels = 2; - const int PeriodicDeep = 3; - const int DevP = 9; - const int NbParticles = 1; + const int PeriodicDeep = 1; + const int DevP = 14; + const int NbParticles = 100; FSphericalCell::Init(DevP); @@ -62,15 +62,18 @@ class TestSphericalDirectPeriodic : public FUTester FReal physicalValue; 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){ FPoint position; loader.fillParticle(&position); + value *= coeff ; + sum += value ; // put in tree - tree.insert(position, idxPart, 0.10); + tree.insert(position, idxPart, value); // get copy - particles[idxPart].position = position; - particles[idxPart].physicalValue = 0.10; + particles[idxPart].position = position; + particles[idxPart].physicalValue = value; particles[idxPart].potential = 0.0; particles[idxPart].forces[0] = 0.0; particles[idxPart].forces[1] = 0.0; @@ -131,7 +134,7 @@ class TestSphericalDirectPeriodic : public FUTester Print("Compute Diff..."); FMath::FAccurater potentialDiff; 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 tree.forEachLeaf([&](LeafClass* leaf){ @@ -149,29 +152,30 @@ class TestSphericalDirectPeriodic : public FUTester fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]); fy.add(particles[indexPartOrig].forces[1],forcesY[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 = "); - printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm()); - printf(" RMSError %e\n",potentialDiff.getRMSError()); + printf(" L2Norm %e\n",potentialDiff.getRelativeL2Norm()); + printf(" RMSError %e\n",potentialDiff.getRMSError()); Print("Fx diff is = "); - printf(" L2Norm %e\n",fx.getRelativeL2Norm()); - printf(" RMSError %e\n",fx.getRMSError()); + printf(" L2Norm %e\n",fx.getRelativeL2Norm()); + printf(" RMSError %e\n",fx.getRMSError()); Print(fx.getRelativeL2Norm()); Print(fx.getRelativeInfNorm()); Print("Fy diff is = "); - printf(" L2Norm %e\n",fy.getRelativeL2Norm()); - printf(" RMSError %e\n",fy.getRMSError()); + printf(" L2Norm %e\n",fy.getRelativeL2Norm()); + printf(" RMSError %e\n",fy.getRMSError()); Print("Fz diff is = "); - printf(" L2Norm %e\n",fz.getRelativeL2Norm()); - printf(" RMSError %e\n",fz.getRMSError()); + printf(" L2Norm %e\n",fz.getRelativeL2Norm()); + printf(" RMSError %e\n",fz.getRMSError()); FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm() + fz.getRelativeL2Norm() *fz.getRelativeL2Norm() ); - printf(" L2 Force Error= %e\n",FMath::Sqrt(L2error)) ; - printf(" Energy = %.12e\n",energy); + 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 -- 2.24.1