Commit bd6e637b authored by COULAUD Olivier's avatar COULAUD Olivier

Add RMS test rather than Infnorm

parent fd6ce2c6
......@@ -281,5 +281,6 @@ INCLUDE(CPack)
#
# OUTPUT
#
MESSAGE(STATUS "CMAKE_CXX_FLAGS = ${CMAKE_CXX_FLAGS}")
MESSAGE(STATUS "ScaLFMM_CXX_FLAGS = ${ScaLFMM_CXX_FLAGS}")
MESSAGE(STATUS "SCALFMM_LIBRARIES = ${SCALFMM_LIBRARIES}")
......@@ -4,4 +4,7 @@
* \section Periodic-1 Periodicity in usual FMM
* \subsection limitation
* support only cubic box
*/
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
// This software is a computer program whose purpose is to compute the FMM.
//
......@@ -203,36 +203,43 @@ struct FMath{
/** A class to compute accuracy */
class FAccurater {
int nbElements;
FReal l2Dot;
FReal l2Diff;
FReal max;
FReal maxDiff;
public:
FAccurater() : l2Dot(0.0), l2Diff(0.0), max(0.0), maxDiff(0.0) {
FAccurater() : nbElements(0),l2Dot(0.0), l2Diff(0.0), max(0.0), maxDiff(0.0) {
}
/** with inital values */
FAccurater(const FReal inGood[], const FReal inBad[], const int nbValues)
: l2Dot(0.0), l2Diff(0.0), max(0.0), maxDiff(0.0) {
: nbElements(0),l2Dot(0.0), l2Diff(0.0), max(0.0), maxDiff(0.0) {
add(inGood, inBad, nbValues);
}
/** Add value to the current list */
void add(const FReal inGood, const FReal inBad){
l2Diff += (inBad - inGood) * (inBad - inGood);
l2Dot += inGood * inGood;
max = Max(max , Abs(inGood));
maxDiff = Max(maxDiff, Abs(inGood-inBad));
l2Diff += (inBad - inGood) * (inBad - inGood);
l2Dot += inGood * inGood;
max = Max(max , Abs(inGood));
maxDiff = Max(maxDiff, Abs(inGood-inBad));
nbElements += 1 ;
}
/** Add array of values */
void add(const FReal inGood[], const FReal inBad[], const int nbValues){
for(int idx = 0 ; idx < nbValues ; ++idx){
add(inGood[idx],inBad[idx]);
}
nbElements += nbValues ;
}
/** Get the L2 norm */
/** Get the root mean squared error*/
FReal getL2Norm() const{
return Sqrt(l2Diff );
}
/** Get the L2 norm */
FReal getRMSError() const{
return Sqrt(l2Diff /static_cast<FReal>(nbElements));
}
/** Get the inf norm */
FReal getInfNorm() const{
return maxDiff;
......@@ -254,10 +261,11 @@ struct FMath{
void reset()
{
l2Dot = FReal(0);
l2Diff = FReal(0);;
max = FReal(0);
maxDiff = FReal(0);
l2Dot = FReal(0);
l2Diff = FReal(0);;
max = FReal(0);
maxDiff = FReal(0);
nbElements = 0 ;
}
};
};
......
......@@ -71,7 +71,7 @@ int main(int argc, char ** argv){
#ifdef ScalFMM_USE_BLAS
// begin Chebyshev kernel
// accuracy
const unsigned int ORDER = 12;
const unsigned int ORDER = 13;
// typedefs
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
......@@ -463,13 +463,16 @@ int main(int argc, char ** argv){
if(FParameters::existParameter(argc, argv, "-verbose")){
std::cout << ">> index " << particles[indexPartOrig].index << std::endl;
std::cout << "Good x " << particles[indexPartOrig].position.getX() << " y " << particles[indexPartOrig].position.getY() << " z " << particles[indexPartOrig].position.getZ() << std::endl;
std::cout << std::fixed << std::setprecision(5) ;
std::cout << "FMM x " << positionsX[idxPart] << " y " << positionsY[idxPart] << " z " << positionsZ[idxPart] << std::endl;
std::cout << "Good fx " <<particles[indexPartOrig].forces[0] << " fy " << particles[indexPartOrig].forces[1] << " fz " << particles[indexPartOrig].forces[2] << std::endl;
std::cout << "FMM fx " << forcesX[idxPart] << " fy " <<forcesY[idxPart] << " fz " << forcesZ[idxPart] << std::endl;
std::cout << "ratio fx " << particles[indexPartOrig].forces[0]/forcesX[idxPart] << " fy " <<particles[indexPartOrig].forces[1]/forcesY[idxPart] << " fz " << particles[indexPartOrig].forces[2]/forcesZ[idxPart] << std::endl;
std::cout << std::scientific << std::setprecision(5) ;
std::cout << "Diff fx " << particles[indexPartOrig].forces[0]-forcesX[idxPart] << " fy " <<particles[indexPartOrig].forces[1]-forcesY[idxPart] << " fz " << particles[indexPartOrig].forces[2]-forcesZ[idxPart] << std::endl;
// std::cout << "GOOD physical value " << particles[indexPartOrig].physicalValue << " potential " << particles[indexPartOrig].potential << std::endl;
// std::cout << "FMM physical value " << physicalValues[idxPart] << " potential " << potentials[idxPart] << " energy cumul " << energy<<std::endl;
std::cout << std::fixed << std::setprecision(5) ;
std::cout << "\n";
}
if(FParameters::existParameter(argc, argv, "-saveError")){
......@@ -531,6 +534,9 @@ int main(int argc, char ** argv){
printf("FmmEwald L2Norm %e\n",fz.getRelativeL2Norm());
printf("FmmEwald InfNorm %e\n",fz.getRelativeInfNorm());
//
double L2error = (fx.getL2Norm()*fx.getL2Norm() + fy.getL2Norm()*fy.getL2Norm() + fz.getL2Norm() *fz.getL2Norm() );
printf("FmmEwald RMS Force Error= %e\n",FMath::Sqrt(L2error/loader->getNumberOfParticles())) ;
//
printf("FmmEwald Energy FMM= %.12e\n",energy);
printf("FmmEwald Energy EWALD= %.12e\n",loader->getEnergy());
printf("FmmEwald Energy EWALD - Energy FMM = %e\n",loader->getEnergy()-energy);
......@@ -543,42 +549,8 @@ int main(int argc, char ** argv){
}
//
// -----------------------------------------------------
// ReGenerate file
// ----------------------------------------------------------------------------
//
if( FParameters::existParameter(argc, argv, "-gen") ){
std::cout << "Generate ewal.out from input file" << std::endl;
std::ofstream fileout("ewal.out",std::ifstream::out);
std::ifstream file(filename,std::ifstream::in);
if(file.is_open()){
const int bufferSize = 512;
char buffer[bufferSize];
file.getline(buffer, bufferSize);
fileout << buffer << '\n';
file.getline(buffer, bufferSize);
fileout << buffer << '\n';
if( !FParameters::existParameter(argc, argv, "-noper") ){
file.getline(buffer, bufferSize);
fileout << buffer << '\n';
file.getline(buffer, bufferSize);
fileout << buffer << '\n';
file.getline(buffer, bufferSize);
fileout << buffer << '\n';
}
for(int idxPart = 0 ; idxPart < loader->getNumberOfParticles() ; ++idxPart){
file.getline(buffer, bufferSize);
fileout << buffer << '\n';
file.getline(buffer, bufferSize);
fileout << buffer << '\n';
file.getline(buffer, bufferSize);
file.getline(buffer, bufferSize);
}
}
}
// end generate
// -----------------------------------------------------
delete[] particles;
delete[] particlesDirect;
......
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, B��renger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
......@@ -148,27 +148,33 @@ class TestRotationDirectPeriodic : public FUTester<TestRotationDirectPeriodic> {
}
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());
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)) ;
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);
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
delete[] particles;
}
......
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