Api updated, testChebInterface will execute both API internal Chebyshev and...

Api updated, testChebInterface will execute both API internal Chebyshev and API user kernel Chebyshev and display some output

Addons/CKernelApi/Src/CScalfmmApi.h

@@ -380,11 +380,20 @@ typedef void (*Callback_M2M)(int level, void* parentCell, int childPosition, voi
  * @param targetCell pointer to cell to be filled
  * @param sourceCellPosition number of source cell (source cell
  * position in the tree can inferred from its number (refer to doc))
- * @param sourceCell array of cell to be read
+ * @param sourceCell cell to be read
  * @param userData datas specific to the user's kernel
  */
 typedef void (*Callback_M2L)(int level, void* targetCell, int sourceCellPosition, void* sourceCell, void* userData);
 
+/**
+ * @brief Function to be filled by user's M2L
+ * @param level current level in the tree
+ * @param targetCell pointer to cell to be filled
+ * @param sourceCell array of cell to be read
+ * @param userData datas specific to the user's kernel
+ */
+typedef void (*Callback_M2LFull)(int level, void* targetCell, void* sourceCell[343], void* userData);
+
 /**
  * @brief Function to be filled by user's L2L
  * @param level current level in the tree
@@ -414,6 +423,19 @@ typedef void (*Callback_L2P)(void* leafCell, int nbParticles,const int* particle
  */
 typedef void (*Callback_P2P)(int nbParticles, const int* particleIndexes, int nbSourceParticles, const int* sourceParticleIndexes, void* userData);
 
+/**
+ * @brief Function to be filled by user's P2P
+ * @param nbParticles number of particle in current leaf
+ * @param particleIndexes indexes of particles currently computed
+ * @param sourceCell array of the neighbors source cells
+ * @param sourceParticleIndexes array of indices of source particles currently computed
+ * @param sourceNbPart array containing the number of part in each neighbors
+ * @param userData datas specific to the user's kernel
+ */
+typedef void (*Callback_P2PFull)(int nbParticles, const int* particleIndexes,
+                                 const int * sourceParticleIndexes[27],int sourceNbPart[27], void* userData);
+
+
 /**
  * @brief Function to be filled by user's P2P inside the leaf
  * @param nbParticles number of particle in current leaf
@@ -433,9 +455,11 @@ typedef struct User_Scalfmm_Kernel_Descriptor {
     Callback_P2M p2m;
     Callback_M2M m2m;
     Callback_M2L m2l;
+    Callback_M2LFull m2l_full;
     Callback_L2L l2l;
     Callback_L2P l2p;
     Callback_P2P p2p;
+    Callback_P2PFull p2p_full;
     Callback_P2PInner p2pinner;
 }Scalfmm_Kernel_Descriptor;
 

Addons/CKernelApi/Src/FInterEngine.hpp

@@ -187,8 +187,6 @@ public:
                     forcesToFill[indexes[idxPart]*3+0] = sources->getForcesX()[idxPart];
                     forcesToFill[indexes[idxPart]*3+1] = sources->getForcesY()[idxPart];
                     forcesToFill[indexes[idxPart]*3+2] = sources->getForcesZ()[idxPart];
-                    printf("forces found : %e,%e,%e\n",
-                           sources->getForcesX()[idxPart],sources->getForcesY()[idxPart],sources->getForcesZ()[idxPart]);
                 }
             });
     }

Addons/CKernelApi/Src/FScalFMMEngine.hpp

@@ -31,10 +31,12 @@
 #include "Components/FSimpleLeaf.hpp"
 #include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
 
-//For lagrange interpolation
-#include "Kernels/Uniform/FUnifCell.hpp"
+//For interpolation
 #include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
-#include "Kernels/Uniform/FUnifKernel.hpp"
+
+//For lagrange interpolation
+// #include "Kernels/Uniform/FUnifCell.hpp"
+// #include "Kernels/Uniform/FUnifKernel.hpp"
 
 //For chebyshev Interpolation
 #include "Kernels/Chebyshev/FChebCell.hpp"

Addons/CKernelApi/Src/FUserKernelEngine.hpp

@@ -72,6 +72,7 @@ public:
     void* getContainer() const {
         return userData;
     }
+
 };
 
 /**
@@ -115,10 +116,25 @@ public:
 
     /** Do nothing */
     virtual void M2L(CellClass* const FRestrict cell, const CellClass* interactions[], const int , const int level) {
-        if(kernel.m2l){
-            for(int idx = 0 ; idx < 343 ; ++idx){
-                if( interactions[idx] ){
-                    kernel.m2l(level, cell->getContainer(), idx, interactions[idx]->getContainer(), userData);
+        if(kernel.m2l_full){//all 343 interactions will be computed directly
+            //First, copy the fmm cell inside an array of user cells
+            void * userCellArray[343];
+            for(int i=0 ; i<343 ; ++i){
+                if(interactions[i] != nullptr){
+                    userCellArray[i] = interactions[i]->getContainer();
+                }
+                else{
+                    userCellArray[i] = nullptr;
+                }
+            }
+            kernel.m2l_full(level,cell->getContainer(),userCellArray,userData);
+        }
+        else{
+            if(kernel.m2l){
+                for(int idx = 0 ; idx < 343 ; ++idx){
+                    if( interactions[idx] ){
+                        kernel.m2l(level, cell->getContainer(), idx, interactions[idx]->getContainer(), userData);
+                    }
                 }
             }
         }
@@ -147,11 +163,27 @@ public:
                      ContainerClass* const neighbors[27], const int ){
         if(kernel.p2pinner) kernel.p2pinner(targets->getNbParticles(), targets->getIndexes().data(), userData);
 
+        if(kernel.p2p_full){
+            //Create the arrays of size and indexes
+            int nbPartPerNeighbors[27];
+            const int * indicesPerNeighbors[27];
+            for(int idx=0 ; idx<27 ; ++idx){
+                if(neighbors[idx]){
+                    nbPartPerNeighbors[idx] = neighbors[idx]->getNbParticles();
+                    indicesPerNeighbors[idx] = neighbors[idx]->getIndexes().data();
+                }
+                else{
+                    nbPartPerNeighbors[idx] = 0;
+                    indicesPerNeighbors[idx] = nullptr;
+                }
+            }
+            kernel.p2p_full(targets->getNbParticles(),targets->getIndexes().data(),indicesPerNeighbors,nbPartPerNeighbors,userData);
+        }
         if(kernel.p2p){
             for(int idx = 0 ; idx < 27 ; ++idx){
                 if( neighbors[idx] ){
                     kernel.p2p(targets->getNbParticles(), targets->getIndexes().data(),
-                                    neighbors[idx]->getNbParticles(), neighbors[idx]->getIndexes().data(), userData);
+                               neighbors[idx]->getNbParticles(), neighbors[idx]->getIndexes().data(), userData);
                 }
             }
         }
@@ -200,12 +232,14 @@ public:
 
     ~FUserKernelEngine(){
         delete octree;
-
+        octree=nullptr;
         if(arranger){
             delete arranger;
+            arranger=nullptr;
         }
         if(kernel){
             delete kernel;
+            kernel=nullptr;
         }
     }
 
@@ -233,6 +267,7 @@ public:
             octree->insert(FPoint(&XYZ[3*idPart]),idPart);
         }
         nbPart += NbPositions;
+        this->init_cell();
     }
 
     /**
@@ -389,8 +424,6 @@ public:
         }
     }
 
-
-
     /*
      * Call the user allocator on userDatas member field of each cell
      */
@@ -422,6 +455,7 @@ public:
         octree->forEachCell([&](CoreCell * currCell){
                 if(currCell->getContainer()){
                     user_cell_deallocator(currCell->getContainer());
+                    currCell->setContainer(nullptr);
                 }
             });
     }

Addons/CKernelApi/Tests/testChebInterface.c

+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+//For timing monitoring
+#include <time.h>
+#include <sys/time.h>
+
+#include "../Src/CScalfmmApi.h"
+
+#include "../../Src/Kernels/Chebyshev/FChebInterface.h"
+
+/**
+ * This file is an example of the user defined kernel API, with link
+ * to our Chebyshev Kernel
+ **/
+
+/**
+ * @brief Wrapper to init internal ChebCell
+ */
+void* cheb_init_cell(int level, long long morton_index, int* tree_position, double* spatial_position){
+    return ChebCellStruct_create(morton_index,tree_position);
+}
+
+/**
+ * @brief Wrapper to free internal ChebCell
+ */
+void cheb_free_cell(void * inCell){
+    ChebCellStruct_free(inCell);
+}
+
+/**
+ * @brief Wrapper to FMM operators (refer to CScalfmmApi.h to get the
+ * detailed descriptions)
+ */
+void cheb_p2m(void* cellData, int nbParticlesInLeaf, const int* particleIndexes, void* userData){
+    ChebKernel_P2M(cellData,nbParticlesInLeaf,particleIndexes,userData);
+}
+void cheb_m2m(int level, void* parentCell, int childPosition, void* childCell, void* userData){
+    ChebKernel_M2M(level,parentCell,childPosition,childCell,userData);
+}
+void cheb_m2l_full(int level, void* targetCell, void* sourceCell[343], void* userData){
+    ChebKernel_M2L(level, targetCell, sourceCell, userData);
+}
+void cheb_l2l(int level, void* parentCell, int childPosition, void* childCell, void* userData){
+    ChebKernel_L2L( level, parentCell, childPosition, childCell,  userData);
+}
+void cheb_l2p(void* leafCell, int nbParticles, const int* particleIndexes, void* userData){
+    ChebKernel_L2P( leafCell, nbParticles, particleIndexes, userData);
+}
+void cheb_p2pFull(int nbParticles, const int* particleIndexes,
+                  const int * sourceParticleIndexes[27], int sourceNbPart[27],void* userData) {
+    ChebKernel_P2P(nbParticles, particleIndexes, sourceParticleIndexes, sourceNbPart, userData);
+}
+
+
+/**
+ * @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
+ * @param number of particle (no default value)
+ */
+int main(int argc, char ** av){
+    omp_set_num_threads(1);
+    printf("Start\n");
+    if(argc<2){
+        printf("Use : %s nb_part\n exiting\n",av[0]);
+        exit(0);
+    }
+    int nbPart= atoi(av[1]);
+    double* particleXYZ = malloc(sizeof(double)*3*nbPart);
+    double* physicalValues = malloc(sizeof(double)*nbPart);
+
+    int treeHeight = 5;
+    double boxWidth = 1.0;
+    double boxCenter[3];
+    boxCenter[0] = boxCenter[1] = boxCenter[2] = 0.0;
+
+    {
+        printf("Creating Particles:\n");
+        int idxPart;
+        for(idxPart = 0 ; idxPart < nbPart ; ++idxPart){
+            particleXYZ[idxPart*3]   = (random()/(double)(RAND_MAX))*boxWidth - boxWidth/2 + boxCenter[0];
+            particleXYZ[idxPart*3+1] = (random()/(double)(RAND_MAX))*boxWidth - boxWidth/2 + boxCenter[1];
+            particleXYZ[idxPart*3+2] = (random()/(double)(RAND_MAX))*boxWidth - boxWidth/2 + boxCenter[2];
+            physicalValues[idxPart] = 1.0;
+
+        }
+
+    }
+
+    {//This part will write generated particles to a file at ScalFMM
+     //format in order to verify numercal results
+        /* FILE * fd = fopen("input.fma","w"); */
+        /* fprintf(fd,"8\t 4\n %d\n %f\t %f\t %f\t %f\n",nbPart,boxWidth/2.0, boxCenter[0],boxCenter[1],boxCenter[2]); */
+        /* int idxPart; */
+        /* for(idxPart=0 ; idxPart<nbPart ; ++idxPart){ */
+        /*     fprintf(fd,"%e\t %e\t %e\t %e \n", */
+        /*             particleXYZ[idxPart*3], */
+        /*             particleXYZ[idxPart*3+1], */
+        /*             particleXYZ[idxPart*3+2], */
+        /*             physicalValues[idxPart]); */
+        /* } */
+        /* fclose(fd); */
+    }
+
+    scalfmm_handle handle = scalfmm_init(user_defined_kernel);
+
+    //For Reference
+    scalfmm_handle handle_ref = scalfmm_init(chebyshev);
+
+    //Struct for user defined kernel
+    struct User_Scalfmm_Cell_Descriptor cellDescriptor;
+    cellDescriptor.user_init_cell = cheb_init_cell;
+    cellDescriptor.user_free_cell = cheb_free_cell;
+
+    //Struct for ref cheb kernel
+    struct User_Scalfmm_Cell_Descriptor user_descr;
+    user_descr.user_init_cell = NULL;
+    user_descr.user_free_cell = NULL;
+
+
+    // Init tree and cell
+    printf("Building the tree and Initizalizing the cells:\n");
+
+    scalfmm_build_tree(handle,treeHeight, boxWidth, boxCenter, cellDescriptor);
+    scalfmm_build_tree(handle_ref,treeHeight, boxWidth, boxCenter, user_descr);
+
+    // Insert particles
+    printf("Inserting particles...\n");
+    scalfmm_tree_insert_particles_xyz(handle, nbPart, particleXYZ);
+    scalfmm_tree_insert_particles_xyz(handle_ref, nbPart, particleXYZ);
+
+    //Set physical values for Cheb_ref
+    scalfmm_set_physical_values(handle_ref,nbPart,physicalValues);
+    //Set our callbacks
+    struct User_Scalfmm_Kernel_Descriptor kernel;
+    kernel.p2m = cheb_p2m;
+    kernel.m2m = cheb_m2m;
+    //init the other to NULL
+    kernel.m2l = NULL;
+    kernel.m2l_full = cheb_m2l_full;
+    kernel.l2l = cheb_l2l;
+    kernel.l2p = cheb_l2p;
+    kernel.p2p_full = cheb_p2pFull;
+    kernel.p2pinner = NULL;
+    kernel.p2p = NULL;
+
+    //Set my datas
+    UserData userDatas;
+    userDatas.kernelStruct = ChebKernelStruct_create(treeHeight,boxWidth,boxCenter); // Set my kernel inside userDatas
+    userDatas.insertedPositions = particleXYZ;                                       // Set the position
+    userDatas.myPhyValues = physicalValues;                                          // Set the physical values
+    double * forcesToStore =  malloc(sizeof(double)*nbPart*3);                       // Create array to store results
+    memset(forcesToStore,0,sizeof(double)*nbPart*3);                                 // memset array
+    userDatas.forcesComputed = forcesToStore;                                        // Set the forces array
+
+    //Give ScalFMM the datas before calling fmm (this will set as well the kernel)
+    scalfmm_user_kernel_config(handle,kernel,&userDatas);
+
+    //Set timers
+    Timer interface_timer,ref_timer;
+
+    //Execute FMM
+    tic(&interface_timer);
+    scalfmm_execute_fmm(handle/*, kernel, &my_data*/);
+    tac(&interface_timer);
+
+    printf("User defined Chebyshev done\n");
+    print_elapsed(&interface_timer);
+
+    tic(&ref_timer);
+    scalfmm_execute_fmm(handle_ref/*, kernel, &my_data*/);
+    tac(&ref_timer);
+
+    printf("Intern Chebyshev done\n");
+    print_elapsed(&ref_timer);
+
+    //Print time results
+    print_difference_elapsed(&interface_timer,&ref_timer);
+
+    //get back the forces for ref_cheb execution
+    double * forcesRef = malloc(sizeof(double)*3*nbPart);
+    memset(forcesRef,0,sizeof(double)*3*nbPart);
+    scalfmm_get_forces_xyz(handle_ref,nbPart,forcesRef);
+    {
+        int idxPart;
+        int nbPartOkay = 0;
+        for(idxPart=0 ; idxPart<nbPart ; ++idxPart ){
+            double diffX,diffY,diffZ;
+            diffX = userDatas.forcesComputed[idxPart*3+0]-forcesRef[idxPart*3+0];
+            diffY = userDatas.forcesComputed[idxPart*3+1]-forcesRef[idxPart*3+1];
+            diffZ = userDatas.forcesComputed[idxPart*3+2]-forcesRef[idxPart*3+2];
+            /* printf("id : %d : %e, %e, %e\n",idxPart,diffX,diffY,diffZ); */
+            if(diffX < 0.00000001 || diffY < 0.00000001 || diffZ < 0.00000001){
+                nbPartOkay++;
+            }
+            //That part is to verify with our usual exec' if everything is alright
+            if(idxPart == 0 || idxPart == nbPart/2 || idxPart == nbPart-1){
+                printf("id : %d : %e, %e, %e\n",idxPart,userDatas.forcesComputed[idxPart*3+0],
+                       userDatas.forcesComputed[idxPart*3+1],
+                       userDatas.forcesComputed[idxPart*3+2]);
+            }
+        }
+        printf("End of simulation \n \t Percentage of good parts : %d/%d (%f %%) \n",
+               nbPartOkay,nbPart,(((double) nbPartOkay)/(double)nbPart)*100);
+    }
+
+    printf("Free the Handle ...\n");
+    scalfmm_dealloc_handle(handle,cheb_free_cell);
+    scalfmm_dealloc_handle(handle_ref,NULL);
+
+    free(particleXYZ);
+    free(physicalValues);
+    free(userDatas.forcesComputed);
+}

Src/Kernels/Chebyshev/FChebInterface.h

+// ===================================================================================
+// Copyright ScalFmm 2011 INRIA,
+// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
+// This software is a computer program whose purpose is to compute the FMM.
+//
+// This software is governed by the CeCILL-C and LGPL licenses and
+// abiding by the rules of distribution of free software.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public and CeCILL-C Licenses for more details.
+// "http://www.cecill.info".
+// "http://www.gnu.org/licenses".
+// ===================================================================================
+
+
+/**
+ * This file provide an interface to the Chebyshev kernel, in order to
+ * call it from C code (and thus use it through API's user defined
+ * kernel feature).
+ */
+
+#ifndef FCHEBINTERFACE_H
+#define FCHEBINTERFACE_H
+
+
+
+//To access a cell
+struct FChebCell_struct;
+typedef struct FChebCell_struct ChebCellStruct;
+ChebCellStruct * ChebCellStruct_create(long long int index,int * tree_position);
+void ChebCellStruct_free(ChebCellStruct * cell);
+
+/* //To access a Leaf */
+/* struct FSimpleLeaf_struct; */
+/* typedef struct FSimpleLeaf_struct ChebLeafStruct; */
+/* ChebLeafStruct * ChebLeafStruct_create(); */
+/* void ChebLeafStruct_free(ChebLeafStruct* leaf); */
+
+//To access the kernel
+struct FChebKernel_struct;
+typedef struct FChebKernel_struct ChebKernelStruct;
+ChebKernelStruct * ChebKernelStruct_create(int inTreeHeight,
+                                           double inBoxWidth,
+                                           double* inBoxCenter);
+
+void ChebKernelStruct_free(void * kernel);
+//To access kernel member function
+
+void ChebKernel_P2M(void * leafCell, int nbParticles,const int* particleIndexes, void* kernel);
+void ChebKernel_M2M(int level, void* parentCell, int childPosition, void* childCell, void* kernel);
+void ChebKernel_M2L(int level, void* targetCell, void* sourceCell[343], void* kernel);
+void ChebKernel_L2L(int level, void* parentCell, int childPosition, void* childCell, void* kernel);
+void ChebKernel_L2P(void* leafCell, int nbParticles, const int* particleIndexes, void* kernel);
+void ChebKernel_P2P(int nbParticles, const int* particleIndexes,
+                    const int * sourceParticleIndexes[27], int sourceNbPart[27],void* userData);
+
+typedef struct myUserDatas{
+    ChebKernelStruct * kernelStruct;
+    double * insertedPositions;
+    double * myPhyValues;
+    double * forcesComputed;
+}UserData;
+
+
+#endif //FCHEBINTERFACE_H