From 2c1cdef1fdb334d89bc45bddc048ba4710b26335 Mon Sep 17 00:00:00 2001
From: piacibel <cyrille.piacibello@inria.fr>
Date: Wed, 6 May 2015 18:33:11 +0200
Subject: [PATCH] Source Target algorithm is now functionnal. An example can be
 found in testSphereElectro.c

---
 Addons/CKernelApi/Src/FInterEngine.hpp      |  18 ++-
 Addons/CKernelApi/Src/FScalfmmApiInit.cpp   |   7 +-
 Addons/CKernelApi/Src/FUserKernelEngine.hpp |   2 +
 Addons/CKernelApi/Tests/Timers.h            |  54 +++++++
 Addons/CKernelApi/Tests/testChebInterface.c |  51 +------
 Addons/CKernelApi/Tests/testSphereElectro.c | 148 ++++++++++++++------
 6 files changed, 183 insertions(+), 97 deletions(-)
 create mode 100644 Addons/CKernelApi/Tests/Timers.h

diff --git a/Addons/CKernelApi/Src/FInterEngine.hpp b/Addons/CKernelApi/Src/FInterEngine.hpp
index b2c75cee..f0fe11d4 100644
--- a/Addons/CKernelApi/Src/FInterEngine.hpp
+++ b/Addons/CKernelApi/Src/FInterEngine.hpp
@@ -31,6 +31,7 @@
 #include "Arranger/FArrangerPeriodic.hpp"
 #include "Arranger/FBasicParticleContainerIndexedMover.hpp"
 #include "Arranger/FParticleTypedIndexedMover.hpp"
+#include "Extensions/FExtendCellType.hpp"
 
 #include "Core/FFmmAlgorithmThread.hpp"
 #include "Core/FFmmAlgorithm.hpp"
@@ -38,6 +39,7 @@
 #include "Core/FFmmAlgorithmThreadTsm.hpp"
 
 
+
 /**
  * @class FInterEngine implements API for Interpolations kernels, its
  * templates can be ChebCell/ChebKernel or UnifCell/UnifKernel
@@ -63,6 +65,7 @@ private:
     // ArrangerClass * arranger;
 
 
+
 public:
     /**
      * @brief Constructor : build the tree and the interpolation
@@ -72,9 +75,10 @@ public:
      * @param BoxCenter double[3] coordinate of the center of the
      * simulation box
      */
-    FInterEngine(scalfmm_kernel_type KernelType) :
+    FInterEngine(scalfmm_kernel_type KernelType, scalfmm_algorithm algo) :
         kernel(nullptr), matrix(nullptr), octree(nullptr)/*,arranger(nullptr)*/{
         FScalFMMEngine<FReal>::kernelType = KernelType;
+        FScalFMMEngine<FReal>::Algorithm = algo;
     }
 
     void build_tree(int TreeHeight, FReal BoxWidth , FReal * BoxCenter,User_Scalfmm_Cell_Descriptor notUsedHere){
@@ -760,11 +764,13 @@ public:
             }
         case 3:
             {
-                // typedef FFmmAlgorithmThreadTsm<OctreeClass,InterCell,ContainerClass,InterKernel,LeafClass> AlgoClassTargetSource;
-                // AlgoClassTargetSource* algoTS = new AlgoClassTargetSource(octree,kernel);
-                // algoTS->execute();
-                // FScalFMMEngine<FReal>::algoTimer = algoTS;
-                // break;
+                // class local : public InterCell, public unExtendedCell{
+                // };
+                typedef FFmmAlgorithmThreadTsm<OctreeClass,InterCell,ContainerClass,InterKernel,LeafClass> AlgoClassTargetSource;
+                AlgoClassTargetSource* algoTS = new AlgoClassTargetSource(octree,kernel);
+                algoTS->execute();
+                FScalFMMEngine<FReal>::algoTimer = algoTS;
+                break;
             }
         default :
             std::cout<< "No algorithm found (probably for strange reasons) : "<< FScalFMMEngine<FReal>::Algorithm <<" exiting" << std::endl;
diff --git a/Addons/CKernelApi/Src/FScalfmmApiInit.cpp b/Addons/CKernelApi/Src/FScalfmmApiInit.cpp
index 19cfb08b..f286a4a7 100644
--- a/Addons/CKernelApi/Src/FScalfmmApiInit.cpp
+++ b/Addons/CKernelApi/Src/FScalfmmApiInit.cpp
@@ -30,7 +30,7 @@ extern "C" scalfmm_handle scalfmm_init(/*int TreeHeight,double BoxWidth,double*
             typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
             typedef FChebSymKernel<FReal,ChebCell,ContainerClass,MatrixKernelClass,7>        ChebKernel;
 
-            handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
+            handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType,algo);
             break;
             // case 2:
             //     //TODO typedefs
@@ -61,13 +61,14 @@ extern "C" scalfmm_handle scalfmm_init(/*int TreeHeight,double BoxWidth,double*
         case 1:
             //TODO typedefs
             typedef FP2PParticleContainerIndexed<FReal>                                 ContainerClass;
-            typedef FChebCell<FReal,7>                                                   ChebCell;
+            //typedef FChebCell<FReal,7>                                                   ChebCell;
+            typedef FTypedChebCell<FReal,7>                                                   ChebCell;
             typedef FSimpleLeaf<FReal,ContainerClass>                                         LeafClass;
 
             typedef FInterpMatrixKernelR<FReal>                                        MatrixKernelClass;
             typedef FChebSymKernel<FReal,ChebCell,ContainerClass,MatrixKernelClass,7>        ChebKernel;
 
-            handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
+            handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType,algo);
             break;
             // case 2:
             //     //TODO typedefs
diff --git a/Addons/CKernelApi/Src/FUserKernelEngine.hpp b/Addons/CKernelApi/Src/FUserKernelEngine.hpp
index 834ee6a4..7004a7fb 100644
--- a/Addons/CKernelApi/Src/FUserKernelEngine.hpp
+++ b/Addons/CKernelApi/Src/FUserKernelEngine.hpp
@@ -108,6 +108,7 @@ public:
         if(kernel.m2m){
             for(int idx = 0 ; idx < 8 ; ++idx){
                 if( children[idx] ){
+                    printf("lvl : %d\n",level);
                     kernel.m2m(level, cell->getContainer(), idx, children[idx]->getContainer(), userData);
                 }
             }
@@ -253,6 +254,7 @@ public:
 
     void build_tree(int TreeHeight,double BoxWidth,double* BoxCenter,Scalfmm_Cell_Descriptor user_cell_descriptor){
         CoreCell::Init(user_cell_descriptor);
+        printf("Tree Height : %d \n",TreeHeight);
         this->octree = new OctreeClass(TreeHeight,FMath::Min(3,TreeHeight-1),BoxWidth,FPoint<FReal>(BoxCenter));
     }
 
diff --git a/Addons/CKernelApi/Tests/Timers.h b/Addons/CKernelApi/Tests/Timers.h
new file mode 100644
index 00000000..3be381a4
--- /dev/null
+++ b/Addons/CKernelApi/Tests/Timers.h
@@ -0,0 +1,54 @@
+#ifndef TIMERS_H
+#define TIMERS_H
+
+#include <time.h>
+#include <sys/time.h>
+
+
+/**
+ * @brief Wrapper on timeval struct
+ */
+typedef struct timer{
+    struct timeval start;
+    struct timeval end;
+}Timer;
+
+/**
+ * @brief monitoring function (init start)
+ */
+void tic(Timer* inTimer){
+    gettimeofday(&inTimer->start, NULL);
+}
+
+/**
+ * @brief monitoring function (init end)
+ */
+void tac(Timer* inTimer){
+    gettimeofday(&inTimer->end, NULL);
+}
+
+/**
+ * @brief monitoring function (return elapsed time in micro second)
+ */
+long int get_elapsed(Timer* inTimer){
+    return ((inTimer->end.tv_sec * 1000000 + inTimer->end.tv_usec)
+            - (inTimer->start.tv_sec * 1000000 + inTimer->start.tv_usec));
+}
+
+/**
+ * @brief monitoring function (print elapsed)
+ */
+void print_elapsed(Timer* inTimer){
+    long int elapsed = get_elapsed(inTimer);
+    printf("Elapsed : %ld us (or %f seconds)\n",elapsed,elapsed/1000000.0);
+}
+
+/**
+ * @brief monitoring function (print difference :: First-Second)
+ */
+void print_difference_elapsed(Timer* inTimer1,Timer*inTimer2){
+    long int diff = get_elapsed(inTimer1)-get_elapsed(inTimer2);
+    printf("Timer Difference : %ld us (or %f seconds)\n",diff,(double)diff/1000000.0);
+}
+
+#endif
diff --git a/Addons/CKernelApi/Tests/testChebInterface.c b/Addons/CKernelApi/Tests/testChebInterface.c
index 0026a20a..1198b37c 100644
--- a/Addons/CKernelApi/Tests/testChebInterface.c
+++ b/Addons/CKernelApi/Tests/testChebInterface.c
@@ -4,8 +4,7 @@
 #include <math.h>
 
 //For timing monitoring
-#include <time.h>
-#include <sys/time.h>
+#include "Timers.h"
 
 #include "../Src/CScalfmmApi.h"
 
@@ -60,52 +59,6 @@ void cheb_resetCell(int level, long long morton_index, int* tree_position, doubl
 }
 
 
-/**
- * @brief Wrapper on timeval struct
- */
-typedef struct timer{
-    struct timeval start;
-    struct timeval end;
-}Timer;
-
-/**
- * @brief monitoring function (init start)
- */
-void tic(Timer* inTimer){
-    gettimeofday(&inTimer->start, NULL);
-}
-
-/**
- * @brief monitoring function (init end)
- */
-void tac(Timer* inTimer){
-    gettimeofday(&inTimer->end, NULL);
-}
-
-/**
- * @brief monitoring function (return elapsed time in micro second)
- */
-long int get_elapsed(Timer* inTimer){
-    return ((inTimer->end.tv_sec * 1000000 + inTimer->end.tv_usec)
-            - (inTimer->start.tv_sec * 1000000 + inTimer->start.tv_usec));
-}
-
-/**
- * @brief monitoring function (print elapsed)
- */
-void print_elapsed(Timer* inTimer){
-    long int elapsed = get_elapsed(inTimer);
-    printf("Elapsed : %ld us (or %f seconds)\n",elapsed,elapsed/1000000.0);
-}
-
-/**
- * @brief monitoring function (print difference :: First-Second)
- */
-void print_difference_elapsed(Timer* inTimer1,Timer*inTimer2){
-    long int diff = get_elapsed(inTimer1)-get_elapsed(inTimer2);
-    printf("Timer Difference : %ld us (or %f seconds)\n",diff,(double)diff/1000000.0);
-}
-
 
 /**
  * @brief Do everything
@@ -253,7 +206,7 @@ int main(int argc, char ** av){
 
     //Set timers
     Timer interface_timer,ref_timer;
-    int ite=0, max_ite=5;
+    int ite=0, max_ite=1;
     while(ite<max_ite){
         //Execute FMM
         tic(&interface_timer);
diff --git a/Addons/CKernelApi/Tests/testSphereElectro.c b/Addons/CKernelApi/Tests/testSphereElectro.c
index 973c4ecc..2ebb5c1b 100644
--- a/Addons/CKernelApi/Tests/testSphereElectro.c
+++ b/Addons/CKernelApi/Tests/testSphereElectro.c
@@ -7,8 +7,7 @@
 
 
 //For timing monitoring
-#include <time.h>
-#include <sys/time.h>
+#include "Timers.h"
 
 #include "../Src/CScalfmmApi.h"
 
@@ -79,10 +78,10 @@ void getNormal(double * positions, double * normeToFill){
     for(i=0 ; i<3 ; ++i){
         normeToFill[i] = positions[i]/norme;
     }
-    printf("Tgt Norme %e - %e - %e\n",
-           normeToFill[0],
-           normeToFill[1],
-           normeToFill[2]);
+    /* printf("Tgt Norme %e - %e - %e\n", */
+    /*        normeToFill[0], */
+    /*        normeToFill[1], */
+    /*        normeToFill[2]); */
 }
 
 void computeNormalXForces(int nbPoints, double * forcesToRead, double * positionsToRead, double * arrayToFill){
@@ -98,15 +97,17 @@ void computeNormalXForces(int nbPoints, double * forcesToRead, double * position
 }
 
 int main(int argc, char ** av){
-    printf("Start\n");
-    if(argc<2){
-        printf("Use : %s nb_part(cible) (optionnal : TreeHeight) \nexiting\n",av[0]);
+    omp_set_num_threads(4);
+    printf("Start %s nb_targets nb_sources \n",av[0]);
+    if(argc<3){
+        printf("Use : %s nb_part(cible) nb_part(source) (optionnal : TreeHeight) \nexiting\n",av[0]);
         exit(0);
     }
     int nbPartTarget= atoi(av[1]);
     int treeHeight = 5 ;
-    if(argc>2){
-        int treeHeight = atoi(av[2]);
+    if(argc>3){
+        int treeHeight = atoi(av[3]);
+        printf("Tree Heigth Input %d\n",treeHeight );
     }
 
     double boxWidth = 2.0;
@@ -122,13 +123,14 @@ int main(int argc, char ** av){
     memset(targetsPhiValues,0,sizeof(double)*nbPartTarget);
     //Fill
     for(i=0 ; i<nbPartTarget ; ++i){
-        targetsPhiValues[i] = -1.0;
+        targetsPhiValues[i] = 1.0;
     }
+
     generateSurfacePoints(1.0,boxCenter,nbPartTarget,targetsXYZ);
-    printf("Surface points generated \n");
+    printf("%d Surface points generated : Targets\n",nbPartTarget);
 
     //Allocation of the sources points
-    int nbPartSource = 10;
+    int nbPartSource = atoi(av[2]);
     double * sourceXYZ = malloc(sizeof(double)* 3*nbPartSource);
     double * sourcePhiValues = malloc(sizeof(double)* nbPartSource);
     //Set to Zero
@@ -138,14 +140,18 @@ int main(int argc, char ** av){
     for(i=0 ; i<nbPartSource ; ++i){
         sourcePhiValues[i] = 1.0;
     }
+
     generateInsidePoints(1.0,boxCenter,nbPartSource,sourceXYZ);
     //displayPoints(nbPartTarget,targetsXYZ);
 
-    printf("Inside points generated \n");
+    printf("%d Inside points generated :  Sources\n",nbPartSource);
+
     //displayPoints(nbPartSource,sourceXYZ);
     //Creation of arrays to store forces
-    double * arrayOfForces = malloc(sizeof(double )* 3 * (nbPartSource+nbPartTarget));
+    double * arrayOfForces = malloc(sizeof(double )* 3 * nbPartTarget);
+    double * arrayOfPot    = malloc(sizeof(double ) * (nbPartTarget));
     memset(arrayOfForces,0,sizeof(double)* 3 * (nbPartTarget));
+    memset(arrayOfPot,0,sizeof(double)* (nbPartTarget));
 
     {//Start of computation
 
@@ -157,7 +163,7 @@ int main(int argc, char ** av){
         user_descr.user_init_cell = NULL;
         user_descr.user_free_cell = NULL;
         //Set algorithm to source target
-        //scalfmm_algorithm_config(handle,source_target);
+        scalfmm_algorithm_config(handle,source_target);
         //Build the tree
         scalfmm_build_tree(handle,treeHeight, boxWidth, boxCenter, user_descr);
 
@@ -166,9 +172,8 @@ int main(int argc, char ** av){
         printf("Sources inserted \n");
         scalfmm_tree_insert_particles_xyz(handle,nbPartTarget,targetsXYZ,TARGET);
         printf("Targets inserted \n");
-        //Since we inserted first the sources, then sources will get
-        //indices from 0 to (nbPartSource-1), and targets from
-        //(nbPartSource) to nbPartSource+nbPartTarget-1).
+
+        //        scalfmm_print_everything(handle);
 
         int * arrayofIndicesSource = malloc(sizeof(int)*nbPartSource);
         int * arrayofIndicesTarget = malloc(sizeof(int)*nbPartTarget);
@@ -189,15 +194,25 @@ int main(int argc, char ** av){
 
 
         }
+        Timer fmm_timer;
+        tic(&fmm_timer);
         //Computation
         scalfmm_execute_fmm(handle/*, kernel, &my_data*/);
+        tac(&fmm_timer);
+        print_elapsed(&fmm_timer);
 
         //Get back the forces
         scalfmm_get_forces_xyz(handle,nbPartTarget,arrayOfForces,TARGET);
-        scalfmm_get_forces_xyz(handle,nbPartSource,&arrayOfForces[nbPartTarget],SOURCE);
+        //Get back the potentials, too
+        scalfmm_get_potentials(handle,nbPartTarget,arrayOfPot,TARGET);
+
+        //No need to get Source forces, since it will be 0 anyway
+        //scalfmm_get_forces_xyz(handle,nbPartSource,&arrayOfForces[nbPartTarget],SOURCE);
+
         printf("Forces computed : \n");
-        displayPoints(nbPartTarget+nbPartSource,arrayOfForces);
-        printf("As expected, Source forces are 0\n \n");
+        //Display array of forces ::
+        //displayPoints(nbPartTarget+nbPartSource,arrayOfForces);
+
         //Release memory used :
         free(arrayofIndicesSource);
         free(arrayofIndicesTarget);
@@ -207,19 +222,24 @@ int main(int argc, char ** av){
     }
 
     {//Let's check the result, we computed fr each target part its forces
-        //Storage of reference forces
+        //Storage of reference forces + potential
         double * arrayRefForces = malloc(sizeof(double)*nbPartTarget*3);
         memset(arrayRefForces,0,sizeof(double)*nbPartTarget*3);
+        double * arrayOfRefPot = malloc(sizeof(double)*nbPartTarget);
+        memset(arrayOfRefPot,0,sizeof(double)*nbPartTarget);
 
         int idTgt;
+        Timer check_timer;
+        tic(&check_timer);
         for(idTgt = 0 ; idTgt<nbPartTarget ; ++idTgt){
             int idSrc;
             double dx,dy,dz;
-            for(idSrc = 0 ; idSrc<nbPartTarget ; ++idSrc){
+
+            for(idSrc = 0 ; idSrc<nbPartSource ; ++idSrc){
                 //First compute dist.
-                dx = sourceXYZ[idSrc+0] - targetsXYZ[idTgt+0];
-                dy = sourceXYZ[idSrc+1] - targetsXYZ[idTgt+1];
-                dz = sourceXYZ[idSrc+2] - targetsXYZ[idTgt+2];
+                dx = sourceXYZ[idSrc*3+0] - targetsXYZ[idTgt*3+0];
+                dy = sourceXYZ[idSrc*3+1] - targetsXYZ[idTgt*3+1];
+                dz = sourceXYZ[idSrc*3+2] - targetsXYZ[idTgt*3+2];
 
                 //Secondly, compute coeff
                 double coeffs = targetsPhiValues[idTgt] * sourcePhiValues[idSrc];
@@ -229,24 +249,70 @@ int main(int argc, char ** av){
                 arrayRefForces[idTgt*3+0] += dx*coeffs*one_over_r3;
                 arrayRefForces[idTgt*3+1] += dy*coeffs*one_over_r3;
                 arrayRefForces[idTgt*3+2] += dz*coeffs*one_over_r3;
-
+                arrayOfRefPot[idTgt] += one_over_r * sourcePhiValues[idSrc];
             }
         }
+        tac(&check_timer);
+        print_elapsed(&check_timer);
 
         {//Then, we compare
-            double errorCumul = 0;
+
+            //For L2 norm
+            double errorCumulXSquared = 0,
+                errorCumulYSquared = 0,
+                errorCumulZSquared = 0,
+                errorCumulPotSquared = 0;
+            //For Inf Norm
+            double maxErrorX = 0,
+                maxErrorY = 0,
+                maxErrorZ = 0,
+                maxErrorPot = 0;
             int idArr;
             for(idArr = 0 ; idArr<nbPartTarget ; ++idArr){
-                errorCumul += fabs(arrayRefForces[idArr+0]-arrayOfForces[idArr+0]);
-                errorCumul += fabs(arrayRefForces[idArr+1]-arrayOfForces[idArr+1]);
-                errorCumul += fabs(arrayRefForces[idArr+2]-arrayOfForces[idArr+2]);
-                printf("Directly Computed %e %e %e\n",
-                       arrayRefForces[idArr+0],
-                       arrayRefForces[idArr+1],
-                       arrayRefForces[idArr+2]);
+                double deltaX = fabs(arrayRefForces[idArr*3+0]-arrayOfForces[idArr*3+0]);
+                double deltaY = fabs(arrayRefForces[idArr*3+1]-arrayOfForces[idArr*3+1]);
+                double deltaZ = fabs(arrayRefForces[idArr*3+2]-arrayOfForces[idArr*3+2]);
+                double deltaPot = fabs(arrayOfRefPot[idArr]-arrayOfPot[idArr]);
+
+                errorCumulXSquared += deltaX*deltaX;
+                errorCumulYSquared += deltaY*deltaY;
+                errorCumulZSquared += deltaZ*deltaZ;
+                errorCumulPotSquared += deltaPot*deltaPot;
+
+                if(maxErrorX < deltaX){ maxErrorX = deltaX ;}
+                if(maxErrorY < deltaY){ maxErrorY = deltaY ;}
+                if(maxErrorZ < deltaZ){ maxErrorZ = deltaZ ;}
+                if(maxErrorPot < deltaPot) {maxErrorPot = deltaPot;}
             }
-            printf("Error cumul : %e\n",errorCumul);
+            //Last check aim to verify the difference between energy
+            //total directly computed and energy total FMM's computed
+            double energy = 0,energyFmm =0;
+            for(idArr = 0; idArr<nbPartTarget ; ++idArr){
+                energy += arrayOfRefPot[idArr] * targetsPhiValues[idArr];
+                energyFmm += arrayOfPot[idArr] * targetsPhiValues[idArr];
+            }
+
+
+
+            printf("\n \t\t X error \t Y error \t Z error \t Pot error\n");
+            printf("Norme Sup :\t %e \t %e \t %e \t %e\n",maxErrorX,maxErrorY,maxErrorZ,maxErrorPot);
+            printf("Norme L2  :\t %e \t %e \t %e \t %e\n",
+                   sqrt(errorCumulXSquared),
+                   sqrt(errorCumulYSquared),
+                   sqrt(errorCumulZSquared),
+                   sqrt(errorCumulPotSquared));
+            printf("Norme Rms :\t %e \t %e \t %e \t %e\n",
+                   sqrt(errorCumulXSquared/((double) nbPartTarget)),
+                   sqrt(errorCumulYSquared/((double) nbPartTarget)),
+                   sqrt(errorCumulZSquared/((double) nbPartTarget)),
+                   sqrt(errorCumulPotSquared/((double) nbPartTarget)));
+            printf(" \n Energy Error  : \t %e \n Energy FMM   : \t %e \n Energy DIRECT : \t %e\n",
+                   fabs(energy-energyFmm),
+                   energyFmm,
+                   energy);
         }
+        free(arrayRefForces);
+        free(arrayOfRefPot);
     }
 
 
@@ -260,14 +326,18 @@ int main(int argc, char ** av){
 
     computeNormalXForces(nbPartTarget,targetsForces,targetsXYZ,normeXForces);
     printf("For each target, we display [Normal To Sphere] . [Force product] \n");
-    displayArray(nbPartTarget,normeXForces);
+    //displayArray(nbPartTarget,normeXForces);
 
 
     //Free memory
+    free(normeXForces);
+    free(targetsForces);
     free(sourceXYZ);
     free(sourcePhiValues);
     free(targetsXYZ);
     free(targetsPhiValues);
     free(arrayOfForces);
+    free(arrayOfPot);
+
     return EXIT_SUCCESS;
 }
-- 
GitLab