testSphereElectro.c 12 KB
Newer Older
1 2 3 4 5 6 7 8 9
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>




//For timing monitoring
10
#include "Timers.h"
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80

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

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

double getRandom(){
    return (random()/(double)(RAND_MAX));
}

void generateSurfacePointOnUnitSphere(int N, double * points){
    double u, v, theta, phi, sinPhi ;
    //
    int j = 0,i=0 ;
    for ( i = 0 ; i< N ; ++i, j+=3)  {
        //
        u = getRandom() ;  v = getRandom() ;
        theta  = 2*M_PI*u ;
        phi     = acos(2*v-1);
        sinPhi = sin(phi);
        //
        points[j]   =     cos(theta)*sinPhi ;
        points[j+1] =   sin(theta)*sinPhi ;
        points[j+2] =   2*v-1 ;
        //
    }
}

void generateSurfacePoints(double rayon, double* centre, int nbDePoints, double* points){

    generateSurfacePointOnUnitSphere(nbDePoints , points) ;
    int j =0,i=0 ;
    for ( i = 0 ; i< nbDePoints ; ++i, j+=3)  {
        points[j]    *= rayon + centre[0];
        points[j+1]  *= rayon + centre[1];
        points[j+2]  *= rayon + centre[2];
    }
}

void generateInsidePoints(double rayon, double*centre, int nbDePoints,  double* points){
    generateSurfacePointOnUnitSphere(nbDePoints, points);
    int j=0;
    double u;
    for(j=0 ; j<nbDePoints ; ++j){
        u = getRandom();
        points[j]    *= (rayon + centre[0])*u;
        points[j+1]  *= (rayon + centre[1])*u;
        points[j+2]  *= (rayon + centre[2])*u;
    }
}

void displayPoints(int nbPoints, double * points){
    int i = 0;
    for(i=0 ; i<nbPoints ; ++i){
        printf("%e %e %e \n",points[i*3],points[i*3+1],points[i*3+2]);
    }
}

void displayArray(int nbValue, double * array){
    int i = 0;
    for(i=0 ; i<nbValue ; ++i){
        printf("%e \n",array[i]);
    }
}

void getNormal(double * positions, double * normeToFill){
    int i;
    double norme = sqrt(positions[0]*positions[0] + positions[1]*positions[1] + positions[2]*positions[2]);
    for(i=0 ; i<3 ; ++i){
        normeToFill[i] = positions[i]/norme;
    }
81 82 83 84
    /* printf("Tgt Norme %e - %e - %e\n", */
    /*        normeToFill[0], */
    /*        normeToFill[1], */
    /*        normeToFill[2]); */
85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
}

void computeNormalXForces(int nbPoints, double * forcesToRead, double * positionsToRead, double * arrayToFill){
    double * currentNormal = malloc(sizeof(double)*3);
    int idxPart,i;
    for(idxPart = 0 ; idxPart<nbPoints ; ++idxPart){
        getNormal(&positionsToRead[idxPart],currentNormal); //get the norme
        for(i=0 ; i<3 ; ++i){
            arrayToFill[idxPart] += currentNormal[i]*forcesToRead[idxPart+i];
        }
    }
    free(currentNormal);
}

int main(int argc, char ** av){
100 101 102 103
    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]);
104 105 106 107
        exit(0);
    }
    int nbPartTarget= atoi(av[1]);
    int treeHeight = 5 ;
108 109 110
    if(argc>3){
        int treeHeight = atoi(av[3]);
        printf("Tree Heigth Input %d\n",treeHeight );
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
    }

    double boxWidth = 2.0;
    double boxCenter[3];
    boxCenter[0] = boxCenter[1] = boxCenter[2] = 0.0;

    int i;
    //Allocation of the target points
    double * targetsXYZ = malloc(sizeof(double)* 3*nbPartTarget);
    double * targetsPhiValues = malloc(sizeof(double)* nbPartTarget);
    //Memset (au cas ou)
    memset(targetsXYZ,0,sizeof(double)*3*nbPartTarget);
    memset(targetsPhiValues,0,sizeof(double)*nbPartTarget);
    //Fill
    for(i=0 ; i<nbPartTarget ; ++i){
126
        targetsPhiValues[i] = 1.0;
127
    }
128

129
    generateSurfacePoints(1.0,boxCenter,nbPartTarget,targetsXYZ);
130
    printf("%d Surface points generated : Targets\n",nbPartTarget);
131 132

    //Allocation of the sources points
133
    int nbPartSource = atoi(av[2]);
134 135 136 137 138 139 140 141 142
    double * sourceXYZ = malloc(sizeof(double)* 3*nbPartSource);
    double * sourcePhiValues = malloc(sizeof(double)* nbPartSource);
    //Set to Zero
    memset(sourceXYZ,0,3*sizeof(double)*nbPartSource);
    memset(sourcePhiValues,0,sizeof(double)*nbPartSource);
    //Fill
    for(i=0 ; i<nbPartSource ; ++i){
        sourcePhiValues[i] = 1.0;
    }
143

144 145 146
    generateInsidePoints(1.0,boxCenter,nbPartSource,sourceXYZ);
    //displayPoints(nbPartTarget,targetsXYZ);

147 148
    printf("%d Inside points generated :  Sources\n",nbPartSource);

149 150
    //displayPoints(nbPartSource,sourceXYZ);
    //Creation of arrays to store forces
151 152
    double * arrayOfForces = malloc(sizeof(double )* 3 * nbPartTarget);
    double * arrayOfPot    = malloc(sizeof(double ) * (nbPartTarget));
153
    memset(arrayOfForces,0,sizeof(double)* 3 * (nbPartTarget));
154
    memset(arrayOfPot,0,sizeof(double)* (nbPartTarget));
155 156 157 158 159 160 161 162 163 164 165

    {//Start of computation

        //For handling the library
        scalfmm_handle handle = scalfmm_init(chebyshev,source_target);

        //Struct for ref cheb kernel
        struct User_Scalfmm_Cell_Descriptor user_descr;
        user_descr.user_init_cell = NULL;
        user_descr.user_free_cell = NULL;
        //Set algorithm to source target
166
        scalfmm_algorithm_config(handle,source_target);
167 168 169 170 171 172 173 174
        //Build the tree
        scalfmm_build_tree(handle,treeHeight, boxWidth, boxCenter, user_descr);

        //Insert Sources and targets
        scalfmm_tree_insert_particles_xyz(handle,nbPartSource,sourceXYZ,SOURCE);
        printf("Sources inserted \n");
        scalfmm_tree_insert_particles_xyz(handle,nbPartTarget,targetsXYZ,TARGET);
        printf("Targets inserted \n");
175 176

        //        scalfmm_print_everything(handle);
177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196

        int * arrayofIndicesSource = malloc(sizeof(int)*nbPartSource);
        int * arrayofIndicesTarget = malloc(sizeof(int)*nbPartTarget);
        {//Set physical values

            //SRC
            int idPart;
            for(idPart = 0 ; idPart<nbPartSource ; ++idPart){
                arrayofIndicesSource[idPart] = idPart;
            }
            scalfmm_set_physical_values_npart(handle,nbPartSource,arrayofIndicesSource,sourcePhiValues,SOURCE);
            //TGT
            for(idPart = 0 ; idPart<nbPartTarget ; ++idPart){
                arrayofIndicesTarget[idPart] = idPart; // here, we add the number of sources previously inserted
            }
            scalfmm_set_physical_values_npart(handle,nbPartTarget,arrayofIndicesTarget,targetsPhiValues,TARGET);



        }
197 198
        Timer fmm_timer;
        tic(&fmm_timer);
199 200
        //Computation
        scalfmm_execute_fmm(handle/*, kernel, &my_data*/);
201 202
        tac(&fmm_timer);
        print_elapsed(&fmm_timer);
203 204 205

        //Get back the forces
        scalfmm_get_forces_xyz(handle,nbPartTarget,arrayOfForces,TARGET);
206 207 208 209 210 211
        //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);

212
        printf("Forces computed : \n");
213 214 215
        //Display array of forces ::
        //displayPoints(nbPartTarget+nbPartSource,arrayOfForces);

216 217 218 219 220 221 222 223 224
        //Release memory used :
        free(arrayofIndicesSource);
        free(arrayofIndicesTarget);

        scalfmm_dealloc_handle(handle,NULL);

    }

    {//Let's check the result, we computed fr each target part its forces
225
        //Storage of reference forces + potential
226 227
        double * arrayRefForces = malloc(sizeof(double)*nbPartTarget*3);
        memset(arrayRefForces,0,sizeof(double)*nbPartTarget*3);
228 229
        double * arrayOfRefPot = malloc(sizeof(double)*nbPartTarget);
        memset(arrayOfRefPot,0,sizeof(double)*nbPartTarget);
230 231

        int idTgt;
232 233
        Timer check_timer;
        tic(&check_timer);
234 235 236
        for(idTgt = 0 ; idTgt<nbPartTarget ; ++idTgt){
            int idSrc;
            double dx,dy,dz;
237 238

            for(idSrc = 0 ; idSrc<nbPartSource ; ++idSrc){
239
                //First compute dist.
240 241 242
                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];
243 244 245 246 247 248 249 250 251

                //Secondly, compute coeff
                double coeffs = targetsPhiValues[idTgt] * sourcePhiValues[idSrc];
                double one_over_r = 1.0/(sqrt(dx*dx+dy*dy+dz*dz));
                double one_over_r3 = one_over_r * one_over_r * one_over_r;

                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;
252
                arrayOfRefPot[idTgt] += one_over_r * sourcePhiValues[idSrc];
253 254
            }
        }
255 256
        tac(&check_timer);
        print_elapsed(&check_timer);
257 258

        {//Then, we compare
259 260 261 262 263 264 265 266 267 268 269

            //For L2 norm
            double errorCumulXSquared = 0,
                errorCumulYSquared = 0,
                errorCumulZSquared = 0,
                errorCumulPotSquared = 0;
            //For Inf Norm
            double maxErrorX = 0,
                maxErrorY = 0,
                maxErrorZ = 0,
                maxErrorPot = 0;
270 271
            int idArr;
            for(idArr = 0 ; idArr<nbPartTarget ; ++idArr){
272 273 274 275 276 277 278 279 280 281 282 283 284 285
                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;}
286
            }
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312
            //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);
313
        }
314 315
        free(arrayRefForces);
        free(arrayOfRefPot);
316 317 318 319 320 321 322 323 324 325 326 327 328
    }


    //Part where we apply normal on target's forces vector
    //Copying each target's parts forces,
    double * targetsForces = malloc(sizeof(double) * 3 * nbPartTarget);
    memcpy(targetsForces,arrayOfForces,sizeof(double)*3*nbPartTarget);

    double * normeXForces =  malloc(sizeof(double) * nbPartTarget);
    memset(normeXForces,0,sizeof(double) * nbPartTarget);

    computeNormalXForces(nbPartTarget,targetsForces,targetsXYZ,normeXForces);
    printf("For each target, we display [Normal To Sphere] . [Force product] \n");
329
    //displayArray(nbPartTarget,normeXForces);
330 331 332


    //Free memory
333 334
    free(normeXForces);
    free(targetsForces);
335 336 337 338 339
    free(sourceXYZ);
    free(sourcePhiValues);
    free(targetsXYZ);
    free(targetsPhiValues);
    free(arrayOfForces);
340 341
    free(arrayOfPot);

342 343
    return EXIT_SUCCESS;
}