Skip to content

S
ScalFMM
Commit 935a7bca authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files
Options

Improvments

parent e226f0b7
Hide whitespace changes
Inline Side-by-side
Showing with 77 additions and 66 deletions
UTests/utestChebyshevDirect.cpp
View file @ 935a7bca
......@@ -93,22 +93,25 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
// put in tree
tree.insert(position, idxPart, physicalValue);
// get copy
particles[idxPart].position = position;
particles[idxPart].position = position;
particles[idxPart].physicalValue = physicalValue;
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
/////////////////////////////////////////////////////////////////////////////////////////////////
// Run FMM computation
/////////////////////////////////////////////////////////////////////////////////////////////////
Print("Fmm...");
KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());
FmmClass algo(&tree,&kernels);
algo.execute();
//0
FReal energy= 0.0 , energyD = 0.0 ;
/////////////////////////////////////////////////////////////////////////////////////////////////
// Run direct computation
/////////////////////////////////////////////////////////////////////////////////////////////////
const MatrixKernelClass MatrixKernel;
Print("Direct...");
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
......@@ -123,12 +126,15 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
&particles[idxOther].forces[2],&particles[idxOther].potential);
}
}
for(int idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){
energyD += particles[idx].potential*particles[idx].physicalValue ;
}
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compare
/////////////////////////////////////////////////////////////////////////////////////////////////
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){
......@@ -147,8 +153,6 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
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;
}
});
}
......@@ -156,31 +160,32 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
delete[] particles;
// Print for information
// Print for information
Print("Potential diff is = ");
printf(" Pot L2Norm %e\n",potentialDiff.getRelativeL2Norm());
printf(" Pot RMSError %e\n",potentialDiff.getRMSError());
printf(" Pot L2Norm %e\n",potentialDiff.getL2Norm());
printf(" Pot RL2Norm %e\n",potentialDiff.getRelativeL2Norm());
printf(" Pot RMSError %e\n",potentialDiff.getRMSError());
Print("Fx diff is = ");
printf(" Fx L2Norm %e\n",fx.getRelativeL2Norm());
printf(" Fx RMSError %e\n",fx.getRMSError());
Print(fx.getRelativeL2Norm());
Print(fx.getRelativeInfNorm());
printf(" Fx L2Norm %e\n",fx.getL2Norm());
printf(" Fx RL2Norm %e\n",fx.getRelativeL2Norm());
printf(" Fx RMSError %e\n",fx.getRMSError());
Print("Fy diff is = ");
printf(" Fy L2Norm %e\n",fy.getRelativeL2Norm());
printf(" Fy RMSError %e\n",fy.getRMSError());
printf(" Fy L2Norm %e\n",fy.getL2Norm());
printf(" Fy RL2Norm %e\n",fy.getRelativeL2Norm());
printf(" Fy RMSError %e\n",fy.getRMSError());
Print("Fz diff is = ");
printf(" Fz L2Norm %e\n",fz.getRelativeL2Norm());
printf(" Fz RMSError %e\n",fz.getRMSError());
printf(" Fz L2Norm %e\n",fz.getL2Norm());
printf(" Fz RL2Norm %e\n",fz.getRelativeL2Norm());
printf(" Fz 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 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-5);
const FReal MaximumDiffForces = FReal(9e-3);
const FReal MaximumDiffPotential = FReal(9e-3);
const FReal MaximumDiffForces = FReal(9e-2);
uassert(potentialDiff.getL2Norm() < MaximumDiffPotential); //1
......@@ -192,7 +197,7 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
uassert(fz.getL2Norm() < MaximumDiffForces); //8
uassert(fz.getRMSError() < MaximumDiffForces); //8
uassert(L2error < MaximumDiffForces); //9 Total Force
uassert(FMath::Abs(energy-energyD) < MaximumDiffPotential); //10 Total Energy
uassert(FMath::Abs(energy-energyD) < 10*MaximumDiffPotential); //10 Total Energy
}
......
UTests/utestChebyshevDirectPeriodic.cpp
View file @ 935a7bca
......@@ -62,7 +62,7 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
const int NbLevels = 4;
const int SizeSubLevels = 2;
const int PeriodicDeep = 2;
const int NbParticles = 100;
const int NbParticles = 500;
FRandomLoader loader(NbParticles);
......@@ -96,19 +96,22 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
particles[idxPart].forces[1] = 0.0;
particles[idxPart].forces[2] = 0.0;
}
// Run FMM
/////////////////////////////////////////////////////////////////////////////////////////////////
// Run FMM computation
/////////////////////////////////////////////////////////////////////////////////////////////////
Print("Fmm...");
FmmClass algo(&tree,PeriodicDeep );
KernelClass kernels(algo.extendedTreeHeight(), algo.extendedBoxWidth(), algo.extendedBoxCenter());
algo.setKernel(&kernels);
algo.execute();
/////////////////////////////////////////////////////////////////////////////////////////////////
// Run direct computation
/////////////////////////////////////////////////////////////////////////////////////////////////
Print("Direct...");
FTreeCoordinate min, max;
algo.repetitionsIntervals(&min, &max);
FReal energy= 0.0 , energyD = 0.0 ;
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
......@@ -148,12 +151,16 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
}
}
}
for(int idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){
energyD += particles[idx].potential*particles[idx].physicalValue ;
}
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compare
/////////////////////////////////////////////////////////////////////////////////////////////////
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
......@@ -172,8 +179,7 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
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;
energy+=potentials[idxPart]*physicalValues[idxPart];
}
});
}
......@@ -182,43 +188,43 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
// Print for information
Print("Potential diff is = ");
printf(" Pot L2Norm %e\n",potentialDiff.getRelativeL2Norm());
printf(" Pot RMSError %e\n",potentialDiff.getRMSError());
Print("Fx diff is = ");
printf(" Fx L2Norm %e\n",fx.getRelativeL2Norm());
printf(" Fx RMSError %e\n",fx.getRMSError());
Print(fx.getRelativeL2Norm());
Print(fx.getRelativeInfNorm());
Print("Fy diff is = ");
printf(" Fy L2Norm %e\n",fy.getRelativeL2Norm());
printf(" Fy RMSError %e\n",fy.getRMSError());
Print("Fz diff is = ");
printf(" Fz L2Norm %e\n",fz.getRelativeL2Norm());
printf(" Fz RMSError %e\n",fz.getRMSError());
FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm() + fz.getRelativeL2Norm() *fz.getRelativeL2Norm() );
Print("Potential diff is = ");
printf(" Pot L2Norm %e\n",potentialDiff.getL2Norm());
printf(" Pot RL2Norm %e\n",potentialDiff.getRelativeL2Norm());
printf(" Pot RMSError %e\n",potentialDiff.getRMSError());
Print("Fx diff is = ");
printf(" Fx L2Norm %e\n",fx.getL2Norm());
printf(" Fx RL2Norm %e\n",fx.getRelativeL2Norm());
printf(" Fx RMSError %e\n",fx.getRMSError());
Print("Fy diff is = ");
printf(" Fy L2Norm %e\n",fy.getL2Norm());
printf(" Fy RL2Norm %e\n",fy.getRelativeL2Norm());
printf(" Fy RMSError %e\n",fy.getRMSError());
Print("Fz diff is = ");
printf(" Fz L2Norm %e\n",fz.getL2Norm());
printf(" Fz RL2Norm %e\n",fz.getRelativeL2Norm());
printf(" Fz 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 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); //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
uassert(fz.getRMSError() < MaximumDiffForces); //8
uassert(L2error < MaximumDiffForces); //9 Total Force
uassert(FMath::Abs(energy-energyD) < MaximumDiffPotential); //10 Total Energy
const FReal MaximumDiffPotential = FReal(9e-3);
const FReal MaximumDiffForces = FReal(9e-2);
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
uassert(L2error < MaximumDiffForces); //9 Total Force
uassert(FMath::Abs(energy-energyD) < 10*MaximumDiffPotential); //10 Total Energy
}
......
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