Commit fc7361fa authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files
parents 950dd3ed c8823edd
......@@ -18,6 +18,7 @@ project(scalfmm)
# Active language
# -----------------------
ENABLE_LANGUAGE(CXX )
MESSAGE(STATUS " CXX ${CMAKE_CXX_COMPILER_ID}" )
# Options
OPTION( SCALFMM_USE_BLAS "Set to ON to build ScaFMM with BLAS" ON )
......@@ -43,8 +44,12 @@ if( SCALFMM_BUILD_DEBUG )
ADD_DEFINITIONS(-O0)
else()
SET(CMAKE_BUILD_TYPE Release)
IF(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
ADD_DEFINITIONS(-ip)
ENDIF()
endif()
# starpu
if( SCALFMM_USE_STARPU )
# starpu
......
......@@ -106,14 +106,17 @@ public:
* @param[out] Interpolator
*/
void assembleInterpolator(const unsigned int NumberOfLocalPoints,
const FPoint *const LocalPoints,
FReal *const Interpolator) const
const FPoint *const LocalPoints,
FReal *const Interpolator) const
{
// values of chebyshev polynomials of source particle: T_o(x_i)
FReal T_of_x[ORDER][3];
// loop: local points (mapped in [-1,1])
for (unsigned int m=0; m<NumberOfLocalPoints; ++m) {
FReal c0, c1, c2 ;
c0 = FReal(0.0) ;
c1 = FReal(1.) / ORDER;
c2 = FReal(2.) *c1 ;
// loop: local points (mapped in [-1,1])
for (unsigned int m=0; m<NumberOfLocalPoints; ++m) {
// evaluate chebyshev polynomials at local points
for (unsigned int o=1; o<ORDER; ++o) {
T_of_x[o][0] = BasisType::T(o, LocalPoints[m].getX());
......@@ -126,9 +129,14 @@ public:
Interpolator[n*nnodes + m] = FReal(1.);
for (unsigned int d=0; d<3; ++d) {
const unsigned int j = node_ids[n][d];
FReal S_d = FReal(1.) / ORDER;
// FReal S_d = FReal(1.) / ORDER;
// for (unsigned int o=1; o<ORDER; ++o)
// S_d += FReal(2.) / ORDER * T_of_x[o][d] * T_of_roots[o][j];
// Interpolator[n*nnodes + m] *= S_d;
FReal S_d = c0 ;
for (unsigned int o=1; o<ORDER; ++o)
S_d += FReal(2.) / ORDER * T_of_x[o][d] * T_of_roots[o][j];
S_d += T_of_x[o][d] * T_of_roots[o][j];
S_d = c1 + c2*S_d ;
Interpolator[n*nnodes + m] *= S_d;
}
}
......@@ -224,8 +232,10 @@ inline void FChebInterpolator<ORDER>::applyP2M(const FPoint& center,
const map_glob_loc map(center, width);
FPoint localPosition;
FReal T_of_x[ORDER][3];
FReal S[3];
FReal S[3],c1;
//
FReal xpx,ypy,zpz ;
c1 = FReal(8.) / nnodes ;
// loop over source particles
typename ContainerClass::ConstBasicIterator iter(*sourceParticles);
while(iter.hasNotFinished()){
......@@ -237,10 +247,17 @@ inline void FChebInterpolator<ORDER>::applyP2M(const FPoint& center,
T_of_x[0][0] = FReal(1.); T_of_x[1][0] = localPosition.getX();
T_of_x[0][1] = FReal(1.); T_of_x[1][1] = localPosition.getY();
T_of_x[0][2] = FReal(1.); T_of_x[1][2] = localPosition.getZ();
xpx = FReal(2.) * localPosition.getX() ;
ypy = FReal(2.) * localPosition.getY() ;
zpz = FReal(2.) * localPosition.getZ() ;
for (unsigned int o=2; o<ORDER; ++o) {
T_of_x[o][0] = FReal(2.) * localPosition.getX() * T_of_x[o-1][0] - T_of_x[o-2][0];
T_of_x[o][1] = FReal(2.) * localPosition.getY() * T_of_x[o-1][1] - T_of_x[o-2][1];
T_of_x[o][2] = FReal(2.) * localPosition.getZ() * T_of_x[o-1][2] - T_of_x[o-2][2];
T_of_x[o][0] = xpx * T_of_x[o-1][0] - T_of_x[o-2][0];
T_of_x[o][1] = ypy * T_of_x[o-1][1] - T_of_x[o-2][1];
T_of_x[o][2] = zpz * T_of_x[o-1][2] - T_of_x[o-2][2];
// T_of_x[o][0] = FReal(2.) * localPosition.getX() * T_of_x[o-1][0] - T_of_x[o-2][0];
// T_of_x[o][1] = FReal(2.) * localPosition.getY() * T_of_x[o-1][1] - T_of_x[o-2][1];
// T_of_x[o][2] = FReal(2.) * localPosition.getZ() * T_of_x[o-1][2] - T_of_x[o-2][2];
}
// anterpolate
......@@ -256,10 +273,15 @@ inline void FChebInterpolator<ORDER>::applyP2M(const FPoint& center,
S[2] += T_of_x[o][2] * T_of_roots[o][j[2]];
}
// gather contributions
S[0] *= FReal(2.); S[0] += FReal(1.);
S[1] *= FReal(2.); S[1] += FReal(1.);
S[2] *= FReal(2.); S[2] += FReal(1.);
multipoleExpansion[n] += FReal(1.) / nnodes * S[0] * S[1] * S[2] * sourceValue;
// Here we consider 1/2 S rather S then we multiply by 8 the results
// S[0] *= FReal(2.); S[0] += FReal(1.);
// S[1] *= FReal(2.); S[1] += FReal(1.);
// S[2] *= FReal(2.); S[2] += FReal(1.);
// multipoleExpansion[n] += c1 * S[0] * S[1] * S[2] * sourceValue;
S[0] += FReal(0.5);
S[1] += FReal(0.5);
S[2] += FReal(0.5);
multipoleExpansion[n] += c1 * S[0] * S[1] * S[2] * sourceValue;
}
// increment source iterator
iter.gotoNext();
......@@ -281,8 +303,10 @@ inline void FChebInterpolator<ORDER>::applyL2P(const FPoint& center,
const map_glob_loc map(center, width);
FPoint localPosition;
FReal T_of_x[ORDER][3];
FReal S[3];
FReal xpx,ypy,zpz ;
FReal S[3],c1;
//
c1 = FReal(8.) / nnodes ;
typename ContainerClass::BasicIterator iter(*localParticles);
while(iter.hasNotFinished()){
......@@ -293,6 +317,9 @@ inline void FChebInterpolator<ORDER>::applyL2P(const FPoint& center,
T_of_x[0][0] = FReal(1.); T_of_x[1][0] = localPosition.getX();
T_of_x[0][1] = FReal(1.); T_of_x[1][1] = localPosition.getY();
T_of_x[0][2] = FReal(1.); T_of_x[1][2] = localPosition.getZ();
xpx = FReal(2.) * localPosition.getX() ;
ypy = FReal(2.) * localPosition.getY() ;
zpz = FReal(2.) * localPosition.getZ() ;
for (unsigned int o=2; o<ORDER; ++o) {
T_of_x[o][0] = FReal(2.) * localPosition.getX() * T_of_x[o-1][0] - T_of_x[o-2][0];
T_of_x[o][1] = FReal(2.) * localPosition.getY() * T_of_x[o-1][1] - T_of_x[o-2][1];
......@@ -312,13 +339,19 @@ inline void FChebInterpolator<ORDER>::applyL2P(const FPoint& center,
S[2] += T_of_x[o][2] * T_of_roots[o][j[2]];
}
// gather contributions
S[0] *= FReal(2.); S[0] += FReal(1.);
S[1] *= FReal(2.); S[1] += FReal(1.);
S[2] *= FReal(2.); S[2] += FReal(1.);
// S[0] *= FReal(2.); S[0] += FReal(1.);
// S[1] *= FReal(2.); S[1] += FReal(1.);
// S[2] *= FReal(2.); S[2] += FReal(1.);
// targetValue += S[0] * S[1] * S[2] * localExpansion[n];
S[0] += FReal(0.5);
S[1] += FReal(0.5);
S[2] += FReal(0.5);
//
targetValue += S[0] * S[1] * S[2] * localExpansion[n];
}
// scale
targetValue /= nnodes;
// targetValue /= nnodes;
targetValue *= c1;
// set potential
iter.data().setPotential(targetValue);
......@@ -472,7 +505,10 @@ inline void FChebInterpolator<ORDER>::applyL2PTotal(const FPoint& center,
FReal T_of_x[ORDER][3];
FReal U_of_x[ORDER][3];
FReal P[3];
//
FReal xpx,ypy,zpz ;
FReal c1 = FReal(8.0) / nnodes ;
//
typename ContainerClass::BasicIterator iter(*localParticles);
while(iter.hasNotFinished()){
......@@ -481,22 +517,36 @@ inline void FChebInterpolator<ORDER>::applyL2PTotal(const FPoint& center,
// evaluate chebyshev polynomials of source particle
// T_0(x_i) and T_1(x_i)
xpx = FReal(2.) * localPosition.getX() ;
ypy = FReal(2.) * localPosition.getY() ;
zpz = FReal(2.) * localPosition.getZ() ;
//
T_of_x[0][0] = FReal(1.); T_of_x[1][0] = localPosition.getX();
T_of_x[0][1] = FReal(1.); T_of_x[1][1] = localPosition.getY();
T_of_x[0][2] = FReal(1.); T_of_x[1][2] = localPosition.getZ();
// U_0(x_i) and U_1(x_i)
U_of_x[0][0] = FReal(1.); U_of_x[1][0] = localPosition.getX() * FReal(2.);
U_of_x[0][1] = FReal(1.); U_of_x[1][1] = localPosition.getY() * FReal(2.);
U_of_x[0][2] = FReal(1.); U_of_x[1][2] = localPosition.getZ() * FReal(2.);
// U_of_x[0][0] = FReal(1.); U_of_x[1][0] = localPosition.getX() * FReal(2.);
// U_of_x[0][1] = FReal(1.); U_of_x[1][1] = localPosition.getY() * FReal(2.);
// U_of_x[0][2] = FReal(1.); U_of_x[1][2] = localPosition.getZ() * FReal(2.);
U_of_x[0][0] = FReal(1.); U_of_x[1][0] = xpx;
U_of_x[0][1] = FReal(1.); U_of_x[1][1] = ypy;
U_of_x[0][2] = FReal(1.); U_of_x[1][2] = zpz;
for (unsigned int o=2; o<ORDER; ++o) {
// T_o(x_i)
T_of_x[o][0] = FReal(2.)*localPosition.getX()*T_of_x[o-1][0] - T_of_x[o-2][0];
T_of_x[o][1] = FReal(2.)*localPosition.getY()*T_of_x[o-1][1] - T_of_x[o-2][1];
T_of_x[o][2] = FReal(2.)*localPosition.getZ()*T_of_x[o-1][2] - T_of_x[o-2][2];
// T_of_x[o][0] = FReal(2.)*localPosition.getX()*T_of_x[o-1][0] - T_of_x[o-2][0];
// T_of_x[o][1] = FReal(2.)*localPosition.getY()*T_of_x[o-1][1] - T_of_x[o-2][1];
// T_of_x[o][2] = FReal(2.)*localPosition.getZ()*T_of_x[o-1][2] - T_of_x[o-2][2];
// // U_o(x_i)
// U_of_x[o][0] = FReal(2.)*localPosition.getX()*U_of_x[o-1][0] - U_of_x[o-2][0];
// U_of_x[o][1] = FReal(2.)*localPosition.getY()*U_of_x[o-1][1] - U_of_x[o-2][1];
// U_of_x[o][2] = FReal(2.)*localPosition.getZ()*U_of_x[o-1][2] - U_of_x[o-2][2];
T_of_x[o][0] = xpx*T_of_x[o-1][0] - T_of_x[o-2][0];
T_of_x[o][1] = ypy*T_of_x[o-1][1] - T_of_x[o-2][1];
T_of_x[o][2] = zpz*T_of_x[o-1][2] - T_of_x[o-2][2];
// U_o(x_i)
U_of_x[o][0] = FReal(2.)*localPosition.getX()*U_of_x[o-1][0] - U_of_x[o-2][0];
U_of_x[o][1] = FReal(2.)*localPosition.getY()*U_of_x[o-1][1] - U_of_x[o-2][1];
U_of_x[o][2] = FReal(2.)*localPosition.getZ()*U_of_x[o-1][2] - U_of_x[o-2][2];
U_of_x[o][0] = xpx*U_of_x[o-1][0] - U_of_x[o-2][0];
U_of_x[o][1] = ypy*U_of_x[o-1][1] - U_of_x[o-2][1];
U_of_x[o][2] = zpz*U_of_x[o-1][2] - U_of_x[o-2][2];
}
// scale, because dT_o/dx = oU_{o-1}
......@@ -509,6 +559,9 @@ inline void FChebInterpolator<ORDER>::applyL2PTotal(const FPoint& center,
// apply P and increment forces
FReal potential = FReal(0.);
FReal forces[3] = {FReal(0.), FReal(0.), FReal(0.)};
//
// Optimization:
// Here we compute 1/2 S and 1/2 P rather S and F like in the paper
for (unsigned int n=0; n<nnodes; ++n) {
// tensor indices of chebyshev nodes
......@@ -527,55 +580,66 @@ inline void FChebInterpolator<ORDER>::applyL2PTotal(const FPoint& center,
// for potential
S0 += T_of_x[o][0] * T_of_roots[o][j[0]];
}
P[0] *= FReal(2.);
P[1] *= FReal(2.); P[1] += FReal(1.);
P[2] *= FReal(2.); P[2] += FReal(1.);
// P[0] *= FReal(2.);
// P[1] *= FReal(2.); P[1] += FReal(1.);
// P[2] *= FReal(2.); P[2] += FReal(1.);
// forces[0] += P[0] * P[1] * P[2] * localExpansion[n];
P[1] += FReal(0.5);
P[2] += FReal(0.5);
forces[0] += P[0] * P[1] * P[2] * localExpansion[n];
// for potential
S0 *= FReal(2.); S0 += FReal(1.);
potential += FReal(1.) / nnodes * S0 * P[1] * P[2] * localExpansion[n];
// S0 *= FReal(2.); S0 += FReal(1.);
S0 += FReal(0.5);
potential += S0 * P[1] * P[2] * localExpansion[n];
// f1 component //////////////////////////////////////
P[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
P[1] = U_of_x[0][1] * T_of_roots[1][j[1]];
P[2] = T_of_x[1][2] * T_of_roots[1][j[2]];
// P[2] = T_of_x[1][2] * T_of_roots[1][j[2]];
for (unsigned int o=2; o<ORDER; ++o) {
P[0] += T_of_x[o ][0] * T_of_roots[o][j[0]];
P[1] += U_of_x[o-1][1] * T_of_roots[o][j[1]];
P[2] += T_of_x[o ][2] * T_of_roots[o][j[2]];
// P[2] += T_of_x[o ][2] * T_of_roots[o][j[2]];
}
P[0] *= FReal(2.); P[0] += FReal(1.);
P[1] *= FReal(2.);
P[2] *= FReal(2.); P[2] += FReal(1.);
// P[0] *= FReal(2.); P[0] += FReal(1.);
// P[1] *= FReal(2.);
// P[2] *= FReal(2.); P[2] += FReal(1.);
P[0] += FReal(0.5);
// P[2] += FReal(0.5);
forces[1] += P[0] * P[1] * P[2] * localExpansion[n];
// f2 component //////////////////////////////////////
P[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
// P[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
P[1] = T_of_x[1][1] * T_of_roots[1][j[1]];
P[2] = U_of_x[0][2] * T_of_roots[1][j[2]];
for (unsigned int o=2; o<ORDER; ++o) {
P[0] += T_of_x[o ][0] * T_of_roots[o][j[0]];
// P[0] += T_of_x[o ][0] * T_of_roots[o][j[0]];
P[1] += T_of_x[o ][1] * T_of_roots[o][j[1]];
P[2] += U_of_x[o-1][2] * T_of_roots[o][j[2]];
}
P[0] *= FReal(2.); P[0] += FReal(1.);
P[1] *= FReal(2.); P[1] += FReal(1.);
P[2] *= FReal(2.);
// P[0] *= FReal(2.); P[0] += FReal(1.);
// P[1] *= FReal(2.); P[1] += FReal(1.);
// P[2] *= FReal(2.);
// P[0] += FReal(0.5);
P[1] += FReal(0.5);
forces[2] += P[0] * P[1] * P[2] * localExpansion[n];
}
// scale forces
forces[0] *= jacobian[0] / nnodes;
forces[1] *= jacobian[1] / nnodes;
forces[2] *= jacobian[2] / nnodes;
// forces[0] *= jacobian[0] / nnodes;
// forces[1] *= jacobian[1] / nnodes;
// forces[2] *= jacobian[2] / nnodes;
forces[0] *= jacobian[0] *c1;
forces[1] *= jacobian[1] *c1;
forces[2] *= jacobian[2] *c1;
potential *= c1 ;
// set computed potential
iter.data().incPotential(potential);
// set computed forces
iter.data().incForces(forces[0] * iter.data().getPhysicalValue(),
forces[1] * iter.data().getPhysicalValue(),
forces[2] * iter.data().getPhysicalValue());
forces[1] * iter.data().getPhysicalValue(),
forces[2] * iter.data().getPhysicalValue());
// increment iterator
iter.gotoNext();
......
......@@ -35,170 +35,204 @@
/*
In this test we compare the spherical fmm results and the direct results.
*/
*/
/** We need to know the position of the particle in the array */
class IndexedParticle : public FChebParticle {
int index;
int index;
public:
IndexedParticle(): index(-1){}
int getIndex() const{
return index;
}
void setIndex( const int inIndex ){
index = inIndex;
}
IndexedParticle(): index(-1){}
int getIndex() const{
return index;
}
void setIndex( const int inIndex ){
index = inIndex;
}
};
/** the test class
*
*/
*
*/
class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
///////////////////////////////////////////////////////////
// The tests!
///////////////////////////////////////////////////////////
/** Classic */
void TestChebyshev(){
const unsigned int ORDER = 5;
// typedefs
typedef IndexedParticle ParticleClass;
typedef FVector<ParticleClass> ContainerClass;
typedef FChebLeaf<ParticleClass,ContainerClass> LeafClass;
typedef FChebMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
typedef FOctree<ParticleClass,CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebKernel<ParticleClass,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
//typedef FChebSymKernel<ParticleClass,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FFmmAlgorithm<OctreeClass,ParticleClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// Warning in make test the exec dir it Build/UTests
// Load particles
FFmaBinLoader<ParticleClass> loader("../../Data/utestSphericalDirect.bin.fma");
if(!loader.isOpen()){
Print("Cannot open particles file.");
uassert(false);
return;
}
Print("Number of particles:");
Print(loader.getNumberOfParticles());
const int NbLevels = 4;
const int SizeSubLevels = 2;
const FReal epsilon = FReal(1e-5);
// Create octree
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
ParticleClass* const particles = new ParticleClass[loader.getNumberOfParticles()];
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
loader.fillParticle(particles[idxPart]);
particles[idxPart].setIndex( idxPart );
tree.insert(particles[idxPart]);
}
// Run FMM
Print("Fmm...");
KernelClass kernels(NbLevels, loader.getCenterOfBox(), loader.getBoxWidth(), epsilon);
FmmClass algo(&tree,&kernels);
algo.execute();
// Run direct computation
const MatrixKernelClass MatrixKernel;
Print("Direct...");
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
//kernels.directInteractionMutual(&particles[idxTarget], &particles[idxOther]);
const FReal wt = particles[idxTarget].getPhysicalValue();
const FReal ws = particles[idxOther ].getPhysicalValue();
const FReal one_over_r = MatrixKernel.evaluate(particles[idxTarget].getPosition(),
particles[idxOther].getPosition());
// potential
particles[idxTarget].incPotential(one_over_r * ws);
particles[idxOther ].incPotential(one_over_r * wt);
// force
F3DPosition force(particles[idxOther].getPosition() - particles[idxTarget].getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
particles[idxTarget].incForces( force.getX(), force.getY(), force.getZ());
particles[idxOther ].incForces( -force.getX(), -force.getY(), -force.getZ());
}
}
// Compare
Print("Compute Diff...");
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
{ // Check that each particle has been summed with all other
OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
ContainerClass::BasicIterator leafIter(*octreeIterator.getCurrentListTargets());
while( leafIter.hasNotFinished() ){
const ParticleClass& other = particles[leafIter.data().getIndex()];
potentialDiff.add(other.getPotential(),leafIter.data().getPotential());
fx.add(other.getForces().getX(),leafIter.data().getForces().getX());
fy.add(other.getForces().getY(),leafIter.data().getForces().getY());
fz.add(other.getForces().getZ(),leafIter.data().getForces().getZ());
leafIter.gotoNext();
}
} while(octreeIterator.moveRight());
}
delete[] particles;
// Print for information
Print("Potential diff is = ");
Print(potentialDiff.getL2Norm());
Print(potentialDiff.getInfNorm());
Print("Fx diff is = ");
Print(fx.getL2Norm());
Print(fx.getInfNorm());
Print("Fy diff is = ");
Print(fy.getL2Norm());
Print(fy.getInfNorm());
Print("Fz diff is = ");
Print(fz.getL2Norm());
Print(fz.getInfNorm());
// Assert
const FReal MaximumDiffPotential = FReal(1e-5);
const FReal MaximumDiffForces = FReal(1e-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);
}
/** If memstas is running print the memory used */
void PostTest() {
if( FMemStats::controler.isUsed() ){
std::cout << "Memory used at the end " << FMemStats::controler.getCurrentAllocated() << " Bytes (" << FMemStats::controler.getCurrentAllocatedMB() << "MB)\n";
std::cout << "Max memory used " << FMemStats::controler.getMaxAllocated() << " Bytes (" << FMemStats::controler.getMaxAllocatedMB() << "MB)\n";
std::cout << "Total memory used " << FMemStats::controler.getTotalAllocated() << " Bytes (" << FMemStats::controler.getTotalAllocatedMB() << "MB)\n";
}
}
///////////////////////////////////////////////////////////
// Set the tests!
///////////////////////////////////////////////////////////
/** set test */
void SetTests(){
AddTest(&TestChebyshevDirect::TestChebyshev,"Test Chebyshev Kernel");
}
///////////////////////////////////////////////////////////
// The tests!
///////////////////////////////////////////////////////////
template <class ParticleClass, class CellClass, class ContainerClass, class KernelClass, class MatrixKernelClass,
class LeafClass, class OctreeClass, class FmmClass>
void RunTest(const FReal epsilon) {
// Warning in make test the exec dir it Build/UTests
// Load particles
FFmaBinLoader<ParticleClass> loader("../../Data/utestSphericalDirect.bin.fma");
if(!loader.isOpen()){
Print("Cannot open particles file.");
uassert(false);
return;
}
Print("Number of particles:");
Print(loader.getNumberOfParticles());
const int NbLevels = 4;
const int SizeSubLevels = 2;
//const FReal epsilon = FReal(1e-5);
// Create octree
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
ParticleClass* const particles = new ParticleClass[loader.getNumberOfParticles()];
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
loader.fillParticle(particles[idxPart]);
particles[idxPart].setIndex( idxPart );
tree.insert(particles[idxPart]);
}
// Run FMM
Print("Fmm...");
KernelClass kernels(NbLevels, loader.getCenterOfBox(), loader.getBoxWidth(), epsilon);
FmmClass algo(&tree,&kernels);
algo.execute();
// Run direct computation
const MatrixKernelClass MatrixKernel;
Print("Direct...");
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
//kernels.directInteractionMutual(&particles[idxTarget], &particles[idxOther]);
const FReal wt = particles[idxTarget].getPhysicalValue();
const FReal ws = particles[idxOther ].getPhysicalValue();
const FReal one_over_r = MatrixKernel.evaluate(particles[idxTarget].getPosition(),
particles[idxOther].getPosition());
// potential
particles[idxTarget].incPotential(one_over_r * ws);
particles[idxOther ].incPotential(one_over_r * wt);
// force
FPoint force(particles[idxOther].getPosition() - particles[idxTarget].getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
particles[idxTarget].incForces( force.getX(), force.getY(), force.getZ());
particles[idxOther ].incForces( -force.getX(), -force.getY(), -force.getZ());
}
}
// Compare
Print("Compute Diff...");
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
{ // Check that each particle has been summed with all other
typename OctreeClass::Iterator octreeIterator(&tree);
octreeIterator.gotoBottomLeft();
do{
typename ContainerClass::BasicIterator leafIter(*octreeIterator.getCurrentListTargets());
while( leafIter.hasNotFinished() ){
const ParticleClass& other = particles[leafIter.data().getIndex()];
potentialDiff.add(other.getPotential(),leafIter.data().getPotential());
fx.add(other.getForces().getX(),leafIter.data().getForces().getX());
fy.add(other.getForces().getY(),leafIter.data().getForces().getY());
fz.add(other.getForces().getZ(),leafIter.data().getForces().getZ());
leafIter.gotoNext();
}
} while(octreeIterator.moveRight());
}
delete[] particles;
// Print for information