diff --git a/CMakeLists.txt b/CMakeLists.txt
index e19f68eeae38e8bd7557cfa270255b08f7b280dc..928896831a695623918273d99d62a2e15d4a9143 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -109,9 +109,9 @@ if (MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/CMakeModules/morse/
option( OPENMP_SUPPORT_PRIORITY "Set to ON to enable tasks priority (KSTAR/StarPU compiler only)" OFF )
option( OPENMP_SUPPORT_TASK_NAME "Set to ON to enable a taskname clause for tasks (KSTAR/StarPU compiler only)" OFF )
if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
- option( SCALFMM_DISABLE_NATIVE_OMP4 "Set to ON to disable the gcc/intel omp4" OFF )
+ option( SCALFMM_USE_OMP4 "Set to ON to disable the gcc/intel omp4" OFF )
else()
- option( SCALFMM_DISABLE_NATIVE_OMP4 "Set to ON to disable the gcc/intel omp4" ON )
+ option( SCALFMM_USE_OMP4 "Set to ON to disable the gcc/intel omp4" ON )
endif()
option( SCALFMM_TIME_OMPTASKS "Set to ON to time omp4 tasks and generate output file" OFF )
# SIMGRID and peformance models options
@@ -616,6 +616,7 @@ if (MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/CMakeModules/morse/
list(APPEND FUSE_LIST "STARPU")
list(APPEND FUSE_LIST "CUDA")
list(APPEND FUSE_LIST "OPENCL")
+ list(APPEND FUSE_LIST "OMP4")
##################################################################
# Use SSE #
@@ -959,6 +960,7 @@ if (MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/CMakeModules/morse/
message(STATUS "SCALFMM_CXX_FLAGS = ${SCALFMM_CXX_FLAGS}")
message(STATUS "SCALFMM_LIBRARIES = ${SCALFMM_LIBRARIES}")
message(STATUS "SCALFMM_INCLUDES = ${SCALFMM_INCLUDES}")
+ message(STATUS "FUSE_LIST = ${FUSE_LIST}")
#
##################################################################
# END #
diff --git a/ScalFMMConfig.cmake.in b/ScalFMMConfig.cmake.in
index bc683cab353236228569d6264bb3f4c1c8e69c9d..bdce4ead7daea776476dfaaa2e28f68de559d373 100644
--- a/ScalFMMConfig.cmake.in
+++ b/ScalFMMConfig.cmake.in
@@ -27,18 +27,18 @@ SET(SCALFMM_BUILD_TYPE "@CMAKE_BUILD_TYPE@")
#
# SCALFMM Options
#
+SET(SCALFMM_USE_ADDONS "@SCALFMM_USE_ADDONS@")
SET(SCALFMM_USE_MPI "@SCALFMM_USE_MPI@")
SET(SCALFMM_USE_BLAS "@SCALFMM_USE_BLAS@")
SET(SCALFMM_USE_FFT "@SCALFMM_USE_FFT@")
SET(SCALFMM_USE_MKL "@SCALFMM_USE_MKL_AS_BLAS@")
-SET(SCALFMM_USE_MEM_STATS "@SCALFMM_USE_MEM_STATS@")
-SET(SCALFMM_USE_SSE "@SCALFMM_USE_SSE@")
SET(SCALFMM_FLAGS "@SCALFMM_FLAGS_OPTI@")
-SET(SCALFMM_USE_ADDONS "@SCALFMM_USE_ADDONS@")
+SET(SCALFMM_USE_MEM_STATS "@SCALFMM_USE_MEM_STATS@")
SET(SCALFMM_USE_LOG "@SCALFMM_USE_LOG@")
+SET(SCALFMM_USE_OMP4 "@SCALFMM_USE_OMP4@")
SET(SCALFMM_USE_STARPU "@SCALFMM_USE_STARPU@")
-SET(SCALFMM_USE_SSE "@SCALFMM_USE_SSE@")
-SET(SCALFMM_USE_AVX "@SCALFMM_USE_AVX@")
+SET(SCALFMM_USE_SSE "@SCALFMM_USE_SSE@")
+SET(SCALFMM_USE_AVX "@SCALFMM_USE_AVX@")
#
SET(SCALFMM_DOC_TAGS "@CMAKE_BINARY_DIR@/Doc/scalfmm.tag")
diff --git a/Tests/Utils/testChebInterpolator.cpp b/Tests/Utils/testChebInterpolator.cpp
index 4bc537b72275f90ba8b466da68067604eb28e150..756e7a7e687ffafe1b2a34f271520c6e620fcaf6 100644
--- a/Tests/Utils/testChebInterpolator.cpp
+++ b/Tests/Utils/testChebInterpolator.cpp
@@ -16,6 +16,7 @@
// ==== CMAKE =====
// @FUSE_BLAS
+// @FUSE_OMP4
// ================
diff --git a/Tests/Utils/testFmmAlgorithmOmp4.cpp b/Tests/Utils/testFmmAlgorithmOmp4.cpp
index 54bdef7049fe43eeb0c293cc5efb7693aefa297d..ba7a8a658bbe7a8be315c56302b4180574177412 100644
--- a/Tests/Utils/testFmmAlgorithmOmp4.cpp
+++ b/Tests/Utils/testFmmAlgorithmOmp4.cpp
@@ -18,7 +18,7 @@
// "http://www.gnu.org/licenses".
// ===================================================================================
-// @SCALFMM_PRIVATE
+// @ SCALFMM_PRIVATE
// @FUSE_OMP4
diff --git a/UTests/FUKernelTester.hpp b/UTests/FUKernelTester.hpp
index 0335d01981e491258327c750ed88aa005bf5bab2..ec98ca93c83e529a10cc9f5bc452a4b1671e0e3b 100644
--- a/UTests/FUKernelTester.hpp
+++ b/UTests/FUKernelTester.hpp
@@ -96,7 +96,7 @@ public:
// Run FMM computation
/////////////////////////////////////////////////////////////////////////////////////////////////
Print("Fmm...");
- std::unique_ptr<KernelClass> kernels(GetKernelFunc(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel));
+ std::unique_ptr<KernelClass> kernels(GetKernelFunc(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel));
FmmClass algo(&tree,kernels.get());
algo.execute();
//
@@ -118,12 +118,12 @@ public:
tree.forEachLeaf([&](LeafClass* leaf){
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();
- const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
- const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
+ 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 FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
+ const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize indexPartOrig = indexes[idxPart];
diff --git a/UTests/utestChebyshev.cpp b/UTests/utestChebyshev.cpp
index 8779ed36448fa558a96ba020234ec998324c01d6..9c5b4476bb18d1b45960940d8de91c81f11b6ec6 100644
--- a/UTests/utestChebyshev.cpp
+++ b/UTests/utestChebyshev.cpp
@@ -34,65 +34,65 @@
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
/*
In this test we compare the Chebyschev fmm results and the direct results.
- */
+*/
/** the test class
*
*/
class TestChebyshevDirect : public FUKernelTester<TestChebyshevDirect> {
- ///////////////////////////////////////////////////////////
- // Set the tests!
- ///////////////////////////////////////////////////////////
-
-
- /** TestChebKernel */
- void TestChebKernel(){
- typedef double FReal;
- const unsigned int ORDER = 6;
- typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
- typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
- typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
- typedef FChebCell<FReal,ORDER> CellClass;
- typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
- typedef FChebKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
- typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
- // run test
+ ///////////////////////////////////////////////////////////
+ // Set the tests!
+ ///////////////////////////////////////////////////////////
+
+
+ /** TestChebKernel */
+ void TestChebKernel(){
+ typedef double FReal;
+ const unsigned int ORDER = 6;
+ typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
+ typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
+ typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
+ typedef FChebCell<FReal,ORDER> CellClass;
+ typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
+ typedef FChebKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
+ typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
+ // run test
RunTest<FReal,CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(
- [&](int NbLevels, FReal boxWidth, FPoint<FReal> centerOfBox, const MatrixKernelClass *const MatrixKernel){
- return std::unique_ptr<KernelClass>(new KernelClass(NbLevels, boxWidth, centerOfBox, MatrixKernel));
- });
- }
-
- /** TestChebSymKernel */
- void TestChebSymKernel(){
- typedef double FReal;
- const unsigned int ORDER = 6;
- typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
- typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
- typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
- typedef FChebCell<FReal,ORDER> CellClass;
- typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
- typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
- typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
- // run test
+ [&](int NbLevels, FReal boxWidth, FPoint<FReal> centerOfBox, const MatrixKernelClass *const MatrixKernel){
+ return std::unique_ptr<KernelClass>(new KernelClass(NbLevels, boxWidth, centerOfBox, MatrixKernel));
+ });
+ }
+
+ /** TestChebSymKernel */
+ void TestChebSymKernel(){
+ typedef double FReal;
+ const unsigned int ORDER = 6;
+ typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
+ typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
+ typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
+ typedef FChebCell<FReal,ORDER> CellClass;
+ typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
+ typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
+ typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
+ // run test
RunTest<FReal,CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>(
- [&](int NbLevels, FReal boxWidth, FPoint<FReal> centerOfBox, const MatrixKernelClass *const MatrixKernel){
- return std::unique_ptr<KernelClass>(new KernelClass(NbLevels, boxWidth, centerOfBox, MatrixKernel));
- });
- }
+ [&](int NbLevels, FReal boxWidth, FPoint<FReal> centerOfBox, const MatrixKernelClass *const MatrixKernel){
+ return std::unique_ptr<KernelClass>(new KernelClass(NbLevels, boxWidth, centerOfBox, MatrixKernel));
+ });
+ }
- ///////////////////////////////////////////////////////////
- // Set the tests!
- ///////////////////////////////////////////////////////////
+ ///////////////////////////////////////////////////////////
+ // Set the tests!
+ ///////////////////////////////////////////////////////////
- /** 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");
- }
+ /** 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");
+ }
};
diff --git a/UTests/utestLagrange.cpp b/UTests/utestLagrange.cpp
index c336bcb7eb9518a29b076453d7c884b380695729..5eb2784262790fd66855f870355e0da22366061a 100644
--- a/UTests/utestLagrange.cpp
+++ b/UTests/utestLagrange.cpp
@@ -45,7 +45,7 @@
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
/*
In this test we compare the TestLagrange fmm results with the direct results.
- */
+*/
/** the test class
@@ -53,191 +53,191 @@
*/
class TestLagrange : public FUTester<TestLagrange> {
- ///////////////////////////////////////////////////////////
- // The tests!
- ///////////////////////////////////////////////////////////
-
- template <class FReal, class CellClass, class ContainerClass, class KernelClass, class MatrixKernelClass,
- class LeafClass, class OctreeClass, class FmmClass>
- void RunTest() {
- //
- // Load particles
- //
- if(sizeof(FReal) == sizeof(float) ) {
- std::cerr << "No input data available for Float "<< std::endl;
- exit(EXIT_FAILURE);
- }
- const std::string parFile( (sizeof(FReal) == sizeof(float))?
- "Test/DirectFloatbfma":
- "UTest/DirectDouble.bfma");
- //
- std::string filename(SCALFMMDataPath+parFile);
- //
- FFmaGenericLoader<FReal> loader(filename);
- Print("Number of particles:");
- Print(loader.getNumberOfParticles());
-
- const int NbLevels = 4;
- const int SizeSubLevels = 2;
+ ///////////////////////////////////////////////////////////
+ // The tests!
+ ///////////////////////////////////////////////////////////
- // Create Matrix Kernel
- const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
+ template <class FReal, class CellClass, class ContainerClass, class KernelClass, class MatrixKernelClass,
+ class LeafClass, class OctreeClass, class FmmClass>
+ void RunTest() {
+ //
+ // Load particles
+ //
+ if(sizeof(FReal) == sizeof(float) ) {
+ std::cerr << "No input data available for Float "<< std::endl;
+ exit(EXIT_FAILURE);
+ }
+ const std::string parFile( (sizeof(FReal) == sizeof(float))?
+ "Test/DirectFloatbfma":
+ "UTest/DirectDouble.bfma");
+ //
+ std::string filename(SCALFMMDataPath+parFile);
+ //
+ FFmaGenericLoader<FReal> loader(filename);
+ Print("Number of particles:");
+ Print(loader.getNumberOfParticles());
+
+ const int NbLevels = 4;
+ const int SizeSubLevels = 2;
// Load particles
- FSize nbParticles = loader.getNumberOfParticles() ;
- FmaRWParticle<FReal, 8,8>* const particles = new FmaRWParticle<FReal, 8,8>[nbParticles];
+ FSize nbParticles = loader.getNumberOfParticles() ;
+ FmaRWParticle<FReal, 8,8>* const particles = new FmaRWParticle<FReal, 8,8>[nbParticles];
- loader.fillParticle(particles,nbParticles);
+ loader.fillParticle(particles,nbParticles);
- // Create octree
- OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+ // Create octree
+ OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
+ //
+ // Insert particle in the tree
+ //
+ for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
+ tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue() );
+ }
+ //
+ // Create Matrix Kernel
+ const MatrixKernelClass MatrixKernel; // FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
//
- // Insert particle in the tree
- //
- for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
- tree.insert(particles[idxPart].getPosition(), idxPart, particles[idxPart].getPhysicalValue() );
- }
- //
- /////////////////////////////////////////////////////////////////////////////////////////////////
- // Run FMM computation
- /////////////////////////////////////////////////////////////////////////////////////////////////
- Print("Fmm...");
- KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
- FmmClass algo(&tree,&kernels);
- algo.execute();
- //0
- FReal energy= 0.0 , energyD = 0.0 ;
- /////////////////////////////////////////////////////////////////////////////////////////////////
- // Compute direct energy
- /////////////////////////////////////////////////////////////////////////////////////////////////
-
- for(FSize idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){
- energyD += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
- }
- /////////////////////////////////////////////////////////////////////////////////////////////////
- // Compare
- /////////////////////////////////////////////////////////////////////////////////////////////////
- Print("Compute Diff...");
- FMath::FAccurater<FReal> potentialDiff;
- FMath::FAccurater<FReal> fx, fy, fz;
- { // Check that each particle has been summed with all other
-
- tree.forEachLeaf([&](LeafClass* leaf){
- 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();
- const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
- const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
-
- for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
- const FSize indexPartOrig = indexes[idxPart];
- potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
- fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
- fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
- fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
- energy += potentials[idxPart]*physicalValues[idxPart];
- }
- });
- }
-
- delete[] particles;
-
- // Print for information
-
- 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 DIRECT = %.12e\n",FMath::Abs(energyD));
-
- // Assert
- const FReal MaximumDiffPotential = FReal(9e-3);
- const FReal MaximumDiffForces = FReal(9e-2);
-
- Print("Test1 - Error Relative L2 norm Potential ");
- uassert(potentialDiff.getRelativeL2Norm() < MaximumDiffPotential); //1
- Print("Test2 - Error RMS L2 norm Potential ");
- uassert(potentialDiff.getRMSError() < MaximumDiffPotential); //2
- Print("Test3 - Error Relative L2 norm FX ");
- uassert(fx.getRelativeL2Norm() < MaximumDiffForces); //3
- Print("Test4 - Error RMS L2 norm FX ");
- uassert(fx.getRMSError() < MaximumDiffForces); //4
- Print("Test5 - Error Relative L2 norm FY ");
- uassert(fy.getRelativeL2Norm() < MaximumDiffForces); //5
- Print("Test6 - Error RMS L2 norm FY ");
- uassert(fy.getRMSError() < MaximumDiffForces); //6
- Print("Test7 - Error Relative L2 norm FZ ");
- uassert(fz.getRelativeL2Norm() < MaximumDiffForces); //8
- Print("Test8 - Error RMS L2 norm FZ ");
- uassert(fz.getRMSError() < MaximumDiffForces); //8
- Print("Test9 - Error Relative L2 norm F ");
- uassert(L2error < MaximumDiffForces); //9 Total Force
- Print("Test10 - Relative error Energy ");
- uassert(FMath::Abs(energy-energyD) /energyD< MaximumDiffPotential); //10 Total Energy
-
-
- }
-
- /** 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!
- ///////////////////////////////////////////////////////////
-
-
- /** TestUnifKernel */
- void TestUnifKernel(){
- typedef double FReal;
- const unsigned int ORDER = 6;
- // typedefs
- typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
- typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
- typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
- typedef FUnifCell<FReal,ORDER> CellClass;
- typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
- typedef FUnifKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
- typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
- // run test
- RunTest<FReal,CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>();
- }
-
- ///////////////////////////////////////////////////////////
- // Set the tests!
- ///////////////////////////////////////////////////////////
-
- /** set test */
- void SetTests(){
- AddTest(&TestLagrange::TestUnifKernel,"Test Lagrange Kernel ");
- }
+ /////////////////////////////////////////////////////////////////////////////////////////////////
+ // Run FMM computation
+ /////////////////////////////////////////////////////////////////////////////////////////////////
+ Print("Fmm...");
+ KernelClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
+ FmmClass algo(&tree,&kernels);
+ algo.execute();
+ //0
+ FReal energy= 0.0 , energyD = 0.0 ;
+ /////////////////////////////////////////////////////////////////////////////////////////////////
+ // Compute direct energy
+ /////////////////////////////////////////////////////////////////////////////////////////////////
+
+ for(FSize idx = 0 ; idx < loader.getNumberOfParticles() ; ++idx){
+ energyD += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
+ }
+ /////////////////////////////////////////////////////////////////////////////////////////////////
+ // Compare
+ /////////////////////////////////////////////////////////////////////////////////////////////////
+ Print("Compute Diff...");
+ FMath::FAccurater<FReal> potentialDiff;
+ FMath::FAccurater<FReal> fx, fy, fz;
+ { // Check that each particle has been summed with all other
+
+ tree.forEachLeaf([&](LeafClass* leaf){
+ 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();
+ const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
+ const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
+
+ for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
+ const FSize indexPartOrig = indexes[idxPart];
+ potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
+ fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
+ fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
+ fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
+ energy += potentials[idxPart]*physicalValues[idxPart];
+ }
+ });
+ }
+
+ delete[] particles;
+
+ // Print for information
+
+ 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 DIRECT = %.12e\n",FMath::Abs(energyD));
+
+ // Assert
+ const FReal MaximumDiffPotential = FReal(9e-3);
+ const FReal MaximumDiffForces = FReal(9e-2);
+
+ Print("Test1 - Error Relative L2 norm Potential ");
+ uassert(potentialDiff.getRelativeL2Norm() < MaximumDiffPotential); //1
+ Print("Test2 - Error RMS L2 norm Potential ");
+ uassert(potentialDiff.getRMSError() < MaximumDiffPotential); //2
+ Print("Test3 - Error Relative L2 norm FX ");
+ uassert(fx.getRelativeL2Norm() < MaximumDiffForces); //3
+ Print("Test4 - Error RMS L2 norm FX ");
+ uassert(fx.getRMSError() < MaximumDiffForces); //4
+ Print("Test5 - Error Relative L2 norm FY ");
+ uassert(fy.getRelativeL2Norm() < MaximumDiffForces); //5
+ Print("Test6 - Error RMS L2 norm FY ");
+ uassert(fy.getRMSError() < MaximumDiffForces); //6
+ Print("Test7 - Error Relative L2 norm FZ ");
+ uassert(fz.getRelativeL2Norm() < MaximumDiffForces); //8
+ Print("Test8 - Error RMS L2 norm FZ ");
+ uassert(fz.getRMSError() < MaximumDiffForces); //8
+ Print("Test9 - Error Relative L2 norm F ");
+ uassert(L2error < MaximumDiffForces); //9 Total Force
+ Print("Test10 - Relative error Energy ");
+ uassert(FMath::Abs(energy-energyD) /energyD< MaximumDiffPotential); //10 Total Energy
+
+
+ }
+
+ /** 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!
+ ///////////////////////////////////////////////////////////
+
+
+ /** TestUnifKernel */
+ void TestUnifKernel(){
+ typedef double FReal;
+ const unsigned int ORDER = 6;
+ // typedefs
+ typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
+ typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
+ typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
+ typedef FUnifCell<FReal,ORDER> CellClass;
+ typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
+ typedef FUnifKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
+ typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
+ // run test
+ RunTest<FReal,CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>();
+ }
+
+ ///////////////////////////////////////////////////////////
+ // Set the tests!
+ ///////////////////////////////////////////////////////////
+
+ /** set test */
+ void SetTests(){
+ AddTest(&TestLagrange::TestUnifKernel,"Test Lagrange Kernel ");
+ }
};