testChebInterface.c 12.1 KB
Newer Older
1 2 3
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
4
#include <math.h>
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
5

6
//For timing monitoring
7
#include "Timers.h"
8 9 10 11 12

#include "../Src/CScalfmmApi.h"

#include "../../Src/Kernels/Chebyshev/FChebInterface.h"

13

14 15 16 17 18 19 20 21
/**
 * This file is an example of the user defined kernel API, with link
 * to our Chebyshev Kernel
 **/

/**
 * @brief Wrapper to init internal ChebCell
 */
22
void* cheb_init_cell(int level, long long morton_index, int* tree_position, double* spatial_position, void * KernelDatas){
23 24 25 26 27 28 29 30 31 32 33 34 35 36
    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)
 */
37
void cheb_p2m(void* cellData, FSize nbParticlesInLeaf, const FSize* particleIndexes, void* userData){
38 39 40 41 42 43 44 45 46 47 48
    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);
}
49
void cheb_l2p(void* leafCell, FSize nbParticles, const FSize* particleIndexes, void* userData){
50 51
    ChebKernel_L2P( leafCell, nbParticles, particleIndexes, userData);
}
52 53
void cheb_p2pFull(FSize nbParticles, const FSize* particleIndexes,
                  const FSize * sourceParticleIndexes[27], FSize sourceNbPart[27],void* userData) {
54 55 56
    ChebKernel_P2P(nbParticles, particleIndexes, sourceParticleIndexes, sourceNbPart, userData);
}

57 58 59 60
void cheb_resetCell(int level, long long morton_index, int* tree_position, double* spatial_position, void * userCell){
    ChebCell_reset(level,morton_index,tree_position,spatial_position,userCell);
}

61 62 63 64 65 66 67


/**
 * @brief Do everything
 * @param number of particle (no default value)
 */
int main(int argc, char ** av){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
68
    //omp_set_num_threads(1);
69 70
    printf("Start\n");
    if(argc<2){
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
71
        printf("Use : %s nb_part (optionnal : TreeHeight) \nexiting\n",av[0]);
72 73 74
        exit(0);
    }
    int nbPart= atoi(av[1]);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
75 76 77 78 79 80
    int treeHeight = 5 ;
    if(argc>2){
        int treeHeight = atoi(av[2]);
    }


81 82 83 84
    double boxWidth = 1.0;
    double boxCenter[3];
    boxCenter[0] = boxCenter[1] = boxCenter[2] = 0.0;

85 86 87 88 89
    //Allocation of the positions and physical values
    double* particleXYZ = malloc(sizeof(double)*3*nbPart);
    double* physicalValues = malloc(sizeof(double)*nbPart);


90 91
    {
        printf("Creating Particles:\n");
92
        FSize idxPart;
93 94 95 96 97 98 99 100 101 102 103 104 105 106
        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]); */
107
        /* FSize idxPart; */
108 109 110 111 112 113 114 115 116 117
        /* 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); */
    }

118
    scalfmm_handle handle = scalfmm_init(user_defined_kernel,multi_thread);
119 120

    //For Reference
121 122
    scalfmm_handle handle_ref = scalfmm_init(chebyshev,multi_thread);

123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140

    //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);

141
    //Once is the tree built, one must set the kernel before inserting particles
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157

    //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;
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
158 159 160 161 162 163 164 165 166 167 168 169

    int nb_threads = omp_get_max_threads();
    int idThreads= 0;

    //Create as many kernels as there are threads in order to void concurrent write
    userDatas.kernelStruct = ChebKernelStruct_create(treeHeight,boxWidth,boxCenter);
    /* malloc(sizeof(void *)*nb_threads); */
    /* for(idThreads=0 ; idThreads<nb_threads ; ++idThreads){ */
    /*     userDatas.kernelStruct[idThreads] = ChebKernelStruct_create(treeHeight,boxWidth,boxCenter); // Set my kernel inside userDatas */
    /* } */

    //Only read, so no split needed
170 171
    userDatas.insertedPositions = particleXYZ;                                       // Set the position
    userDatas.myPhyValues = physicalValues;                                          // Set the physical values
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
172 173 174 175 176 177 178 179 180

    //Create as many array of forces as there are threads in order to void concurrent write
    double ** forcesToStore = malloc(sizeof(double*)*nb_threads);
    //For each array, initialisation
    for(idThreads=0 ; idThreads<nb_threads ; ++idThreads){
        forcesToStore[idThreads] =  malloc(sizeof(double)*nbPart*3); //allocate memory
        memset(forcesToStore[idThreads],0,sizeof(double)*nbPart*3);  //set to zero (IMPORTANT, as operators usually "+=" on forces)
    }
    userDatas.forcesComputed = forcesToStore;
181

182 183 184 185 186 187 188 189 190 191
    //Same idea with potential -> could be merge with loop just above,
    //but it's kept apart in order to simply
    double ** potentialToStore = malloc(sizeof(double*)*nb_threads);
    //For each array, initialisation
    for(idThreads=0 ; idThreads<nb_threads ; ++idThreads){
        potentialToStore[idThreads] = malloc(sizeof(double)* nbPart);
        memset(potentialToStore[idThreads],0,sizeof(double)* nbPart);
    }
    userDatas.potentials = potentialToStore;

192 193 194
    //Give ScalFMM the datas before calling fmm (this will set as well the kernel)
    scalfmm_user_kernel_config(handle,kernel,&userDatas);

195 196 197 198 199 200 201 202 203 204 205 206

    // Insert particles
    printf("Inserting particles...\n");
    scalfmm_tree_insert_particles_xyz(handle, nbPart, particleXYZ,BOTH);
    scalfmm_tree_insert_particles_xyz(handle_ref, nbPart, particleXYZ,BOTH);


    //Set physical values for Cheb_ref
    scalfmm_set_physical_values(handle_ref,nbPart,physicalValues,BOTH);



207 208
    //Set timers
    Timer interface_timer,ref_timer;
209
    int ite=0, max_ite=1;
210 211 212 213 214 215
    while(ite<max_ite){
        //Execute FMM
        tic(&interface_timer);
        scalfmm_execute_fmm(handle/*, kernel, &my_data*/);
        tac(&interface_timer);

216

217 218 219 220 221 222 223 224 225 226 227 228
        //Reduction on forces & potential arrays
        {
            FSize idxPart;
            for(idThreads=1 ; idThreads<nb_threads ; ++idThreads){
                for(idxPart=0 ; idxPart<nbPart ; ++idxPart){
                    //Everything is stored in first array

                    forcesToStore[0][3*idxPart+0] += forcesToStore[idThreads][3*idxPart+0];
                    forcesToStore[0][3*idxPart+1] += forcesToStore[idThreads][3*idxPart+1];
                    forcesToStore[0][3*idxPart+2] += forcesToStore[idThreads][3*idxPart+2];
                    potentialToStore[0][idxPart] += potentialToStore[idThreads][idxPart];
                }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
229 230
            }
        }
231

232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251
        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 & potentials for ref_cheb execution
        double * forcesRef     = malloc(sizeof(double)*3*nbPart);
        double * potentialsRef = malloc(sizeof(double)*nbPart);

        memset(forcesRef,0,sizeof(double)*3*nbPart);
        memset(potentialsRef,0,sizeof(double)*nbPart);

252 253 254
        scalfmm_get_forces_xyz(handle_ref,nbPart,forcesRef,BOTH);
        scalfmm_get_potentials(handle_ref,nbPart,potentialsRef,BOTH);
        //scalfmm_print_everything(handle_ref);
255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293

        {//Comparison part
            FSize idxPart;
            int nbPartOkay = 0;
            for(idxPart=0 ; idxPart<nbPart ; ++idxPart ){
                double diffX,diffY,diffZ,diffPot;
                diffX = fabs( forcesToStore[0][idxPart*3+0]-forcesRef[idxPart*3+0] );
                diffY = fabs( forcesToStore[0][idxPart*3+1]-forcesRef[idxPart*3+1] );
                diffZ = fabs( forcesToStore[0][idxPart*3+2]-forcesRef[idxPart*3+2] );
                diffPot = fabs( potentialToStore[0][idxPart]-potentialsRef[idxPart] );

                //THERE

                if(diffX < 0.00000001 && diffY < 0.00000001 && diffZ < 0.00000001 && diffPot < 0.00000001){
                    nbPartOkay++;
                }
                else{
                    printf("id : %lld : %e, %e, %e, %e, ChebInterf Pot : %e  Cheb Pot : %e \n",
                           idxPart,diffX,diffY,diffZ,diffPot,
                           potentialToStore[0][idxPart],
                           potentialsRef[idxPart]);
                }
                //That part is to verify with our usual exec' if everything is alright
                if(idxPart == 0 || idxPart == nbPart/2 || idxPart == nbPart-1){
                    printf("User one's id : %lld : %e, %e, %e, %e\n",idxPart,
                           forcesToStore[0][idxPart*3+0],
                           forcesToStore[0][idxPart*3+1],
                           forcesToStore[0][idxPart*3+2],
                           potentialToStore[0][idxPart]);
                    printf("Chebyshev one's id : %lld : %e, %e, %e, %e\n",idxPart,
                           forcesRef[idxPart*3+0],
                           forcesRef[idxPart*3+1],
                           forcesRef[idxPart*3+2],
                           potentialsRef[idxPart]);
                }
            }
            printf("End of simulation -- \t %d\n \t Percentage of good parts : %d/%d (%f %%) \n",ite,
                   nbPartOkay,nbPart,(((double) nbPartOkay)/(double)nbPart)*100);
        }
294

295 296
        free(forcesRef);
        free(potentialsRef);
297

298 299 300
        //Reset
        scalfmm_reset_tree(handle,cheb_resetCell);
        scalfmm_reset_tree(handle_ref,NULL);
301

302
        printf("Internal resets done \n");
303

304 305 306 307 308
        {//Reset User's datas
            FSize idThreads;
            for(idThreads=0;idThreads<nb_threads;++idThreads){
                memset(potentialToStore[idThreads],0,sizeof(double)*nbPart);
                memset(forcesToStore[idThreads],0,sizeof(double)*nbPart*3);
309 310
            }
        }
311 312 313
        printf("External resets done ...\n");

        ite++;
314
    }
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
315
    printf("Free the kernels\n");
316

PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
317
    printf("Free the Handles ...\n");
318 319 320 321 322
    scalfmm_dealloc_handle(handle,cheb_free_cell);
    scalfmm_dealloc_handle(handle_ref,NULL);

    free(particleXYZ);
    free(physicalValues);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
323 324 325
    //free the thread' specific datas
    for(idThreads=0 ; idThreads<nb_threads ; ++idThreads){
        free(forcesToStore[idThreads]);
326
        free(potentialToStore[idThreads]);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
327
    }
328 329
    free(forcesToStore);
    free(potentialToStore);
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
330

331 332
    ChebKernelStruct_free(userDatas.kernelStruct);
    return EXIT_SUCCESS;
333
}